src/align.c: Raw import align_element method from ALE.

[libale.git] / src / align.c

/*  
 *  Copyright 2008, 2009 David Hilvert <dhilvert@gmail.com>
 *
 *  This file is part of libale.
 *
 *  libale is free software: you can redistribute it and/or modify it under the
 *  terms of the GNU Affero General Public License as published by the Free
 *  Software Foundation, either version 3 of the License, or (at your option)
 *  any later version.
 *
 *  libale is distributed in the hope that it will be useful, but WITHOUT ANY
 *  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 *  FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for
 *  more details.
 *
 *  You should have received a copy of the GNU Affero General Public License
 *  along with libale.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "libale.h"

/*
 * API Implementation.
 */

#warning raw imported code should be revised for Libale.
#if 0

/* XXX: Raw import of code from ALE.
 */

/*
 * Original ALE header copyright notice gives '2002, 2004, 2007 David Hilvert'
 * for the following, but git-blame (correctly) gives 2008.
 */

/*
 * Alignment state
 *
 * This structure contains alignment state information.  The change
 * between the non-default old initial alignment and old final
 * alignment is used to adjust the non-default current initial
 * alignment.  If either the old or new initial alignment is a default
 * alignment, the old --follow semantics are preserved.
 */

class astate_t {
        ale_trans old_initial_alignment;
        ale_trans old_final_alignment;
        ale_trans default_initial_alignment;
        int old_is_default;
        std::vector<int> is_default;
        ale_image input_frame;

public:
        astate_t() : 
                        is_default(1) {

                old_initial_alignment = ale_new_trans(accel::context(), NULL);
                old_final_alignment = ale_new_trans(accel::context(), NULL);
                default_initial_alignment = ale_new_trans(accel::context(), NULL);

                input_frame = NULL;
                is_default[0] = 1;
                old_is_default = 1;
        }

        ale_image get_input_frame() const {
                return input_frame;
        }

        void set_is_default(unsigned int index, int value) {

                /*
                 * Expand the array, if necessary.
                 */
                if (index == is_default.size());
                        is_default.resize(index + 1);

                assert (index < is_default.size());
                is_default[index] = value;
        }

        int get_is_default(unsigned int index) {
                assert (index < is_default.size());
                return is_default[index];
        }

        ale_trans get_default() {
                return default_initial_alignment;
        }

        void set_default(ale_trans t) {
                default_initial_alignment = t;
        }
        
        void default_set_original_bounds(ale_image i) {
                ale_trans_set_original_bounds(default_initial_alignment, i);
        }

        void set_final(ale_trans t) {
                old_final_alignment = t;
        }

        void set_input_frame(ale_image i) {
                input_frame = i;
        }

        /*
         * Implement new delta --follow semantics.
         *
         * If we have a transformation T such that
         *
         *      prev_final == T(prev_init)
         *
         * Then we also have
         *
         *      current_init_follow == T(current_init)
         *
         * We can calculate T as follows:
         *
         *      T == prev_final(prev_init^-1)
         *
         * Where ^-1 is the inverse operator.
         */
        static trans_single follow(trans_single a, trans_single b, trans_single c) {
                trans_single cc = c;

                if (alignment_class == 0) {
                        /*
                         * Translational transformations
                         */

                        ale_pos t0 = -a.eu_get(0) +  b.eu_get(0);
                        ale_pos t1 = -a.eu_get(1) +  b.eu_get(1);

                        cc.eu_modify(0, t0);
                        cc.eu_modify(1, t1);

                } else if (alignment_class == 1) {
                        /*
                         * Euclidean transformations
                         */

                        ale_pos t2 = -a.eu_get(2) +  b.eu_get(2);

                        cc.eu_modify(2, t2);

                        point p( c.scaled_height()/2 + c.eu_get(0) - a.eu_get(0),
                                 c.scaled_width()/2 + c.eu_get(1) - a.eu_get(1) );

                        p = b.transform_scaled(p);

                        cc.eu_modify(0, p[0] - c.scaled_height()/2 - c.eu_get(0));
                        cc.eu_modify(1, p[1] - c.scaled_width()/2 - c.eu_get(1));
                        
                } else if (alignment_class == 2) {
                        /*
                         * Projective transformations
                         */

                        point p[4];

                        p[0] = b.transform_scaled(a
                             . scaled_inverse_transform(c.transform_scaled(point(      0        ,       0       ))));
                        p[1] = b.transform_scaled(a
                             . scaled_inverse_transform(c.transform_scaled(point(c.scaled_height(),       0       ))));
                        p[2] = b.transform_scaled(a
                             . scaled_inverse_transform(c.transform_scaled(point(c.scaled_height(), c.scaled_width()))));
                        p[3] = b.transform_scaled(a
                             . scaled_inverse_transform(c.transform_scaled(point(      0        , c.scaled_width()))));

                        cc.gpt_set(p);
                }

                return cc;
        }

