1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
/*
* hrs-scaling.c
*
* Intensity scaling using generalised HRS target function
*
* (c) 2006-2010 Thomas White <taw@physics.org>
*
* Part of CrystFEL - crystallography with a FEL
*
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <stdlib.h>
#include <assert.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include "image.h"
#include "peaks.h"
#include "symmetry.h"
#include "geometry.h"
#include "cell.h"
static void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r)
{
int i, j;
for ( i=0; i<r; i++ ) {
STATUS("[ ");
for ( j=0; j<r; j++ ) {
STATUS("%+9.3e ", gsl_matrix_get(M, i, j));
}
STATUS("][ a%2i ] = [ %+9.3e ]\n", i, gsl_vector_get(v, i));
}
}
/* Scale the stack of images */
void scale_intensities(struct image *images, int n, const char *sym)
{
#if 0
gsl_matrix *M;
gsl_vector *v;
gsl_vector *shifts;
int j;
M = gsl_matrix_calloc(n, n);
v = gsl_vector_calloc(n);
for ( j=0; j<n; j++ ) {
signed int hind, kind, lind;
signed int ha, ka, la;
double I_full, delta_I;
float I_partial;
float xc, yc;
int h;
struct image *image = &images[j];
struct cpeak *spots = image->cpeaks;
for ( h=0; h<image->n_cpeaks; h++ ) {
int g;
double v_c, gr, I_full;
hind = spots[h].h;
kind = spots[h].k;
lind = spots[h].l;
/* Don't attempt to use spots with very small
* partialities, since it won't be accurate. */
if ( spots[h].p < 0.1 ) continue;
if ( integrate_peak(image, spots[h].x, spots[h].y,
&xc, &yc, &I_partial, NULL, NULL, 1, 1) ) {
continue;
}
get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
I_full = lookup_intensity(i_full, ha, ka, la);
delta_I = I_partial - (spots[h].p * I_full / image->osf);
for ( g=0; g<NUM_PARAMS; g++ ) {
double M_curr, M_c;
M_curr = gsl_matrix_get(M, g, k);
M_c = gradient(image, g, spots[h],
image->profile_radius)
* gradient(image, k, spots[h],
image->profile_radius);
M_c *= pow(I_full, 2.0);
gsl_matrix_set(M, g, k, M_curr + M_c);
}
gr = gradient(image, k, spots[h],
image->profile_radius);
v_c = delta_I * I_full * gr;
gsl_vector_set(v, k, v_c);
}
}
show_matrix_eqn(M, v, NUM_PARAMS);
shifts = gsl_vector_alloc(NUM_PARAMS);
gsl_linalg_HH_solve(M, v, shifts);
for ( param=0; param<NUM_PARAMS; param++ ) {
double shift = gsl_vector_get(shifts, param);
apply_shift(image, param, shift);
}
gsl_matrix_free(M);
gsl_vector_free(v);
gsl_vector_free(shifts);
free(spots);
spots = find_intersections(image, image->indexed_cell, n, 0);
*pspots = spots;
return mean_partial_dev(image, spots, *n, sym, i_full, NULL);
#endif
}
|
