2 * Copyright 2002, 2004, 2005, 2006, 2007, 2008, 2009 David Hilvert <dhilvert@gmail.com>
4 * This file is part of libale.
6 * libale is free software: you can redistribute it and/or modify it under the
7 * terms of the GNU Affero General Public License as published by the Free
8 * Software Foundation, either version 3 of the License, or (at your option)
11 * libale is distributed in the hope that it will be useful, but WITHOUT ANY
12 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13 * FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
16 * You should have received a copy of the GNU Affero General Public License
17 * along with libale. If not, see <http://www.gnu.org/licenses/>.
27 * Types for scale clusters.
33 * Experimental non-linear scaling approach.
36 struct nl_scale_cluster
{
37 const image
*accum_max
;
38 const image
*accum_min
;
39 const image
*certainty_max
;
40 const image
*certainty_min
;
41 const image
*aweight_max
;
42 const image
*aweight_min
;
43 exclusion
*ax_parameters
;
46 const image
*input_certainty_max
;
47 const image
*input_certainty_min
;
48 const image
*input_max
;
49 const image
*input_min
;
53 struct scale_cluster
{
57 ale_exclusion_list ax_parameters
;
60 ale_image input_certainty
;
64 nl_scale_cluster
*nl_scale_clusters
;
69 static struct scale_cluster
*init_clusters(int frame
, ale_pos scale_factor
,
70 const image
*input_frame
, unsigned int steps
,
71 int *local_ax_count
) {
74 * Allocate memory for the array.
77 struct scale_cluster
*scale_clusters
=
78 (struct scale_cluster
*) malloc(steps
* sizeof(struct scale_cluster
));
80 assert (scale_clusters
);
83 ui::get()->memory_error("alignment");
86 * Prepare images and exclusion regions for the highest level
90 scale_clusters
[0].accum
= reference_image
;
92 ui::get()->constructing_lod_clusters(0.0);
93 scale_clusters
[0].input_scale
= scale_factor
;
94 if (scale_factor
< 1.0 && interpolant
== NULL
)
95 scale_clusters
[0].input
= input_frame
->scale(scale_factor
, "alignment");
97 scale_clusters
[0].input
= input_frame
;
99 scale_clusters
[0].certainty
= reference_defined
;
100 scale_clusters
[0].aweight
= alignment_weights
;
101 scale_clusters
[0].ax_parameters
= filter_ax_parameters(frame
, local_ax_count
);
104 * Allocate and determine input frame certainty.
107 if (scale_clusters
[0].input
->get_bayer() != IMAGE_BAYER_NONE
) {
108 scale_clusters
[0].input_certainty
= new_image_bayer_ale_real(
109 scale_clusters
[0].input
->height(),
110 scale_clusters
[0].input
->width(),
111 scale_clusters
[0].input
->depth(),
112 scale_clusters
[0].input
->get_bayer());
114 scale_clusters
[0].input_certainty
= scale_clusters
[0].input
->clone("certainty");
117 for (unsigned int i
= 0; i
< scale_clusters
[0].input_certainty
->height(); i
++)
118 for (unsigned int j
= 0; j
< scale_clusters
[0].input_certainty
->width(); j
++)
119 for (unsigned int k
= 0; k
< 3; k
++)
120 if (scale_clusters
[0].input
->get_channels(i
, j
) & (1 << k
))
121 ((image
*) scale_clusters
[0].input_certainty
)->set_chan(i
, j
, k
,
122 scale_clusters
[0].input
->
123 exp().confidence(scale_clusters
[0].input
->get_pixel(i
, j
))[k
]);
125 scale_ax_parameters(*local_ax_count
, scale_clusters
[0].ax_parameters
, scale_factor
,
126 (scale_factor
< 1.0 && interpolant
== NULL
) ? scale_factor
: (ale_pos
) 1);
128 init_nl_cluster(&(scale_clusters
[0]));
131 * Prepare reduced-detail images and exclusion
135 for (unsigned int step
= 1; step
< steps
; step
++) {
136 ui::get()->constructing_lod_clusters(step
);
137 scale_clusters
[step
].accum
= prepare_lod(scale_clusters
[step
- 1].accum
);
138 scale_clusters
[step
].certainty
= prepare_lod_def(scale_clusters
[step
- 1].certainty
);
139 scale_clusters
[step
].aweight
= prepare_lod_def(scale_clusters
[step
- 1].aweight
);
140 scale_clusters
[step
].ax_parameters
141 = copy_ax_parameters(*local_ax_count
, scale_clusters
[step
- 1].ax_parameters
);
143 double sf
= scale_clusters
[step
- 1].input_scale
/ 2;
144 scale_clusters
[step
].input_scale
= sf
;
146 if (sf
>= 1.0 || interpolant
!= NULL
) {
147 scale_clusters
[step
].input
= scale_clusters
[step
- 1].input
;
148 scale_clusters
[step
].input_certainty
= scale_clusters
[step
- 1].input_certainty
;
149 scale_ax_parameters(*local_ax_count
, scale_clusters
[step
].ax_parameters
, 0.5, 1);
150 } else if (sf
> 0.5) {
151 scale_clusters
[step
].input
= scale_clusters
[step
- 1].input
->scale(sf
, "alignment");
152 scale_clusters
[step
].input_certainty
= scale_clusters
[step
- 1].input
->scale(sf
, "alignment", 1);
153 scale_ax_parameters(*local_ax_count
, scale_clusters
[step
].ax_parameters
, 0.5, sf
);
155 scale_clusters
[step
].input
= scale_clusters
[step
- 1].input
->scale(0.5, "alignment");
156 scale_clusters
[step
].input_certainty
= scale_clusters
[step
- 1].input_certainty
->scale(0.5, "alignment", 1);
157 scale_ax_parameters(*local_ax_count
, scale_clusters
[step
].ax_parameters
, 0.5, 0.5);
160 init_nl_cluster(&(scale_clusters
[step
]));
163 return scale_clusters
;
166 static diff_stat_t
_align_element(ale_pos perturb
, ale_pos local_lower
,
167 scale_cluster
*scale_clusters
, diff_stat_t here
,
168 ale_pos adj_p
, ale_pos adj_o
, ale_pos adj_b
,
169 ale_pos
*current_bd
, ale_pos
*modified_bd
,
170 astate_t
*astate
, int lod
, scale_cluster si
) {
173 * Run initial tests to get perturbation multipliers and error
177 std::vector
<d2::trans_single
> t_set
;
179 here
.get_perturb_set(&t_set
, adj_p
, adj_o
, adj_b
, current_bd
, modified_bd
);
181 int stable_count
= 0;
183 while (perturb
>= local_lower
) {
185 ui::get()->alignment_dims(scale_clusters
[lod
].accum
->height(), scale_clusters
[lod
].accum
->width(),
186 scale_clusters
[lod
].input
->height(), scale_clusters
[lod
].input
->width());
189 * Orientational adjustment value in degrees.
191 * Since rotational perturbation is now specified as an
192 * arclength, we have to convert.
195 ale_pos adj_o
= 2 * (double) perturb
196 / sqrt(pow(scale_clusters
[0].input
->height(), 2)
197 + pow(scale_clusters
[0].input
->width(), 2))
202 * Barrel distortion adjustment value
205 ale_pos adj_b
= perturb
* bda_mult
;
207 trans_single old_offset
= here
.get_offset();
209 here
.perturb_test(perturb
, adj_p
, adj_o
, adj_b
, current_bd
, modified_bd
,
212 if (here
.get_offset() == old_offset
)
217 if (stable_count
== 3) {
224 && lod
> lrint(log(perturb
) / log(2)) - lod_preferred
) {
227 * Work with images twice as large
231 si
= scale_clusters
[lod
];
234 * Rescale the transforms.
237 ale_pos rescale_factor
= (double) scale_factor
238 / (double) pow(2, lod
)
239 / (double) here
.get_offset().scale();
241 here
.rescale(rescale_factor
, si
);
244 adj_p
= perturb
/ pow(2, lod
);
251 ui::get()->alignment_perturbation_level(perturb
, lod
);
254 ui::get()->set_match(here
.get_error());
255 ui::get()->set_offset(here
.get_offset());
259 ale_pos rescale_factor
= (double) scale_factor
260 / (double) here
.get_offset().scale();
262 here
.rescale(rescale_factor
, scale_clusters
[0]);
269 * Align frame m against the reference.
271 * XXX: the transformation class currently combines ordinary
272 * transformations with scaling. This is somewhat convenient for
273 * some things, but can also be confusing. This method, _align(), is
274 * one case where special care must be taken to ensure that the scale
275 * is always set correctly (by using the 'rescale' method).
277 static diff_stat_multi
_align(int m
, int local_gs
, astate_t
*astate
) {
279 const image
*input_frame
= astate
->get_input_frame();
282 * Local upper/lower data, possibly dependent on image
286 ale_pos local_lower
, local_upper
;
287 ale_accum local_gs_mo
;
290 * Select the minimum dimension as the reference.
293 ale_pos reference_size
= input_frame
->height();
294 if (input_frame
->width() < reference_size
)
295 reference_size
= input_frame
->width();
296 ale_accum reference_area
= input_frame
->height()
297 * input_frame
->width();
299 if (perturb_lower_percent
)
300 local_lower
= (double) perturb_lower
301 * (double) reference_size
303 * (double) scale_factor
;
305 local_lower
= perturb_lower
;
307 if (perturb_upper_percent
)
308 local_upper
= (double) perturb_upper
309 * (double) reference_size
311 * (double) scale_factor
;
313 local_upper
= perturb_upper
;
315 local_upper
= pow(2, floor(log(local_upper
) / log(2)));
318 local_gs_mo
= (double) _gs_mo
319 * (double) reference_area
321 * (double) scale_factor
;
323 local_gs_mo
= _gs_mo
;
326 * Logarithms aren't exact, so we divide repeatedly to discover
327 * how many steps will occur, and pass this information to the
332 double step_variable
= local_upper
;
333 while (step_variable
>= local_lower
) {
338 ale_pos perturb
= local_upper
;
341 kept_t
[latest
] = latest_t
;
342 kept_ok
[latest
] = latest_ok
;
346 * Determine how many levels of detail should be prepared, by
347 * calculating the initial (largest) value for the
348 * level-of-detail variable.
351 int lod
= lrint(log(perturb
) / log(2)) - lod_preferred
;
356 while (lod
> 0 && (reference_image
->width() < pow(2, lod
) * min_dimension
357 || reference_image
->height() < pow(2, lod
) * min_dimension
))
360 unsigned int steps
= (unsigned int) lod
+ 1;
363 * Prepare multiple levels of detail.
367 struct scale_cluster
*scale_clusters
= init_clusters(m
,
368 scale_factor
, input_frame
, steps
,
372 * Initialize the default initial transform
375 astate
->init_default();
378 * Set the default transformation.
381 transformation offset
= astate
->get_default();
384 * Establish boundaries
387 offset
.set_current_bounds(reference_image
);
389 ui::get()->alignment_degree_max(offset
.get_coordinate(offset
.stack_depth() - 1).degree
);
391 if (offset
.stack_depth() == 1) {
392 ui::get()->set_steps(step_count
, 0);
394 ui::get()->set_steps(offset
.get_coordinate(offset
.stack_depth() - 1).degree
+ 1, 1);
398 * Load any file-specified transformations
401 for (unsigned int index
= 0; index
< offset
.stack_depth(); index
++) {
403 unsigned int index_2
;
404 offset
.set_current_index(index
);
406 offset
= tload_next(tload
, alignment_class
== 2,
408 &is_default
, offset
.get_current_index() == 0);
410 index_2
= offset
.get_current_index();
412 if (index_2
> index
) {
413 for (unsigned int index_3
= index
; index_3
< index_2
; index_3
++)
414 astate
->set_is_default(index_3
, 1);
419 astate
->set_is_default(index
, is_default
);
422 offset
.set_current_index(0);
424 astate
->init_frame_alignment_primary(&offset
, lod
, perturb
);
427 * Control point alignment
432 transformation o
= offset
;
435 * Determine centroid data
438 point current
, previous
;
439 centroids(m
, ¤t
, &previous
);
441 if (current
.defined() && previous
.defined()) {
443 o
.set_dimensions(input_frame
);
444 o
.translate((previous
- current
) * o
.scale());
449 * Determine rotation for alignment classes other than translation.
452 ale_pos lowest_error
= cp_rms_error(m
, o
);
454 ale_pos rot_lower
= 2 * (double) local_lower
455 / sqrt(pow(scale_clusters
[0].input
->height(), 2)
456 + pow(scale_clusters
[0].input
->width(), 2))
460 if (alignment_class
> 0)
461 for (double rot
= 30; rot
> rot_lower
; rot
/= 2)
462 for (double srot
= -rot
; srot
< rot
* 1.5; srot
+= rot
* 2) {
464 while (is_improved
) {
466 transformation test_t
= o
;
468 * XXX: is this right?
470 test_t
.rotate(current
* o
.scale(), srot
);
471 ale_pos test_v
= cp_rms_error(m
, test_t
);
473 if (test_v
< lowest_error
) {
474 lowest_error
= test_v
;
483 * Determine projective parameters through a local
487 if (alignment_class
== 2) {
488 ale_pos adj_p
= lowest_error
;
490 if (adj_p
< local_lower
)
493 while (adj_p
>= local_lower
) {
494 transformation test_t
= o
;
499 while (is_improved
) {
502 for (int i
= 0; i
< 4; i
++)
503 for (int j
= 0; j
< 2; j
++)
504 for (adj_s
= -adj_p
; adj_s
<= adj_p
; adj_s
+= 2 * adj_p
) {
508 if (perturb_type
== 0)
509 test_t
.gpt_modify(j
, i
, adj_s
);
510 else if (perturb_type
== 1)
511 test_t
.gr_modify(j
, i
, adj_s
);
515 test_v
= cp_rms_error(m
, test_t
);
517 if (test_v
< lowest_error
) {
518 lowest_error
= test_v
;
531 * Pre-alignment exposure adjustment
535 ui::get()->exposure_1();
536 set_exposure_ratio(m
, scale_clusters
[0], offset
, local_ax_count
, 0);
540 * Scale transform for lod
543 for (int lod_
= 0; lod_
< lod
; lod_
++) {
544 transformation s
= offset
;
545 transformation t
= offset
;
547 t
.rescale(1 / (double) 2);
549 if (!(t
.scaled_height() > 0 && t
.scaled_height() < s
.scaled_height())
550 || !(t
.scaled_width() > 0 && t
.scaled_width() < s
.scaled_width())) {
551 perturb
/= pow(2, lod
- lod_
);
559 ui::get()->set_offset(offset
);
561 struct scale_cluster si
= scale_clusters
[lod
];
564 * Projective adjustment value
567 ale_pos adj_p
= perturb
/ pow(2, lod
);
570 * Orientational adjustment value in degrees.
572 * Since rotational perturbation is now specified as an
573 * arclength, we have to convert.
576 ale_pos adj_o
= (double) 2 * (double) perturb
577 / sqrt(pow((double) scale_clusters
[0].input
->height(), (double) 2)
578 + pow((double) scale_clusters
[0].input
->width(), (double) 2))
583 * Barrel distortion adjustment value
586 ale_pos adj_b
= perturb
* bda_mult
;
589 * Global search overlap requirements.
592 local_gs_mo
= (double) local_gs_mo
/ pow(pow(2, lod
), 2);
595 * Alignment statistics.
598 diff_stat_t
here(offset
.elem_bounds());
601 * Current difference (error) value
604 ui::get()->prematching();
605 here
.diff(si
, offset
.get_current_element(), local_ax_count
, m
);
606 ui::get()->set_match(here
.get_error());
609 * Current and modified barrel distortion parameters
612 ale_pos current_bd
[BARREL_DEGREE
];
613 ale_pos modified_bd
[BARREL_DEGREE
];
614 offset
.bd_get(current_bd
);
615 offset
.bd_get(modified_bd
);
618 * Translational global search step
621 if (perturb
>= local_lower
&& local_gs
!= 0 && local_gs
!= 5
622 && (local_gs
!= 6 || astate
->get_is_default(0))) {
624 ui::get()->global_alignment(perturb
, lod
);
625 ui::get()->gs_mo(local_gs_mo
);
627 test_globals(&here
, si
, offset
, local_gs
, adj_p
,
628 local_ax_count
, m
, local_gs_mo
, perturb
);
630 ui::get()->set_match(here
.get_error());
631 ui::get()->set_offset(here
.get_offset());
635 * Perturbation adjustment loop.
638 offset
.set_current_element(here
.get_offset());
640 for (unsigned int i
= 0; i
< offset
.stack_depth(); i
++) {
642 ui::get()->aligning_element(i
, offset
.stack_depth());
644 offset
.set_current_index(i
);
646 ui::get()->start_multi_alignment_element(offset
);
648 ui::get()->set_offset(offset
);
651 astate
->init_frame_alignment_nonprimary(&offset
, lod
, perturb
, i
);
653 if (_exp_register
== 1) {
654 ui::get()->exposure_1();
655 pixel_accum
asum(0, 0, 0), bsum(0, 0, 0);
656 exposure_ratio_iterate
eri(&asum
, &bsum
, scale_clusters
[0], offset
, local_ax_count
, 0,
657 offset
.elem_bounds().scale_to_bounds(scale_clusters
[0].accum
->height(),
658 scale_clusters
[0].accum
->width()));
661 pixel_accum tr
= asum
/ bsum
;
662 ui::get()->exp_multiplier(tr
[0], tr
[1], tr
[2]);
663 offset
.set_tonal_multiplier(tr
);
668 int e_div
= offset
.get_current_coordinate().degree
;
669 ale_pos e_perturb
= perturb
;
670 ale_pos e_adj_p
= adj_p
;
671 ale_pos e_adj_b
= adj_b
;
673 for (int d
= 0; d
< e_div
; d
++) {
685 d2::trans_multi::elem_bounds_t b
= offset
.elem_bounds();
687 for (int dim_satisfied
= 0; e_lod
> 0 && !dim_satisfied
; ) {
688 int height
= scale_clusters
[e_lod
].accum
->height();
689 int width
= scale_clusters
[e_lod
].accum
->width();
691 d2::trans_multi::elem_bounds_int_t bi
= b
.scale_to_bounds(height
, width
);
693 dim_satisfied
= bi
.satisfies_min_dim(min_dimension
);
695 if (!dim_satisfied
) {
702 * Scale transform for lod
705 for (int lod_
= 0; lod_
< e_lod
; lod_
++) {
706 trans_single s
= offset
.get_element(i
);
707 trans_single t
= offset
.get_element(i
);
709 t
.rescale(1 / (double) 2);
711 if (!(t
.scaled_height() > 0 && t
.scaled_height() < s
.scaled_height())
712 || !(t
.scaled_width() > 0 && t
.scaled_width() < s
.scaled_width())) {
713 e_perturb
/= pow(2, e_lod
- lod_
);
717 offset
.set_element(i
, t
);
721 ui::get()->set_offset(offset
);
725 * Announce perturbation size
728 ui::get()->aligning(e_perturb
, e_lod
);
730 si
= scale_clusters
[e_lod
];
732 here
.set_elem_bounds(offset
.elem_bounds());
734 here
.diff(si
, offset
.get_current_element(), local_ax_count
, m
);
738 here
= check_ancestor_path(offset
, si
, here
, local_ax_count
, m
);
740 here
= _align_element(e_perturb
, local_lower
, scale_clusters
,
741 here
, e_adj_p
, adj_o
, e_adj_b
, current_bd
, modified_bd
,
744 offset
.rescale(here
.get_offset().scale() / offset
.scale());
746 offset
.set_current_element(here
.get_offset());
748 if (i
> 0 && _exp_register
== 1) {
749 if (ma_cert_satisfied(scale_clusters
[0], offset
, i
)) {
750 ui::get()->exposure_2();
751 pixel_accum
asum(0, 0, 0), bsum(0, 0, 0);
752 exposure_ratio_iterate
eri(&asum
, &bsum
, scale_clusters
[0], offset
, local_ax_count
, 1,
753 offset
.elem_bounds().scale_to_bounds(scale_clusters
[0].accum
->height(),
754 scale_clusters
[0].accum
->width()));
757 pixel_accum tr
= asum
/ bsum
;
758 ui::get()->exp_multiplier(tr
[0], tr
[1], tr
[2]);
759 offset
.set_tonal_multiplier(tr
);
761 offset
.set_tonal_multiplier(offset
.get_element(offset
.parent_index(i
)).get_tonal_multiplier(point(0, 0)));
763 } else if (_exp_register
== 1) {
764 ui::get()->exposure_2();
765 set_exposure_ratio(m
, scale_clusters
[0], offset
, local_ax_count
, 1);
768 ui::get()->set_offset(offset
);
770 if (i
+ 1 == offset
.stack_depth())
771 ui::get()->alignment_degree_complete(offset
.get_coordinate(i
).degree
);
772 else if (offset
.get_coordinate(i
).degree
!= offset
.get_coordinate(i
+ 1).degree
)
773 ui::get()->alignment_degree_complete(offset
.get_coordinate(i
+ 1).degree
);
776 offset
.set_current_index(0);
779 offset
.set_multi(reference_image
, input_frame
);
782 * Recalculate error on whole frame.
785 ui::get()->postmatching();
786 diff_stat_generic
<transformation
> multi_here(offset
.elem_bounds());
787 multi_here
.diff(scale_clusters
[0], offset
, local_ax_count
, m
);
788 ui::get()->set_match(multi_here
.get_error());
791 * Free the level-of-detail structures
794 final_clusters(scale_clusters
, scale_factor
, steps
);
797 * Ensure that the match meets the threshold.
800 if (threshold_ok(multi_here
.get_error())) {
802 * Update alignment variables
805 astate
->set_default(offset
);
806 astate
->set_final(offset
);
807 ui::get()->alignment_match_ok();
808 } else if (local_gs
== 4) {
811 * Align with outer starting points.
815 * XXX: This probably isn't exactly the right thing to do,
816 * since variables like old_initial_value have been overwritten.
819 diff_stat_multi nested_result
= _align(m
, -1, astate
);
821 if (threshold_ok(nested_result
.get_error())) {
822 return nested_result
;
823 } else if (nested_result
.get_error() < multi_here
.get_error()) {
824 multi_here
= nested_result
;
828 offset
= astate
->get_default();
830 ui::get()->set_match(multi_here
.get_error());
831 ui::get()->alignment_no_match();
833 } else if (local_gs
== -1) {
841 offset
= astate
->get_default();
843 ui::get()->alignment_no_match();
847 * Write the tonal registration multiplier as a comment.
850 pixel trm
= image_rw::exp(m
).get_multiplier();
851 tsave_trm(tsave
, trm
[0], trm
[1], trm
[2]);
854 * Save the transformation information
857 for (unsigned int index
= 0; index
< offset
.stack_depth(); index
++) {
858 offset
.set_current_index(index
);
860 tsave_next(tsave
, offset
, alignment_class
== 2,
861 offset
.get_current_index() == 0);
864 offset
.set_current_index(0);
867 * Update match statistics.
870 match_sum
+= (1 - multi_here
.get_error()) * (ale_accum
) 100;
878 int ale_align(ale_image a
, ale_image b
, ale_trans start
,
879 ale_align_properties align_properties
, ale_trans result
) {
880 #warning function unfinished.