        /* 
         * For multi-alignment following, we use the following approach, not
         * guaranteed to work with large changes in scene or perspective, but
         * which should be somewhat flexible:
         *
         * For 
         *
         *      t[][] calculated final alignments 
         *      s[][] alignments as loaded from file
         *      previous frame n
         *      current frame n+1
         *      fundamental (primary) 0
         *      non-fundamental (non-primary) m!=0
         *      parent element m'
         *      follow(a, b, c) applying the (a, b) delta T=b(a^-1) to c
         *
         * following in the case of missing file data might be generated by 
         *
         *      t[n+1][0] = t[n][0]
         *      t[n+1][m!=0] = follow(t[n][m'], t[n+1][m'], t[n][m])
         *
         * cases with all noted file data present might be generated by 
         *
         *      t[n+1][0] = follow(s[n][0], t[n][0], s[n+1][0])
         *      t[n+1][m!=0] = follow(s[n+1][m'], t[n+1][m'], s[n+1][m])
         *
         * For non-following behavior, or where assigning the above is
         * impossible, we assign the following default
         *
         *      t[n+1][0] = Identity
         *      t[n+1][m!=0] = t[n+1][m']
         */

        void init_frame_alignment_primary(transformation *offset, int lod, ale_pos perturb) {

                if (perturb > 0 && !old_is_default && !get_is_default(0)
                 && default_initial_alignment_type == 1) {

                        /*
                         * Apply following logic for the primary element.
                         */

                        ui::get()->following();

                        trans_single new_offset = follow(old_initial_alignment.get_element(0),
                                                        old_final_alignment.get_element(0), 
                                                        offset->get_element(0));

                        old_initial_alignment = *offset;

                        offset->set_element(0, new_offset);

                        ui::get()->set_offset(new_offset);
                } else {
                        old_initial_alignment = *offset;
                }

                is_default.resize(old_initial_alignment.stack_depth());
        }

        void init_frame_alignment_nonprimary(transformation *offset, 
                        int lod, ale_pos perturb, unsigned int index) {

                assert (index > 0);

                unsigned int parent_index = offset->parent_index(index);

                if (perturb > 0 
                 && !get_is_default(parent_index)
                 && !get_is_default(index)
                 && default_initial_alignment_type == 1) {

                        /*
                         * Apply file-based following logic for the
                         * given element.
                         */

                        ui::get()->following();

                        trans_single new_offset = follow(old_initial_alignment.get_element(parent_index),
                                                        offset->get_element(parent_index), 
                                                        offset->get_element(index));

                        old_initial_alignment.set_element(index, offset->get_element(index));
                        offset->set_element(index, new_offset);

                        ui::get()->set_offset(new_offset);

                        return;
                }

                offset->get_coordinate(parent_index);


                if (perturb > 0
                 && old_final_alignment.exists(offset->get_coordinate(parent_index))
                 && old_final_alignment.exists(offset->get_current_coordinate())
                 && default_initial_alignment_type == 1) {

                        /*
                         * Apply nonfile-based following logic for
                         * the given element.
                         */

                        ui::get()->following();

                        /*
                         * XXX: Although it is different, the below
                         * should be equivalent to the comment
                         * description.
                         */

                        trans_single a = old_final_alignment.get_element(offset->get_coordinate(parent_index));
                        trans_single b = old_final_alignment.get_element(offset->get_current_coordinate());
                        trans_single c = offset->get_element(parent_index);
                        
                        trans_single new_offset = follow(a, b, c);

                        offset->set_element(index, new_offset);
                        ui::get()->set_offset(new_offset);

                        return;
                }

                /*
                 * Handle other cases.
                 */

                if (get_is_default(index)) {
                        offset->set_element(index, offset->get_element(parent_index));
                        ui::get()->set_offset(offset->get_element(index));
                }
        }

