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