src/align.c

   1 /*
   2  *  Copyright 2002, 2004, 2005, 2006, 2007, 2008, 2009 David Hilvert <dhilvert@gmail.com>
   3  *
   4  *  This file is part of libale.
   5  *
   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)
   9  *  any later version.
  10  *
  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
  14  *  more details.
  15  *
  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/>.
  18 */
  19
  20 #include <math.h>
  21
  22 #include "libale.h"
  23
  24 /*
  25  * API Implementation.
  26  */
  27
  28 /*
  29  * Types for scale clusters.
  30  */
  31
  32 #if 0
  33
  34 /*
  35  * Experimental non-linear scaling approach.
  36  */
  37
  38 struct nl_scale_cluster {
  39         ale_image accum_max;
  40         ale_image accum_min;
  41         ale_image aweight_max;
  42         ale_image aweight_min;
  43
  44         ale_pos input_scale;
  45         ale_image input_certainty_max;
  46         ale_image input_certainty_min;
  47         ale_image input_max;
  48         ale_image input_min;
  49 };
  50 #endif
  51
  52 struct scale_cluster {
  53         ale_image accum;
  54         ale_image aweight;
  55
  56         ale_pos input_scale;
  57         ale_image input_certainty;
  58         ale_image input;
  59
  60 #if 0
  61         nl_scale_cluster *nl_scale_clusters;
  62 #endif
  63 };
  64
  65
  66 static struct scale_cluster *init_clusters(ale_context ac, int frame, ale_pos scale_factor,
  67                 ale_image input_frame, ale_trans input_trans, ale_image reference_image,
  68                 ale_image alignment_weights, unsigned int steps, ale_filter interpolant) {
  69
  70         int input_bayer = ale_trans_get_bayer(input_trans);
  71         ale_real input_gamma = ale_trans_get_gamma(input_trans);
  72         ale_certainty input_certainty = ale_trans_retain_certainty(input_trans);
  73         ale_real input_black = ale_trans_get_black(input_trans);
  74
  75         unsigned int step;
  76
  77         /*
  78          * Allocate memory for the array.
  79          */
  80
  81         struct scale_cluster *scale_clusters =
  82                 (struct scale_cluster *) malloc(steps * sizeof(struct scale_cluster));
  83
  84         if (!scale_clusters)
  85                 return NULL;
  86
  87         /*
  88          * Prepare images and exclusion regions for the highest level
  89          * of detail.
  90          */
  91
  92         scale_clusters[0].accum = reference_image;
  93
  94         ale_context_signal(ac, "ui_constructing_lod_clusters", 0.0);
  95
  96         scale_clusters[0].input_scale = scale_factor;
  97         if (scale_factor < 1.0 && interpolant == NULL)
  98                 scale_clusters[0].input = ale_scale_image(input_frame, scale_factor, ale_context_get_color_type(ac), input_gamma, input_bayer);
  99         else
 100                 scale_clusters[0].input = input_frame;
 101
 102         scale_clusters[0].aweight = alignment_weights;
 103
 104         /*
 105          * Allocate and determine input frame certainty.
 106          */
 107
 108         scale_clusters[0].input_certainty = ale_clone_image(scale_clusters[0].input, ale_context_get_certainty_type(ac));
 109
 110         ale_eval("MAP_PIXEL(%0I, p, CERTAINTY(%1t, p, %2f, GET_AFFINE_PIXEL(%3i, p, %4f)))",
 111                         scale_clusters[0].input_certainty, input_black, input_certainty,
 112                         scale_clusters[0].input, (scale_factor < 1.0 && interpolant == NULL) ? 1 : input_gamma);
 113
 114 #if 0
 115         init_nl_cluster(&(scale_clusters[0]));
 116 #endif
 117
 118         /*
 119          * Prepare reduced-detail images and exclusion
 120          * regions.
 121          */
 122
 123         for (step = 1; step < steps; step++) {
 124                 ale_context_signal(ac, "ui_constructing_lod_clusters", step);
 125                 scale_clusters[step].accum = ale_scale_image(scale_clusters[step - 1].accum, 0.5, ale_context_get_color_type(ac), 1, ALE_BAYER_NONE);
 126                 scale_clusters[step].aweight = ale_scale_image_wt(scale_clusters[step - 1].aweight, 0.5, ale_context_get_weight_type(ac), ALE_BAYER_NONE);
 127
 128                 ale_pos sf = scale_clusters[step - 1].input_scale;
 129
 130                 scale_clusters[step].input_scale = sf / 2;
 131
 132                 if (sf >= 2.0 || interpolant != NULL) {
 133                         scale_clusters[step].input = scale_clusters[step - 1].input;
 134                         scale_clusters[step].input_certainty = scale_clusters[step - 1].input_certainty;
 135                 } else if (sf >= 1.0) {
 136                         scale_clusters[step].input = ale_scale_image(scale_clusters[step - 1].input, sf / 2, ale_context_get_color_type(ac), input_gamma, input_bayer);
 137                         scale_clusters[step].input_certainty = ale_scale_image_wt(scale_clusters[step - 1].input_certainty, sf / 2, ale_context_get_certainty_type(ac), input_bayer);
 138                 } else {
 139                         scale_clusters[step].input = ale_scale_image(scale_clusters[step - 1].input, 0.5, ale_context_get_color_type(ac), 1, ALE_BAYER_NONE);
 140                         scale_clusters[step].input_certainty = ale_scale_image_wt(scale_clusters[step - 1].input_certainty, 0.5, ale_context_get_certainty_type(ac), ALE_BAYER_NONE);
 141                 }
 142
 143 #if 0
 144                 init_nl_cluster(&(scale_clusters[step]));
 145 #endif
 146         }
 147
 148         ale_certainty_release(input_certainty);
 149
 150         return scale_clusters;
 151 }
 152
 153 /*
 154  * Not-quite-symmetric difference function.  Determines the difference in areas
 155  * where the arrays overlap.  Uses the reference array's notion of pixel positions.
 156  *
 157  * For the purposes of determining the difference, this function divides each
 158  * pixel value by the corresponding image's average pixel magnitude, unless we
 159  * are:
 160  *
 161  * a) Extending the boundaries of the image, or
 162  *
 163  * b) following the previous frame's transform
 164  *
 165  * If we are doing monte-carlo pixel sampling for alignment, we
 166  * typically sample a subset of available pixels; otherwise, we sample
 167  * all pixels.
 168  *
 169  */
 170
 171 ALE_STRUCT(run, \
 172         ale_trans offset; \
 173         ale_accum result; \
 174         ale_accum divisor; \
 175         point max, min; \
 176         ale_accum centroid[2], centroid_divisor; \
 177         ale_accum de_centroid[2], de_centroid_v, de_sum; \
 178         )
 179
 180 static void run_init(run *r, ale_trans offset) {
 181
 182         r->result = 0;
 183         r->divisor = 0;
 184
 185         r->min = point_posinf(2);
 186         r->max = point_neginf(2);
 187
 188         r->centroid[0] = 0;
 189         r->centroid[1] = 0;
 190         r->centroid_divisor = 0;
 191
 192         r->de_centroid[0] = 0;
 193         r->de_centroid[1] = 0;
 194
 195         r->de_centroid_v = 0;
 196
 197         r->de_sum = 0;
 198
 199         r->offset = offset;
 200 }
 201
 202 static void run_add(run *r, const run *s) {
 203         int d;
 204
 205         r->result += s->result;
 206         r->divisor += s->divisor;
 207
 208         for (d = 0; d < 2; d++) {
 209                 if (r->min.x[d] > s->min.x[d])
 210                         r->min.x[d] = s->min.x[d];
 211                 if (r->max.x[d] < s->max.x[d])
 212                         r->max.x[d] = s->max.x[d];
 213                 r->centroid[d] += s->centroid[d];
 214                 r->de_centroid[d] += s->de_centroid[d];
 215         }
 216
 217         r->centroid_divisor += s->centroid_divisor;
 218         r->de_centroid_v += s->de_centroid_v;
 219         r->de_sum += s->de_sum;
 220 }
 221
 222 static ale_accum run_get_error(run *r, ale_real metric_exponent) {
 223         return pow(r->result / r->divisor, 1/(ale_accum) metric_exponent);
 224 }
 225
 226 static void run_sample(int f, scale_cluster c, int i, int j, point t, point u,
 227                 const run &comparison) {
 228
 229         pixel pa = c.accum->get_pixel(i, j);
 230
 231         ale_real this_result[2] = { 0, 0 };
 232         ale_real this_divisor[2] = { 0, 0 };
 233
 234         pixel p[2];
 235         pixel weight[2];
 236         weight[0] = pixel(1, 1, 1);
 237         weight[1] = pixel(1, 1, 1);
 238
 239         pixel tm = offset.get_tonal_multiplier(point(i, j) + c.accum->offset());
 240
 241         if (interpolant != NULL) {
 242                 interpolant->filtered(i, j, &p[0], &weight[0], 1, f);
 243
 244                 if (weight[0].min_norm() > ale_real_weight_floor) {
 245                         p[0] /= weight[0];
 246                 } else {
 247                         return;
 248                 }
 249
 250         } else {
 251                 p[0] = c.input->get_bl(t);
 252         }
 253
 254         p[0] *= tm;
 255
 256         if (u.defined()) {
 257                 p[1] = c.input->get_bl(u);
 258                 p[1] *= tm;
 259         }
 260
 261
 262         /*
 263          * Handle certainty.
 264          */
 265
 266         if (certainty_weights == 1) {
 267
 268                 /*
 269                  * For speed, use arithmetic interpolation (get_bl(.))
 270                  * instead of geometric (get_bl(., 1))
 271                  */
 272
 273                 weight[0] *= c.input_certainty->get_bl(t);
 274                 if (u.defined())
 275                         weight[1] *= c.input_certainty->get_bl(u);
 276                 weight[0] *= c.certainty->get_pixel(i, j);
 277                 weight[1] *= c.certainty->get_pixel(i, j);
 278         }
 279
 280         if (c.aweight != NULL) {
 281                 weight[0] *= c.aweight->get_pixel(i, j);
 282                 weight[1] *= c.aweight->get_pixel(i, j);
 283         }
 284
 285         /*
 286          * Update sampling area statistics
 287          */
 288
 289         if (min[0] > i)
 290                 min[0] = i;
 291         if (min[1] > j)
 292                 min[1] = j;
 293         if (max[0] < i)
 294                 max[0] = i;
 295         if (max[1] < j)
 296                 max[1] = j;
 297
 298         centroid[0] += (weight[0][0] + weight[0][1] + weight[0][2]) * i;
 299         centroid[1] += (weight[0][0] + weight[0][1] + weight[0][2]) * j;
 300         centroid_divisor += (weight[0][0] + weight[0][1] + weight[0][2]);
 301
 302         /*
 303          * Determine alignment type.
 304          */
 305
 306         for (int m = 0; m < (u.defined() ? 2 : 1); m++)
 307         if (channel_alignment_type == 0) {
 308                 /*
 309                  * Align based on all channels.
 310                  */
 311
 312
 313                 for (int k = 0; k < 3; k++) {
 314                         ale_real achan = pa[k];
 315                         ale_real bchan = p[m][k];
 316
 317                         this_result[m] += weight[m][k] * pow(fabs(achan - bchan), metric_exponent);
 318                         this_divisor[m] += weight[m][k] * pow(achan > bchan ? achan : bchan, metric_exponent);
 319                 }
 320         } else if (channel_alignment_type == 1) {
 321                 /*
 322                  * Align based on the green channel.
 323                  */
 324
 325                 ale_real achan = pa[1];
 326                 ale_real bchan = p[m][1];
 327
 328                 this_result[m] = weight[m][1] * pow(fabs(achan - bchan), metric_exponent);
 329                 this_divisor[m] = weight[m][1] * pow(achan > bchan ? achan : bchan, metric_exponent);
 330         } else if (channel_alignment_type == 2) {
 331                 /*
 332                  * Align based on the sum of all channels.
 333                  */
 334
 335                 ale_real asum = 0;
 336                 ale_real bsum = 0;
 337                 ale_real wsum = 0;
 338
 339                 for (int k = 0; k < 3; k++) {
 340                         asum += pa[k];
 341                         bsum += p[m][k];
 342                         wsum += weight[m][k] / 3;
 343                 }
 344
 345                 this_result[m] = wsum * pow(fabs(asum - bsum), metric_exponent);
 346                 this_divisor[m] = wsum * pow(asum > bsum ? asum : bsum, metric_exponent);
 347         }
 348
 349         if (u.defined()) {
 350 //                                      ale_real de = fabs(this_result[0] / this_divisor[0]
 351 //                                                       - this_result[1] / this_divisor[1]);
 352                 ale_real de = fabs(this_result[0] - this_result[1]);
 353
 354                 de_centroid[0] += de * (ale_real) i;
 355                 de_centroid[1] += de * (ale_real) j;
 356
 357                 de_centroid_v += de * (ale_real) t.lengthto(u);
 358
 359                 de_sum += de;
 360         }
 361
 362         result += (this_result[0]);
 363         divisor += (this_divisor[0]);
 364 }
 365
 366 static void run_rescale(ale_pos scale) {
 367         offset.rescale(scale);
 368
 369         de_centroid[0] *= scale;
 370         de_centroid[1] *= scale;
 371         de_centroid_v *= scale;
 372 }
 373
 374 static point run_get_centroid() {
 375         point result = point(centroid[0] / centroid_divisor, centroid[1] / centroid_divisor);
 376
 377         assert (finite(centroid[0])
 378              && finite(centroid[1])
 379              && (result.defined() || centroid_divisor == 0));
 380
 381         return result;
 382 }
 383
 384 static point run_get_error_centroid() {
 385         point result = point(de_centroid[0] / de_sum, de_centroid[1] / de_sum);
 386         return result;
 387 }
 388
 389
 390 static ale_pos run_get_error_perturb() {
 391         ale_pos result = de_centroid_v / de_sum;
 392
 393         return result;
 394 }
 395
 396 template<class diff_trans>
 397 class diff_stat_generic {
 398
 399         transformation::elem_bounds_t elem_bounds;
 400
 401         /*
 402          * When non-empty, runs.front() is best, runs.back() is
 403          * testing.
 404          */
 405
 406         std::vector<run> runs;
 407
 408         /*
 409          * old_runs stores the latest available perturbation set for
 410          * each multi-alignment element.
 411          */
 412
 413         typedef int run_index;
 414         std::map<run_index, run> old_runs;
 415
 416         static void *diff_subdomain(void *args);
 417
 418         struct subdomain_args {
 419                 struct scale_cluster c;
 420                 std::vector<run> runs;
 421                 int ax_count;
 422                 int f;
 423                 int i_min, i_max, j_min, j_max;
 424                 int subdomain;
 425         };
 426
 427         struct scale_cluster si;
 428         int ax_count;
 429         int frame;
 430
 431         std::vector<ale_pos> perturb_multipliers;
 432
 433 public:
 434         void diff(struct scale_cluster c, const diff_trans &t,
 435                         int _ax_count, int f) {
 436
 437                 if (runs.size() == 2)
 438                         runs.pop_back();
 439
 440                 runs.push_back(run(t));
 441
 442                 si = c;
 443                 ax_count = _ax_count;
 444                 frame = f;
 445
 446                 ui::get()->d2_align_sample_start();
 447
 448                 if (interpolant != NULL) {
 449
 450                         /*
 451                          * XXX: This has not been tested, and may be completely
 452                          * wrong.
 453                          */
 454
 455                         transformation tt = transformation::eu_identity();
 456                         tt.set_current_element(t);
 457                         interpolant->set_parameters(tt, c.input, c.accum->offset());
 458                 }
 459
 460                 int N;
 461 #ifdef USE_PTHREAD
 462                 N = thread::count();
 463
 464                 pthread_t *threads = (pthread_t *) malloc(sizeof(pthread_t) * N);
 465                 pthread_attr_t *thread_attr = (pthread_attr_t *) malloc(sizeof(pthread_attr_t) * N);
 466
 467 #else
 468                 N = 1;
 469 #endif
 470
 471                 subdomain_args *args = new subdomain_args[N];
 472
 473                 transformation::elem_bounds_int_t b = elem_bounds.scale_to_bounds(c.accum->height(), c.accum->width());
 474
 475 //                      fprintf(stdout, "[%d %d] [%d %d]\n",
 476 //                              global_i_min, global_i_max, global_j_min, global_j_max);
 477
 478                 for (int ti = 0; ti < N; ti++) {
 479                         args[ti].c = c;
 480                         args[ti].runs = runs;
 481                         args[ti].ax_count = ax_count;
 482                         args[ti].f = f;
 483                         args[ti].i_min = b.imin + ((b.imax - b.imin) * ti) / N;
 484                         args[ti].i_max = b.imin + ((b.imax - b.imin) * (ti + 1)) / N;
 485                         args[ti].j_min = b.jmin;
 486                         args[ti].j_max = b.jmax;
 487                         args[ti].subdomain = ti;
 488 #ifdef USE_PTHREAD
 489                         pthread_attr_init(&thread_attr[ti]);
 490                         pthread_attr_setdetachstate(&thread_attr[ti], PTHREAD_CREATE_JOINABLE);
 491                         pthread_create(&threads[ti], &thread_attr[ti], diff_subdomain, &args[ti]);
 492 #else
 493                         diff_subdomain(&args[ti]);
 494 #endif
 495                 }
 496
 497                 for (int ti = 0; ti < N; ti++) {
 498 #ifdef USE_PTHREAD
 499                         pthread_join(threads[ti], NULL);
 500 #endif
 501                         runs.back().add(args[ti].runs.back());
 502                 }
 503
 504 #ifdef USE_PTHREAD
 505                 free(threads);
 506                 free(thread_attr);
 507 #endif
 508
 509                 delete[] args;
 510
 511                 ui::get()->d2_align_sample_stop();
 512
 513         }
 514
 515 private:
 516         void rediff() {
 517                 std::vector<diff_trans> t_array;
 518
 519                 for (unsigned int r = 0; r < runs.size(); r++) {
 520                         t_array.push_back(runs[r].offset);
 521                 }
 522
 523                 runs.clear();
 524
 525                 for (unsigned int r = 0; r < t_array.size(); r++)
 526                         diff(si, t_array[r], ax_count, frame);
 527         }
 528
 529
 530 public:
 531         int better() {
 532                 assert(runs.size() >= 2);
 533                 assert(runs[0].offset.scale() == runs[1].offset.scale());
 534
 535                 return (runs[1].get_error() < runs[0].get_error()
 536                      || (!finite(runs[0].get_error()) && finite(runs[1].get_error())));
 537         }
 538
 539         int better_defined() {
 540                 assert(runs.size() >= 2);
 541                 assert(runs[0].offset.scale() == runs[1].offset.scale());
 542
 543                 return (runs[1].get_error() < runs[0].get_error());
 544         }
 545
 546         diff_stat_generic(transformation::elem_bounds_t e)
 547                         : runs(), old_runs(), perturb_multipliers() {
 548                 elem_bounds = e;
 549         }
 550
 551         run_index get_run_index(unsigned int perturb_index) {
 552                 return perturb_index;
 553         }
 554
 555         run &get_run(unsigned int perturb_index) {
 556                 run_index index = get_run_index(perturb_index);
 557
 558                 assert(old_runs.count(index));
 559                 return old_runs[index];
 560         }
 561
 562         void rescale(ale_pos scale, scale_cluster _si) {
 563                 assert(runs.size() == 1);
 564
 565                 si = _si;
 566
 567                 runs[0].rescale(scale);
 568
 569                 rediff();
 570         }
 571
 572         ~diff_stat_generic() {
 573         }
 574
 575         diff_stat_generic &operator=(const diff_stat_generic &dst) {
 576                 /*
 577                  * Copy run information.
 578                  */
 579                 runs = dst.runs;
 580                 old_runs = dst.old_runs;
 581
 582                 /*
 583                  * Copy diff variables
 584                  */
 585                 si = dst.si;
 586                 ax_count = dst.ax_count;
 587                 frame = dst.frame;
 588                 perturb_multipliers = dst.perturb_multipliers;
 589                 elem_bounds = dst.elem_bounds;
 590
 591                 return *this;
 592         }
 593
 594         diff_stat_generic(const diff_stat_generic &dst) : runs(), old_runs(),
 595                                               perturb_multipliers() {
 596                 operator=(dst);
 597         }
 598
 599         void set_elem_bounds(transformation::elem_bounds_t e) {
 600                 elem_bounds = e;
 601         }
 602
 603         ale_accum get_result() {
 604                 assert(runs.size() == 1);
 605                 return runs[0].result;
 606         }
 607
 608         ale_accum get_divisor() {
 609                 assert(runs.size() == 1);
 610                 return runs[0].divisor;
 611         }
 612
 613         diff_trans get_offset() {
 614                 assert(runs.size() == 1);
 615                 return runs[0].offset;
 616         }
 617
 618         int operator!=(diff_stat_generic &param) {
 619                 return (get_error() != param.get_error());
 620         }
 621
 622         int operator==(diff_stat_generic &param) {
 623                 return !(operator!=(param));
 624         }
 625
 626         ale_pos get_error_perturb() {
 627                 assert(runs.size() == 1);
 628                 return runs[0].get_error_perturb();
 629         }
 630
 631         ale_accum get_error() const {
 632                 assert(runs.size() == 1);
 633                 return runs[0].get_error();
 634         }
 635
 636 public:
 637         /*
 638          * Get the set of transformations produced by a given perturbation
 639          */
 640         void get_perturb_set(std::vector<diff_trans> *set,
 641                         ale_pos adj_p, ale_pos adj_o, ale_pos adj_b,
 642                         ale_pos *current_bd, ale_pos *modified_bd,
 643                         std::vector<ale_pos> multipliers = std::vector<ale_pos>()) {
 644
 645                 assert(runs.size() == 1);
 646
 647                 diff_trans test_t(diff_trans::eu_identity());
 648
 649                 /*
 650                  * Translational or euclidean transformation
 651                  */
 652
 653                 for (unsigned int i = 0; i < 2; i++)
 654                 for (unsigned int s = 0; s < 2; s++) {
 655
 656                         if (!multipliers.size())
 657                                 multipliers.push_back(1);
 658
 659                         assert(finite(multipliers[0]));
 660
 661                         test_t = get_offset();
 662
 663                         // test_t.eu_modify(i, (s ? -adj_p : adj_p) * multipliers[0]);
 664                         test_t.translate((i ? point(1, 0) : point(0, 1))
 665                                        * (s ? -adj_p : adj_p)
 666                                        * multipliers[0]);
 667
 668                         test_t.snap(adj_p / 2);
 669
 670                         set->push_back(test_t);
 671                         multipliers.erase(multipliers.begin());
 672
 673                 }
 674
 675                 if (alignment_class > 0)
 676                 for (unsigned int s = 0; s < 2; s++) {
 677
 678                         if (!multipliers.size())
 679                                 multipliers.push_back(1);
 680
 681                         assert(multipliers.size());
 682                         assert(finite(multipliers[0]));
 683
 684                         if (!(adj_o * multipliers[0] < rot_max))
 685                                 return;
 686
 687                         ale_pos adj_s = (s ? 1 : -1) * adj_o * multipliers[0];
 688
 689                         test_t = get_offset();
 690
 691                         run_index ori = get_run_index(set->size());
 692                         point centroid = point::undefined();
 693
 694                         if (!old_runs.count(ori))
 695                                 ori = get_run_index(0);
 696
 697                         if (!centroid.finite() && old_runs.count(ori)) {
 698                                 centroid = old_runs[ori].get_error_centroid();
 699
 700                                 if (!centroid.finite())
 701                                         centroid = old_runs[ori].get_centroid();
 702
 703                                 centroid *= test_t.scale()
 704                                           / old_runs[ori].offset.scale();
 705                         }
 706
 707                         if (!centroid.finite() && !test_t.is_projective()) {
 708                                 test_t.eu_modify(2, adj_s);
 709                         } else if (!centroid.finite()) {
 710                                 centroid = point(si.input->height() / 2,
 711                                                  si.input->width() / 2);
 712
 713                                 test_t.rotate(centroid + si.accum->offset(),
 714                                                 adj_s);
 715                         } else {
 716                                 test_t.rotate(centroid + si.accum->offset(),
 717                                                 adj_s);
 718                         }
 719
 720                         test_t.snap(adj_p / 2);
 721
 722                         set->push_back(test_t);
 723                         multipliers.erase(multipliers.begin());
 724                 }
 725
 726                 if (alignment_class == 2) {
 727
 728                         /*
 729                          * Projective transformation
 730                          */
 731
 732                         for (unsigned int i = 0; i < 4; i++)
 733                         for (unsigned int j = 0; j < 2; j++)
 734                         for (unsigned int s = 0; s < 2; s++) {
 735
 736                                 if (!multipliers.size())
 737                                         multipliers.push_back(1);
 738
 739                                 assert(multipliers.size());
 740                                 assert(finite(multipliers[0]));
 741
 742                                 ale_pos adj_s = (s ? -1 : 1) * adj_p * multipliers [0];
 743
 744                                 test_t = get_offset();
 745
 746                                 if (perturb_type == 0)
 747                                         test_t.gpt_modify_bounded(j, i, adj_s, elem_bounds.scale_to_bounds(si.accum->height(), si.accum->width()));
 748                                 else if (perturb_type == 1)
 749                                         test_t.gr_modify(j, i, adj_s);
 750                                 else
 751                                         assert(0);
 752
 753                                 test_t.snap(adj_p / 2);
 754
 755                                 set->push_back(test_t);
 756                                 multipliers.erase(multipliers.begin());
 757                         }
 758
 759                 }
 760
 761                 /*
 762                  * Barrel distortion
 763                  */
 764
 765                 if (bda_mult != 0 && adj_b != 0) {
 766
 767                         for (unsigned int d = 0; d < get_offset().bd_count(); d++)
 768                         for (unsigned int s = 0; s < 2; s++) {
 769
 770                                 if (!multipliers.size())
 771                                         multipliers.push_back(1);
 772
 773                                 assert (multipliers.size());
 774                                 assert (finite(multipliers[0]));
 775
 776                                 ale_pos adj_s = (s ? -1 : 1) * adj_b * multipliers[0];
 777
 778                                 if (bda_rate > 0 && fabs(modified_bd[d] + adj_s - current_bd[d]) > bda_rate)
 779                                         continue;
 780
 781                                 diff_trans test_t = get_offset();
 782
 783                                 test_t.bd_modify(d, adj_s);
 784
 785                                 set->push_back(test_t);
 786                         }
 787                 }
 788         }
 789
 790         void confirm() {
 791                 assert(runs.size() == 2);
 792                 runs[0] = runs[1];
 793                 runs.pop_back();
 794         }
 795
 796         void discard() {
 797                 assert(runs.size() == 2);
 798                 runs.pop_back();
 799         }
 800
 801         void perturb_test(ale_pos perturb, ale_pos adj_p, ale_pos adj_o, ale_pos adj_b,
 802                         ale_pos *current_bd, ale_pos *modified_bd, int stable) {
 803
 804                 assert(runs.size() == 1);
 805
 806                 std::vector<diff_trans> t_set;
 807
 808                 if (perturb_multipliers.size() == 0) {
 809                         get_perturb_set(&t_set, adj_p, adj_o, adj_b,
 810                                         current_bd, modified_bd);
 811
 812                         for (unsigned int i = 0; i < t_set.size(); i++) {
 813                                 diff_stat_generic test = *this;
 814
 815                                 test.diff(si, t_set[i], ax_count, frame);
 816
 817                                 test.confirm();
 818
 819                                 if (finite(test.get_error_perturb()))
 820                                         perturb_multipliers.push_back(adj_p / test.get_error_perturb());
 821                                 else
 822                                         perturb_multipliers.push_back(1);
 823
 824                         }
 825
 826                         t_set.clear();
 827                 }
 828
 829                 get_perturb_set(&t_set, adj_p, adj_o, adj_b, current_bd, modified_bd,
 830                                 perturb_multipliers);
 831
 832                 int found_unreliable = 1;
 833                 std::vector<int> tested(t_set.size(), 0);
 834
 835                 for (unsigned int i = 0; i < t_set.size(); i++) {
 836                         run_index ori = get_run_index(i);
 837
 838                         /*
 839                          * Check for stability
 840                          */
 841                         if (stable
 842                          && old_runs.count(ori)
 843                          && old_runs[ori].offset == t_set[i])
 844                                 tested[i] = 1;
 845                 }
 846
 847                 std::vector<ale_pos> perturb_multipliers_original = perturb_multipliers;
 848
 849                 while (found_unreliable) {
 850
 851                         found_unreliable = 0;
 852
 853                         for (unsigned int i = 0; i < t_set.size(); i++) {
 854
 855                                 if (tested[i])
 856                                         continue;
 857
 858                                 diff(si, t_set[i], ax_count, frame);
 859
 860                                 if (!(i < perturb_multipliers.size())
 861                                  || !finite(perturb_multipliers[i])) {
 862
 863                                         perturb_multipliers.resize(i + 1);
 864
 865                                         if (finite(perturb_multipliers[i])
 866                                          && finite(adj_p)
 867                                          && finite(adj_p / runs[1].get_error_perturb())) {
 868                                                 perturb_multipliers[i] =
 869                                                         adj_p / runs[1].get_error_perturb();
 870
 871                                                 found_unreliable = 1;
 872                                         } else
 873                                                 perturb_multipliers[i] = 1;
 874
 875                                         continue;
 876                                 }
 877
 878                                 run_index ori = get_run_index(i);
 879
 880                                 if (old_runs.count(ori) == 0)
 881                                         old_runs.insert(std::pair<run_index, run>(ori, runs[1]));
 882                                 else
 883                                         old_runs[ori] = runs[1];
 884
 885                                 if (finite(perturb_multipliers_original[i])
 886                                  && finite(runs[1].get_error_perturb())
 887                                  && finite(adj_p)
 888                                  && finite(perturb_multipliers_original[i] * adj_p / runs[1].get_error_perturb()))
 889                                         perturb_multipliers[i] = perturb_multipliers_original[i]
 890                                                 * adj_p / runs[1].get_error_perturb();
 891                                 else
 892                                         perturb_multipliers[i] = 1;
 893
 894                                 tested[i] = 1;
 895
 896                                 if (better()
 897                                  && runs[1].get_error() < runs[0].get_error()
 898                                  && perturb_multipliers[i]
 899                                   / perturb_multipliers_original[i] < 2) {
 900                                         runs[0] = runs[1];
 901                                         runs.pop_back();
 902                                         return;
 903                                 }
 904
 905                         }
 906
 907                         if (runs.size() > 1)
 908                                 runs.pop_back();
 909
 910                         if (!found_unreliable)
 911                                 return;
 912                 }
 913         }
 914
 915 };
 916
 917 typedef diff_stat_generic<trans_single> diff_stat_t;
 918 typedef diff_stat_generic<trans_multi> diff_stat_multi;
 919
 920
 921 /*
 922  * Adjust exposure for an aligned frame B against reference A.
 923  *
 924  * Expects full-LOD images.
 925  *
 926  * Note: This method does not use any weighting, by certainty or
 927  * otherwise, in the first exposure registration pass, as any bias of
 928  * weighting according to color may also bias the exposure registration
 929  * result; it does use weighting, including weighting by certainty
 930  * (even if certainty weighting is not specified), in the second pass,
 931  * under the assumption that weighting by certainty improves handling
 932  * of out-of-range highlights, and that bias of exposure measurements
 933  * according to color may generally be less harmful after spatial
 934  * registration has been performed.
 935  */
 936 class exposure_ratio_iterate : public thread::decompose_domain {
 937         pixel_accum *asums;
 938         pixel_accum *bsums;
 939         pixel_accum *asum;
 940         pixel_accum *bsum;
 941         struct scale_cluster c;
 942         const transformation &t;
 943         int ax_count;
 944         int pass_number;
 945 protected:
 946         void prepare_subdomains(unsigned int N) {
 947                 asums = new pixel_accum[N];
 948                 bsums = new pixel_accum[N];
 949         }
 950         void subdomain_algorithm(unsigned int thread,
 951                         int i_min, int i_max, int j_min, int j_max) {
 952
 953                 point offset = c.accum->offset();
 954
 955                 for (mc_iterate m(i_min, i_max, j_min, j_max, thread); !m.done(); m++) {
 956
 957                         unsigned int i = (unsigned int) m.get_i();
 958                         unsigned int j = (unsigned int) m.get_j();
 959
 960                         if (ref_excluded(i, j, offset, c.ax_parameters, ax_count))
 961                                 continue;
 962
 963                         /*
 964                          * Transform
 965                          */
 966
 967                         struct point q;
 968
 969                         q = (c.input_scale < 1.0 && interpolant == NULL)
 970                           ? t.scaled_inverse_transform(
 971                                 point(i + offset[0], j + offset[1]))
 972                           : t.unscaled_inverse_transform(
 973                                 point(i + offset[0], j + offset[1]));
 974
 975                         /*
 976                          * Check that the transformed coordinates are within
 977                          * the boundaries of array c.input, that they are not
 978                          * subject to exclusion, and that the weight value in
 979                          * the accumulated array is nonzero.
 980                          */
 981
 982                         if (input_excluded(q[0], q[1], c.ax_parameters, ax_count))
 983                                 continue;
 984
 985                         if (q[0] >= 0
 986                          && q[0] <= c.input->height() - 1
 987                          && q[1] >= 0
 988                          && q[1] <= c.input->width() - 1
 989                          && ((pixel) c.certainty->get_pixel(i, j)).minabs_norm() != 0) {
 990                                 pixel a = c.accum->get_pixel(i, j);
 991                                 pixel b;
 992
 993                                 b = c.input->get_bl(q);
 994
 995                                 pixel weight = ((c.aweight && pass_number)
 996                                               ? (pixel) c.aweight->get_pixel(i, j)
 997                                               : pixel(1, 1, 1))
 998                                              * (pass_number
 999                                               ? psqrt(c.certainty->get_pixel(i, j)
1000                                                     * c.input_certainty->get_bl(q, 1))
1001                                               : pixel(1, 1, 1));
1002
1003                                 asums[thread] += a * weight;
1004                                 bsums[thread] += b * weight;
1005                         }
1006                 }
1007         }
1008
1009         void finish_subdomains(unsigned int N) {
1010                 for (unsigned int n = 0; n < N; n++) {
1011                         *asum += asums[n];
1012                         *bsum += bsums[n];
1013                 }
1014                 delete[] asums;
1015                 delete[] bsums;
1016         }
1017 public:
1018         exposure_ratio_iterate(pixel_accum *_asum,
1019                                pixel_accum *_bsum,
1020                                struct scale_cluster _c,
1021                                const transformation &_t,
1022                                int _ax_count,
1023                                int _pass_number) : decompose_domain(0, _c.accum->height(),
1024                                                                     0, _c.accum->width()),
1025                                                    t(_t) {
1026
1027                 asum = _asum;
1028                 bsum = _bsum;
1029                 c = _c;
1030                 ax_count = _ax_count;
1031                 pass_number = _pass_number;
1032         }
1033
1034         exposure_ratio_iterate(pixel_accum *_asum,
1035                                pixel_accum *_bsum,
1036                                struct scale_cluster _c,
1037                                const transformation &_t,
1038                                int _ax_count,
1039                                int _pass_number,
1040                                transformation::elem_bounds_int_t b) : decompose_domain(b.imin, b.imax,
1041                                                                        b.jmin, b.jmax),
1042                                                       t(_t) {
1043
1044                 asum = _asum;
1045                 bsum = _bsum;
1046                 c = _c;
1047                 ax_count = _ax_count;
1048                 pass_number = _pass_number;
1049         }
1050 };
1051
1052 static void set_exposure_ratio(unsigned int m, struct scale_cluster c,
1053                 const transformation &t, int ax_count, int pass_number) {
1054
1055         if (_exp_register == 2) {
1056                 /*
1057                  * Use metadata only.
1058                  */
1059                 ale_real gain_multiplier = image_rw::exp(m).get_gain_multiplier();
1060                 pixel multiplier = pixel(gain_multiplier, gain_multiplier, gain_multiplier);
1061
1062                 image_rw::exp(m).set_multiplier(multiplier);
1063                 ui::get()->exp_multiplier(multiplier[0],
1064                                           multiplier[1],
1065                                           multiplier[2]);
1066
1067                 return;
1068         }
1069
1070         pixel_accum asum(0, 0, 0), bsum(0, 0, 0);
1071
1072         exposure_ratio_iterate eri(&asum, &bsum, c, t, ax_count, pass_number);
1073         eri.run();
1074
1075         // std::cerr << (asum / bsum) << " ";
1076
1077         pixel_accum new_multiplier;
1078
1079         new_multiplier = asum / bsum * image_rw::exp(m).get_multiplier();
1080
1081         if (finite(new_multiplier[0])
1082          && finite(new_multiplier[1])
1083          && finite(new_multiplier[2])) {
1084                 image_rw::exp(m).set_multiplier(new_multiplier);
1085                 ui::get()->exp_multiplier(new_multiplier[0],
1086                                           new_multiplier[1],
1087                                           new_multiplier[2]);
1088         }
1089 }
1090
1091 static diff_stat_t _align_element(ale_pos perturb, ale_pos local_lower,
1092                 scale_cluster *scale_clusters, diff_stat_t here,
1093                 ale_pos adj_p, ale_pos adj_o, ale_pos adj_b,
1094                 ale_pos *current_bd, ale_pos *modified_bd,
1095                 astate_t *astate, int lod, scale_cluster si) {
1096
1097         /*
1098          * Run initial tests to get perturbation multipliers and error
1099          * centroids.
1100          */
1101
1102         std::vector<d2::trans_single> t_set;
1103
1104         here.get_perturb_set(&t_set, adj_p, adj_o, adj_b, current_bd, modified_bd);
1105
1106         int stable_count = 0;
1107
1108         while (perturb >= local_lower) {
1109
1110                 ui::get()->alignment_dims(scale_clusters[lod].accum->height(), scale_clusters[lod].accum->width(),
1111                                           scale_clusters[lod].input->height(), scale_clusters[lod].input->width());
1112
1113                 /*
1114                  * Orientational adjustment value in degrees.
1115                  *
1116                  * Since rotational perturbation is now specified as an
1117                  * arclength, we have to convert.
1118                  */
1119
1120                 ale_pos adj_o = 2 * (double) perturb
1121                                   / sqrt(pow(scale_clusters[0].input->height(), 2)
1122                                        + pow(scale_clusters[0].input->width(),  2))
1123                                   * 180
1124                                   / M_PI;
1125
1126                 /*
1127                  * Barrel distortion adjustment value
1128                  */
1129
1130                 ale_pos adj_b = perturb * bda_mult;
1131
1132                 trans_single old_offset = here.get_offset();
1133
1134                 here.perturb_test(perturb, adj_p, adj_o, adj_b, current_bd, modified_bd,
1135                         stable_count);
1136
1137                 if (here.get_offset() == old_offset)
1138                         stable_count++;
1139                 else
1140                         stable_count = 0;
1141
1142                 if (stable_count == 3) {
1143
1144                         stable_count = 0;
1145
1146                         perturb *= 0.5;
1147
1148                         if (lod > 0
1149                          && lod > lrint(log(perturb) / log(2)) - lod_preferred) {
1150
1151                                 /*
1152                                  * Work with images twice as large
1153                                  */
1154
1155                                 lod--;
1156                                 si = scale_clusters[lod];
1157
1158                                 /*
1159                                  * Rescale the transforms.
1160                                  */
1161
1162                                 ale_pos rescale_factor = (double) scale_factor
1163                                                        / (double) pow(2, lod)
1164                                                        / (double) here.get_offset().scale();
1165
1166                                 here.rescale(rescale_factor, si);
1167
1168                         } else {
1169                                 adj_p = perturb / pow(2, lod);
1170                         }
1171
1172                         /*
1173                          * Announce changes
1174                          */
1175
1176                         ui::get()->alignment_perturbation_level(perturb, lod);
1177                 }
1178
1179                 ui::get()->set_match(here.get_error());
1180                 ui::get()->set_offset(here.get_offset());
1181         }
1182
1183         if (lod > 0) {
1184                 ale_pos rescale_factor = (double) scale_factor
1185                                        / (double) here.get_offset().scale();
1186
1187                 here.rescale(rescale_factor, scale_clusters[0]);
1188         }
1189
1190         return here;
1191 }
1192
1193 /*
1194  * Align frame m against the reference.
1195  *
1196  * XXX: the transformation class currently combines ordinary
1197  * transformations with scaling.  This is somewhat convenient for
1198  * some things, but can also be confusing.  This method, _align(), is
1199  * one case where special care must be taken to ensure that the scale
1200  * is always set correctly (by using the 'rescale' method).
1201  */
1202 static diff_stat_multi _align(int m, int local_gs, astate_t *astate) {
1203
1204         ale_image input_frame = astate->get_input_frame();
1205
1206         /*
1207          * Local upper/lower data, possibly dependent on image
1208          * dimensions.
1209          */
1210
1211         ale_pos local_lower, local_upper;
1212         ale_accum local_gs_mo;
1213
1214         /*
1215          * Select the minimum dimension as the reference.
1216          */
1217
1218         ale_pos reference_size = input_frame->height();
1219         if (input_frame->width() < reference_size)
1220                 reference_size = input_frame->width();
1221         ale_accum reference_area = input_frame->height()
1222                                  * input_frame->width();
1223
1224         if (perturb_lower_percent)
1225                 local_lower = (double) perturb_lower
1226                             * (double) reference_size
1227                             * (double) 0.01
1228                             * (double) scale_factor;
1229         else
1230                 local_lower = perturb_lower;
1231
1232         if (perturb_upper_percent)
1233                 local_upper = (double) perturb_upper
1234                             * (double) reference_size
1235                             * (double) 0.01
1236                             * (double) scale_factor;
1237         else
1238                 local_upper = perturb_upper;
1239
1240         local_upper = pow(2, floor(log(local_upper) / log(2)));
1241
1242         if (gs_mo_percent)
1243                 local_gs_mo = (double) _gs_mo
1244                             * (double) reference_area
1245                             * (double) 0.01
1246                             * (double) scale_factor;
1247         else
1248                 local_gs_mo = _gs_mo;
1249
1250         /*
1251          * Logarithms aren't exact, so we divide repeatedly to discover
1252          * how many steps will occur, and pass this information to the
1253          * user interface.
1254          */
1255
1256         int step_count = 0;
1257         double step_variable = local_upper;
1258         while (step_variable >= local_lower) {
1259                 step_variable /= 2;
1260                 step_count++;
1261         }
1262
1263         ale_pos perturb = local_upper;
1264
1265         if (_keep) {
1266                 kept_t[latest] = latest_t;
1267                 kept_ok[latest] = latest_ok;
1268         }
1269
1270         /*
1271          * Determine how many levels of detail should be prepared, by
1272          * calculating the initial (largest) value for the
1273          * level-of-detail variable.
1274          */
1275
1276         int lod = lrint(log(perturb) / log(2)) - lod_preferred;
1277
1278         if (lod < 0)
1279                 lod = 0;
1280
1281         while (lod > 0 && (reference_image->width() < pow(2, lod) * min_dimension
1282                         || reference_image->height() < pow(2, lod) * min_dimension))
1283                 lod--;
1284
1285         unsigned int steps = (unsigned int) lod + 1;
1286
1287         /*
1288          * Prepare multiple levels of detail.
1289          */
1290
1291         int local_ax_count;
1292         struct scale_cluster *scale_clusters = init_clusters(m,
1293                         scale_factor, input_frame, steps,
1294                         &local_ax_count);
1295
1296 #error add check for non-NULL return
1297
1298         /*
1299          * Initialize the default initial transform
1300          */
1301
1302         astate->init_default();
1303
1304         /*
1305          * Set the default transformation.
1306          */
1307
1308         transformation offset = astate->get_default();
1309
1310         /*
1311          * Establish boundaries
1312          */
1313
1314         offset.set_current_bounds(reference_image);
1315
1316         ui::get()->alignment_degree_max(offset.get_coordinate(offset.stack_depth() - 1).degree);
1317
1318         if (offset.stack_depth() == 1) {
1319                 ui::get()->set_steps(step_count, 0);
1320         } else {
1321                 ui::get()->set_steps(offset.get_coordinate(offset.stack_depth() - 1).degree + 1, 1);
1322         }
1323
1324         /*
1325          * Load any file-specified transformations
1326          */
1327
1328         for (unsigned int index = 0; index < offset.stack_depth(); index++) {
1329                 int is_default = 1;
1330                 unsigned int index_2;
1331                 offset.set_current_index(index);
1332
1333                 offset = tload_next(tload, alignment_class == 2,
1334                                 offset,
1335                                 &is_default, offset.get_current_index() == 0);
1336
1337                 index_2 = offset.get_current_index();
1338
1339                 if (index_2 > index) {
1340                         for (unsigned int index_3 = index; index_3 < index_2; index_3++)
1341                                 astate->set_is_default(index_3, 1);
1342
1343                         index = index_2;
1344                 }
1345
1346                 astate->set_is_default(index, is_default);
1347         }
1348
1349         offset.set_current_index(0);
1350
1351         astate->init_frame_alignment_primary(&offset, lod, perturb);
1352
1353         /*
1354          * Control point alignment
1355          */
1356
1357         if (local_gs == 5) {
1358
1359                 transformation o = offset;
1360
1361                 /*
1362                  * Determine centroid data
1363                  */
1364
1365                 point current, previous;
1366                 centroids(m, &current, &previous);
1367
1368                 if (current.defined() && previous.defined()) {
1369                         o = orig_t;
1370                         o.set_dimensions(input_frame);
1371                         o.translate((previous - current) * o.scale());
1372                         current = previous;
1373                 }
1374
1375                 /*
1376                  * Determine rotation for alignment classes other than translation.
1377                  */
1378
1379                 ale_pos lowest_error = cp_rms_error(m, o);
1380
1381                 ale_pos rot_lower = 2 * (double) local_lower
1382                                   / sqrt(pow(scale_clusters[0].input->height(), 2)
1383                                        + pow(scale_clusters[0].input->width(),  2))
1384                                   * 180
1385                                   / M_PI;
1386
1387                 if  (alignment_class > 0)
1388                 for (double rot = 30; rot > rot_lower; rot /= 2)
1389                 for (double srot = -rot; srot < rot * 1.5; srot += rot * 2) {
1390                         int is_improved = 1;
1391                         while (is_improved) {
1392                                 is_improved = 0;
1393                                 transformation test_t = o;
1394                                 /*
1395                                  * XXX: is this right?
1396                                  */
1397                                 test_t.rotate(current * o.scale(), srot);
1398                                 ale_pos test_v = cp_rms_error(m, test_t);
1399
1400                                 if (test_v < lowest_error) {
1401                                         lowest_error = test_v;
1402                                         o = test_t;
1403                                         srot += 3 * rot;
1404                                         is_improved = 1;
1405                                 }
1406                         }
1407                 }
1408
1409                 /*
1410                  * Determine projective parameters through a local
1411                  * minimum search.
1412                  */
1413
1414                 if (alignment_class == 2) {
1415                         ale_pos adj_p = lowest_error;
1416
1417                         if (adj_p < local_lower)
1418                                 adj_p = local_lower;
1419
1420                         while (adj_p >= local_lower) {
1421                                 transformation test_t = o;
1422                                 int is_improved = 1;
1423                                 ale_pos test_v;
1424                                 ale_pos adj_s;
1425
1426                                 while (is_improved) {
1427                                         is_improved = 0;
1428
1429                                         for (int i = 0; i < 4; i++)
1430                                         for (int j = 0; j < 2; j++)
1431                                         for (adj_s = -adj_p; adj_s <= adj_p; adj_s += 2 * adj_p) {
1432
1433                                                 test_t = o;
1434
1435                                                 if (perturb_type == 0)
1436                                                         test_t.gpt_modify(j, i, adj_s);
1437                                                 else if (perturb_type == 1)
1438                                                         test_t.gr_modify(j, i, adj_s);
1439                                                 else
1440                                                         assert(0);
1441
1442                                                 test_v = cp_rms_error(m, test_t);
1443
1444                                                 if (test_v < lowest_error) {
1445                                                         lowest_error = test_v;
1446                                                         o = test_t;
1447                                                         adj_s += 3 * adj_p;
1448                                                         is_improved = 1;
1449                                                 }
1450                                         }
1451                                 }
1452                                 adj_p /= 2;
1453                         }
1454                 }
1455         }
1456
1457         /*
1458          * Pre-alignment exposure adjustment
1459          */
1460
1461         if (_exp_register) {
1462                 ui::get()->exposure_1();
1463                 set_exposure_ratio(m, scale_clusters[0], offset, local_ax_count, 0);
1464         }
1465
1466         /*
1467          * Scale transform for lod
1468          */
1469
1470         for (int lod_ = 0; lod_ < lod; lod_++) {
1471                 transformation s = offset;
1472                 transformation t = offset;
1473
1474                 t.rescale(1 / (double) 2);
1475
1476                 if (!(t.scaled_height() > 0 && t.scaled_height() < s.scaled_height())
1477                  || !(t.scaled_width() > 0  && t.scaled_width() < s.scaled_width())) {
1478                         perturb /= pow(2, lod - lod_);
1479                         lod = lod_;
1480                         break;
1481                 } else {
1482                         offset = t;
1483                 }
1484         }
1485
1486         ui::get()->set_offset(offset);
1487
1488         struct scale_cluster si = scale_clusters[lod];
1489
1490         /*
1491          * Projective adjustment value
1492          */
1493
1494         ale_pos adj_p = perturb / pow(2, lod);
1495
1496         /*
1497          * Orientational adjustment value in degrees.
1498          *
1499          * Since rotational perturbation is now specified as an
1500          * arclength, we have to convert.
1501          */
1502
1503         ale_pos adj_o = (double) 2 * (double) perturb
1504                           / sqrt(pow((double) scale_clusters[0].input->height(), (double) 2)
1505                                + pow((double) scale_clusters[0].input->width(),  (double) 2))
1506                           * (double) 180
1507                           / M_PI;
1508
1509         /*
1510          * Barrel distortion adjustment value
1511          */
1512
1513         ale_pos adj_b = perturb * bda_mult;
1514
1515         /*
1516          * Global search overlap requirements.
1517          */
1518
1519         local_gs_mo = (double) local_gs_mo / pow(pow(2, lod), 2);
1520
1521         /*
1522          * Alignment statistics.
1523          */
1524
1525         diff_stat_t here(offset.elem_bounds());
1526
1527         /*
1528          * Current difference (error) value
1529          */
1530
1531         ui::get()->prematching();
1532         here.diff(si, offset.get_current_element(), local_ax_count, m);
1533         ui::get()->set_match(here.get_error());
1534
1535         /*
1536          * Current and modified barrel distortion parameters
1537          */
1538
1539         ale_pos current_bd[BARREL_DEGREE];
1540         ale_pos modified_bd[BARREL_DEGREE];
1541         offset.bd_get(current_bd);
1542         offset.bd_get(modified_bd);
1543
1544         /*
1545          * Translational global search step
1546          */
1547
1548         if (perturb >= local_lower && local_gs != 0 && local_gs != 5
1549          && (local_gs != 6 || astate->get_is_default(0))) {
1550
1551                 ui::get()->global_alignment(perturb, lod);
1552                 ui::get()->gs_mo(local_gs_mo);
1553
1554                 test_globals(&here, si, offset, local_gs, adj_p,
1555                                 local_ax_count, m, local_gs_mo, perturb);
1556
1557                 ui::get()->set_match(here.get_error());
1558                 ui::get()->set_offset(here.get_offset());
1559         }
1560
1561         /*
1562          * Perturbation adjustment loop.
1563          */
1564
1565         offset.set_current_element(here.get_offset());
1566
1567         for (unsigned int i = 0; i < offset.stack_depth(); i++) {
1568
1569                 ui::get()->aligning_element(i, offset.stack_depth());
1570
1571                 offset.set_current_index(i);
1572
1573                 ui::get()->start_multi_alignment_element(offset);
1574
1575                 ui::get()->set_offset(offset);
1576
1577                 if (i > 0) {
1578                         astate->init_frame_alignment_nonprimary(&offset, lod, perturb, i);
1579
1580                         if (_exp_register == 1) {
1581                                 ui::get()->exposure_1();
1582                                 pixel_accum asum(0, 0, 0), bsum(0, 0, 0);
1583                                 exposure_ratio_iterate eri(&asum, &bsum, scale_clusters[0], offset, local_ax_count, 0,
1584                                         offset.elem_bounds().scale_to_bounds(scale_clusters[0].accum->height(),
1585                                                                              scale_clusters[0].accum->width()));
1586
1587                                 eri.run();
1588                                 pixel_accum tr = asum / bsum;
1589                                 ui::get()->exp_multiplier(tr[0], tr[1], tr[2]);
1590                                 offset.set_tonal_multiplier(tr);
1591                         }
1592                 }
1593
1594                 int e_lod = lod;
1595                 int e_div = offset.get_current_coordinate().degree;
1596                 ale_pos e_perturb = perturb;
1597                 ale_pos e_adj_p = adj_p;
1598                 ale_pos e_adj_b = adj_b;
1599
1600                 for (int d = 0; d < e_div; d++) {
1601                         e_adj_b = 0;
1602                         e_perturb *= 0.5;
1603                         if (e_lod > 0) {
1604                                 e_lod--;
1605                         } else {
1606                                 e_adj_p *= 0.5;
1607                         }
1608                 }
1609
1610                 if (i > 0) {
1611
1612                         d2::trans_multi::elem_bounds_t b = offset.elem_bounds();
1613
1614                         for (int dim_satisfied = 0; e_lod > 0 && !dim_satisfied; ) {
1615                                 int height = scale_clusters[e_lod].accum->height();
1616                                 int width = scale_clusters[e_lod].accum->width();
1617
1618                                 d2::trans_multi::elem_bounds_int_t bi = b.scale_to_bounds(height, width);
1619
1620                                 dim_satisfied = bi.satisfies_min_dim(min_dimension);
1621
1622                                 if (!dim_satisfied) {
1623                                         e_lod--;
1624                                         e_adj_p *= 2;
1625                                 }
1626                         }
1627
1628                         /*
1629                          * Scale transform for lod
1630                          */
1631
1632                         for (int lod_ = 0; lod_ < e_lod; lod_++) {
1633                                 trans_single s = offset.get_element(i);
1634                                 trans_single t = offset.get_element(i);
1635
1636                                 t.rescale(1 / (double) 2);
1637
1638                                 if (!(t.scaled_height() > 0 && t.scaled_height() < s.scaled_height())
1639                                  || !(t.scaled_width() > 0  && t.scaled_width() < s.scaled_width())) {
1640                                         e_perturb /= pow(2, e_lod - lod_);
1641                                         e_lod = lod_;
1642                                         break;
1643                                 } else {
1644                                         offset.set_element(i, t);
1645                                 }
1646                         }
1647
1648                         ui::get()->set_offset(offset);
1649                 }
1650
1651                 /*
1652                  * Announce perturbation size
1653                  */
1654
1655                 ui::get()->aligning(e_perturb, e_lod);
1656
1657                 si = scale_clusters[e_lod];
1658
1659                 here.set_elem_bounds(offset.elem_bounds());
1660
1661                 here.diff(si, offset.get_current_element(), local_ax_count, m);
1662
1663                 here.confirm();
1664
1665                 here = check_ancestor_path(offset, si, here, local_ax_count, m);
1666
1667                 here = _align_element(e_perturb, local_lower, scale_clusters,
1668                                 here, e_adj_p, adj_o, e_adj_b, current_bd, modified_bd,
1669                                 astate, e_lod, si);
1670
1671                 offset.rescale(here.get_offset().scale() / offset.scale());
1672
1673                 offset.set_current_element(here.get_offset());
1674
1675                 if (i > 0 && _exp_register == 1) {
1676                         if (ma_cert_satisfied(scale_clusters[0], offset, i)) {
1677                                 ui::get()->exposure_2();
1678                                 pixel_accum asum(0, 0, 0), bsum(0, 0, 0);
1679                                 exposure_ratio_iterate eri(&asum, &bsum, scale_clusters[0], offset, local_ax_count, 1,
1680                                         offset.elem_bounds().scale_to_bounds(scale_clusters[0].accum->height(),
1681                                                                              scale_clusters[0].accum->width()));
1682
1683                                 eri.run();
1684                                 pixel_accum tr = asum / bsum;
1685                                 ui::get()->exp_multiplier(tr[0], tr[1], tr[2]);
1686                                 offset.set_tonal_multiplier(tr);
1687                         } else {
1688                                 offset.set_tonal_multiplier(offset.get_element(offset.parent_index(i)).get_tonal_multiplier(point(0, 0)));
1689                         }
1690                 } else if (_exp_register == 1) {
1691                         ui::get()->exposure_2();
1692                         set_exposure_ratio(m, scale_clusters[0], offset, local_ax_count, 1);
1693                 }
1694
1695                 ui::get()->set_offset(offset);
1696
1697                 if (i + 1 == offset.stack_depth())
1698                         ui::get()->alignment_degree_complete(offset.get_coordinate(i).degree);
1699                 else if (offset.get_coordinate(i).degree != offset.get_coordinate(i + 1).degree)
1700                         ui::get()->alignment_degree_complete(offset.get_coordinate(i + 1).degree);
1701         }
1702
1703         offset.set_current_index(0);
1704
1705         ui::get()->multi();
1706         offset.set_multi(reference_image, input_frame);
1707
1708         /*
1709          * Recalculate error on whole frame.
1710          */
1711
1712         ui::get()->postmatching();
1713         diff_stat_generic<transformation> multi_here(offset.elem_bounds());
1714         multi_here.diff(scale_clusters[0], offset, local_ax_count, m);
1715         ui::get()->set_match(multi_here.get_error());
1716
1717         /*
1718          * Free the level-of-detail structures
1719          */
1720
1721         final_clusters(scale_clusters, scale_factor, steps);
1722
1723         /*
1724          * Ensure that the match meets the threshold.
1725          */
1726
1727         if (threshold_ok(multi_here.get_error())) {
1728                 /*
1729                  * Update alignment variables
1730                  */
1731                 latest_ok = 1;
1732                 astate->set_default(offset);
1733                 astate->set_final(offset);
1734                 ui::get()->alignment_match_ok();
1735         } else if (local_gs == 4) {
1736
1737                 /*
1738                  * Align with outer starting points.
1739                  */
1740
1741                 /*
1742                  * XXX: This probably isn't exactly the right thing to do,
1743                  * since variables like old_initial_value have been overwritten.
1744                  */
1745
1746                 diff_stat_multi nested_result = _align(m, -1, astate);
1747
1748                 if (threshold_ok(nested_result.get_error())) {
1749                         return nested_result;
1750                 } else if (nested_result.get_error() < multi_here.get_error()) {
1751                         multi_here = nested_result;
1752                 }
1753
1754                 if (is_fail_default)
1755                         offset = astate->get_default();
1756
1757                 ui::get()->set_match(multi_here.get_error());
1758                 ui::get()->alignment_no_match();
1759
1760         } else if (local_gs == -1) {
1761
1762                 latest_ok = 0;
1763                 latest_t = offset;
1764                 return multi_here;
1765
1766         } else {
1767                 if (is_fail_default)
1768                         offset = astate->get_default();
1769                 latest_ok = 0;
1770                 ui::get()->alignment_no_match();
1771         }
1772
1773         /*
1774          * Write the tonal registration multiplier as a comment.
1775          */
1776
1777         pixel trm = image_rw::exp(m).get_multiplier();
1778         tsave_trm(tsave, trm[0], trm[1], trm[2]);
1779
1780         /*
1781          * Save the transformation information
1782          */
1783
1784         for (unsigned int index = 0; index < offset.stack_depth(); index++) {
1785                 offset.set_current_index(index);
1786
1787                 tsave_next(tsave, offset, alignment_class == 2,
1788                                   offset.get_current_index() == 0);
1789         }
1790
1791         offset.set_current_index(0);
1792
1793         /*
1794          * Update match statistics.
1795          */
1796
1797         match_sum += (1 - multi_here.get_error()) * (ale_accum) 100;
1798         match_count++;
1799         latest = m;
1800         latest_t = offset;
1801
1802         return multi_here;
1803 }
1804
1805 int ale_align(ale_context ac, ale_image a, ale_image b, ale_trans start,
1806                 ale_align_properties align_properties, ale_trans result) {
1807 #warning function unfinished.
1808 }
1809