        void init_default() {

                if (default_initial_alignment_type == 0) {
                        
                        /*
                         * Follow the transformation of the original frame,
                         * setting new image dimensions.
                         */

                        // astate->default_initial_alignment = orig_t;
                        default_initial_alignment.set_current_element(orig_t.get_element(0));
                        default_initial_alignment.set_dimensions(input_frame);

                } else if (default_initial_alignment_type == 1)

                        /*
                         * Follow previous transformation, setting new image
                         * dimensions.
                         */

                        default_initial_alignment.set_dimensions(input_frame);

                else
                        assert(0);

                old_is_default = get_is_default(0);
        }
};


/*
 * Original ALE header copyright notice gives '2002, 2004, 2007 David Hilvert'
 * git-blame (which does not recall prior to 2005 at the moment) additionally
 * gives changes for the following several methods in years 2005, 2006, 2007,
 * 2008, and 2009.
 */

static struct scale_cluster *init_clusters(int frame, ale_pos scale_factor,
                const image *input_frame, unsigned int steps,
                int *local_ax_count) {

        /*
         * Allocate memory for the array.
         */

        struct scale_cluster *scale_clusters = 
                (struct scale_cluster *) malloc(steps * sizeof(struct scale_cluster));

        assert (scale_clusters);

        if (!scale_clusters)
                ui::get()->memory_error("alignment");

        /*
         * Prepare images and exclusion regions for the highest level
         * of detail.  
         */

        scale_clusters[0].accum = reference_image;

        ui::get()->constructing_lod_clusters(0.0);
        scale_clusters[0].input_scale = scale_factor;
        if (scale_factor < 1.0 && interpolant == NULL)
                scale_clusters[0].input = input_frame->scale(scale_factor, "alignment");
        else
                scale_clusters[0].input = input_frame;

        scale_clusters[0].certainty = reference_defined;
        scale_clusters[0].aweight = alignment_weights;
        scale_clusters[0].ax_parameters = filter_ax_parameters(frame, local_ax_count);

        /*
         * Allocate and determine input frame certainty.
         */

        if (scale_clusters[0].input->get_bayer() != IMAGE_BAYER_NONE) {
                scale_clusters[0].input_certainty = new_image_bayer_ale_real(
                                scale_clusters[0].input->height(),
                                scale_clusters[0].input->width(),
                                scale_clusters[0].input->depth(),
                                scale_clusters[0].input->get_bayer());
        } else {
                scale_clusters[0].input_certainty = scale_clusters[0].input->clone("certainty");
        }

        for (unsigned int i = 0; i < scale_clusters[0].input_certainty->height(); i++)
        for (unsigned int j = 0; j < scale_clusters[0].input_certainty->width(); j++)
        for (unsigned int k = 0; k < 3; k++)
        if (scale_clusters[0].input->get_channels(i, j) & (1 << k))
                ((image *) scale_clusters[0].input_certainty)->set_chan(i, j, k,
                        scale_clusters[0].input->
                                exp().confidence(scale_clusters[0].input->get_pixel(i, j))[k]);

        scale_ax_parameters(*local_ax_count, scale_clusters[0].ax_parameters, scale_factor, 
                        (scale_factor < 1.0 && interpolant == NULL) ? scale_factor : (ale_pos) 1);

        init_nl_cluster(&(scale_clusters[0]));

        /*
         * Prepare reduced-detail images and exclusion
         * regions.
         */

        for (unsigned int step = 1; step < steps; step++) {
                ui::get()->constructing_lod_clusters(step);
                scale_clusters[step].accum = prepare_lod(scale_clusters[step - 1].accum);
                scale_clusters[step].certainty = prepare_lod_def(scale_clusters[step - 1].certainty);
                scale_clusters[step].aweight = prepare_lod_def(scale_clusters[step - 1].aweight);
                scale_clusters[step].ax_parameters 
                        = copy_ax_parameters(*local_ax_count, scale_clusters[step - 1].ax_parameters);

                double sf = scale_clusters[step - 1].input_scale / 2;
                scale_clusters[step].input_scale = sf;

                if (sf >= 1.0 || interpolant != NULL) {
                        scale_clusters[step].input = scale_clusters[step - 1].input;
                        scale_clusters[step].input_certainty = scale_clusters[step - 1].input_certainty;
                        scale_ax_parameters(*local_ax_count, scale_clusters[step].ax_parameters, 0.5, 1);
                } else if (sf > 0.5) {
                        scale_clusters[step].input = scale_clusters[step - 1].input->scale(sf, "alignment");
                        scale_clusters[step].input_certainty = scale_clusters[step - 1].input->scale(sf, "alignment", 1);
                        scale_ax_parameters(*local_ax_count, scale_clusters[step].ax_parameters, 0.5, sf);
                } else {
                        scale_clusters[step].input = scale_clusters[step - 1].input->scale(0.5, "alignment");
                        scale_clusters[step].input_certainty = scale_clusters[step - 1].input_certainty->scale(0.5, "alignment", 1);
                        scale_ax_parameters(*local_ax_count, scale_clusters[step].ax_parameters, 0.5, 0.5);
                }

                init_nl_cluster(&(scale_clusters[step]));
        }

        return scale_clusters;
}

