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
|
/*
* intensities.c
*
* Extract integrated intensities by relrod estimation
*
* (c) 2007 Thomas White <taw27@cam.ac.uk>
*
* dtr - Diffraction Tomography Reconstruction
*
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include "control.h"
#include "reflections.h"
#include "image.h"
#include "reproject.h"
/* Extract integrated reflection intensities by estimating the spike function
* based on the observed intensity and the calculated excitation error from
* the lattice refinement. Easy. */
void intensities_extract(ControlContext *ctx) {
int i, j;
int n_meas, n_dupl, n_notf;
double max;
Reflection *reflection;
/* Free previous analysis if required */
if ( ctx->integrated != NULL ) {
reflectionlist_free(ctx->integrated);
}
ctx->integrated = reflectionlist_new();
n_meas = 0;
n_dupl = 0;
n_notf = 0;
max = 0;
for ( i=0; i<ctx->images->n_images; i++ ) {
ImageRecord *image;
image = &ctx->images->images[i];
if ( image->rflist == NULL ) image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
for ( j=0; j<image->rflist->n_features; j++ ) {
ImageFeature *feature;
signed int h, k, l;
feature = &image->rflist->features[j];
h = feature->reflection->h;
k = feature->reflection->k;
l = feature->reflection->l;
if ( feature->partner != NULL ) {
if ( (h!=0) || (k!=0) || (l!=0) ) {
double intensity;
Reflection *new;
/* Perform relrod calculation of doom here.
* TODO: Figure out if this is even possible. */
intensity = feature->partner->intensity;
new = reflection_add(ctx->integrated,
feature->reflection->x, feature->reflection->y, feature->reflection->z,
intensity, REFLECTION_GENERATED);
if ( new != NULL ) {
new->h = h;
new->k = k;
new->l = l;
//printf("IN: Adding %3i %3i %3i, intensity=%f\n", h, k, l, intensity);
if ( intensity > max ) max = intensity;
n_meas++;
} else {
printf("IN: Duplicate measurement for %3i %3i %3i\n", h, k, l);
n_dupl++;
}
}
} else {
//printf("IN: %3i %3i %3i not found\n", h, k, l);
n_notf++;
}
}
}
/* Normalise all reflections to the most intense reflection */
reflection = ctx->integrated->reflections;
while ( reflection ) {
reflection->intensity /= max;
reflection = reflection->next;
}
printf("IN: %5i intensities measured\n", n_meas);
printf("IN: %5i duplicated measurements\n", n_dupl);
printf("IN: %5i predicted reflections not found\n", n_notf);
}
|
