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