static diff_stat_t _align_element(ale_pos perturb, ale_pos local_lower,
                scale_cluster *scale_clusters, diff_stat_t here,
                ale_pos adj_p, ale_pos adj_o, ale_pos adj_b, 
                ale_pos *current_bd, ale_pos *modified_bd,
                astate_t *astate, int lod, scale_cluster si) {

        /*
         * Run initial tests to get perturbation multipliers and error
         * centroids.
         */

        std::vector<d2::trans_single> t_set;

        here.get_perturb_set(&t_set, adj_p, adj_o, adj_b, current_bd, modified_bd);

        int stable_count = 0;

        while (perturb >= local_lower) {
                
                ui::get()->alignment_dims(scale_clusters[lod].accum->height(), scale_clusters[lod].accum->width(),
                                          scale_clusters[lod].input->height(), scale_clusters[lod].input->width());

                /*
                 * Orientational adjustment value in degrees.
                 *
                 * Since rotational perturbation is now specified as an
                 * arclength, we have to convert.
                 */

                ale_pos adj_o = 2 * (double) perturb 
                                  / sqrt(pow(scale_clusters[0].input->height(), 2)
                                       + pow(scale_clusters[0].input->width(),  2))
                                  * 180
                                  / M_PI;

                /*
                 * Barrel distortion adjustment value
                 */

                ale_pos adj_b = perturb * bda_mult;

                trans_single old_offset = here.get_offset();

                here.perturb_test(perturb, adj_p, adj_o, adj_b, current_bd, modified_bd,
                        stable_count);

                if (here.get_offset() == old_offset)
                        stable_count++;
                else
                        stable_count = 0;

                if (stable_count == 3) {

                        stable_count = 0;

                        perturb *= 0.5;

                        if (lod > 0 
                         && lod > lrint(log(perturb) / log(2)) - lod_preferred) {

                                /*
                                 * Work with images twice as large
                                 */

                                lod--;
                                si = scale_clusters[lod];

                                /* 
                                 * Rescale the transforms.
                                 */

                                ale_pos rescale_factor = (double) scale_factor
                                                       / (double) pow(2, lod)
                                                       / (double) here.get_offset().scale();

                                here.rescale(rescale_factor, si);

                        } else {
                                adj_p = perturb / pow(2, lod);
                        }

                        /*
                         * Announce changes
                         */

                        ui::get()->alignment_perturbation_level(perturb, lod);
                }

                ui::get()->set_match(here.get_error());
                ui::get()->set_offset(here.get_offset());
        }

        if (lod > 0) {
                ale_pos rescale_factor = (double) scale_factor
                                       / (double) here.get_offset().scale();

                here.rescale(rescale_factor, scale_clusters[0]);
        }

        return here;
}

/*
 * Align frame m against the reference.
 *
 * XXX: the transformation class currently combines ordinary
 * transformations with scaling.  This is somewhat convenient for
 * some things, but can also be confusing.  This method, _align(), is
 * one case where special care must be taken to ensure that the scale
 * is always set correctly (by using the 'rescale' method).
 */
