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