?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
d2/align: Revert last patch until the corresponding code in Libale (currently
[Ale.git] / d2 / align.h
blob296df91e4a9c65cd6f36ab8eb8ce827b7ce80b40
1 // Copyright 2002, 2004, 2007 David Hilvert <dhilvert@auricle.dyndns.org>,
2 // <dhilvert@ugcs.caltech.edu>
3
4 /* This file is part of the Anti-Lamenessing Engine.
5
6 The Anti-Lamenessing Engine is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 3 of the License, or
9 (at your option) any later version.
10
11 The Anti-Lamenessing Engine is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with the Anti-Lamenessing Engine; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 */
20
21 /*
22 * align.h: Handle alignment of frames.
23 */
24
25 #ifndef __d2align_h__
26 #define __d2align_h__
27
28 #include "filter.h"
29 #include "transformation.h"
30 #include "image.h"
31 #include "point.h"
32 #include "render.h"
33 #include "tfile.h"
34 #include "image_rw.h"
35
36 class align {
37 private:
38
39 /*
40 * Alignment properties
41 */
42
43 static ale_align_properties align_properties() {
44 static ale_align_properties data = NULL;
45
46 if (data == NULL)
47 data = ale_new_align_properties();
48
49 assert(data);
50
51 return data;
52 }
53
54 /*
55 * Private data members
56 */
57
58 static ale_pos scale_factor;
59
60 /*
61 * Original frame transformation
62 */
63 static transformation orig_t;
64
65 /*
66 * Keep data older than latest
67 */
68 static int _keep;
69 static transformation *kept_t;
70 static int *kept_ok;
71
72 /*
73 * Transformation file handlers
74 */
75
76 static tload_t *tload;
77 static tsave_t *tsave;
78
79 /*
80 * Control point variables
81 */
82
83 static const point **cp_array;
84 static unsigned int cp_count;
85
86 /*
87 * Reference rendering to align against
88 */
89
90 static render *reference;
91 static filter::scaled_filter *interpolant;
92 static ale_image reference_image;
93
94 /*
95 * Per-pixel alignment weight map
96 */
97
98 static ale_image weight_map;
99
100 /*
101 * Frequency-dependent alignment weights
102 */
103
104 static double horiz_freq_cut;
105 static double vert_freq_cut;
106 static double avg_freq_cut;
107 static const char *fw_output;
108
109 /*
110 * Algorithmic alignment weighting
111 */
112
113 static const char *wmx_exec;
114 static const char *wmx_file;
115 static const char *wmx_defs;
116
117 /*
118 * Non-certainty alignment weights
119 */
120
121 static ale_image alignment_weights;
122
123 /*
124 * Latest transformation.
125 */
126
127 static transformation latest_t;
128
129 /*
130 * Flag indicating whether the latest transformation
131 * resulted in a match.
132 */
133
134 static int latest_ok;
135
136 /*
137 * Frame number most recently aligned.
138 */
139
140 static int latest;
141
142 /*
143 * Exposure registration
144 *
145 * 0. Preserve the original exposure of images.
146 *
147 * 1. Match exposure between images.
148 *
149 * 2. Use only image metadata for registering exposure.
150 */
151
152 static int _exp_register;
153
154 /*
155 * Alignment class.
156 *
157 * 0. Translation only. Only adjust the x and y position of images.
158 * Do not rotate input images or perform projective transformations.
159 *
160 * 1. Euclidean transformations only. Adjust the x and y position
161 * of images and the orientation of the image about the image center.
162 *
163 * 2. Perform general projective transformations. See the file gpt.h
164 * for more information about general projective transformations.
165 */
166
167 static int alignment_class;
168
169 /*
170 * Default initial alignment type.
171 *
172 * 0. Identity transformation.
173 *
174 * 1. Most recently accepted frame's final transformation.
175 */
176
177 static int default_initial_alignment_type;
178
179 /*
180 * Projective group behavior
181 *
182 * 0. Perturb in output coordinates.
183 *
184 * 1. Perturb in source coordinates
185 */
186
187 static int perturb_type;
188
189 /*
190 * Alignment for failed frames -- default or optimal?
191 *
192 * A frame that does not meet the match threshold can be assigned the
193 * best alignment found, or can be assigned its alignment default.
194 */
195
196 static int is_fail_default;
197
198 /*
199 * Alignment code.
200 *
201 * 0. Align images with an error contribution from each color channel.
202 *
203 * 1. Align images with an error contribution only from the green channel.
204 * Other color channels do not affect alignment.
205 *
206 * 2. Align images using a summation of channels. May be useful when dealing
207 * with images that have high frequency color ripples due to color aliasing.
208 */
209
210 static int channel_alignment_type;
211
212 /*
213 * Error metric exponent
214 */
215
216 static ale_real metric_exponent;
217
218 /*
219 * Match threshold
220 */
221
222 static float match_threshold;
223
224 /*
225 * Perturbation lower and upper bounds.
226 */
227
228 static ale_pos perturb_lower;
229 static int perturb_lower_percent;
230 static ale_pos perturb_upper;
231 static int perturb_upper_percent;
232
233 /*
234 * Preferred level-of-detail scale factor is 2^lod_preferred/perturb.
235 */
236
237 static int lod_preferred;
238
239 /*
240 * Minimum dimension for reduced LOD.
241 */
242
243 static int min_dimension;
244
245 /*
246 * Maximum rotational perturbation
247 */
248
249 static ale_pos rot_max;
250
251 /*
252 * Barrel distortion alignment multiplier
253 */
254
255 static ale_pos bda_mult;
256
257 /*
258 * Barrel distortion maximum adjustment rate
259 */
260
261 static ale_pos bda_rate;
262
263 /*
264 * Alignment match sum
265 */
266
267 static ale_accum match_sum;
268
269 /*
270 * Alignment match count.
271 */
272
273 static int match_count;
274
275 /*
276 * Monte Carlo parameter
277 */
278
279 static ale_pos _mc;
280
281 /*
282 * Certainty weight flag
283 *
284 * 0. Don't use certainty weights for alignment.
285 *
286 * 1. Use certainty weights for alignment.
287 */
288
289 static int certainty_weights;
290
291 /*
292 * Global search parameter
293 *
294 * 0. Local: Local search only.
295 * 1. Inner: Alignment reference image inner region
296 * 2. Outer: Alignment reference image outer region
297 * 3. All: Alignment reference image inner and outer regions.
298 * 4. Central: Inner if possible; else, best of inner and outer.
299 * 5. Points: Align by control points.
300 */
301
302 static int _gs;
303
304 /*
305 * Minimum overlap for global searches
306 */
307
308 static ale_accum _gs_mo;
309 static int gs_mo_percent;
310
311 /*
312 * Minimum certainty for multi-alignment element registration.
313 */
314
315 static ale_real _ma_cert;
316
317 /*
318 * Exclusion regions
319 */
320
321 static exclusion *ax_parameters;
322 static int ax_count;
323
324 /*
325 * Types for scale clusters.
326 */
327
328 struct nl_scale_cluster {
329 const image *accum_max;
330 const image *accum_min;
331 const image *certainty_max;
332 const image *certainty_min;
333 const image *aweight_max;
334 const image *aweight_min;
335 exclusion *ax_parameters;
336
337 ale_pos input_scale;
338 const image *input_certainty_max;
339 const image *input_certainty_min;
340 const image *input_max;
341 const image *input_min;
342 };
343
344 struct scale_cluster {
345 const image *accum;
346 const image *certainty;
347 const image *aweight;
348 exclusion *ax_parameters;
349
350 ale_pos input_scale;
351 const image *input_certainty;
352 const image *input;
353
354 nl_scale_cluster *nl_scale_clusters;
355 };
356
357 /*
358 * Check for exclusion region coverage in the reference
359 * array.
360 */
361 static int ref_excluded(int i, int j, point offset, exclusion *params, int param_count) {
362 for (int idx = 0; idx < param_count; idx++)
363 if (params[idx].type == exclusion::RENDER
364 && i + offset[0] >= params[idx].x[0]
365 && i + offset[0] <= params[idx].x[1]
366 && j + offset[1] >= params[idx].x[2]
367 && j + offset[1] <= params[idx].x[3])
368 return 1;
369
370 return 0;
371 }
372
373 /*
374 * Check for exclusion region coverage in the input
375 * array.
376 */
377 static int input_excluded(ale_pos ti, ale_pos tj, exclusion *params, int param_count) {
378 for (int idx = 0; idx < param_count; idx++)
379 if (params[idx].type == exclusion::FRAME
380 && ti >= params[idx].x[0]
381 && ti <= params[idx].x[1]
382 && tj >= params[idx].x[2]
383 && tj <= params[idx].x[3])
384 return 1;
385
386 return 0;
387 }
388
389 /*
390 * Overlap function. Determines the number of pixels in areas where
391 * the arrays overlap. Uses the reference array's notion of pixel
392 * positions.
393 */
394 static unsigned int overlap(struct scale_cluster c, transformation t, int ax_count) {
395 assert (reference_image);
396
397 unsigned int result = 0;
398
399 point offset = c.accum->offset();
400
401 for (unsigned int i = 0; i < c.accum->height(); i++)
402 for (unsigned int j = 0; j < c.accum->width(); j++) {
403
404 if (ref_excluded(i, j, offset, c.ax_parameters, ax_count))
405 continue;
406
407 /*
408 * Transform
409 */
410
411 struct point q;
412
413 q = (c.input_scale < 1.0 && interpolant == NULL)
414 ? t.scaled_inverse_transform(
415 point(i + offset[0], j + offset[1]))
416 : t.unscaled_inverse_transform(
417 point(i + offset[0], j + offset[1]));
418
419 ale_pos ti = q[0];
420 ale_pos tj = q[1];
421
422 /*
423 * Check that the transformed coordinates are within
424 * the boundaries of array c.input, and check that the
425 * weight value in the accumulated array is nonzero,
426 * unless we know it is nonzero by virtue of the fact
427 * that it falls within the region of the original
428 * frame (e.g. when we're not increasing image
429 * extents). Also check for frame exclusion.
430 */
431
432 if (input_excluded(ti, tj, c.ax_parameters, ax_count))
433 continue;
434
435 if (ti >= 0
436 && ti <= c.input->height() - 1
437 && tj >= 0
438 && tj <= c.input->width() - 1
439 && c.certainty->get_pixel(i, j)[0] != 0)
440 result++;
441 }
442
443 return result;
444 }
445
446 /*
447 * Monte carlo iteration class.
448 *
449 * Monte Carlo alignment has been used for statistical comparisons in
450 * spatial registration, and is now used for tonal registration
451 * and final match calculation.
452 */
453
454 /*
455 * We use a random process for which the expected number of sampled
456 * pixels is +/- .000003 from the coverage in the range [.005,.995] for
457 * an image with 100,000 pixels. (The actual number may still deviate
458 * from the expected number by more than this amount, however.) The
459 * method is as follows:
460 *
461 * We have coverage == USE/ALL, or (expected # pixels to use)/(# total
462 * pixels). We derive from this SKIP/USE.
463 *
464 * SKIP/USE == (SKIP/ALL)/(USE/ALL) == (1 - (USE/ALL))/(USE/ALL)
465 *
466 * Once we have SKIP/USE, we know the expected number of pixels to skip
467 * in each iteration. We use a random selection process that provides
468 * SKIP/USE close to this calculated value.
469 *
470 * If we can draw uniformly to select the number of pixels to skip, we
471 * do. In this case, the maximum number of pixels to skip is twice the
472 * expected number.
473 *
474 * If we cannot draw uniformly, we still assign equal probability to
475 * each of the integer values in the interval [0, 2 * (SKIP/USE)], but
476 * assign an unequal amount to the integer value ceil(2 * SKIP/USE) +
477 * 1.
478 */
479
480 /*
481 * When reseeding the random number generator, we want the same set of
482 * pixels to be used in cases where two alignment options are compared.
483 * If we wanted to avoid bias from repeatedly utilizing the same seed,
484 * we could seed with the number of the frame most recently aligned:
485 *
486 * srand(latest);
487 *
488 * However, in cursory tests, it seems okay to just use the default
489 * seed of 1, and so we do this, since it is simpler; both of these
490 * approaches to reseeding achieve better results than not reseeding.
491 * (1 is the default seed according to the GNU Manual Page for
492 * rand(3).)
493 *
494 * For subdomain calculations, we vary the seed by adding the subdomain
495 * index.
496 */
497
498 class mc_iterate {
499 ale_pos mc_max;
500 unsigned int index;
501 unsigned int index_max;
502 int i_min;
503 int i_max;
504 int j_min;
505 int j_max;
506
507 rng_t rng;
508
509 public:
510 mc_iterate(int _i_min, int _i_max, int _j_min, int _j_max, unsigned int subdomain)
511 : rng() {
512
513 ale_pos coverage;
514
515 i_min = _i_min;
516 i_max = _i_max;
517 j_min = _j_min;
518 j_max = _j_max;
519
520 if (i_max < i_min || j_max < j_min)
521 index_max = 0;
522 else
523 index_max = (i_max - i_min) * (j_max - j_min);
524
525 if (index_max < 500 || _mc > 100 || _mc <= 0)
526 coverage = 1;
527 else
528 coverage = _mc / 100;
529
530 double su = (1 - coverage) / coverage;
531
532 mc_max = (floor(2*su) * (1 + floor(2*su)) + 2*su)
533 / (2 + 2 * floor(2*su) - 2*su);
534
535 rng.seed(1 + subdomain);
536
537 #define FIXED16 3
538 #if ALE_COORDINATES == FIXED16
539 /*
540 * XXX: This calculation might not yield the correct
541 * expected value.
542 */
543 index = -1 + (int) ceil(((ale_pos) mc_max+1)
544 / (ale_pos) ( (1 + 0xffffff)
545 / (1 + (rng.get() & 0xffffff))));
546 #else
547 index = -1 + (int) ceil((ale_accum) (mc_max+1)
548 * ( (1 + ((ale_accum) (rng.get())) )
549 / (1 + ((ale_accum) RAND_MAX)) ));
550 #endif
551 #undef FIXED16
552 }
553
554 int get_i() {
555 return index / (j_max - j_min) + i_min;
556 }
557
558 int get_j() {
559 return index % (j_max - j_min) + j_min;
560 }
561
562 void operator++(int whats_this_for) {
563 #define FIXED16 3
564 #if ALE_COORDINATES == FIXED16
565 index += (int) ceil(((ale_pos) mc_max+1)
566 / (ale_pos) ( (1 + 0xffffff)
567 / (1 + (rng.get() & 0xffffff))));
568 #else
569 index += (int) ceil((ale_accum) (mc_max+1)
570 * ( (1 + ((ale_accum) (rng.get())) )
571 / (1 + ((ale_accum) RAND_MAX)) ));
572 #endif
573 #undef FIXED16
574 }
575
576 int done() {
577 return (index >= index_max);
578 }
579 };
580
581 /*
582 * Not-quite-symmetric difference function. Determines the difference in areas
583 * where the arrays overlap. Uses the reference array's notion of pixel positions.
584 *
585 * For the purposes of determining the difference, this function divides each
586 * pixel value by the corresponding image's average pixel magnitude, unless we
587 * are:
588 *
589 * a) Extending the boundaries of the image, or
590 *
591 * b) following the previous frame's transform
592 *
593 * If we are doing monte-carlo pixel sampling for alignment, we
594 * typically sample a subset of available pixels; otherwise, we sample
595 * all pixels.
596 *
597 */
598
599 template<class diff_trans>
600 class diff_stat_generic {
601
602 transformation::elem_bounds_t elem_bounds;
603
604 struct run {
605
606 diff_trans offset;
607
608 ale_accum result;
609 ale_accum divisor;
610
611 point max, min;
612 ale_accum centroid[2], centroid_divisor;
613 ale_accum de_centroid[2], de_centroid_v, de_sum;
614
615 void init() {
616
617 result = 0;
618 divisor = 0;
619
620 min = point::posinf();
621 max = point::neginf();
622
623 centroid[0] = 0;
624 centroid[1] = 0;
625 centroid_divisor = 0;
626
627 de_centroid[0] = 0;
628 de_centroid[1] = 0;
629
630 de_centroid_v = 0;
631
632 de_sum = 0;
633 }
634
635 void init(diff_trans _offset) {
636 offset = _offset;
637 init();
638 }
639
640 /*
641 * Required for STL sanity.
642 */
643 run() : offset(diff_trans::eu_identity()) {
644 init();
645 }
646
647 run(diff_trans _offset) : offset(_offset) {
648 init(_offset);
649 }
650
651 void add(const run &_run) {
652 result += _run.result;
653 divisor += _run.divisor;
654
655 for (int d = 0; d < 2; d++) {
656 if (min[d] > _run.min[d])
657 min[d] = _run.min[d];
658 if (max[d] < _run.max[d])
659 max[d] = _run.max[d];
660 centroid[d] += _run.centroid[d];
661 de_centroid[d] += _run.de_centroid[d];
662 }
663
664 centroid_divisor += _run.centroid_divisor;
665 de_centroid_v += _run.de_centroid_v;
666 de_sum += _run.de_sum;
667 }
668
669 run(const run &_run) : offset(_run.offset) {
670
671 /*
672 * Initialize
673 */
674 init(_run.offset);
675
676 /*
677 * Add
678 */
679 add(_run);
680 }
681
682 run &operator=(const run &_run) {
683
684 /*
685 * Initialize
686 */
687 init(_run.offset);
688
689 /*
690 * Add
691 */
692 add(_run);
693
694 return *this;
695 }
696
697 ~run() {
698 }
699
700 ale_accum get_error() const {
701 return pow(result / divisor, 1/(ale_accum) metric_exponent);
702 }
703
704 void sample(int f, scale_cluster c, int i, int j, point t, point u,
705 const run &comparison) {
706
707 pixel pa = c.accum->get_pixel(i, j);
708
709 ale_real this_result[2] = { 0, 0 };
710 ale_real this_divisor[2] = { 0, 0 };
711
712 pixel p[2];
713 pixel weight[2];
714 weight[0] = pixel(1, 1, 1);
715 weight[1] = pixel(1, 1, 1);
716
717 pixel tm = offset.get_tonal_multiplier(point(i, j) + c.accum->offset());
718
719 if (interpolant != NULL) {
720 interpolant->filtered(i, j, &p[0], &weight[0], 1, f);
721
722 if (weight[0].min_norm() > ale_real_weight_floor) {
723 p[0] /= weight[0];
724 } else {
725 return;
726 }
727
728 } else {
729 p[0] = c.input->get_bl(t);
730 }
731
732 p[0] *= tm;
733
734 if (u.defined()) {
735 p[1] = c.input->get_bl(u);
736 p[1] *= tm;
737 }
738
739
740 /*
741 * Handle certainty.
742 */
743
744 if (certainty_weights == 1) {
745
746 /*
747 * For speed, use arithmetic interpolation (get_bl(.))
748 * instead of geometric (get_bl(., 1))
749 */
750
751 weight[0] *= c.input_certainty->get_bl(t);
752 if (u.defined())
753 weight[1] *= c.input_certainty->get_bl(u);
754 weight[0] *= c.certainty->get_pixel(i, j);
755 weight[1] *= c.certainty->get_pixel(i, j);
756 }
757
758 if (c.aweight != NULL) {
759 weight[0] *= c.aweight->get_pixel(i, j);
760 weight[1] *= c.aweight->get_pixel(i, j);
761 }
762
763 /*
764 * Update sampling area statistics
765 */
766
767 if (min[0] > i)
768 min[0] = i;
769 if (min[1] > j)
770 min[1] = j;
771 if (max[0] < i)
772 max[0] = i;
773 if (max[1] < j)
774 max[1] = j;
775
776 centroid[0] += (weight[0][0] + weight[0][1] + weight[0][2]) * i;
777 centroid[1] += (weight[0][0] + weight[0][1] + weight[0][2]) * j;
778 centroid_divisor += (weight[0][0] + weight[0][1] + weight[0][2]);
779
780 /*
781 * Determine alignment type.
782 */
783
784 for (int m = 0; m < (u.defined() ? 2 : 1); m++)
785 if (channel_alignment_type == 0) {
786 /*
787 * Align based on all channels.
788 */
789
790
791 for (int k = 0; k < 3; k++) {
792 ale_real achan = pa[k];
793 ale_real bchan = p[m][k];
794
795 this_result[m] += weight[m][k] * pow(fabs(achan - bchan), metric_exponent);
796 this_divisor[m] += weight[m][k] * pow(achan > bchan ? achan : bchan, metric_exponent);
797 }
798 } else if (channel_alignment_type == 1) {
799 /*
800 * Align based on the green channel.
801 */
802
803 ale_real achan = pa[1];
804 ale_real bchan = p[m][1];
805
806 this_result[m] = weight[m][1] * pow(fabs(achan - bchan), metric_exponent);
807 this_divisor[m] = weight[m][1] * pow(achan > bchan ? achan : bchan, metric_exponent);
808 } else if (channel_alignment_type == 2) {
809 /*
810 * Align based on the sum of all channels.
811 */
812
813 ale_real asum = 0;
814 ale_real bsum = 0;
815 ale_real wsum = 0;
816
817 for (int k = 0; k < 3; k++) {
818 asum += pa[k];
819 bsum += p[m][k];
820 wsum += weight[m][k] / 3;
821 }
822
823 this_result[m] = wsum * pow(fabs(asum - bsum), metric_exponent);
824 this_divisor[m] = wsum * pow(asum > bsum ? asum : bsum, metric_exponent);
825 }
826
827 if (u.defined()) {
828 // ale_real de = fabs(this_result[0] / this_divisor[0]
829 // - this_result[1] / this_divisor[1]);
830 ale_real de = fabs(this_result[0] - this_result[1]);
831
832 de_centroid[0] += de * (ale_real) i;
833 de_centroid[1] += de * (ale_real) j;
834
835 de_centroid_v += de * (ale_real) t.lengthto(u);
836
837 de_sum += de;
838 }
839
840 result += (this_result[0]);
841 divisor += (this_divisor[0]);
842 }
843
844 void rescale(ale_pos scale) {
845 offset.rescale(scale);
846
847 de_centroid[0] *= scale;
848 de_centroid[1] *= scale;
849 de_centroid_v *= scale;
850 }
851
852 point get_centroid() {
853 point result = point(centroid[0] / centroid_divisor, centroid[1] / centroid_divisor);
854
855 assert (finite(centroid[0])
856 && finite(centroid[1])
857 && (result.defined() || centroid_divisor == 0));
858
859 return result;
860 }
861
862 point get_error_centroid() {
863 point result = point(de_centroid[0] / de_sum, de_centroid[1] / de_sum);
864 return result;
865 }
866
867
868 ale_pos get_error_perturb() {
869 ale_pos result = de_centroid_v / de_sum;
870
871 return result;
872 }
873
874 };
875
876 /*
877 * When non-empty, runs.front() is best, runs.back() is
878 * testing.
879 */
880
881 std::vector<run> runs;
882
883 /*
884 * old_runs stores the latest available perturbation set for
885 * each multi-alignment element.
886 */
887
888 typedef int run_index;
889 std::map<run_index, run> old_runs;
890
891 static void *diff_subdomain(void *args);
892
893 struct subdomain_args {
894 struct scale_cluster c;
895 std::vector<run> runs;
896 int ax_count;
897 int f;
898 int i_min, i_max, j_min, j_max;
899 int subdomain;
900 };
901
902 struct scale_cluster si;
903 int ax_count;
904 int frame;
905
906 std::vector<ale_pos> perturb_multipliers;
907
908 public:
909 void diff(struct scale_cluster c, const diff_trans &t,
910 int _ax_count, int f) {
911
912 if (runs.size() == 2)
913 runs.pop_back();
914
915 runs.push_back(run(t));
916
917 si = c;
918 ax_count = _ax_count;
919 frame = f;
920
921 ui::get()->d2_align_sample_start();
922
923 if (interpolant != NULL) {
924
925 /*
926 * XXX: This has not been tested, and may be completely
927 * wrong.
928 */
929
930 transformation tt = transformation::eu_identity();
931 tt.set_current_element(t);
932 interpolant->set_parameters(tt, c.input, c.accum->offset());
933 }
934
935 int N;
936 #ifdef USE_PTHREAD
937 N = thread::count();
938
939 pthread_t *threads = (pthread_t *) malloc(sizeof(pthread_t) * N);
940 pthread_attr_t *thread_attr = (pthread_attr_t *) malloc(sizeof(pthread_attr_t) * N);
941
942 #else
943 N = 1;
944 #endif
945
946 subdomain_args *args = new subdomain_args[N];
947
948 transformation::elem_bounds_int_t b = elem_bounds.scale_to_bounds(c.accum->height(), c.accum->width());
949
950 // fprintf(stdout, "[%d %d] [%d %d]\n",
951 // global_i_min, global_i_max, global_j_min, global_j_max);
952
953 for (int ti = 0; ti < N; ti++) {
954 args[ti].c = c;
955 args[ti].runs = runs;
956 args[ti].ax_count = ax_count;
957 args[ti].f = f;
958 args[ti].i_min = b.imin + ((b.imax - b.imin) * ti) / N;
959 args[ti].i_max = b.imin + ((b.imax - b.imin) * (ti + 1)) / N;
960 args[ti].j_min = b.jmin;
961 args[ti].j_max = b.jmax;
962 args[ti].subdomain = ti;
963 #ifdef USE_PTHREAD
964 pthread_attr_init(&thread_attr[ti]);
965 pthread_attr_setdetachstate(&thread_attr[ti], PTHREAD_CREATE_JOINABLE);
966 pthread_create(&threads[ti], &thread_attr[ti], diff_subdomain, &args[ti]);
967 #else
968 diff_subdomain(&args[ti]);
969 #endif
970 }
971
972 for (int ti = 0; ti < N; ti++) {
973 #ifdef USE_PTHREAD
974 pthread_join(threads[ti], NULL);
975 #endif
976 runs.back().add(args[ti].runs.back());
977 }
978
979 #ifdef USE_PTHREAD
980 free(threads);
981 free(thread_attr);
982 #endif
983
984 delete[] args;
985
986 ui::get()->d2_align_sample_stop();
987
988 }
989
990 private:
991 void rediff() {
992 std::vector<diff_trans> t_array;
993
994 for (unsigned int r = 0; r < runs.size(); r++) {
995 t_array.push_back(runs[r].offset);
996 }
997
998 runs.clear();
999
1000 for (unsigned int r = 0; r < t_array.size(); r++)
1001 diff(si, t_array[r], ax_count, frame);
1002 }
1003
1004
1005 public:
1006 int better() {
1007 assert(runs.size() >= 2);
1008 assert(runs[0].offset.scale() == runs[1].offset.scale());
1009
1010 return (runs[1].get_error() < runs[0].get_error()
1011 || (!finite(runs[0].get_error()) && finite(runs[1].get_error())));
1012 }
1013
1014 int better_defined() {
1015 assert(runs.size() >= 2);
1016 assert(runs[0].offset.scale() == runs[1].offset.scale());
1017
1018 return (runs[1].get_error() < runs[0].get_error());
1019 }
1020
1021 diff_stat_generic(transformation::elem_bounds_t e)
1022 : runs(), old_runs(), perturb_multipliers() {
1023 elem_bounds = e;
1024 }
1025
1026 run_index get_run_index(unsigned int perturb_index) {
1027 return perturb_index;
1028 }
1029
1030 run &get_run(unsigned int perturb_index) {
1031 run_index index = get_run_index(perturb_index);
1032
1033 assert(old_runs.count(index));
1034 return old_runs[index];
1035 }
1036
1037 void rescale(ale_pos scale, scale_cluster _si) {
1038 assert(runs.size() == 1);
1039
1040 si = _si;
1041
1042 runs[0].rescale(scale);
1043
1044 rediff();
1045 }
1046
1047 ~diff_stat_generic() {
1048 }
1049
1050 diff_stat_generic &operator=(const diff_stat_generic &dst) {
1051 /*
1052 * Copy run information.
1053 */
1054 runs = dst.runs;
1055 old_runs = dst.old_runs;
1056
1057 /*
1058 * Copy diff variables
1059 */
1060 si = dst.si;
1061 ax_count = dst.ax_count;
1062 frame = dst.frame;
1063 perturb_multipliers = dst.perturb_multipliers;
1064 elem_bounds = dst.elem_bounds;
1065
1066 return *this;
1067 }
1068
1069 diff_stat_generic(const diff_stat_generic &dst) : runs(), old_runs(),
1070 perturb_multipliers() {
1071 operator=(dst);
1072 }
1073
1074 void set_elem_bounds(transformation::elem_bounds_t e) {
1075 elem_bounds = e;
1076 }
1077
1078 ale_accum get_result() {
1079 assert(runs.size() == 1);
1080 return runs[0].result;
1081 }
1082
1083 ale_accum get_divisor() {
1084 assert(runs.size() == 1);
1085 return runs[0].divisor;
1086 }
1087
1088 diff_trans get_offset() {
1089 assert(runs.size() == 1);
1090 return runs[0].offset;
1091 }
1092
1093 int operator!=(diff_stat_generic &param) {
1094 return (get_error() != param.get_error());
1095 }
1096
1097 int operator==(diff_stat_generic &param) {
1098 return !(operator!=(param));
1099 }
1100
1101 ale_pos get_error_perturb() {
1102 assert(runs.size() == 1);
1103 return runs[0].get_error_perturb();
1104 }
1105
1106 ale_accum get_error() const {
1107 assert(runs.size() == 1);
1108 return runs[0].get_error();
1109 }
1110
1111 public:
1112 /*
1113 * Get the set of transformations produced by a given perturbation
1114 */
1115 void get_perturb_set(std::vector<diff_trans> *set,
1116 ale_pos adj_p, ale_pos adj_o, ale_pos adj_b,
1117 ale_pos *current_bd, ale_pos *modified_bd,
1118 std::vector<ale_pos> multipliers = std::vector<ale_pos>()) {
1119
1120 assert(runs.size() == 1);
1121
1122 diff_trans test_t(diff_trans::eu_identity());
1123
1124 /*
1125 * Translational or euclidean transformation
1126 */
1127
1128 for (unsigned int i = 0; i < 2; i++)
1129 for (unsigned int s = 0; s < 2; s++) {
1130
1131 if (!multipliers.size())
1132 multipliers.push_back(1);
1133
1134 assert(finite(multipliers[0]));
1135
1136 test_t = get_offset();
1137
1138 // test_t.eu_modify(i, (s ? -adj_p : adj_p) * multipliers[0]);
1139 test_t.translate((i ? point(1, 0) : point(0, 1))
1140 * (s ? -adj_p : adj_p)
1141 * multipliers[0]);
1142
1143 test_t.snap(adj_p / 2);
1144
1145 set->push_back(test_t);
1146 multipliers.erase(multipliers.begin());
1147
1148 }
1149
1150 if (alignment_class > 0)
1151 for (unsigned int s = 0; s < 2; s++) {
1152
1153 if (!multipliers.size())
1154 multipliers.push_back(1);
1155
1156 assert(multipliers.size());
1157 assert(finite(multipliers[0]));
1158
1159 if (!(adj_o * multipliers[0] < rot_max))
1160 return;
1161
1162 ale_pos adj_s = (s ? 1 : -1) * adj_o * multipliers[0];
1163
1164 test_t = get_offset();
1165
1166 run_index ori = get_run_index(set->size());
1167 point centroid = point::undefined();
1168
1169 if (!old_runs.count(ori))
1170 ori = get_run_index(0);
1171
1172 if (!centroid.finite() && old_runs.count(ori)) {
1173 centroid = old_runs[ori].get_error_centroid();
1174
1175 if (!centroid.finite())
1176 centroid = old_runs[ori].get_centroid();
1177
1178 centroid *= test_t.scale()
1179 / old_runs[ori].offset.scale();
1180 }
1181
1182 if (!centroid.finite() && !test_t.is_projective()) {
1183 test_t.eu_modify(2, adj_s);
1184 } else if (!centroid.finite()) {
1185 centroid = point(si.input->height() / 2,
1186 si.input->width() / 2);
1187
1188 test_t.rotate(centroid + si.accum->offset(),
1189 adj_s);
1190 } else {
1191 test_t.rotate(centroid + si.accum->offset(),
1192 adj_s);
1193 }
1194
1195 test_t.snap(adj_p / 2);
1196
1197 set->push_back(test_t);
1198 multipliers.erase(multipliers.begin());
1199 }
1200
1201 if (alignment_class == 2) {
1202
1203 /*
1204 * Projective transformation
1205 */
1206
1207 for (unsigned int i = 0; i < 4; i++)
1208 for (unsigned int j = 0; j < 2; j++)
1209 for (unsigned int s = 0; s < 2; s++) {
1210
1211 if (!multipliers.size())
1212 multipliers.push_back(1);
1213
1214 assert(multipliers.size());
1215 assert(finite(multipliers[0]));
1216
1217 ale_pos adj_s = (s ? -1 : 1) * adj_p * multipliers [0];
1218
1219 test_t = get_offset();
1220
1221 if (perturb_type == 0)
1222 test_t.gpt_modify_bounded(j, i, adj_s, elem_bounds.scale_to_bounds(si.accum->height(), si.accum->width()));
1223 else if (perturb_type == 1)
1224 test_t.gr_modify(j, i, adj_s);
1225 else
1226 assert(0);
1227
1228 test_t.snap(adj_p / 2);
1229
1230 set->push_back(test_t);
1231 multipliers.erase(multipliers.begin());
1232 }
1233
1234 }
1235
1236 /*
1237 * Barrel distortion
1238 */
1239
1240 if (bda_mult != 0 && adj_b != 0) {
1241
1242 for (unsigned int d = 0; d < get_offset().bd_count(); d++)
1243 for (unsigned int s = 0; s < 2; s++) {
1244
1245 if (!multipliers.size())
1246 multipliers.push_back(1);
1247
1248 assert (multipliers.size());
1249 assert (finite(multipliers[0]));
1250
1251 ale_pos adj_s = (s ? -1 : 1) * adj_b * multipliers[0];
1252
1253 if (bda_rate > 0 && fabs(modified_bd[d] + adj_s - current_bd[d]) > bda_rate)
1254 continue;
1255
1256 diff_trans test_t = get_offset();
1257
1258 test_t.bd_modify(d, adj_s);
1259
1260 set->push_back(test_t);
1261 }
1262 }
1263 }
1264
1265 void confirm() {
1266 assert(runs.size() == 2);
1267 runs[0] = runs[1];
1268 runs.pop_back();
1269 }
1270
1271 void discard() {
1272 assert(runs.size() == 2);
1273 runs.pop_back();
1274 }
1275
1276 void perturb_test(ale_pos perturb, ale_pos adj_p, ale_pos adj_o, ale_pos adj_b,
1277 ale_pos *current_bd, ale_pos *modified_bd, int stable) {
1278
1279 assert(runs.size() == 1);
1280
1281 std::vector<diff_trans> t_set;
1282
1283 if (perturb_multipliers.size() == 0) {
1284 get_perturb_set(&t_set, adj_p, adj_o, adj_b,
1285 current_bd, modified_bd);
1286
1287 for (unsigned int i = 0; i < t_set.size(); i++) {
1288 diff_stat_generic test = *this;
1289
1290 test.diff(si, t_set[i], ax_count, frame);
1291
1292 test.confirm();
1293
1294 if (finite(test.get_error_perturb()))
1295 perturb_multipliers.push_back(adj_p / test.get_error_perturb());
1296 else
1297 perturb_multipliers.push_back(1);
1298
1299 }
1300
1301 t_set.clear();
1302 }
1303
1304 get_perturb_set(&t_set, adj_p, adj_o, adj_b, current_bd, modified_bd,
1305 perturb_multipliers);
1306
1307 int found_unreliable = 1;
1308 std::vector<int> tested(t_set.size(), 0);
1309
1310 for (unsigned int i = 0; i < t_set.size(); i++) {
1311 run_index ori = get_run_index(i);
1312
1313 /*
1314 * Check for stability
1315 */
1316 if (stable
1317 && old_runs.count(ori)
1318 && old_runs[ori].offset == t_set[i])
1319 tested[i] = 1;
1320 }
1321
1322 std::vector<ale_pos> perturb_multipliers_original = perturb_multipliers;
1323
1324 while (found_unreliable) {
1325
1326 found_unreliable = 0;
1327
1328 for (unsigned int i = 0; i < t_set.size(); i++) {
1329
1330 if (tested[i])
1331 continue;
1332
1333 diff(si, t_set[i], ax_count, frame);
1334
1335 if (!(i < perturb_multipliers.size())
1336 || !finite(perturb_multipliers[i])) {
1337
1338 perturb_multipliers.resize(i + 1);
1339
1340 if (finite(perturb_multipliers[i])
1341 && finite(adj_p)
1342 && finite(adj_p / runs[1].get_error_perturb())) {
1343 perturb_multipliers[i] =
1344 adj_p / runs[1].get_error_perturb();
1345
1346 found_unreliable = 1;
1347 } else
1348 perturb_multipliers[i] = 1;
1349
1350 continue;
1351 }
1352
1353 run_index ori = get_run_index(i);
1354
1355 if (old_runs.count(ori) == 0)
1356 old_runs.insert(std::pair<run_index, run>(ori, runs[1]));
1357 else
1358 old_runs[ori] = runs[1];
1359
1360 if (finite(perturb_multipliers_original[i])
1361 && finite(runs[1].get_error_perturb())
1362 && finite(adj_p)
1363 && finite(perturb_multipliers_original[i] * adj_p / runs[1].get_error_perturb()))
1364 perturb_multipliers[i] = perturb_multipliers_original[i]
1365 * adj_p / runs[1].get_error_perturb();
1366 else
1367 perturb_multipliers[i] = 1;
1368
1369 tested[i] = 1;
1370
1371 if (better()
1372 && runs[1].get_error() < runs[0].get_error()
1373 && perturb_multipliers[i]
1374 / perturb_multipliers_original[i] < 2) {
1375 runs[0] = runs[1];
1376 runs.pop_back();
1377 return;
1378 }
1379
1380 }
1381
1382 if (runs.size() > 1)
1383 runs.pop_back();
1384
1385 if (!found_unreliable)
1386 return;
1387 }
1388 }
1389
1390 };
1391
1392 typedef diff_stat_generic<trans_single> diff_stat_t;
1393 typedef diff_stat_generic<trans_multi> diff_stat_multi;
1394
1395
1396 /*
1397 * Adjust exposure for an aligned frame B against reference A.
1398 *
1399 * Expects full-LOD images.
1400 *
1401 * Note: This method does not use any weighting, by certainty or
1402 * otherwise, in the first exposure registration pass, as any bias of
1403 * weighting according to color may also bias the exposure registration
1404 * result; it does use weighting, including weighting by certainty
1405 * (even if certainty weighting is not specified), in the second pass,
1406 * under the assumption that weighting by certainty improves handling
1407 * of out-of-range highlights, and that bias of exposure measurements
1408 * according to color may generally be less harmful after spatial
1409 * registration has been performed.
1410 */
1411 class exposure_ratio_iterate : public thread::decompose_domain {
1412 pixel_accum *asums;
1413 pixel_accum *bsums;
1414 pixel_accum *asum;
1415 pixel_accum *bsum;
1416 struct scale_cluster c;
1417 const transformation &t;
1418 int ax_count;
1419 int pass_number;
1420 protected:
1421 void prepare_subdomains(unsigned int N) {
1422 asums = new pixel_accum[N];
1423 bsums = new pixel_accum[N];
1424 }
1425 void subdomain_algorithm(unsigned int thread,
1426 int i_min, int i_max, int j_min, int j_max) {
1427
1428 point offset = c.accum->offset();
1429
1430 for (mc_iterate m(i_min, i_max, j_min, j_max, thread); !m.done(); m++) {
1431
1432 unsigned int i = (unsigned int) m.get_i();
1433 unsigned int j = (unsigned int) m.get_j();
1434
1435 if (ref_excluded(i, j, offset, c.ax_parameters, ax_count))
1436 continue;
1437
1438 /*
1439 * Transform
1440 */
1441
1442 struct point q;
1443
1444 q = (c.input_scale < 1.0 && interpolant == NULL)
1445 ? t.scaled_inverse_transform(
1446 point(i + offset[0], j + offset[1]))
1447 : t.unscaled_inverse_transform(
1448 point(i + offset[0], j + offset[1]));
1449
1450 /*
1451 * Check that the transformed coordinates are within
1452 * the boundaries of array c.input, that they are not
1453 * subject to exclusion, and that the weight value in
1454 * the accumulated array is nonzero.
1455 */
1456
1457 if (input_excluded(q[0], q[1], c.ax_parameters, ax_count))
1458 continue;
1459
1460 if (q[0] >= 0
1461 && q[0] <= c.input->height() - 1
1462 && q[1] >= 0
1463 && q[1] <= c.input->width() - 1
1464 && ((pixel) c.certainty->get_pixel(i, j)).minabs_norm() != 0) {
1465 pixel a = c.accum->get_pixel(i, j);
1466 pixel b;
1467
1468 b = c.input->get_bl(q);
1469
1470 pixel weight = ((c.aweight && pass_number)
1471 ? (pixel) c.aweight->get_pixel(i, j)
1472 : pixel(1, 1, 1))
1473 * (pass_number
1474 ? psqrt(c.certainty->get_pixel(i, j)
1475 * c.input_certainty->get_bl(q, 1))
1476 : pixel(1, 1, 1));
1477
1478 asums[thread] += a * weight;
1479 bsums[thread] += b * weight;
1480 }
1481 }
1482 }
1483
1484 void finish_subdomains(unsigned int N) {
1485 for (unsigned int n = 0; n < N; n++) {
1486 *asum += asums[n];
1487 *bsum += bsums[n];
1488 }
1489 delete[] asums;
1490 delete[] bsums;
1491 }
1492 public:
1493 exposure_ratio_iterate(pixel_accum *_asum,
1494 pixel_accum *_bsum,
1495 struct scale_cluster _c,
1496 const transformation &_t,
1497 int _ax_count,
1498 int _pass_number) : decompose_domain(0, _c.accum->height(),
1499 0, _c.accum->width()),
1500 t(_t) {
1501
1502 asum = _asum;
1503 bsum = _bsum;
1504 c = _c;
1505 ax_count = _ax_count;
1506 pass_number = _pass_number;
1507 }
1508
1509 exposure_ratio_iterate(pixel_accum *_asum,
1510 pixel_accum *_bsum,
1511 struct scale_cluster _c,
1512 const transformation &_t,
1513 int _ax_count,
1514 int _pass_number,
1515 transformation::elem_bounds_int_t b) : decompose_domain(b.imin, b.imax,
1516 b.jmin, b.jmax),
1517 t(_t) {
1518
1519 asum = _asum;
1520 bsum = _bsum;
1521 c = _c;
1522 ax_count = _ax_count;
1523 pass_number = _pass_number;
1524 }
1525 };
1526
1527 static void set_exposure_ratio(unsigned int m, struct scale_cluster c,
1528 const transformation &t, int ax_count, int pass_number) {
1529
1530 if (_exp_register == 2) {
1531 /*
1532 * Use metadata only.
1533 */
1534 ale_real gain_multiplier = image_rw::exp(m).get_gain_multiplier();
1535 pixel multiplier = pixel(gain_multiplier, gain_multiplier, gain_multiplier);
1536
1537 image_rw::exp(m).set_multiplier(multiplier);
1538 ui::get()->exp_multiplier(multiplier[0],
1539 multiplier[1],
1540 multiplier[2]);
1541
1542 return;
1543 }
1544
1545 pixel_accum asum(0, 0, 0), bsum(0, 0, 0);
1546
1547 exposure_ratio_iterate eri(&asum, &bsum, c, t, ax_count, pass_number);
1548 eri.run();
1549
1550 // std::cerr << (asum / bsum) << " ";
1551
1552 pixel_accum new_multiplier;
1553
1554 new_multiplier = asum / bsum * image_rw::exp(m).get_multiplier();
1555
1556 if (finite(new_multiplier[0])
1557 && finite(new_multiplier[1])
1558 && finite(new_multiplier[2])) {
1559 image_rw::exp(m).set_multiplier(new_multiplier);
1560 ui::get()->exp_multiplier(new_multiplier[0],
1561 new_multiplier[1],
1562 new_multiplier[2]);
1563 }
1564 }
1565
1566 /*
1567 * Copy all ax parameters.
1568 */
1569 static exclusion *copy_ax_parameters(int local_ax_count, exclusion *source) {
1570
1571 exclusion *dest = (exclusion *) malloc(local_ax_count * sizeof(exclusion));
1572
1573 assert (dest);
1574
1575 if (!dest)
1576 ui::get()->memory_error("exclusion regions");
1577
1578 for (int idx = 0; idx < local_ax_count; idx++)
1579 dest[idx] = source[idx];
1580
1581 return dest;
1582 }
1583
1584 /*
1585 * Copy ax parameters according to frame.
1586 */
1587 static exclusion *filter_ax_parameters(int frame, int *local_ax_count) {
1588
1589 exclusion *dest = (exclusion *) malloc(ax_count * sizeof(exclusion));
1590
1591 assert (dest);
1592
1593 if (!dest)
1594 ui::get()->memory_error("exclusion regions");
1595
1596 *local_ax_count = 0;
1597
1598 for (int idx = 0; idx < ax_count; idx++) {
1599 if (ax_parameters[idx].x[4] > frame
1600 || ax_parameters[idx].x[5] < frame)
1601 continue;
1602
1603 dest[*local_ax_count] = ax_parameters[idx];
1604
1605 (*local_ax_count)++;
1606 }
1607
1608 return dest;
1609 }
1610
1611 static void scale_ax_parameters(int local_ax_count, exclusion *ax_parameters,
1612 ale_pos ref_scale, ale_pos input_scale) {
1613 for (int i = 0; i < local_ax_count; i++) {
1614 ale_pos scale = (ax_parameters[i].type == exclusion::RENDER)
1615 ? ref_scale
1616 : input_scale;
1617
1618 for (int n = 0; n < 6; n++) {
1619 ax_parameters[i].x[n] = ax_parameters[i].x[n] * scale;
1620 }
1621 }
1622 }
1623
1624 /*
1625 * Prepare the next level of detail for ordinary images.
1626 */
1627 static const image *prepare_lod(const image *current) {
1628 if (current == NULL)
1629 return NULL;
1630
1631 return current->scale_by_half("prepare_lod");
1632 }
1633
1634 /*
1635 * Prepare the next level of detail for definition maps.
1636 */
1637 static const image *prepare_lod_def(const image *current) {
1638 if (current == NULL)
1639 return NULL;
1640
1641 return current->defined_scale_by_half("prepare_lod_def");
1642 }
1643
1644 /*
1645 * Initialize scale cluster data structures.
1646 */
1647
1648 static void init_nl_cluster(struct scale_cluster *sc) {
1649 }
1650
1651 /*
1652 * Destroy the first element in the scale cluster data structure.
1653 */
1654 static void final_clusters(struct scale_cluster *scale_clusters, ale_pos scale_factor,
1655 unsigned int steps) {
1656
1657 if (scale_clusters[0].input_scale < 1.0) {
1658 delete scale_clusters[0].input;
1659 }
1660
1661 delete scale_clusters[0].input_certainty;
1662
1663 free((void *)scale_clusters[0].ax_parameters);
1664
1665 for (unsigned int step = 1; step < steps; step++) {
1666 delete scale_clusters[step].accum;
1667 delete scale_clusters[step].certainty;
1668 delete scale_clusters[step].aweight;
1669
1670 if (scale_clusters[step].input_scale < 1.0) {
1671 delete scale_clusters[step].input;
1672 delete scale_clusters[step].input_certainty;
1673 }
1674
1675 free((void *)scale_clusters[step].ax_parameters);
1676 }
1677
1678 free(scale_clusters);
1679 }
1680
1681 /*
1682 * Calculate the centroid of a control point for the set of frames
1683 * having index lower than m. Divide by any scaling of the output.
1684 */
1685 static point unscaled_centroid(unsigned int m, unsigned int p) {
1686 assert(_keep);
1687
1688 point point_sum(0, 0);
1689 ale_accum divisor = 0;
1690
1691 for(unsigned int j = 0; j < m; j++) {
1692 point pp = cp_array[p][j];
1693
1694 if (pp.defined()) {
1695 point_sum += kept_t[j].transform_unscaled(pp)
1696 / kept_t[j].scale();
1697 divisor += 1;
1698 }
1699 }
1700
1701 if (divisor == 0)
1702 return point::undefined();
1703
1704 return point_sum / divisor;
1705 }
1706
1707 /*
1708 * Calculate centroid of this frame, and of all previous frames,
1709 * from points common to both sets.
1710 */
1711 static void centroids(unsigned int m, point *current, point *previous) {
1712 /*
1713 * Calculate the translation
1714 */
1715 point other_centroid(0, 0);
1716 point this_centroid(0, 0);
1717 ale_pos divisor = 0;
1718
1719 for (unsigned int i = 0; i < cp_count; i++) {
1720 point other_c = unscaled_centroid(m, i);
1721 point this_c = cp_array[i][m];
1722
1723 if (!other_c.defined() || !this_c.defined())
1724 continue;
1725
1726 other_centroid += other_c;
1727 this_centroid += this_c;
1728 divisor += 1;
1729
1730 }
1731
1732 if (divisor == 0) {
1733 *current = point::undefined();
1734 *previous = point::undefined();
1735 return;
1736 }
1737
1738 *current = this_centroid / divisor;
1739 *previous = other_centroid / divisor;
1740 }
1741
1742 /*
1743 * Calculate the RMS error of control points for frame m, with
1744 * transformation t, against control points for earlier frames.
1745 */
1746 static ale_pos cp_rms_error(unsigned int m, transformation t) {
1747 assert (_keep);
1748
1749 ale_accum err = 0;
1750 ale_accum divisor = 0;
1751
1752 for (unsigned int i = 0; i < cp_count; i++)
1753 for (unsigned int j = 0; j < m; j++) {
1754 const point *p = cp_array[i];
1755 point p_ref = kept_t[j].transform_unscaled(p[j]);
1756 point p_cur = t.transform_unscaled(p[m]);
1757
1758 if (!p_ref.defined() || !p_cur.defined())
1759 continue;
1760
1761 err += p_ref.lengthtosq(p_cur);
1762 divisor += 1;
1763 }
1764
1765 return (ale_pos) sqrt(err / divisor);
1766 }
1767
1768
1769 static void test_global(diff_stat_t *here, scale_cluster si, transformation t,
1770 int local_ax_count, int m, ale_accum local_gs_mo, ale_pos perturb) {
1771
1772 diff_stat_t test(*here);
1773
1774 test.diff(si, t.get_current_element(), local_ax_count, m);
1775
1776 unsigned int ovl = overlap(si, t, local_ax_count);
1777
1778 if (ovl >= local_gs_mo && test.better()) {
1779 test.confirm();
1780 *here = test;
1781 ui::get()->set_match(here->get_error());
1782 ui::get()->set_offset(here->get_offset());
1783 } else {
1784 test.discard();
1785 }
1786
1787 }
1788
1789 /*
1790 * Get the set of global transformations for a given density
1791 */
1792 static void test_globals(diff_stat_t *here,
1793 scale_cluster si, transformation t, int local_gs, ale_pos adj_p,
1794 int local_ax_count, int m, ale_accum local_gs_mo, ale_pos perturb) {
1795
1796 transformation offset = t;
1797
1798 point min, max;
1799
1800 transformation offset_p = offset;
1801
1802 if (!offset_p.is_projective())
1803 offset_p.eu_to_gpt();
1804
1805 min = max = offset_p.gpt_get(0);
1806 for (int p_index = 1; p_index < 4; p_index++) {
1807 point p = offset_p.gpt_get(p_index);
1808 if (p[0] < min[0])
1809 min[0] = p[0];
1810 if (p[1] < min[1])
1811 min[1] = p[1];
1812 if (p[0] > max[0])
1813 max[0] = p[0];
1814 if (p[1] > max[1])
1815 max[1] = p[1];
1816 }
1817
1818 point inner_min_t = -min;
1819 point inner_max_t = -max + point(si.accum->height(), si.accum->width());
1820 point outer_min_t = -max + point(adj_p - 1, adj_p - 1);
1821 point outer_max_t = point(si.accum->height(), si.accum->width()) - point(adj_p, adj_p);
1822
1823 if (local_gs == 1 || local_gs == 3 || local_gs == 4 || local_gs == 6) {
1824
1825 /*
1826 * Inner
1827 */
1828
1829 for (ale_pos i = inner_min_t[0]; i <= inner_max_t[0]; i += adj_p)
1830 for (ale_pos j = inner_min_t[1]; j <= inner_max_t[1]; j += adj_p) {
1831 transformation test_t = offset;
1832 test_t.translate(point(i, j));
1833 test_global(here, si, test_t, local_ax_count, m, local_gs_mo, perturb);
1834 }
1835 }
1836
1837 if (local_gs == 2 || local_gs == 3 || local_gs == -1 || local_gs == 6) {
1838
1839 /*
1840 * Outer
1841 */
1842
1843 for (ale_pos i = outer_min_t[0]; i <= outer_max_t[0]; i += adj_p)
1844 for (ale_pos j = outer_min_t[1]; j < inner_min_t[1]; j += adj_p) {
1845 transformation test_t = offset;
1846 test_t.translate(point(i, j));
1847 test_global(here, si, test_t, local_ax_count, m, local_gs_mo, perturb);
1848 }
1849 for (ale_pos i = outer_min_t[0]; i <= outer_max_t[0]; i += adj_p)
1850 for (ale_pos j = outer_max_t[1]; j > inner_max_t[1]; j -= adj_p) {
1851 transformation test_t = offset;
1852 test_t.translate(point(i, j));
1853 test_global(here, si, test_t, local_ax_count, m, local_gs_mo, perturb);
1854 }
1855 for (ale_pos i = outer_min_t[0]; i < inner_min_t[0]; i += adj_p)
1856 for (ale_pos j = outer_min_t[1]; j <= outer_max_t[1]; j += adj_p) {
1857 transformation test_t = offset;
1858 test_t.translate(point(i, j));
1859 test_global(here, si, test_t, local_ax_count, m, local_gs_mo, perturb);
1860 }
1861 for (ale_pos i = outer_max_t[0]; i > inner_max_t[0]; i -= adj_p)
1862 for (ale_pos j = outer_min_t[1]; j <= outer_max_t[1]; j += adj_p) {
1863 transformation test_t = offset;
1864 test_t.translate(point(i, j));
1865 test_global(here, si, test_t, local_ax_count, m, local_gs_mo, perturb);
1866 }
1867 }
1868 }
1869
1870 static void get_translational_set(std::vector<transformation> *set,
1871 transformation t, ale_pos adj_p) {
1872
1873 ale_pos adj_s;
1874
1875 transformation offset = t;
1876 transformation test_t(transformation::eu_identity());
1877
1878 for (int i = 0; i < 2; i++)
1879 for (adj_s = -adj_p; adj_s <= adj_p; adj_s += 2 * adj_p) {
1880
1881 test_t = offset;
1882
1883 test_t.translate(i ? point(adj_s, 0) : point(0, adj_s));
1884
1885 set->push_back(test_t);
1886 }
1887 }
1888
1889 static int threshold_ok(ale_accum error) {
1890 if ((1 - error) * (ale_accum) 100 >= match_threshold)
1891 return 1;
1892
1893 if (!(match_threshold >= 0))
1894 return 1;
1895
1896 return 0;
1897 }
1898
1899 /*
1900 * Check for satisfaction of the certainty threshold.
1901 */
1902 static int ma_cert_satisfied(const scale_cluster &c, const transformation &t, unsigned int i) {
1903 transformation::elem_bounds_int_t b = t.elem_bounds().scale_to_bounds(c.accum->height(), c.accum->width());
1904 ale_accum sum[3] = {0, 0, 0};
1905
1906 for (unsigned int ii = b.imin; ii < b.imax; ii++)
1907 for (unsigned int jj = b.jmin; jj < b.jmax; jj++) {
1908 pixel p = c.accum->get_pixel(ii, jj);
1909 sum[0] += p[0];
1910 sum[1] += p[1];
1911 sum[2] += p[2];
1912 }
1913
1914 unsigned int count = (b.jmax - b.jmin) * (b.imax - b.imin);
1915
1916 for (int k = 0; k < 3; k++)
1917 if (sum[k] / count < _ma_cert)
1918 return 0;
1919
1920 return 1;
1921 }
1922
1923 static diff_stat_t check_ancestor_path(const trans_multi &offset, const scale_cluster &si, diff_stat_t here, int local_ax_count, int frame) {
1924
1925 if (offset.get_current_index() > 0)
1926 for (int i = offset.parent_index(offset.get_current_index()); i >= 0; i = (i > 0) ? (int) offset.parent_index(i) : -1) {
1927 trans_single t = offset.get_element(i);
1928 t.rescale(offset.get_current_element().scale());
1929
1930 here.diff(si, t, local_ax_count, frame);
1931
1932 if (here.better_defined())
1933 here.confirm();
1934 else
1935 here.discard();
1936 }
1937
1938 return here;
1939 }
1940
1941 #ifdef USE_FFTW
1942 /*
1943 * High-pass filter for frequency weights
1944 */
1945 static void hipass(int rows, int cols, fftw_complex *inout) {
1946 for (int i = 0; i < rows * vert_freq_cut; i++)
1947 for (int j = 0; j < cols; j++)
1948 for (int k = 0; k < 2; k++)
1949 inout[i * cols + j][k] = 0;
1950 for (int i = 0; i < rows; i++)
1951 for (int j = 0; j < cols * horiz_freq_cut; j++)
1952 for (int k = 0; k < 2; k++)
1953 inout[i * cols + j][k] = 0;
1954 for (int i = 0; i < rows; i++)
1955 for (int j = 0; j < cols; j++)
1956 for (int k = 0; k < 2; k++)
1957 if (i / (double) rows + j / (double) cols < 2 * avg_freq_cut)
1958 inout[i * cols + j][k] = 0;
1959 }
1960 #endif
1961
1962
1963 /*
1964 * Reset alignment weights
1965 */
1966 static void reset_weights() {
1967 if (alignment_weights != NULL)
1968 ale_image_release(alignment_weights);
1969
1970 alignment_weights = NULL;
1971 }
1972
1973 /*
1974 * Initialize alignment weights
1975 */
1976 static void init_weights() {
1977 if (alignment_weights != NULL)
1978 return;
1979
1980 alignment_weights = ale_new_image(accel::context(), ALE_IMAGE_RGB, ale_image_get_type(reference_image));
1981
1982 assert (alignment_weights);
1983
1984 assert (!ale_resize_image(alignment_weights, 0, 0, ale_image_get_width(reference_image), ale_image_get_height(reference_image)));
1985
1986 ale_image_map_0(alignment_weights, "SET_PIXEL(p, PIXEL(1, 1, 1))");
1987 }
1988
1989 /*
1990 * Update alignment weights with weight map
1991 */
1992 static void map_update() {
1993
1994 if (weight_map == NULL)
1995 return;
1996
1997 init_weights();
1998
1999 ale_image_map_2(alignment_weights, alignment_weights, weight_map, " \
2000 SET_PIXEL(p, GET_PIXEL(0, p) * GET_PIXEL_BG(1, p))");
2001 }
2002
2003 /*
2004 * Update alignment weights with algorithmic weights
2005 */
2006 static void wmx_update() {
2007 #ifdef USE_UNIX
2008
2009 static exposure *exp_def = new exposure_default();
2010 static exposure *exp_bool = new exposure_boolean();
2011
2012 if (wmx_file == NULL || wmx_exec == NULL || wmx_defs == NULL)
2013 return;
2014
2015 unsigned int rows = ale_image_get_height(reference_image);
2016 unsigned int cols = ale_image_get_width(reference_image);
2017
2018 image_rw::write_image(wmx_file, reference_image);
2019 image_rw::write_image(wmx_defs, reference_image, exp_bool->get_gamma(), 0, 0, 1);
2020
2021 /* execute ... */
2022 int exit_status = 1;
2023 if (!fork()) {
2024 execlp(wmx_exec, wmx_exec, wmx_file, wmx_defs, NULL);
2025 ui::get()->exec_failure(wmx_exec, wmx_file, wmx_defs);
2026 }
2027
2028 wait(&exit_status);
2029
2030 if (exit_status)
2031 ui::get()->fork_failure("d2::align");
2032
2033 ale_image wmx_weights = image_rw::read_image(wmx_file, exp_def);
2034
2035 ale_image_set_x_offset(wmx_weights, ale_image_get_x_offset(reference_image));
2036 ale_image_set_y_offset(wmx_weights, ale_image_get_y_offset(reference_image));
2037
2038 if (ale_image_get_height(wmx_weights) != rows || ale_image_get_width(wmx_weights) != cols)
2039 ui::get()->error("algorithmic weighting must not change image size");
2040
2041 if (alignment_weights == NULL)
2042 alignment_weights = wmx_weights;
2043 else
2044 ale_image_map_2(alignment_weights, alignment_weights, wmx_weights, "\
2045 SET_PIXEL(p, GET_PIXEL(0, p) * GET_PIXEL(1, p))");
2046 #endif
2047 }
2048
2049 /*
2050 * Update alignment weights with frequency weights
2051 */
2052 static void fw_update() {
2053 #ifdef USE_FFTW
2054 if (horiz_freq_cut == 0
2055 && vert_freq_cut == 0
2056 && avg_freq_cut == 0)
2057 return;
2058
2059 /*
2060 * Required for correct operation of --fwshow
2061 */
2062
2063 assert (alignment_weights == NULL);
2064
2065 int rows = reference_image->height();
2066 int cols = reference_image->width();
2067 int colors = reference_image->depth();
2068
2069 alignment_weights = new_image_ale_real(rows, cols,
2070 colors, "alignment_weights");
2071
2072 fftw_complex *inout = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * rows * cols);
2073
2074 assert (inout);
2075
2076 fftw_plan pf = fftw_plan_dft_2d(rows, cols,
2077 inout, inout,
2078 FFTW_FORWARD, FFTW_ESTIMATE);
2079
2080 fftw_plan pb = fftw_plan_dft_2d(rows, cols,
2081 inout, inout,
2082 FFTW_BACKWARD, FFTW_ESTIMATE);
2083
2084 for (int k = 0; k < colors; k++) {
2085 for (int i = 0; i < rows * cols; i++) {
2086 inout[i][0] = reference_image->get_pixel(i / cols, i % cols)[k];
2087 inout[i][1] = 0;
2088 }
2089
2090 fftw_execute(pf);
2091 hipass(rows, cols, inout);
2092 fftw_execute(pb);
2093
2094 for (int i = 0; i < rows * cols; i++) {
2095 #if 0
2096 alignment_weights->pix(i / cols, i % cols)[k] = fabs(inout[i][0] / (rows * cols));
2097 #else
2098 alignment_weights->set_chan(i / cols, i % cols, k,
2099 sqrt(pow(inout[i][0] / (rows * cols), 2)
2100 + pow(inout[i][1] / (rows * cols), 2)));
2101 #endif
2102 }
2103 }
2104
2105 fftw_destroy_plan(pf);
2106 fftw_destroy_plan(pb);
2107 fftw_free(inout);
2108
2109 if (fw_output != NULL)
2110 image_rw::write_image(fw_output, alignment_weights);
2111 #endif
2112 }
2113
2114 /*
2115 * Update alignment to frame N.
2116 */
2117 static void update_to(int n) {
2118
2119 assert (n <= latest + 1);
2120 assert (n >= 0);
2121
2122 ale_align_properties astate = align_properties();
2123
2124 ui::get()->set_frame_num(n);
2125
2126 if (latest < 0) {
2127
2128 /*
2129 * Handle the initial frame
2130 */
2131
2132 astate.set_input_frame(image_rw::open(n));
2133
2134 const image *i = astate.get_input_frame();
2135 int is_default;
2136 transformation result = alignment_class == 2
2137 ? transformation::gpt_identity(i, scale_factor)
2138 : transformation::eu_identity(i, scale_factor);
2139 result = tload_first(tload, alignment_class == 2, result, &is_default);
2140 tsave_first(tsave, result, alignment_class == 2);
2141
2142 if (_keep > 0) {
2143 kept_t = transformation::new_eu_identity_array(image_rw::count());
2144 kept_ok = (int *) malloc(image_rw::count()
2145 * sizeof(int));
2146 assert (kept_t);
2147 assert (kept_ok);
2148
2149 if (!kept_t || !kept_ok)
2150 ui::get()->memory_error("alignment");
2151
2152 kept_ok[0] = 1;
2153 kept_t[0] = result;
2154 }
2155
2156 latest = 0;
2157 latest_ok = 1;
2158 latest_t = result;
2159
2160 astate.set_default(result);
2161 orig_t = result;
2162
2163 image_rw::close(n);
2164 }
2165
2166 for (int i = latest + 1; i <= n; i++) {
2167
2168 /*
2169 * Handle supplemental frames.
2170 */
2171
2172 assert (reference != NULL);
2173
2174 ui::get()->set_arender_current();
2175 reference->sync(i - 1);
2176 ui::get()->clear_arender_current();
2177 reference_image = reference->get_image();
2178 reference_defined = reference->get_defined();
2179
2180 if (i == 1)
2181 astate.default_set_original_bounds(reference_image);
2182
2183 reset_weights();
2184 fw_update();
2185 wmx_update();
2186 map_update();
2187
2188 assert (reference_image != NULL);
2189 assert (reference_defined != NULL);
2190
2191 astate.set_input_frame(image_rw::open(i));
2192
2193 _align(i, _gs, &astate);
2194
2195 image_rw::close(n);
2196 }
2197 }
2198
2199 public:
2200
2201 /*
2202 * Set the control point count
2203 */
2204 static void set_cp_count(unsigned int n) {
2205 assert (cp_array == NULL);
2206
2207 cp_count = n;
2208 cp_array = (const point **) malloc(n * sizeof(const point *));
2209 }
2210
2211 /*
2212 * Set control points.
2213 */
2214 static void set_cp(unsigned int i, const point *p) {
2215 cp_array[i] = p;
2216 }
2217
2218 /*
2219 * Register exposure
2220 */
2221 static void exp_register() {
2222 _exp_register = 1;
2223 }
2224
2225 /*
2226 * Register exposure only based on metadata
2227 */
2228 static void exp_meta_only() {
2229 _exp_register = 2;
2230 }
2231
2232 /*
2233 * Don't register exposure
2234 */
2235 static void exp_noregister() {
2236 _exp_register = 0;
2237 }
2238
2239 /*
2240 * Set alignment class to translation only. Only adjust the x and y
2241 * position of images. Do not rotate input images or perform
2242 * projective transformations.
2243 */
2244 static void class_translation() {
2245 alignment_class = 0;
2246 }
2247
2248 /*
2249 * Set alignment class to Euclidean transformations only. Adjust the x
2250 * and y position of images and the orientation of the image about the
2251 * image center.
2252 */
2253 static void class_euclidean() {
2254 alignment_class = 1;
2255 }
2256
2257 /*
2258 * Set alignment class to perform general projective transformations.
2259 * See the file gpt.h for more information about general projective
2260 * transformations.
2261 */
2262 static void class_projective() {
2263 alignment_class = 2;
2264 }
2265
2266 /*
2267 * Set the default initial alignment to the identity transformation.
2268 */
2269 static void initial_default_identity() {
2270 default_initial_alignment_type = 0;
2271 }
2272
2273 /*
2274 * Set the default initial alignment to the most recently matched
2275 * frame's final transformation.
2276 */
2277 static void initial_default_follow() {
2278 default_initial_alignment_type = 1;
2279 }
2280
2281 /*
2282 * Perturb output coordinates.
2283 */
2284 static void perturb_output() {
2285 perturb_type = 0;
2286 }
2287
2288 /*
2289 * Perturb source coordinates.
2290 */
2291 static void perturb_source() {
2292 perturb_type = 1;
2293 }
2294
2295 /*
2296 * Frames under threshold align optimally
2297 */
2298 static void fail_optimal() {
2299 is_fail_default = 0;
2300 }
2301
2302 /*
2303 * Frames under threshold keep their default alignments.
2304 */
2305 static void fail_default() {
2306 is_fail_default = 1;
2307 }
2308
2309 /*
2310 * Align images with an error contribution from each color channel.
2311 */
2312 static void all() {
2313 channel_alignment_type = 0;
2314 }
2315
2316 /*
2317 * Align images with an error contribution only from the green channel.
2318 * Other color channels do not affect alignment.
2319 */
2320 static void green() {
2321 channel_alignment_type = 1;
2322 }
2323
2324 /*
2325 * Align images using a summation of channels. May be useful when
2326 * dealing with images that have high frequency color ripples due to
2327 * color aliasing.
2328 */
2329 static void sum() {
2330 channel_alignment_type = 2;
2331 }
2332
2333 /*
2334 * Error metric exponent
2335 */
2336
2337 static void set_metric_exponent(float me) {
2338 metric_exponent = me;
2339 }
2340
2341 /*
2342 * Match threshold
2343 */
2344
2345 static void set_match_threshold(float mt) {
2346 match_threshold = mt;
2347 }
2348
2349 /*
2350 * Perturbation lower and upper bounds.
2351 */
2352
2353 static void set_perturb_lower(ale_pos pl, int plp) {
2354 perturb_lower = pl;
2355 perturb_lower_percent = plp;
2356 }
2357
2358 static void set_perturb_upper(ale_pos pu, int pup) {
2359 perturb_upper = pu;
2360 perturb_upper_percent = pup;
2361 }
2362
2363 /*
2364 * Maximum rotational perturbation.
2365 */
2366
2367 static void set_rot_max(int rm) {
2368
2369 /*
2370 * Obtain the largest power of two not larger than rm.
2371 */
2372
2373 rot_max = pow(2, floor(log(rm) / log(2)));
2374 }
2375
2376 /*
2377 * Barrel distortion adjustment multiplier
2378 */
2379
2380 static void set_bda_mult(ale_pos m) {
2381 bda_mult = m;
2382 }
2383
2384 /*
2385 * Barrel distortion maximum rate of change
2386 */
2387
2388 static void set_bda_rate(ale_pos m) {
2389 bda_rate = m;
2390 }
2391
2392 /*
2393 * Level-of-detail
2394 */
2395
2396 static void set_lod_preferred(int lm) {
2397 lod_preferred = lm;
2398 }
2399
2400 /*
2401 * Minimum dimension for reduced level-of-detail.
2402 */
2403
2404 static void set_min_dimension(int md) {
2405 min_dimension = md;
2406 }
2407
2408 /*
2409 * Set the scale factor
2410 */
2411 static void set_scale(ale_pos s) {
2412 scale_factor = s;
2413 }
2414
2415 /*
2416 * Set reference rendering to align against
2417 */
2418 static void set_reference(render *r) {
2419 reference = r;
2420 }
2421
2422 /*
2423 * Set the interpolant
2424 */
2425 static void set_interpolant(filter::scaled_filter *f) {
2426 interpolant = f;
2427 }
2428
2429 /*
2430 * Set alignment weights image
2431 */
2432 static void set_weight_map(const image *i) {
2433 weight_map = i;
2434 }
2435
2436 /*
2437 * Set frequency cuts
2438 */
2439 static void set_frequency_cut(double h, double v, double a) {
2440 horiz_freq_cut = h;
2441 vert_freq_cut = v;
2442 avg_freq_cut = a;
2443 }
2444
2445 /*
2446 * Set algorithmic alignment weighting
2447 */
2448 static void set_wmx(const char *e, const char *f, const char *d) {
2449 wmx_exec = e;
2450 wmx_file = f;
2451 wmx_defs = d;
2452 }
2453
2454 /*
2455 * Show frequency weights
2456 */
2457 static void set_fl_show(const char *filename) {
2458 fw_output = filename;
2459 }
2460
2461 /*
2462 * Set transformation file loader.
2463 */
2464 static void set_tload(tload_t *tl) {
2465 tload = tl;
2466 }
2467
2468 /*
2469 * Set transformation file saver.
2470 */
2471 static void set_tsave(tsave_t *ts) {
2472 tsave = ts;
2473 }
2474
2475 /*
2476 * Get match statistics for frame N.
2477 */
2478 static int match(int n) {
2479 update_to(n);
2480
2481 if (n == latest)
2482 return latest_ok;
2483 else if (_keep)
2484 return kept_ok[n];
2485 else {
2486 assert(0);
2487 exit(1);
2488 }
2489 }
2490
2491 /*
2492 * Message that old alignment data should be kept.
2493 */
2494 static void keep() {
2495 assert (latest == -1);
2496 _keep = 1;
2497 }
2498
2499 /*
2500 * Get alignment for frame N.
2501 */
2502 static transformation of(int n) {
2503 update_to(n);
2504 if (n == latest)
2505 return latest_t;
2506 else if (_keep)
2507 return kept_t[n];
2508 else {
2509 assert(0);
2510 exit(1);
2511 }
2512 }
2513
2514 /*
2515 * Use Monte Carlo alignment sampling with argument N.
2516 */
2517 static void mc(ale_pos n) {
2518 _mc = n;
2519 }
2520
2521 /*
2522 * Set the certainty-weighted flag.
2523 */
2524 static void certainty_weighted(int flag) {
2525 certainty_weights = flag;
2526 }
2527
2528 /*
2529 * Set the global search type.
2530 */
2531 static void gs(const char *type) {
2532 if (!strcmp(type, "local")) {
2533 _gs = 0;
2534 } else if (!strcmp(type, "inner")) {
2535 _gs = 1;
2536 } else if (!strcmp(type, "outer")) {
2537 _gs = 2;
2538 } else if (!strcmp(type, "all")) {
2539 _gs = 3;
2540 } else if (!strcmp(type, "central")) {
2541 _gs = 4;
2542 } else if (!strcmp(type, "defaults")) {
2543 _gs = 6;
2544 } else if (!strcmp(type, "points")) {
2545 _gs = 5;
2546 keep();
2547 } else {
2548 ui::get()->error("bad global search type");
2549 }
2550 }
2551
2552 /*
2553 * Set the minimum overlap for global searching
2554 */
2555 static void gs_mo(ale_pos value, int _gs_mo_percent) {
2556 _gs_mo = value;
2557 gs_mo_percent = _gs_mo_percent;
2558 }
2559
2560 /*
2561 * Set mutli-alignment certainty lower bound.
2562 */
2563 static void set_ma_cert(ale_real value) {
2564 _ma_cert = value;
2565 }
2566
2567 /*
2568 * Set alignment exclusion regions
2569 */
2570 static void set_exclusion(exclusion *_ax_parameters, int _ax_count) {
2571 ax_count = _ax_count;
2572 ax_parameters = _ax_parameters;
2573 }
2574
2575 /*
2576 * Get match summary statistics.
2577 */
2578 static ale_accum match_summary() {
2579 return match_sum / (ale_accum) match_count;
2580 }
2581 };
2582
2583 template<class diff_trans>
2584 void *d2::align::diff_stat_generic<diff_trans>::diff_subdomain(void *args) {
2585
2586 subdomain_args *sargs = (subdomain_args *) args;
2587
2588 struct scale_cluster c = sargs->c;
2589 std::vector<run> runs = sargs->runs;
2590 int ax_count = sargs->ax_count;
2591 int f = sargs->f;
2592 int i_min = sargs->i_min;
2593 int i_max = sargs->i_max;
2594 int j_min = sargs->j_min;
2595 int j_max = sargs->j_max;
2596 int subdomain = sargs->subdomain;
2597
2598 assert (reference_image);
2599
2600 point offset = c.accum->offset();
2601
2602 for (mc_iterate m(i_min, i_max, j_min, j_max, subdomain); !m.done(); m++) {
2603
2604 int i = m.get_i();
2605 int j = m.get_j();
2606
2607 /*
2608 * Check reference frame definition.
2609 */
2610
2611 if (!((pixel) c.accum->get_pixel(i, j)).finite()
2612 || !(((pixel) c.certainty->get_pixel(i, j)).minabs_norm() > 0))
2613 continue;
2614
2615 /*
2616 * Check for exclusion in render coordinates.
2617 */
2618
2619 if (ref_excluded(i, j, offset, c.ax_parameters, ax_count))
2620 continue;
2621
2622 /*
2623 * Transform
2624 */
2625
2626 struct point q, r = point::undefined();
2627
2628 q = (c.input_scale < 1.0 && interpolant == NULL)
2629 ? runs.back().offset.scaled_inverse_transform(
2630 point(i + offset[0], j + offset[1]))
2631 : runs.back().offset.unscaled_inverse_transform(
2632 point(i + offset[0], j + offset[1]));
2633
2634 if (runs.size() == 2) {
2635 r = (c.input_scale < 1.0)
2636 ? runs.front().offset.scaled_inverse_transform(
2637 point(i + offset[0], j + offset[1]))
2638 : runs.front().offset.unscaled_inverse_transform(
2639 point(i + offset[0], j + offset[1]));
2640 }
2641
2642 ale_pos ti = q[0];
2643 ale_pos tj = q[1];
2644
2645 /*
2646 * Check that the transformed coordinates are within
2647 * the boundaries of array c.input and that they
2648 * are not subject to exclusion.
2649 *
2650 * Also, check that the weight value in the accumulated array
2651 * is nonzero, unless we know it is nonzero by virtue of the
2652 * fact that it falls within the region of the original frame
2653 * (e.g. when we're not increasing image extents).
2654 */
2655
2656 if (input_excluded(ti, tj, c.ax_parameters, ax_count))
2657 continue;
2658
2659 if (input_excluded(r[0], r[1], c.ax_parameters, ax_count))
2660 r = point::undefined();
2661
2662 /*
2663 * Check the boundaries of the input frame
2664 */
2665
2666 if (!(ti >= 0
2667 && ti <= c.input->height() - 1
2668 && tj >= 0
2669 && tj <= c.input->width() - 1))
2670 continue;
2671
2672 if (!(r[0] >= 0
2673 && r[0] <= c.input->height() - 1
2674 && r[1] >= 0
2675 && r[1] <= c.input->width() - 1))
2676 r = point::undefined();
2677
2678 sargs->runs.back().sample(f, c, i, j, q, r, runs.front());
2679 }
2680
2681 return NULL;
2682 }
2683
2684 #endif
synthetic capture engine and renderer