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