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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
/*
* ipr.c
*
* Iterative Prediction-Refinement Reconstruction
*
* (c) 2007 Thomas White <taw27@cam.ac.uk>
* Gordon Ball <gfb21@cam.ac.uk>
*
* dtr - Diffraction Tomography Reconstruction
*
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "control.h"
#include "reflections.h"
#include "itrans.h"
#include "utils.h"
#include "imagedisplay.h"
#include "reproject.h"
#include "ipr.h"
#include "displaywindow.h"
#include "basis.h"
/* Generate list of reflections to analyse */
static ReflectionList *ipr_generate(ControlContext *ctx, Basis *basis) {
ReflectionList *ordered;
double max_res;
signed int h, k, l;
int max_order_a, max_order_b, max_order_c;
/* NB This assumes the diffraction pattern is "vaguely" centered... */
if ( ctx->inputfiletype != INPUT_CACHE) {
int max_width, max_height;
max_width = ctx->images[0].width;
max_height = ctx->images[0].height;
if ( ctx->fmode == FORMULATION_PIXELSIZE ) {
max_res = sqrt(max_width*max_width + max_height*max_height) * ctx->images[0].pixel_size;
} else {
max_res = sqrt(max_width*max_width + max_height*max_height) / (ctx->images[0].lambda * ctx->images[0].camera_len);
}
max_res = max_res / 2;
} else {
max_res = 1e10; //works until I put some more in the reflect.cache header
}
ordered = reflectionlist_new();
max_order_a = max_res/modulus(basis->a.x, basis->a.y, basis->a.z);
max_order_b = max_res/modulus(basis->b.x, basis->b.y, basis->b.z);
max_order_c = max_res/modulus(basis->c.x, basis->c.y, basis->c.z);
for ( h=-max_order_a; h<=max_order_a; h++ ) {
for ( k=-max_order_b; k<=max_order_b; k++ ) {
for ( l=-max_order_c; l<=max_order_c; l++ ) {
double x, y, z;
x = h*basis->a.x + k*basis->b.x + l*basis->c.x;
y = h*basis->a.y + k*basis->b.y + l*basis->c.y;
z = h*basis->a.z + k*basis->b.z + l*basis->c.z;
if ( ( x*x + y*y + z*z ) <= max_res*max_res ) {
Reflection *ref;
ref = reflection_add(ordered, x, y, z, 1, REFLECTION_GENERATED);
if ( ref ) { /* Check it's not NULL */
ref->h = h; ref->k = k; ref->l = l;
}
if ( (h!=0) || (k!=0) || (l!=0) ) {
// reflection_add(ctx->reflectionlist, x, y, z, 1, REFLECTION_GENERATED);
reflection_add(ordered, x, y, z, 1, REFLECTION_GENERATED);
}
}
}
}
}
return ordered;
}
static gint ipr_clicked(GtkWidget *widget, GdkEventButton *event, ControlContext *ctx) {
ImageReflection *refl;
size_t n, j;
ctx->ipr_cur_image++;
if ( ctx->ipr_cur_image == ctx->n_images ) ctx->ipr_cur_image = 0;
imagedisplay_clear_circles(ctx->ipr_id);
reflection_clear_markers(ctx->reflectionlist);
refl = reproject_get_reflections(ctx->images[ctx->ipr_cur_image], &n, ctx->ipr_lat, ctx);
for ( j=0; j<n; j++ ) {
imagedisplay_mark_circle(ctx->ipr_id, refl[j].x, refl[j].y);
}
imagedisplay_put_data(ctx->ipr_id, ctx->images[ctx->ipr_cur_image]);
displaywindow_update(ctx->dw);
return 0;
}
int ipr_refine(ControlContext *ctx) {
Basis *basis;
int finished;
basis = basis_find(ctx);
if ( !basis ) {
printf("IP: Unable to find basis\n");
return -1;
}
ctx->ipr_basis = basis;
ctx->ipr_lat = ipr_generate(ctx, basis);
printf("IP: Performing refinement...\n");
finished = 0;
do {
size_t n, j;
ImageReflection *refl;
/* Select an image */
ctx->ipr_cur_image = 0;
ctx->ipr_id = imagedisplay_open_with_message(ctx->images[ctx->ipr_cur_image], "Current Image", "Click to change image",
IMAGEDISPLAY_SHOW_CENTRE | IMAGEDISPLAY_SHOW_TILT_AXIS, G_CALLBACK(ipr_clicked), ctx);
refl = reproject_get_reflections(ctx->images[ctx->ipr_cur_image], &n, ctx->ipr_lat, ctx);
for ( j=0; j<n; j++ ) {
imagedisplay_mark_circle(ctx->ipr_id, refl[j].x, refl[j].y);
}
finished = 1;
} while ( !finished );
return 0;
}
|
