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