1 // Copyright 2002, 2004, 2007 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 3 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 type.
164 * 0. Identity transformation.
166 * 1. Most recently accepted frame's final transformation.
169 static int default_initial_alignment_type
;
172 * Projective group behavior
174 * 0. Perturb in output coordinates.
176 * 1. Perturb in source coordinates
179 static int perturb_type
;
184 * This structure contains variables necessary for handling a
185 * multi-alignment element. The change between the non-default old
186 * initial alignment and old final alignment is used to adjust the
187 * non-default current initial alignment. If either the old or new
188 * initial alignment is a default alignment, the old --follow semantics
193 int is_default
, old_is_default
;
196 transformation old_initial_alignment
;
197 transformation old_final_alignment
;
198 transformation default_initial_alignment
;
199 const image
*input_frame
;
209 * Alignment for failed frames -- default or optimal?
211 * A frame that does not meet the match threshold can be assigned the
212 * best alignment found, or can be assigned its alignment default.
215 static int is_fail_default
;
220 * 0. Align images with an error contribution from each color channel.
222 * 1. Align images with an error contribution only from the green channel.
223 * Other color channels do not affect alignment.
225 * 2. Align images using a summation of channels. May be useful when dealing
226 * with images that have high frequency color ripples due to color aliasing.
229 static int channel_alignment_type
;
232 * Error metric exponent
235 static float metric_exponent
;
241 static float match_threshold
;
244 * Perturbation lower and upper bounds.
247 static ale_pos perturb_lower
;
248 static int perturb_lower_percent
;
249 static ale_pos perturb_upper
;
250 static int perturb_upper_percent
;
253 * Maximum level-of-detail scale factor is 2^lod_max/perturb.
259 * Maximum rotational perturbation
262 static ale_pos rot_max
;
265 * Barrel distortion alignment multiplier
268 static ale_pos bda_mult
;
271 * Barrel distortion maximum adjustment rate
274 static ale_pos bda_rate
;
277 * Alignment match sum
280 static ale_accum match_sum
;
283 * Alignment match count.
286 static int match_count
;
289 * Monte Carlo parameter
291 * 0. Don't use monte carlo alignment sampling.
293 * (0,1]. Select, on average, a number of pixels which is the
294 * specified fraction of the number of pixels in the accumulated image.
298 static unsigned int _mcd_limit
;
301 * Certainty weight flag
303 * 0. Don't use certainty weights for alignment.
305 * 1. Use certainty weights for alignment.
308 static int certainty_weights
;
311 * Global search parameter
313 * 0. Local: Local search only.
314 * 1. Inner: Alignment reference image inner region
315 * 2. Outer: Alignment reference image outer region
316 * 3. All: Alignment reference image inner and outer regions.
317 * 4. Central: Inner if possible; else, best of inner and outer.
318 * 5. Points: Align by control points.
324 * Multi-alignment cardinality.
327 static unsigned int _ma_card
;
330 * Multi-alignment contiguity.
333 static double _ma_cont
;
336 * Minimum overlap for global searches
339 static ale_pos _gs_mo
;
340 static int gs_mo_percent
;
346 static exclusion
*ax_parameters
;
350 * Type for scale cluster.
353 struct scale_cluster
{
355 const image
*defined
;
356 const image
*aweight
;
357 exclusion
*ax_parameters
;
364 * Check for exclusion region coverage in the reference
367 static int ref_excluded(int i
, int j
, point offset
, exclusion
*params
, int param_count
) {
368 for (int idx
= 0; idx
< param_count
; idx
++)
369 if (params
[idx
].type
== exclusion::RENDER
370 && i
+ offset
[0] >= params
[idx
].x
[0]
371 && i
+ offset
[0] <= params
[idx
].x
[1]
372 && j
+ offset
[1] >= params
[idx
].x
[2]
373 && j
+ offset
[1] <= params
[idx
].x
[3])
380 * Check for exclusion region coverage in the input
383 static int input_excluded(ale_pos ti
, ale_pos tj
, exclusion
*params
, int param_count
) {
384 for (int idx
= 0; idx
< param_count
; idx
++)
385 if (params
[idx
].type
== exclusion::FRAME
386 && ti
>= params
[idx
].x
[0]
387 && ti
<= params
[idx
].x
[1]
388 && tj
>= params
[idx
].x
[2]
389 && tj
<= params
[idx
].x
[3])
396 * Overlap function. Determines the number of pixels in areas where
397 * the arrays overlap. Uses the reference array's notion of pixel
400 static unsigned int overlap(struct scale_cluster c
, transformation t
, int ax_count
) {
401 assert (reference_image
);
403 unsigned int result
= 0;
405 point offset
= c
.accum
->offset();
407 for (unsigned int i
= 0; i
< c
.accum
->height(); i
++)
408 for (unsigned int j
= 0; j
< c
.accum
->width(); j
++) {
410 if (ref_excluded(i
, j
, offset
, c
.ax_parameters
, ax_count
))
419 q
= (c
.input_scale
< 1.0 && interpolant
== NULL
)
420 ? t
.scaled_inverse_transform(
421 point(i
+ offset
[0], j
+ offset
[1]))
422 : t
.unscaled_inverse_transform(
423 point(i
+ offset
[0], j
+ offset
[1]));
429 * Check that the transformed coordinates are within
430 * the boundaries of array c.input, and check that the
431 * weight value in the accumulated array is nonzero,
432 * unless we know it is nonzero by virtue of the fact
433 * that it falls within the region of the original
434 * frame (e.g. when we're not increasing image
435 * extents). Also check for frame exclusion.
438 if (input_excluded(ti
, tj
, c
.ax_parameters
, ax_count
))
442 && ti
<= c
.input
->height() - 1
444 && tj
<= c
.input
->width() - 1
445 && c
.defined
->get_pixel(i
, j
)[0] != 0)
453 * Calculate the region associated with the current multi-alignment
456 static void calculate_element_region(transformation
*t
, scale_cluster si
,
457 int local_ax_count
) {
459 unsigned int i_max
= si
.accum
->height();
460 unsigned int j_max
= si
.accum
->width();
461 point offset
= si
.accum
->offset();
463 if (si
.input_scale
< 1.0 && interpolant
== NULL
)
464 t
->begin_calculate_scaled_region(i_max
, j_max
, offset
);
466 t
->begin_calculate_unscaled_region(i_max
, j_max
, offset
);
468 for (unsigned int i
= 0; i
< i_max
; i
++)
469 for (unsigned int j
= 0; j
< j_max
; j
++) {
471 if (ref_excluded(i
, j
, offset
, si
.ax_parameters
, local_ax_count
))
476 while ((q
= t
->get_query_point((int) (i
+ offset
[0]),
477 (int) (j
+ offset
[1]))).defined()) {
482 if (input_excluded(ti
, tj
, si
.ax_parameters
, ax_count
))
486 && ti
<= si
.input
->height() - 1
488 && tj
<= si
.input
->width() - 1
489 && si
.defined
->get_pixel(i
, j
)[0] != 0) {
496 t
->end_calculate_region();
500 * Not-quite-symmetric difference function. Determines the difference in areas
501 * where the arrays overlap. Uses the reference array's notion of pixel positions.
503 * For the purposes of determining the difference, this function divides each
504 * pixel value by the corresponding image's average pixel magnitude, unless we
507 * a) Extending the boundaries of the image, or
509 * b) following the previous frame's transform
511 * If we are doing monte-carlo pixel sampling for alignment, we
512 * typically sample a subset of available pixels; otherwise, we sample
519 struct removal_pair
{
522 ale_pos zero_favorability_gradient
;
527 transformation offset
;
535 ale_accum centroid
[2], centroid_divisor
;
536 ale_accum de_centroid
[2], de_centroid_v
, de_sum
;
538 std::vector
<removal_pair
> errors_favoring
[2];
540 typedef unsigned int hist_bin
;
554 errors_favoring
[0].clear();
555 errors_favoring
[1].clear();
557 hist_min_r
= INT_MIN
;
559 hist_min_d
= INT_MIN
;
562 histogram
= (hist_bin
*) calloc(hist_dim
* hist_dim
, sizeof(hist_bin
));
565 min
= point::posinf();
566 max
= point::neginf();
570 centroid_divisor
= 0;
580 void init(transformation _offset
, ale_pos _perturb
, ale_pos _mc
) {
588 * Required for STL sanity.
594 run(transformation _offset
, ale_pos _perturb
, ale_pos _mc
) : offset() {
595 init(_offset
, _perturb
, _mc
);
598 void add_removal(const removal_pair
&rp
) {
599 for (int r
= 0; r
< 2; r
++) {
600 for (unsigned int i
= errors_favoring
[r
].size();; i
--) {
603 || (r
&& rp
.zero_favorability_gradient
604 >= errors_favoring
[r
][i
- 1].zero_favorability_gradient
)
605 || (!r
&& rp
.zero_favorability_gradient
606 <= errors_favoring
[r
][i
- 1].zero_favorability_gradient
)) {
611 std::vector
<removal_pair
>::iterator ii
= errors_favoring
[r
].begin();
614 errors_favoring
[r
].insert(ii
, rp
);
619 if (errors_favoring
[r
].size() > _mcd_limit
)
620 errors_favoring
[r
].pop_back();
622 assert(errors_favoring
[r
].size() <= _mcd_limit
);
626 void add_hist(int r
, int d
, int count
) {
630 int r_shift
= 0, d_shift
= 0;
632 if (r
- hist_min_r
>= hist_dim
) {
633 r_shift
= (r
- hist_min_r
) - hist_dim
+ 1;
634 hist_min_r
+= r_shift
;
637 if (d
- hist_min_d
>= hist_dim
) {
638 d_shift
= (d
- hist_min_d
) - hist_dim
+ 1;
639 hist_min_d
+= d_shift
;
642 assert (r_shift
>= 0);
643 assert (d_shift
>= 0);
645 if (r_shift
|| d_shift
) {
646 for (int rr
= 0; rr
< hist_dim
; rr
++)
647 for (int dd
= 0; dd
< hist_dim
; dd
++) {
649 hist_bin value
= histogram
[rr
* hist_dim
+ dd
];
651 histogram
[rr
* hist_dim
+ dd
] = 0;
653 int rrr
= rr
- r_shift
;
654 int ddd
= dd
- d_shift
;
661 histogram
[rrr
* hist_dim
+ ddd
] += value
;
673 histogram
[r
* hist_dim
+ d
] += count
;
676 void add_hist(ale_accum result
, ale_accum divisor
) {
677 ale_accum rbin
= log(result
) / log(2);
678 ale_accum dbin
= log(divisor
) / log(2);
680 if (!(rbin
> INT_MIN
))
682 if (!(dbin
> INT_MIN
))
685 add_hist((int) floor(rbin
), (int) floor(dbin
), 1);
688 void add(const run
&_run
) {
689 result
+= _run
.result
;
690 divisor
+= _run
.divisor
;
692 for (int r
= 0; r
< 2; r
++)
693 for (unsigned int i
= 0; i
< _run
.errors_favoring
[r
].size(); i
++)
694 add_removal(_run
.errors_favoring
[r
][i
]);
696 for (int r
= 0; r
< _run
.hist_dim
; r
++)
697 for (int d
= 0; d
< _run
.hist_dim
; d
++)
698 add_hist(r
+ _run
.hist_min_r
, d
+ _run
.hist_min_d
,
699 _run
.histogram
[r
* hist_dim
+ d
]);
701 for (int d
= 0; d
< 2; d
++) {
702 if (min
[d
] > _run
.min
[d
])
703 min
[d
] = _run
.min
[d
];
704 if (max
[d
] < _run
.max
[d
])
705 max
[d
] = _run
.max
[d
];
706 centroid
[d
] += _run
.centroid
[d
];
707 de_centroid
[d
] += _run
.de_centroid
[d
];
710 centroid_divisor
+= _run
.centroid_divisor
;
711 de_centroid_v
+= _run
.de_centroid_v
;
712 de_sum
+= _run
.de_sum
;
715 run(const run
&_run
) : offset() {
720 init(_run
.offset
, _run
.perturb
, _run
.mc
);
728 run
&operator=(const run
&_run
) {
731 * Free old structures
738 init(_run
.offset
, _run
.perturb
, _run
.mc
);
752 ale_accum
get_error() const {
753 return pow(result
/ divisor
, 1/metric_exponent
);
756 int check_synchronized_removal(const run
*with
) const {
758 if (hist_total
< 100 || with
->hist_total
< 100)
763 if (get_error() > with
->get_error())
766 if ((sense
== 0 && !(get_error() < with
->get_error()))
767 || (sense
== 1 && !(get_error() > with
->get_error())))
770 ale_accum error
[2] = { with
->result
, this->result
};
771 ale_accum divisor
[2] = { with
->divisor
, this->divisor
};
773 for (unsigned int i
= 0; i
< errors_favoring
[sense
].size(); i
++)
774 for (unsigned int d
= 0; d
< 2; d
++) {
775 error
[d
] -= errors_favoring
[sense
][i
].error
[d
];
776 divisor
[d
] -= errors_favoring
[sense
][i
].divisor
[d
];
779 for (int d
= 0; d
< 2; d
++)
781 || !(divisor
[d
] >= 0))
784 if ((sense
== 0 && !(error
[0] / divisor
[0] > error
[1] / divisor
[1]))
785 || (sense
== 1 && !(error
[0] / divisor
[0] < error
[1] / divisor
[1])))
792 * XXX: A better removal technique would probably match
793 * removals between compared transformations (see
794 * check_synchronized_removal). This simpler approach
795 * is likely to under-estimate reliability for any
796 * given removal count, due to its less constrained
797 * method of performing removals.
799 int check_removal(const run
*with
) const {
801 if (hist_total
< 100 || with
->hist_total
< 100)
804 if (get_error() > with
->get_error())
805 return with
->check_removal(this);
807 if (!(get_error() < with
->get_error()))
810 ale_accum bresult
, bdivisor
, wresult
, wdivisor
;
811 hist_bin
*bhist
, *whist
;
813 ui::get()->d2_align_sim_start();
815 bhist
= (hist_bin
*)malloc(sizeof(hist_bin
) * hist_dim
* hist_dim
);
816 whist
= (hist_bin
*)malloc(sizeof(hist_bin
) * with
->hist_dim
* with
->hist_dim
);
820 wresult
= with
->result
;
821 wdivisor
= with
->divisor
;
823 for (int i
= 0; i
< hist_dim
* hist_dim
; i
++)
824 bhist
[i
] = histogram
[i
];
826 for (int i
= 0; i
< with
->hist_dim
* with
->hist_dim
; i
++)
827 whist
[i
] = with
->histogram
[i
];
829 for (unsigned int r
= 0; r
< _mcd_limit
; r
++) {
830 int max_gradient_bin
= -1;
831 int max_gradient_hist
= -1;
832 ale_accum max_gradient
= 0;
834 for (int i
= 0; i
< hist_dim
* hist_dim
; i
++) {
838 ale_accum br
= pow(2, hist_min_r
+ i
/ hist_dim
);
839 ale_accum bd
= pow(2, hist_min_d
+ i
% hist_dim
);
840 ale_accum b_test_gradient
=
841 (bresult
- br
) / (bdivisor
- bd
) - bresult
/ bdivisor
;
843 if (!(b_test_gradient
< max_gradient
)) {
844 max_gradient_bin
= i
;
845 max_gradient_hist
= 0;
846 max_gradient
= b_test_gradient
;
850 for (int i
= 0; i
< with
->hist_dim
* with
->hist_dim
; i
++) {
854 ale_accum wr
= pow(2, with
->hist_min_r
+ i
/ with
->hist_dim
);
855 ale_accum wd
= pow(2, with
->hist_min_d
+ i
% with
->hist_dim
);
856 ale_accum w_test_gradient
=
857 wresult
/ wdivisor
- (wresult
- wr
) / (wdivisor
- wd
);
859 if (!(w_test_gradient
< max_gradient
)) {
860 max_gradient_bin
= i
;
861 max_gradient_hist
= 1;
862 max_gradient
= w_test_gradient
;
866 if (max_gradient_hist
== -1)
869 if (max_gradient_hist
== 0) {
870 bhist
[max_gradient_bin
]--;
871 bresult
-= pow(2, hist_min_r
+ max_gradient_bin
/ hist_dim
);
872 bdivisor
-= pow(2, hist_min_d
+ max_gradient_bin
% hist_dim
);
873 } else if (max_gradient_hist
== 1) {
874 whist
[max_gradient_bin
]--;
875 wresult
-= pow(2, with
->hist_min_r
+ max_gradient_bin
/ with
->hist_dim
);
876 wdivisor
-= pow(2, with
->hist_min_d
+ max_gradient_bin
% with
->hist_dim
);
883 ui::get()->d2_align_sim_stop();
885 if (finite(bresult
/ bdivisor
)
886 && finite(wresult
/ wdivisor
)
887 && bresult
/ bdivisor
< wresult
/ wdivisor
)
893 void sample(int f
, scale_cluster c
, int i
, int j
, point t
, point u
,
894 const run
&comparison
) {
896 pixel pa
= c
.accum
->get_pixel(i
, j
);
900 ale_accum this_result
[2] = { 0, 0 };
901 ale_accum this_divisor
[2] = { 0, 0 };
903 if (interpolant
!= NULL
)
904 interpolant
->filtered(i
, j
, &p
[0], &weight
[0], 1, f
);
908 c
.input
->get_bl(t
, result
);
910 weight
[0] = result
[1];
915 c
.input
->get_bl(u
, result
);
917 weight
[1] = result
[1];
925 if (certainty_weights
== 0) {
926 weight
[0] = pixel(1, 1, 1);
927 weight
[1] = pixel(1, 1, 1);
930 if (c
.aweight
!= NULL
) {
931 weight
[0] *= c
.aweight
->get_pixel(i
, j
);
932 weight
[1] *= c
.aweight
->get_pixel(i
, j
);
936 * Update sampling area statistics
948 centroid
[0] += (weight
[0][0] + weight
[0][1] + weight
[0][2]) * i
;
949 centroid
[1] += (weight
[0][0] + weight
[0][1] + weight
[0][2]) * j
;
950 centroid_divisor
+= (weight
[0][0] + weight
[0][1] + weight
[0][2]);
953 * Determine alignment type.
956 for (int m
= 0; m
< (u
.defined() ? 2 : 1); m
++)
957 if (channel_alignment_type
== 0) {
959 * Align based on all channels.
963 for (int k
= 0; k
< 3; k
++) {
964 ale_real achan
= pa
[k
];
965 ale_real bchan
= p
[m
][k
];
967 this_result
[m
] += weight
[m
][k
] * pow(fabs(achan
- bchan
), metric_exponent
);
968 this_divisor
[m
] += weight
[m
][k
] * pow(achan
> bchan
? achan
: bchan
, metric_exponent
);
970 } else if (channel_alignment_type
== 1) {
972 * Align based on the green channel.
975 ale_real achan
= pa
[1];
976 ale_real bchan
= p
[m
][1];
978 this_result
[m
] = weight
[m
][1] * pow(fabs(achan
- bchan
), metric_exponent
);
979 this_divisor
[m
] = weight
[m
][1] * pow(achan
> bchan
? achan
: bchan
, metric_exponent
);
980 } else if (channel_alignment_type
== 2) {
982 * Align based on the sum of all channels.
989 for (int k
= 0; k
< 3; k
++) {
992 wsum
+= weight
[m
][k
] / 3;
995 this_result
[m
] = wsum
* pow(fabs(asum
- bsum
), metric_exponent
);
996 this_divisor
[m
] = wsum
* pow(asum
> bsum
? asum
: bsum
, metric_exponent
);
1000 // ale_accum de = fabs(this_result[0] / this_divisor[0]
1001 // - this_result[1] / this_divisor[1]);
1002 ale_accum de
= fabs(this_result
[0] - this_result
[1]);
1004 de_centroid
[0] += de
* i
;
1005 de_centroid
[1] += de
* j
;
1007 de_centroid_v
+= de
* t
.lengthto(u
);
1013 rp
.error
[0] = this_result
[0];
1014 rp
.error
[1] = this_result
[1];
1015 rp
.divisor
[0] = this_divisor
[0];
1016 rp
.divisor
[1] = this_divisor
[1];
1018 ale_accum result_approx
= comparison
.result
;
1019 ale_accum divisor_approx
= comparison
.divisor
;
1021 rp
.zero_favorability_gradient
=
1022 (result_approx
- rp
.error
[1])
1023 / (divisor_approx
- rp
.error
[1])
1024 - (result_approx
- rp
.error
[0])
1025 / (divisor_approx
- rp
.error
[0]);
1030 result
+= (this_result
[0]);
1031 divisor
+= (this_divisor
[0]);
1033 add_hist(this_result
[0], this_divisor
[0]);
1037 fprintf(stderr
, "\n");
1038 fprintf(stderr
, "hist_min_r = %d\n", hist_min_r
);
1039 fprintf(stderr
, "hist_min_d = %d\n", hist_min_d
);
1040 fprintf(stderr
, "hist_dim = %d\n", hist_dim
);
1041 fprintf(stderr
, "hist_total = %d\n", hist_total
);
1042 fprintf(stderr
, "\n");
1044 hist_bin recalc_total
= 0;
1046 for (int r
= 0; r
< hist_dim
; r
++) {
1047 for (int d
= 0; d
< hist_dim
; d
++) {
1048 recalc_total
+= histogram
[r
* hist_dim
+ d
];
1049 fprintf(stderr
, "\t%d", histogram
[r
* hist_dim
+ d
]);
1051 fprintf(stderr
, "\n");
1054 fprintf(stderr
, "\n");
1055 fprintf(stderr
, "recalc_total = %d\n", recalc_total
);
1058 void rescale(ale_pos scale
) {
1059 offset
.rescale(scale
);
1061 de_centroid
[0] *= scale
;
1062 de_centroid
[1] *= scale
;
1063 de_centroid_v
*= scale
;
1066 point
get_centroid() {
1067 point result
= point(centroid
[0] / centroid_divisor
, centroid
[1] / centroid_divisor
);
1069 assert (finite(centroid
[0])
1070 && finite(centroid
[1])
1071 && (result
.defined() || centroid_divisor
== 0));
1076 point
get_error_centroid() {
1077 point result
= point(de_centroid
[0] / de_sum
, de_centroid
[1] / de_sum
);
1082 ale_pos
get_error_perturb() {
1083 ale_pos result
= de_centroid_v
/ de_sum
;
1091 * When non-empty, runs.front() is best, runs.back() is
1095 std::vector
<run
> runs
;
1098 * old_runs stores the latest available perturbation set for
1099 * each multi-alignment element.
1102 typedef std::pair
<unsigned int, unsigned int> run_index
;
1103 std::map
<run_index
, run
> old_runs
;
1104 std::map
<run_index
, ale_pos
> current_mc
;
1106 static void *diff_subdomain(void *args
);
1108 struct subdomain_args
{
1109 struct scale_cluster c
;
1110 std::vector
<run
> runs
;
1114 int i_min
, i_max
, j_min
, j_max
;
1122 transformation offset
[2];
1125 memo_entry(int _frame
, ale_pos _perturb
, transformation dest
, transformation origin
, ale_pos _mc
) {
1133 int subequals(const memo_entry
&m
, int o
, int u
) const {
1135 if (frame
!= m
.frame
)
1138 if (perturb
!= m
.perturb
)
1141 if (offset
[o
] != m
.offset
[u
])
1147 int equals(const memo_entry
&m
) const {
1148 return subequals(m
, 0, 0) && subequals(m
, 1, 1);
1151 int conflicts(const memo_entry
&m
) const {
1152 memo_entry
mr(m
.frame
, m
.perturb
, m
.offset
[1], m
.offset
[0], mc
);
1158 static std::vector
<memo_entry
> memos
;
1160 static rng_t loop_handling_randomness
;
1162 static int is_loop(std::vector
<memo_entry
> mf
) {
1164 assert (mf
.size() > 0);
1167 * Loop containing me
1169 if (mf
.back().subequals(mf
.front(), 1, 0))
1173 * Loop not containing me
1175 for (unsigned int i
= 1; i
< mf
.size(); i
++) {
1176 if (mf
.back().subequals(mf
[i
], 1, 0))
1180 for (unsigned int i
= 0; i
< memos
.size(); i
++) {
1181 if (memos
[i
].subequals(mf
.back(), 0, 1)) {
1182 mf
.push_back(memos
[i
]);
1191 if (loop_handling_randomness
.get() % 10 == 0)
1198 int get_current_index() const {
1199 assert(runs
.size());
1200 return runs
[0].offset
.get_current_index();
1203 struct scale_cluster si
;
1207 ale_accum mc_default
;
1209 std::vector
<ale_pos
> perturb_multipliers
;
1212 ale_pos
get_mc() const {
1213 assert(runs
.size() >= 1);
1217 static void clear_memos() {
1221 void diff(struct scale_cluster c
, ale_pos _mc_arg
, ale_pos perturb
,
1223 int _ax_count
, int f
) {
1225 if (runs
.size() == 2)
1228 if (runs
.size() > 0 && runs
[0].mc
!= _mc_arg
) {
1229 runs
[0].mc
= _mc_arg
;
1233 runs
.push_back(run(t
, perturb
, _mc_arg
));
1236 ax_count
= _ax_count
;
1239 assert (_mc_arg
> 0);
1241 ui::get()->d2_align_sample_start();
1243 if (interpolant
!= NULL
)
1244 interpolant
->set_parameters(t
, c
.input
, c
.accum
->offset());
1248 N
= thread::count();
1250 pthread_t
*threads
= (pthread_t
*) malloc(sizeof(pthread_t
) * N
);
1251 pthread_attr_t
*thread_attr
= (pthread_attr_t
*) malloc(sizeof(pthread_attr_t
) * N
);
1257 subdomain_args
*args
= new subdomain_args
[N
];
1259 for (int ti
= 0; ti
< N
; ti
++) {
1261 args
[ti
].runs
= runs
;
1262 args
[ti
]._mc_arg
= _mc_arg
;
1263 args
[ti
].ax_count
= ax_count
;
1265 args
[ti
].i_min
= (c
.accum
->height() * ti
) / N
;
1266 args
[ti
].i_max
= (c
.accum
->height() * (ti
+ 1)) / N
;
1268 args
[ti
].j_max
= c
.accum
->width();
1269 args
[ti
].subdomain
= ti
;
1271 pthread_attr_init(&thread_attr
[ti
]);
1272 pthread_attr_setdetachstate(&thread_attr
[ti
], PTHREAD_CREATE_JOINABLE
);
1273 pthread_create(&threads
[ti
], &thread_attr
[ti
], diff_subdomain
, &args
[ti
]);
1275 diff_subdomain(&args
[ti
]);
1279 for (int ti
= 0; ti
< N
; ti
++) {
1281 pthread_join(threads
[ti
], NULL
);
1283 runs
.back().add(args
[ti
].runs
.back());
1288 ui::get()->d2_align_sample_stop();
1294 std::vector
<transformation
> t_array
;
1295 std::vector
<ale_pos
> p_array
;
1296 std::vector
<ale_pos
> mc_array
;
1298 for (unsigned int r
= 0; r
< runs
.size(); r
++) {
1299 t_array
.push_back(runs
[r
].offset
);
1300 p_array
.push_back(runs
[r
].perturb
);
1301 mc_array
.push_back(runs
[r
].mc
);
1306 for (unsigned int r
= 0; r
< t_array
.size(); r
++)
1307 diff(si
, mc_array
[r
], p_array
[r
], t_array
[r
], ax_count
, frame
);
1310 int reliable() const {
1311 assert(runs
.size() >= 2);
1319 if (!runs
[1].check_synchronized_removal(&runs
[0]))
1322 memo_entry
me(frame
, runs
[1].perturb
, runs
[1].offset
,
1323 runs
[0].offset
, get_mc());
1327 if (runs
[1].get_error() > runs
[0].get_error()) {
1329 me
.offset
[0] = runs
[0].offset
;
1330 me
.offset
[1] = runs
[1].offset
;
1333 for (unsigned int i
= 0; i
< memos
.size(); i
++) {
1335 if (better
&& me
.subequals(memos
[i
], 1, 1)) {
1336 if (me
.mc
> memos
[i
].mc
) {
1343 if (0 && me
.conflicts(memos
[i
])) {
1344 if (me
.mc
> memos
[i
].mc
) {
1355 memos
.push_back(me
);
1363 assert(runs
.size() >= 2);
1364 assert(runs
[0].offset
.scale() == runs
[1].offset
.scale());
1368 while (!reliable()) {
1378 return (runs
[1].get_error() < runs
[0].get_error()
1379 || (!finite(runs
[0].get_error()) && finite(runs
[1].get_error())));
1382 diff_stat_t(ale_pos _mc_arg
) : runs(), old_runs(), current_mc(),
1383 perturb_multipliers() {
1385 mc_default
= _mc_arg
;
1387 // ui::get()->alignment_monte_carlo_parameter(_mc_arg);
1390 run_index
get_run_index(unsigned int perturb_index
) {
1391 return run_index(get_current_index(), perturb_index
);
1394 run
&get_run(unsigned int perturb_index
) {
1395 run_index index
= get_run_index(perturb_index
);
1397 assert(old_runs
.count(index
));
1398 return old_runs
[index
];
1401 void rescale(ale_pos scale
, scale_cluster _si
) {
1402 assert(runs
.size() == 1);
1406 runs
[0].rescale(scale
);
1409 for (unsigned int i
= 0;
1410 current_mc
.count(run_index(i
, 0)); i
++)
1411 for (unsigned int j
= 0;
1412 current_mc
.count(run_index(i
, j
)); j
++)
1413 current_mc
[run_index(i
, j
)] /= pow(scale
, 2);
1415 runs
[0].mc
/= pow(scale
, 2);
1417 mc_default
/= pow(scale
, 2);
1421 ui::get()->alignment_monte_carlo_parameter(get_mc());
1427 assert(runs
.size() == 1);
1432 for (unsigned int i
= 0;
1433 current_mc
.count(run_index(i
, 0)); i
++)
1434 for (unsigned int j
= 0;
1435 current_mc
.count(run_index(i
, j
)); j
++)
1436 current_mc
[run_index(i
, j
)] /= 2;
1440 ui::get()->alignment_monte_carlo_parameter(get_mc());
1445 void push_element() {
1446 assert(runs
.size() == 1);
1448 runs
[0].offset
.push_element();
1453 unsigned int get_current_index() {
1454 assert (runs
.size() > 0);
1456 return runs
[0].offset
.get_current_index();
1459 void set_current_index(unsigned int i
) {
1460 assert(runs
.size() == 1);
1461 runs
[0].offset
.set_current_index(i
);
1465 void calculate_element_region() {
1466 assert(runs
.size() == 1);
1468 if (get_offset().get_current_index() > 0
1469 && get_offset().is_nontrivial())
1470 align::calculate_element_region(&runs
[0].offset
, si
, ax_count
);
1476 diff_stat_t
&operator=(const diff_stat_t
&dst
) {
1478 * Copy run information.
1481 old_runs
= dst
.old_runs
;
1482 current_mc
= dst
.current_mc
;
1485 * Copy diff variables
1487 mc_default
= dst
.mc_default
;
1489 ax_count
= dst
.ax_count
;
1491 perturb_multipliers
= dst
.perturb_multipliers
;
1496 diff_stat_t(const diff_stat_t
&dst
) : runs(), old_runs(),
1497 perturb_multipliers() {
1501 ale_accum
get_result() {
1502 assert(runs
.size() == 1);
1503 return runs
[0].result
;
1506 ale_accum
get_divisor() {
1507 assert(runs
.size() == 1);
1508 return runs
[0].divisor
;
1511 transformation
get_offset() {
1512 assert(runs
.size() == 1);
1513 return runs
[0].offset
;
1516 int operator!=(diff_stat_t
¶m
) {
1517 return (get_mc() != param
.get_mc()
1518 || get_error() != param
.get_error());
1521 int operator==(diff_stat_t
¶m
) {
1522 return !(operator!=(param
));
1525 ale_pos
get_error_perturb() {
1526 assert(runs
.size() == 1);
1527 return runs
[0].get_error_perturb();
1530 ale_accum
get_error() const {
1531 assert(runs
.size() == 1);
1532 return runs
[0].get_error();
1537 void mcd_increase() {
1539 assert (runs
.size() > 0);
1541 for (unsigned int r
= 0; r
< runs
.size(); r
++)
1544 // ui::get()->alignment_monte_carlo_parameter(get_mc());
1549 void mcd_decrease() {
1551 assert (runs
.size() > 0);
1553 for (unsigned int r
= 0; r
< runs
.size(); r
++)
1556 // ui::get()->alignment_monte_carlo_parameter(get_mc());
1561 void try_decrease();
1565 * Get the set of transformations produced by a given perturbation
1567 void get_perturb_set(std::vector
<transformation
> *set
,
1568 ale_pos adj_p
, ale_pos adj_o
, ale_pos adj_b
,
1569 ale_pos
*current_bd
, ale_pos
*modified_bd
,
1570 std::vector
<ale_pos
> multipliers
= std::vector
<ale_pos
>()) {
1572 assert(runs
.size() == 1);
1574 transformation test_t
;
1577 * Translational or euclidean transformation
1580 for (unsigned int i
= 0; i
< 2; i
++)
1581 for (unsigned int s
= 0; s
< 2; s
++) {
1583 if (!multipliers
.size())
1584 multipliers
.push_back(1);
1586 assert(finite(multipliers
[0]));
1588 test_t
= get_offset();
1590 // test_t.eu_modify(i, (s ? -adj_p : adj_p) * multipliers[0]);
1591 test_t
.translate((i
? point(1, 0) : point(0, 1))
1592 * (s
? -adj_p
: adj_p
)
1595 test_t
.snap(adj_p
/ 2);
1597 set
->push_back(test_t
);
1598 multipliers
.erase(multipliers
.begin());
1602 if (alignment_class
> 0)
1603 for (unsigned int s
= 0; s
< 2; s
++) {
1605 if (!multipliers
.size())
1606 multipliers
.push_back(1);
1608 assert(multipliers
.size());
1609 assert(finite(multipliers
[0]));
1611 if (!(adj_o
* multipliers
[0] < rot_max
))
1614 ale_pos adj_s
= (s
? 1 : -1) * adj_o
* multipliers
[0];
1616 test_t
= get_offset();
1618 run_index ori
= get_run_index(set
->size());
1619 point centroid
= point::undefined();
1621 if (!old_runs
.count(ori
))
1622 ori
= get_run_index(0);
1624 if (!centroid
.finite() && old_runs
.count(ori
)) {
1625 centroid
= old_runs
[ori
].get_error_centroid();
1627 if (!centroid
.finite())
1628 centroid
= old_runs
[ori
].get_centroid();
1630 centroid
*= test_t
.scale()
1631 / old_runs
[ori
].offset
.scale();
1634 if (!centroid
.finite() && !test_t
.is_projective()) {
1635 test_t
.eu_modify(2, adj_s
);
1636 } else if (!centroid
.finite()) {
1637 centroid
= point(si
.input
->height() / 2,
1638 si
.input
->width() / 2);
1640 test_t
.rotate(centroid
+ si
.accum
->offset(),
1643 test_t
.rotate(centroid
+ si
.accum
->offset(),
1647 test_t
.snap(adj_p
/ 2);
1649 set
->push_back(test_t
);
1650 multipliers
.erase(multipliers
.begin());
1653 if (alignment_class
== 2) {
1656 * Projective transformation
1659 for (unsigned int i
= 0; i
< 4; i
++)
1660 for (unsigned int j
= 0; j
< 2; j
++)
1661 for (unsigned int s
= 0; s
< 2; s
++) {
1663 if (!multipliers
.size())
1664 multipliers
.push_back(1);
1666 assert(multipliers
.size());
1667 assert(finite(multipliers
[0]));
1669 ale_pos adj_s
= (s
? -1 : 1) * adj_p
* multipliers
[0];
1671 test_t
= get_offset();
1673 if (perturb_type
== 0)
1674 test_t
.gpt_modify(j
, i
, adj_s
);
1675 else if (perturb_type
== 1)
1676 test_t
.gr_modify(j
, i
, adj_s
);
1680 test_t
.snap(adj_p
/ 2);
1682 set
->push_back(test_t
);
1683 multipliers
.erase(multipliers
.begin());
1692 if (bda_mult
!= 0 && adj_b
!= 0) {
1694 for (unsigned int d
= 0; d
< get_offset().bd_count(); d
++)
1695 for (unsigned int s
= 0; s
< 2; s
++) {
1697 if (!multipliers
.size())
1698 multipliers
.push_back(1);
1700 assert (multipliers
.size());
1701 assert (finite(multipliers
[0]));
1703 ale_pos adj_s
= (s
? -1 : 1) * adj_b
* multipliers
[0];
1705 if (bda_rate
> 0 && fabs(modified_bd
[d
] + adj_s
- current_bd
[d
]) > bda_rate
)
1708 transformation test_t
= get_offset();
1710 test_t
.bd_modify(d
, adj_s
);
1712 set
->push_back(test_t
);
1718 assert(runs
.size() == 2);
1724 assert(runs
.size() == 2);
1728 void perturb_test(ale_pos perturb
, ale_pos adj_p
, ale_pos adj_o
, ale_pos adj_b
,
1729 ale_pos
*current_bd
, ale_pos
*modified_bd
, int stable
) {
1731 assert(runs
.size() == 1);
1733 std::vector
<transformation
> t_set
;
1735 if (perturb_multipliers
.size() == 0) {
1736 get_perturb_set(&t_set
, adj_p
, adj_o
, adj_b
,
1737 current_bd
, modified_bd
);
1739 for (unsigned int i
= 0; i
< t_set
.size(); i
++) {
1740 diff_stat_t test
= *this;
1742 test
.diff(si
, mc_default
, perturb
, t_set
[i
], ax_count
, frame
);
1746 if (finite(adj_p
/ test
.get_error_perturb()))
1747 perturb_multipliers
.push_back(adj_p
/ test
.get_error_perturb());
1749 perturb_multipliers
.push_back(1);
1756 get_perturb_set(&t_set
, adj_p
, adj_o
, adj_b
, current_bd
, modified_bd
,
1757 perturb_multipliers
);
1759 int found_unreliable
= 1;
1760 std::vector
<int> tested(t_set
.size(), 0);
1762 ale_pos current_mc_arg
= mc_default
;
1764 for (unsigned int i
= 0; i
< t_set
.size(); i
++) {
1765 run_index ori
= get_run_index(i
);
1768 * Check for stability
1771 && old_runs
.count(ori
)
1772 && old_runs
[ori
].offset
== t_set
[i
])
1776 * Update starting mc arg.
1778 if (current_mc
.count(ori
) && current_mc
[ori
] < current_mc_arg
)
1779 current_mc_arg
= current_mc
[ori
];
1782 std::vector
<ale_pos
> perturb_multipliers_original
= perturb_multipliers
;
1784 while (found_unreliable
) {
1786 ui::get()->alignment_monte_carlo_parameter(current_mc_arg
);
1788 found_unreliable
= 0;
1790 for (unsigned int i
= 0; i
< t_set
.size(); i
++) {
1795 diff(si
, current_mc_arg
, perturb
, t_set
[i
], ax_count
, frame
);
1797 if (!(i
< perturb_multipliers
.size())
1798 || !finite(perturb_multipliers
[i
])) {
1800 perturb_multipliers
.resize(i
+ 1);
1802 perturb_multipliers
[i
] =
1803 adj_p
/ runs
[1].get_error_perturb();
1805 if (finite(perturb_multipliers
[i
]))
1806 found_unreliable
= 1;
1811 run_index ori
= get_run_index(i
);
1813 if (old_runs
.count(ori
) == 0)
1814 old_runs
.insert(std::pair
<run_index
, run
>(ori
, runs
[1]));
1816 old_runs
[ori
] = runs
[1];
1818 perturb_multipliers
[i
] = perturb_multipliers_original
[i
]
1819 * adj_p
/ runs
[1].get_error_perturb();
1821 if (!finite(perturb_multipliers
[i
]))
1822 perturb_multipliers
[i
] = 1;
1824 if ((!current_mc
.count(ori
)
1825 && current_mc_arg
< 1)
1826 || (current_mc
.count(ori
)
1827 && current_mc
[ori
] > current_mc_arg
1828 && current_mc_arg
< 1)
1830 found_unreliable
= 1;
1834 current_mc
[ori
] = runs
[1].mc
;
1838 && runs
[1].get_error() < runs
[0].get_error()
1839 && perturb_multipliers
[i
]
1840 / perturb_multipliers_original
[i
] < 2) {
1841 current_mc
[ori
] = runs
[1].mc
;
1842 memos
.push_back(memo_entry(frame
, runs
[1].perturb
,
1843 runs
[1].offset
, runs
[0].offset
, get_mc()));
1851 if (runs
.size() > 1)
1854 if (!found_unreliable
)
1857 current_mc_arg
*= 2;
1862 * Attempt to make the current element non-trivial, by finding a nearby
1863 * alignment admitting a non-empty element region.
1865 void make_element_nontrivial(ale_pos adj_p
, ale_pos adj_o
) {
1866 assert(runs
.size() == 1);
1868 transformation
*t
= &runs
[0].offset
;
1870 if (t
->is_nontrivial())
1873 calculate_element_region();
1875 if (t
->is_nontrivial())
1878 std::vector
<transformation
> t_set
;
1879 get_perturb_set(&t_set
, adj_p
, adj_o
, 0, NULL
, NULL
);
1881 for (unsigned int i
= 0; i
< t_set
.size(); i
++) {
1883 align::calculate_element_region(&t_set
[i
], si
, ax_count
);
1885 if (t_set
[i
].is_nontrivial()) {
1896 * Adjust exposure for an aligned frame B against reference A.
1898 * Expects full-LOD images.
1900 * This function is a bit of a mess, as it reflects rather ad-hoc rules
1901 * regarding what seems to work w.r.t. certainty. Using certainty in the
1902 * first pass seems to result in worse alignment, while not using certainty
1903 * in the second pass results in incorrect determination of exposure.
1905 * [Note that this may have been due to a bug in certainty determination
1906 * within this function.]
1908 static void set_exposure_ratio(unsigned int m
, struct scale_cluster c
,
1909 transformation t
, int ax_count
, int pass_number
) {
1911 if (_exp_register
== 2) {
1913 * Use metadata only.
1915 ale_real gain_multiplier
= image_rw::exp(m
).get_gain_multiplier();
1916 pixel multiplier
= pixel(gain_multiplier
, gain_multiplier
, gain_multiplier
);
1918 image_rw::exp(m
).set_multiplier(multiplier
);
1919 ui::get()->exp_multiplier(multiplier
[0],
1926 pixel_accum asum
, bsum
;
1928 point offset
= c
.accum
->offset();
1930 for (unsigned int i
= 0; i
< c
.accum
->height(); i
++)
1931 for (unsigned int j
= 0; j
< c
.accum
->width(); j
++) {
1933 if (ref_excluded(i
, j
, offset
, c
.ax_parameters
, ax_count
))
1942 q
= (c
.input_scale
< 1.0 && interpolant
== NULL
)
1943 ? t
.scaled_inverse_transform(
1944 point(i
+ offset
[0], j
+ offset
[1]))
1945 : t
.unscaled_inverse_transform(
1946 point(i
+ offset
[0], j
+ offset
[1]));
1949 * Check that the transformed coordinates are within
1950 * the boundaries of array c.input, that they are not
1951 * subject to exclusion, and that the weight value in
1952 * the accumulated array is nonzero.
1955 if (input_excluded(q
[0], q
[1], c
.ax_parameters
, ax_count
))
1959 && q
[0] <= c
.input
->height() - 1
1961 && q
[1] <= c
.input
->width() - 1
1962 && c
.defined
->get_pixel(i
, j
).minabs_norm() != 0) {
1963 pixel a
= c
.accum
->get_pixel(i
, j
);
1966 c
.input
->get_bl(q
, b
);
1969 pixel weight
= (c
.aweight
1970 ? c
.aweight
->get_pixel(i
, j
)
1972 * ((!certainty_weights
&& pass_number
)
1973 ? c
.defined
->get_pixel(i
, j
)
1979 pixel weight
= pixel(1, 1, 1);
1983 bsum
+= b
[0] * weight
;
1987 // std::cerr << (asum / bsum) << " ";
1989 pixel_accum new_multiplier
;
1991 new_multiplier
= asum
/ bsum
* image_rw::exp(m
).get_multiplier();
1993 if (finite(new_multiplier
[0])
1994 && finite(new_multiplier
[1])
1995 && finite(new_multiplier
[2])) {
1996 image_rw::exp(m
).set_multiplier(new_multiplier
);
1997 ui::get()->exp_multiplier(new_multiplier
[0],
2004 * Copy all ax parameters.
2006 static exclusion
*copy_ax_parameters(int local_ax_count
, exclusion
*source
) {
2008 exclusion
*dest
= (exclusion
*) malloc(local_ax_count
* sizeof(exclusion
));
2013 ui::get()->memory_error("exclusion regions");
2015 for (int idx
= 0; idx
< local_ax_count
; idx
++)
2016 dest
[idx
] = source
[idx
];
2022 * Copy ax parameters according to frame.
2024 static exclusion
*filter_ax_parameters(int frame
, int *local_ax_count
) {
2026 exclusion
*dest
= (exclusion
*) malloc(ax_count
* sizeof(exclusion
));
2031 ui::get()->memory_error("exclusion regions");
2033 *local_ax_count
= 0;
2035 for (int idx
= 0; idx
< ax_count
; idx
++) {
2036 if (ax_parameters
[idx
].x
[4] > frame
2037 || ax_parameters
[idx
].x
[5] < frame
)
2040 dest
[*local_ax_count
] = ax_parameters
[idx
];
2042 (*local_ax_count
)++;
2048 static void scale_ax_parameters(int local_ax_count
, exclusion
*ax_parameters
,
2049 ale_pos ref_scale
, ale_pos input_scale
) {
2050 for (int i
= 0; i
< local_ax_count
; i
++) {
2051 ale_pos scale
= (ax_parameters
[i
].type
== exclusion::RENDER
)
2055 for (int n
= 0; n
< 6; n
++) {
2056 ax_parameters
[i
].x
[n
] = (int) floor(ax_parameters
[i
].x
[n
]
2063 * Prepare the next level of detail for ordinary images.
2065 static const image
*prepare_lod(const image
*current
) {
2066 if (current
== NULL
)
2069 return current
->scale_by_half("prepare_lod");
2073 * Prepare the next level of detail for weighted images.
2075 static const image
*prepare_lod(const image
*current
, const image
*weights
) {
2076 if (current
== NULL
)
2079 return current
->scale_by_half(weights
, "prepare_lod");
2083 * Prepare the next level of detail for definition maps.
2085 static const image
*prepare_lod_def(const image
*current
) {
2088 return current
->defined_scale_by_half("prepare_lod_def");
2092 * Initialize the scale cluster data structure.
2094 static struct scale_cluster
*init_clusters(int frame
, ale_real scale_factor
,
2095 const image
*input_frame
, unsigned int steps
,
2096 int *local_ax_count
) {
2099 * Allocate memory for the array.
2102 struct scale_cluster
*scale_clusters
=
2103 (struct scale_cluster
*) malloc(steps
* sizeof(struct scale_cluster
));
2105 assert (scale_clusters
);
2107 if (!scale_clusters
)
2108 ui::get()->memory_error("alignment");
2111 * Prepare images and exclusion regions for the highest level
2115 scale_clusters
[0].accum
= reference_image
;
2117 ui::get()->constructing_lod_clusters(0.0);
2118 scale_clusters
[0].input_scale
= scale_factor
;
2119 if (scale_factor
< 1.0 && interpolant
== NULL
)
2120 scale_clusters
[0].input
= input_frame
->scale(scale_factor
, "alignment");
2122 scale_clusters
[0].input
= input_frame
;
2124 scale_clusters
[0].defined
= reference_defined
;
2125 scale_clusters
[0].aweight
= alignment_weights_const
;
2126 scale_clusters
[0].ax_parameters
= filter_ax_parameters(frame
, local_ax_count
);
2128 scale_ax_parameters(*local_ax_count
, scale_clusters
[0].ax_parameters
, scale_factor
,
2129 (scale_factor
< 1.0 && interpolant
== NULL
) ? scale_factor
: 1);
2132 * Prepare reduced-detail images and exclusion
2136 for (unsigned int step
= 1; step
< steps
; step
++) {
2137 ui::get()->constructing_lod_clusters(step
);
2138 scale_clusters
[step
].accum
= prepare_lod(scale_clusters
[step
- 1].accum
, scale_clusters
[step
- 1].aweight
);
2139 scale_clusters
[step
].defined
= prepare_lod_def(scale_clusters
[step
- 1].defined
);
2140 scale_clusters
[step
].aweight
= prepare_lod(scale_clusters
[step
- 1].aweight
);
2141 scale_clusters
[step
].ax_parameters
2142 = copy_ax_parameters(*local_ax_count
, scale_clusters
[step
- 1].ax_parameters
);
2144 double sf
= scale_clusters
[step
- 1].input_scale
/ 2;
2145 scale_clusters
[step
].input_scale
= sf
;
2147 if (sf
>= 1.0 || interpolant
!= NULL
) {
2148 scale_clusters
[step
].input
= scale_clusters
[step
- 1].input
;
2149 scale_ax_parameters(*local_ax_count
, scale_clusters
[step
].ax_parameters
, 0.5, 1);
2150 } else if (sf
> 0.5) {
2151 scale_clusters
[step
].input
= scale_clusters
[step
- 1].input
->scale(sf
, "alignment");
2152 scale_ax_parameters(*local_ax_count
, scale_clusters
[step
].ax_parameters
, 0.5, sf
);
2154 scale_clusters
[step
].input
= scale_clusters
[step
- 1].input
->scale(0.5, "alignment");
2155 scale_ax_parameters(*local_ax_count
, scale_clusters
[step
].ax_parameters
, 0.5, 0.5);
2159 return scale_clusters
;
2163 * Destroy the first element in the scale cluster data structure.
2165 static void final_clusters(struct scale_cluster
*scale_clusters
, ale_real scale_factor
,
2166 unsigned int steps
) {
2168 if (scale_clusters
[0].input_scale
< 1.0)
2169 delete scale_clusters
[0].input
;
2171 free((void *)scale_clusters
[0].ax_parameters
);
2173 for (unsigned int step
= 1; step
< steps
; step
++) {
2174 delete scale_clusters
[step
].accum
;
2175 delete scale_clusters
[step
].defined
;
2176 delete scale_clusters
[step
].aweight
;
2178 if (scale_clusters
[step
].input_scale
< 1.0)
2179 delete scale_clusters
[step
].input
;
2181 free((void *)scale_clusters
[step
].ax_parameters
);
2184 free(scale_clusters
);
2188 * Calculate the centroid of a control point for the set of frames
2189 * having index lower than m. Divide by any scaling of the output.
2191 static point
unscaled_centroid(unsigned int m
, unsigned int p
) {
2194 point
point_sum(0, 0);
2195 ale_accum divisor
= 0;
2197 for(unsigned int j
= 0; j
< m
; j
++) {
2198 point pp
= cp_array
[p
][j
];
2201 point_sum
+= kept_t
[j
].transform_unscaled(pp
)
2202 / kept_t
[j
].scale();
2208 return point::undefined();
2210 return point_sum
/ divisor
;
2214 * Calculate centroid of this frame, and of all previous frames,
2215 * from points common to both sets.
2217 static void centroids(unsigned int m
, point
*current
, point
*previous
) {
2219 * Calculate the translation
2221 point
other_centroid(0, 0);
2222 point
this_centroid(0, 0);
2223 ale_pos divisor
= 0;
2225 for (unsigned int i
= 0; i
< cp_count
; i
++) {
2226 point other_c
= unscaled_centroid(m
, i
);
2227 point this_c
= cp_array
[i
][m
];
2229 if (!other_c
.defined() || !this_c
.defined())
2232 other_centroid
+= other_c
;
2233 this_centroid
+= this_c
;
2239 *current
= point::undefined();
2240 *previous
= point::undefined();
2244 *current
= this_centroid
/ divisor
;
2245 *previous
= other_centroid
/ divisor
;
2249 * Calculate the RMS error of control points for frame m, with
2250 * transformation t, against control points for earlier frames.
2252 static ale_accum
cp_rms_error(unsigned int m
, transformation t
) {
2256 ale_accum divisor
= 0;
2258 for (unsigned int i
= 0; i
< cp_count
; i
++)
2259 for (unsigned int j
= 0; j
< m
; j
++) {
2260 const point
*p
= cp_array
[i
];
2261 point p_ref
= kept_t
[j
].transform_unscaled(p
[j
]);
2262 point p_cur
= t
.transform_unscaled(p
[m
]);
2264 if (!p_ref
.defined() || !p_cur
.defined())
2267 err
+= p_ref
.lengthtosq(p_cur
);
2271 return sqrt(err
/ divisor
);
2275 * Implement new delta --follow semantics.
2277 * If we have a transformation T such that
2279 * prev_final == T(prev_init)
2283 * current_init_follow == T(current_init)
2285 * We can calculate T as follows:
2287 * T == prev_final(prev_init^-1)
2289 * Where ^-1 is the inverse operator.
2291 static transformation
follow(element_t
*element
, transformation offset
, int lod
) {
2293 transformation new_offset
= offset
;
2296 * Criteria for using following.
2299 if (!element
->old_is_default
&& !element
->is_default
&&
2300 default_initial_alignment_type
== 1) {
2302 * Ensure that the lod for the old initial and final
2303 * alignments are equal to the lod for the new initial
2307 ui::get()->following();
2309 element
->old_final_alignment
.rescale (1 / pow(2, lod
));
2310 element
->old_initial_alignment
.rescale(1 / pow(2, lod
- element
->old_lod
));
2312 for (offset
.set_current_index(0),
2313 element
->old_initial_alignment
.set_current_index(0),
2314 element
->old_final_alignment
.set_current_index(0),
2315 new_offset
.set_current_index(0);
2317 offset
.get_current_index() < _ma_card
;
2319 offset
.push_element(),
2320 new_offset
.push_element()) {
2322 if (alignment_class
== 0) {
2324 * Translational transformations
2327 ale_pos t0
= -element
->old_initial_alignment
.eu_get(0)
2328 + element
->old_final_alignment
.eu_get(0);
2329 ale_pos t1
= -element
->old_initial_alignment
.eu_get(1)
2330 + element
->old_final_alignment
.eu_get(1);
2332 new_offset
.eu_modify(0, t0
);
2333 new_offset
.eu_modify(1, t1
);
2335 } else if (alignment_class
== 1) {
2337 * Euclidean transformations
2340 ale_pos t2
= -element
->old_initial_alignment
.eu_get(2)
2341 + element
->old_final_alignment
.eu_get(2);
2343 new_offset
.eu_modify(2, t2
);
2345 point
p( offset
.scaled_height()/2 + offset
.eu_get(0) - element
->old_initial_alignment
.eu_get(0),
2346 offset
.scaled_width()/2 + offset
.eu_get(1) - element
->old_initial_alignment
.eu_get(1) );
2348 p
= element
->old_final_alignment
.transform_scaled(p
);
2350 new_offset
.eu_modify(0, p
[0] - offset
.scaled_height()/2 - offset
.eu_get(0));
2351 new_offset
.eu_modify(1, p
[1] - offset
.scaled_width()/2 - offset
.eu_get(1));
2353 } else if (alignment_class
== 2) {
2355 * Projective transformations
2360 p
[0] = element
->old_final_alignment
.transform_scaled(element
->old_initial_alignment
2361 . scaled_inverse_transform(offset
.get_current_element().transform_scaled(point( 0 , 0 ))));
2362 p
[1] = element
->old_final_alignment
.transform_scaled(element
->old_initial_alignment
2363 . scaled_inverse_transform(offset
.get_current_element().transform_scaled(point(offset
.scaled_height(), 0 ))));
2364 p
[2] = element
->old_final_alignment
.transform_scaled(element
->old_initial_alignment
2365 . scaled_inverse_transform(offset
.get_current_element().transform_scaled(point(offset
.scaled_height(), offset
.scaled_width()))));
2366 p
[3] = element
->old_final_alignment
.transform_scaled(element
->old_initial_alignment
2367 . scaled_inverse_transform(offset
.get_current_element().transform_scaled(point( 0 , offset
.scaled_width()))));
2369 new_offset
.gpt_set(p
);
2373 ui::get()->set_offset(offset
);
2379 static void test_global(diff_stat_t
*here
, scale_cluster si
, transformation t
,
2380 int local_ax_count
, int m
, ale_pos local_gs_mo
, ale_pos perturb
) {
2382 diff_stat_t
test(*here
);
2384 test
.diff(si
, test
.get_mc(), perturb
, t
, local_ax_count
, m
);
2386 unsigned int ovl
= overlap(si
, t
, local_ax_count
);
2388 if (ovl
>= local_gs_mo
&& test
.better()) {
2391 ui::get()->set_match(here
->get_error());
2392 ui::get()->set_offset(here
->get_offset());
2400 * Get the set of global transformations for a given density
2402 static void test_globals(diff_stat_t
*here
,
2403 scale_cluster si
, transformation t
, int local_gs
, ale_pos adj_p
,
2404 int local_ax_count
, int m
, ale_pos local_gs_mo
, ale_pos perturb
) {
2406 transformation offset
= t
;
2410 transformation offset_p
= offset
;
2412 if (!offset_p
.is_projective())
2413 offset_p
.eu_to_gpt();
2415 min
= max
= offset_p
.gpt_get(0);
2416 for (int p_index
= 1; p_index
< 4; p_index
++) {
2417 point p
= offset_p
.gpt_get(p_index
);
2428 point inner_min_t
= -min
;
2429 point inner_max_t
= -max
+ point(si
.accum
->height(), si
.accum
->width());
2430 point outer_min_t
= -max
+ point(adj_p
- 1, adj_p
- 1);
2431 point outer_max_t
= point(si
.accum
->height(), si
.accum
->width()) - point(adj_p
, adj_p
);
2433 if (local_gs
== 1 || local_gs
== 3 || local_gs
== 4 || local_gs
== 6) {
2439 for (ale_pos i
= inner_min_t
[0]; i
<= inner_max_t
[0]; i
+= adj_p
)
2440 for (ale_pos j
= inner_min_t
[1]; j
<= inner_max_t
[1]; j
+= adj_p
) {
2441 transformation test_t
= offset
;
2442 test_t
.translate(point(i
, j
));
2443 test_global(here
, si
, test_t
, local_ax_count
, m
, local_gs_mo
, perturb
);
2447 if (local_gs
== 2 || local_gs
== 3 || local_gs
== -1 || local_gs
== 6) {
2453 for (ale_pos i
= outer_min_t
[0]; i
<= outer_max_t
[0]; i
+= adj_p
)
2454 for (ale_pos j
= outer_min_t
[1]; j
< inner_min_t
[1]; j
+= adj_p
) {
2455 transformation test_t
= offset
;
2456 test_t
.translate(point(i
, j
));
2457 test_global(here
, si
, test_t
, local_ax_count
, m
, local_gs_mo
, perturb
);
2459 for (ale_pos i
= outer_min_t
[0]; i
<= outer_max_t
[0]; i
+= adj_p
)
2460 for (ale_pos j
= outer_max_t
[1]; j
> inner_max_t
[1]; j
-= adj_p
) {
2461 transformation test_t
= offset
;
2462 test_t
.translate(point(i
, j
));
2463 test_global(here
, si
, test_t
, local_ax_count
, m
, local_gs_mo
, perturb
);
2465 for (ale_pos i
= outer_min_t
[0]; i
< inner_min_t
[0]; i
+= adj_p
)
2466 for (ale_pos j
= outer_min_t
[1]; j
<= outer_max_t
[1]; j
+= adj_p
) {
2467 transformation test_t
= offset
;
2468 test_t
.translate(point(i
, j
));
2469 test_global(here
, si
, test_t
, local_ax_count
, m
, local_gs_mo
, perturb
);
2471 for (ale_pos i
= outer_max_t
[0]; i
> inner_max_t
[0]; i
-= adj_p
)
2472 for (ale_pos j
= outer_min_t
[1]; j
<= outer_max_t
[1]; j
+= adj_p
) {
2473 transformation test_t
= offset
;
2474 test_t
.translate(point(i
, j
));
2475 test_global(here
, si
, test_t
, local_ax_count
, m
, local_gs_mo
, perturb
);
2480 static void get_translational_set(std::vector
<transformation
> *set
,
2481 transformation t
, ale_pos adj_p
) {
2485 transformation offset
= t
;
2486 transformation test_t
;
2488 for (int i
= 0; i
< 2; i
++)
2489 for (adj_s
= -adj_p
; adj_s
<= adj_p
; adj_s
+= 2 * adj_p
) {
2493 test_t
.translate(i
? point(adj_s
, 0) : point(0, adj_s
));
2495 set
->push_back(test_t
);
2499 static int threshold_ok(ale_accum error
) {
2500 if ((1 - error
) * 100 >= match_threshold
)
2503 if (!(match_threshold
>= 0))
2510 * Align frame m against the reference.
2512 * XXX: the transformation class currently combines ordinary
2513 * transformations with scaling. This is somewhat convenient for
2514 * some things, but can also be confusing. This method, _align(), is
2515 * one case where special care must be taken to ensure that the scale
2516 * is always set correctly (by using the 'rescale' method).
2518 static diff_stat_t
_align(int m
, int local_gs
, element_t
*element
) {
2520 const image
*input_frame
= element
->input_frame
;
2523 * Local upper/lower data, possibly dependent on image
2527 ale_pos local_lower
, local_upper
, local_gs_mo
;
2530 * Select the minimum dimension as the reference.
2533 ale_pos reference_size
= input_frame
->height();
2534 if (input_frame
->width() < reference_size
)
2535 reference_size
= input_frame
->width();
2536 ale_pos reference_area
= input_frame
->height()
2537 * input_frame
->width();
2539 if (perturb_lower_percent
)
2540 local_lower
= perturb_lower
2545 local_lower
= perturb_lower
;
2547 if (perturb_upper_percent
)
2548 local_upper
= perturb_upper
2553 local_upper
= perturb_upper
;
2555 local_upper
= pow(2, floor(log(local_upper
) / log(2)));
2558 local_gs_mo
= _gs_mo
2563 local_gs_mo
= _gs_mo
;
2566 * Logarithms aren't exact, so we divide repeatedly to discover
2567 * how many steps will occur, and pass this information to the
2572 double step_variable
= local_upper
;
2573 while (step_variable
>= local_lower
) {
2578 ui::get()->set_steps(step_count
);
2580 ale_pos perturb
= local_upper
;
2584 kept_t
[latest
] = latest_t
;
2585 kept_ok
[latest
] = latest_ok
;
2589 * Maximum level-of-detail. Use a level of detail at most
2590 * 2^lod_diff finer than the adjustment resolution. lod_diff
2591 * is a synonym for lod_max.
2594 const int lod_diff
= lod_max
;
2597 * Determine how many levels of detail should be prepared.
2600 unsigned int steps
= (perturb
> pow(2, lod_max
))
2601 ? (unsigned int) (log(perturb
) / log(2)) - lod_max
+ 1 : 1;
2605 * Prepare multiple levels of detail.
2609 struct scale_cluster
*scale_clusters
= init_clusters(m
,
2610 scale_factor
, input_frame
, steps
,
2614 * Initialize variables used in the main loop.
2620 * Initialize the default initial transform
2623 if (default_initial_alignment_type
== 0) {
2626 * Follow the transformation of the original frame,
2627 * setting new image dimensions.
2630 // element->default_initial_alignment = orig_t;
2631 element
->default_initial_alignment
.set_current_element(orig_t
.get_element(0));
2632 element
->default_initial_alignment
.set_dimensions(input_frame
);
2634 } else if (default_initial_alignment_type
== 1)
2637 * Follow previous transformation, setting new image
2641 element
->default_initial_alignment
.set_dimensions(input_frame
);
2646 element
->old_is_default
= element
->is_default
;
2649 * Scale default initial transform for lod
2652 element
->default_initial_alignment
.rescale(1 / pow(2, lod
));
2655 * Set the default transformation.
2658 transformation offset
= element
->default_initial_alignment
;
2661 * Load any file-specified transformations
2664 for (offset
.set_current_index(0);
2665 offset
.get_current_index() < _ma_card
;
2666 offset
.push_element()) {
2668 offset
= tload_next(tload
, alignment_class
== 2,
2670 &element
->is_default
,
2671 offset
.get_current_index() == 0);
2675 offset
.set_current_index(0);
2677 ui::get()->set_offset(offset
);
2682 * Apply following logic
2685 transformation new_offset
= follow(element
, offset
, lod
);
2687 new_offset
.set_current_index(0);
2689 element
->old_initial_alignment
= offset
;
2690 element
->old_lod
= lod
;
2691 offset
= new_offset
;
2694 element
->old_initial_alignment
= offset
;
2695 element
->old_lod
= lod
;
2698 struct scale_cluster si
= scale_clusters
[lod
];
2699 ale_pos _mc_arg
= (_mc
> 0) ? (_mc
* pow(2, 2 * lod
))
2702 ui::get()->alignment_monte_carlo_parameter(_mc_arg
);
2705 * Projective adjustment value
2708 ale_pos adj_p
= (perturb
>= pow(2, lod_diff
))
2709 ? pow(2, lod_diff
) : (double) perturb
;
2712 * Orientational adjustment value in degrees.
2714 * Since rotational perturbation is now specified as an
2715 * arclength, we have to convert.
2718 ale_pos adj_o
= 2 * perturb
2719 / sqrt(pow(scale_clusters
[0].input
->height(), 2)
2720 + pow(scale_clusters
[0].input
->width(), 2))
2725 * Barrel distortion adjustment value
2728 ale_pos adj_b
= perturb
* bda_mult
;
2731 * Global search overlap requirements.
2734 local_gs_mo
/= pow(pow(2, lod
), 2);
2737 * Pre-alignment exposure adjustment
2740 if (_exp_register
) {
2741 ui::get()->exposure_1();
2742 transformation o
= offset
;
2743 for (int k
= lod
; k
> 0; k
--)
2745 set_exposure_ratio(m
, scale_clusters
[0], o
, local_ax_count
, 0);
2749 * Alignment statistics.
2752 diff_stat_t
here(_mc_arg
);
2755 * Current difference (error) value
2758 ui::get()->prematching();
2759 here
.diff(si
, _mc_arg
, perturb
, offset
, local_ax_count
, m
);
2760 ui::get()->set_match(here
.get_error());
2763 * Current and modified barrel distortion parameters
2766 ale_pos current_bd
[BARREL_DEGREE
];
2767 ale_pos modified_bd
[BARREL_DEGREE
];
2768 offset
.bd_get(current_bd
);
2769 offset
.bd_get(modified_bd
);
2772 * Translational global search step
2775 if (perturb
>= local_lower
&& local_gs
!= 0 && local_gs
!= 5
2776 && (local_gs
!= 6 || element
->is_default
)) {
2778 ui::get()->global_alignment(perturb
, lod
);
2779 ui::get()->gs_mo(local_gs_mo
);
2781 test_globals(&here
, si
, here
.get_offset(), local_gs
, adj_p
,
2782 local_ax_count
, m
, local_gs_mo
, perturb
);
2784 ui::get()->set_match(here
.get_error());
2785 ui::get()->set_offset(here
.get_offset());
2789 * Control point alignment
2792 if (local_gs
== 5) {
2794 transformation o
= here
.get_offset();
2796 for (int k
= lod
; k
> 0; k
--)
2800 * Determine centroid data
2803 point current
, previous
;
2804 centroids(m
, ¤t
, &previous
);
2806 if (current
.defined() && previous
.defined()) {
2808 o
.set_dimensions(input_frame
);
2809 o
.translate((previous
- current
) * o
.scale());
2814 * Determine rotation for alignment classes other than translation.
2817 ale_accum lowest_error
= cp_rms_error(m
, o
);
2819 ale_pos rot_lower
= 2 * local_lower
2820 / sqrt(pow(scale_clusters
[0].input
->height(), 2)
2821 + pow(scale_clusters
[0].input
->width(), 2))
2825 if (alignment_class
> 0)
2826 for (ale_pos rot
= 30; rot
> rot_lower
; rot
/= 2)
2827 for (ale_pos srot
= -rot
; srot
< rot
* 1.5; srot
+= rot
* 2) {
2828 int is_improved
= 1;
2829 while (is_improved
) {
2831 transformation test_t
= o
;
2833 * XXX: is this right?
2835 test_t
.rotate(current
* o
.scale(), srot
);
2836 ale_pos test_v
= cp_rms_error(m
, test_t
);
2838 if (test_v
< lowest_error
) {
2839 lowest_error
= test_v
;
2848 * Determine projective parameters through a local
2852 if (alignment_class
== 2) {
2853 ale_accum adj_p
= lowest_error
;
2855 if (adj_p
< local_lower
)
2856 adj_p
= local_lower
;
2858 while (adj_p
>= local_lower
) {
2859 transformation test_t
= o
;
2860 int is_improved
= 1;
2864 while (is_improved
) {
2867 for (int i
= 0; i
< 4; i
++)
2868 for (int j
= 0; j
< 2; j
++)
2869 for (adj_s
= -adj_p
; adj_s
<= adj_p
; adj_s
+= 2 * adj_p
) {
2873 if (perturb_type
== 0)
2874 test_t
.gpt_modify(j
, i
, adj_s
);
2875 else if (perturb_type
== 1)
2876 test_t
.gr_modify(j
, i
, adj_s
);
2880 test_v
= cp_rms_error(m
, test_t
);
2882 if (test_v
< lowest_error
) {
2883 lowest_error
= test_v
;
2895 set_exposure_ratio(m
, scale_clusters
[0], o
, local_ax_count
, 0);
2897 for (int k
= lod
; k
> 0; k
--)
2900 here
.diff(si
, here
.get_mc(), perturb
, o
, local_ax_count
, m
);
2902 ui::get()->set_match(here
.get_error());
2903 ui::get()->set_offset(here
.get_offset());
2907 * Announce perturbation size
2910 ui::get()->aligning(perturb
, lod
);
2913 * Run initial tests to get perturbation multipliers and error
2917 std::vector
<transformation
> t_set
;
2919 here
.get_perturb_set(&t_set
, adj_p
, adj_o
, adj_b
, current_bd
, modified_bd
);
2922 * Perturbation adjustment loop.
2925 int stable_count
= 0;
2927 while (perturb
>= local_lower
) {
2930 * Orientational adjustment value in degrees.
2932 * Since rotational perturbation is now specified as an
2933 * arclength, we have to convert.
2936 ale_pos adj_o
= 2 * perturb
2937 / sqrt(pow(scale_clusters
[0].input
->height(), 2)
2938 + pow(scale_clusters
[0].input
->width(), 2))
2943 * Barrel distortion adjustment value
2946 ale_pos adj_b
= perturb
* bda_mult
;
2948 diff_stat_t old_here
= here
;
2950 here
.perturb_test(perturb
, adj_p
, adj_o
, adj_b
, current_bd
, modified_bd
,
2953 if (here
.get_offset() == old_here
.get_offset())
2958 if (stable_count
== 3) {
2963 here
.calculate_element_region();
2965 if (here
.get_current_index() + 1 < _ma_card
) {
2966 here
.push_element();
2967 here
.make_element_nontrivial(adj_p
, adj_o
);
2968 element
->is_primary
= 0;
2971 here
.set_current_index(0);
2973 element
->is_primary
= 1;
2982 * Work with images twice as large
2986 si
= scale_clusters
[lod
];
2989 * Rescale the transforms.
2992 here
.rescale(2, si
);
2993 element
->default_initial_alignment
.rescale(2);
3003 ui::get()->alignment_perturbation_level(perturb
, lod
);
3008 ui::get()->set_match(here
.get_error());
3009 ui::get()->set_offset(here
.get_offset());
3012 here
.set_current_index(0);
3014 offset
= here
.get_offset();
3017 here
.rescale(pow(2, lod
), scale_clusters
[0]);
3018 element
->default_initial_alignment
.rescale(pow(2, lod
));
3022 * Post-alignment exposure adjustment
3025 if (_exp_register
== 1) {
3026 ui::get()->exposure_2();
3027 set_exposure_ratio(m
, scale_clusters
[0], offset
, local_ax_count
, 1);
3034 ui::get()->postmatching();
3035 offset
.use_full_support();
3036 here
.diff(scale_clusters
[0], here
.get_mc(), perturb
, offset
, local_ax_count
, m
);
3038 offset
.use_restricted_support();
3039 ui::get()->set_match(here
.get_error());
3042 * Free the level-of-detail structures
3045 final_clusters(scale_clusters
, scale_factor
, steps
);
3048 * Ensure that the match meets the threshold.
3051 if (threshold_ok(here
.get_error()) {
3053 * Update alignment variables
3056 element
->default_initial_alignment
= offset
;
3057 element
->old_final_alignment
= offset
;
3058 ui::get()->alignment_match_ok();
3059 } else if (local_gs
== 4) {
3062 * Align with outer starting points.
3066 * XXX: This probably isn't exactly the right thing to do,
3067 * since variables like old_initial_value have been overwritten.
3070 diff_stat_t nested_result
= _align(m
, -1, element
);
3072 if (threshold_ok(nested_result
.get_error())) {
3073 return nested_result
;
3074 } else if (nested_result
.get_error() < here
.get_error()) {
3075 here
= nested_result
;
3078 if (is_fail_default
)
3079 offset
= element
->default_initial_alignment
;
3081 ui::get()->set_match(here
.get_error());
3082 ui::get()->alignment_no_match();
3084 } else if (local_gs
== -1) {
3091 if (is_fail_default
)
3092 offset
= element
->default_initial_alignment
;
3094 ui::get()->alignment_no_match();
3098 * Write the tonal registration multiplier as a comment.
3101 pixel trm
= image_rw::exp(m
).get_multiplier();
3102 tsave_trm(tsave
, trm
[0], trm
[1], trm
[2]);
3105 * Save the transformation information
3108 for (offset
.set_current_index(0);
3109 offset
.get_current_index() < _ma_card
;
3110 offset
.push_element()) {
3112 tsave_next(tsave
, offset
, alignment_class
== 2,
3113 offset
.get_current_index() == 0);
3116 offset
.set_current_index(0);
3121 * Update match statistics.
3124 match_sum
+= (1 - here
.get_error()) * 100;
3133 * High-pass filter for frequency weights
3135 static void hipass(int rows
, int cols
, fftw_complex
*inout
) {
3136 for (int i
= 0; i
< rows
* vert_freq_cut
; i
++)
3137 for (int j
= 0; j
< cols
; j
++)
3138 for (int k
= 0; k
< 2; k
++)
3139 inout
[i
* cols
+ j
][k
] = 0;
3140 for (int i
= 0; i
< rows
; i
++)
3141 for (int j
= 0; j
< cols
* horiz_freq_cut
; j
++)
3142 for (int k
= 0; k
< 2; k
++)
3143 inout
[i
* cols
+ j
][k
] = 0;
3144 for (int i
= 0; i
< rows
; i
++)
3145 for (int j
= 0; j
< cols
; j
++)
3146 for (int k
= 0; k
< 2; k
++)
3147 if (i
/ (double) rows
+ j
/ (double) cols
< 2 * avg_freq_cut
)
3148 inout
[i
* cols
+ j
][k
] = 0;
3154 * Reset alignment weights
3156 static void reset_weights() {
3158 alignment_weights_const
= NULL
;
3160 if (alignment_weights
!= NULL
)
3161 delete alignment_weights
;
3163 alignment_weights
= NULL
;
3167 * Initialize alignment weights
3169 static void init_weights() {
3170 if (alignment_weights
!= NULL
)
3173 int rows
= reference_image
->height();
3174 int cols
= reference_image
->width();
3175 int colors
= reference_image
->depth();
3177 alignment_weights
= new image_ale_real(rows
, cols
,
3178 colors
, "alignment_weights");
3180 assert (alignment_weights
);
3182 for (int i
= 0; i
< rows
; i
++)
3183 for (int j
= 0; j
< cols
; j
++)
3184 alignment_weights
->set_pixel(i
, j
, pixel(1, 1, 1));
3188 * Update alignment weights with weight map
3190 static void map_update() {
3192 if (weight_map
== NULL
)
3197 point map_offset
= reference_image
->offset() - weight_map
->offset();
3199 int rows
= reference_image
->height();
3200 int cols
= reference_image
->width();
3202 for (int i
= 0; i
< rows
; i
++)
3203 for (int j
= 0; j
< cols
; j
++) {
3204 point map_weight_position
= map_offset
+ point(i
, j
);
3205 if (map_weight_position
[0] >= 0
3206 && map_weight_position
[1] >= 0
3207 && map_weight_position
[0] <= weight_map
->height() - 1
3208 && map_weight_position
[1] <= weight_map
->width() - 1)
3209 alignment_weights
->pix(i
, j
) *= weight_map
->get_bl(map_weight_position
);
3214 * Update alignment weights with an internal weight map, reflecting a
3215 * summation of certainty values. Use existing memory structures if
3216 * possible. This function updates alignment_weights_const; hence, it
3217 * should not be called prior to any functions that modify the
3218 * alignment_weights structure.
3220 static void imap_update() {
3221 if (alignment_weights
== NULL
) {
3222 alignment_weights_const
= reference_defined
;
3224 int rows
= reference_image
->height();
3225 int cols
= reference_image
->width();
3227 for (int i
= 0; i
< rows
; i
++)
3228 for (int j
= 0; j
< cols
; j
++)
3229 alignment_weights
->pix(i
, j
) *= reference_defined
->get_pixel(i
, j
);
3231 alignment_weights_const
= alignment_weights
;
3236 * Update alignment weights with algorithmic weights
3238 static void wmx_update() {
3241 static exposure
*exp_def
= new exposure_default();
3242 static exposure
*exp_bool
= new exposure_boolean();
3244 if (wmx_file
== NULL
|| wmx_exec
== NULL
|| wmx_defs
== NULL
)
3247 unsigned int rows
= reference_image
->height();
3248 unsigned int cols
= reference_image
->width();
3250 image_rw::write_image(wmx_file
, reference_image
);
3251 image_rw::write_image(wmx_defs
, reference_defined
, exp_bool
);
3254 int exit_status
= 1;
3256 execlp(wmx_exec
, wmx_exec
, wmx_file
, wmx_defs
, NULL
);
3257 ui::get()->exec_failure(wmx_exec
, wmx_file
, wmx_defs
);
3263 ui::get()->fork_failure("d2::align");
3265 image
*wmx_weights
= image_rw::read_image(wmx_file
, exp_def
);
3267 if (wmx_weights
->height() != rows
|| wmx_weights
->width() != cols
)
3268 ui::get()->error("algorithmic weighting must not change image size");
3270 if (alignment_weights
== NULL
)
3271 alignment_weights
= wmx_weights
;
3273 for (unsigned int i
= 0; i
< rows
; i
++)
3274 for (unsigned int j
= 0; j
< cols
; j
++)
3275 alignment_weights
->pix(i
, j
) *= wmx_weights
->pix(i
, j
);
3280 * Update alignment weights with frequency weights
3282 static void fw_update() {
3284 if (horiz_freq_cut
== 0
3285 && vert_freq_cut
== 0
3286 && avg_freq_cut
== 0)
3290 * Required for correct operation of --fwshow
3293 assert (alignment_weights
== NULL
);
3295 int rows
= reference_image
->height();
3296 int cols
= reference_image
->width();
3297 int colors
= reference_image
->depth();
3299 alignment_weights
= new image_ale_real(rows
, cols
,
3300 colors
, "alignment_weights");
3302 fftw_complex
*inout
= (fftw_complex
*) fftw_malloc(sizeof(fftw_complex
) * rows
* cols
);
3306 fftw_plan pf
= fftw_plan_dft_2d(rows
, cols
,
3308 FFTW_FORWARD
, FFTW_ESTIMATE
);
3310 fftw_plan pb
= fftw_plan_dft_2d(rows
, cols
,
3312 FFTW_BACKWARD
, FFTW_ESTIMATE
);
3314 for (int k
= 0; k
< colors
; k
++) {
3315 for (int i
= 0; i
< rows
* cols
; i
++) {
3316 inout
[i
][0] = reference_image
->get_pixel(i
/ cols
, i
% cols
)[k
];
3321 hipass(rows
, cols
, inout
);
3324 for (int i
= 0; i
< rows
* cols
; i
++) {
3326 alignment_weights
->pix(i
/ cols
, i
% cols
)[k
] = fabs(inout
[i
][0] / (rows
* cols
));
3328 alignment_weights
->pix(i
/ cols
, i
% cols
)[k
] =
3329 sqrt(pow(inout
[i
][0] / (rows
* cols
), 2)
3330 + pow(inout
[i
][1] / (rows
* cols
), 2));
3335 fftw_destroy_plan(pf
);
3336 fftw_destroy_plan(pb
);
3339 if (fw_output
!= NULL
)
3340 image_rw::write_image(fw_output
, alignment_weights
);
3345 * Update alignment to frame N.
3347 static void update_to(int n
) {
3349 assert (n
<= latest
+ 1);
3352 static std::vector
<element_t
> elements
;
3359 * Handle the initial frame
3362 elements
[0].input_frame
= image_rw::open(n
);
3364 const image
*i
= elements
[0].input_frame
;
3366 transformation result
= alignment_class
== 2
3367 ? transformation::gpt_identity(i
, scale_factor
)
3368 : transformation::eu_identity(i
, scale_factor
);
3369 result
= tload_first(tload
, alignment_class
== 2, result
, &is_default
);
3370 tsave_first(tsave
, result
, alignment_class
== 2);
3373 kept_t
= new transformation
[image_rw::count()];
3374 kept_ok
= (int *) malloc(image_rw::count()
3379 if (!kept_t
|| !kept_ok
)
3380 ui::get()->memory_error("alignment");
3390 elements
[0].default_initial_alignment
= result
;
3396 for (int i
= latest
+ 1; i
<= n
; i
++) {
3400 * Handle supplemental frames.
3403 assert (reference
!= NULL
);
3405 ui::get()->set_arender_current();
3406 reference
->sync(i
- 1);
3407 ui::get()->clear_arender_current();
3408 reference_image
= reference
->get_image();
3409 reference_defined
= reference
->get_defined();
3416 if (certainty_weights
)
3417 imap_update(); /* Must be called after all other _updates */
3419 assert (reference_image
!= NULL
);
3420 assert (reference_defined
!= NULL
);
3422 elements
[j
].input_frame
= image_rw::open(i
);
3423 elements
[j
].is_primary
= 1;
3425 _align(i
, _gs
, &elements
[j
]);
3430 if (elements
.size() > _ma_card
)
3431 elements
.resize(_ma_card
);
3437 * Set the control point count
3439 static void set_cp_count(unsigned int n
) {
3440 assert (cp_array
== NULL
);
3443 cp_array
= (const point
**) malloc(n
* sizeof(const point
*));
3447 * Set control points.
3449 static void set_cp(unsigned int i
, const point
*p
) {
3456 static void exp_register() {
3461 * Register exposure only based on metadata
3463 static void exp_meta_only() {
3468 * Don't register exposure
3470 static void exp_noregister() {
3475 * Set alignment class to translation only. Only adjust the x and y
3476 * position of images. Do not rotate input images or perform
3477 * projective transformations.
3479 static void class_translation() {
3480 alignment_class
= 0;
3484 * Set alignment class to Euclidean transformations only. Adjust the x
3485 * and y position of images and the orientation of the image about the
3488 static void class_euclidean() {
3489 alignment_class
= 1;
3493 * Set alignment class to perform general projective transformations.
3494 * See the file gpt.h for more information about general projective
3497 static void class_projective() {
3498 alignment_class
= 2;
3502 * Set the default initial alignment to the identity transformation.
3504 static void initial_default_identity() {
3505 default_initial_alignment_type
= 0;
3509 * Set the default initial alignment to the most recently matched
3510 * frame's final transformation.
3512 static void initial_default_follow() {
3513 default_initial_alignment_type
= 1;
3517 * Perturb output coordinates.
3519 static void perturb_output() {
3524 * Perturb source coordinates.
3526 static void perturb_source() {
3531 * Frames under threshold align optimally
3533 static void fail_optimal() {
3534 is_fail_default
= 0;
3538 * Frames under threshold keep their default alignments.
3540 static void fail_default() {
3541 is_fail_default
= 1;
3545 * Align images with an error contribution from each color channel.
3548 channel_alignment_type
= 0;
3552 * Align images with an error contribution only from the green channel.
3553 * Other color channels do not affect alignment.
3555 static void green() {
3556 channel_alignment_type
= 1;
3560 * Align images using a summation of channels. May be useful when
3561 * dealing with images that have high frequency color ripples due to
3565 channel_alignment_type
= 2;
3569 * Error metric exponent
3572 static void set_metric_exponent(float me
) {
3573 metric_exponent
= me
;
3580 static void set_match_threshold(float mt
) {
3581 match_threshold
= mt
;
3585 * Perturbation lower and upper bounds.
3588 static void set_perturb_lower(ale_pos pl
, int plp
) {
3590 perturb_lower_percent
= plp
;
3593 static void set_perturb_upper(ale_pos pu
, int pup
) {
3595 perturb_upper_percent
= pup
;
3599 * Maximum rotational perturbation.
3602 static void set_rot_max(int rm
) {
3605 * Obtain the largest power of two not larger than rm.
3608 rot_max
= pow(2, floor(log(rm
) / log(2)));
3612 * Barrel distortion adjustment multiplier
3615 static void set_bda_mult(ale_pos m
) {
3620 * Barrel distortion maximum rate of change
3623 static void set_bda_rate(ale_pos m
) {
3631 static void set_lod_max(int lm
) {
3636 * Set the scale factor
3638 static void set_scale(ale_pos s
) {
3643 * Set reference rendering to align against
3645 static void set_reference(render
*r
) {
3650 * Set the interpolant
3652 static void set_interpolant(filter::scaled_filter
*f
) {
3657 * Set alignment weights image
3659 static void set_weight_map(const image
*i
) {
3664 * Set frequency cuts
3666 static void set_frequency_cut(double h
, double v
, double a
) {
3673 * Set algorithmic alignment weighting
3675 static void set_wmx(const char *e
, const char *f
, const char *d
) {
3682 * Show frequency weights
3684 static void set_fl_show(const char *filename
) {
3685 fw_output
= filename
;
3689 * Set transformation file loader.
3691 static void set_tload(tload_t
*tl
) {
3696 * Set transformation file saver.
3698 static void set_tsave(tsave_t
*ts
) {
3703 * Get match statistics for frame N.
3705 static int match(int n
) {
3719 * Message that old alignment data should be kept.
3721 static void keep() {
3722 assert (latest
== -1);
3727 * Get alignment for frame N.
3729 static transformation
of(int n
) {
3742 * Use Monte Carlo alignment sampling with argument N.
3744 static void mc(ale_pos n
) {
3748 static void mcd_limit(int n
) {
3753 * Set the certainty-weighted flag.
3755 static void certainty_weighted(int flag
) {
3756 certainty_weights
= flag
;
3760 * Set the global search type.
3762 static void gs(const char *type
) {
3763 if (!strcmp(type
, "local")) {
3765 } else if (!strcmp(type
, "inner")) {
3767 } else if (!strcmp(type
, "outer")) {
3769 } else if (!strcmp(type
, "all")) {
3771 } else if (!strcmp(type
, "central")) {
3773 } else if (!strcmp(type
, "defaults")) {
3775 } else if (!strcmp(type
, "points")) {
3779 ui::get()->error("bad global search type");
3784 * Multi-alignment contiguity
3786 static void ma_cont(double value
) {
3791 * Multi-alignment cardinality
3793 static void ma_card(unsigned int value
) {
3794 assert (value
>= 1);
3799 * Set the minimum overlap for global searching
3801 static void gs_mo(ale_pos value
, int _gs_mo_percent
) {
3803 gs_mo_percent
= _gs_mo_percent
;
3807 * Don't use Monte Carlo alignment sampling.
3809 static void no_mc() {
3814 * Set alignment exclusion regions
3816 static void set_exclusion(exclusion
*_ax_parameters
, int _ax_count
) {
3817 ax_count
= _ax_count
;
3818 ax_parameters
= _ax_parameters
;
3822 * Get match summary statistics.
3824 static ale_accum
match_summary() {
3825 return match_sum
/ match_count
;