?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
[haskell] use the right version of the hkl library >= 5.0.0.3381-1~
[hkl.git] / binoculars / hkl-binoculars.c
blob9cbfa86d322aa9e763850c8710ae1de8a3246064
1 /* This file is part of the hkl library.
2 *
3 * The hkl library is free software: you can redistribute it and/or modify
4 * it under the terms of the GNU General Public License as published by
5 * the Free Software Foundation, either version 3 of the License, or
6 * (at your option) any later version.
7 *
8 * The hkl library is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * GNU General Public License for more details.
12 *
13 * You should have received a copy of the GNU General Public License
14 * along with the hkl library. If not, see <http://www.gnu.org/licenses/>.
15 *
16 * Copyright (C) 2003-2023 Synchrotron SOLEIL
17 * L'Orme des Merisiers Saint-Aubin
18 * BP 48 91192 GIF-sur-YVETTE CEDEX
19 *
20 * Authors: Picca Frédéric-Emmanuel <picca@synchrotron-soleil.fr>
21 */
22
23 /* #define DEBUG */
24
25 #ifdef DEBUG
26 # define CGLM_DEFINE_PRINTS
27 # define CGLM_PRINT_PRECISION 7
28 #endif
29
30 #include <stdio.h>
31 #include <time.h>
32
33 #include "ccan/array_size/array_size.h"
34 #include "hkl-binoculars-private.h"
35 #include "hkl-matrix-private.h"
36 #include "hkl-quaternion-private.h"
37 #include "hkl-sample-private.h"
38 #include "hkl-vector-private.h"
39
40 /* mark the masked pixels with this value */
41 #define MASKED PTRDIFF_MAX
42 #define REMOVED 0
43
44 /* Math */
45
46 static inline ptrdiff_t min(ptrdiff_t x, ptrdiff_t y)
47 {
48 return y ^ ((x ^ y) & -(x < y));
49 }
50
51 static inline ptrdiff_t max(ptrdiff_t x, ptrdiff_t y)
52 {
53 return x ^ ((x ^ y) & -(x < y));
54 }
55
56 /* Axis Limits */
57
58 datatype(
59 AxisLimitType,
60 (NoLimit),
61 (Limit, ptrdiff_t)
62 );
63
64 struct _HklBinocularsAxisLimits
65 {
66 AxisLimitType imin;
67 AxisLimitType imax;
68 };
69
70 void hkl_binoculars_axis_limits_free(HklBinocularsAxisLimits *self)
71 {
72 free(self);
73 }
74
75 HklBinocularsAxisLimits *hkl_binoculars_axis_limits_new(const ptrdiff_t *imin,
76 const ptrdiff_t *imax)
77 {
78 HklBinocularsAxisLimits *self = g_new(HklBinocularsAxisLimits, 1);
79
80 if (NULL == imin){
81 self->imin = NoLimit();
82 } else {
83 self->imin = Limit(*imin);
84 }
85
86 if (NULL == imax){
87 self->imax = NoLimit();
88 } else {
89 self->imax = Limit(*imax);
90 }
91
92 return self;
93 }
94
95 /* Axis */
96
97 static inline double axis_min(const HklBinocularsAxis *self)
98 {
99 return self->imin * self->resolution;
100 }
101
102 static inline double axis_max(const HklBinocularsAxis *self)
103 {
104 return self->imax * self->resolution;
105 }
106
107 static inline void hkl_binoculars_axis_init(HklBinocularsAxis *self,
108 const char *name,
109 size_t index,
110 ptrdiff_t imin,
111 ptrdiff_t imax,
112 double resolution)
113 {
114 self->name = name;
115 self->index = index;
116 self->resolution = resolution;
117 self->imin = imin;
118 self->imax = imax;
119 }
120
121 static inline int hkl_binoculars_axis_cmp(const HklBinocularsAxis *self,
122 const HklBinocularsAxis *other)
123 {
124 int res = 0;
125
126 res |= strcmp(self->name, other->name);
127 res |= self->index != other->index;
128 res |= self->resolution != other->resolution;
129 res |= self->imin != other->imin;
130 res |= self->imax != other->imax;
131
132 return res;
133 }
134
135 double *hkl_binoculars_axis_array(const HklBinocularsAxis *self)
136 {
137 double *arr = g_new0(double, 6);
138
139 arr[0] = self->index;
140 arr[1] = axis_min(self);
141 arr[2] = axis_max(self);
142 arr[3] = self->resolution;
143 arr[4] = self->imin;
144 arr[5] = self->imax;
145
146 return arr;
147 }
148
149 /* check if *self contains *other */
150 static inline int hkl_binoculars_axis_contains_axis(const HklBinocularsAxis *self,
151 const HklBinocularsAxis *other)
152 {
153 return self->imin <= other->imin && self->imax >= other->imax;
154 }
155
156 static inline void hkl_binoculars_axis_merge(HklBinocularsAxis *self, const HklBinocularsAxis *other)
157 {
158 self->imin = min(self->imin, other->imin);
159 self->imax = max(self->imax, other->imax);
160 }
161
162 void hkl_binoculars_axis_fprintf(FILE *f, const HklBinocularsAxis *self)
163 {
164 fprintf(f, "%s : %ld min: %f(%ld) max: %f(%ld) res: %f size: %ld",
165 self->name, self->index,
166 axis_min(self), self->imin,
167 axis_max(self), self->imax,
168 self->resolution, axis_size(self));
169 }
170
171
172 /* darray_axis */
173
174 static inline void merge_axes(darray_axis *axes,
175 const darray_axis *others)
176 {
177 size_t i;
178
179 for(i=0; i<darray_size(*axes); ++i)
180 hkl_binoculars_axis_merge(&darray_item(*axes, i),
181 &darray_item(*others, i));
182 }
183
184 static inline int does_not_include(const darray_axis *axes,
185 const darray_axis *others)
186 {
187 size_t i;
188 int res = 0;
189
190 if (darray_size(*axes) == darray_size(*others)){
191 for(i=0; i<darray_size(*axes); ++i){
192 if (0 == hkl_binoculars_axis_contains_axis(&darray_item(*axes, i),
193 &darray_item(*others, i))){
194 res = 1;
195 break;
196 }
197 }
198 } else {
199 res = 1;
200 }
201
202 return res;
203 }
204
205
206 /* Space */
207
208 static inline void hkl_binoculars_space_item_fprintf(FILE *f, const HklBinocularsSpaceItem *self)
209 {
210 fprintf(f, "item->indexes(%p) v: %ld %ld %ld, intensity: %d", &self->indexes_0[0],
211 self->indexes_0[0], self->indexes_0[1], self->indexes_0[2], self->intensity);
212 }
213
214 static inline int space_is_empty(const HklBinocularsSpace *space)
215 {
216 return 0 == darray_size(space->items);
217 }
218
219 static inline void space_update_axes(HklBinocularsSpace *space,
220 const char *names[],
221 size_t n_pixels,
222 const double resolutions[])
223 {
224 size_t i;
225 HklBinocularsSpaceItem *item;
226 HklBinocularsSpaceItem minimum;
227 HklBinocularsSpaceItem maximum;
228
229 if (space_is_empty(space))
230 return;
231
232 minimum = maximum = darray_item(space->items, 0);
233
234 for(i=0; i<darray_size(space->axes); ++i){
235 darray_foreach(item, space->items){
236 minimum.indexes_0[i] = min(minimum.indexes_0[i], item->indexes_0[i]);
237 maximum.indexes_0[i] = max(maximum.indexes_0[i], item->indexes_0[i]);
238 }
239 }
240
241 for(i=0; i<darray_size(space->axes); ++i){
242 HklBinocularsAxis *axis = &darray_item(space->axes, i);
243 hkl_binoculars_axis_init(axis, names[i], i,
244 minimum.indexes_0[i], maximum.indexes_0[i],
245 resolutions[i]);
246 }
247 }
248
249 static inline int item_in_the_limits(const HklBinocularsSpaceItem *item,
250 const HklBinocularsAxisLimits **limits,
251 size_t n_limits)
252 {
253 int res = TRUE;
254
255 if (NULL != limits){
256 for(size_t i=0; i<n_limits; ++i){
257 ptrdiff_t v = item->indexes_0[i];
258
259 match(limits[i]->imin){
260 of(NoLimit){
261 }
262 of(Limit, imin){
263 if (v < *imin) res = FALSE;
264 break;
265 }
266 }
267
268 match(limits[i]->imax){
269 of(NoLimit){
270 }
271 of(Limit, imax){
272 if(v > *imax) res = FALSE;
273 break;
274 }
275 }
276 }
277 }
278
279 return res;
280 }
281
282 HklBinocularsSpace *hkl_binoculars_space_new(size_t max_items, size_t n_axes)
283 {
284 HklBinocularsSpace *self = g_new(HklBinocularsSpace, 1);
285
286 self->max_items = max_items;
287 darray_init(self->items);
288 darray_resize(self->items, max_items);
289 darray_init(self->axes);
290 darray_resize(self->axes, n_axes);
291
292 return self;
293 }
294
295 void hkl_binoculars_space_free(HklBinocularsSpace *self)
296 {
297 darray_free(self->axes);
298 darray_free(self->items);
299 free(self);
300 }
301
302 void hkl_binoculars_space_fprintf(FILE *f, const HklBinocularsSpace *self)
303 {
304 size_t masked;
305 HklBinocularsAxis *axis;
306
307 fprintf(f, "\nHklBinocularsSpace: %p", self);
308 fprintf(f, "\nn_indexes_0: %ld", darray_size(self->items));
309 if(space_is_empty(self)){
310 fprintf(f, "\nempty");
311 } else {
312 fprintf(f, "\nn_axes: %ld", darray_size(self->axes));
313 darray_foreach(axis, self->axes){
314 fprintf(f, "\n");
315 hkl_binoculars_axis_fprintf(f, axis);
316 }
317
318 masked = self->max_items - darray_size(self->items);
319 fprintf(f, "\nmasked pixels: %ld (%f%%)", masked, (double)masked / self->max_items * 100);
320 }
321 }
322
323 /* angles */
324
325 #define HKL_BINOCULARS_SPACE_ANGLES_IMPL(image_t) \
326 HKL_BINOCULARS_SPACE_ANGLES_DECL(image_t) \
327 { \
328 size_t i, j; \
329 const char * names[] = {"delta_lab", "gamma_lab", "tth"}; \
330 double delta0, gamma0, tth; \
331 \
332 assert(ARRAY_SIZE(names) == darray_size(space->axes)); \
333 assert(ARRAY_SIZE(names) == n_resolutions); \
334 assert(n_pixels == space->max_items); \
335 \
336 darray_size(space->items) = 0; \
337 \
338 const double *p_x = &pixels_coordinates[0 * n_pixels]; \
339 const double *p_y = &pixels_coordinates[1 * n_pixels]; \
340 const double *p_z = &pixels_coordinates[2 * n_pixels]; \
341 \
342 HklDetector *detector = hkl_detector_factory_new(HKL_DETECTOR_TYPE_0D); \
343 const HklQuaternion q = hkl_geometry_detector_rotation_get(geometry, detector); \
344 \
345 for(i=0;i<n_pixels;++i){ \
346 if(NULL == masked || 0 == masked[i]){ \
347 HklBinocularsSpaceItem item; \
348 HklVector v = {{p_x[i], p_y[i], p_z[i]}}; \
349 \
350 hkl_vector_rotated_quaternion(&v, &q); \
351 delta0 = atan2(v.data[2], v.data[0]); \
352 gamma0 = M_PI_2 - atan2(sqrt(v.data[2] * v.data[2] + v.data[0] * v.data[0]), v.data[1]); \
353 tth = acos(v.data[0]); \
354 \
355 v.data[0] = delta0 / M_PI * 180.0; \
356 v.data[1] = gamma0 / M_PI * 180.0; \
357 v.data[2] = tth / M_PI * 180.0; \
358 \
359 for(j=0; j<ARRAY_SIZE(names); ++j){ \
360 item.indexes_0[j] = rint(v.data[j] / resolutions[j]); \
361 } \
362 item.intensity = rint((double)image[i] * weight); \
363 \
364 \
365 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
366 darray_append(space->items, item); \
367 } \
368 } \
369 \
370 space_update_axes(space, names, n_pixels, resolutions); \
371 \
372 hkl_detector_free(detector); \
373 }
374
375 HKL_BINOCULARS_SPACE_ANGLES_IMPL(int32_t);
376 HKL_BINOCULARS_SPACE_ANGLES_IMPL(uint16_t);
377 HKL_BINOCULARS_SPACE_ANGLES_IMPL(uint32_t);
378
379 /* qcustom */
380
381 static const char **axis_name_from_subprojection(HklBinocularsQCustomSubProjectionEnum subprojection,
382 HklBinocularsSpace *space,
383 int n_resolutions)
384 {
385 const char **names;
386
387 switch(subprojection){
388 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QX_QY_QZ:
389 case HKL_BINOCULARS_QCUSTOM_NUM_SUBPROJECTIONS:
390 {
391 static const char *names_qx_qy_qz[] = {"qx", "qy", "qz"};
392 assert(ARRAY_SIZE(names_qx_qy_qz) == darray_size(space->axes));
393 assert(ARRAY_SIZE(names_qx_qy_qz) == n_resolutions);
394 names = names_qx_qy_qz;
395 break;
396 }
397 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_TTH_TIMESTAMP:
398 {
399 static const char *names_q_tth_timestamp[] = {"q", "tth", "timestamp"};
400 assert(ARRAY_SIZE(names_q_tth_timestamp) == darray_size(space->axes));
401 assert(ARRAY_SIZE(names_q_tth_timestamp) == n_resolutions);
402 names = names_q_tth_timestamp;
403 break;
404 }
405 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_TIMESTAMP:
406 {
407 static const char *names_q_timestamp[] = {"q", "timestamp"};
408 assert(ARRAY_SIZE(names_q_timestamp) <= darray_size(space->axes));
409 assert(ARRAY_SIZE(names_q_timestamp) <= n_resolutions);
410 names = names_q_timestamp;
411 break;
412 }
413 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QPAR_QPER_TIMESTAMP:
414 {
415 static const char *names_qpar_qper_timestamp[] = {"qpar", "qper", "timestamp"};
416 assert(ARRAY_SIZE(names_qpar_qper_timestamp) == darray_size(space->axes));
417 assert(ARRAY_SIZE(names_qpar_qper_timestamp) == n_resolutions);
418 names = names_qpar_qper_timestamp;
419 break;
420 }
421 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QPAR_QPER:
422 {
423 static const char *names_qpar_qper[] = {"qpar", "qper"};
424 assert(ARRAY_SIZE(names_qpar_qper) <= darray_size(space->axes));
425 assert(ARRAY_SIZE(names_qpar_qper) <= n_resolutions);
426 names = names_qpar_qper;
427 break;
428 }
429 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_PHI_QX:
430 {
431 static const char *names_q_phi_qx[] = {"q", "phi", "qx"};
432 assert(ARRAY_SIZE(names_q_phi_qx) == darray_size(space->axes));
433 assert(ARRAY_SIZE(names_q_phi_qx) == n_resolutions);
434 names = names_q_phi_qx;
435 break;
436 }
437 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_PHI_QY:
438 {
439 static const char *names_q_phi_qy[] = {"q", "phi", "qy"};
440 assert(ARRAY_SIZE(names_q_phi_qy) == darray_size(space->axes));
441 assert(ARRAY_SIZE(names_q_phi_qy) == n_resolutions);
442 names = names_q_phi_qy;
443 break;
444 }
445 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_PHI_QZ:
446 {
447 static const char *names_q_phi_qz[] = {"q", "phi", "qz"};
448 assert(ARRAY_SIZE(names_q_phi_qz) == darray_size(space->axes));
449 assert(ARRAY_SIZE(names_q_phi_qz) == n_resolutions);
450 names = names_q_phi_qz;
451 break;
452 }
453 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_STEREO:
454 {
455 static const char *names_q_stereo[] = {"q", "xp", "yp"};
456 assert(ARRAY_SIZE(names_q_stereo) == darray_size(space->axes));
457 assert(ARRAY_SIZE(names_q_stereo) == n_resolutions);
458 names = names_q_stereo;
459 break;
460 }
461 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_DELTALAB_GAMMALAB_SAMPLEAXIS:
462 {
463 static const char *names_deltalab_gammalab_sampleaxis[] = {"deltalab", "gammalab", "sampleaxis"};
464 assert(ARRAY_SIZE(names_deltalab_gammalab_sampleaxis) == darray_size(space->axes));
465 assert(ARRAY_SIZE(names_deltalab_gammalab_sampleaxis) == n_resolutions);
466 names = names_deltalab_gammalab_sampleaxis;
467 break;
468 }
469 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_X_Y_Z:
470 {
471 static const char *names_x_y_z[] = {"x", "y", "z"};
472 assert(ARRAY_SIZE(names_x_y_z) == darray_size(space->axes));
473 assert(ARRAY_SIZE(names_x_y_z) == n_resolutions);
474 names = names_x_y_z;
475 break;
476 }
477 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Y_Z_TIMESTAMP:
478 {
479 static const char *names_y_z_timestamp[] = {"y", "z", "timestamp"};
480 assert(ARRAY_SIZE(names_y_z_timestamp) == darray_size(space->axes));
481 assert(ARRAY_SIZE(names_y_z_timestamp) == n_resolutions);
482 names = names_y_z_timestamp;
483 break;
484 }
485 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_QPAR_QPER:
486 {
487 static const char *names_q_qpar_qper[] = {"q", "qpar", "qper"};
488 assert(ARRAY_SIZE(names_q_qpar_qper) == darray_size(space->axes));
489 assert(ARRAY_SIZE(names_q_qpar_qper) == n_resolutions);
490 names = names_q_qpar_qper;
491 break;
492 }
493 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QPARS_QPER_TIMESTAMP:
494 {
495 static const char *names_qpars_qper_timestamp[] = {"qpars", "qper", "timestamp"};
496 assert(ARRAY_SIZE(names_qpars_qper_timestamp) == darray_size(space->axes));
497 assert(ARRAY_SIZE(names_qpars_qper_timestamp) == n_resolutions);
498 names = names_qpars_qper_timestamp;
499 break;
500 }
501 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QPAR_QPER_SAMPLEAXIS:
502 {
503 static const char *names_qpar_qper_sampleaxis[] = {"qpar", "qper", "sample axis"};
504 assert(ARRAY_SIZE(names_qpar_qper_sampleaxis) == darray_size(space->axes));
505 assert(ARRAY_SIZE(names_qpar_qper_sampleaxis) == n_resolutions);
506 names = names_qpar_qper_sampleaxis;
507 break;
508 }
509 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_SAMPLEAXIS_TTH:
510 {
511 static const char *names_q_sampleaxis_tth[] = {"q", "sample axis", "tth"};
512 assert(ARRAY_SIZE(names_q_sampleaxis_tth) == darray_size(space->axes));
513 assert(ARRAY_SIZE(names_q_sampleaxis_tth) == n_resolutions);
514 names = names_q_sampleaxis_tth;
515 break;
516 }
517 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_SAMPLEAXIS_TIMESTAMP:
518 {
519 static const char *names_q_sampleaxis_timestamp[] = {"q", "sample axis", "timestamp"};
520 assert(ARRAY_SIZE(names_q_sampleaxis_timestamp) == darray_size(space->axes));
521 assert(ARRAY_SIZE(names_q_sampleaxis_timestamp) == n_resolutions);
522 names = names_q_sampleaxis_timestamp;
523 break;
524 }
525 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QX_QY_TIMESTAMP:
526 {
527 static const char *names_[] = {"qx", "qy", "timestamp"};
528 assert(ARRAY_SIZE(names_) == darray_size(space->axes));
529 assert(ARRAY_SIZE(names_) == n_resolutions);
530 names = names_;
531 break;
532 }
533 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QX_QZ_TIMESTAMP:
534 {
535 static const char *names_[] = {"qx", "qz", "timestamp"};
536 assert(ARRAY_SIZE(names_) == darray_size(space->axes));
537 assert(ARRAY_SIZE(names_) == n_resolutions);
538 names = names_;
539 break;
540 }
541 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QY_QZ_TIMESTAMP:
542 {
543 static const char *names_[] = {"qy", "qz", "timestamp"};
544 assert(ARRAY_SIZE(names_) == darray_size(space->axes));
545 assert(ARRAY_SIZE(names_) == n_resolutions);
546 names = names_;
547 break;
548 }
549 default:
550 {
551 static const char *names_qx_qy_qz[] = {"qx", "qy", "qz"};
552 assert(ARRAY_SIZE(names_qx_qy_qz) == darray_size(space->axes));
553 assert(ARRAY_SIZE(names_qx_qy_qz) == n_resolutions);
554 names = names_qx_qy_qz;
555 break;
556 }
557 }
558 return names;
559 }
560
561 static inline int not_masked(const uint8_t *masked, size_t idx)
562 {
563 return NULL == masked || 0 == masked[idx];
564 }
565
566 static inline double polarisation(vec3s kf, double weight, int do_polarisation)
567 {
568 if (do_polarisation){
569 CGLM_ALIGN_MAT vec3s epsilon = {{0, 1, 0}};
570 float p = glms_vec3_dot(epsilon, kf) / glms_vec3_norm2(kf);
571 weight = weight / (1 - p*p);
572 }
573 return weight;
574 }
575
576 #define HKL_BINOCULARS_SPACE_QCUSTOM_IMPL(image_t) \
577 HKL_BINOCULARS_SPACE_QCUSTOM_DECL(image_t) \
578 { \
579 size_t i; \
580 HklBinocularsSpaceItem item; \
581 double correction; \
582 const char **names = axis_name_from_subprojection(subprojection, space, n_resolutions); \
583 assert(n_pixels == space->max_items); \
584 \
585 const double *q_x = &pixels_coordinates[0 * n_pixels]; \
586 const double *q_y = &pixels_coordinates[1 * n_pixels]; \
587 const double *q_z = &pixels_coordinates[2 * n_pixels]; \
588 \
589 HklSample *sample = hkl_sample_new("test"); \
590 HklDetector *detector = hkl_detector_factory_new(HKL_DETECTOR_TYPE_0D); \
591 HklHolder *holder_d = hkl_geometry_detector_holder_get(geometry, detector); \
592 CGLM_ALIGN_MAT mat4s m_holder_d = hkl_binoculars_holder_transformation_get(holder_d); \
593 const HklVector ki_v = hkl_geometry_ki_get(geometry); \
594 CGLM_ALIGN_MAT vec3s ki = {{ki_v.data[0], ki_v.data[1], ki_v.data[2]}}; \
595 float k = glms_vec3_norm(ki); \
596 HklHolder *holder_s = hkl_geometry_sample_holder_get(geometry,sample); \
597 CGLM_ALIGN_MAT mat4s m_holder_s = hkl_binoculars_holder_transformation_get(holder_s); \
598 \
599 switch(surf){ \
600 case HKL_BINOCULARS_SURFACE_ORIENTATION_VERTICAL: \
601 { \
602 CGLM_ALIGN_MAT vec3s axis = GLMS_XUP; \
603 CGLM_ALIGN_MAT mat4s m_q_ub = glms_rotate_make(-M_PI_2, axis); \
604 \
605 m_holder_s = glms_mat4_mul(m_holder_s, m_q_ub); \
606 break; \
607 } \
608 case HKL_BINOCULARS_SURFACE_ORIENTATION_HORIZONTAL: \
609 case HKL_BINOCULARS_SURFACE_ORIENTATION_NUM_ORIENTATION: \
610 break; \
611 } \
612 \
613 CGLM_ALIGN_MAT vec3s euler_xyz = {{uqx, uqy, uqz}}; \
614 m_holder_s = glms_mat4_mul(m_holder_s, glms_euler_xyz(euler_xyz)); \
615 \
616 m_holder_s = glms_mat4_inv(m_holder_s); \
617 \
618 darray_size(space->items) = 0; \
619 \
620 glms_mat4_print(m_holder_s, stdout); \
621 glms_mat4_print(m_holder_d, stdout); \
622 \
623 switch(subprojection){ \
624 case HKL_BINOCULARS_QCUSTOM_NUM_SUBPROJECTIONS: \
625 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QX_QY_QZ: \
626 { \
627 for(i=0;i<n_pixels;++i){ \
628 if(not_masked(masked, i)){ \
629 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
630 \
631 v = glms_mat4_mulv3(m_holder_d, v, 1); \
632 v = glms_vec3_scale_as(v, k); \
633 correction = polarisation(v, weight, do_polarisation_correction); \
634 v = glms_vec3_sub(v, ki); \
635 v = glms_mat4_mulv3(m_holder_s, v, 0); \
636 \
637 item.indexes_0[0] = rint(v.raw[0] / resolutions[0]); \
638 item.indexes_0[1] = rint(v.raw[1] / resolutions[1]); \
639 item.indexes_0[2] = rint(v.raw[2] / resolutions[2]); \
640 item.intensity = rint((double)image[i] * correction); \
641 \
642 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
643 darray_append(space->items, item); \
644 } \
645 } \
646 break; \
647 } \
648 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_TTH_TIMESTAMP: \
649 { \
650 for(i=0;i<n_pixels;++i){ \
651 if(not_masked(masked, i)){ \
652 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
653 \
654 v = glms_mat4_mulv3(m_holder_d, v, 1); \
655 v = glms_vec3_scale_as(v, k); \
656 correction = polarisation(v, weight, do_polarisation_correction); \
657 v = glms_vec3_sub(v , ki); \
658 v = glms_mat4_mulv3(m_holder_s, v, 0); \
659 \
660 item.indexes_0[0] = rint(glms_vec3_norm(v) / resolutions[0]); \
661 item.indexes_0[1] = rint(asin(glms_vec3_norm(v) / 2 / k) * 2 / M_PI * 180 / resolutions[1]); \
662 item.indexes_0[2] = rint(timestamp / resolutions[2]); \
663 item.intensity = rint((double)image[i] * correction); \
664 \
665 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
666 darray_append(space->items, item); \
667 } \
668 } \
669 break; \
670 } \
671 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_TIMESTAMP: \
672 { \
673 for(i=0;i<n_pixels;++i){ \
674 if(not_masked(masked, i)){ \
675 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
676 \
677 v = glms_mat4_mulv3(m_holder_d, v, 1); \
678 v = glms_vec3_scale_as(v, k); \
679 correction = polarisation(v, weight, do_polarisation_correction); \
680 v = glms_vec3_sub(v , ki); \
681 v = glms_mat4_mulv3(m_holder_s, v, 0); \
682 \
683 item.indexes_0[0] = rint(glms_vec3_norm(v) / resolutions[0]); \
684 item.indexes_0[1] = rint(timestamp / resolutions[1]); \
685 item.indexes_0[2] = REMOVED; \
686 item.intensity = rint((double)image[i] * correction); \
687 \
688 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
689 darray_append(space->items, item); \
690 } \
691 } \
692 break; \
693 } \
694 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QPAR_QPER_TIMESTAMP: \
695 { \
696 for(i=0;i<n_pixels;++i){ \
697 if(not_masked(masked, i)){ \
698 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
699 \
700 v = glms_mat4_mulv3(m_holder_d, v, 1); \
701 v = glms_vec3_scale_as(v, k); \
702 correction = polarisation(v, weight, do_polarisation_correction); \
703 v = glms_vec3_sub(v , ki); \
704 v = glms_mat4_mulv3(m_holder_s, v, 0); \
705 \
706 item.indexes_0[0] = rint(sqrt(v.raw[0] * v.raw[0] + v.raw[1] * v.raw[1]) / resolutions[0]); \
707 item.indexes_0[1] = rint(v.raw[2] / resolutions[1]); \
708 item.indexes_0[2] = rint(timestamp / resolutions[2]); \
709 item.intensity = rint((double)image[i] * correction); \
710 \
711 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
712 darray_append(space->items, item); \
713 } \
714 } \
715 break; \
716 } \
717 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QPAR_QPER: \
718 { \
719 for(i=0;i<n_pixels;++i){ \
720 if(not_masked(masked, i)){ \
721 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
722 \
723 v = glms_mat4_mulv3(m_holder_d, v, 1); \
724 v = glms_vec3_scale_as(v, k); \
725 correction = polarisation(v, weight, do_polarisation_correction); \
726 v = glms_vec3_sub(v , ki); \
727 v = glms_mat4_mulv3(m_holder_s, v, 0); \
728 \
729 item.indexes_0[0] = rint(sqrt(v.raw[0] * v.raw[0] + v.raw[1] * v.raw[1]) / resolutions[0]); \
730 item.indexes_0[1] = rint(v.raw[2] / resolutions[1]); \
731 item.indexes_0[2] = REMOVED; \
732 item.intensity = rint((double)image[i] * correction); \
733 \
734 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
735 darray_append(space->items, item); \
736 } \
737 } \
738 break; \
739 } \
740 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_PHI_QX: \
741 { \
742 for(i=0;i<n_pixels;++i){ \
743 if(not_masked(masked, i)){ \
744 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
745 \
746 v = glms_mat4_mulv3(m_holder_d, v, 1); \
747 v = glms_vec3_scale_as(v, k); \
748 correction = polarisation(v, weight, do_polarisation_correction); \
749 v = glms_vec3_sub(v , ki); \
750 v = glms_mat4_mulv3(m_holder_s, v, 0); \
751 \
752 item.indexes_0[0] = rint(glms_vec3_norm(v) / resolutions[0]); \
753 item.indexes_0[1] = rint((atan2(v.raw[2], -v.raw[1])) / M_PI * 180 / resolutions[1]); \
754 item.indexes_0[2] = rint(v.raw[0] / resolutions[2]); \
755 item.intensity = rint((double)image[i] * correction); \
756 \
757 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
758 darray_append(space->items, item); \
759 } \
760 } \
761 break; \
762 } \
763 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_PHI_QY: \
764 { \
765 for(i=0;i<n_pixels;++i){ \
766 if(not_masked(masked, i)){ \
767 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
768 \
769 v = glms_mat4_mulv3(m_holder_d, v, 1); \
770 v = glms_vec3_scale_as(v, k); \
771 correction = polarisation(v, weight, do_polarisation_correction); \
772 v = glms_vec3_sub(v , ki); \
773 v = glms_mat4_mulv3(m_holder_s, v, 0); \
774 \
775 item.indexes_0[0] = rint(glms_vec3_norm(v) / resolutions[0]); \
776 item.indexes_0[1] = rint((atan2(v.raw[2], v.raw[0])) / M_PI * 180 / resolutions[1]); \
777 item.indexes_0[2] = rint(v.raw[1] / resolutions[2]); \
778 item.intensity = rint((double)image[i] * correction); \
779 \
780 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
781 darray_append(space->items, item); \
782 } \
783 } \
784 break; \
785 } \
786 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_PHI_QZ: \
787 { \
788 for(i=0;i<n_pixels;++i){ \
789 if(not_masked(masked, i)){ \
790 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
791 \
792 v = glms_mat4_mulv3(m_holder_d, v, 1); \
793 v = glms_vec3_scale_as(v, k); \
794 correction = polarisation(v, weight, do_polarisation_correction); \
795 v = glms_vec3_sub(v , ki); \
796 v = glms_mat4_mulv3(m_holder_s, v, 0); \
797 \
798 item.indexes_0[0] = rint(glms_vec3_norm(v) / resolutions[0]); \
799 item.indexes_0[1] = rint((atan2(v.raw[0], v.raw[1])) / M_PI * 180 / resolutions[1]); \
800 item.indexes_0[2] = rint(v.raw[2] / resolutions[2]); \
801 item.intensity = rint((double)image[i] * correction); \
802 \
803 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
804 darray_append(space->items, item); \
805 } \
806 } \
807 break; \
808 } \
809 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_STEREO: \
810 { \
811 for(i=0;i<n_pixels;++i){ \
812 if(not_masked(masked, i)){ \
813 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
814 \
815 v = glms_mat4_mulv3(m_holder_d, v, 1); \
816 v = glms_vec3_scale_as(v, k); \
817 correction = polarisation(v, weight, do_polarisation_correction); \
818 v = glms_vec3_sub(v , ki); \
819 v = glms_mat4_mulv3(m_holder_s, v, 0); \
820 \
821 item.indexes_0[0] = rint(glms_vec3_norm(v) / resolutions[0]); \
822 double ratio = v.raw[2] + item.indexes_0[0]; \
823 item.indexes_0[1] = rint(v.raw[0] / ratio / resolutions[1]); \
824 item.indexes_0[2] = rint(v.raw[1] / ratio / resolutions[2]); \
825 item.intensity = rint((double)image[i] * correction); \
826 \
827 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
828 darray_append(space->items, item); \
829 } \
830 } \
831 break; \
832 } \
833 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_DELTALAB_GAMMALAB_SAMPLEAXIS: \
834 { \
835 const HklParameter *p = hkl_geometry_axis_get(geometry, sample_axis, NULL); \
836 if (NULL != p){ \
837 int axis = rint(hkl_parameter_value_get(p, HKL_UNIT_USER) / resolutions[2]); \
838 \
839 for(i=0;i<n_pixels;++i){ \
840 if(not_masked(masked, i)){ \
841 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
842 \
843 v = glms_mat4_mulv3(m_holder_d, v, 1); \
844 correction = polarisation(v, weight, do_polarisation_correction); \
845 \
846 item.indexes_0[0] = rint(atan2(v.raw[2], sqrt(v.raw[0] * v.raw[0] + v.raw[1] * v.raw[1])) / M_PI * 180 / resolutions[0]); \
847 item.indexes_0[1] = rint(atan2(v.raw[1], v.raw[0]) / M_PI * 180 / resolutions[1]); \
848 item.indexes_0[2] = axis; \
849 item.intensity = rint((double)image[i] * correction); \
850 \
851 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
852 darray_append(space->items, item); \
853 } \
854 } \
855 } \
856 break; \
857 } \
858 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_X_Y_Z: \
859 { \
860 for(i=0;i<n_pixels;++i){ \
861 if(not_masked(masked, i)){ \
862 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
863 \
864 glms_vec3_print(v, stdout); \
865 v = glms_mat4_mulv3(m_holder_d, v, 1); \
866 correction = polarisation(v, weight, do_polarisation_correction); \
867 glms_mat4_print(m_holder_d, stdout); \
868 glms_vec3_print(v, stdout); \
869 \
870 item.indexes_0[0] = rint(v.raw[0] / resolutions[0]); \
871 item.indexes_0[1] = rint(v.raw[1] / resolutions[1]); \
872 item.indexes_0[2] = rint(v.raw[2] / resolutions[2]); \
873 item.intensity = rint((double)image[i] * correction); \
874 \
875 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
876 darray_append(space->items, item); \
877 } \
878 } \
879 break; \
880 } \
881 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Y_Z_TIMESTAMP: \
882 { \
883 for(i=0;i<n_pixels;++i){ \
884 if(not_masked(masked, i)){ \
885 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
886 \
887 v = glms_mat4_mulv3(m_holder_d, v, 1); \
888 correction = polarisation(v, weight, do_polarisation_correction); \
889 \
890 item.indexes_0[0] = rint(v.raw[1] / resolutions[0]); \
891 item.indexes_0[1] = rint(v.raw[2] / resolutions[1]); \
892 item.indexes_0[2] = rint(timestamp / resolutions[2]); \
893 item.intensity = rint((double)image[i] * correction); \
894 \
895 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
896 darray_append(space->items, item); \
897 } \
898 } \
899 break; \
900 } \
901 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_QPAR_QPER: \
902 { \
903 for(i=0;i<n_pixels;++i){ \
904 if(not_masked(masked, i)){ \
905 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
906 \
907 v = glms_mat4_mulv3(m_holder_d, v, 1); \
908 v = glms_vec3_scale_as(v, k); \
909 correction = polarisation(v, weight, do_polarisation_correction); \
910 v = glms_vec3_sub(v , ki); \
911 v = glms_mat4_mulv3(m_holder_s, v, 0); \
912 \
913 item.indexes_0[0] = rint(glms_vec3_norm(v) / resolutions[0]); \
914 item.indexes_0[1] = rint(sqrt(v.raw[0] * v.raw[0] + v.raw[1] * v.raw[1]) / resolutions[1]); \
915 item.indexes_0[2] = rint(v.raw[2] / resolutions[2]); \
916 item.intensity = rint((double)image[i] * correction); \
917 \
918 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
919 darray_append(space->items, item); \
920 } \
921 } \
922 break; \
923 } \
924 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QPARS_QPER_TIMESTAMP: \
925 { \
926 for(i=0;i<n_pixels;++i){ \
927 if(not_masked(masked, i)){ \
928 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
929 \
930 v = glms_mat4_mulv3(m_holder_d, v, 1); \
931 v = glms_vec3_scale_as(v, k); \
932 correction = polarisation(v, weight, do_polarisation_correction); \
933 v = glms_vec3_sub(v , ki); \
934 v = glms_mat4_mulv3(m_holder_s, v, 0); \
935 \
936 item.indexes_0[0] = rint(sqrt(v.raw[0] * v.raw[0] + v.raw[1] * v.raw[1]) / resolutions[0]); \
937 item.indexes_0[1] = rint(v.raw[2] / resolutions[1]); \
938 item.indexes_0[2] = rint(timestamp / resolutions[2]); \
939 item.intensity = rint((double)image[i] * correction); \
940 \
941 if(v.raw[1] != 0.0){ \
942 if(signbit(v.raw[1]) != 0){ \
943 item.indexes_0[0] = -item.indexes_0[0]; \
944 } \
945 }else{ \
946 if(signbit(v.raw[0]) != 0){ \
947 item.indexes_0[0] = -item.indexes_0[0]; \
948 } \
949 } \
950 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
951 darray_append(space->items, item); \
952 } \
953 } \
954 break; \
955 } \
956 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QPAR_QPER_SAMPLEAXIS: \
957 { \
958 const HklParameter *p = hkl_geometry_axis_get(geometry, sample_axis, NULL); \
959 if (NULL != p){ \
960 int axis = rint(hkl_parameter_value_get(p, HKL_UNIT_USER) / resolutions[2]); \
961 \
962 for(i=0;i<n_pixels;++i){ \
963 if(not_masked(masked, i)){ \
964 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
965 \
966 v = glms_mat4_mulv3(m_holder_d, v, 1); \
967 v = glms_vec3_scale_as(v, k); \
968 correction = polarisation(v, weight, do_polarisation_correction); \
969 v = glms_vec3_sub(v , ki); \
970 v = glms_mat4_mulv3(m_holder_s, v, 0); \
971 \
972 item.indexes_0[0] = rint(sqrt(v.raw[0] * v.raw[0] + v.raw[1] * v.raw[1]) / resolutions[0]); \
973 item.indexes_0[1] = rint(v.raw[2] / resolutions[1]); \
974 item.indexes_0[2] = axis; \
975 item.intensity = rint((double)image[i] * correction); \
976 \
977 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
978 darray_append(space->items, item); \
979 } \
980 } \
981 break; \
982 } \
983 } \
984 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_SAMPLEAXIS_TTH: \
985 { \
986 const HklParameter *p = hkl_geometry_axis_get(geometry, sample_axis, NULL); \
987 if (NULL != p){ \
988 int axis = rint(hkl_parameter_value_get(p, HKL_UNIT_USER) / resolutions[1]); \
989 \
990 for(i=0;i<n_pixels;++i){ \
991 if(not_masked(masked, i)){ \
992 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
993 \
994 v = glms_mat4_mulv3(m_holder_d, v, 1); \
995 v = glms_vec3_scale_as(v, k); \
996 correction = polarisation(v, weight, do_polarisation_correction); \
997 v = glms_vec3_sub(v , ki); \
998 v = glms_mat4_mulv3(m_holder_s, v, 0); \
999 \
1000 item.indexes_0[0] = rint(glms_vec3_norm(v) / resolutions[0]); \
1001 item.indexes_0[1] = axis; \
1002 item.indexes_0[2] = rint(asin(glms_vec3_norm(v) / 2 / k) * 2 / M_PI * 180 / resolutions[2]); \
1003 item.intensity = rint((double)image[i] * correction); \
1004 \
1005 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
1006 darray_append(space->items, item); \
1007 } \
1008 } \
1009 break; \
1010 } \
1011 } \
1012 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_Q_SAMPLEAXIS_TIMESTAMP: \
1013 { \
1014 const HklParameter *p = hkl_geometry_axis_get(geometry, sample_axis, NULL); \
1015 if (NULL != p){ \
1016 int axis = rint(hkl_parameter_value_get(p, HKL_UNIT_USER) / resolutions[1]); \
1017 \
1018 for(i=0;i<n_pixels;++i){ \
1019 if(not_masked(masked, i)){ \
1020 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
1021 \
1022 v = glms_mat4_mulv3(m_holder_d, v, 1); \
1023 v = glms_vec3_scale_as(v, k); \
1024 correction = polarisation(v, weight, do_polarisation_correction); \
1025 v = glms_vec3_sub(v , ki); \
1026 v = glms_mat4_mulv3(m_holder_s, v, 0); \
1027 \
1028 item.indexes_0[0] = rint(glms_vec3_norm(v) / resolutions[0]); \
1029 item.indexes_0[1] = axis; \
1030 item.indexes_0[2] = rint(timestamp / resolutions[2]); \
1031 item.intensity = rint((double)image[i] * correction); \
1032 \
1033 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
1034 darray_append(space->items, item); \
1035 } \
1036 } \
1037 break; \
1038 } \
1039 } \
1040 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QX_QY_TIMESTAMP: \
1041 { \
1042 for(i=0;i<n_pixels;++i){ \
1043 if(not_masked(masked, i)){ \
1044 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
1045 \
1046 v = glms_mat4_mulv3(m_holder_d, v, 1); \
1047 v = glms_vec3_scale_as(v, k); \
1048 correction = polarisation(v, weight, do_polarisation_correction); \
1049 v = glms_vec3_sub(v, ki); \
1050 v = glms_mat4_mulv3(m_holder_s, v, 0); \
1051 \
1052 item.indexes_0[0] = rint(v.raw[0] / resolutions[0]); \
1053 item.indexes_0[1] = rint(v.raw[1] / resolutions[1]); \
1054 item.indexes_0[2] = rint(timestamp / resolutions[2]); \
1055 item.intensity = rint((double)image[i] * correction); \
1056 \
1057 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
1058 darray_append(space->items, item); \
1059 } \
1060 } \
1061 break; \
1062 } \
1063 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QX_QZ_TIMESTAMP: \
1064 { \
1065 for(i=0;i<n_pixels;++i){ \
1066 if(not_masked(masked, i)){ \
1067 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
1068 \
1069 v = glms_mat4_mulv3(m_holder_d, v, 1); \
1070 v = glms_vec3_scale_as(v, k); \
1071 correction = polarisation(v, weight, do_polarisation_correction); \
1072 v = glms_vec3_sub(v, ki); \
1073 v = glms_mat4_mulv3(m_holder_s, v, 0); \
1074 \
1075 item.indexes_0[0] = rint(v.raw[0] / resolutions[0]); \
1076 item.indexes_0[1] = rint(v.raw[2] / resolutions[1]); \
1077 item.indexes_0[2] = rint(timestamp / resolutions[2]); \
1078 item.intensity = rint((double)image[i] * correction); \
1079 \
1080 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
1081 darray_append(space->items, item); \
1082 } \
1083 } \
1084 break; \
1085 } \
1086 case HKL_BINOCULARS_QCUSTOM_SUB_PROJECTION_QY_QZ_TIMESTAMP: \
1087 { \
1088 for(i=0;i<n_pixels;++i){ \
1089 if(not_masked(masked, i)){ \
1090 CGLM_ALIGN_MAT vec3s v = {{q_x[i], q_y[i], q_z[i]}}; \
1091 \
1092 v = glms_mat4_mulv3(m_holder_d, v, 1); \
1093 v = glms_vec3_scale_as(v, k); \
1094 correction = polarisation(v, weight, do_polarisation_correction); \
1095 v = glms_vec3_sub(v, ki); \
1096 v = glms_mat4_mulv3(m_holder_s, v, 0); \
1097 \
1098 item.indexes_0[0] = rint(v.raw[1] / resolutions[0]); \
1099 item.indexes_0[1] = rint(v.raw[2] / resolutions[1]); \
1100 item.indexes_0[2] = rint(timestamp / resolutions[2]); \
1101 item.intensity = rint((double)image[i] * correction); \
1102 \
1103 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
1104 darray_append(space->items, item); \
1105 } \
1106 } \
1107 break; \
1108 } \
1109 } \
1110 space_update_axes(space, names, n_pixels, resolutions); \
1111 hkl_detector_free(detector); \
1112 hkl_sample_free(sample); \
1113 }
1114
1115 HKL_BINOCULARS_SPACE_QCUSTOM_IMPL(int32_t);
1116 HKL_BINOCULARS_SPACE_QCUSTOM_IMPL(uint16_t);
1117 HKL_BINOCULARS_SPACE_QCUSTOM_IMPL(uint32_t);
1118
1119 /* hkl */
1120
1121 #define HKL_BINOCULARS_SPACE_HKL_IMPL(image_t) \
1122 HKL_BINOCULARS_SPACE_HKL_DECL(image_t) \
1123 { \
1124 size_t i; \
1125 double correction; \
1126 const char * names[] = {"H", "K", "L"}; \
1127 HklBinocularsSpaceItem item; \
1128 \
1129 assert(ARRAY_SIZE(names) == darray_size(space->axes)); \
1130 assert(ARRAY_SIZE(names) == n_resolutions); \
1131 assert(n_pixels == space->max_items); \
1132 \
1133 const double *h = &pixels_coordinates[0 * n_pixels]; \
1134 const double *k = &pixels_coordinates[1 * n_pixels]; \
1135 const double *l = &pixels_coordinates[2 * n_pixels]; \
1136 \
1137 HklDetector *detector = hkl_detector_factory_new(HKL_DETECTOR_TYPE_0D); \
1138 HklHolder *holder_d = hkl_geometry_detector_holder_get(geometry, detector); \
1139 CGLM_ALIGN_MAT mat4s m_holder_d = hkl_binoculars_holder_transformation_get(holder_d); \
1140 const HklVector ki_v = hkl_geometry_ki_get(geometry); \
1141 CGLM_ALIGN_MAT vec3s ki = {{ki_v.data[0], ki_v.data[1], ki_v.data[2]}}; \
1142 float K = glms_vec3_norm(ki); \
1143 HklHolder *holder_s = hkl_geometry_sample_holder_get(geometry,sample); \
1144 CGLM_ALIGN_MAT mat4s m_holder_s = hkl_binoculars_holder_transformation_get(holder_s); \
1145 \
1146 const HklMatrix *UB = hkl_sample_UB_get(sample); \
1147 CGLM_ALIGN_MAT mat4s ub = {{{UB->data[0][0], UB->data[1][0], UB->data[2][0], 0}, \
1148 {UB->data[0][1], UB->data[1][1], UB->data[2][1], 0}, \
1149 {UB->data[0][2], UB->data[1][2], UB->data[2][2], 0}, \
1150 {0, 0, 0, 1}}}; \
1151 glms_mat4_print(ub, stdout); \
1152 m_holder_s = glms_mat4_mul(m_holder_s, ub); \
1153 glms_mat4_print(m_holder_s, stdout); \
1154 m_holder_s = glms_mat4_inv(m_holder_s); \
1155 \
1156 darray_size(space->items) = 0; \
1157 \
1158 glms_mat4_print(m_holder_s, stdout); \
1159 glms_mat4_print(m_holder_d, stdout); \
1160 \
1161 for(i=0;i<n_pixels;++i){ \
1162 if(not_masked(masked, i)){ \
1163 CGLM_ALIGN_MAT vec3s v = {{h[i], k[i], l[i]}}; \
1164 \
1165 v = glms_mat4_mulv3(m_holder_d, v, 1); \
1166 v = glms_vec3_scale_as(v, K); \
1167 correction = polarisation(v, weight, do_polarisation_correction); \
1168 v = glms_vec3_sub(v, ki); \
1169 v = glms_mat4_mulv3(m_holder_s, v, 0); \
1170 \
1171 item.indexes_0[0] = rint(v.raw[0] / resolutions[0]); \
1172 item.indexes_0[1] = rint(v.raw[1] / resolutions[1]); \
1173 item.indexes_0[2] = rint(v.raw[2] / resolutions[2]); \
1174 item.intensity = rint((double)image[i] * correction); \
1175 \
1176 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
1177 darray_append(space->items, item); \
1178 } \
1179 } \
1180 \
1181 space_update_axes(space, names, n_pixels, resolutions); \
1182 \
1183 hkl_detector_free(detector); \
1184 }
1185
1186 HKL_BINOCULARS_SPACE_HKL_IMPL(int32_t);
1187 HKL_BINOCULARS_SPACE_HKL_IMPL(uint16_t);
1188 HKL_BINOCULARS_SPACE_HKL_IMPL(uint32_t);
1189
1190 /* test */
1191
1192 #define HKL_BINOCULARS_SPACE_TEST_IMPL(image_t) \
1193 HKL_BINOCULARS_SPACE_TEST_DECL(image_t) \
1194 { \
1195 size_t i; \
1196 const char * names[] = {"H", "K", "L"}; \
1197 double correction; \
1198 HklBinocularsSpaceItem item; \
1199 \
1200 assert(ARRAY_SIZE(names) == darray_size(space->axes)); \
1201 assert(ARRAY_SIZE(names) == n_resolutions); \
1202 assert(n_pixels == space->max_items); \
1203 \
1204 const double *h = &pixels_coordinates[0 * n_pixels]; \
1205 const double *k = &pixels_coordinates[1 * n_pixels]; \
1206 const double *l = &pixels_coordinates[2 * n_pixels]; \
1207 \
1208 HklDetector *detector = hkl_detector_factory_new(HKL_DETECTOR_TYPE_0D); \
1209 HklHolder *holder_d = hkl_geometry_detector_holder_get(geometry, detector); \
1210 CGLM_ALIGN_MAT mat4s m_holder_d = hkl_binoculars_holder_transformation_get(holder_d); \
1211 const HklVector ki_v = hkl_geometry_ki_get(geometry); \
1212 CGLM_ALIGN_MAT vec3s ki = {{ki_v.data[0], ki_v.data[1], ki_v.data[2]}}; \
1213 float K = glms_vec3_norm(ki); \
1214 HklHolder *holder_s = hkl_geometry_sample_holder_get(geometry,sample); \
1215 CGLM_ALIGN_MAT mat4s m_holder_s = hkl_binoculars_holder_transformation_get(holder_s); \
1216 \
1217 const HklMatrix *UB = hkl_sample_UB_get(sample); \
1218 CGLM_ALIGN_MAT mat4s ub = {{{UB->data[0][0], UB->data[1][0], UB->data[2][0], 0}, \
1219 {UB->data[0][1], UB->data[1][1], UB->data[2][1], 0}, \
1220 {UB->data[0][2], UB->data[1][2], UB->data[2][2], 0}, \
1221 {0, 0, 0, 1}}}; \
1222 glms_mat4_print(ub, stdout); \
1223 m_holder_s = glms_mat4_mul(m_holder_s, ub); \
1224 glms_mat4_print(m_holder_s, stdout); \
1225 m_holder_s = glms_mat4_inv(m_holder_s); \
1226 \
1227 darray_size(space->items) = 0; \
1228 \
1229 glms_mat4_print(m_holder_s, stdout); \
1230 glms_mat4_print(m_holder_d, stdout); \
1231 \
1232 for(i=0;i<n_pixels;++i){ \
1233 if(not_masked(masked, i)){ \
1234 CGLM_ALIGN_MAT vec3s v = {{h[i], k[i], l[i]}}; \
1235 \
1236 v = glms_mat4_mulv3(m_holder_d, v, 1); \
1237 v = glms_vec3_scale_as(v, K); \
1238 correction = polarisation(v, weight, do_polarisation_correction); \
1239 v = glms_vec3_sub(v, ki); \
1240 v = glms_mat4_mulv3(m_holder_s, v, 0); \
1241 \
1242 item.indexes_0[0] = rint(v.raw[0] / resolutions[0]); \
1243 item.indexes_0[1] = rint(v.raw[1] / resolutions[1]); \
1244 item.indexes_0[2] = rint(v.raw[2] / resolutions[2]); \
1245 item.intensity = rint((double)image[i] * correction); \
1246 \
1247 if(TRUE == item_in_the_limits(&item, limits, n_limits)) \
1248 darray_append(space->items, item); \
1249 } \
1250 } \
1251 \
1252 space_update_axes(space, names, n_pixels, resolutions); \
1253 \
1254 hkl_detector_free(detector); \
1255 }
1256
1257 HKL_BINOCULARS_SPACE_TEST_IMPL(int32_t);
1258 HKL_BINOCULARS_SPACE_TEST_IMPL(uint16_t);
1259 HKL_BINOCULARS_SPACE_TEST_IMPL(uint32_t);
1260
1261
1262 /* Cube */
1263
1264 /* this method compute the linear coordinates of the first element in
1265 * the absolute coordinates of the bin */
1266 static inline ptrdiff_t compute_offset0(const darray_axis *axes)
1267 {
1268 size_t i;
1269 ptrdiff_t len = 1;
1270 ptrdiff_t offset = 0;
1271 size_t n_axes = darray_size(*axes);
1272
1273 for(i=0; i<n_axes; ++i){
1274 const HklBinocularsAxis *axis = &darray_item(*axes, n_axes-1-i);
1275 offset += len * axis->imin;
1276 len *= axis_size(axis);
1277 }
1278 return offset;
1279 }
1280
1281 static inline HklBinocularsCube *empty_cube_from_axes(const darray_axis *axes)
1282 {
1283 HklBinocularsCube *self = NULL;
1284
1285 if (0 != darray_size(*axes)){
1286 HklBinocularsAxis *axis;
1287
1288 self = g_new0(HklBinocularsCube, 1);
1289 darray_foreach(axis, *axes){
1290 darray_append(self->axes, *axis);
1291 }
1292 self->offset0 = compute_offset0(&self->axes);
1293 self->photons = NULL;
1294 self->contributions = NULL;
1295 }
1296
1297 return self;
1298 }
1299
1300 static inline size_t cube_size(const HklBinocularsCube *self)
1301 {
1302 size_t n = 1;
1303 HklBinocularsAxis *axis;
1304
1305 darray_foreach(axis, self->axes){
1306 n *= axis_size(axis);
1307 }
1308
1309 return n;
1310 }
1311
1312 static inline size_t malloc_cube(HklBinocularsCube *self)
1313 {
1314 size_t n = cube_size(self);
1315
1316 self->photons = malloc(n * sizeof(*self->photons));
1317 self->contributions = malloc(n * sizeof(*self->contributions));
1318
1319 return n;
1320 }
1321
1322 static inline size_t calloc_cube(HklBinocularsCube *self)
1323 {
1324 size_t n = cube_size(self);
1325
1326 self->photons = calloc(n, sizeof(*self->photons));
1327 self->contributions = calloc(n, sizeof(*self->contributions));
1328
1329 return n;
1330 }
1331
1332 /* Using this method the Cube has already the right dimensions, we
1333 * just add the Space data into it. */
1334 static inline void add_non_empty_space(HklBinocularsCube *cube,
1335 const HklBinocularsSpace *space)
1336 {
1337 size_t i;
1338 size_t n_axes = darray_size(cube->axes);
1339 ptrdiff_t lens[n_axes];
1340 HklBinocularsSpaceItem *item;
1341
1342 assert(n_axes == darray_size(space->axes));
1343
1344 /* compute the lens */
1345 lens[0] = 1;
1346
1347 for(i=1; i<n_axes; ++i){
1348 lens[i] = lens[i - 1] * axis_size(&darray_item(cube->axes, n_axes - i));
1349 }
1350
1351 darray_foreach(item, space->items){
1352 ptrdiff_t w = -cube->offset0;
1353
1354 for(i=0; i<n_axes; ++i){
1355 w += lens[i] * item->indexes_0[n_axes - 1 - i];
1356 /* fprintf(stdout, " %ld %ld", lens[i], item->indexes_0[n_axes - 1 - i]); */
1357 }
1358
1359 /* fprintf(stdout, " w: %ld %ld\n", w, cube_size(cube)); */
1360 cube->photons[w] += item->intensity;
1361 cube->contributions[w] += 1;
1362 }
1363 }
1364
1365 static inline int cube_is_empty(const HklBinocularsCube *self)
1366 {
1367 return 0 == darray_size(self->axes);
1368 }
1369
1370 HklBinocularsCube *hkl_binoculars_cube_new_empty(void)
1371 {
1372 HklBinocularsCube *self = g_new(HklBinocularsCube, 1);
1373
1374 darray_init(self->axes);
1375 self->offset0 = 0;
1376 self->photons = NULL;
1377 self->contributions = NULL;
1378
1379 return self;
1380 }
1381
1382 void hkl_binoculars_cube_free(HklBinocularsCube *self)
1383 {
1384 free(self->contributions);
1385 free(self->photons);
1386 darray_free(self->axes);
1387 free(self);
1388 }
1389
1390 unsigned int hkl_binoculars_cube_cmp(const HklBinocularsCube *self, const HklBinocularsCube *other)
1391 {
1392 int i;
1393 unsigned int res = 0;
1394
1395 res |= darray_size(self->axes) != darray_size(other->axes);
1396 res |= self->offset0 != other->offset0;
1397 for(i=0; i<darray_size(self->axes); ++i){
1398 res |= hkl_binoculars_axis_cmp(&darray_item(self->axes, i),
1399 &darray_item(other->axes, i));
1400 }
1401
1402 /* for(i=0; i<cube_size(self); ++i){ */
1403 /* res |= self->photons[i] != other->photons[i]; */
1404 /* res |= self->contributions[i] != other->contributions[i]; */
1405 /* if(res){ */
1406 /* fprintf(stdout, "\nphotons: %d %u %u ", i, self->photons[i], other->photons[i]); */
1407 /* fprintf(stdout, "\ncontributions: %d %u %u", i, self->contributions[i], other->contributions[i]); */
1408 /* break; */
1409 /* } */
1410 /* } */
1411
1412 if(res){
1413 fprintf(stdout, "cube: self\n");
1414 hkl_binoculars_cube_fprintf(stdout, self);
1415 fprintf(stdout, "cube: other\n");
1416 hkl_binoculars_cube_fprintf(stdout, other);
1417 }
1418
1419 return res;
1420 }
1421
1422 void hkl_binoculars_cube_fprintf(FILE *f, const HklBinocularsCube *self)
1423 {
1424 HklBinocularsAxis *axis;
1425
1426 fprintf(f, "HklBinocularsCube: %p", self);
1427
1428 fprintf(f, "\nn_axes: %ld", darray_size(self->axes));
1429 darray_foreach(axis, self->axes){
1430 fprintf(f, "\n");
1431 hkl_binoculars_axis_fprintf(f, axis);
1432 }
1433 fprintf(f, "\nphotons: %p", self->photons);
1434 fprintf(f, "\ncontributions: %p", self->contributions);
1435 }
1436
1437 void hkl_binoculars_cube_dims(const HklBinocularsCube *self, size_t ndims, size_t *dims)
1438 {
1439 HklBinocularsAxis *axis;
1440
1441 assert(ndims == darray_size(self->axes));
1442
1443 darray_foreach(axis, self->axes){
1444 *(dims++) = axis_size(axis);
1445 }
1446 }
1447
1448 static inline size_t find_first_non_empty_space_index(size_t n_spaces,
1449 const HklBinocularsSpace *const *spaces)
1450 {
1451 size_t i;
1452
1453 for(i=0; i<n_spaces; ++i)
1454 if(!space_is_empty(spaces[i]))
1455 return i;
1456
1457 return n_spaces;
1458 }
1459
1460 HklBinocularsCube *hkl_binoculars_cube_new(size_t n_spaces,
1461 const HklBinocularsSpace *const *spaces)
1462 {
1463 size_t i, i0;
1464 HklBinocularsCube *self = NULL;
1465
1466 i0 = find_first_non_empty_space_index(n_spaces, spaces);
1467
1468 if(i0 < n_spaces){
1469 self = empty_cube_from_axes(&spaces[i0]->axes);
1470
1471 if(NULL != self){
1472 /* compute the final cube dimensions and the index offset */
1473 for(i=i0; i<n_spaces; ++i){
1474 const HklBinocularsSpace *space = spaces[i];
1475
1476 if(space_is_empty(space))
1477 continue;
1478
1479 merge_axes(&self->axes, &space->axes);
1480 }
1481 self->offset0 = compute_offset0(&self->axes);
1482
1483 /* allocated the final cube photons and contributions */
1484 calloc_cube(self);
1485
1486 /* add all the spaces */
1487 for(i=i0; i<n_spaces; ++i){
1488 const HklBinocularsSpace *space = spaces[i];
1489
1490 if(space_is_empty(space))
1491 continue;
1492
1493 add_non_empty_space(self, space);
1494 }
1495 }
1496 }
1497
1498 return self;
1499 }
1500
1501 HklBinocularsCube *hkl_binoculars_cube_new_empty_from_cube(const HklBinocularsCube *cube)
1502 {
1503 HklBinocularsCube *self = empty_cube_from_axes(&cube->axes);
1504 if(NULL != self)
1505 calloc_cube(self);
1506 else{
1507 self = hkl_binoculars_cube_new_empty();
1508 }
1509
1510
1511 return self;
1512 }
1513
1514 HklBinocularsCube *hkl_binoculars_cube_new_from_space(const HklBinocularsSpace *space)
1515 {
1516 HklBinocularsCube *self = NULL;
1517
1518 if(space_is_empty(space)){
1519 self = hkl_binoculars_cube_new_empty();
1520 } else {
1521 self = empty_cube_from_axes(&space->axes);
1522 if(NULL != self){
1523 calloc_cube(self);
1524
1525 add_non_empty_space(self, space);
1526 }
1527 }
1528
1529 return self;
1530 }
1531
1532 HklBinocularsCube *hkl_binoculars_cube_new_copy(const HklBinocularsCube *src)
1533 {
1534 size_t n;
1535 HklBinocularsCube *self = empty_cube_from_axes(&src->axes);
1536
1537 if(NULL != self){
1538 /* allocate the final cube */
1539 n = malloc_cube(self);
1540
1541 /* copy the data */
1542 if(self->photons)
1543 memcpy(self->photons, src->photons, n * sizeof(*self->photons));
1544 if(self->contributions)
1545 memcpy(self->contributions, src->contributions, n * sizeof(*self->contributions));
1546 }
1547
1548 return self;
1549 }
1550
1551 static inline void cube_add_cube_2(HklBinocularsCube *self,
1552 const HklBinocularsCube *other)
1553 {
1554 size_t i, j;
1555 size_t i_offset, j_offset;
1556
1557 size_t stride_i = 1;
1558 size_t stride_j = stride_i * axis_size(&darray_item(self->axes, 1));
1559
1560 /* fill the values of other */
1561 size_t stride_i_other = 1;
1562 size_t stride_j_other = stride_i_other * axis_size(&darray_item(other->axes, 1));
1563
1564 i_offset = darray_item(other->axes, 1).imin - darray_item(self->axes, 1).imin;
1565 j_offset = darray_item(other->axes, 0).imin - darray_item(self->axes, 0).imin;
1566
1567 for(j=0; j<axis_size(&darray_item(other->axes, 0)); ++j){
1568 for(i=0; i<axis_size(&darray_item(other->axes, 1)); ++i){
1569 size_t w = (i + i_offset) * stride_i + (j + j_offset) * stride_j;
1570 size_t w1 = i * stride_i_other + j * stride_j_other;
1571
1572 self->photons[w] += other->photons[w1];
1573 self->contributions[w] += other->contributions[w1];
1574 }
1575 }
1576 }
1577
1578 static inline void cube_add_cube_3(HklBinocularsCube *self,
1579 const HklBinocularsCube *other)
1580 {
1581 size_t i, j, k;
1582 size_t i_offset, j_offset, k_offset;
1583
1584 size_t stride_i = 1;
1585 size_t stride_j = stride_i * axis_size(&darray_item(self->axes, 2));
1586 size_t stride_k = stride_j * axis_size(&darray_item(self->axes, 1));
1587
1588 /* fill the values of other */
1589 size_t stride_i_other = 1;
1590 size_t stride_j_other = stride_i_other * axis_size(&darray_item(other->axes, 2));
1591 size_t stride_k_other = stride_j_other * axis_size(&darray_item(other->axes, 1));
1592
1593 i_offset = darray_item(other->axes, 2).imin - darray_item(self->axes, 2).imin;
1594 j_offset = darray_item(other->axes, 1).imin - darray_item(self->axes, 1).imin;
1595 k_offset = darray_item(other->axes, 0).imin - darray_item(self->axes, 0).imin;
1596
1597 for(k=0; k<axis_size(&darray_item(other->axes, 0)); ++k){
1598 for(j=0; j<axis_size(&darray_item(other->axes, 1)); ++j){
1599 for(i=0; i<axis_size(&darray_item(other->axes, 2)); ++i){
1600 size_t w = (i + i_offset) * stride_i + (j + j_offset) * stride_j + (k + k_offset) * stride_k;
1601 size_t w1 = i * stride_i_other + j * stride_j_other + k * stride_k_other;
1602
1603 self->photons[w] += other->photons[w1];
1604 self->contributions[w] += other->contributions[w1];
1605 }
1606 }
1607 }
1608 }
1609
1610 static inline void cube_add_cube(HklBinocularsCube *self,
1611 const HklBinocularsCube *other)
1612 {
1613 assert(darray_size(self->axes) == darray_size(other->axes));
1614
1615 switch(darray_size(self->axes)){
1616 case 2: cube_add_cube_2(self, other);
1617 break;
1618 case 3: cube_add_cube_3(self, other);
1619 break;
1620 default: assert(0);
1621 }
1622 }
1623
1624 static inline void compute_strides(const darray_axis *axes, size_t strides[], size_t n_strides)
1625 {
1626 size_t i;
1627
1628 strides[0] = 1;
1629 for(i=1; i<n_strides; ++i)
1630 strides[i] = strides[i-1] * axis_size(&darray_item(*axes, n_strides - i + 1));
1631 }
1632
1633 HklBinocularsCube *hkl_binoculars_cube_new_merge(const HklBinocularsCube *cube1,
1634 const HklBinocularsCube *cube2)
1635 {
1636 HklBinocularsCube *self = empty_cube_from_axes(&cube1->axes);
1637 if(NULL != self){
1638 merge_axes(&self->axes, &cube2->axes);
1639 calloc_cube(self);
1640
1641 cube_add_cube(self, cube1);
1642 cube_add_cube(self, cube2);
1643 }
1644 return self;
1645 }
1646
1647 static inline void switch_content(HklBinocularsCube *self,
1648 HklBinocularsCube *other)
1649 {
1650 unsigned int *ptr;
1651 darray_axis tmp;
1652 ptrdiff_t offset0;
1653
1654 tmp = self->axes;
1655 self->axes = other->axes;
1656 other->axes = tmp;
1657
1658 offset0 = self->offset0;
1659 self->offset0 = other->offset0;
1660 other->offset0 = offset0;
1661
1662 ptr = self->photons;
1663 self->photons = other->photons;
1664 other->photons = ptr;
1665
1666 ptr = self->contributions;
1667 self->contributions = other->contributions;
1668 other->contributions = ptr;
1669 }
1670
1671 void hkl_binoculars_cube_add_space(HklBinocularsCube *self,
1672 const HklBinocularsSpace *space)
1673 {
1674 #ifdef DEBUG
1675 fprintf(stdout, "\nENTERING hkl_binoculars_cube_add_space:\n");
1676 hkl_binoculars_cube_fprintf(stdout, self);
1677 hkl_binoculars_space_fprintf(stdout, space);
1678 #endif
1679 /* check the compatibility of the cube and the space. */
1680 if (1 != space_is_empty(space)){
1681 if (does_not_include(&self->axes, &space->axes)){
1682 #ifdef DEBUG
1683 fprintf(stdout, "\nthe Cube does not contain the space, so create a new cube.");
1684 #endif
1685 if(0 != darray_size(self->axes)){ /* self cube is not empty */
1686 HklBinocularsCube *cube = empty_cube_from_axes(&self->axes);
1687 if(NULL != cube){
1688 merge_axes(&cube->axes, &space->axes); /* circonscript */
1689 cube->offset0 = compute_offset0(&cube->axes);
1690 calloc_cube(cube);
1691 cube_add_cube(cube, self);
1692 switch_content(self, cube);
1693 hkl_binoculars_cube_free(cube);
1694 }
1695 } else { /* self cube is empty */
1696 HklBinocularsCube *cube = empty_cube_from_axes(&space->axes);
1697 if(NULL != cube){
1698 cube->offset0 = compute_offset0(&cube->axes);
1699 calloc_cube(cube);
1700 switch_content(self, cube);
1701 hkl_binoculars_cube_free(cube);
1702 }
1703 }
1704 }
1705 add_non_empty_space(self, space);
1706 }
1707
1708 #ifdef DEBUG
1709 fprintf(stdout, "\n");
1710 hkl_binoculars_cube_fprintf(stdout, self);
1711 fprintf(stdout, "\nLEAVING hkl_binoculars_cube_add_space:\n");
1712 #endif
1713 }
Diffractometer computation library