?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
d2::align: Fix active channel determination in scale cluster initialization.
[Ale.git] / d2 / align.h
blob2553b7cb60b60e506205656955e9655af24e5923
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 #define D2_ALIGN_IMPRECISE_COVERAGE 0.3
37
38 class align {
39 private:
40
41 /*
42 * Private data members
43 */
44
45 static ale_pos scale_factor;
46
47 /*
48 * Original frame transformation
49 */
50 static transformation orig_t;
51
52 /*
53 * Keep data older than latest
54 */
55 static int _keep;
56 static transformation *kept_t;
57 static int *kept_ok;
58
59 /*
60 * Transformation file handlers
61 */
62
63 static tload_t *tload;
64 static tsave_t *tsave;
65
66 /*
67 * Control point variables
68 */
69
70 static const point **cp_array;
71 static unsigned int cp_count;
72
73 /*
74 * Reference rendering to align against
75 */
76
77 static render *reference;
78 static filter::scaled_filter *interpolant;
79 static const image *reference_image;
80 static const image *reference_defined;
81
82 /*
83 * Per-pixel alignment weight map
84 */
85
86 static const image *weight_map;
87
88 /*
89 * Frequency-dependent alignment weights
90 */
91
92 static double horiz_freq_cut;
93 static double vert_freq_cut;
94 static double avg_freq_cut;
95 static const char *fw_output;
96
97 /*
98 * Algorithmic alignment weighting
99 */
100
101 static const char *wmx_exec;
102 static const char *wmx_file;
103 static const char *wmx_defs;
104
105 /*
106 * Non-certainty alignment weights
107 */
108
109 static image *alignment_weights;
110
111 /*
112 * Latest transformation.
113 */
114
115 static transformation latest_t;
116
117 /*
118 * Flag indicating whether the latest transformation
119 * resulted in a match.
120 */
121
122 static int latest_ok;
123
124 /*
125 * Frame number most recently aligned.
126 */
127
128 static int latest;
129
130 /*
131 * Exposure registration
132 *
133 * 0. Preserve the original exposure of images.
134 *
135 * 1. Match exposure between images.
136 *
137 * 2. Use only image metadata for registering exposure.
138 */
139
140 static int _exp_register;
141
142 /*
143 * Alignment class.
144 *
145 * 0. Translation only. Only adjust the x and y position of images.
146 * Do not rotate input images or perform projective transformations.
147 *
148 * 1. Euclidean transformations only. Adjust the x and y position
149 * of images and the orientation of the image about the image center.
150 *
151 * 2. Perform general projective transformations. See the file gpt.h
152 * for more information about general projective transformations.
153 */
154
155 static int alignment_class;
156
157 /*
158 * Default initial alignment type.
159 *
160 * 0. Identity transformation.
161 *
162 * 1. Most recently accepted frame's final transformation.
163 */
164
165 static int default_initial_alignment_type;
166
167 /*
168 * Projective group behavior
169 *
170 * 0. Perturb in output coordinates.
171 *
172 * 1. Perturb in source coordinates
173 */
174
175 static int perturb_type;
176
177 /*
178 * Alignment element
179 *
180 * This structure contains variables necessary for handling a
181 * multi-alignment element. The change between the non-default old
182 * initial alignment and old final alignment is used to adjust the
183 * non-default current initial alignment. If either the old or new
184 * initial alignment is a default alignment, the old --follow semantics
185 * are preserved.
186 */
187
188 struct element_t {
189 int is_default, old_is_default;
190 int is_primary;
191 int old_lod;
192 transformation old_initial_alignment;
193 transformation old_final_alignment;
194 transformation default_initial_alignment;
195 const image *input_frame;
196
197 public:
198 element_t() {
199 is_default = 1;
200 input_frame = NULL;
201 }
202 };
203
204 /*
205 * Alignment for failed frames -- default or optimal?
206 *
207 * A frame that does not meet the match threshold can be assigned the
208 * best alignment found, or can be assigned its alignment default.
209 */
210
211 static int is_fail_default;
212
213 /*
214 * Alignment code.
215 *
216 * 0. Align images with an error contribution from each color channel.
217 *
218 * 1. Align images with an error contribution only from the green channel.
219 * Other color channels do not affect alignment.
220 *
221 * 2. Align images using a summation of channels. May be useful when dealing
222 * with images that have high frequency color ripples due to color aliasing.
223 */
224
225 static int channel_alignment_type;
226
227 /*
228 * Error metric exponent
229 */
230
231 static float metric_exponent;
232
233 /*
234 * Match threshold
235 */
236
237 static float match_threshold;
238
239 /*
240 * Perturbation lower and upper bounds.
241 */
242
243 static ale_pos perturb_lower;
244 static int perturb_lower_percent;
245 static ale_pos perturb_upper;
246 static int perturb_upper_percent;
247
248 /*
249 * Maximum level-of-detail scale factor is 2^lod_max/perturb.
250 */
251
252 static int lod_max;
253
254 /*
255 * Maximum rotational perturbation
256 */
257
258 static ale_pos rot_max;
259
260 /*
261 * Barrel distortion alignment multiplier
262 */
263
264 static ale_pos bda_mult;
265
266 /*
267 * Barrel distortion maximum adjustment rate
268 */
269
270 static ale_pos bda_rate;
271
272 /*
273 * Alignment match sum
274 */
275
276 static ale_accum match_sum;
277
278 /*
279 * Alignment match count.
280 */
281
282 static int match_count;
283
284 /*
285 * Exposure registration and final match precision flag
286 *
287 * 0. Do not require precise calculation.
288 *
289 * 1. Require precise exposure registration and final match.
290 */
291
292 static int precise_calculation;
293
294 /*
295 * Certainty weight flag
296 *
297 * 0. Don't use certainty weights for alignment.
298 *
299 * 1. Use certainty weights for alignment.
300 */
301
302 static int certainty_weights;
303
304 /*
305 * Global search parameter
306 *
307 * 0. Local: Local search only.
308 * 1. Inner: Alignment reference image inner region
309 * 2. Outer: Alignment reference image outer region
310 * 3. All: Alignment reference image inner and outer regions.
311 * 4. Central: Inner if possible; else, best of inner and outer.
312 * 5. Points: Align by control points.
313 */
314
315 static int _gs;
316
317 /*
318 * Multi-alignment cardinality.
319 */
320
321 static unsigned int _ma_card;
322
323 /*
324 * Multi-alignment contiguity.
325 */
326
327 static double _ma_cont;
328
329 /*
330 * Minimum overlap for global searches
331 */
332
333 static ale_pos _gs_mo;
334 static int gs_mo_percent;
335
336 /*
337 * Exclusion regions
338 */
339
340 static exclusion *ax_parameters;
341 static int ax_count;
342
343 /*
344 * Types for scale clusters.
345 */
346
347 struct nl_scale_cluster {
348 const image *accum_max;
349 const image *accum_min;
350 const image *certainty_max;
351 const image *certainty_min;
352 const image *aweight_max;
353 const image *aweight_min;
354 exclusion *ax_parameters;
355
356 ale_pos input_scale;
357 const image *input_certainty_max;
358 const image *input_certainty_min;
359 const image *input_max;
360 const image *input_min;
361 };
362
363 struct scale_cluster {
364 const image *accum;
365 const image *certainty;
366 const image *aweight;
367 exclusion *ax_parameters;
368
369 ale_pos input_scale;
370 const image *input_certainty;
371 const image *input;
372
373 nl_scale_cluster *nl_scale_clusters;
374 };
375
376 /*
377 * Check for exclusion region coverage in the reference
378 * array.
379 */
380 static int ref_excluded(int i, int j, point offset, exclusion *params, int param_count) {
381 for (int idx = 0; idx < param_count; idx++)
382 if (params[idx].type == exclusion::RENDER
383 && i + offset[0] >= params[idx].x[0]
384 && i + offset[0] <= params[idx].x[1]
385 && j + offset[1] >= params[idx].x[2]
386 && j + offset[1] <= params[idx].x[3])
387 return 1;
388
389 return 0;
390 }
391
392 /*
393 * Check for exclusion region coverage in the input
394 * array.
395 */
396 static int input_excluded(ale_pos ti, ale_pos tj, exclusion *params, int param_count) {
397 for (int idx = 0; idx < param_count; idx++)
398 if (params[idx].type == exclusion::FRAME
399 && ti >= params[idx].x[0]
400 && ti <= params[idx].x[1]
401 && tj >= params[idx].x[2]
402 && tj <= params[idx].x[3])
403 return 1;
404
405 return 0;
406 }
407
408 /*
409 * Overlap function. Determines the number of pixels in areas where
410 * the arrays overlap. Uses the reference array's notion of pixel
411 * positions.
412 */
413 static unsigned int overlap(struct scale_cluster c, transformation t, int ax_count) {
414 assert (reference_image);
415
416 unsigned int result = 0;
417
418 point offset = c.accum->offset();
419
420 for (unsigned int i = 0; i < c.accum->height(); i++)
421 for (unsigned int j = 0; j < c.accum->width(); j++) {
422
423 if (ref_excluded(i, j, offset, c.ax_parameters, ax_count))
424 continue;
425
426 /*
427 * Transform
428 */
429
430 struct point q;
431
432 q = (c.input_scale < 1.0 && interpolant == NULL)
433 ? t.scaled_inverse_transform(
434 point(i + offset[0], j + offset[1]))
435 : t.unscaled_inverse_transform(
436 point(i + offset[0], j + offset[1]));
437
438 ale_pos ti = q[0];
439 ale_pos tj = q[1];
440
441 /*
442 * Check that the transformed coordinates are within
443 * the boundaries of array c.input, and check that the
444 * weight value in the accumulated array is nonzero,
445 * unless we know it is nonzero by virtue of the fact
446 * that it falls within the region of the original
447 * frame (e.g. when we're not increasing image
448 * extents). Also check for frame exclusion.
449 */
450
451 if (input_excluded(ti, tj, c.ax_parameters, ax_count))
452 continue;
453
454 if (ti >= 0
455 && ti <= c.input->height() - 1
456 && tj >= 0
457 && tj <= c.input->width() - 1
458 && c.certainty->get_pixel(i, j)[0] != 0)
459 result++;
460 }
461
462 return result;
463 }
464
465 /*
466 * Calculate the region associated with the current multi-alignment
467 * element.
468 */
469 static void calculate_element_region(transformation *t, scale_cluster si,
470 int local_ax_count) {
471
472 unsigned int i_max = si.accum->height();
473 unsigned int j_max = si.accum->width();
474 point offset = si.accum->offset();
475
476 if (si.input_scale < 1.0 && interpolant == NULL)
477 t->begin_calculate_scaled_region(i_max, j_max, offset);
478 else
479 t->begin_calculate_unscaled_region(i_max, j_max, offset);
480
481 for (unsigned int i = 0; i < i_max; i++)
482 for (unsigned int j = 0; j < j_max; j++) {
483
484 if (ref_excluded(i, j, offset, si.ax_parameters, local_ax_count))
485 continue;
486
487 point q;
488
489 while ((q = t->get_query_point((int) (i + offset[0]),
490 (int) (j + offset[1]))).defined()) {
491
492 ale_pos ti = q[0];
493 ale_pos tj = q[1];
494
495 if (input_excluded(ti, tj, si.ax_parameters, ax_count))
496 continue;
497
498 if (ti >= 0
499 && ti <= si.input->height() - 1
500 && tj >= 0
501 && tj <= si.input->width() - 1
502 && si.certainty->get_pixel(i, j)[0] != 0) {
503
504 assert(0);
505 }
506 }
507 }
508
509 t->end_calculate_region();
510 }
511
512 /*
513 * Monte carlo iteration class.
514 *
515 * Monte Carlo alignment has been used for statistical comparisons in
516 * spatial registration, and is now used for tonal registration
517 * and final match calculation.
518 */
519
520 /*
521 * We use a random process for which the expected number of sampled
522 * pixels is +/- .000003 from the coverage in the range [.005,.995] for
523 * an image with 100,000 pixels. (The actual number may still deviate
524 * from the expected number by more than this amount, however.) The
525 * method is as follows:
526 *
527 * We have coverage == USE/ALL, or (expected # pixels to use)/(# total
528 * pixels). We derive from this SKIP/USE.
529 *
530 * SKIP/USE == (SKIP/ALL)/(USE/ALL) == (1 - (USE/ALL))/(USE/ALL)
531 *
532 * Once we have SKIP/USE, we know the expected number of pixels to skip
533 * in each iteration. We use a random selection process that provides
534 * SKIP/USE close to this calculated value.
535 *
536 * If we can draw uniformly to select the number of pixels to skip, we
537 * do. In this case, the maximum number of pixels to skip is twice the
538 * expected number.
539 *
540 * If we cannot draw uniformly, we still assign equal probability to
541 * each of the integer values in the interval [0, 2 * (SKIP/USE)], but
542 * assign an unequal amount to the integer value ceil(2 * SKIP/USE) +
543 * 1.
544 */
545
546 /*
547 * When reseeding the random number generator, we want the same set of
548 * pixels to be used in cases where two alignment options are compared.
549 * If we wanted to avoid bias from repeatedly utilizing the same seed,
550 * we could seed with the number of the frame most recently aligned:
551 *
552 * srand(latest);
553 *
554 * However, in cursory tests, it seems okay to just use the default
555 * seed of 1, and so we do this, since it is simpler; both of these
556 * approaches to reseeding achieve better results than not reseeding.
557 * (1 is the default seed according to the GNU Manual Page for
558 * rand(3).)
559 *
560 * For subdomain calculations, we vary the seed by adding the subdomain
561 * index.
562 */
563
564 class mc_iterate {
565 ale_pos mc_max;
566 unsigned int index;
567 unsigned int index_max;
568 int i_min;
569 int i_max;
570 int j_min;
571 int j_max;
572
573 rng_t rng;
574
575 public:
576 mc_iterate(int _i_min, int _i_max, int _j_min, int _j_max, unsigned int subdomain)
577 : rng() {
578
579 ale_pos coverage;
580
581 i_min = _i_min;
582 i_max = _i_max;
583 j_min = _j_min;
584 j_max = _j_max;
585
586 index_max = (i_max - i_min) * (j_max - j_min);
587
588 if (index_max < 500 || precise_calculation)
589 coverage = 1;
590 else
591 coverage = D2_ALIGN_IMPRECISE_COVERAGE;
592
593 ale_pos su = (1 - coverage) / coverage;
594
595 mc_max = (floor(2*su) * (1 + floor(2*su)) + 2*su)
596 / (2 + 2 * floor(2*su) - 2*su);
597
598 rng.seed(1 + subdomain);
599
600 index = -1 + (int) ceil((mc_max+1)
601 * ( (1 + ((ale_pos) (rng.get())) )
602 / (1 + ((ale_pos) RAND_MAX)) ));
603 }
604
605 int get_i() {
606 return index / (j_max - j_min) + i_min;
607 }
608
609 int get_j() {
610 return index % (j_max - j_min) + j_min;
611 }
612
613 void operator++(int whats_this_for) {
614 index += (int) ceil((mc_max+1)
615 * ( (1 + ((ale_pos) (rng.get())) )
616 / (1 + ((ale_pos) RAND_MAX)) ));
617 }
618
619 int done() {
620 return (index >= index_max);
621 }
622 };
623
624 /*
625 * Not-quite-symmetric difference function. Determines the difference in areas
626 * where the arrays overlap. Uses the reference array's notion of pixel positions.
627 *
628 * For the purposes of determining the difference, this function divides each
629 * pixel value by the corresponding image's average pixel magnitude, unless we
630 * are:
631 *
632 * a) Extending the boundaries of the image, or
633 *
634 * b) following the previous frame's transform
635 *
636 * If we are doing monte-carlo pixel sampling for alignment, we
637 * typically sample a subset of available pixels; otherwise, we sample
638 * all pixels.
639 *
640 */
641
642 class diff_stat_t {
643
644 struct run {
645
646 transformation offset;
647 ale_pos perturb;
648
649 ale_accum result;
650 ale_accum divisor;
651
652 point max, min;
653 ale_accum centroid[2], centroid_divisor;
654 ale_accum de_centroid[2], de_centroid_v, de_sum;
655
656 void init() {
657
658 result = 0;
659 divisor = 0;
660
661 min = point::posinf();
662 max = point::neginf();
663
664 centroid[0] = 0;
665 centroid[1] = 0;
666 centroid_divisor = 0;
667
668 de_centroid[0] = 0;
669 de_centroid[1] = 0;
670
671 de_centroid_v = 0;
672
673 de_sum = 0;
674 }
675
676 void init(transformation _offset, ale_pos _perturb) {
677 offset = _offset;
678 perturb = _perturb;
679 init();
680 }
681
682 /*
683 * Required for STL sanity.
684 */
685 run() : offset() {
686 init();
687 }
688
689 run(transformation _offset, ale_pos _perturb) : offset() {
690 init(_offset, _perturb);
691 }
692
693 void add(const run &_run) {
694 result += _run.result;
695 divisor += _run.divisor;
696
697 for (int d = 0; d < 2; d++) {
698 if (min[d] > _run.min[d])
699 min[d] = _run.min[d];
700 if (max[d] < _run.max[d])
701 max[d] = _run.max[d];
702 centroid[d] += _run.centroid[d];
703 de_centroid[d] += _run.de_centroid[d];
704 }
705
706 centroid_divisor += _run.centroid_divisor;
707 de_centroid_v += _run.de_centroid_v;
708 de_sum += _run.de_sum;
709 }
710
711 run(const run &_run) : offset() {
712
713 /*
714 * Initialize
715 */
716 init(_run.offset, _run.perturb);
717
718 /*
719 * Add
720 */
721 add(_run);
722 }
723
724 run &operator=(const run &_run) {
725
726 /*
727 * Initialize
728 */
729 init(_run.offset, _run.perturb);
730
731 /*
732 * Add
733 */
734 add(_run);
735
736 return *this;
737 }
738
739 ~run() {
740 }
741
742 ale_accum get_error() const {
743 return pow(result / divisor, 1/metric_exponent);
744 }
745
746 void sample(int f, scale_cluster c, int i, int j, point t, point u,
747 const run &comparison) {
748
749 pixel pa = c.accum->get_pixel(i, j);
750
751 ale_accum this_result[2] = { 0, 0 };
752 ale_accum this_divisor[2] = { 0, 0 };
753
754 pixel p[2];
755 pixel weight[2];
756 weight[0] = pixel(1, 1, 1);
757 weight[1] = pixel(1, 1, 1);
758
759 if (interpolant != NULL) {
760 interpolant->filtered(i, j, &p[0], &weight[1], 1, f);
761 } else {
762 p[0] = c.input->get_bl(t);
763 }
764
765 if (u.defined()) {
766 p[1] = c.input->get_bl(u);
767 }
768
769
770 /*
771 * Handle certainty.
772 */
773
774 if (certainty_weights == 1) {
775
776 /*
777 * For speed, use arithmetic interpolation (get_bl(.))
778 * instead of geometric (get_bl(., 1))
779 */
780
781 weight[0] *= c.input_certainty->get_bl(t);
782 if (u.defined())
783 weight[1] *= c.input_certainty->get_bl(u);
784 weight[0] *= c.certainty->get_pixel(i, j);
785 weight[1] *= c.certainty->get_pixel(i, j);
786 }
787
788 if (c.aweight != NULL) {
789 weight[0] *= c.aweight->get_pixel(i, j);
790 weight[1] *= c.aweight->get_pixel(i, j);
791 }
792
793 /*
794 * Update sampling area statistics
795 */
796
797 if (min[0] > i)
798 min[0] = i;
799 if (min[1] > j)
800 min[1] = j;
801 if (max[0] < i)
802 max[0] = i;
803 if (max[1] < j)
804 max[1] = j;
805
806 centroid[0] += (weight[0][0] + weight[0][1] + weight[0][2]) * i;
807 centroid[1] += (weight[0][0] + weight[0][1] + weight[0][2]) * j;
808 centroid_divisor += (weight[0][0] + weight[0][1] + weight[0][2]);
809
810 /*
811 * Determine alignment type.
812 */
813
814 for (int m = 0; m < (u.defined() ? 2 : 1); m++)
815 if (channel_alignment_type == 0) {
816 /*
817 * Align based on all channels.
818 */
819
820
821 for (int k = 0; k < 3; k++) {
822 ale_real achan = pa[k];
823 ale_real bchan = p[m][k];
824
825 this_result[m] += weight[m][k] * pow(fabs(achan - bchan), metric_exponent);
826 this_divisor[m] += weight[m][k] * pow(achan > bchan ? achan : bchan, metric_exponent);
827 }
828 } else if (channel_alignment_type == 1) {
829 /*
830 * Align based on the green channel.
831 */
832
833 ale_real achan = pa[1];
834 ale_real bchan = p[m][1];
835
836 this_result[m] = weight[m][1] * pow(fabs(achan - bchan), metric_exponent);
837 this_divisor[m] = weight[m][1] * pow(achan > bchan ? achan : bchan, metric_exponent);
838 } else if (channel_alignment_type == 2) {
839 /*
840 * Align based on the sum of all channels.
841 */
842
843 ale_real asum = 0;
844 ale_real bsum = 0;
845 ale_real wsum = 0;
846
847 for (int k = 0; k < 3; k++) {
848 asum += pa[k];
849 bsum += p[m][k];
850 wsum += weight[m][k] / 3;
851 }
852
853 this_result[m] = wsum * pow(fabs(asum - bsum), metric_exponent);
854 this_divisor[m] = wsum * pow(asum > bsum ? asum : bsum, metric_exponent);
855 }
856
857 if (u.defined()) {
858 // ale_accum de = fabs(this_result[0] / this_divisor[0]
859 // - this_result[1] / this_divisor[1]);
860 ale_accum de = fabs(this_result[0] - this_result[1]);
861
862 de_centroid[0] += de * i;
863 de_centroid[1] += de * j;
864
865 de_centroid_v += de * t.lengthto(u);
866
867 de_sum += de;
868 }
869
870 result += (this_result[0]);
871 divisor += (this_divisor[0]);
872 }
873
874 void rescale(ale_pos scale) {
875 offset.rescale(scale);
876
877 de_centroid[0] *= scale;
878 de_centroid[1] *= scale;
879 de_centroid_v *= scale;
880 }
881
882 point get_centroid() {
883 point result = point(centroid[0] / centroid_divisor, centroid[1] / centroid_divisor);
884
885 assert (finite(centroid[0])
886 && finite(centroid[1])
887 && (result.defined() || centroid_divisor == 0));
888
889 return result;
890 }
891
892 point get_error_centroid() {
893 point result = point(de_centroid[0] / de_sum, de_centroid[1] / de_sum);
894 return result;
895 }
896
897
898 ale_pos get_error_perturb() {
899 ale_pos result = de_centroid_v / de_sum;
900
901 return result;
902 }
903
904 };
905
906 /*
907 * When non-empty, runs.front() is best, runs.back() is
908 * testing.
909 */
910
911 std::vector<run> runs;
912
913 /*
914 * old_runs stores the latest available perturbation set for
915 * each multi-alignment element.
916 */
917
918 typedef std::pair<unsigned int, unsigned int> run_index;
919 std::map<run_index, run> old_runs;
920
921 static void *diff_subdomain(void *args);
922
923 struct subdomain_args {
924 struct scale_cluster c;
925 std::vector<run> runs;
926 int ax_count;
927 int f;
928 int i_min, i_max, j_min, j_max;
929 int subdomain;
930 };
931
932 int get_current_index() const {
933 assert(runs.size());
934 return runs[0].offset.get_current_index();
935 }
936
937 struct scale_cluster si;
938 int ax_count;
939 int frame;
940
941 std::vector<ale_pos> perturb_multipliers;
942
943 public:
944 void diff(struct scale_cluster c, ale_pos perturb,
945 transformation t,
946 int _ax_count, int f) {
947
948 if (runs.size() == 2)
949 runs.pop_back();
950
951 runs.push_back(run(t, perturb));
952
953 si = c;
954 ax_count = _ax_count;
955 frame = f;
956
957 ui::get()->d2_align_sample_start();
958
959 if (interpolant != NULL)
960 interpolant->set_parameters(t, c.input, c.accum->offset());
961
962 int N;
963 #ifdef USE_PTHREAD
964 N = thread::count();
965
966 pthread_t *threads = (pthread_t *) malloc(sizeof(pthread_t) * N);
967 pthread_attr_t *thread_attr = (pthread_attr_t *) malloc(sizeof(pthread_attr_t) * N);
968
969 #else
970 N = 1;
971 #endif
972
973 subdomain_args *args = new subdomain_args[N];
974
975 for (int ti = 0; ti < N; ti++) {
976 args[ti].c = c;
977 args[ti].runs = runs;
978 args[ti].ax_count = ax_count;
979 args[ti].f = f;
980 args[ti].i_min = (c.accum->height() * ti) / N;
981 args[ti].i_max = (c.accum->height() * (ti + 1)) / N;
982 args[ti].j_min = 0;
983 args[ti].j_max = c.accum->width();
984 args[ti].subdomain = ti;
985 #ifdef USE_PTHREAD
986 pthread_attr_init(&thread_attr[ti]);
987 pthread_attr_setdetachstate(&thread_attr[ti], PTHREAD_CREATE_JOINABLE);
988 pthread_create(&threads[ti], &thread_attr[ti], diff_subdomain, &args[ti]);
989 #else
990 diff_subdomain(&args[ti]);
991 #endif
992 }
993
994 for (int ti = 0; ti < N; ti++) {
995 #ifdef USE_PTHREAD
996 pthread_join(threads[ti], NULL);
997 #endif
998 runs.back().add(args[ti].runs.back());
999 }
1000
1001 delete[] args;
1002
1003 ui::get()->d2_align_sample_stop();
1004
1005 }
1006
1007 private:
1008 void rediff() {
1009 std::vector<transformation> t_array;
1010 std::vector<ale_pos> p_array;
1011
1012 for (unsigned int r = 0; r < runs.size(); r++) {
1013 t_array.push_back(runs[r].offset);
1014 p_array.push_back(runs[r].perturb);
1015 }
1016
1017 runs.clear();
1018
1019 for (unsigned int r = 0; r < t_array.size(); r++)
1020 diff(si, p_array[r], t_array[r], ax_count, frame);
1021 }
1022
1023
1024 public:
1025 int better() {
1026 assert(runs.size() >= 2);
1027 assert(runs[0].offset.scale() == runs[1].offset.scale());
1028
1029 return (runs[1].get_error() < runs[0].get_error()
1030 || (!finite(runs[0].get_error()) && finite(runs[1].get_error())));
1031 }
1032
1033 diff_stat_t() : runs(), old_runs(), perturb_multipliers() {
1034 }
1035
1036 run_index get_run_index(unsigned int perturb_index) {
1037 return run_index(get_current_index(), perturb_index);
1038 }
1039
1040 run &get_run(unsigned int perturb_index) {
1041 run_index index = get_run_index(perturb_index);
1042
1043 assert(old_runs.count(index));
1044 return old_runs[index];
1045 }
1046
1047 void rescale(ale_pos scale, scale_cluster _si) {
1048 assert(runs.size() == 1);
1049
1050 si = _si;
1051
1052 runs[0].rescale(scale);
1053
1054 rediff();
1055 }
1056
1057 void push_element() {
1058 assert(runs.size() == 1);
1059
1060 runs[0].offset.push_element();
1061
1062 rediff();
1063 }
1064
1065 unsigned int get_current_index() {
1066 assert (runs.size() > 0);
1067
1068 return runs[0].offset.get_current_index();
1069 }
1070
1071 void set_current_index(unsigned int i) {
1072 assert(runs.size() == 1);
1073 runs[0].offset.set_current_index(i);
1074 rediff();
1075 }
1076
1077 void calculate_element_region() {
1078 assert(runs.size() == 1);
1079
1080 if (get_offset().get_current_index() > 0
1081 && get_offset().is_nontrivial())
1082 align::calculate_element_region(&runs[0].offset, si, ax_count);
1083 }
1084
1085 ~diff_stat_t() {
1086 }
1087
1088 diff_stat_t &operator=(const diff_stat_t &dst) {
1089 /*
1090 * Copy run information.
1091 */
1092 runs = dst.runs;
1093 old_runs = dst.old_runs;
1094
1095 /*
1096 * Copy diff variables
1097 */
1098 si = dst.si;
1099 ax_count = dst.ax_count;
1100 frame = dst.frame;
1101 perturb_multipliers = dst.perturb_multipliers;
1102
1103 return *this;
1104 }
1105
1106 diff_stat_t(const diff_stat_t &dst) : runs(), old_runs(),
1107 perturb_multipliers() {
1108 operator=(dst);
1109 }
1110
1111 ale_accum get_result() {
1112 assert(runs.size() == 1);
1113 return runs[0].result;
1114 }
1115
1116 ale_accum get_divisor() {
1117 assert(runs.size() == 1);
1118 return runs[0].divisor;
1119 }
1120
1121 transformation get_offset() {
1122 assert(runs.size() == 1);
1123 return runs[0].offset;
1124 }
1125
1126 int operator!=(diff_stat_t &param) {
1127 return (get_error() != param.get_error());
1128 }
1129
1130 int operator==(diff_stat_t &param) {
1131 return !(operator!=(param));
1132 }
1133
1134 ale_pos get_error_perturb() {
1135 assert(runs.size() == 1);
1136 return runs[0].get_error_perturb();
1137 }
1138
1139 ale_accum get_error() const {
1140 assert(runs.size() == 1);
1141 return runs[0].get_error();
1142 }
1143
1144 public:
1145 /*
1146 * Get the set of transformations produced by a given perturbation
1147 */
1148 void get_perturb_set(std::vector<transformation> *set,
1149 ale_pos adj_p, ale_pos adj_o, ale_pos adj_b,
1150 ale_pos *current_bd, ale_pos *modified_bd,
1151 std::vector<ale_pos> multipliers = std::vector<ale_pos>()) {
1152
1153 assert(runs.size() == 1);
1154
1155 transformation test_t;
1156
1157 /*
1158 * Translational or euclidean transformation
1159 */
1160
1161 for (unsigned int i = 0; i < 2; i++)
1162 for (unsigned int s = 0; s < 2; s++) {
1163
1164 if (!multipliers.size())
1165 multipliers.push_back(1);
1166
1167 assert(finite(multipliers[0]));
1168
1169 test_t = get_offset();
1170
1171 // test_t.eu_modify(i, (s ? -adj_p : adj_p) * multipliers[0]);
1172 test_t.translate((i ? point(1, 0) : point(0, 1))
1173 * (s ? -adj_p : adj_p)
1174 * multipliers[0]);
1175
1176 test_t.snap(adj_p / 2);
1177
1178 set->push_back(test_t);
1179 multipliers.erase(multipliers.begin());
1180
1181 }
1182
1183 if (alignment_class > 0)
1184 for (unsigned int s = 0; s < 2; s++) {
1185
1186 if (!multipliers.size())
1187 multipliers.push_back(1);
1188
1189 assert(multipliers.size());
1190 assert(finite(multipliers[0]));
1191
1192 if (!(adj_o * multipliers[0] < rot_max))
1193 return;
1194
1195 ale_pos adj_s = (s ? 1 : -1) * adj_o * multipliers[0];
1196
1197 test_t = get_offset();
1198
1199 run_index ori = get_run_index(set->size());
1200 point centroid = point::undefined();
1201
1202 if (!old_runs.count(ori))
1203 ori = get_run_index(0);
1204
1205 if (!centroid.finite() && old_runs.count(ori)) {
1206 centroid = old_runs[ori].get_error_centroid();
1207
1208 if (!centroid.finite())
1209 centroid = old_runs[ori].get_centroid();
1210
1211 centroid *= test_t.scale()
1212 / old_runs[ori].offset.scale();
1213 }
1214
1215 if (!centroid.finite() && !test_t.is_projective()) {
1216 test_t.eu_modify(2, adj_s);
1217 } else if (!centroid.finite()) {
1218 centroid = point(si.input->height() / 2,
1219 si.input->width() / 2);
1220
1221 test_t.rotate(centroid + si.accum->offset(),
1222 adj_s);
1223 } else {
1224 test_t.rotate(centroid + si.accum->offset(),
1225 adj_s);
1226 }
1227
1228 test_t.snap(adj_p / 2);
1229
1230 set->push_back(test_t);
1231 multipliers.erase(multipliers.begin());
1232 }
1233
1234 if (alignment_class == 2) {
1235
1236 /*
1237 * Projective transformation
1238 */
1239
1240 for (unsigned int i = 0; i < 4; i++)
1241 for (unsigned int j = 0; j < 2; j++)
1242 for (unsigned int s = 0; s < 2; s++) {
1243
1244 if (!multipliers.size())
1245 multipliers.push_back(1);
1246
1247 assert(multipliers.size());
1248 assert(finite(multipliers[0]));
1249
1250 ale_pos adj_s = (s ? -1 : 1) * adj_p * multipliers [0];
1251
1252 test_t = get_offset();
1253
1254 if (perturb_type == 0)
1255 test_t.gpt_modify(j, i, adj_s);
1256 else if (perturb_type == 1)
1257 test_t.gr_modify(j, i, adj_s);
1258 else
1259 assert(0);
1260
1261 test_t.snap(adj_p / 2);
1262
1263 set->push_back(test_t);
1264 multipliers.erase(multipliers.begin());
1265 }
1266
1267 }
1268
1269 /*
1270 * Barrel distortion
1271 */
1272
1273 if (bda_mult != 0 && adj_b != 0) {
1274
1275 for (unsigned int d = 0; d < get_offset().bd_count(); d++)
1276 for (unsigned int s = 0; s < 2; s++) {
1277
1278 if (!multipliers.size())
1279 multipliers.push_back(1);
1280
1281 assert (multipliers.size());
1282 assert (finite(multipliers[0]));
1283
1284 ale_pos adj_s = (s ? -1 : 1) * adj_b * multipliers[0];
1285
1286 if (bda_rate > 0 && fabs(modified_bd[d] + adj_s - current_bd[d]) > bda_rate)
1287 continue;
1288
1289 transformation test_t = get_offset();
1290
1291 test_t.bd_modify(d, adj_s);
1292
1293 set->push_back(test_t);
1294 }
1295 }
1296 }
1297
1298 void confirm() {
1299 assert(runs.size() == 2);
1300 runs[0] = runs[1];
1301 runs.pop_back();
1302 }
1303
1304 void discard() {
1305 assert(runs.size() == 2);
1306 runs.pop_back();
1307 }
1308
1309 void perturb_test(ale_pos perturb, ale_pos adj_p, ale_pos adj_o, ale_pos adj_b,
1310 ale_pos *current_bd, ale_pos *modified_bd, int stable) {
1311
1312 assert(runs.size() == 1);
1313
1314 std::vector<transformation> t_set;
1315
1316 if (perturb_multipliers.size() == 0) {
1317 get_perturb_set(&t_set, adj_p, adj_o, adj_b,
1318 current_bd, modified_bd);
1319
1320 for (unsigned int i = 0; i < t_set.size(); i++) {
1321 diff_stat_t test = *this;
1322
1323 test.diff(si, perturb, t_set[i], ax_count, frame);
1324
1325 test.confirm();
1326
1327 if (finite(adj_p / test.get_error_perturb()))
1328 perturb_multipliers.push_back(adj_p / test.get_error_perturb());
1329 else
1330 perturb_multipliers.push_back(1);
1331
1332 }
1333
1334 t_set.clear();
1335 }
1336
1337 get_perturb_set(&t_set, adj_p, adj_o, adj_b, current_bd, modified_bd,
1338 perturb_multipliers);
1339
1340 int found_unreliable = 1;
1341 std::vector<int> tested(t_set.size(), 0);
1342
1343 for (unsigned int i = 0; i < t_set.size(); i++) {
1344 run_index ori = get_run_index(i);
1345
1346 /*
1347 * Check for stability
1348 */
1349 if (stable
1350 && old_runs.count(ori)
1351 && old_runs[ori].offset == t_set[i])
1352 tested[i] = 1;
1353 }
1354
1355 std::vector<ale_pos> perturb_multipliers_original = perturb_multipliers;
1356
1357 while (found_unreliable) {
1358
1359 found_unreliable = 0;
1360
1361 for (unsigned int i = 0; i < t_set.size(); i++) {
1362
1363 if (tested[i])
1364 continue;
1365
1366 diff(si, perturb, t_set[i], ax_count, frame);
1367
1368 if (!(i < perturb_multipliers.size())
1369 || !finite(perturb_multipliers[i])) {
1370
1371 perturb_multipliers.resize(i + 1);
1372
1373 perturb_multipliers[i] =
1374 adj_p / runs[1].get_error_perturb();
1375
1376 if (finite(perturb_multipliers[i]))
1377 found_unreliable = 1;
1378
1379 continue;
1380 }
1381
1382 run_index ori = get_run_index(i);
1383
1384 if (old_runs.count(ori) == 0)
1385 old_runs.insert(std::pair<run_index, run>(ori, runs[1]));
1386 else
1387 old_runs[ori] = runs[1];
1388
1389 perturb_multipliers[i] = perturb_multipliers_original[i]
1390 * adj_p / runs[1].get_error_perturb();
1391
1392 if (!finite(perturb_multipliers[i]))
1393 perturb_multipliers[i] = 1;
1394
1395 tested[i] = 1;
1396
1397 if (better()
1398 && runs[1].get_error() < runs[0].get_error()
1399 && perturb_multipliers[i]
1400 / perturb_multipliers_original[i] < 2) {
1401 runs[0] = runs[1];
1402 runs.pop_back();
1403 return;
1404 }
1405
1406 }
1407
1408 if (runs.size() > 1)
1409 runs.pop_back();
1410
1411 if (!found_unreliable)
1412 return;
1413 }
1414 }
1415
1416 /*
1417 * Attempt to make the current element non-trivial, by finding a nearby
1418 * alignment admitting a non-empty element region.
1419 */
1420 void make_element_nontrivial(ale_pos adj_p, ale_pos adj_o) {
1421 assert(runs.size() == 1);
1422
1423 transformation *t = &runs[0].offset;
1424
1425 if (t->is_nontrivial())
1426 return;
1427
1428 calculate_element_region();
1429
1430 if (t->is_nontrivial())
1431 return;
1432
1433 std::vector<transformation> t_set;
1434 get_perturb_set(&t_set, adj_p, adj_o, 0, NULL, NULL);
1435
1436 for (unsigned int i = 0; i < t_set.size(); i++) {
1437
1438 align::calculate_element_region(&t_set[i], si, ax_count);
1439
1440 if (t_set[i].is_nontrivial()) {
1441 *t = t_set[i];
1442 return;
1443 }
1444 }
1445 }
1446
1447 };
1448
1449
1450 /*
1451 * Adjust exposure for an aligned frame B against reference A.
1452 *
1453 * Expects full-LOD images.
1454 *
1455 * Note: This method does not use any weighting, by certainty or
1456 * otherwise, in the first exposure registration pass, as any bias of
1457 * weighting according to color may also bias the exposure registration
1458 * result; it does use weighting, including weighting by certainty
1459 * (even if certainty weighting is not specified), in the second pass,
1460 * under the assumption that weighting by certainty improves handling
1461 * of out-of-range highlights, and that bias of exposure measurements
1462 * according to color may generally be less harmful after spatial
1463 * registration has been performed.
1464 */
1465 class exposure_ratio_iterate : public thread::decompose_domain {
1466 pixel_accum *asums;
1467 pixel_accum *bsums;
1468 pixel_accum *asum;
1469 pixel_accum *bsum;
1470 struct scale_cluster c;
1471 transformation t;
1472 int ax_count;
1473 int pass_number;
1474 protected:
1475 void prepare_subdomains(unsigned int N) {
1476 asums = new pixel_accum[N];
1477 bsums = new pixel_accum[N];
1478 }
1479 void subdomain_algorithm(unsigned int thread,
1480 int i_min, int i_max, int j_min, int j_max) {
1481
1482 point offset = c.accum->offset();
1483
1484 for (mc_iterate m(i_min, i_max, j_min, j_max, thread); !m.done(); m++) {
1485
1486 unsigned int i = (unsigned int) m.get_i();
1487 unsigned int j = (unsigned int) m.get_j();
1488
1489 if (ref_excluded(i, j, offset, c.ax_parameters, ax_count))
1490 continue;
1491
1492 /*
1493 * Transform
1494 */
1495
1496 struct point q;
1497
1498 q = (c.input_scale < 1.0 && interpolant == NULL)
1499 ? t.scaled_inverse_transform(
1500 point(i + offset[0], j + offset[1]))
1501 : t.unscaled_inverse_transform(
1502 point(i + offset[0], j + offset[1]));
1503
1504 /*
1505 * Check that the transformed coordinates are within
1506 * the boundaries of array c.input, that they are not
1507 * subject to exclusion, and that the weight value in
1508 * the accumulated array is nonzero.
1509 */
1510
1511 if (input_excluded(q[0], q[1], c.ax_parameters, ax_count))
1512 continue;
1513
1514 if (q[0] >= 0
1515 && q[0] <= c.input->height() - 1
1516 && q[1] >= 0
1517 && q[1] <= c.input->width() - 1
1518 && c.certainty->get_pixel(i, j).minabs_norm() != 0) {
1519 pixel a = c.accum->get_pixel(i, j);
1520 pixel b;
1521
1522 b = c.input->get_bl(q);
1523
1524 pixel weight = ((c.aweight && pass_number)
1525 ? c.aweight->get_pixel(i, j)
1526 : pixel(1, 1, 1))
1527 * (pass_number
1528 ? ppow(c.certainty->get_pixel(i, j)
1529 * c.input_certainty->get_bl(q, 1), 0.5)
1530 : pixel(1, 1, 1));
1531
1532 asums[thread] += a * weight;
1533 bsums[thread] += b * weight;
1534 }
1535 }
1536 }
1537
1538 void finish_subdomains(unsigned int N) {
1539 for (unsigned int n = 0; n < N; n++) {
1540 *asum += asums[n];
1541 *bsum += bsums[n];
1542 }
1543 delete asums;
1544 delete bsums;
1545 }
1546 public:
1547 exposure_ratio_iterate(pixel_accum *_asum,
1548 pixel_accum *_bsum,
1549 struct scale_cluster _c,
1550 transformation _t,
1551 int _ax_count,
1552 int _pass_number) : decompose_domain(0, _c.accum->height(),
1553 0, _c.accum->width()){
1554
1555 asum = _asum;
1556 bsum = _bsum;
1557 c = _c;
1558 t = _t;
1559 ax_count = _ax_count;
1560 pass_number = _pass_number;
1561 }
1562 };
1563
1564 static void set_exposure_ratio(unsigned int m, struct scale_cluster c,
1565 transformation t, int ax_count, int pass_number) {
1566
1567 if (_exp_register == 2) {
1568 /*
1569 * Use metadata only.
1570 */
1571 ale_real gain_multiplier = image_rw::exp(m).get_gain_multiplier();
1572 pixel multiplier = pixel(gain_multiplier, gain_multiplier, gain_multiplier);
1573
1574 image_rw::exp(m).set_multiplier(multiplier);
1575 ui::get()->exp_multiplier(multiplier[0],
1576 multiplier[1],
1577 multiplier[2]);
1578
1579 return;
1580 }
1581
1582 pixel_accum asum(0, 0, 0), bsum(0, 0, 0);
1583
1584 exposure_ratio_iterate eri(&asum, &bsum, c, t, ax_count, pass_number);
1585 eri.run();
1586
1587 // std::cerr << (asum / bsum) << " ";
1588
1589 pixel_accum new_multiplier;
1590
1591 new_multiplier = asum / bsum * image_rw::exp(m).get_multiplier();
1592
1593 if (finite(new_multiplier[0])
1594 && finite(new_multiplier[1])
1595 && finite(new_multiplier[2])) {
1596 image_rw::exp(m).set_multiplier(new_multiplier);
1597 ui::get()->exp_multiplier(new_multiplier[0],
1598 new_multiplier[1],
1599 new_multiplier[2]);
1600 }
1601 }
1602
1603 /*
1604 * Copy all ax parameters.
1605 */
1606 static exclusion *copy_ax_parameters(int local_ax_count, exclusion *source) {
1607
1608 exclusion *dest = (exclusion *) malloc(local_ax_count * sizeof(exclusion));
1609
1610 assert (dest);
1611
1612 if (!dest)
1613 ui::get()->memory_error("exclusion regions");
1614
1615 for (int idx = 0; idx < local_ax_count; idx++)
1616 dest[idx] = source[idx];
1617
1618 return dest;
1619 }
1620
1621 /*
1622 * Copy ax parameters according to frame.
1623 */
1624 static exclusion *filter_ax_parameters(int frame, int *local_ax_count) {
1625
1626 exclusion *dest = (exclusion *) malloc(ax_count * sizeof(exclusion));
1627
1628 assert (dest);
1629
1630 if (!dest)
1631 ui::get()->memory_error("exclusion regions");
1632
1633 *local_ax_count = 0;
1634
1635 for (int idx = 0; idx < ax_count; idx++) {
1636 if (ax_parameters[idx].x[4] > frame
1637 || ax_parameters[idx].x[5] < frame)
1638 continue;
1639
1640 dest[*local_ax_count] = ax_parameters[idx];
1641
1642 (*local_ax_count)++;
1643 }
1644
1645 return dest;
1646 }
1647
1648 static void scale_ax_parameters(int local_ax_count, exclusion *ax_parameters,
1649 ale_pos ref_scale, ale_pos input_scale) {
1650 for (int i = 0; i < local_ax_count; i++) {
1651 ale_pos scale = (ax_parameters[i].type == exclusion::RENDER)
1652 ? ref_scale
1653 : input_scale;
1654
1655 for (int n = 0; n < 6; n++) {
1656 ax_parameters[i].x[n] = ax_parameters[i].x[n] * scale;
1657 }
1658 }
1659 }
1660
1661 /*
1662 * Prepare the next level of detail for ordinary images.
1663 */
1664 static const image *prepare_lod(const image *current) {
1665 if (current == NULL)
1666 return NULL;
1667
1668 return current->scale_by_half("prepare_lod");
1669 }
1670
1671 /*
1672 * Prepare the next level of detail for definition maps.
1673 */
1674 static const image *prepare_lod_def(const image *current) {
1675 if (current == NULL)
1676 return NULL;
1677
1678 return current->defined_scale_by_half("prepare_lod_def");
1679 }
1680
1681 /*
1682 * Initialize scale cluster data structures.
1683 */
1684
1685 static void init_nl_cluster(struct scale_cluster *sc) {
1686 }
1687
1688 static struct scale_cluster *init_clusters(int frame, ale_real scale_factor,
1689 const image *input_frame, unsigned int steps,
1690 int *local_ax_count) {
1691
1692 /*
1693 * Allocate memory for the array.
1694 */
1695
1696 struct scale_cluster *scale_clusters =
1697 (struct scale_cluster *) malloc(steps * sizeof(struct scale_cluster));
1698
1699 assert (scale_clusters);
1700
1701 if (!scale_clusters)
1702 ui::get()->memory_error("alignment");
1703
1704 /*
1705 * Prepare images and exclusion regions for the highest level
1706 * of detail.
1707 */
1708
1709 scale_clusters[0].accum = reference_image;
1710
1711 ui::get()->constructing_lod_clusters(0.0);
1712 scale_clusters[0].input_scale = scale_factor;
1713 if (scale_factor < 1.0 && interpolant == NULL)
1714 scale_clusters[0].input = input_frame->scale(scale_factor, "alignment");
1715 else
1716 scale_clusters[0].input = input_frame;
1717
1718 scale_clusters[0].certainty = reference_defined->clone("certainty");
1719 scale_clusters[0].aweight = alignment_weights;
1720 scale_clusters[0].ax_parameters = filter_ax_parameters(frame, local_ax_count);
1721
1722 /*
1723 * Allocate and determine input frame certainty.
1724 */
1725
1726 if (scale_clusters[0].input->get_bayer() != IMAGE_BAYER_NONE) {
1727 scale_clusters[0].input_certainty = new image_bayer_ale_real(
1728 scale_clusters[0].input->height(),
1729 scale_clusters[0].input->width(),
1730 scale_clusters[0].input->depth(),
1731 scale_clusters[0].input->get_bayer());
1732 } else {
1733 scale_clusters[0].input_certainty = scale_clusters[0].input->clone("certainty");
1734 }
1735
1736 for (unsigned int i = 0; i < scale_clusters[0].input_certainty->height(); i++)
1737 for (unsigned int j = 0; j < scale_clusters[0].input_certainty->width(); j++)
1738 for (unsigned int k = 0; k < 3; k++)
1739 if (scale_clusters[0].input->get_channels(i, j) & (1 << k))
1740 ((image *) scale_clusters[0].input_certainty)->chan(i, j, k) =
1741 scale_clusters[0].input->
1742 exp().confidence(scale_clusters[0].input->get_pixel(i, j))[k];
1743
1744 scale_ax_parameters(*local_ax_count, scale_clusters[0].ax_parameters, scale_factor,
1745 (scale_factor < 1.0 && interpolant == NULL) ? scale_factor : 1);
1746
1747 init_nl_cluster(&(scale_clusters[0]));
1748
1749 /*
1750 * Prepare reduced-detail images and exclusion
1751 * regions.
1752 */
1753
1754 for (unsigned int step = 1; step < steps; step++) {
1755 ui::get()->constructing_lod_clusters(step);
1756 scale_clusters[step].accum = prepare_lod(scale_clusters[step - 1].accum);
1757 scale_clusters[step].certainty = prepare_lod_def(scale_clusters[step - 1].certainty);
1758 scale_clusters[step].aweight = prepare_lod_def(scale_clusters[step - 1].aweight);
1759 scale_clusters[step].ax_parameters
1760 = copy_ax_parameters(*local_ax_count, scale_clusters[step - 1].ax_parameters);
1761
1762 double sf = scale_clusters[step - 1].input_scale / 2;
1763 scale_clusters[step].input_scale = sf;
1764
1765 if (sf >= 1.0 || interpolant != NULL) {
1766 scale_clusters[step].input = scale_clusters[step - 1].input;
1767 scale_clusters[step].input_certainty = scale_clusters[step - 1].input_certainty;
1768 scale_ax_parameters(*local_ax_count, scale_clusters[step].ax_parameters, 0.5, 1);
1769 } else if (sf > 0.5) {
1770 scale_clusters[step].input = scale_clusters[step - 1].input->scale(sf, "alignment");
1771 scale_clusters[step].input_certainty = scale_clusters[step - 1].input->scale(sf, "alignment", 1);
1772 scale_ax_parameters(*local_ax_count, scale_clusters[step].ax_parameters, 0.5, sf);
1773 } else {
1774 scale_clusters[step].input = scale_clusters[step - 1].input->scale(0.5, "alignment");
1775 scale_clusters[step].input_certainty = scale_clusters[step - 1].input_certainty->scale(0.5, "alignment", 1);
1776 scale_ax_parameters(*local_ax_count, scale_clusters[step].ax_parameters, 0.5, 0.5);
1777 }
1778
1779 init_nl_cluster(&(scale_clusters[step]));
1780 }
1781
1782 return scale_clusters;
1783 }
1784
1785 /*
1786 * Destroy the first element in the scale cluster data structure.
1787 */
1788 static void final_clusters(struct scale_cluster *scale_clusters, ale_real scale_factor,
1789 unsigned int steps) {
1790
1791 if (scale_clusters[0].input_scale < 1.0)
1792 delete scale_clusters[0].input;
1793
1794 free((void *)scale_clusters[0].ax_parameters);
1795
1796 delete scale_clusters[0].certainty;
1797
1798 for (unsigned int step = 1; step < steps; step++) {
1799 delete scale_clusters[step].accum;
1800 delete scale_clusters[step].certainty;
1801 delete scale_clusters[step].aweight;
1802
1803 if (scale_clusters[step].input_scale < 1.0)
1804 delete scale_clusters[step].input;
1805
1806 free((void *)scale_clusters[step].ax_parameters);
1807 }
1808
1809 free(scale_clusters);
1810 }
1811
1812 /*
1813 * Calculate the centroid of a control point for the set of frames
1814 * having index lower than m. Divide by any scaling of the output.
1815 */
1816 static point unscaled_centroid(unsigned int m, unsigned int p) {
1817 assert(_keep);
1818
1819 point point_sum(0, 0);
1820 ale_accum divisor = 0;
1821
1822 for(unsigned int j = 0; j < m; j++) {
1823 point pp = cp_array[p][j];
1824
1825 if (pp.defined()) {
1826 point_sum += kept_t[j].transform_unscaled(pp)
1827 / kept_t[j].scale();
1828 divisor += 1;
1829 }
1830 }
1831
1832 if (divisor == 0)
1833 return point::undefined();
1834
1835 return point_sum / divisor;
1836 }
1837
1838 /*
1839 * Calculate centroid of this frame, and of all previous frames,
1840 * from points common to both sets.
1841 */
1842 static void centroids(unsigned int m, point *current, point *previous) {
1843 /*
1844 * Calculate the translation
1845 */
1846 point other_centroid(0, 0);
1847 point this_centroid(0, 0);
1848 ale_pos divisor = 0;
1849
1850 for (unsigned int i = 0; i < cp_count; i++) {
1851 point other_c = unscaled_centroid(m, i);
1852 point this_c = cp_array[i][m];
1853
1854 if (!other_c.defined() || !this_c.defined())
1855 continue;
1856
1857 other_centroid += other_c;
1858 this_centroid += this_c;
1859 divisor += 1;
1860
1861 }
1862
1863 if (divisor == 0) {
1864 *current = point::undefined();
1865 *previous = point::undefined();
1866 return;
1867 }
1868
1869 *current = this_centroid / divisor;
1870 *previous = other_centroid / divisor;
1871 }
1872
1873 /*
1874 * Calculate the RMS error of control points for frame m, with
1875 * transformation t, against control points for earlier frames.
1876 */
1877 static ale_accum cp_rms_error(unsigned int m, transformation t) {
1878 assert (_keep);
1879
1880 ale_accum err = 0;
1881 ale_accum divisor = 0;
1882
1883 for (unsigned int i = 0; i < cp_count; i++)
1884 for (unsigned int j = 0; j < m; j++) {
1885 const point *p = cp_array[i];
1886 point p_ref = kept_t[j].transform_unscaled(p[j]);
1887 point p_cur = t.transform_unscaled(p[m]);
1888
1889 if (!p_ref.defined() || !p_cur.defined())
1890 continue;
1891
1892 err += p_ref.lengthtosq(p_cur);
1893 divisor += 1;
1894 }
1895
1896 return sqrt(err / divisor);
1897 }
1898
1899 /*
1900 * Implement new delta --follow semantics.
1901 *
1902 * If we have a transformation T such that
1903 *
1904 * prev_final == T(prev_init)
1905 *
1906 * Then we also have
1907 *
1908 * current_init_follow == T(current_init)
1909 *
1910 * We can calculate T as follows:
1911 *
1912 * T == prev_final(prev_init^-1)
1913 *
1914 * Where ^-1 is the inverse operator.
1915 */
1916 static transformation follow(element_t *element, transformation offset, int lod) {
1917
1918 transformation new_offset = offset;
1919
1920 /*
1921 * Criteria for using following.
1922 */
1923
1924 if (!element->old_is_default && !element->is_default &&
1925 default_initial_alignment_type == 1) {
1926 /*
1927 * Ensure that the lod for the old initial and final
1928 * alignments are equal to the lod for the new initial
1929 * alignment.
1930 */
1931
1932 ui::get()->following();
1933
1934 element->old_final_alignment.rescale (1 / pow(2, lod));
1935 element->old_initial_alignment.rescale(1 / pow(2, lod - element->old_lod));
1936
1937 for (offset.set_current_index(0),
1938 element->old_initial_alignment.set_current_index(0),
1939 element->old_final_alignment.set_current_index(0),
1940 new_offset.set_current_index(0);
1941
1942 offset.get_current_index() < _ma_card;
1943
1944 offset.push_element(),
1945 new_offset.push_element()) {
1946
1947 if (alignment_class == 0) {
1948 /*
1949 * Translational transformations
1950 */
1951
1952 ale_pos t0 = -element->old_initial_alignment.eu_get(0)
1953 + element->old_final_alignment.eu_get(0);
1954 ale_pos t1 = -element->old_initial_alignment.eu_get(1)
1955 + element->old_final_alignment.eu_get(1);
1956
1957 new_offset.eu_modify(0, t0);
1958 new_offset.eu_modify(1, t1);
1959
1960 } else if (alignment_class == 1) {
1961 /*
1962 * Euclidean transformations
1963 */
1964
1965 ale_pos t2 = -element->old_initial_alignment.eu_get(2)
1966 + element->old_final_alignment.eu_get(2);
1967
1968 new_offset.eu_modify(2, t2);
1969
1970 point p( offset.scaled_height()/2 + offset.eu_get(0) - element->old_initial_alignment.eu_get(0),
1971 offset.scaled_width()/2 + offset.eu_get(1) - element->old_initial_alignment.eu_get(1) );
1972
1973 p = element->old_final_alignment.transform_scaled(p);
1974
1975 new_offset.eu_modify(0, p[0] - offset.scaled_height()/2 - offset.eu_get(0));
1976 new_offset.eu_modify(1, p[1] - offset.scaled_width()/2 - offset.eu_get(1));
1977
1978 } else if (alignment_class == 2) {
1979 /*
1980 * Projective transformations
1981 */
1982
1983 point p[4];
1984
1985 p[0] = element->old_final_alignment.transform_scaled(element->old_initial_alignment
1986 . scaled_inverse_transform(offset.get_current_element().transform_scaled(point( 0 , 0 ))));
1987 p[1] = element->old_final_alignment.transform_scaled(element->old_initial_alignment
1988 . scaled_inverse_transform(offset.get_current_element().transform_scaled(point(offset.scaled_height(), 0 ))));
1989 p[2] = element->old_final_alignment.transform_scaled(element->old_initial_alignment
1990 . scaled_inverse_transform(offset.get_current_element().transform_scaled(point(offset.scaled_height(), offset.scaled_width()))));
1991 p[3] = element->old_final_alignment.transform_scaled(element->old_initial_alignment
1992 . scaled_inverse_transform(offset.get_current_element().transform_scaled(point( 0 , offset.scaled_width()))));
1993
1994 new_offset.gpt_set(p);
1995 }
1996 }
1997
1998 ui::get()->set_offset(offset);
1999 }
2000
2001 return new_offset;
2002 }
2003
2004 static void test_global(diff_stat_t *here, scale_cluster si, transformation t,
2005 int local_ax_count, int m, ale_pos local_gs_mo, ale_pos perturb) {
2006
2007 diff_stat_t test(*here);
2008
2009 test.diff(si, perturb, t, local_ax_count, m);
2010
2011 unsigned int ovl = overlap(si, t, local_ax_count);
2012
2013 if (ovl >= local_gs_mo && test.better()) {
2014 test.confirm();
2015 *here = test;
2016 ui::get()->set_match(here->get_error());
2017 ui::get()->set_offset(here->get_offset());
2018 } else {
2019 test.discard();
2020 }
2021
2022 }
2023
2024 /*
2025 * Get the set of global transformations for a given density
2026 */
2027 static void test_globals(diff_stat_t *here,
2028 scale_cluster si, transformation t, int local_gs, ale_pos adj_p,
2029 int local_ax_count, int m, ale_pos local_gs_mo, ale_pos perturb) {
2030
2031 transformation offset = t;
2032
2033 point min, max;
2034
2035 transformation offset_p = offset;
2036
2037 if (!offset_p.is_projective())
2038 offset_p.eu_to_gpt();
2039
2040 min = max = offset_p.gpt_get(0);
2041 for (int p_index = 1; p_index < 4; p_index++) {
2042 point p = offset_p.gpt_get(p_index);
2043 if (p[0] < min[0])
2044 min[0] = p[0];
2045 if (p[1] < min[1])
2046 min[1] = p[1];
2047 if (p[0] > max[0])
2048 max[0] = p[0];
2049 if (p[1] > max[1])
2050 max[1] = p[1];
2051 }
2052
2053 point inner_min_t = -min;
2054 point inner_max_t = -max + point(si.accum->height(), si.accum->width());
2055 point outer_min_t = -max + point(adj_p - 1, adj_p - 1);
2056 point outer_max_t = point(si.accum->height(), si.accum->width()) - point(adj_p, adj_p);
2057
2058 if (local_gs == 1 || local_gs == 3 || local_gs == 4 || local_gs == 6) {
2059
2060 /*
2061 * Inner
2062 */
2063
2064 for (ale_pos i = inner_min_t[0]; i <= inner_max_t[0]; i += adj_p)
2065 for (ale_pos j = inner_min_t[1]; j <= inner_max_t[1]; j += adj_p) {
2066 transformation test_t = offset;
2067 test_t.translate(point(i, j));
2068 test_global(here, si, test_t, local_ax_count, m, local_gs_mo, perturb);
2069 }
2070 }
2071
2072 if (local_gs == 2 || local_gs == 3 || local_gs == -1 || local_gs == 6) {
2073
2074 /*
2075 * Outer
2076 */
2077
2078 for (ale_pos i = outer_min_t[0]; i <= outer_max_t[0]; i += adj_p)
2079 for (ale_pos j = outer_min_t[1]; j < inner_min_t[1]; j += adj_p) {
2080 transformation test_t = offset;
2081 test_t.translate(point(i, j));
2082 test_global(here, si, test_t, local_ax_count, m, local_gs_mo, perturb);
2083 }
2084 for (ale_pos i = outer_min_t[0]; i <= outer_max_t[0]; i += adj_p)
2085 for (ale_pos j = outer_max_t[1]; j > inner_max_t[1]; j -= adj_p) {
2086 transformation test_t = offset;
2087 test_t.translate(point(i, j));
2088 test_global(here, si, test_t, local_ax_count, m, local_gs_mo, perturb);
2089 }
2090 for (ale_pos i = outer_min_t[0]; i < inner_min_t[0]; i += adj_p)
2091 for (ale_pos j = outer_min_t[1]; j <= outer_max_t[1]; j += adj_p) {
2092 transformation test_t = offset;
2093 test_t.translate(point(i, j));
2094 test_global(here, si, test_t, local_ax_count, m, local_gs_mo, perturb);
2095 }
2096 for (ale_pos i = outer_max_t[0]; i > inner_max_t[0]; i -= adj_p)
2097 for (ale_pos j = outer_min_t[1]; j <= outer_max_t[1]; j += adj_p) {
2098 transformation test_t = offset;
2099 test_t.translate(point(i, j));
2100 test_global(here, si, test_t, local_ax_count, m, local_gs_mo, perturb);
2101 }
2102 }
2103 }
2104
2105 static void get_translational_set(std::vector<transformation> *set,
2106 transformation t, ale_pos adj_p) {
2107
2108 ale_pos adj_s;
2109
2110 transformation offset = t;
2111 transformation test_t;
2112
2113 for (int i = 0; i < 2; i++)
2114 for (adj_s = -adj_p; adj_s <= adj_p; adj_s += 2 * adj_p) {
2115
2116 test_t = offset;
2117
2118 test_t.translate(i ? point(adj_s, 0) : point(0, adj_s));
2119
2120 set->push_back(test_t);
2121 }
2122 }
2123
2124 static int threshold_ok(ale_accum error) {
2125 if ((1 - error) * 100 >= match_threshold)
2126 return 1;
2127
2128 if (!(match_threshold >= 0))
2129 return 1;
2130
2131 return 0;
2132 }
2133
2134 /*
2135 * Align frame m against the reference.
2136 *
2137 * XXX: the transformation class currently combines ordinary
2138 * transformations with scaling. This is somewhat convenient for
2139 * some things, but can also be confusing. This method, _align(), is
2140 * one case where special care must be taken to ensure that the scale
2141 * is always set correctly (by using the 'rescale' method).
2142 */
2143 static diff_stat_t _align(int m, int local_gs, element_t *element) {
2144
2145 const image *input_frame = element->input_frame;
2146
2147 /*
2148 * Local upper/lower data, possibly dependent on image
2149 * dimensions.
2150 */
2151
2152 ale_pos local_lower, local_upper, local_gs_mo;
2153
2154 /*
2155 * Select the minimum dimension as the reference.
2156 */
2157
2158 ale_pos reference_size = input_frame->height();
2159 if (input_frame->width() < reference_size)
2160 reference_size = input_frame->width();
2161 ale_pos reference_area = input_frame->height()
2162 * input_frame->width();
2163
2164 if (perturb_lower_percent)
2165 local_lower = perturb_lower
2166 * reference_size
2167 * 0.01
2168 * scale_factor;
2169 else
2170 local_lower = perturb_lower;
2171
2172 if (perturb_upper_percent)
2173 local_upper = perturb_upper
2174 * reference_size
2175 * 0.01
2176 * scale_factor;
2177 else
2178 local_upper = perturb_upper;
2179
2180 local_upper = pow(2, floor(log(local_upper) / log(2)));
2181
2182 if (gs_mo_percent)
2183 local_gs_mo = _gs_mo
2184 * reference_area
2185 * 0.01
2186 * scale_factor;
2187 else
2188 local_gs_mo = _gs_mo;
2189
2190 /*
2191 * Logarithms aren't exact, so we divide repeatedly to discover
2192 * how many steps will occur, and pass this information to the
2193 * user interface.
2194 */
2195
2196 int step_count = 0;
2197 double step_variable = local_upper;
2198 while (step_variable >= local_lower) {
2199 step_variable /= 2;
2200 step_count++;
2201 }
2202
2203 ui::get()->set_steps(step_count);
2204
2205 ale_pos perturb = local_upper;
2206 int lod;
2207
2208 if (_keep) {
2209 kept_t[latest] = latest_t;
2210 kept_ok[latest] = latest_ok;
2211 }
2212
2213 /*
2214 * Maximum level-of-detail. Use a level of detail at most
2215 * 2^lod_diff finer than the adjustment resolution. lod_diff
2216 * is a synonym for lod_max.
2217 */
2218
2219 const int lod_diff = lod_max;
2220
2221 /*
2222 * Determine how many levels of detail should be prepared.
2223 */
2224
2225 /*
2226 * Plain (unsigned int) casting seems to be broken in some cases.
2227 */
2228
2229 unsigned int steps = (perturb > pow(2, lod_max))
2230 ? (unsigned int) lrint(log(perturb) / log(2)) - lod_max + 1 : 1;
2231
2232 /*
2233 * Prepare multiple levels of detail.
2234 */
2235
2236 int local_ax_count;
2237 struct scale_cluster *scale_clusters = init_clusters(m,
2238 scale_factor, input_frame, steps,
2239 &local_ax_count);
2240
2241 /*
2242 * Initialize variables used in the main loop.
2243 */
2244
2245 lod = (steps - 1);
2246
2247 /*
2248 * Initialize the default initial transform
2249 */
2250
2251 if (default_initial_alignment_type == 0) {
2252
2253 /*
2254 * Follow the transformation of the original frame,
2255 * setting new image dimensions.
2256 */
2257
2258 // element->default_initial_alignment = orig_t;
2259 element->default_initial_alignment.set_current_element(orig_t.get_element(0));
2260 element->default_initial_alignment.set_dimensions(input_frame);
2261
2262 } else if (default_initial_alignment_type == 1)
2263
2264 /*
2265 * Follow previous transformation, setting new image
2266 * dimensions.
2267 */
2268
2269 element->default_initial_alignment.set_dimensions(input_frame);
2270
2271 else
2272 assert(0);
2273
2274 element->old_is_default = element->is_default;
2275
2276 /*
2277 * Scale default initial transform for lod
2278 */
2279
2280 element->default_initial_alignment.rescale(1 / pow(2, lod));
2281
2282 /*
2283 * Set the default transformation.
2284 */
2285
2286 transformation offset = element->default_initial_alignment;
2287
2288 /*
2289 * Load any file-specified transformations
2290 */
2291
2292 for (offset.set_current_index(0);
2293 offset.get_current_index() < _ma_card;
2294 offset.push_element()) {
2295
2296 offset = tload_next(tload, alignment_class == 2,
2297 offset,
2298 &element->is_default,
2299 offset.get_current_index() == 0);
2300
2301 }
2302
2303 offset.set_current_index(0);
2304
2305 ui::get()->set_offset(offset);
2306
2307 if (perturb > 0) {
2308
2309 /*
2310 * Apply following logic
2311 */
2312
2313 transformation new_offset = follow(element, offset, lod);
2314
2315 new_offset.set_current_index(0);
2316
2317 element->old_initial_alignment = offset;
2318 element->old_lod = lod;
2319 offset = new_offset;
2320
2321 } else {
2322 element->old_initial_alignment = offset;
2323 element->old_lod = lod;
2324 }
2325
2326 struct scale_cluster si = scale_clusters[lod];
2327
2328 /*
2329 * Projective adjustment value
2330 */
2331
2332 ale_pos adj_p = (perturb >= pow(2, lod_diff))
2333 ? pow(2, lod_diff) : (double) perturb;
2334
2335 /*
2336 * Orientational adjustment value in degrees.
2337 *
2338 * Since rotational perturbation is now specified as an
2339 * arclength, we have to convert.
2340 */
2341
2342 ale_pos adj_o = 2 * perturb
2343 / sqrt(pow(scale_clusters[0].input->height(), 2)
2344 + pow(scale_clusters[0].input->width(), 2))
2345 * 180
2346 / M_PI;
2347
2348 /*
2349 * Barrel distortion adjustment value
2350 */
2351
2352 ale_pos adj_b = perturb * bda_mult;
2353
2354 /*
2355 * Global search overlap requirements.
2356 */
2357
2358 local_gs_mo /= pow(pow(2, lod), 2);
2359
2360 /*
2361 * Pre-alignment exposure adjustment
2362 */
2363
2364 if (_exp_register) {
2365 ui::get()->exposure_1();
2366 transformation o = offset;
2367 for (int k = lod; k > 0; k--)
2368 o.rescale(2);
2369 set_exposure_ratio(m, scale_clusters[0], o, local_ax_count, 0);
2370 }
2371
2372 /*
2373 * Alignment statistics.
2374 */
2375
2376 diff_stat_t here;
2377
2378 /*
2379 * Current difference (error) value
2380 */
2381
2382 ui::get()->prematching();
2383 here.diff(si, perturb, offset, local_ax_count, m);
2384 ui::get()->set_match(here.get_error());
2385
2386 /*
2387 * Current and modified barrel distortion parameters
2388 */
2389
2390 ale_pos current_bd[BARREL_DEGREE];
2391 ale_pos modified_bd[BARREL_DEGREE];
2392 offset.bd_get(current_bd);
2393 offset.bd_get(modified_bd);
2394
2395 /*
2396 * Translational global search step
2397 */
2398
2399 if (perturb >= local_lower && local_gs != 0 && local_gs != 5
2400 && (local_gs != 6 || element->is_default)) {
2401
2402 ui::get()->global_alignment(perturb, lod);
2403 ui::get()->gs_mo(local_gs_mo);
2404
2405 test_globals(&here, si, here.get_offset(), local_gs, adj_p,
2406 local_ax_count, m, local_gs_mo, perturb);
2407
2408 ui::get()->set_match(here.get_error());
2409 ui::get()->set_offset(here.get_offset());
2410 }
2411
2412 /*
2413 * Control point alignment
2414 */
2415
2416 if (local_gs == 5) {
2417
2418 transformation o = here.get_offset();
2419
2420 for (int k = lod; k > 0; k--)
2421 o.rescale(2);
2422
2423 /*
2424 * Determine centroid data
2425 */
2426
2427 point current, previous;
2428 centroids(m, &current, &previous);
2429
2430 if (current.defined() && previous.defined()) {
2431 o = orig_t;
2432 o.set_dimensions(input_frame);
2433 o.translate((previous - current) * o.scale());
2434 current = previous;
2435 }
2436
2437 /*
2438 * Determine rotation for alignment classes other than translation.
2439 */
2440
2441 ale_accum lowest_error = cp_rms_error(m, o);
2442
2443 ale_pos rot_lower = 2 * local_lower
2444 / sqrt(pow(scale_clusters[0].input->height(), 2)
2445 + pow(scale_clusters[0].input->width(), 2))
2446 * 180
2447 / M_PI;
2448
2449 if (alignment_class > 0)
2450 for (ale_pos rot = 30; rot > rot_lower; rot /= 2)
2451 for (ale_pos srot = -rot; srot < rot * 1.5; srot += rot * 2) {
2452 int is_improved = 1;
2453 while (is_improved) {
2454 is_improved = 0;
2455 transformation test_t = o;
2456 /*
2457 * XXX: is this right?
2458 */
2459 test_t.rotate(current * o.scale(), srot);
2460 ale_pos test_v = cp_rms_error(m, test_t);
2461
2462 if (test_v < lowest_error) {
2463 lowest_error = test_v;
2464 o = test_t;
2465 srot += 3 * rot;
2466 is_improved = 1;
2467 }
2468 }
2469 }
2470
2471 /*
2472 * Determine projective parameters through a local
2473 * minimum search.
2474 */
2475
2476 if (alignment_class == 2) {
2477 ale_accum adj_p = lowest_error;
2478
2479 if (adj_p < local_lower)
2480 adj_p = local_lower;
2481
2482 while (adj_p >= local_lower) {
2483 transformation test_t = o;
2484 int is_improved = 1;
2485 ale_accum test_v;
2486 ale_accum adj_s;
2487
2488 while (is_improved) {
2489 is_improved = 0;
2490
2491 for (int i = 0; i < 4; i++)
2492 for (int j = 0; j < 2; j++)
2493 for (adj_s = -adj_p; adj_s <= adj_p; adj_s += 2 * adj_p) {
2494
2495 test_t = o;
2496
2497 if (perturb_type == 0)
2498 test_t.gpt_modify(j, i, adj_s);
2499 else if (perturb_type == 1)
2500 test_t.gr_modify(j, i, adj_s);
2501 else
2502 assert(0);
2503
2504 test_v = cp_rms_error(m, test_t);
2505
2506 if (test_v < lowest_error) {
2507 lowest_error = test_v;
2508 o = test_t;
2509 adj_s += 3 * adj_p;
2510 is_improved = 1;
2511 }
2512 }
2513 }
2514 adj_p /= 2;
2515 }
2516 }
2517
2518 if (_exp_register)
2519 set_exposure_ratio(m, scale_clusters[0], o, local_ax_count, 0);
2520
2521 for (int k = lod; k > 0; k--)
2522 o.rescale(0.5);
2523
2524 here.diff(si, perturb, o, local_ax_count, m);
2525 here.confirm();
2526 ui::get()->set_match(here.get_error());
2527 ui::get()->set_offset(here.get_offset());
2528 }
2529
2530 /*
2531 * Announce perturbation size
2532 */
2533
2534 ui::get()->aligning(perturb, lod);
2535
2536 /*
2537 * Run initial tests to get perturbation multipliers and error
2538 * centroids.
2539 */
2540
2541 std::vector<transformation> t_set;
2542
2543 here.get_perturb_set(&t_set, adj_p, adj_o, adj_b, current_bd, modified_bd);
2544
2545 /*
2546 * Perturbation adjustment loop.
2547 */
2548
2549 int stable_count = 0;
2550
2551 while (perturb >= local_lower) {
2552
2553 /*
2554 * Orientational adjustment value in degrees.
2555 *
2556 * Since rotational perturbation is now specified as an
2557 * arclength, we have to convert.
2558 */
2559
2560 ale_pos adj_o = 2 * perturb
2561 / sqrt(pow(scale_clusters[0].input->height(), 2)
2562 + pow(scale_clusters[0].input->width(), 2))
2563 * 180
2564 / M_PI;
2565
2566 /*
2567 * Barrel distortion adjustment value
2568 */
2569
2570 ale_pos adj_b = perturb * bda_mult;
2571
2572 diff_stat_t old_here = here;
2573
2574 here.perturb_test(perturb, adj_p, adj_o, adj_b, current_bd, modified_bd,
2575 stable_count);
2576
2577 if (here.get_offset() == old_here.get_offset())
2578 stable_count++;
2579 else
2580 stable_count = 0;
2581
2582 if (stable_count == 3) {
2583
2584 stable_count = 0;
2585
2586 here.calculate_element_region();
2587
2588 if (here.get_current_index() + 1 < _ma_card) {
2589 here.push_element();
2590 here.make_element_nontrivial(adj_p, adj_o);
2591 element->is_primary = 0;
2592 } else {
2593
2594 here.set_current_index(0);
2595
2596 element->is_primary = 1;
2597
2598 perturb *= 0.5;
2599
2600 if (lod > 0) {
2601
2602 /*
2603 * Work with images twice as large
2604 */
2605
2606 lod--;
2607 si = scale_clusters[lod];
2608
2609 /*
2610 * Rescale the transforms.
2611 */
2612
2613 here.rescale(2, si);
2614 element->default_initial_alignment.rescale(2);
2615
2616 } else {
2617 adj_p = perturb;
2618 }
2619
2620 /*
2621 * Announce changes
2622 */
2623
2624 ui::get()->alignment_perturbation_level(perturb, lod);
2625
2626 }
2627 }
2628
2629 ui::get()->set_match(here.get_error());
2630 ui::get()->set_offset(here.get_offset());
2631 }
2632
2633 here.set_current_index(0);
2634
2635 offset = here.get_offset();
2636
2637 if (lod > 0) {
2638 here.rescale(pow(2, lod), scale_clusters[0]);
2639 element->default_initial_alignment.rescale(pow(2, lod));
2640 }
2641
2642 /*
2643 * Post-alignment exposure adjustment
2644 */
2645
2646 if (_exp_register == 1) {
2647 ui::get()->exposure_2();
2648 set_exposure_ratio(m, scale_clusters[0], offset, local_ax_count, 1);
2649 }
2650
2651 /*
2652 * Recalculate error
2653 */
2654
2655 ui::get()->postmatching();
2656 offset.use_full_support();
2657 here.diff(scale_clusters[0], perturb, offset, local_ax_count, m);
2658 here.confirm();
2659 offset.use_restricted_support();
2660 ui::get()->set_match(here.get_error());
2661
2662 /*
2663 * Free the level-of-detail structures
2664 */
2665
2666 final_clusters(scale_clusters, scale_factor, steps);
2667
2668 /*
2669 * Ensure that the match meets the threshold.
2670 */
2671
2672 if (threshold_ok(here.get_error())) {
2673 /*
2674 * Update alignment variables
2675 */
2676 latest_ok = 1;
2677 element->default_initial_alignment = offset;
2678 element->old_final_alignment = offset;
2679 ui::get()->alignment_match_ok();
2680 } else if (local_gs == 4) {
2681
2682 /*
2683 * Align with outer starting points.
2684 */
2685
2686 /*
2687 * XXX: This probably isn't exactly the right thing to do,
2688 * since variables like old_initial_value have been overwritten.
2689 */
2690
2691 diff_stat_t nested_result = _align(m, -1, element);
2692
2693 if (threshold_ok(nested_result.get_error())) {
2694 return nested_result;
2695 } else if (nested_result.get_error() < here.get_error()) {
2696 here = nested_result;
2697 }
2698
2699 if (is_fail_default)
2700 offset = element->default_initial_alignment;
2701
2702 ui::get()->set_match(here.get_error());
2703 ui::get()->alignment_no_match();
2704
2705 } else if (local_gs == -1) {
2706
2707 latest_ok = 0;
2708 latest_t = offset;
2709 return here;
2710
2711 } else {
2712 if (is_fail_default)
2713 offset = element->default_initial_alignment;
2714 latest_ok = 0;
2715 ui::get()->alignment_no_match();
2716 }
2717
2718 /*
2719 * Write the tonal registration multiplier as a comment.
2720 */
2721
2722 pixel trm = image_rw::exp(m).get_multiplier();
2723 tsave_trm(tsave, trm[0], trm[1], trm[2]);
2724
2725 /*
2726 * Save the transformation information
2727 */
2728
2729 for (offset.set_current_index(0);
2730 offset.get_current_index() < _ma_card;
2731 offset.push_element()) {
2732
2733 tsave_next(tsave, offset, alignment_class == 2,
2734 offset.get_current_index() == 0);
2735 }
2736
2737 offset.set_current_index(0);
2738
2739 latest_t = offset;
2740
2741 /*
2742 * Update match statistics.
2743 */
2744
2745 match_sum += (1 - here.get_error()) * 100;
2746 match_count++;
2747 latest = m;
2748
2749 return here;
2750 }
2751
2752 #ifdef USE_FFTW
2753 /*
2754 * High-pass filter for frequency weights
2755 */
2756 static void hipass(int rows, int cols, fftw_complex *inout) {
2757 for (int i = 0; i < rows * vert_freq_cut; i++)
2758 for (int j = 0; j < cols; j++)
2759 for (int k = 0; k < 2; k++)
2760 inout[i * cols + j][k] = 0;
2761 for (int i = 0; i < rows; i++)
2762 for (int j = 0; j < cols * horiz_freq_cut; j++)
2763 for (int k = 0; k < 2; k++)
2764 inout[i * cols + j][k] = 0;
2765 for (int i = 0; i < rows; i++)
2766 for (int j = 0; j < cols; j++)
2767 for (int k = 0; k < 2; k++)
2768 if (i / (double) rows + j / (double) cols < 2 * avg_freq_cut)
2769 inout[i * cols + j][k] = 0;
2770 }
2771 #endif
2772
2773
2774 /*
2775 * Reset alignment weights
2776 */
2777 static void reset_weights() {
2778 if (alignment_weights != NULL)
2779 delete alignment_weights;
2780
2781 alignment_weights = NULL;
2782 }
2783
2784 /*
2785 * Initialize alignment weights
2786 */
2787 static void init_weights() {
2788 if (alignment_weights != NULL)
2789 return;
2790
2791 int rows = reference_image->height();
2792 int cols = reference_image->width();
2793 int colors = reference_image->depth();
2794
2795 alignment_weights = new image_ale_real(rows, cols,
2796 colors, "alignment_weights");
2797
2798 assert (alignment_weights);
2799
2800 for (int i = 0; i < rows; i++)
2801 for (int j = 0; j < cols; j++)
2802 alignment_weights->set_pixel(i, j, pixel(1, 1, 1));
2803 }
2804
2805 /*
2806 * Update alignment weights with weight map
2807 */
2808 static void map_update() {
2809
2810 if (weight_map == NULL)
2811 return;
2812
2813 init_weights();
2814
2815 point map_offset = reference_image->offset() - weight_map->offset();
2816
2817 int rows = reference_image->height();
2818 int cols = reference_image->width();
2819
2820 for (int i = 0; i < rows; i++)
2821 for (int j = 0; j < cols; j++) {
2822 point map_weight_position = map_offset + point(i, j);
2823 if (map_weight_position[0] >= 0
2824 && map_weight_position[1] >= 0
2825 && map_weight_position[0] <= weight_map->height() - 1
2826 && map_weight_position[1] <= weight_map->width() - 1)
2827 alignment_weights->pix(i, j) *= weight_map->get_bl(map_weight_position);
2828 }
2829 }
2830
2831 /*
2832 * Update alignment weights with algorithmic weights
2833 */
2834 static void wmx_update() {
2835 #ifdef USE_UNIX
2836
2837 static exposure *exp_def = new exposure_default();
2838 static exposure *exp_bool = new exposure_boolean();
2839
2840 if (wmx_file == NULL || wmx_exec == NULL || wmx_defs == NULL)
2841 return;
2842
2843 unsigned int rows = reference_image->height();
2844 unsigned int cols = reference_image->width();
2845
2846 image_rw::write_image(wmx_file, reference_image);
2847 image_rw::write_image(wmx_defs, reference_defined, exp_bool);
2848
2849 /* execute ... */
2850 int exit_status = 1;
2851 if (!fork()) {
2852 execlp(wmx_exec, wmx_exec, wmx_file, wmx_defs, NULL);
2853 ui::get()->exec_failure(wmx_exec, wmx_file, wmx_defs);
2854 }
2855
2856 wait(&exit_status);
2857
2858 if (exit_status)
2859 ui::get()->fork_failure("d2::align");
2860
2861 image *wmx_weights = image_rw::read_image(wmx_file, exp_def);
2862
2863 if (wmx_weights->height() != rows || wmx_weights->width() != cols)
2864 ui::get()->error("algorithmic weighting must not change image size");
2865
2866 if (alignment_weights == NULL)
2867 alignment_weights = wmx_weights;
2868 else
2869 for (unsigned int i = 0; i < rows; i++)
2870 for (unsigned int j = 0; j < cols; j++)
2871 alignment_weights->pix(i, j) *= wmx_weights->pix(i, j);
2872 #endif
2873 }
2874
2875 /*
2876 * Update alignment weights with frequency weights
2877 */
2878 static void fw_update() {
2879 #ifdef USE_FFTW
2880 if (horiz_freq_cut == 0
2881 && vert_freq_cut == 0
2882 && avg_freq_cut == 0)
2883 return;
2884
2885 /*
2886 * Required for correct operation of --fwshow
2887 */
2888
2889 assert (alignment_weights == NULL);
2890
2891 int rows = reference_image->height();
2892 int cols = reference_image->width();
2893 int colors = reference_image->depth();
2894
2895 alignment_weights = new image_ale_real(rows, cols,
2896 colors, "alignment_weights");
2897
2898 fftw_complex *inout = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * rows * cols);
2899
2900 assert (inout);
2901
2902 fftw_plan pf = fftw_plan_dft_2d(rows, cols,
2903 inout, inout,
2904 FFTW_FORWARD, FFTW_ESTIMATE);
2905
2906 fftw_plan pb = fftw_plan_dft_2d(rows, cols,
2907 inout, inout,
2908 FFTW_BACKWARD, FFTW_ESTIMATE);
2909
2910 for (int k = 0; k < colors; k++) {
2911 for (int i = 0; i < rows * cols; i++) {
2912 inout[i][0] = reference_image->get_pixel(i / cols, i % cols)[k];
2913 inout[i][1] = 0;
2914 }
2915
2916 fftw_execute(pf);
2917 hipass(rows, cols, inout);
2918 fftw_execute(pb);
2919
2920 for (int i = 0; i < rows * cols; i++) {
2921 #if 0
2922 alignment_weights->pix(i / cols, i % cols)[k] = fabs(inout[i][0] / (rows * cols));
2923 #else
2924 alignment_weights->pix(i / cols, i % cols)[k] =
2925 sqrt(pow(inout[i][0] / (rows * cols), 2)
2926 + pow(inout[i][1] / (rows * cols), 2));
2927 #endif
2928 }
2929 }
2930
2931 fftw_destroy_plan(pf);
2932 fftw_destroy_plan(pb);
2933 fftw_free(inout);
2934
2935 if (fw_output != NULL)
2936 image_rw::write_image(fw_output, alignment_weights);
2937 #endif
2938 }
2939
2940 /*
2941 * Update alignment to frame N.
2942 */
2943 static void update_to(int n) {
2944
2945 assert (n <= latest + 1);
2946 assert (n >= 0);
2947
2948 static std::vector<element_t> elements;
2949
2950 if (latest < 0) {
2951
2952 elements.resize(1);
2953
2954 /*
2955 * Handle the initial frame
2956 */
2957
2958 elements[0].input_frame = image_rw::open(n);
2959
2960 const image *i = elements[0].input_frame;
2961 int is_default;
2962 transformation result = alignment_class == 2
2963 ? transformation::gpt_identity(i, scale_factor)
2964 : transformation::eu_identity(i, scale_factor);
2965 result = tload_first(tload, alignment_class == 2, result, &is_default);
2966 tsave_first(tsave, result, alignment_class == 2);
2967
2968 if (_keep > 0) {
2969 kept_t = new transformation[image_rw::count()];
2970 kept_ok = (int *) malloc(image_rw::count()
2971 * sizeof(int));
2972 assert (kept_t);
2973 assert (kept_ok);
2974
2975 if (!kept_t || !kept_ok)
2976 ui::get()->memory_error("alignment");
2977
2978 kept_ok[0] = 1;
2979 kept_t[0] = result;
2980 }
2981
2982 latest = 0;
2983 latest_ok = 1;
2984 latest_t = result;
2985
2986 elements[0].default_initial_alignment = result;
2987 orig_t = result;
2988
2989 image_rw::close(n);
2990 }
2991
2992 for (int i = latest + 1; i <= n; i++) {
2993 int j = 0;
2994
2995 /*
2996 * Handle supplemental frames.
2997 */
2998
2999 assert (reference != NULL);
3000
3001 ui::get()->set_arender_current();
3002 reference->sync(i - 1);
3003 ui::get()->clear_arender_current();
3004 reference_image = reference->get_image();
3005 reference_defined = reference->get_defined();
3006
3007 reset_weights();
3008 fw_update();
3009 wmx_update();
3010 map_update();
3011
3012 assert (reference_image != NULL);
3013 assert (reference_defined != NULL);
3014
3015 elements[j].input_frame = image_rw::open(i);
3016 elements[j].is_primary = 1;
3017
3018 _align(i, _gs, &elements[j]);
3019
3020 image_rw::close(n);
3021 }
3022
3023 if (elements.size() > _ma_card)
3024 elements.resize(_ma_card);
3025 }
3026
3027 public:
3028
3029 /*
3030 * Set the control point count
3031 */
3032 static void set_cp_count(unsigned int n) {
3033 assert (cp_array == NULL);
3034
3035 cp_count = n;
3036 cp_array = (const point **) malloc(n * sizeof(const point *));
3037 }
3038
3039 /*
3040 * Set control points.
3041 */
3042 static void set_cp(unsigned int i, const point *p) {
3043 cp_array[i] = p;
3044 }
3045
3046 /*
3047 * Register exposure
3048 */
3049 static void exp_register() {
3050 _exp_register = 1;
3051 }
3052
3053 /*
3054 * Register exposure only based on metadata
3055 */
3056 static void exp_meta_only() {
3057 _exp_register = 2;
3058 }
3059
3060 /*
3061 * Don't register exposure
3062 */
3063 static void exp_noregister() {
3064 _exp_register = 0;
3065 }
3066
3067 /*
3068 * Set alignment class to translation only. Only adjust the x and y
3069 * position of images. Do not rotate input images or perform
3070 * projective transformations.
3071 */
3072 static void class_translation() {
3073 alignment_class = 0;
3074 }
3075
3076 /*
3077 * Set alignment class to Euclidean transformations only. Adjust the x
3078 * and y position of images and the orientation of the image about the
3079 * image center.
3080 */
3081 static void class_euclidean() {
3082 alignment_class = 1;
3083 }
3084
3085 /*
3086 * Set alignment class to perform general projective transformations.
3087 * See the file gpt.h for more information about general projective
3088 * transformations.
3089 */
3090 static void class_projective() {
3091 alignment_class = 2;
3092 }
3093
3094 /*
3095 * Set the default initial alignment to the identity transformation.
3096 */
3097 static void initial_default_identity() {
3098 default_initial_alignment_type = 0;
3099 }
3100
3101 /*
3102 * Set the default initial alignment to the most recently matched
3103 * frame's final transformation.
3104 */
3105 static void initial_default_follow() {
3106 default_initial_alignment_type = 1;
3107 }
3108
3109 /*
3110 * Perturb output coordinates.
3111 */
3112 static void perturb_output() {
3113 perturb_type = 0;
3114 }
3115
3116 /*
3117 * Perturb source coordinates.
3118 */
3119 static void perturb_source() {
3120 perturb_type = 1;
3121 }
3122
3123 /*
3124 * Frames under threshold align optimally
3125 */
3126 static void fail_optimal() {
3127 is_fail_default = 0;
3128 }
3129
3130 /*
3131 * Frames under threshold keep their default alignments.
3132 */
3133 static void fail_default() {
3134 is_fail_default = 1;
3135 }
3136
3137 /*
3138 * Align images with an error contribution from each color channel.
3139 */
3140 static void all() {
3141 channel_alignment_type = 0;
3142 }
3143
3144 /*
3145 * Align images with an error contribution only from the green channel.
3146 * Other color channels do not affect alignment.
3147 */
3148 static void green() {
3149 channel_alignment_type = 1;
3150 }
3151
3152 /*
3153 * Align images using a summation of channels. May be useful when
3154 * dealing with images that have high frequency color ripples due to
3155 * color aliasing.
3156 */
3157 static void sum() {
3158 channel_alignment_type = 2;
3159 }
3160
3161 /*
3162 * Error metric exponent
3163 */
3164
3165 static void set_metric_exponent(float me) {
3166 metric_exponent = me;
3167 }
3168
3169 /*
3170 * Match threshold
3171 */
3172
3173 static void set_match_threshold(float mt) {
3174 match_threshold = mt;
3175 }
3176
3177 /*
3178 * Perturbation lower and upper bounds.
3179 */
3180
3181 static void set_perturb_lower(ale_pos pl, int plp) {
3182 perturb_lower = pl;
3183 perturb_lower_percent = plp;
3184 }
3185
3186 static void set_perturb_upper(ale_pos pu, int pup) {
3187 perturb_upper = pu;
3188 perturb_upper_percent = pup;
3189 }
3190
3191 /*
3192 * Maximum rotational perturbation.
3193 */
3194
3195 static void set_rot_max(int rm) {
3196
3197 /*
3198 * Obtain the largest power of two not larger than rm.
3199 */
3200
3201 rot_max = pow(2, floor(log(rm) / log(2)));
3202 }
3203
3204 /*
3205 * Barrel distortion adjustment multiplier
3206 */
3207
3208 static void set_bda_mult(ale_pos m) {
3209 bda_mult = m;
3210 }
3211
3212 /*
3213 * Barrel distortion maximum rate of change
3214 */
3215
3216 static void set_bda_rate(ale_pos m) {
3217 bda_rate = m;
3218 }
3219
3220 /*
3221 * Level-of-detail
3222 */
3223
3224 static void set_lod_max(int lm) {
3225 lod_max = lm;
3226 }
3227
3228 /*
3229 * Set the scale factor
3230 */
3231 static void set_scale(ale_pos s) {
3232 scale_factor = s;
3233 }
3234
3235 /*
3236 * Set reference rendering to align against
3237 */
3238 static void set_reference(render *r) {
3239 reference = r;
3240 }
3241
3242 /*
3243 * Set the interpolant
3244 */
3245 static void set_interpolant(filter::scaled_filter *f) {
3246 interpolant = f;
3247 }
3248
3249 /*
3250 * Set alignment weights image
3251 */
3252 static void set_weight_map(const image *i) {
3253 weight_map = i;
3254 }
3255
3256 /*
3257 * Set frequency cuts
3258 */
3259 static void set_frequency_cut(double h, double v, double a) {
3260 horiz_freq_cut = h;
3261 vert_freq_cut = v;
3262 avg_freq_cut = a;
3263 }
3264
3265 /*
3266 * Set algorithmic alignment weighting
3267 */
3268 static void set_wmx(const char *e, const char *f, const char *d) {
3269 wmx_exec = e;
3270 wmx_file = f;
3271 wmx_defs = d;
3272 }
3273
3274 /*
3275 * Show frequency weights
3276 */
3277 static void set_fl_show(const char *filename) {
3278 fw_output = filename;
3279 }
3280
3281 /*
3282 * Set transformation file loader.
3283 */
3284 static void set_tload(tload_t *tl) {
3285 tload = tl;
3286 }
3287
3288 /*
3289 * Set transformation file saver.
3290 */
3291 static void set_tsave(tsave_t *ts) {
3292 tsave = ts;
3293 }
3294
3295 /*
3296 * Get match statistics for frame N.
3297 */
3298 static int match(int n) {
3299 update_to(n);
3300
3301 if (n == latest)
3302 return latest_ok;
3303 else if (_keep)
3304 return kept_ok[n];
3305 else {
3306 assert(0);
3307 exit(1);
3308 }
3309 }
3310
3311 /*
3312 * Message that old alignment data should be kept.
3313 */
3314 static void keep() {
3315 assert (latest == -1);
3316 _keep = 1;
3317 }
3318
3319 /*
3320 * Get alignment for frame N.
3321 */
3322 static transformation of(int n) {
3323 update_to(n);
3324 if (n == latest)
3325 return latest_t;
3326 else if (_keep)
3327 return kept_t[n];
3328 else {
3329 assert(0);
3330 exit(1);
3331 }
3332 }
3333
3334 /*
3335 * Set calculation precision
3336 */
3337 static void set_precise(int p) {
3338 precise_calculation = p;
3339 }
3340
3341 /*
3342 * Set the certainty-weighted flag.
3343 */
3344 static void certainty_weighted(int flag) {
3345 certainty_weights = flag;
3346 }
3347
3348 /*
3349 * Set the global search type.
3350 */
3351 static void gs(const char *type) {
3352 if (!strcmp(type, "local")) {
3353 _gs = 0;
3354 } else if (!strcmp(type, "inner")) {
3355 _gs = 1;
3356 } else if (!strcmp(type, "outer")) {
3357 _gs = 2;
3358 } else if (!strcmp(type, "all")) {
3359 _gs = 3;
3360 } else if (!strcmp(type, "central")) {
3361 _gs = 4;
3362 } else if (!strcmp(type, "defaults")) {
3363 _gs = 6;
3364 } else if (!strcmp(type, "points")) {
3365 _gs = 5;
3366 keep();
3367 } else {
3368 ui::get()->error("bad global search type");
3369 }
3370 }
3371
3372 /*
3373 * Multi-alignment contiguity
3374 */
3375 static void ma_cont(double value) {
3376 _ma_cont = value;
3377 }
3378
3379 /*
3380 * Multi-alignment cardinality
3381 */
3382 static void ma_card(unsigned int value) {
3383 assert (value >= 1);
3384 _ma_card = value;
3385 }
3386
3387 /*
3388 * Set the minimum overlap for global searching
3389 */
3390 static void gs_mo(ale_pos value, int _gs_mo_percent) {
3391 _gs_mo = value;
3392 gs_mo_percent = _gs_mo_percent;
3393 }
3394
3395 /*
3396 * Set alignment exclusion regions
3397 */
3398 static void set_exclusion(exclusion *_ax_parameters, int _ax_count) {
3399 ax_count = _ax_count;
3400 ax_parameters = _ax_parameters;
3401 }
3402
3403 /*
3404 * Get match summary statistics.
3405 */
3406 static ale_accum match_summary() {
3407 return match_sum / match_count;
3408 }
3409 };
3410
3411 #endif
synthetic capture engine and renderer