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
152
153
|
/*
* mapping.c
*
* 3D Mapping
*
* (c) 2007 Thomas White <taw27@cam.ac.uk>
*
* dtr - Diffraction Tomography Reconstruction
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "reflections.h"
#include "control.h"
#include "itrans.h"
#include "image.h"
#include "displaywindow.h"
static int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {
/* "Input" space */
double d, x, y;
ImageRecord *imagerecord;
/* Angular description of reflection */
double theta, psi, k;
/* Reciprocal space */
double tilt;
double omega;
double x_temp, y_temp, z_temp;
double nx, ny, nz;
imagerecord = refl->parent;
x = refl->x; y = refl->y;
tilt = 2*M_PI*(imagerecord->tilt/360); /* Convert to Radians */
omega = 2*M_PI*(imagerecord->omega/360); /* Likewise */
k = 1/imagerecord->lambda;
/* Calculate an angular description of the reflection */
if ( imagerecord->fmode == FORMULATION_CLEN ) {
x /= imagerecord->resolution;
y /= imagerecord->resolution; /* Convert pixels to metres */
d = sqrt((x*x) + (y*y));
theta = atan2(d, imagerecord->camera_len);
} else if (imagerecord->fmode == FORMULATION_PIXELSIZE ) {
x *= imagerecord->pixel_size;
y *= imagerecord->pixel_size; /* Convert pixels to metres^-1 */
d = sqrt((x*x) + (y*y));
theta = atan2(d, k);
} else {
fprintf(stderr, "Unrecognised formulation mode in reflection_add_from_dp\n");
return -1;
}
psi = atan2(y, x);
x_temp = k*sin(theta)*cos(psi);
y_temp = k*sin(theta)*sin(psi);
z_temp = k - k*cos(theta);
/* Apply the rotations...
First: rotate image clockwise until tilt axis is aligned horizontally. */
nx = x_temp*cos(omega) + y_temp*sin(omega);
ny = -x_temp*sin(omega) + y_temp*cos(omega);
nz = z_temp;
/* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way.
This is because the crystal is rotated in the experiment, not the Ewald sphere. */
x_temp = nx; y_temp = ny; z_temp = nz;
nx = x_temp;
ny = cos(tilt)*y_temp + sin(tilt)*z_temp;
nz = -sin(tilt)*y_temp + cos(tilt)*z_temp;
/* Finally, reverse the omega rotation to restore the location of the image in 3D space */
x_temp = nx; y_temp = ny; z_temp = nz;
nx = x_temp*cos(-omega) + y_temp*sin(-omega);
ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
nz = z_temp;
*ddx = nx;
*ddy = ny;
*ddz = nz;
*twotheta = theta; /* Radians. I've used the "wrong" nomenclature above */
return 0;
}
static void mapping_map_features(ControlContext *ctx) {
int i;
/* Create reflection list for measured reflections */
if ( ctx->reflectionlist ) reflectionlist_free(ctx->reflectionlist);
ctx->reflectionlist = reflectionlist_new();
printf("MP: Mapping to 3D..."); fflush(stdout);
for ( i=0; i<ctx->images->n_images; i++ ) {
int j;
/* Iterate over the features in this image */
for ( j=0; j<ctx->images->images[i].features->n_features; j++ ) {
double nx, ny, nz, twotheta;
if ( !mapping_map_to_space(&ctx->images->images[i].features->features[j], &nx, &ny, &nz, &twotheta) ) {
reflection_add(ctx->reflectionlist, nx, ny, nz, ctx->images->images[i].features->features[j].intensity, REFLECTION_NORMAL);
}
}
}
printf("done.\n");
}
ReflectionList *mapping_create(ControlContext *ctx) {
int i;
/* Find all the features */
printf("MP: Analysing images..."); fflush(stdout);
for ( i=0; i<ctx->images->n_images; i++ ) {
ctx->images->images[i].features = itrans_process_image(&ctx->images->images[i], ctx->psmode);
}
printf("done.\n");
mapping_map_features(ctx);
return ctx->reflectionlist;
}
void mapping_adjust_axis(ControlContext *ctx, double offset) {
int i;
for ( i=0; i<ctx->images->n_images; i++ ) {
printf("Image #%3i: old omega=%f deg, new omega=%f deg\n", i, ctx->images->images[i].omega, ctx->images->images[i].omega+offset);
ctx->images->images[i].omega += offset;
}
mapping_map_features(ctx);
if ( ctx->dw ) displaywindow_update(ctx->dw);
}
|