static diff_stat_multi _align(int m, int local_gs, astate_t *astate) {

        const image *input_frame = astate->get_input_frame();

        /*
         * Local upper/lower data, possibly dependent on image
         * dimensions.
         */

        ale_pos local_lower, local_upper;
        ale_accum local_gs_mo;

        /*
         * Select the minimum dimension as the reference.
         */

        ale_pos reference_size = input_frame->height();
        if (input_frame->width() < reference_size)
                reference_size = input_frame->width();
        ale_accum reference_area = input_frame->height()
                                 * input_frame->width();

        if (perturb_lower_percent)
                local_lower = (double) perturb_lower
                            * (double) reference_size
                            * (double) 0.01
                            * (double) scale_factor;
        else
                local_lower = perturb_lower;

        if (perturb_upper_percent)
                local_upper = (double) perturb_upper
                            * (double) reference_size
                            * (double) 0.01
                            * (double) scale_factor;
        else
                local_upper = perturb_upper;

        local_upper = pow(2, floor(log(local_upper) / log(2)));

        if (gs_mo_percent)
                local_gs_mo = (double) _gs_mo
                            * (double) reference_area
                            * (double) 0.01
                            * (double) scale_factor;
        else
                local_gs_mo = _gs_mo;

        /*
         * Logarithms aren't exact, so we divide repeatedly to discover
         * how many steps will occur, and pass this information to the
         * user interface.
         */

        int step_count = 0;
        double step_variable = local_upper;
        while (step_variable >= local_lower) {
                step_variable /= 2;
                step_count++;
        }

        ale_pos perturb = local_upper;

        if (_keep) {
                kept_t[latest] = latest_t;
                kept_ok[latest] = latest_ok;
        }

        /*
         * Determine how many levels of detail should be prepared, by
         * calculating the initial (largest) value for the
         * level-of-detail variable.
         */

        int lod = lrint(log(perturb) / log(2)) - lod_preferred;

        if (lod < 0)
                lod = 0;

        while (lod > 0 && (reference_image->width() < pow(2, lod) * min_dimension
                        || reference_image->height() < pow(2, lod) * min_dimension))
                lod--;

        unsigned int steps = (unsigned int) lod + 1;

        /*
         * Prepare multiple levels of detail.
         */

        int local_ax_count;
        struct scale_cluster *scale_clusters = init_clusters(m,
                        scale_factor, input_frame, steps,
                        &local_ax_count);

        /*
         * Initialize the default initial transform
         */

        astate->init_default();

        /*
         * Set the default transformation.
         */

        transformation offset = astate->get_default();

        /*
         * Establish boundaries
         */

        offset.set_current_bounds(reference_image);

        ui::get()->alignment_degree_max(offset.get_coordinate(offset.stack_depth() - 1).degree);

        if (offset.stack_depth() == 1) {
                ui::get()->set_steps(step_count, 0);
        } else {
                ui::get()->set_steps(offset.get_coordinate(offset.stack_depth() - 1).degree + 1, 1);
        }

        /*
         * Load any file-specified transformations
         */

        for (unsigned int index = 0; index < offset.stack_depth(); index++) {
                int is_default = 1;
                unsigned int index_2;
                offset.set_current_index(index);

                offset = tload_next(tload, alignment_class == 2, 
                                offset, 
                                &is_default, offset.get_current_index() == 0);

                index_2 = offset.get_current_index();

                if (index_2 > index) {
                        for (unsigned int index_3 = index; index_3 < index_2; index_3++)
                                astate->set_is_default(index_3, 1);

                        index = index_2;
                }

                astate->set_is_default(index, is_default);
        }

        offset.set_current_index(0);

        astate->init_frame_alignment_primary(&offset, lod, perturb);

        /*
         * Control point alignment
         */

        if (local_gs == 5) {

                transformation o = offset;

                /*
                 * Determine centroid data
                 */

                point current, previous;
                centroids(m, &current, &previous);

                if (current.defined() && previous.defined()) {
                        o = orig_t;
                        o.set_dimensions(input_frame);
                        o.translate((previous - current) * o.scale());
                        current = previous;
                }

                /*
                 * Determine rotation for alignment classes other than translation.
                 */

                ale_pos lowest_error = cp_rms_error(m, o);

                ale_pos rot_lower = 2 * (double) local_lower
                                  / sqrt(pow(scale_clusters[0].input->height(), 2)
                                       + pow(scale_clusters[0].input->width(),  2))
                                  * 180
                                  / M_PI;

                if  (alignment_class > 0)
                for (double rot = 30; rot > rot_lower; rot /= 2) 
                for (double srot = -rot; srot < rot * 1.5; srot += rot * 2) {
                        int is_improved = 1;
                        while (is_improved) {
                                is_improved = 0;
                                transformation test_t = o;
                                /*
                                 * XXX: is this right?
                                 */
                                test_t.rotate(current * o.scale(), srot);
                                ale_pos test_v = cp_rms_error(m, test_t);

                                if (test_v < lowest_error) {
                                        lowest_error = test_v;
                                        o = test_t;
                                        srot += 3 * rot;
                                        is_improved = 1;
                                }
                        }
                }

                /*
                 * Determine projective parameters through a local
                 * minimum search.
                 */

                if (alignment_class == 2) {
                        ale_pos adj_p = lowest_error;

                        if (adj_p < local_lower)
                                adj_p = local_lower;

                        while (adj_p >= local_lower) {
                                transformation test_t = o;
                                int is_improved = 1;
                                ale_pos test_v;
                                ale_pos adj_s;

                                while (is_improved) {
                                        is_improved = 0;

                                        for (int i = 0; i < 4; i++)
                                        for (int j = 0; j < 2; j++)
                                        for (adj_s = -adj_p; adj_s <= adj_p; adj_s += 2 * adj_p) {

                                                test_t = o;

                                                if (perturb_type == 0)
                                                        test_t.gpt_modify(j, i, adj_s);
                                                else if (perturb_type == 1)
                                                        test_t.gr_modify(j, i, adj_s);
                                                else
                                                        assert(0);

                                                test_v = cp_rms_error(m, test_t);

                                                if (test_v < lowest_error) {
                                                        lowest_error = test_v;
                                                        o = test_t;
                                                        adj_s += 3 * adj_p;
                                                        is_improved = 1;
                                                }
                                        }
                                }
                                adj_p /= 2;
                        }
                }
        }

        /*
         * Pre-alignment exposure adjustment
         */

        if (_exp_register) {
                ui::get()->exposure_1();
                set_exposure_ratio(m, scale_clusters[0], offset, local_ax_count, 0);
        }

        /*
         * Scale transform for lod
         */

        for (int lod_ = 0; lod_ < lod; lod_++) {
                transformation s = offset;
                transformation t = offset;

                t.rescale(1 / (double) 2);

                if (!(t.scaled_height() > 0 && t.scaled_height() < s.scaled_height())
                 || !(t.scaled_width() > 0  && t.scaled_width() < s.scaled_width())) {
                        perturb /= pow(2, lod - lod_);
                        lod = lod_;
                        break;
                } else {
                        offset = t;
                }
        }

        ui::get()->set_offset(offset);

        struct scale_cluster si = scale_clusters[lod];

        /*
         * Projective adjustment value
         */

        ale_pos adj_p = perturb / pow(2, lod);

        /*
         * Orientational adjustment value in degrees.
         *
         * Since rotational perturbation is now specified as an
         * arclength, we have to convert.
         */

        ale_pos adj_o = (double) 2 * (double) perturb 
                          / sqrt(pow((double) scale_clusters[0].input->height(), (double) 2)
                               + pow((double) scale_clusters[0].input->width(),  (double) 2))
                          * (double) 180
                          / M_PI;

        /*
         * Barrel distortion adjustment value
         */

        ale_pos adj_b = perturb * bda_mult;

        /*
         * Global search overlap requirements.
         */

        local_gs_mo = (double) local_gs_mo / pow(pow(2, lod), 2);

        /*
         * Alignment statistics.
         */

        diff_stat_t here(offset.elem_bounds());

        /*
         * Current difference (error) value
         */

        ui::get()->prematching();
        here.diff(si, offset.get_current_element(), local_ax_count, m);
        ui::get()->set_match(here.get_error());

        /*
         * Current and modified barrel distortion parameters
         */

        ale_pos current_bd[BARREL_DEGREE];
        ale_pos modified_bd[BARREL_DEGREE];
        offset.bd_get(current_bd);
        offset.bd_get(modified_bd);

        /*
         * Translational global search step
         */

        if (perturb >= local_lower && local_gs != 0 && local_gs != 5
         && (local_gs != 6 || astate->get_is_default(0))) {
                
                ui::get()->global_alignment(perturb, lod);
                ui::get()->gs_mo(local_gs_mo);

                test_globals(&here, si, offset, local_gs, adj_p,
                                local_ax_count, m, local_gs_mo, perturb);

                ui::get()->set_match(here.get_error());
                ui::get()->set_offset(here.get_offset());
        }

        /*
         * Perturbation adjustment loop.  
         */

        offset.set_current_element(here.get_offset());

        for (unsigned int i = 0; i < offset.stack_depth(); i++) {

                ui::get()->aligning_element(i, offset.stack_depth());
        
                offset.set_current_index(i);

                ui::get()->start_multi_alignment_element(offset);

                ui::get()->set_offset(offset);

                if (i > 0) {
                        astate->init_frame_alignment_nonprimary(&offset, lod, perturb, i);

                        if (_exp_register == 1) {
                                ui::get()->exposure_1();
                                pixel_accum asum(0, 0, 0), bsum(0, 0, 0);
                                exposure_ratio_iterate eri(&asum, &bsum, scale_clusters[0], offset, local_ax_count, 0, 
                                        offset.elem_bounds().scale_to_bounds(scale_clusters[0].accum->height(),
                                                                             scale_clusters[0].accum->width()));

                                eri.run();
                                pixel_accum tr = asum / bsum;
                                ui::get()->exp_multiplier(tr[0], tr[1], tr[2]);
                                offset.set_tonal_multiplier(tr);
                        }
                }

                int e_lod = lod;
                int e_div = offset.get_current_coordinate().degree;
                ale_pos e_perturb = perturb;
                ale_pos e_adj_p = adj_p;
                ale_pos e_adj_b = adj_b;

                for (int d = 0; d < e_div; d++) {
                        e_adj_b = 0;
                        e_perturb *= 0.5;
                        if (e_lod > 0) {
                                e_lod--;
                        } else {
                                e_adj_p *= 0.5;
                        }
                }

                if (i > 0) {

                        d2::trans_multi::elem_bounds_t b = offset.elem_bounds();

                        for (int dim_satisfied = 0; e_lod > 0 && !dim_satisfied; ) {
                                int height = scale_clusters[e_lod].accum->height();
                                int width = scale_clusters[e_lod].accum->width();

                                d2::trans_multi::elem_bounds_int_t bi = b.scale_to_bounds(height, width);

                                dim_satisfied = bi.satisfies_min_dim(min_dimension);
                                
                                if (!dim_satisfied) {
                                        e_lod--;
                                        e_adj_p *= 2;
                                }
                        }

                        /*
                         * Scale transform for lod
                         */

                        for (int lod_ = 0; lod_ < e_lod; lod_++) {
                                trans_single s = offset.get_element(i);
                                trans_single t = offset.get_element(i);

                                t.rescale(1 / (double) 2);

                                if (!(t.scaled_height() > 0 && t.scaled_height() < s.scaled_height())
                                 || !(t.scaled_width() > 0  && t.scaled_width() < s.scaled_width())) {
                                        e_perturb /= pow(2, e_lod - lod_);
                                        e_lod = lod_;
                                        break;
                                } else {
                                        offset.set_element(i, t);
                                }
                        }

                        ui::get()->set_offset(offset);
                }

                /*
                 * Announce perturbation size
                 */

                ui::get()->aligning(e_perturb, e_lod);

                si = scale_clusters[e_lod];

                here.set_elem_bounds(offset.elem_bounds());

                here.diff(si, offset.get_current_element(), local_ax_count, m);

                here.confirm();

                here = check_ancestor_path(offset, si, here, local_ax_count, m);

                here = _align_element(e_perturb, local_lower, scale_clusters, 
                                here, e_adj_p, adj_o, e_adj_b, current_bd, modified_bd,
                                astate, e_lod, si);

                offset.rescale(here.get_offset().scale() / offset.scale());

                offset.set_current_element(here.get_offset());

                if (i > 0 && _exp_register == 1) {
                        if (ma_cert_satisfied(scale_clusters[0], offset, i)) {
                                ui::get()->exposure_2();
                                pixel_accum asum(0, 0, 0), bsum(0, 0, 0);
                                exposure_ratio_iterate eri(&asum, &bsum, scale_clusters[0], offset, local_ax_count, 1, 
                                        offset.elem_bounds().scale_to_bounds(scale_clusters[0].accum->height(),
                                                                             scale_clusters[0].accum->width()));

                                eri.run();
                                pixel_accum tr = asum / bsum;
                                ui::get()->exp_multiplier(tr[0], tr[1], tr[2]);
                                offset.set_tonal_multiplier(tr);
                        } else {
                                offset.set_tonal_multiplier(offset.get_element(offset.parent_index(i)).get_tonal_multiplier(point(0, 0)));
                        }
                } else if (_exp_register == 1) {
                        ui::get()->exposure_2();
                        set_exposure_ratio(m, scale_clusters[0], offset, local_ax_count, 1);
                }

                ui::get()->set_offset(offset);

                if (i + 1 == offset.stack_depth())
                        ui::get()->alignment_degree_complete(offset.get_coordinate(i).degree);
                else if (offset.get_coordinate(i).degree != offset.get_coordinate(i + 1).degree)
                        ui::get()->alignment_degree_complete(offset.get_coordinate(i + 1).degree);
        }

        offset.set_current_index(0);

        ui::get()->multi();
        offset.set_multi(reference_image, input_frame);

        /*
         * Recalculate error on whole frame.
         */

        ui::get()->postmatching();
        diff_stat_generic<transformation> multi_here(offset.elem_bounds());
        multi_here.diff(scale_clusters[0], offset, local_ax_count, m);
        ui::get()->set_match(multi_here.get_error());

        /*
         * Free the level-of-detail structures
         */

        final_clusters(scale_clusters, scale_factor, steps);

        /*
         * Ensure that the match meets the threshold.
         */

        if (threshold_ok(multi_here.get_error())) {
                /*
                 * Update alignment variables
                 */
                latest_ok = 1;
                astate->set_default(offset);
                astate->set_final(offset);
                ui::get()->alignment_match_ok();
        } else if (local_gs == 4) {

                /*
                 * Align with outer starting points.
                 */

                /*
                 * XXX: This probably isn't exactly the right thing to do,
                 * since variables like old_initial_value have been overwritten.
                 */

                diff_stat_multi nested_result = _align(m, -1, astate);

                if (threshold_ok(nested_result.get_error())) {
                        return nested_result;
                } else if (nested_result.get_error() < multi_here.get_error()) {
                        multi_here = nested_result;
                }

                if (is_fail_default)
                        offset = astate->get_default();

                ui::get()->set_match(multi_here.get_error());
                ui::get()->alignment_no_match();
        
        } else if (local_gs == -1) {
                
                latest_ok = 0;
                latest_t = offset;
                return multi_here;

        } else {
                if (is_fail_default)
                        offset = astate->get_default();
                latest_ok = 0;
                ui::get()->alignment_no_match();
        }

        /*
         * Write the tonal registration multiplier as a comment.
         */

        pixel trm = image_rw::exp(m).get_multiplier();
        tsave_trm(tsave, trm[0], trm[1], trm[2]);

        /*
         * Save the transformation information
         */

        for (unsigned int index = 0; index < offset.stack_depth(); index++) {
                offset.set_current_index(index);

                tsave_next(tsave, offset, alignment_class == 2, 
                                  offset.get_current_index() == 0);
        }

        offset.set_current_index(0);

        /*
         * Update match statistics.
         */

        match_sum += (1 - multi_here.get_error()) * (ale_accum) 100;
        match_count++;
        latest = m;
        latest_t = offset;

        return multi_here;
}

#endif

int ale_align(ale_image a, ale_image b, ale_trans start,
                ale_align_properties align_properties, ale_trans result) {
#warning function unfinished.
}