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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
|
from __future__ import annotations
import math
import random
from dxtbx.serialize import load
from scitbx import matrix
from scitbx.array_family import flex
from dials.algorithms.profile_model.gaussian_rs import (
BBoxCalculator3D,
CoordinateSystem,
)
from dials.algorithms.profile_model.gaussian_rs.transform import (
MapFramesForward,
MapFramesReverse,
)
def test_map_frames_forward(dials_data):
sequence = load.imageset(
dials_data("centroid_test_data", pathlib=True) / "sweep.json"
)
# Get the models
beam = sequence.get_beam()
detector = sequence.get_detector()
gonio = sequence.get_goniometer()
scan = sequence.get_scan()
# Set the delta_divergence/mosaicity
n_sigma = 3
sigma_divergence = 0.060 * math.pi / 180
mosaicity = 0.154 * math.pi / 180
delta_divergence = n_sigma * sigma_divergence
delta_mosaicity = n_sigma * mosaicity
# Set the grid size
grid_size = (4, 4, 4)
# Create the E3 fraction object
transform = MapFramesForward(
scan.get_array_range()[0],
scan.get_oscillation(deg=False)[0],
scan.get_oscillation(deg=False)[1],
mosaicity,
n_sigma,
grid_size[2],
)
# Create the bounding box calculator
calculate_bbox = BBoxCalculator3D(
beam, detector, gonio, scan, delta_divergence, delta_mosaicity
)
assert len(detector) == 1
s0 = beam.get_s0()
m2 = gonio.get_rotation_axis()
s0_length = matrix.col(beam.get_s0()).length()
for i in range(100):
# Get random x, y, z
x = random.uniform(0, 2000)
y = random.uniform(0, 2000)
z = random.uniform(0, 9)
# Get random s1, phi, panel
s1 = matrix.col(detector[0].get_pixel_lab_coord((x, y))).normalize() * s0_length
phi = scan.get_angle_from_array_index(z, deg=False)
panel = 0
# Calculate the bounding box
bbox = calculate_bbox(s1, z, panel)
# Create the XDS coordinate system
xcs = CoordinateSystem(m2, s0, s1, phi)
# Calculate the transform fraction
fraction = transform(bbox[4:], phi, xcs.zeta())
# Ensure the minimum and maximum are 0 < 1
fmax = flex.max(fraction)
fmin = flex.min(fraction)
assert fmax <= (1.0 + 5e-15) and fmax > 0.0, f"{fmax:.16f} not between 0 and 1"
assert fmin >= 0.0 and fmin <= 1.0
# Ensure the fraction for each image frame adds up to 1.0 for
# all those frames completely within the grid
for j in range(1, fraction.all()[0] - 1):
tot = flex.sum(fraction[j : j + 1, :])
assert abs(tot - 1.0) < 1e-7
# Ensure the frames follow a progression through the grid. I.e,
# check that values increase then decrease and don't jump around
for j in range(fraction.all()[0]):
f = fraction[j : j + 1, :]
last = f[0]
rev = False
for i in range(1, len(f)):
curr = f[1]
if rev is False:
if curr < last:
rev = True
else:
assert curr <= last
last = curr
def test_map_frames_reverse(dials_data):
sequence = load.imageset(
dials_data("centroid_test_data", pathlib=True) / "sweep.json"
)
# Get the models
beam = sequence.get_beam()
detector = sequence.get_detector()
gonio = sequence.get_goniometer()
scan = sequence.get_scan()
# Set the delta_divergence/mosaicity
n_sigma = 3
sigma_divergence = 0.060 * math.pi / 180
mosaicity = 0.154 * math.pi / 180
delta_divergence = n_sigma * sigma_divergence
delta_mosaicity = n_sigma * mosaicity
# Set the grid size
grid_size = (4, 4, 4)
# Create the E3 fraction object
transform = MapFramesReverse(
scan.get_array_range()[0],
scan.get_oscillation(deg=False)[0],
scan.get_oscillation(deg=False)[1],
mosaicity,
n_sigma,
grid_size[2],
)
# Create the bounding box calculator
calculate_bbox = BBoxCalculator3D(
beam, detector, gonio, scan, delta_divergence, delta_mosaicity
)
s0 = beam.get_s0()
m2 = gonio.get_rotation_axis()
s0_length = matrix.col(beam.get_s0()).length()
for i in range(100):
# Get random x, y, z
x = random.uniform(0, 2000)
y = random.uniform(0, 2000)
z = random.uniform(0, 9)
# Get random s1, phi, panel
s1 = matrix.col(detector[0].get_pixel_lab_coord((x, y))).normalize() * s0_length
phi = scan.get_angle_from_array_index(z, deg=False)
panel = 0
# Calculate the bounding box
bbox = calculate_bbox(s1, phi, panel)
x1, x2 = bbox[0], bbox[1]
y1, y2 = bbox[2], bbox[3]
z1, z2 = bbox[4], bbox[5]
if x1 == 0 or y1 == 0 or z1 == 0:
continue
if x2 == 2000 or y2 == 2000 or z2 == 9:
continue
# Create the XDS coordinate system
xcs = CoordinateSystem(m2, s0, s1, phi)
# Calculate the transform fraction
fraction = transform(bbox[4:], phi, xcs.zeta())
# Ensure the minimum and maximum are 0 < 1
fmax = flex.max(fraction)
fmin = flex.min(fraction)
assert fmax <= 1.0 and fmax > 0.0
assert fmin >= 0.0 and fmin <= 1.0
# Ensure the fraction for image adds up to 1.0 for
# all those images completely within the image
for v3 in range(fraction.all()[0]):
tot = flex.sum(fraction[v3 : v3 + 1, :])
assert abs(tot - 1.0) < 1e-7
# Ensure the frames follow a progression through the grid. I.e,
# check that values increase then decrease and don't jump around
for v3 in range(fraction.all()[0]):
f = fraction[v3 : v3 + 1, :]
last = f[0]
rev = False
for i in range(1, len(f)):
curr = f[1]
if rev is False:
if curr < last:
rev = True
else:
assert curr <= last
last = curr
def test_map_forward_reverse(dials_data):
sequence = load.imageset(
dials_data("centroid_test_data", pathlib=True) / "sweep.json"
)
# Get the models
beam = sequence.get_beam()
detector = sequence.get_detector()
gonio = sequence.get_goniometer()
scan = sequence.get_scan()
# Set the delta_divergence/mosaicity
n_sigma = 3
sigma_divergence = 0.060 * math.pi / 180
mosaicity = 0.154 * math.pi / 180
delta_divergence = n_sigma * sigma_divergence
delta_mosaicity = n_sigma * mosaicity
# Set the grid size
grid_size = (4, 4, 4)
# Create the E3 fraction object
transform_forward = MapFramesForward(
scan.get_array_range()[0],
scan.get_oscillation(deg=False)[0],
scan.get_oscillation(deg=False)[1],
mosaicity,
n_sigma,
grid_size[2],
)
# Create the E3 fraction object
transform_reverse = MapFramesReverse(
scan.get_array_range()[0],
scan.get_oscillation(deg=False)[0],
scan.get_oscillation(deg=False)[1],
mosaicity,
n_sigma,
grid_size[2],
)
# Create the bounding box calculator
calculate_bbox = BBoxCalculator3D(
beam, detector, gonio, scan, delta_divergence, delta_mosaicity
)
s0 = beam.get_s0()
m2 = gonio.get_rotation_axis()
s0_length = matrix.col(beam.get_s0()).length()
for i in range(100):
# Get random x, y, z
x = random.uniform(0, 2000)
y = random.uniform(0, 2000)
z = random.uniform(0, 9)
# Get random s1, phi, panel
s1 = matrix.col(detector[0].get_pixel_lab_coord((x, y))).normalize() * s0_length
phi = scan.get_angle_from_array_index(z, deg=False)
panel = 0
# Calculate the bounding box
bbox = calculate_bbox(s1, phi, panel)
# Create the XDS coordinate system
xcs = CoordinateSystem(m2, s0, s1, phi)
# Calculate the transform fraction
forward_fraction = transform_forward(bbox[4:], phi, xcs.zeta())
# Calculate the transform fraction
reverse_fraction = transform_reverse(bbox[4:], phi, xcs.zeta())
# Check the same points are non-zero
eps = 1e-7
for j in range(forward_fraction.all()[0]):
for i in range(forward_fraction.all()[1]):
if forward_fraction[j, i] > 0.0:
assert reverse_fraction[i, j] > 0.0
else:
assert reverse_fraction[i, j] < eps
|
