1 // Copyright 2002, 2004 David Hilvert <dhilvert@auricle.dyndns.org>,
2 // <dhilvert@ugcs.caltech.edu>
4 /* This file is part of the Anti-Lamenessing Engine.
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.
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.
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
22 * align.h: Handle alignment of frames.
29 #include "transformation.h"
40 * Private data members
43 static ale_pos scale_factor
;
46 * Original frame transformation
48 static transformation orig_t
;
51 * Keep data older than latest
54 static transformation
*kept_t
;
58 * Transformation file handlers
61 static tload_t
*tload
;
62 static tsave_t
*tsave
;
65 * Control point variables
68 static const point
**cp_array
;
69 static unsigned int cp_count
;
72 * Reference rendering to align against
75 static render
*reference
;
76 static filter::scaled_filter
*interpolant
;
77 static const image
*reference_image
;
78 static const image
*reference_defined
;
81 * Per-pixel alignment weight map
84 static const image
*weight_map
;
87 * Frequency-dependent alignment weights
90 static double horiz_freq_cut
;
91 static double vert_freq_cut
;
92 static double avg_freq_cut
;
93 static const char *fw_output
;
96 * Algorithmic alignment weighting
99 static const char *wmx_exec
;
100 static const char *wmx_file
;
101 static const char *wmx_defs
;
104 * Consolidated alignment weights
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.
112 static image
*alignment_weights
;
113 static const image
*alignment_weights_const
;
116 * Latest transformation.
119 static transformation latest_t
;
122 * Flag indicating whether the latest transformation
123 * resulted in a match.
126 static int latest_ok
;
129 * Frame number most recently aligned.
135 * Exposure registration
137 * 0. Preserve the original exposure of images.
139 * 1. Match exposure between images.
141 * 2. Use only image metadata for registering exposure.
144 static int _exp_register
;
149 * 0. Translation only. Only adjust the x and y position of images.
150 * Do not rotate input images or perform projective transformations.
152 * 1. Euclidean transformations only. Adjust the x and y position
153 * of images and the orientation of the image about the image center.
155 * 2. Perform general projective transformations. See the file gpt.h
156 * for more information about general projective transformations.
159 static int alignment_class
;
162 * Default initial alignment.
164 * 0. Identity transformation.
166 * 1. Most recently accepted frame's final transformation.
169 static int default_initial_alignment_type
;
170 static transformation default_initial_alignment
;
173 * Projective group behavior
175 * 0. Perturb in output coordinates.
177 * 1. Perturb in source coordinates
180 static int perturb_type
;
183 * Helper variables for new --follow semantics (as of 0.5.0).
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.
192 static int is_default
, old_is_default
;
194 static transformation old_initial_alignment
;
195 static transformation old_final_alignment
;
198 * Alignment for failed frames -- default or optimal?
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.
204 static int is_fail_default
;
209 * 0. Align images with an error contribution from each color channel.
211 * 1. Align images with an error contribution only from the green channel.
212 * Other color channels do not affect alignment.
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.
218 static int channel_alignment_type
;
221 * Error metric exponent
224 static float metric_exponent
;
230 static float match_threshold
;
233 * Perturbation lower and upper bounds.
236 static ale_pos perturb_lower
;
237 static int perturb_lower_percent
;
238 static ale_pos perturb_upper
;
239 static int perturb_upper_percent
;
242 * Maximum level-of-detail scale factor is 2^lod_max/perturb.
248 * Maximum rotational perturbation
251 static ale_pos rot_max
;
254 * Barrel distortion alignment multiplier
257 static ale_pos bda_mult
;
260 * Barrel distortion maximum adjustment rate
263 static ale_pos bda_rate
;
266 * Alignment match sum
269 static ale_accum match_sum
;
272 * Alignment match count.
275 static int match_count
;
278 * Monte Carlo parameter
280 * 0. Don't use monte carlo alignment sampling.
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.
287 static int _mcd_limit
;
290 * Certainty weight flag
292 * 0. Don't use certainty weights for alignment.
294 * 1. Use certainty weights for alignment.
297 static int certainty_weights
;
300 * Global search parameter
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.
313 * Minimum overlap for global searches
316 static unsigned int _gs_mo
;
322 static exclusion
*ax_parameters
;
326 * Type for scale cluster.
329 struct scale_cluster
{
331 const image
*defined
;
332 const image
*aweight
;
333 exclusion
*ax_parameters
;
340 * Check for exclusion region coverage in the reference
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])
356 * Check for exclusion region coverage in the input
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])
372 * Overlap function. Determines the number of pixels in areas where
373 * the arrays overlap. Uses the reference array's notion of pixel
376 static unsigned int overlap(struct scale_cluster c
, transformation t
, int ax_count
) {
377 assert (reference_image
);
379 unsigned int result
= 0;
381 point offset
= c
.accum
->offset();
383 for (unsigned int i
= 0; i
< c
.accum
->height(); i
++)
384 for (unsigned int j
= 0; j
< c
.accum
->width(); j
++) {
386 if (ref_excluded(i
, j
, offset
, c
.ax_parameters
, ax_count
))
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]));
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.
414 if (input_excluded(ti
, tj
, c
.ax_parameters
, ax_count
))
418 && ti
<= c
.input
->height() - 1
420 && tj
<= c
.input
->width() - 1
421 && c
.defined
->get_pixel(i
, j
)[0] != 0)
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.
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
436 * a) Extending the boundaries of the image, or
438 * b) following the previous frame's transform
440 * If we are doing monte-carlo pixel sampling for alignment, we
441 * typically sample a subset of available pixels; otherwise, we sample
451 ale_accum centroid
[2], centroid_divisor
;
453 typedef unsigned int hist_bin
;
459 hist_bin
*histogram_integral
;
463 void add_hist(int r
, int d
, int count
) {
467 int r_shift
= 0, d_shift
= 0;
469 if (r
- hist_min_r
>= hist_dim
) {
470 r_shift
= (r
- hist_min_r
) - hist_dim
+ 1;
471 hist_min_r
+= r_shift
;
474 if (d
- hist_min_d
>= hist_dim
) {
475 d_shift
= (d
- hist_min_d
) - hist_dim
+ 1;
476 hist_min_d
+= d_shift
;
479 assert (r_shift
>= 0);
480 assert (d_shift
>= 0);
482 if (r_shift
|| d_shift
) {
483 for (int rr
= 0; rr
< hist_dim
; rr
++)
484 for (int dd
= 0; dd
< hist_dim
; dd
++) {
486 hist_bin value
= histogram
[rr
* hist_dim
+ dd
];
488 histogram
[rr
* hist_dim
+ dd
] = 0;
490 int rrr
= rr
- r_shift
;
491 int ddd
= dd
- d_shift
;
498 histogram
[rrr
* hist_dim
+ ddd
] += value
;
510 histogram
[r
* hist_dim
+ d
] += count
;
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);
517 if (!(rbin
> INT_MIN
))
519 if (!(dbin
> INT_MIN
))
522 add_hist((int) floor(rbin
), (int) floor(dbin
), 1);
529 hist_min_r
= INT_MIN
;
530 hist_min_d
= INT_MIN
;
533 histogram
= (hist_bin
*) calloc(hist_dim
* hist_dim
, sizeof(hist_bin
));
535 min
= point::posinf();
536 max
= point::neginf();
540 centroid_divisor
= 0;
548 hist_min_r
= INT_MIN
;
549 hist_min_d
= INT_MIN
;
551 histogram
= (hist_bin
*) calloc(hist_dim
* hist_dim
, sizeof(hist_bin
));
558 int check_removal(diff_stat_t
*with
) {
559 ale_accum bresult
, bdivisor
, wresult
, wdivisor
;
560 hist_bin
*bhist
, *whist
;
562 ui::get()->d2_align_sim_start();
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
);
569 wresult
= with
->result
;
570 wdivisor
= with
->divisor
;
572 for (int i
= 0; i
< hist_dim
* hist_dim
; i
++)
573 bhist
[i
] = histogram
[i
];
575 for (int i
= 0; i
< with
->hist_dim
* with
->hist_dim
; i
++)
576 whist
[i
] = with
->histogram
[i
];
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;
583 for (int i
= 0; i
< hist_dim
* hist_dim
; i
++) {
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
;
592 if (b_test_gradient
> max_gradient
) {
593 max_gradient_bin
= i
;
594 max_gradient_hist
= 0;
595 max_gradient
= b_test_gradient
;
599 for (int i
= 0; i
< with
->hist_dim
* with
->hist_dim
; i
++) {
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
);
608 if (w_test_gradient
> max_gradient
) {
609 max_gradient_bin
= i
;
610 max_gradient_hist
= 1;
611 max_gradient
= w_test_gradient
;
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
);
629 ui::get()->d2_align_sim_stop();
631 if (bresult
/ bdivisor
< wresult
/ wdivisor
)
637 int reliable(diff_stat_t
*with
, ale_pos mc
) {
638 return check_removal(with
);
641 void add(const diff_stat_t
*ds
) {
642 result
+= ds
->result
;
643 divisor
+= ds
->divisor
;
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
]);
650 for (int d
= 0; d
< 2; d
++) {
651 if (min
[d
] > ds
->min
[d
])
653 if (max
[d
] < ds
->max
[d
])
655 centroid
[d
] += ds
->centroid
[d
];
658 centroid_divisor
+= ds
->centroid_divisor
;
661 ale_accum
get_result() {
665 ale_accum
get_divisor() {
669 point
get_centroid() {
670 return point(centroid
[0] / centroid_divisor
, centroid
[1] / centroid_divisor
);
673 ale_accum
get_error() {
674 return pow(result
/ divisor
, 1/metric_exponent
);
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
);
681 ale_accum this_result
= 0;
682 ale_accum this_divisor
= 0;
685 if (interpolant
!= NULL
)
686 interpolant
->filtered(i
, j
, &pb
, &weight
, 1, f
);
689 c
.input
->get_bl(point(ti
, tj
), result
);
698 if (certainty_weights
== 0)
699 weight
= pixel(1, 1, 1);
701 if (c
.aweight
!= NULL
)
702 weight
*= c
.aweight
->get_pixel(i
, j
);
705 * Update sampling area statistics
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]);
722 * Determine alignment type.
725 if (channel_alignment_type
== 0) {
727 * Align based on all channels.
731 for (int k
= 0; k
< 3; k
++) {
732 ale_real achan
= pa
[k
];
733 ale_real bchan
= pb
[k
];
735 this_result
+= weight
[k
] * pow(fabs(achan
- bchan
), metric_exponent
);
736 this_divisor
+= weight
[k
] * pow(achan
> bchan
? achan
: bchan
, metric_exponent
);
738 } else if (channel_alignment_type
== 1) {
740 * Align based on the green channel.
743 ale_real achan
= pa
[1];
744 ale_real bchan
= pb
[1];
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) {
750 * Align based on the sum of all channels.
757 for (int k
= 0; k
< 3; k
++) {
760 wsum
+= weight
[k
] / 3;
763 this_result
= wsum
* pow(fabs(asum
- bsum
), metric_exponent
);
764 this_divisor
= wsum
* pow(asum
> bsum
? asum
: bsum
, metric_exponent
);
767 result
+= this_result
;
768 divisor
+= this_divisor
;
770 add_hist(this_result
, this_divisor
);
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");
781 hist_bin recalc_total
= 0;
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
]);
788 fprintf(stderr
, "\n");
791 fprintf(stderr
, "\n");
792 fprintf(stderr
, "recalc_total = %d\n", recalc_total
);
796 struct subdomain_args
{
797 struct scale_cluster c
;
802 diff_stat_t
* diff_stat
;
803 int i_min
, i_max
, j_min
, j_max
;
807 static void *diff_subdomain(void *args
) {
809 subdomain_args
*sargs
= (subdomain_args
*) args
;
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
;
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
;
822 assert (reference_image
);
824 point offset
= c
.accum
->offset();
829 * We always use the same code for exhaustive and Monte Carlo
830 * pixel sampling, setting _mc_arg = 1 when all pixels are to
834 if (_mc_arg
<= 0 || _mc_arg
>= 1)
839 int index_max
= (i_max
- i_min
) * (j_max
- j_min
);
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:
849 * We have mc_arg == USE/ALL, or (expected # pixels to
850 * use)/(# total pixels). We derive from this
853 * SKIP/USE == (SKIP/ALL)/(USE/ALL) == (1 - (USE/ALL))/(USE/ALL)
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.
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
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.
871 ale_pos u
= (1 - _mc_arg
) / _mc_arg
;
873 ale_pos mc_max
= (floor(2*u
) * (1 + floor(2*u
)) + 2*u
)
874 / (2 + 2 * floor(2*u
) - 2*u
);
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:
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
892 * For subdomain calculations, we vary the seed by subdomain.
897 rng
.seed(1 + subdomain
);
899 for(index
= -1 + (int) ceil((mc_max
+1)
900 * ( (1 + ((ale_pos
) (rng
.get())) )
901 / (1 + ((ale_pos
) RAND_MAX
)) ));
903 index
+= (int) ceil((mc_max
+1)
904 * ( (1 + ((ale_pos
) (rng
.get())) )
905 / (1 + ((ale_pos
) RAND_MAX
)) ))){
907 i
= index
/ (j_max
- j_min
) + i_min
;
908 j
= index
% (j_max
- j_min
) + j_min
;
911 * Check for exclusion in render coordinates.
914 if (ref_excluded(i
, j
, offset
, c
.ax_parameters
, ax_count
))
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]));
933 * Check that the transformed coordinates are within
934 * the boundaries of array c.input and that they
935 * are not subject to exclusion.
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).
943 if (input_excluded(ti
, tj
, c
.ax_parameters
, ax_count
))
947 && ti
<= c
.input
->height() - 1
949 && tj
<= c
.input
->width() - 1
950 && c
.defined
->get_pixel(i
, j
)[0] != 0)
952 sargs
->diff_stat
->sample(f
, c
, i
, j
, ti
, tj
);
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
) {
962 diff_stat_t
*diff_stat
= diff_stat_in
;
964 if (diff_stat
== NULL
)
965 diff_stat
= new diff_stat_t();
969 ui::get()->d2_align_sample_start();
971 if (interpolant
!= NULL
)
972 interpolant
->set_parameters(t
, c
.input
, c
.accum
->offset());
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
);
985 subdomain_args
*args
= new subdomain_args
[N
];
987 for (int ti
= 0; ti
< N
; ti
++) {
990 args
[ti
]._mc_arg
= _mc_arg
;
991 args
[ti
].ax_count
= ax_count
;
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
;
997 args
[ti
].j_max
= c
.accum
->width();
998 args
[ti
].subdomain
= ti
;
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
]);
1005 diff_subdomain(&args
[ti
]);
1009 for (int ti
= 0; ti
< N
; ti
++) {
1011 pthread_join(threads
[ti
], NULL
);
1013 diff_stat
->add(args
[ti
].diff_stat
);
1015 delete args
[ti
].diff_stat
;
1020 ui::get()->d2_align_sample_stop();
1022 ale_accum result
= diff_stat
->get_error();
1024 if (diff_stat_in
== NULL
)
1032 * Adjust exposure for an aligned frame B against reference A.
1034 * Expects full-LOD images.
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.
1041 * [Note that this may have been due to a bug in certainty determination
1042 * within this function.]
1044 static void set_exposure_ratio(unsigned int m
, struct scale_cluster c
,
1045 transformation t
, int ax_count
, int pass_number
) {
1047 if (_exp_register
== 2) {
1049 * Use metadata only.
1051 ale_real gain_multiplier
= image_rw::exp(m
).get_gain_multiplier();
1052 pixel multiplier
= pixel(gain_multiplier
, gain_multiplier
, gain_multiplier
);
1054 image_rw::exp(m
).set_multiplier(multiplier
);
1055 ui::get()->exp_multiplier(multiplier
[0],
1062 pixel_accum asum
, bsum
;
1064 point offset
= c
.accum
->offset();
1066 for (unsigned int i
= 0; i
< c
.accum
->height(); i
++)
1067 for (unsigned int j
= 0; j
< c
.accum
->width(); j
++) {
1069 if (ref_excluded(i
, j
, offset
, c
.ax_parameters
, ax_count
))
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]));
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.
1091 if (input_excluded(q
[0], q
[1], c
.ax_parameters
, ax_count
))
1095 && q
[0] <= c
.input
->height() - 1
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
);
1102 c
.input
->get_bl(q
, b
);
1105 pixel weight
= (c
.aweight
1106 ? c
.aweight
->get_pixel(i
, j
)
1108 * ((!certainty_weights
&& pass_number
)
1109 ? c
.defined
->get_pixel(i
, j
)
1115 pixel weight
= pixel(1, 1, 1);
1119 bsum
+= b
[0] * weight
;
1123 // std::cerr << (asum / bsum) << " ";
1125 pixel_accum new_multiplier
;
1127 new_multiplier
= asum
/ bsum
* image_rw::exp(m
).get_multiplier();
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],
1140 * Copy all ax parameters.
1142 static exclusion
*copy_ax_parameters(int local_ax_count
, exclusion
*source
) {
1144 exclusion
*dest
= (exclusion
*) malloc(local_ax_count
* sizeof(exclusion
));
1149 ui::get()->memory_error("exclusion regions");
1151 for (int idx
= 0; idx
< local_ax_count
; idx
++)
1152 dest
[idx
] = source
[idx
];
1158 * Copy ax parameters according to frame.
1160 static exclusion
*filter_ax_parameters(int frame
, int *local_ax_count
) {
1162 exclusion
*dest
= (exclusion
*) malloc(ax_count
* sizeof(exclusion
));
1167 ui::get()->memory_error("exclusion regions");
1169 *local_ax_count
= 0;
1171 for (int idx
= 0; idx
< ax_count
; idx
++) {
1172 if (ax_parameters
[idx
].x
[4] > frame
1173 || ax_parameters
[idx
].x
[5] < frame
)
1176 dest
[*local_ax_count
] = ax_parameters
[idx
];
1178 (*local_ax_count
)++;
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
)
1191 for (int n
= 0; n
< 6; n
++) {
1192 ax_parameters
[i
].x
[n
] = (int) floor(ax_parameters
[i
].x
[n
]
1199 * Prepare the next level of detail for ordinary images.
1201 static const image
*prepare_lod(const image
*current
) {
1202 if (current
== NULL
)
1205 return current
->scale_by_half("prepare_lod");
1209 * Prepare the next level of detail for weighted images.
1211 static const image
*prepare_lod(const image
*current
, const image
*weights
) {
1212 if (current
== NULL
)
1215 return current
->scale_by_half(weights
, "prepare_lod");
1219 * Prepare the next level of detail for definition maps.
1221 static const image
*prepare_lod_def(const image
*current
) {
1224 return current
->defined_scale_by_half("prepare_lod_def");
1228 * Initialize the scale cluster data structure.
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
) {
1235 * Allocate memory for the array.
1238 struct scale_cluster
*scale_clusters
=
1239 (struct scale_cluster
*) malloc(steps
* sizeof(struct scale_cluster
));
1241 assert (scale_clusters
);
1243 if (!scale_clusters
)
1244 ui::get()->memory_error("alignment");
1247 * Prepare images and exclusion regions for the highest level
1251 scale_clusters
[0].accum
= reference_image
;
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");
1258 scale_clusters
[0].input
= input_frame
;
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
);
1264 scale_ax_parameters(*local_ax_count
, scale_clusters
[0].ax_parameters
, scale_factor
,
1265 (scale_factor
< 1.0 && interpolant
== NULL
) ? scale_factor
: 1);
1268 * Prepare reduced-detail images and exclusion
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
);
1280 double sf
= scale_clusters
[step
- 1].input_scale
/ 2;
1281 scale_clusters
[step
].input_scale
= sf
;
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
);
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);
1295 return scale_clusters
;
1299 * Destroy the first element in the scale cluster data structure.
1301 static void final_clusters(struct scale_cluster
*scale_clusters
, ale_real scale_factor
,
1302 unsigned int steps
) {
1304 if (scale_clusters
[0].input_scale
< 1.0)
1305 delete scale_clusters
[0].input
;
1307 free((void *)scale_clusters
[0].ax_parameters
);
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
;
1314 if (scale_clusters
[step
].input_scale
< 1.0)
1315 delete scale_clusters
[step
].input
;
1317 free((void *)scale_clusters
[step
].ax_parameters
);
1320 free(scale_clusters
);
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.
1327 static point
unscaled_centroid(unsigned int m
, unsigned int p
) {
1330 point
point_sum(0, 0);
1331 ale_accum divisor
= 0;
1333 for(unsigned int j
= 0; j
< m
; j
++) {
1334 point pp
= cp_array
[p
][j
];
1337 point_sum
+= kept_t
[j
].transform_unscaled(pp
)
1338 / kept_t
[j
].scale();
1344 return point::undefined();
1346 return point_sum
/ divisor
;
1350 * Calculate centroid of this frame, and of all previous frames,
1351 * from points common to both sets.
1353 static void centroids(unsigned int m
, point
*current
, point
*previous
) {
1355 * Calculate the translation
1357 point
other_centroid(0, 0);
1358 point
this_centroid(0, 0);
1359 ale_pos divisor
= 0;
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
];
1365 if (!other_c
.defined() || !this_c
.defined())
1368 other_centroid
+= other_c
;
1369 this_centroid
+= this_c
;
1375 *current
= point::undefined();
1376 *previous
= point::undefined();
1380 *current
= this_centroid
/ divisor
;
1381 *previous
= other_centroid
/ divisor
;
1385 * Calculate the RMS error of control points for frame m, with
1386 * transformation t, against control points for earlier frames.
1388 static ale_accum
cp_rms_error(unsigned int m
, transformation t
) {
1392 ale_accum divisor
= 0;
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
]);
1400 if (!p_ref
.defined() || !p_cur
.defined())
1403 err
+= p_ref
.lengthtosq(p_cur
);
1407 return sqrt(err
/ divisor
);
1411 * Align frame m against the reference.
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).
1419 static ale_accum
_align(int m
, int local_gs
) {
1422 * Open the input frame.
1426 ui::get()->loading_file();
1427 const image
*input_frame
= image_rw::open(m
);
1430 * Local upper/lower data, possibly dependent on image
1434 ale_pos local_lower
, local_upper
;
1437 * Select the minimum dimension as the reference.
1440 ale_pos reference_size
= input_frame
->height();
1441 if (input_frame
->width() < reference_size
)
1442 reference_size
= input_frame
->width();
1444 if (perturb_lower_percent
)
1445 local_lower
= perturb_lower
1450 local_lower
= perturb_lower
;
1452 if (perturb_upper_percent
)
1453 local_upper
= perturb_upper
1458 local_upper
= perturb_upper
;
1460 local_upper
= pow(2, floor(log(local_upper
) / log(2)));
1463 * Logarithms aren't exact, so we divide repeatedly to discover
1464 * how many steps will occur, and pass this information to the
1469 double step_variable
= local_upper
;
1470 while (step_variable
>= local_lower
) {
1475 ui::get()->set_steps(step_count
);
1477 ale_pos perturb
= local_upper
;
1478 transformation offset
;
1480 diff_stat_t
*here_diff_stat
= new diff_stat_t();
1484 kept_t
[latest
] = latest_t
;
1485 kept_ok
[latest
] = latest_ok
;
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.
1494 const int lod_diff
= lod_max
;
1497 * Determine how many levels of detail should be prepared.
1500 unsigned int steps
= (perturb
> pow(2, lod_max
))
1501 ? (unsigned int) (log(perturb
) / log(2)) - lod_max
+ 1 : 1;
1505 * Prepare multiple levels of detail.
1509 struct scale_cluster
*scale_clusters
= init_clusters(m
,
1510 scale_factor
, input_frame
, steps
,
1514 * Initialize variables used in the main loop.
1520 * Initialize the default initial transform
1523 if (default_initial_alignment_type
== 0) {
1526 * Follow the transformation of the original frame,
1527 * setting new image dimensions.
1530 default_initial_alignment
= orig_t
;
1531 default_initial_alignment
.set_dimensions(input_frame
);
1533 } else if (default_initial_alignment_type
== 1)
1536 * Follow previous transformation, setting new image
1540 default_initial_alignment
.set_dimensions(input_frame
);
1545 old_is_default
= is_default
;
1548 * Scale default initial transform for lod
1551 default_initial_alignment
.rescale(1 / pow(2, lod
));
1554 * Load any file-specified transformation
1557 offset
= tload_next(tload
, alignment_class
== 2, default_initial_alignment
, &is_default
);
1559 transformation new_offset
= offset
;
1561 if (!old_is_default
&& !is_default
&& default_initial_alignment_type
== 1) {
1564 * Implement new delta --follow semantics.
1566 * If we have a transformation T such that
1568 * prev_final == T(prev_init)
1572 * current_init_follow == T(current_init)
1574 * We can calculate T as follows:
1576 * T == prev_final(prev_init^-1)
1578 * Where ^-1 is the inverse operator.
1582 * Ensure that the lod for the old initial and final
1583 * alignments are equal to the lod for the new initial
1587 old_final_alignment
.rescale (1 / pow(2, lod
));
1588 old_initial_alignment
.rescale(1 / pow(2, lod
- old_lod
));
1590 if (alignment_class
== 0) {
1592 * Translational transformations
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);
1598 new_offset
.eu_modify(0, t0
);
1599 new_offset
.eu_modify(1, t1
);
1601 } else if (alignment_class
== 1) {
1603 * Euclidean transformations
1606 ale_pos t2
= -old_initial_alignment
.eu_get(2) + old_final_alignment
.eu_get(2);
1608 new_offset
.eu_modify(2, t2
);
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) );
1613 p
= old_final_alignment
.transform_scaled(p
);
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));
1618 } else if (alignment_class
== 2) {
1620 * Projective transformations
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()))));
1634 new_offset
.gpt_set(p
);
1638 old_initial_alignment
= offset
;
1640 offset
= new_offset
;
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())); */
1647 ui::get()->alignment_monte_carlo_parameter(_mc_arg
);
1650 * Projective adjustment value
1653 ale_pos adj_p
= (perturb
>= pow(2, lod_diff
))
1654 ? pow(2, lod_diff
) : (double) perturb
;
1657 * Pre-alignment exposure adjustment
1660 if (_exp_register
) {
1661 ui::get()->exposure_1();
1662 transformation o
= offset
;
1663 for (int k
= lod
; k
> 0; k
--)
1665 set_exposure_ratio(m
, scale_clusters
[0], o
, local_ax_count
, 0);
1669 * Current difference (error) value
1672 ui::get()->prematching();
1673 here
= diff(si
, offset
, _mc_arg
, local_ax_count
, m
, here_diff_stat
);
1674 ui::get()->set_match(here
);
1677 * Current and modified barrel distortion parameters
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
);
1686 * Translational global search step
1689 if (perturb
>= local_lower
&& local_gs
!= 0 && local_gs
!= 5) {
1691 ui::get()->aligning(perturb
, lod
);
1693 transformation lowest_t
= offset
;
1694 ale_accum lowest_v
= here
;
1698 transformation offset_p
= offset
;
1700 if (!offset_p
.is_projective())
1701 offset_p
.eu_to_gpt();
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
);
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
);
1721 if (local_gs
== 1 || local_gs
== 3 || local_gs
== 4) {
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
);
1734 if ((v
< lowest_v
&& ovl
>= _gs_mo
)
1735 || (!finite(lowest_v
) && finite(v
))) {
1742 if (local_gs
== 2 || local_gs
== 3 || local_gs
== -1) {
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
);
1755 if ((v
< lowest_v
&& ovl
>= _gs_mo
)
1756 || (!finite(lowest_v
) && finite(v
))) {
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
);
1768 if ((v
< lowest_v
&& ovl
>= _gs_mo
)
1769 || (!finite(lowest_v
) && finite(v
))) {
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
);
1781 if ((v
< lowest_v
&& ovl
>= _gs_mo
)
1782 || (!finite(lowest_v
) && finite(v
))) {
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
);
1794 if ((v
< lowest_v
&& ovl
>= _gs_mo
)
1795 || (!finite(lowest_v
) && finite(v
))) {
1805 ui::get()->set_match(here
);
1809 * Control point alignment
1812 if (local_gs
== 5) {
1814 transformation o
= offset
;
1815 for (int k
= lod
; k
> 0; k
--)
1819 * Determine centroid data
1822 point current
, previous
;
1823 centroids(m
, ¤t
, &previous
);
1825 if (current
.defined() && previous
.defined()) {
1827 o
.set_dimensions(input_frame
);
1828 o
.translate((previous
- current
) * o
.scale());
1833 * Determine rotation for alignment classes other than translation.
1836 ale_accum lowest_error
= cp_rms_error(m
, o
);
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))
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
) {
1850 transformation test_t
= o
;
1851 test_t
.rotate(current
* o
.scale(), srot
);
1852 ale_pos test_v
= cp_rms_error(m
, test_t
);
1854 if (test_v
< lowest_error
) {
1855 lowest_error
= test_v
;
1864 * Determine projective parameters through a local
1868 if (alignment_class
== 2) {
1869 ale_accum adj_p
= lowest_error
;
1871 if (adj_p
< local_lower
)
1872 adj_p
= local_lower
;
1874 while (adj_p
>= local_lower
) {
1875 transformation test_t
= o
;
1876 int is_improved
= 1;
1880 while (is_improved
) {
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
) {
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
);
1896 test_v
= cp_rms_error(m
, test_t
);
1898 if (test_v
< lowest_error
) {
1899 lowest_error
= test_v
;
1911 set_exposure_ratio(m
, scale_clusters
[0], o
, local_ax_count
, 0);
1913 for (int k
= lod
; k
> 0; k
--)
1920 * Perturbation adjustment loop.
1923 while (perturb
>= local_lower
) {
1925 ui::get()->aligning(perturb
, lod
);
1930 * Orientational adjustment value in degrees.
1932 * Since rotational perturbation is now specified as an
1933 * arclength, we have to convert.
1936 ale_pos adj_o
= 2 * perturb
1937 / sqrt(pow(scale_clusters
[0].input
->height(), 2)
1938 + pow(scale_clusters
[0].input
->width(), 2))
1943 * Barrel distortion adjustment value
1946 ale_pos adj_b
= perturb
* bda_mult
;
1948 transformation test_t
;
1949 transformation old_offset
= offset
;
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;
1959 if (alignment_class
< 2 && alignment_class
>= 0) {
1962 * Translational or euclidean transformation
1965 for (int i
= 0; i
< 2; i
++)
1966 for (adj_s
= -adj_p
; adj_s
<= adj_p
; adj_s
+= 2 * adj_p
) {
1970 test_t
.eu_modify(i
, adj_s
);
1972 test_d
= diff(si
, test_t
, _mc_arg
, local_ax_count
, m
, here_diff_stat
);
1974 if (test_d
< here
|| (!finite(here
) && finite(test_d
))) {
1978 || here_diff_stat
->reliable(old_here_diff_stat
, _mc_arg
)) {
1979 found_reliable_better
= 1;
1986 if (_mc
<= 0 && test_d
> here
) {
1988 || old_here_diff_stat
->reliable(here_diff_stat
, _mc_arg
))
1989 found_reliable_worse
= 1;
1991 found_unreliable_worse
= 1;
1995 if (alignment_class
== 1 && adj_o
< rot_max
)
1996 for (adj_s
= -adj_o
; adj_s
<= adj_o
; adj_s
+= 2 * adj_o
) {
1998 point sample_centroid
= old_here_diff_stat
->get_centroid() + si
.accum
->offset();
2002 // test_t.eu_modify(2, adj_s);
2005 test_t
.eu_rotate_about_scaled(sample_centroid
, adj_s
);
2007 test_d
= diff(si
, test_t
, _mc_arg
, local_ax_count
, m
, here_diff_stat
);
2009 if (test_d
< here
|| (!finite(here
) && finite(test_d
))) {
2013 || here_diff_stat
->reliable(old_here_diff_stat
, _mc_arg
)) {
2014 found_reliable_better
= 1;
2021 if (_mc
<= 0 && test_d
> here
) {
2023 || old_here_diff_stat
->reliable(here_diff_stat
, _mc_arg
))
2024 found_reliable_worse
= 1;
2026 found_unreliable_worse
= 1;
2030 } else if (alignment_class
== 2) {
2033 * Projective transformation
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
) {
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
);
2049 test_d
= diff(si
, test_t
, _mc_arg
, local_ax_count
, m
, here_diff_stat
);
2052 fprintf(stderr
, "old\n");
2053 old_here_diff_stat
->print_hist();
2054 fprintf(stderr
, "new\n");
2055 here_diff_stat
->print_hist();
2058 if (test_d
< here
|| (!finite(here
) && finite(test_d
))) {
2059 // fprintf(stderr, "found better\n");
2063 || here_diff_stat
->reliable(old_here_diff_stat
, _mc_arg
)) {
2064 // fprintf(stderr, "found reliable better\n");
2065 found_reliable_better
= 1;
2072 if (_mc
<= 0 && test_d
> here
) {
2074 || old_here_diff_stat
->reliable(here_diff_stat
, _mc_arg
))
2075 found_reliable_worse
= 1;
2077 found_unreliable_worse
= 1;
2087 if (bda_mult
!= 0) {
2089 static unsigned int d_rotation
= 0;
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
) {
2094 unsigned int rd
= (d
+ d_rotation
) % offset
.bd_count();
2095 d_rotation
= (d_rotation
+ 1) % offset
.bd_count();
2097 if (bda_rate
> 0 && fabs(modified_bd
[rd
] + adj_s
- current_bd
[rd
]) > bda_rate
)
2102 test_t
.bd_modify(rd
, adj_s
);
2104 test_d
= diff(si
, test_t
, _mc_arg
, local_ax_count
, m
, here_diff_stat
);
2106 if (test_d
< here
|| (!finite(here
) && finite(test_d
))) {
2110 || here_diff_stat
->reliable(old_here_diff_stat
, _mc_arg
)) {
2111 found_reliable_better
= 1;
2114 modified_bd
[rd
] += adj_s
;
2119 if (_mc
<= 0 && test_d
> here
) {
2121 || old_here_diff_stat
->reliable(here_diff_stat
, _mc_arg
))
2122 found_reliable_worse
= 1;
2124 found_unreliable_worse
= 1;
2131 if (_mc_arg
< 1 && _mc
<= 0
2132 && (!found_reliable_worse
|| found_better
)
2133 /* && found_unreliable_worse */
2134 && !found_reliable_better
) {
2136 ui::get()->alignment_monte_carlo_parameter(_mc_arg
);
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();
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
);
2148 if (here_diff_stat_half
->reliable(old_here_diff_stat_half
, _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
;
2156 delete here_diff_stat_half
;
2157 delete old_here_diff_stat_half
;
2161 if (!(here
< old_here
) && !(!finite(old_here
) && finite(here
))) {
2170 * We're about to work with images
2171 * twice as large, so rescale the
2176 default_initial_alignment
.rescale(2);
2179 * Work with images twice as large
2183 si
= scale_clusters
[lod
];
2188 here
= diff(si
, offset
, _mc_arg
, local_ax_count
, m
, here_diff_stat
);
2189 delete old_here_diff_stat
;
2192 delete here_diff_stat
;
2193 here_diff_stat
= old_here_diff_stat
;
2201 ui::get()->alignment_monte_carlo_parameter(_mc_arg
);
2202 ui::get()->alignment_perturbation_level(perturb
, lod
);
2205 delete old_here_diff_stat
;
2208 ui::get()->set_match(here
);
2212 offset
.rescale(pow(2, lod
));
2213 default_initial_alignment
.rescale(pow(2, lod
));
2217 * Post-alignment exposure adjustment
2220 if (_exp_register
== 1) {
2221 ui::get()->exposure_2();
2222 set_exposure_ratio(m
, scale_clusters
[0], offset
, local_ax_count
, 1);
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
);
2236 * Free the level-of-detail structures
2239 final_clusters(scale_clusters
, scale_factor
, steps
);
2242 * Ensure that the match meets the threshold.
2245 if ((1 - here
) * 100 > match_threshold
) {
2247 * Update alignment variables
2250 default_initial_alignment
= offset
;
2251 old_final_alignment
= offset
;
2252 ui::get()->alignment_match_ok();
2253 } else if (local_gs
== 4) {
2256 * Align with outer starting points.
2260 * XXX: This probably isn't exactly the right thing to do,
2261 * since variables like old_initial_value have been overwritten.
2264 ale_accum nested_result
= _align(m
, -1);
2266 if ((1 - nested_result
) * 100 > match_threshold
) {
2267 return nested_result
;
2268 } else if (nested_result
< here
) {
2269 here
= nested_result
;
2273 if (is_fail_default
)
2274 offset
= default_initial_alignment
;
2276 ui::get()->set_match(here
);
2277 ui::get()->alignment_no_match();
2279 } else if (local_gs
== -1) {
2286 if (is_fail_default
)
2287 offset
= default_initial_alignment
;
2289 ui::get()->alignment_no_match();
2293 * Write the tonal registration multiplier as a comment.
2296 pixel trm
= image_rw::exp(m
).get_multiplier();
2297 tsave_trm(tsave
, trm
[0], trm
[1], trm
[2]);
2300 * Save the transformation information
2303 tsave_next(tsave
, offset
, alignment_class
== 2);
2308 * Update match statistics.
2311 match_sum
+= (1 - here
) * 100;
2316 * Close the image file.
2326 * High-pass filter for frequency weights
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;
2347 * Reset alignment weights
2349 static void reset_weights() {
2351 alignment_weights_const
= NULL
;
2353 if (alignment_weights
!= NULL
)
2354 delete alignment_weights
;
2356 alignment_weights
= NULL
;
2360 * Initialize alignment weights
2362 static void init_weights() {
2363 if (alignment_weights
!= NULL
)
2366 int rows
= reference_image
->height();
2367 int cols
= reference_image
->width();
2368 int colors
= reference_image
->depth();
2370 alignment_weights
= new image_ale_real(rows
, cols
,
2371 colors
, "alignment_weights");
2373 assert (alignment_weights
);
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));
2381 * Update alignment weights with weight map
2383 static void map_update() {
2385 if (weight_map
== NULL
)
2390 point map_offset
= reference_image
->offset() - weight_map
->offset();
2392 int rows
= reference_image
->height();
2393 int cols
= reference_image
->width();
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
);
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.
2413 static void imap_update() {
2414 if (alignment_weights
== NULL
) {
2415 alignment_weights_const
= reference_defined
;
2417 int rows
= reference_image
->height();
2418 int cols
= reference_image
->width();
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
);
2424 alignment_weights_const
= alignment_weights
;
2429 * Update alignment weights with algorithmic weights
2431 static void wmx_update() {
2434 static exposure
*exp_def
= new exposure_default();
2435 static exposure
*exp_bool
= new exposure_boolean();
2437 if (wmx_file
== NULL
|| wmx_exec
== NULL
|| wmx_defs
== NULL
)
2440 unsigned int rows
= reference_image
->height();
2441 unsigned int cols
= reference_image
->width();
2443 image_rw::write_image(wmx_file
, reference_image
);
2444 image_rw::write_image(wmx_defs
, reference_defined
, exp_bool
);
2447 int exit_status
= 1;
2449 execlp(wmx_exec
, wmx_exec
, wmx_file
, wmx_defs
, NULL
);
2450 ui::get()->exec_failure(wmx_exec
, wmx_file
, wmx_defs
);
2456 ui::get()->fork_failure("d2::align");
2458 image
*wmx_weights
= image_rw::read_image(wmx_file
, exp_def
);
2460 if (wmx_weights
->height() != rows
|| wmx_weights
->width() != cols
)
2461 ui::get()->error("algorithmic weighting must not change image size");
2463 if (alignment_weights
== NULL
)
2464 alignment_weights
= wmx_weights
;
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
);
2473 * Update alignment weights with frequency weights
2475 static void fw_update() {
2477 if (horiz_freq_cut
== 0
2478 && vert_freq_cut
== 0
2479 && avg_freq_cut
== 0)
2483 * Required for correct operation of --fwshow
2486 assert (alignment_weights
== NULL
);
2488 int rows
= reference_image
->height();
2489 int cols
= reference_image
->width();
2490 int colors
= reference_image
->depth();
2492 alignment_weights
= new image_ale_real(rows
, cols
,
2493 colors
, "alignment_weights");
2495 fftw_complex
*inout
= (fftw_complex
*) fftw_malloc(sizeof(fftw_complex
) * rows
* cols
);
2499 fftw_plan pf
= fftw_plan_dft_2d(rows
, cols
,
2501 FFTW_FORWARD
, FFTW_ESTIMATE
);
2503 fftw_plan pb
= fftw_plan_dft_2d(rows
, cols
,
2505 FFTW_BACKWARD
, FFTW_ESTIMATE
);
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
];
2514 hipass(rows
, cols
, inout
);
2517 for (int i
= 0; i
< rows
* cols
; i
++) {
2519 alignment_weights
->pix(i
/ cols
, i
% cols
)[k
] = fabs(inout
[i
][0] / (rows
* cols
));
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));
2528 fftw_destroy_plan(pf
);
2529 fftw_destroy_plan(pb
);
2532 if (fw_output
!= NULL
)
2533 image_rw::write_image(fw_output
, alignment_weights
);
2538 * Update alignment to frame N.
2540 static void update_to(int n
) {
2541 assert (n
<= latest
+ 1);
2544 const image
*i
= image_rw::open(n
);
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);
2554 kept_t
= (transformation
*) malloc(image_rw::count()
2555 * sizeof(transformation
));
2556 kept_ok
= (int *) malloc(image_rw::count()
2561 if (!kept_t
|| !kept_ok
)
2562 ui::get()->memory_error("alignment");
2572 default_initial_alignment
= result
;
2576 for (int i
= latest
+ 1; i
<= n
; i
++) {
2577 assert (reference
!= NULL
);
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();
2590 if (certainty_weights
)
2591 imap_update(); /* Must be called after all other _updates */
2593 assert (reference_image
!= NULL
);
2594 assert (reference_defined
!= NULL
);
2603 * Set the control point count
2605 static void set_cp_count(unsigned int n
) {
2606 assert (cp_array
== NULL
);
2609 cp_array
= (const point
**) malloc(n
* sizeof(const point
*));
2613 * Set control points.
2615 static void set_cp(unsigned int i
, const point
*p
) {
2622 static void exp_register() {
2627 * Register exposure only based on metadata
2629 static void exp_meta_only() {
2634 * Don't register exposure
2636 static void exp_noregister() {
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.
2645 static void class_translation() {
2646 alignment_class
= 0;
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
2654 static void class_euclidean() {
2655 alignment_class
= 1;
2659 * Set alignment class to perform general projective transformations.
2660 * See the file gpt.h for more information about general projective
2663 static void class_projective() {
2664 alignment_class
= 2;
2668 * Set the default initial alignment to the identity transformation.
2670 static void initial_default_identity() {
2671 default_initial_alignment_type
= 0;
2675 * Set the default initial alignment to the most recently matched
2676 * frame's final transformation.
2678 static void initial_default_follow() {
2679 default_initial_alignment_type
= 1;
2683 * Perturb output coordinates.
2685 static void perturb_output() {
2690 * Perturb source coordinates.
2692 static void perturb_source() {
2697 * Frames under threshold align optimally
2699 static void fail_optimal() {
2700 is_fail_default
= 0;
2704 * Frames under threshold keep their default alignments.
2706 static void fail_default() {
2707 is_fail_default
= 1;
2711 * Align images with an error contribution from each color channel.
2714 channel_alignment_type
= 0;
2718 * Align images with an error contribution only from the green channel.
2719 * Other color channels do not affect alignment.
2721 static void green() {
2722 channel_alignment_type
= 1;
2726 * Align images using a summation of channels. May be useful when
2727 * dealing with images that have high frequency color ripples due to
2731 channel_alignment_type
= 2;
2735 * Error metric exponent
2738 static void set_metric_exponent(float me
) {
2739 metric_exponent
= me
;
2746 static void set_match_threshold(float mt
) {
2747 match_threshold
= mt
;
2751 * Perturbation lower and upper bounds.
2754 static void set_perturb_lower(ale_pos pl
, int plp
) {
2756 perturb_lower_percent
= plp
;
2759 static void set_perturb_upper(ale_pos pu
, int pup
) {
2761 perturb_upper_percent
= pup
;
2765 * Maximum rotational perturbation.
2768 static void set_rot_max(int rm
) {
2771 * Obtain the largest power of two not larger than rm.
2774 rot_max
= pow(2, floor(log(rm
) / log(2)));
2778 * Barrel distortion adjustment multiplier
2781 static void set_bda_mult(ale_pos m
) {
2786 * Barrel distortion maximum rate of change
2789 static void set_bda_rate(ale_pos m
) {
2797 static void set_lod_max(int lm
) {
2802 * Set the scale factor
2804 static void set_scale(ale_pos s
) {
2809 * Set reference rendering to align against
2811 static void set_reference(render
*r
) {
2816 * Set the interpolant
2818 static void set_interpolant(filter::scaled_filter
*f
) {
2823 * Set alignment weights image
2825 static void set_weight_map(const image
*i
) {
2830 * Set frequency cuts
2832 static void set_frequency_cut(double h
, double v
, double a
) {
2839 * Set algorithmic alignment weighting
2841 static void set_wmx(const char *e
, const char *f
, const char *d
) {
2848 * Show frequency weights
2850 static void set_fl_show(const char *filename
) {
2851 fw_output
= filename
;
2855 * Set transformation file loader.
2857 static void set_tload(tload_t
*tl
) {
2862 * Set transformation file saver.
2864 static void set_tsave(tsave_t
*ts
) {
2869 * Get match statistics for frame N.
2871 static int match(int n
) {
2885 * Message that old alignment data should be kept.
2887 static void keep() {
2888 assert (latest
== -1);
2893 * Get alignment for frame N.
2895 static transformation
of(int n
) {
2908 * Use Monte Carlo alignment sampling with argument N.
2910 static void mc(ale_pos n
) {
2914 static void mcd_limit(int n
) {
2919 * Set the certainty-weighted flag.
2921 static void certainty_weighted(int flag
) {
2922 certainty_weights
= flag
;
2926 * Set the global search type.
2928 static void gs(const char *type
) {
2929 if (!strcmp(type
, "local")) {
2931 } else if (!strcmp(type
, "inner")) {
2933 } else if (!strcmp(type
, "outer")) {
2935 } else if (!strcmp(type
, "all")) {
2937 } else if (!strcmp(type
, "central")) {
2939 } else if (!strcmp(type
, "points")) {
2943 ui::get()->error("bad global search type");
2948 * Set the minimum overlap for global searching
2950 static void gs_mo(unsigned int value
) {
2955 * Don't use Monte Carlo alignment sampling.
2957 static void no_mc() {
2962 * Set alignment exclusion regions
2964 static void set_exclusion(exclusion
*_ax_parameters
, int _ax_count
) {
2965 ax_count
= _ax_count
;
2966 ax_parameters
= _ax_parameters
;
2970 * Get match summary statistics.
2972 static ale_accum
match_summary() {
2973 return match_sum
/ match_count
;