d3/scene.h

   1 // Copyright 2003, 2004, 2005, 2006 David Hilvert <dhilvert@auricle.dyndns.org>,
   2 //                                                <dhilvert@ugcs.caltech.edu>
   3
   4 /*  This file is part of the Anti-Lamenessing Engine.
   5
   6     The Anti-Lamenessing Engine is free software; you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation; either version 2 of the License, or
   9     (at your option) any later version.
  10
  11     The Anti-Lamenessing Engine is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with the Anti-Lamenessing Engine; if not, write to the Free Software
  18     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  19 */
  20
  21 /*
  22  * d3/scene.h: Representation of a 3D scene.
  23  */
  24
  25 #ifndef __scene_h__
  26 #define __scene_h__
  27
  28 #include "point.h"
  29
  30 /*
  31  * View angle multiplier.
  32  *
  33  * Setting this to a value larger than one can be useful for debugging.
  34  */
  35
  36 #define VIEW_ANGLE_MULTIPLIER 1
  37
  38 class scene {
  39
  40         /*
  41          * Clipping planes
  42          */
  43         static ale_pos front_clip;
  44         static ale_pos rear_clip;
  45
  46         /*
  47          * Decimation exponents for geometry
  48          */
  49         static int primary_decimation_upper;
  50         static int input_decimation_lower;
  51         static int output_decimation_preferred;
  52
  53         /*
  54          * Output clipping
  55          */
  56         static int output_clip;
  57
  58         /*
  59          * Model files
  60          */
  61         static const char *load_model_name;
  62         static const char *save_model_name;
  63
  64         /*
  65          * Occupancy attenuation
  66          */
  67
  68         static double occ_att;
  69
  70         /*
  71          * Normalization of output by weight
  72          */
  73
  74         static int normalize_weights;
  75
  76         /*
  77          * Filtering data
  78          */
  79
  80         static int use_filter;
  81         static const char *d3chain_type;
  82
  83         /*
  84          * Falloff exponent
  85          */
  86
  87         static double falloff_exponent;
  88
  89         /*
  90          * Third-camera error multiplier
  91          */
  92         static double tc_multiplier;
  93
  94         /*
  95          * Occupancy update iterations
  96          */
  97         static unsigned int ou_iterations;
  98
  99         /*
 100          * Pairwise ambiguity
 101          */
 102         static unsigned int pairwise_ambiguity;
 103
 104         /*
 105          * Pairwise comparisons
 106          */
 107         static const char *pairwise_comparisons;
 108
 109         /*
 110          * 3D Post-exclusion
 111          */
 112         static int d3px_count;
 113         static double *d3px_parameters;
 114
 115         /*
 116          * Nearness threshold
 117          */
 118         static const ale_real nearness;
 119
 120         /*
 121          * Encounter threshold for defined pixels.
 122          */
 123         static double encounter_threshold;
 124
 125         /*
 126          * Median calculation radii.
 127          */
 128         static double depth_median_radius;
 129         static double diff_median_radius;
 130
 131         /*
 132          * Flag for subspace traversal.
 133          */
 134         static int subspace_traverse;
 135
 136         /*
 137          * Structure to hold input frame information at levels of
 138          * detail between full detail and full decimation.
 139          */
 140         class lod_image {
 141                 unsigned int f;
 142                 unsigned int entries;
 143                 std::vector<const d2::image *> im;
 144                 std::vector<pt> transformation;
 145
 146         public:
 147                 /*
 148                  * Constructor
 149                  */
 150                 lod_image(unsigned int _f) {
 151
 152                         pt _pt;
 153
 154                         f = _f;
 155                         im.push_back(d2::image_rw::copy(f, "3D reference image"));
 156                         assert(im.back());
 157                         entries = 1;
 158                         _pt = d3::align::projective(f);
 159                         _pt.scale(1 / _pt.scale_2d());
 160                         transformation.push_back(_pt);
 161
 162                         while (im.back()->height() > 4
 163                             && im.back()->width() > 4) {
 164
 165                                 im.push_back(im.back()->scale_by_half("3D, reduced LOD"));
 166                                 assert(im.back());
 167
 168                                 _pt.scale(1 / _pt.scale_2d() / pow(2, entries));
 169                                 transformation.push_back(_pt);
 170
 171                                 entries += 1;
 172                         }
 173                 }
 174
 175                 /*
 176                  * Get the number of scales
 177                  */
 178                 unsigned int count() {
 179                         return entries;
 180                 }
 181
 182                 /*
 183                  * Get an image.
 184                  */
 185                 const d2::image *get_image(unsigned int i) {
 186                         assert(i >= 0);
 187                         assert(i < entries);
 188                         return im[i];
 189                 }
 190
 191                 int in_bounds(d2::point p) {
 192                         return im[0]->in_bounds(p);
 193                 }
 194
 195                 /*
 196                  * Get a 'trilinear' color.  We currently don't do interpolation
 197                  * between levels of detail; hence, it's discontinuous in tl_coord.
 198                  */
 199                 d2::pixel get_tl(d2::point p, ale_pos tl_coord) {
 200
 201                         assert(in_bounds(p));
 202
 203                         tl_coord = round(tl_coord);
 204
 205                         if (tl_coord >= entries)
 206                                 tl_coord = entries;
 207                         if (tl_coord < 0)
 208                                 tl_coord = 0;
 209
 210                         p = p / (ale_pos) pow(2, tl_coord);
 211
 212                         unsigned int itlc = (unsigned int) tl_coord;
 213
 214                         if (p[0] > im[itlc]->height() - 1)
 215                                 p[0] = im[itlc]->height() - 1;
 216                         if (p[1] > im[itlc]->width() - 1)
 217                                 p[1] = im[itlc]->width() - 1;
 218
 219                         assert(p[0] >= 0);
 220                         assert(p[1] >= 0);
 221
 222                         return im[itlc]->get_bl(p);
 223                 }
 224
 225                 d2::pixel get_max_diff(d2::point p, ale_pos tl_coord) {
 226                         assert(in_bounds(p));
 227
 228                         tl_coord = round(tl_coord);
 229
 230                         if (tl_coord >= entries)
 231                                 tl_coord = entries;
 232                         if (tl_coord < 0)
 233                                 tl_coord = 0;
 234
 235                         p = p / (ale_pos) pow(2, tl_coord);
 236
 237                         unsigned int itlc = (unsigned int) tl_coord;
 238
 239                         if (p[0] > im[itlc]->height() - 1)
 240                                 p[0] = im[itlc]->height() - 1;
 241                         if (p[1] > im[itlc]->width() - 1)
 242                                 p[1] = im[itlc]->width() - 1;
 243
 244                         assert(p[0] >= 0);
 245                         assert(p[1] >= 0);
 246
 247                         return im[itlc]->get_max_diff(p);
 248                 }
 249
 250                 /*
 251                  * Get the transformation
 252                  */
 253                 pt get_t(unsigned int i) {
 254                         assert(i >= 0);
 255                         assert(i < entries);
 256                         return transformation[i];
 257                 }
 258
 259                 /*
 260                  * Get the camera origin in world coordinates
 261                  */
 262                 point origin() {
 263                         return transformation[0].origin();
 264                 }
 265
 266                 /*
 267                  * Destructor
 268                  */
 269                 ~lod_image() {
 270                         for (unsigned int i = 0; i < entries; i++) {
 271                                 delete im[i];
 272                         }
 273                 }
 274         };
 275
 276         /*
 277          * Structure to hold weight information for reference images.
 278          */
 279         class ref_weights {
 280                 unsigned int f;
 281                 std::vector<d2::image *> weights;
 282                 pt transformation;
 283                 int tc_low;
 284                 int tc_high;
 285                 int initialized;
 286
 287                 void set_image(d2::image *im, ale_real value) {
 288                         assert(im);
 289                         for (unsigned int i = 0; i < im->height(); i++)
 290                         for (unsigned int j = 0; j < im->width(); j++)
 291                                 im->pix(i, j) = d2::pixel(value, value, value);
 292                 }
 293
 294                 d2::image *make_image(ale_pos sf, ale_real init_value = 0) {
 295                         d2::image *result = new d2::image_ale_real(
 296                                         (unsigned int) ceil(transformation.unscaled_height() * sf),
 297                                         (unsigned int) ceil(transformation.unscaled_width() * sf), 3);
 298                         assert(result);
 299
 300                         if (init_value != 0)
 301                                 set_image(result, init_value);
 302
 303                         return result;
 304                 }
 305
 306         public:
 307
 308                 /*
 309                  * Explicit weight subtree
 310                  */
 311                 struct subtree {
 312                         ale_real node_value;
 313                         subtree *children[2][2];
 314
 315                         subtree(ale_real nv, subtree *a, subtree *b, subtree *c, subtree *d) {
 316                                 node_value = nv;
 317                                 children[0][0] = a;
 318                                 children[0][1] = b;
 319                                 children[1][0] = c;
 320                                 children[1][1] = d;
 321                         }
 322
 323                         ~subtree() {
 324                                 for (int i = 0; i < 2; i++)
 325                                 for (int j = 0; j < 2; j++)
 326                                         delete children[i][j];
 327                         }
 328                 };
 329
 330                 /*
 331                  * Constructor
 332                  */
 333                 ref_weights(unsigned int _f) {
 334                         f = _f;
 335                         transformation = d3::align::projective(f);
 336                         initialized = 0;
 337                 }
 338
 339                 /*
 340                  * Check spatial bounds.
 341                  */
 342                 int in_spatial_bounds(point p) {
 343
 344                         if (!p.defined())
 345                                 return 0;
 346
 347                         if (p[0] < 0)
 348                                 return 0;
 349                         if (p[1] < 0)
 350                                 return 0;
 351                         if (p[0] > transformation.unscaled_height() - 1)
 352                                 return 0;
 353                         if (p[1] > transformation.unscaled_width() - 1)
 354                                 return 0;
 355                         if (p[2] >= 0)
 356                                 return 0;
 357
 358                         return 1;
 359                 }
 360
 361                 int in_spatial_bounds(const space::traverse &t) {
 362                         point p = transformation.centroid(t);
 363                         return in_spatial_bounds(p);
 364                 }
 365
 366                 /*
 367                  * Increase resolution to the given level.
 368                  */
 369                 void increase_resolution(int tc, unsigned int i, unsigned int j) {
 370                         d2::image *im = weights[tc - tc_low];
 371                         assert(im);
 372                         assert(i <= im->height() - 1);
 373                         assert(j <= im->width() - 1);
 374
 375                         /*
 376                          * Check for the cases known to have no lower level of detail.
 377                          */
 378
 379                         if (im->pix(i, j)[0] == -1)
 380                                 return;
 381
 382                         if (tc == tc_high)
 383                                 return;
 384
 385                         increase_resolution(tc + 1, i / 2, j / 2);
 386
 387                         /*
 388                          * Load the lower-level image structure.
 389                          */
 390
 391                         d2::image *iim = weights[tc + 1 - tc_low];
 392
 393                         assert(iim);
 394                         assert(i / 2 <= iim->height() - 1);
 395                         assert(j / 2 <= iim->width() - 1);
 396
 397                         /*
 398                          * Check for the case where no lower level of detail is
 399                          * available.
 400                          */
 401
 402                         if (iim->pix(i / 2, j / 2)[0] == -1)
 403                                 return;
 404
 405                         /*
 406                          * Spread out the lower level of detail among (uninitialized)
 407                          * peer values.
 408                          */
 409
 410                         for (unsigned int ii = (i / 2) * 2; ii < (i / 2) * 2 + 1; ii++)
 411                         for (unsigned int jj = (j / 2) * 2; jj < (j / 2) * 2 + 1; jj++) {
 412                                 assert(ii <= im->height() - 1);
 413                                 assert(jj <= im->width() - 1);
 414                                 assert(im->pix(ii, jj)[0] == 0);
 415
 416                                 im->pix(ii, jj) = iim->pix(i / 2, j / 2);
 417                         }
 418
 419                         /*
 420                          * Set the lower level of detail to point here.
 421                          */
 422
 423                         iim->pix(i / 2, j / 2)[0] = -1;
 424                 }
 425
 426                 /*
 427                  * Add weights to positive higher-resolution pixels where
 428                  * found when their current values match the given subtree
 429                  * values; set negative pixels to zero and return 0 if no
 430                  * positive higher-resolution pixels are found.
 431                  */
 432                 int add_partial(int tc, unsigned int i, unsigned int j, ale_real weight, subtree *st) {
 433                         d2::image *im = weights[tc - tc_low];
 434                         assert(im);
 435
 436                         if (i == im->height() - 1
 437                          || j == im->width() - 1) {
 438                                 return 1;
 439                         }
 440
 441                         assert(i <= im->height() - 1);
 442                         assert(j <= im->width() - 1);
 443
 444                         /*
 445                          * Check for positive values.
 446                          */
 447
 448                         if (im->pix(i, j)[0] > 0) {
 449                                 if (st && st->node_value == im->pix(i, j)[0])
 450                                         im->pix(i, j)[0] += weight * (1 - im->pix(i, j)[0]);
 451                                 return 1;
 452                         }
 453
 454                         /*
 455                          * Handle the case where there are no higher levels of detail.
 456                          */
 457
 458                         if (tc == tc_low) {
 459                                 if (im->pix(i, j)[0] != 0) {
 460                                         fprintf(stderr, "failing assertion: im[%d]->pix(%d, %d)[0] == %g\n", tc, i, j,
 461                                                         im->pix(i, j)[0]);
 462                                 }
 463                                 assert(im->pix(i, j)[0] == 0);
 464                                 return 0;
 465                         }
 466
 467                         /*
 468                          * Handle the case where higher levels of detail are available.
 469                          */
 470
 471                         int success[2][2];
 472
 473                         for (int ii = 0; ii < 2; ii++)
 474                         for (int jj = 0; jj < 2; jj++)
 475                                 success[ii][jj] = add_partial(tc - 1, i * 2 + ii, j * 2 + jj, weight,
 476                                                 st ? st->children[ii][jj] : NULL);
 477
 478                         if (!success[0][0]
 479                          && !success[0][1]
 480                          && !success[1][0]
 481                          && !success[1][1]) {
 482                                 im->pix(i, j)[0] = 0;
 483                                 return 0;
 484                         }
 485
 486                         d2::image *iim = weights[tc - 1 - tc_low];
 487                         assert(iim);
 488
 489                         for (int ii = 0; ii < 2; ii++)
 490                         for (int jj = 0; jj < 2; jj++)
 491                                 if (success[ii][jj] == 0) {
 492                                         assert(i * 2 + ii < iim->height());
 493                                         assert(j * 2 + jj < iim->width());
 494
 495                                         iim->pix(i * 2 + ii, j * 2 + jj)[0] = weight;
 496                                 }
 497
 498                         im->pix(i, j)[0] = -1;
 499
 500                         return 1;
 501                 }
 502
 503                 /*
 504                  * Add weight.
 505                  */
 506                 void add_weight(int tc, unsigned int i, unsigned int j, ale_real weight, subtree *st) {
 507
 508                         assert (weight >= 0);
 509
 510                         d2::image *im = weights[tc - tc_low];
 511                         assert(im);
 512
 513 //                      fprintf(stderr, "[aw tc=%d i=%d j=%d imax=%d jmax=%d]\n",
 514 //                                      tc, i, j, im->height(), im->width());
 515
 516                         assert(i <= im->height() - 1);
 517                         assert(j <= im->width() - 1);
 518
 519                         /*
 520                          * Increase resolution, if necessary
 521                          */
 522
 523                         increase_resolution(tc, i, j);
 524
 525                         /*
 526                          * Attempt to add the weight at levels of detail
 527                          * where weight is defined.
 528                          */
 529
 530                         if (add_partial(tc, i, j, weight, st))
 531                                 return;
 532
 533                         /*
 534                          * If no weights are defined at any level of detail,
 535                          * then set the weight here.
 536                          */
 537
 538                         im->pix(i, j)[0] = weight;
 539                 }
 540
 541                 void add_weight(int tc, d2::point p, ale_real weight, subtree *st) {
 542
 543                         assert (weight >= 0);
 544
 545                         p *= pow(2, -tc);
 546
 547                         unsigned int i = (unsigned int) floor(p[0]);
 548                         unsigned int j = (unsigned int) floor(p[1]);
 549
 550                         add_weight(tc, i, j, weight, st);
 551                 }
 552
 553                 void add_weight(const space::traverse &t, ale_real weight, subtree *st) {
 554
 555                         assert (weight >= 0);
 556
 557                         if (weight == 0)
 558                                 return;
 559
 560                         ale_pos tc = transformation.trilinear_coordinate(t);
 561                         point p = transformation.centroid(t);
 562                         assert(in_spatial_bounds(p));
 563
 564                         tc = round(tc);
 565
 566                         /*
 567                          * Establish a reasonable (?) upper bound on resolution.
 568                          */
 569
 570                         if (tc < input_decimation_lower) {
 571                                 weight /= pow(4, (input_decimation_lower - tc));
 572                                 tc = input_decimation_lower;
 573                         }
 574
 575                         /*
 576                          * Initialize, if necessary.
 577                          */
 578
 579                         if (!initialized) {
 580                                 tc_low = tc_high = (int) tc;
 581
 582                                 ale_real sf = pow(2, -tc);
 583
 584                                 weights.push_back(make_image(sf));
 585
 586                                 initialized = 1;
 587                         }
 588
 589                         /*
 590                          * Check resolution bounds
 591                          */
 592
 593                         assert (tc_low <= tc_high);
 594
 595                         /*
 596                          * Generate additional levels of detail, if necessary.
 597                          */
 598
 599                         while (tc < tc_low) {
 600                                 tc_low--;
 601
 602                                 ale_real sf = pow(2, -tc_low);
 603
 604                                 weights.insert(weights.begin(), make_image(sf));
 605                         }
 606
 607                         while (tc > tc_high) {
 608                                 tc_high++;
 609
 610                                 ale_real sf = pow(2, -tc_high);
 611
 612                                 weights.push_back(make_image(sf, -1));
 613                         }
 614
 615                         add_weight((int) tc, p.xy(), weight, st);
 616                 }
 617
 618                 /*
 619                  * Get weight
 620                  */
 621                 ale_real get_weight(int tc, unsigned int i, unsigned int j) {
 622
 623 //                      fprintf(stderr, "[gw tc=%d i=%u j=%u tclow=%d tchigh=%d]\n",
 624 //                                      tc, i, j, tc_low, tc_high);
 625
 626                         if (tc < tc_low || !initialized)
 627                                 return 0;
 628
 629                         if (tc > tc_high) {
 630                                 return (get_weight(tc - 1, i * 2 + 0, j * 2 + 0)
 631                                       + get_weight(tc - 1, i * 2 + 1, j * 2 + 0)
 632                                       + get_weight(tc - 1, i * 2 + 1, j * 2 + 1)
 633                                       + get_weight(tc - 1, i * 2 + 0, j * 2 + 1)) / 4;
 634                         }
 635
 636                         assert(weights.size() > (unsigned int) (tc - tc_low));
 637
 638                         d2::image *im = weights[tc - tc_low];
 639                         assert(im);
 640
 641                         if (i == im->height())
 642                                 return 1;
 643                         if (j == im->width())
 644                                 return 1;
 645
 646                         assert(i < im->height());
 647                         assert(j < im->width());
 648
 649                         if (im->pix(i, j)[0] == -1) {
 650                                 return (get_weight(tc - 1, i * 2 + 0, j * 2 + 0)
 651                                       + get_weight(tc - 1, i * 2 + 1, j * 2 + 0)
 652                                       + get_weight(tc - 1, i * 2 + 1, j * 2 + 1)
 653                                       + get_weight(tc - 1, i * 2 + 0, j * 2 + 1)) / 4;
 654                         }
 655
 656                         if (im->pix(i, j)[0] == 0) {
 657                                 if (tc == tc_high)
 658                                         return 0;
 659                                 if (weights[tc - tc_low + 1]->pix(i / 2, j / 2)[0] == -1)
 660                                         return 0;
 661                                 return get_weight(tc + 1, i / 2, j / 2);
 662                         }
 663
 664                         return im->pix(i, j)[0];
 665                 }
 666
 667                 ale_real get_weight(int tc, d2::point p) {
 668
 669                         p *= pow(2, -tc);
 670
 671                         unsigned int i = (unsigned int) floor(p[0]);
 672                         unsigned int j = (unsigned int) floor(p[1]);
 673
 674                         return get_weight(tc, i, j);
 675                 }
 676
 677                 ale_real get_weight(const space::traverse &t) {
 678                         ale_pos tc = transformation.trilinear_coordinate(t);
 679                         point p = transformation.centroid(t);
 680                         assert(in_spatial_bounds(p));
 681
 682                         if (!initialized)
 683                                 return 0;
 684
 685                         tc = round(tc);
 686
 687                         if (tc < tc_low) {
 688                                 tc = tc_low;
 689                         }
 690
 691                         return get_weight((int) tc, p.xy());
 692                 }
 693
 694                 /*
 695                  * Check whether a subtree is simple.
 696                  */
 697                 int is_simple(subtree *s) {
 698                         assert (s);
 699
 700                         if (s->node_value == -1
 701                          && s->children[0][0] == NULL
 702                          && s->children[0][1] == NULL
 703                          && s->children[1][0] == NULL
 704                          && s->children[1][1] == NULL)
 705                                 return 1;
 706
 707                         return 0;
 708                 }
 709
 710                 /*
 711                  * Get a weight subtree.
 712                  */
 713                 subtree *get_subtree(int tc, unsigned int i, unsigned int j) {
 714
 715                         /*
 716                          * tc > tc_high is handled recursively.
 717                          */
 718
 719                         if (tc > tc_high) {
 720                                 subtree *result = new subtree(-1,
 721                                                 get_subtree(tc - 1, i * 2 + 0, j * 2 + 0),
 722                                                 get_subtree(tc - 1, i * 2 + 0, j * 2 + 1),
 723                                                 get_subtree(tc - 1, i * 2 + 1, j * 2 + 0),
 724                                                 get_subtree(tc - 1, i * 2 + 1, j * 2 + 1));
 725
 726                                 if (is_simple(result)) {
 727                                         delete result;
 728                                         return NULL;
 729                                 }
 730
 731                                 return result;
 732                         }
 733
 734                         assert(tc >= tc_low);
 735                         assert(weights[tc - tc_low]);
 736
 737                         d2::image *im = weights[tc - tc_low];
 738
 739                         /*
 740                          * Rectangular images will, in general, have
 741                          * out-of-bounds tree sections.  Handle this case.
 742                          */
 743
 744                         if (i >= im->height())
 745                                 return NULL;
 746                         if (j >= im->width())
 747                                 return NULL;
 748
 749                         /*
 750                          * -1 weights are handled recursively
 751                          */
 752
 753                         if (im->pix(i, j)[0] == -1) {
 754                                 subtree *result = new subtree(-1,
 755                                                 get_subtree(tc - 1, i * 2 + 0, j * 2 + 0),
 756                                                 get_subtree(tc - 1, i * 2 + 0, j * 2 + 1),
 757                                                 get_subtree(tc - 1, i * 2 + 1, j * 2 + 0),
 758                                                 get_subtree(tc - 1, i * 2 + 1, j * 2 + 1));
 759
 760                                 if (is_simple(result)) {
 761                                         im->pix(i, j)[0] = 0;
 762                                         delete result;
 763                                         return NULL;
 764                                 }
 765
 766                                 return result;
 767                         }
 768
 769                         /*
 770                          * Zero weights have NULL subtrees.
 771                          */
 772
 773                         if (im->pix(i, j)[0] == 0)
 774                                 return NULL;
 775
 776                         /*
 777                          * Handle the remaining case.
 778                          */
 779
 780                         return new subtree(im->pix(i, j)[0], NULL, NULL, NULL, NULL);
 781                 }
 782
 783                 subtree *get_subtree(int tc, d2::point p) {
 784                         p *= pow(2, -tc);
 785
 786                         unsigned int i = (unsigned int) floor(p[0]);
 787                         unsigned int j = (unsigned int) floor(p[1]);
 788
 789                         return get_subtree(tc, i, j);
 790                 }
 791
 792                 subtree *get_subtree(const space::traverse &t) {
 793                         ale_pos tc = transformation.trilinear_coordinate(t);
 794                         point p = transformation.centroid(t);
 795                         assert(in_spatial_bounds(p));
 796
 797                         if (!initialized)
 798                                 return NULL;
 799
 800                         if (tc < input_decimation_lower)
 801                                 tc = input_decimation_lower;
 802
 803                         tc = round(tc);
 804
 805                         if (tc < tc_low)
 806                                 return NULL;
 807
 808                         return get_subtree((int) tc, p.xy());
 809                 }
 810
 811                 /*
 812                  * Destructor
 813                  */
 814                 ~ref_weights() {
 815                         for (unsigned int i = 0; i < weights.size(); i++) {
 816                                 delete weights[i];
 817                         }
 818                 }
 819         };
 820
 821         /*
 822          * Resolution check.
 823          */
 824         static int resolution_ok(pt transformation, ale_pos tc) {
 825
 826                 if (pow(2, tc) > transformation.unscaled_height()
 827                  || pow(2, tc) > transformation.unscaled_width())
 828                         return 0;
 829
 830                 if (tc < input_decimation_lower - 1.5)
 831                         return 0;
 832
 833                 return 1;
 834         }
 835
 836         /*
 837          * Structure to hold input frame information at all levels of detail.
 838          */
 839         class lod_images {
 840
 841                 /*
 842                  * All images.
 843                  */
 844
 845                 std::vector<lod_image *> images;
 846
 847         public:
 848
 849                 lod_images() {
 850                         images.resize(d2::image_rw::count(), NULL);
 851                 }
 852
 853                 unsigned int count() {
 854                         return d2::image_rw::count();
 855                 }
 856
 857                 void open(unsigned int f) {
 858                         assert (images[f] == NULL);
 859
 860                         if (images[f] == NULL)
 861                                 images[f] = new lod_image(f);
 862                 }
 863
 864                 void open_all() {
 865                         for (unsigned int f = 0; f < d2::image_rw::count(); f++)
 866                                 open(f);
 867                 }
 868
 869                 lod_image *get(unsigned int f) {
 870                         assert (images[f] != NULL);
 871                         return images[f];
 872                 }
 873
 874                 void close(unsigned int f) {
 875                         assert (images[f] != NULL);
 876                         delete images[f];
 877                         images[f] = NULL;
 878                 }
 879
 880                 void close_all() {
 881                         for (unsigned int f = 0; f < d2::image_rw::count(); f++)
 882                                 close(f);
 883                 }
 884
 885                 ~lod_images() {
 886                         close_all();
 887                 }
 888         };
 889
 890         /*
 891          * All levels-of-detail
 892          */
 893
 894         static struct lod_images *al;
 895
 896         /*
 897          * Data structure for storing best encountered subspace candidates.
 898          */
 899         class candidates {
 900                 std::vector<std::vector<std::pair<ale_pos, ale_real> > > levels;
 901                 int image_index;
 902                 unsigned int height;
 903                 unsigned int width;
 904
 905                 /*
 906                  * Point p is in world coordinates.
 907                  */
 908                 void generate_subspace(point iw, ale_pos diagonal) {
 909
 910 //                      fprintf(stderr, "[gs iw=%f %f %f d=%f]\n",
 911 //                                      iw[0], iw[1], iw[2], diagonal);
 912
 913                         space::traverse st = space::traverse::root();
 914
 915                         if (!st.includes(iw)) {
 916                                 assert(0);
 917                                 return;
 918                         }
 919
 920                         int highres = 0;
 921                         int lowres = 0;
 922
 923                         /*
 924                          * Loop until resolutions of interest have been generated.
 925                          */
 926
 927                         for(;;) {
 928
 929                                 ale_pos current_diagonal = (st.get_max() - st.get_min()).norm();
 930
 931                                 assert(!isnan(current_diagonal));
 932
 933                                 /*
 934                                  * Generate any new desired spatial registers.
 935                                  */
 936
 937                                 /*
 938                                  * Inputs
 939                                  */
 940
 941                                 for (int f = 0; f < 2; f++) {
 942
 943                                         /*
 944                                          * Low resolution
 945                                          */
 946
 947                                         if (current_diagonal < 2 * diagonal
 948                                          && lowres == 0) {
 949                                                 if (spatial_info_map.find(st.get_node()) == spatial_info_map.end()) {
 950                                                         spatial_info_map[st.get_node()];
 951                                                         ui::get()->d3_increment_spaces();
 952                                                 }
 953                                                 lowres = 1;
 954                                         }
 955
 956                                         /*
 957                                          * High resolution.
 958                                          */
 959
 960                                         if (current_diagonal < diagonal
 961                                          && highres == 0) {
 962                                                 if (spatial_info_map.find(st.get_node()) == spatial_info_map.end()) {
 963                                                         spatial_info_map[st.get_node()];
 964                                                         ui::get()->d3_increment_spaces();
 965                                                 }
 966                                                 highres = 1;
 967                                         }
 968                                 }
 969
 970                                 /*
 971                                  * Check for completion
 972                                  */
 973
 974                                 if (highres && lowres)
 975                                         return;
 976
 977                                 /*
 978                                  * Check precision before analyzing space further.
 979                                  */
 980
 981                                 if (st.precision_wall()) {
 982                                         fprintf(stderr, "\n\n*** Error: reached subspace precision wall ***\n\n");
 983                                         assert(0);
 984                                         return;
 985                                 }
 986
 987                                 if (st.positive().includes(iw)) {
 988                                         st = st.positive();
 989                                         total_tsteps++;
 990                                 } else if (st.negative().includes(iw)) {
 991                                         st = st.negative();
 992                                         total_tsteps++;
 993                                 } else {
 994                                         fprintf(stderr, "failed iw = (%f, %f, %f)\n",
 995                                                         iw[0], iw[1], iw[2]);
 996                                         assert(0);
 997                                 }
 998                         }
 999                 }
1000
1001         public:
1002                 candidates(int f) {
1003
1004                         image_index = f;
1005                         height = (unsigned int) al->get(f)->get_t(0).unscaled_height();
1006                         width = (unsigned int) al->get(f)->get_t(0).unscaled_width();
1007
1008                         /*
1009                          * Is this necessary?
1010                          */
1011
1012                         levels.resize(primary_decimation_upper - input_decimation_lower + 1);
1013                         for (int l = input_decimation_lower; l <= primary_decimation_upper; l++) {
1014                                 levels[l - input_decimation_lower].resize((unsigned int) (floor(height / pow(2, l))
1015                                                                * floor(width / pow(2, l))
1016                                                                * pairwise_ambiguity),
1017                                                  std::pair<ale_pos, ale_real>(0, 0));
1018                         }
1019                 }
1020
1021                 /*
1022                  * Point p is expected to be in local projective coordinates.
1023                  */
1024
1025                 void add_candidate(point p, int tc, ale_real score) {
1026                         assert(tc <= primary_decimation_upper);
1027                         assert(tc >= input_decimation_lower);
1028                         assert(p[2] < 0);
1029                         assert(score >= 0);
1030
1031                         int i = (unsigned int) floor(p[0] / pow(2, tc));
1032                         int j = (unsigned int) floor(p[1] / pow(2, tc));
1033
1034                         int swidth  = (int) floor(width / pow(2, tc));
1035
1036                         assert(j < swidth);
1037                         assert(i < (int) floor(height / pow(2, tc)));
1038
1039                         for (unsigned int k = 0; k < pairwise_ambiguity; k++) {
1040                                 std::pair<ale_pos, ale_real> *pk =
1041                                         &(levels[tc - input_decimation_lower][i * swidth * pairwise_ambiguity + j * pairwise_ambiguity + k]);
1042
1043                                 if (pk->first != 0 && score >= pk->second)
1044                                         continue;
1045
1046                                 if (i == 1 && j == 1 && tc == 4)
1047                                         fprintf(stderr, "[ac p2=%f score=%f]\n", p[2], score);
1048
1049                                 ale_pos tp = pk->first;
1050                                 ale_real tr = pk->second;
1051
1052                                 pk->first = p[2];
1053                                 pk->second = score;
1054
1055                                 p[2] = tp;
1056                                 score = tr;
1057
1058                                 if (p[2] == 0)
1059                                         break;
1060                         }
1061                 }
1062
1063                 /*
1064                  * Generate subspaces for candidates.
1065                  */
1066
1067                 void generate_subspaces() {
1068
1069                         fprintf(stderr, "+");
1070                         for (int l = input_decimation_lower; l <= primary_decimation_upper; l++) {
1071                                 unsigned int sheight = (unsigned int) floor(height / pow(2, l));
1072                                 unsigned int swidth  = (unsigned int) floor(width  / pow(2, l));
1073
1074                                 for (unsigned int i = 0; i < sheight; i++)
1075                                 for (unsigned int j = 0; j < swidth; j++)
1076                                 for (unsigned int k = 0; k < pairwise_ambiguity; k++) {
1077                                         std::pair<ale_pos, ale_real> *pk =
1078                                                 &(levels[l - input_decimation_lower]
1079                                                         [i * swidth * pairwise_ambiguity + j * pairwise_ambiguity + k]);
1080
1081                                         if (pk->first == 0) {
1082                                                 fprintf(stderr, "o");
1083                                                 continue;
1084                                         } else {
1085                                                 fprintf(stderr, "|");
1086                                         }
1087
1088                                         ale_pos si = i * pow(2, l) + ((l > 0) ? pow(2, l - 1) : 0);
1089                                         ale_pos sj = j * pow(2, l) + ((l > 0) ? pow(2, l - 1) : 0);
1090
1091 //                                      fprintf(stderr, "[gss l=%d i=%d j=%d d=%g]\n", l, i, j, pk->first);
1092
1093                                         point p = al->get(image_index)->get_t(0).pw_unscaled(point(si, sj, pk->first));
1094
1095                                         generate_subspace(p,
1096                                                         al->get(image_index)->get_t(0).diagonal_distance_3d(pk->first, l));
1097                                 }
1098                         }
1099                 }
1100         };
1101
1102         /*
1103          * List for calculating weighted median.
1104          */
1105         class wml {
1106                 ale_real *data;
1107                 unsigned int size;
1108                 unsigned int used;
1109
1110                 ale_real &_w(unsigned int i) {
1111                         assert(i <= used);
1112                         return data[i * 2];
1113                 }
1114
1115                 ale_real &_d(unsigned int i) {
1116                         assert(i <= used);
1117                         return data[i * 2 + 1];
1118                 }
1119
1120                 void increase_capacity() {
1121
1122                         if (size > 0)
1123                                 size *= 2;
1124                         else
1125                                 size = 1;
1126
1127                         data = (ale_real *) realloc(data, sizeof(ale_real) * 2 * (size * 1));
1128
1129                         assert(data);
1130                         assert (size > used);
1131
1132                         if (!data) {
1133                                 fprintf(stderr, "Unable to allocate %d bytes of memory\n",
1134                                                 sizeof(ale_real) * 2 * (size * 1));
1135                                 exit(1);
1136                         }
1137                 }
1138
1139                 void insert_weight(unsigned int i, ale_real v, ale_real w) {
1140                         assert(used < size);
1141                         assert(used >= i);
1142                         for (unsigned int j = used; j > i; j--) {
1143                                 _w(j) = _w(j - 1);
1144                                 _d(j) = _d(j - 1);
1145                         }
1146
1147                         _w(i) = w;
1148                         _d(i) = v;
1149
1150                         used++;
1151                 }
1152
1153         public:
1154
1155                 unsigned int get_size() {
1156                         return size;
1157                 }
1158
1159                 unsigned int get_used() {
1160                         return used;
1161                 }
1162
1163                 void print_info() {
1164                         fprintf(stderr, "[st %p sz %u el", this, size);
1165                         for (unsigned int i = 0; i < used; i++)
1166                                 fprintf(stderr, " (%f, %f)", _d(i), _w(i));
1167                         fprintf(stderr, "]\n");
1168                 }
1169
1170                 void clear() {
1171                         used = 0;
1172                 }
1173
1174                 void insert_weight(ale_real v, ale_real w) {
1175                         for (unsigned int i = 0; i < used; i++) {
1176                                 if (_d(i) == v) {
1177                                         _w(i) += w;
1178                                         return;
1179                                 }
1180                                 if (_d(i) > v) {
1181                                         if (used == size)
1182                                                 increase_capacity();
1183                                         insert_weight(i, v, w);
1184                                         return;
1185                                 }
1186                         }
1187                         if (used == size)
1188                                 increase_capacity();
1189                         insert_weight(used, v, w);
1190                 }
1191
1192                 /*
1193                  * Finds the median at half-weight, or between half-weight
1194                  * and zero-weight, depending on the attenuation value.
1195                  */
1196
1197                 ale_real find_median(double attenuation = 0) {
1198
1199                         assert(attenuation >= 0);
1200                         assert(attenuation <= 1);
1201
1202                         ale_real zero1 = 0;
1203                         ale_real zero2 = 0;
1204                         ale_real undefined = zero1 / zero2;
1205
1206                         ale_accum weight_sum = 0;
1207
1208                         for (unsigned int i = 0; i < used; i++)
1209                                 weight_sum += _w(i);
1210
1211 //                      if (weight_sum == 0)
1212 //                              return undefined;
1213
1214                         if (used == 0 || used == 1)
1215                                 return undefined;
1216
1217                         if (weight_sum == 0) {
1218                                 ale_accum data_sum = 0;
1219                                 for (unsigned int i = 0; i < used; i++)
1220                                         data_sum += _d(i);
1221                                 return data_sum / used;
1222                         }
1223
1224
1225                         ale_accum midpoint = weight_sum * (0.5 - 0.5 * attenuation);
1226
1227                         ale_accum weight_sum_2 = 0;
1228
1229                         for (unsigned int i = 0; i < used && weight_sum_2 < midpoint; i++) {
1230                                 weight_sum_2 += _w(i);
1231
1232                                 if (weight_sum_2 > midpoint) {
1233                                         return _d(i);
1234                                 } else if (weight_sum_2 == midpoint) {
1235                                         assert (i + 1 < used);
1236                                         return (_d(i) + _d(i + 1)) / 2;
1237                                 }
1238                         }
1239
1240                         return undefined;
1241                 }
1242
1243                 wml(int initial_size = 0) {
1244
1245 //                      if (initial_size == 0) {
1246 //                              initial_size = (int) (d2::image_rw::count() * 1.5);
1247 //                      }
1248
1249                         size = initial_size;
1250                         used = 0;
1251
1252                         if (size > 0) {
1253                                 data = (ale_real *) malloc(size * sizeof(ale_real) * 2);
1254                                 assert(data);
1255                         } else {
1256                                 data = NULL;
1257                         }
1258                 }
1259
1260                 /*
1261                  * copy constructor.  This is required to avoid undesired frees.
1262                  */
1263
1264                 wml(const wml &w) {
1265                         size = w.size;
1266                         used = w.used;
1267                         data = (ale_real *) malloc(size * sizeof(ale_real) * 2);
1268                         assert(data);
1269
1270                         memcpy(data, w.data, size * sizeof(ale_real) * 2);
1271                 }
1272
1273                 ~wml() {
1274                         free(data);
1275                 }
1276         };
1277
1278         /*
1279          * Class for information regarding spatial regions of interest.
1280          *
1281          * This class is configured for convenience in cases where sampling is
1282          * performed using an approximation of the fine:box:1,triangle:2 chain.
1283          * In this case, the *_1 variables would store the fine data and the
1284          * *_2 variables would store the coarse data.  Other uses are also
1285          * possible.
1286          */
1287
1288         class spatial_info {
1289                 /*
1290                  * Map channel value --> weight.
1291                  */
1292                 wml color_weights_1[3];
1293                 wml color_weights_2[3];
1294
1295                 /*
1296                  * Current color.
1297                  */
1298                 d2::pixel color;
1299
1300                 /*
1301                  * Map occupancy value --> weight.
1302                  */
1303                 wml occupancy_weights_1;
1304                 wml occupancy_weights_2;
1305
1306                 /*
1307                  * Current occupancy value.
1308                  */
1309                 ale_real occupancy;
1310
1311                 /*
1312                  * pocc/socc density
1313                  */
1314
1315                 unsigned int pocc_density;
1316                 unsigned int socc_density;
1317
1318                 /*
1319                  * Insert a weight into a list.
1320                  */
1321                 void insert_weight(wml *m, ale_real v, ale_real w) {
1322                         m->insert_weight(v, w);
1323                 }
1324
1325                 /*
1326                  * Find the median of a weighted list.  Uses NaN for undefined.
1327                  */
1328                 ale_real find_median(wml *m, double attenuation = 0) {
1329                         return m->find_median(attenuation);
1330                 }
1331
1332         public:
1333                 /*
1334                  * Constructor.
1335                  */
1336                 spatial_info() {
1337                         color = d2::pixel::zero();
1338                         occupancy = 0;
1339                         pocc_density = 0;
1340                         socc_density = 0;
1341                 }
1342
1343                 /*
1344                  * Accumulate color; primary data set.
1345                  */
1346                 void accumulate_color_1(int f, d2::pixel color, d2::pixel weight) {
1347                         for (int k = 0; k < 3; k++)
1348                                 insert_weight(&color_weights_1[k], color[k], weight[k]);
1349                 }
1350
1351                 /*
1352                  * Accumulate color; secondary data set.
1353                  */
1354                 void accumulate_color_2(d2::pixel color, d2::pixel weight) {
1355                         for (int k = 0; k < 3; k++)
1356                                 insert_weight(&color_weights_2[k], color[k], weight[k]);
1357                 }
1358
1359                 /*
1360                  * Accumulate occupancy; primary data set.
1361                  */
1362                 void accumulate_occupancy_1(int f, ale_real occupancy, ale_real weight) {
1363                         insert_weight(&occupancy_weights_1, occupancy, weight);
1364                 }
1365
1366                 /*
1367                  * Accumulate occupancy; secondary data set.
1368                  */
1369                 void accumulate_occupancy_2(ale_real occupancy, ale_real weight) {
1370                         insert_weight(&occupancy_weights_2, occupancy, weight);
1371
1372                         if (occupancy == 0 || occupancy_weights_2.get_size() < 96)
1373                                 return;
1374
1375                         // fprintf(stderr, "%p updated socc with: %f %f\n", this, occupancy, weight);
1376                         // occupancy_weights_2.print_info();
1377                 }
1378
1379                 /*
1380                  * Update color (and clear accumulation structures).
1381                  */
1382                 void update_color() {
1383                         for (int d = 0; d < 3; d++) {
1384                                 ale_real c = find_median(&color_weights_1[d]);
1385                                 if (isnan(c))
1386                                         c = find_median(&color_weights_2[d]);
1387                                 if (isnan(c))
1388                                         c = 0;
1389
1390                                 color[d] = c;
1391
1392                                 color_weights_1[d].clear();
1393                                 color_weights_2[d].clear();
1394                         }
1395                 }
1396
1397                 /*
1398                  * Update occupancy (and clear accumulation structures).
1399                  */
1400                 void update_occupancy() {
1401                         ale_real o = find_median(&occupancy_weights_1, occ_att);
1402                         if (isnan(o))
1403                                 o = find_median(&occupancy_weights_2, occ_att);
1404                         if (isnan(o))
1405                                 o = 0;
1406
1407                         occupancy = o;
1408
1409                         pocc_density = occupancy_weights_1.get_used();
1410                         socc_density = occupancy_weights_2.get_used();
1411
1412                         occupancy_weights_1.clear();
1413                         occupancy_weights_2.clear();
1414
1415                 }
1416
1417                 /*
1418                  * Get current color.
1419                  */
1420                 d2::pixel get_color() {
1421                         return color;
1422                 }
1423
1424                 /*
1425                  * Get current occupancy.
1426                  */
1427                 ale_real get_occupancy() {
1428                         assert (finite(occupancy));
1429                         return occupancy;
1430                 }
1431
1432                 /*
1433                  * Get primary color density.
1434                  */
1435
1436                 unsigned int get_pocc_density() {
1437                         return pocc_density;
1438                 }
1439
1440                 unsigned int get_socc_density() {
1441                         return socc_density;
1442                 }
1443         };
1444
1445         /*
1446          * Map spatial regions of interest to spatial info structures.  XXX:
1447          * This may get very poor cache behavior in comparison with, say, an
1448          * array.  Unfortunately, there is no immediately obvious array
1449          * representation.  If some kind of array representation were adopted,
1450          * it would probably cluster regions of similar depth from the
1451          * perspective of the typical camera.  In particular, for a
1452          * stereoscopic view, depth ordering for two random points tends to be
1453          * similar between cameras, I think.  Unfortunately, it is never
1454          * identical for all points (unless cameras are co-located).  One
1455          * possible approach would be to order based on, say, camera 0's idea
1456          * of depth.
1457          */
1458
1459 #if !defined(HASH_MAP_GNU) && !defined(HASH_MAP_STD)
1460         typedef std::map<struct space::node *, spatial_info> spatial_info_map_t;
1461 #elif defined(HASH_MAP_GNU)
1462         struct node_hash
1463         {
1464                 size_t operator()(struct space::node *n) const
1465                 {
1466                         return __gnu_cxx::hash<long>()((long) n);
1467                 }
1468         };
1469         typedef __gnu_cxx::hash_map<struct space::node *, spatial_info, node_hash > spatial_info_map_t;
1470 #elif defined(HASH_MAP_STD)
1471         typedef std::hash_map<struct space::node *, spatial_info> spatial_info_map_t;
1472 #endif
1473
1474         static spatial_info_map_t spatial_info_map;
1475
1476 public:
1477
1478         /*
1479          * Debugging variables.
1480          */
1481
1482         static unsigned long total_ambiguity;
1483         static unsigned long total_pixels;
1484         static unsigned long total_divisions;
1485         static unsigned long total_tsteps;
1486
1487         /*
1488          * Member functions
1489          */
1490
1491         static void et(double et_parameter) {
1492                 encounter_threshold = et_parameter;
1493         }
1494
1495         static void dmr(double dmr_parameter) {
1496                 depth_median_radius = dmr_parameter;
1497         }
1498
1499         static void fmr(double fmr_parameter) {
1500                 diff_median_radius = fmr_parameter;
1501         }
1502
1503         static void load_model(const char *name) {
1504                 load_model_name = name;
1505         }
1506
1507         static void save_model(const char *name) {
1508                 save_model_name = name;
1509         }
1510
1511         static void fc(ale_pos fc) {
1512                 front_clip = fc;
1513         }
1514
1515         static void di_upper(ale_pos _dgi) {
1516                 primary_decimation_upper = (int) round(_dgi);
1517         }
1518
1519         static void do_try(ale_pos _dgo) {
1520                 output_decimation_preferred = (int) round(_dgo);
1521         }
1522
1523         static void di_lower(ale_pos _idiv) {
1524                 input_decimation_lower = (int) round(_idiv);
1525         }
1526
1527         static void oc() {
1528                 output_clip = 1;
1529         }
1530
1531         static void no_oc() {
1532                 output_clip = 0;
1533         }
1534
1535         static void rc(ale_pos rc) {
1536                 rear_clip = rc;
1537         }
1538
1539         /*
1540          * Initialize 3D scene from 2D scene, using 2D and 3D alignment
1541          * information.
1542          */
1543         static void init_from_d2() {
1544
1545                 /*
1546                  * Rear clip value of 0 is converted to infinity.
1547                  */
1548
1549                 if (rear_clip == 0) {
1550                         ale_pos one = +1;
1551                         ale_pos zero = +0;
1552
1553                         rear_clip = one / zero;
1554                         assert(isinf(rear_clip) == +1);
1555                 }
1556
1557                 /*
1558                  * Scale and translate clipping plane depths.
1559                  */
1560
1561                 ale_pos cp_scalar = d3::align::projective(0).wc(point(0, 0, 0))[2];
1562
1563                 front_clip = front_clip * cp_scalar - cp_scalar;
1564                 rear_clip = rear_clip * cp_scalar - cp_scalar;
1565
1566                 /*
1567                  * Allocate image structures.
1568                  */
1569
1570                 al = new lod_images;
1571
1572                 if (tc_multiplier != 0) {
1573                         al->open_all();
1574                 }
1575         }
1576
1577         /*
1578          * Perform spatial_info updating on a given subspace, for given
1579          * parameters.
1580          */
1581         static void subspace_info_update(space::iterate si, int f, ref_weights *weights) {
1582
1583                 while(!si.done()) {
1584
1585                         space::traverse st = si.get();
1586
1587                         /*
1588                          * Skip spaces with no color information.
1589                          *
1590                          * XXX: This could be more efficient, perhaps.
1591                          */
1592
1593                         if (spatial_info_map.count(st.get_node()) == 0) {
1594                                 si.next();
1595                                 continue;
1596                         }
1597
1598                         ui::get()->d3_increment_space_num();
1599
1600
1601                         /*
1602                          * Get in-bounds centroid, if one exists.
1603                          */
1604
1605                         point p = al->get(f)->get_t(0).centroid(st);
1606
1607                         if (!p.defined()) {
1608                                 si.next();
1609                                 continue;
1610                         }
1611
1612                         /*
1613                          * Get information on the subspace.
1614                          */
1615
1616                         spatial_info *sn = &spatial_info_map[st.get_node()];
1617                         d2::pixel color = sn->get_color();
1618                         ale_real occupancy = sn->get_occupancy();
1619
1620                         /*
1621                          * Store current weight so we can later check for
1622                          * modification by higher-resolution subspaces.
1623                          */
1624
1625                         ref_weights::subtree *tree = weights->get_subtree(st);
1626
1627                         /*
1628                          * Check for higher resolution subspaces, and
1629                          * update the space iterator.
1630                          */
1631
1632                         if (st.get_node()->positive
1633                          || st.get_node()->negative) {
1634
1635                                 /*
1636                                  * Cleave space for the higher-resolution pass,
1637                                  * skipping the current space, since we will
1638                                  * process that later.
1639                                  */
1640
1641                                 space::iterate cleaved_space = si.cleave();
1642
1643                                 cleaved_space.next();
1644
1645                                 subspace_info_update(cleaved_space, f, weights);
1646
1647                         } else {
1648                                 si.next();
1649                         }
1650
1651                         /*
1652                          * Add new data on the subspace and update weights.
1653                          */
1654
1655                         ale_pos tc = al->get(f)->get_t(0).trilinear_coordinate(st);
1656                         d2::pixel pcolor = al->get(f)->get_tl(p.xy(), tc);
1657                         d2::pixel colordiff = (color - pcolor) * (ale_real) 256;
1658
1659                         if (falloff_exponent != 0) {
1660                                 d2::pixel max_diff = al->get(f)->get_max_diff(p.xy(), tc) * (ale_real) 256;
1661
1662                                 for (int k = 0; k < 3; k++)
1663                                 if (max_diff[k] > 1)
1664                                         colordiff[k] /= pow(max_diff[k], falloff_exponent);
1665                         }
1666
1667                         /*
1668                          * Determine the probability of encounter.
1669                          */
1670
1671                         d2::pixel encounter = d2::pixel(1, 1, 1) * (1 - weights->get_weight(st));
1672
1673                         /*
1674                          * Update weights
1675                          */
1676
1677                         weights->add_weight(st, occupancy, tree);
1678
1679                         /*
1680                          * Delete the subtree, if necessary.
1681                          */
1682
1683                         delete tree;
1684
1685                         /*
1686                          * Check for cases in which the subspace should not be
1687                          * updated.
1688                          */
1689
1690                         if (!resolution_ok(al->get(f)->get_t(0), tc))
1691                                 continue;
1692
1693                         if (d2::render::is_excluded_f(p.xy(), f))
1694                                 continue;
1695
1696                         /*
1697                          * Update subspace.
1698                          */
1699
1700                         sn->accumulate_color_1(f, pcolor, encounter);
1701                         d2::pixel channel_occ = pexp(-colordiff * colordiff);
1702
1703                         ale_accum occ = channel_occ[0];
1704
1705                         for (int k = 1; k < 3; k++)
1706                                 if (channel_occ[k] < occ)
1707                                         occ = channel_occ[k];
1708
1709                         sn->accumulate_occupancy_1(f, occ, encounter[0]);
1710
1711                 }
1712         }
1713
1714         /*
1715          * Run a single iteration of the spatial_info update cycle.
1716          */
1717         static void spatial_info_update() {
1718                 /*
1719                  * Iterate through each frame.
1720                  */
1721                 for (unsigned int f = 0; f < d2::image_rw::count(); f++) {
1722
1723                         ui::get()->d3_occupancy_status(f);
1724
1725                         /*
1726                          * Open the frame and transformation.
1727                          */
1728
1729                         if (tc_multiplier == 0)
1730                                 al->open(f);
1731
1732                         /*
1733                          * Allocate weights data structure for storing encounter
1734                          * probabilities.
1735                          */
1736
1737                         ref_weights *weights = new ref_weights(f);
1738
1739                         /*
1740                          * Call subspace_info_update for the root space.
1741                          */
1742
1743                         subspace_info_update(space::iterate(al->get(f)->origin()), f, weights);
1744
1745                         /*
1746                          * Free weights.
1747                          */
1748
1749                         delete weights;
1750
1751                         /*
1752                          * Close the frame and transformation.
1753                          */
1754
1755                         if (tc_multiplier == 0)
1756                                 al->close(f);
1757                 }
1758
1759                 /*
1760                  * Update all spatial_info structures.
1761                  */
1762                 for (spatial_info_map_t::iterator i = spatial_info_map.begin(); i != spatial_info_map.end(); i++) {
1763                         i->second.update_color();
1764                         i->second.update_occupancy();
1765
1766 //                      d2::pixel color = i->second.get_color();
1767
1768 //                      fprintf(stderr, "space p=%p updated to c=[%f %f %f] o=%f\n",
1769 //                                      i->first, color[0], color[1], color[2],
1770 //                                      i->second.get_occupancy());
1771                 }
1772         }
1773
1774         /*
1775          * Support function for view() and depth().  This function
1776          * always performs exclusion.
1777          */
1778
1779         static const void view_recurse(int type, d2::image *im, d2::image *weights, space::iterate si, pt _pt,
1780                         int prune = 0, d2::point pl = d2::point(0, 0), d2::point ph = d2::point(0, 0)) {
1781                 while (!si.done()) {
1782                         space::traverse st = si.get();
1783
1784                         /*
1785                          * Remove excluded regions.
1786                          */
1787
1788                         if (excluded(st)) {
1789                                 si.cleave();
1790                                 continue;
1791                         }
1792
1793                         /*
1794                          * Prune.
1795                          */
1796
1797                         if (prune && !_pt.check_inclusion_scaled(st, pl, ph)) {
1798                                 si.cleave();
1799                                 continue;
1800                         }
1801
1802                         /*
1803                          * XXX: This could be more efficient, perhaps.
1804                          */
1805
1806                         if (spatial_info_map.count(st.get_node()) == 0) {
1807                                 si.next();
1808                                 continue;
1809                         }
1810
1811                         ui::get()->d3_increment_space_num();
1812
1813                         spatial_info sn = spatial_info_map[st.get_node()];
1814
1815                         /*
1816                          * Get information on the subspace.
1817                          */
1818
1819                         d2::pixel color = sn.get_color();
1820                         // d2::pixel color = d2::pixel(1, 1, 1) * (double) (((unsigned int) (st.get_node()) / sizeof(space)) % 65535);
1821                         ale_real occupancy = sn.get_occupancy();
1822
1823                         /*
1824                          * Determine the view-local bounding box for the
1825                          * subspace.
1826                          */
1827
1828                         point bb[2];
1829
1830                         _pt.get_view_local_bb_scaled(st, bb);
1831
1832                         point min = bb[0];
1833                         point max = bb[1];
1834
1835                         if (prune) {
1836                                 if (min[0] > ph[0]
1837                                  || min[1] > ph[1]
1838                                  || max[0] < pl[0]
1839                                  || max[1] < pl[1]) {
1840                                         si.next();
1841                                         continue;
1842                                 }
1843
1844                                 if (min[0] < pl[0])
1845                                         min[0] = pl[0];
1846                                 if (min[1] < pl[1])
1847                                         min[1] = pl[1];
1848                                 if (max[0] > ph[0])
1849                                         max[0] = ph[0];
1850                                 if (max[1] > ph[1])
1851                                         max[1] = ph[1];
1852
1853                                 min[0] -= pl[0];
1854                                 min[1] -= pl[1];
1855                                 max[0] -= pl[0];
1856                                 max[1] -= pl[1];
1857                         }
1858
1859                         /*
1860                          * Data structure to check modification of weights by
1861                          * higher-resolution subspaces.
1862                          */
1863
1864                         std::queue<d2::pixel> weight_queue;
1865
1866                         /*
1867                          * Check for higher resolution subspaces, and
1868                          * update the space iterator.
1869                          */
1870
1871                         if (st.get_node()->positive
1872                          || st.get_node()->negative) {
1873
1874                                 /*
1875                                  * Store information about current weights,
1876                                  * so we will know which areas have been
1877                                  * covered by higher-resolution subspaces.
1878                                  */
1879
1880                                 for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
1881                                 for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++)
1882                                         weight_queue.push(weights->get_pixel(i, j));
1883
1884                                 /*
1885                                  * Cleave space for the higher-resolution pass,
1886                                  * skipping the current space, since we will
1887                                  * process that afterward.
1888                                  */
1889
1890                                 space::iterate cleaved_space = si.cleave();
1891
1892                                 cleaved_space.next();
1893
1894                                 view_recurse(type, im, weights, cleaved_space, _pt, prune, pl, ph);
1895
1896                         } else {
1897                                 si.next();
1898                         }
1899
1900
1901                         /*
1902                          * Iterate over pixels in the bounding box, finding
1903                          * pixels that intersect the subspace.  XXX: assume
1904                          * for now that all pixels in the bounding box
1905                          * intersect the subspace.
1906                          */
1907
1908                         for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
1909                         for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++) {
1910
1911                                 /*
1912                                  * Check for higher-resolution updates.
1913                                  */
1914
1915                                 if (weight_queue.size()) {
1916                                         if (weight_queue.front() != weights->get_pixel(i, j)) {
1917                                                 weight_queue.pop();
1918                                                 continue;
1919                                         }
1920                                         weight_queue.pop();
1921                                 }
1922
1923                                 /*
1924                                  * Determine the probability of encounter.
1925                                  */
1926
1927                                 d2::pixel encounter = (d2::pixel(1, 1, 1)
1928                                                      - weights->get_pixel(i, j))
1929                                                     * occupancy;
1930
1931                                 /*
1932                                  * Update images.
1933                                  */
1934
1935                                 if (type == 0) {
1936
1937                                         /*
1938                                          * Color view
1939                                          */
1940
1941                                         weights->pix(i, j) += encounter;
1942                                         im->pix(i, j)      += encounter * color;
1943
1944                                 } else if (type == 1) {
1945
1946                                         /*
1947                                          * Weighted (transparent) depth display
1948                                          */
1949
1950                                         ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
1951                                         weights->pix(i, j) += encounter;
1952                                         im->pix(i, j)      += encounter * depth_value;
1953
1954                                 } else if (type == 2) {
1955
1956                                         /*
1957                                          * Ambiguity (ambivalence) measure.
1958                                          */
1959
1960                                         weights->pix(i, j) = d2::pixel(1, 1, 1);
1961                                         im->pix(i, j) += 0.1 * d2::pixel(1, 1, 1);
1962
1963                                 } else if (type == 3) {
1964
1965                                         /*
1966                                          * Closeness measure.
1967                                          */
1968
1969                                         ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
1970                                         if (weights->pix(i, j)[0] == 0) {
1971                                                 weights->pix(i, j) = d2::pixel(1, 1, 1);
1972                                                 im->pix(i, j) = d2::pixel(1, 1, 1) * depth_value;
1973                                         } else if (im->pix(i, j)[2] < depth_value) {
1974                                                 im->pix(i, j) = d2::pixel(1, 1, 1) * depth_value;
1975                                         } else {
1976                                                 continue;
1977                                         }
1978
1979                                 } else if (type == 4) {
1980
1981                                         /*
1982                                          * Weighted (transparent) contribution display
1983                                          */
1984
1985                                         ale_pos contribution_value = sn.get_pocc_density() /* + sn.get_socc_density() */;
1986                                         weights->pix(i, j) += encounter;
1987                                         im->pix(i, j)      += encounter * contribution_value;
1988
1989                                         assert (finite(encounter[0]));
1990                                         assert (finite(contribution_value));
1991
1992                                 } else if (type == 5) {
1993
1994                                         /*
1995                                          * Weighted (transparent) occupancy display
1996                                          */
1997
1998                                         ale_pos contribution_value = occupancy;
1999                                         weights->pix(i, j) += encounter;
2000                                         im->pix(i, j)      += encounter * contribution_value;
2001
2002                                 } else if (type == 6) {
2003
2004                                         /*
2005                                          * (Depth, xres, yres) triple
2006                                          */
2007
2008                                         ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
2009                                         weights->pix(i, j)[0] += encounter[0];
2010                                         if (weights->pix(i, j)[1] < encounter[0]) {
2011                                                 weights->pix(i, j)[1] = encounter[0];
2012                                                 im->pix(i, j)[0] = weights->pix(i, j)[1] * depth_value;
2013                                                 im->pix(i, j)[1] = max[0] - min[0];
2014                                                 im->pix(i, j)[2] = max[1] - min[1];
2015                                         }
2016
2017                                 } else if (type == 7) {
2018
2019                                         /*
2020                                          * (xoff, yoff, 0) triple
2021                                          */
2022
2023                                         weights->pix(i, j)[0] += encounter[0];
2024                                         if (weights->pix(i, j)[1] < encounter[0]) {
2025                                                 weights->pix(i, j)[1] = encounter[0];
2026                                                 im->pix(i, j)[0] = i - min[0];
2027                                                 im->pix(i, j)[1] = j - min[1];
2028                                                 im->pix(i, j)[2] = 0;
2029                                         }
2030
2031                                 } else if (type == 8) {
2032
2033                                         /*
2034                                          * Value = 1 for any intersected space.
2035                                          */
2036
2037                                         weights->pix(i, j) = d2::pixel(1, 1, 1);
2038                                         im->pix(i, j) = d2::pixel(1, 1, 1);
2039
2040                                 } else if (type == 9) {
2041
2042                                         /*
2043                                          * Number of contributions for the nearest space.
2044                                          */
2045
2046                                         if (weights->pix(i, j)[0] == 1)
2047                                                 continue;
2048
2049                                         weights->pix(i, j) = d2::pixel(1, 1, 1);
2050                                         im->pix(i, j) = d2::pixel(1, 1, 1) * (sn.get_pocc_density() * 0.1);
2051
2052                                 } else
2053                                         assert(0);
2054                         }
2055                 }
2056         }
2057
2058         /*
2059          * Generate an depth image from a specified view.
2060          */
2061         static const d2::image *depth(pt _pt, int n = -1, int prune = 0,
2062                         d2::point pl = d2::point(0, 0), d2::point ph = d2::point(0, 0)) {
2063                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
2064
2065                 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
2066
2067                 if (n >= 0) {
2068                         assert((int) floor(d2::align::of(n).scaled_height())
2069                              == (int) floor(_pt.scaled_height()));
2070                         assert((int) floor(d2::align::of(n).scaled_width())
2071                              == (int) floor(_pt.scaled_width()));
2072                 }
2073
2074                 d2::image *im1, *im2, *im3, *weights;;
2075
2076                 if (prune) {
2077
2078                         im1 = new d2::image_ale_real((int) floor(ph[0] - pl[0]) + 1,
2079                                         (int) floor(ph[1] - pl[1]) + 1, 3);
2080
2081                         im2 = new d2::image_ale_real((int) floor(ph[0] - pl[0]) + 1,
2082                                         (int) floor(ph[1] - pl[1]) + 1, 3);
2083
2084                         im3 = new d2::image_ale_real((int) floor(ph[0] - pl[0]) + 1,
2085                                         (int) floor(ph[1] - pl[1]) + 1, 3);
2086
2087                         weights = new d2::image_ale_real((int) floor(ph[0] - pl[0]) + 1,
2088                                         (int) floor(ph[1] - pl[1]) + 1, 3);
2089
2090                 } else {
2091
2092                         im1 = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2093                                        (int) floor(_pt.scaled_width()), 3);
2094
2095                         im2 = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2096                                        (int) floor(_pt.scaled_width()), 3);
2097
2098                         im3 = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2099                                        (int) floor(_pt.scaled_width()), 3);
2100
2101                         weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2102                                                         (int) floor(_pt.scaled_width()), 3);
2103                 }
2104
2105                 /*
2106                  * Iterate through subspaces.
2107                  */
2108
2109                 space::iterate si(_pt.origin());
2110
2111                 view_recurse(6, im1, weights, si, _pt, prune, pl, ph);
2112
2113                 delete weights;
2114
2115                 if (prune) {
2116                         weights = new d2::image_ale_real((int) floor(ph[0] - pl[0]) + 1,
2117                                         (int) floor(ph[1] - pl[1]) + 1, 3);
2118                 } else {
2119                         weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2120                                                         (int) floor(_pt.scaled_width()), 3);
2121                 }
2122
2123 #if 1
2124                 view_recurse(7, im2, weights, si, _pt, prune, pl, ph);
2125 #else
2126                 view_recurse(8, im2, weights, si, _pt, prune, pl, ph);
2127                 return im2;
2128 #endif
2129
2130                 /*
2131                  * Normalize depths by weights
2132                  */
2133
2134                 if (normalize_weights)
2135                 for (unsigned int i = 0; i < im1->height(); i++)
2136                 for (unsigned int j = 0; j < im1->width();  j++)
2137                         im1->pix(i, j)[0] /= weights->pix(i, j)[1];
2138
2139
2140                 for (unsigned int i = 0; i < im1->height(); i++)
2141                 for (unsigned int j = 0; j < im1->width();  j++) {
2142
2143                         /*
2144                          * Handle interpolation.
2145                          */
2146
2147                         d2::point x;
2148                         d2::point blx;
2149                         d2::point res(im1->pix(i, j)[1], im1->pix(i, j)[2]);
2150
2151                         for (int d = 0; d < 2; d++) {
2152
2153                                 if (im2->pix(i, j)[d] < res[d] / 2)
2154                                         x[d] = (ale_pos) (d?j:i) - res[d] / 2 - im2->pix(i, j)[d];
2155                                 else
2156                                         x[d] = (ale_pos) (d?j:i) + res[d] / 2 - im2->pix(i, j)[d];
2157
2158                                 blx[d] = 1 - ((d?j:i) - x[d]) / res[d];
2159                         }
2160
2161                         ale_real depth_val = 0;
2162                         ale_real depth_weight = 0;
2163
2164                         for (int ii = 0; ii < 2; ii++)
2165                         for (int jj = 0; jj < 2; jj++) {
2166                                 d2::point p = x + d2::point(ii, jj) * res;
2167                                 if (im1->in_bounds(p)) {
2168
2169                                         ale_real d = im1->get_bl(p)[0];
2170
2171                                         if (isnan(d))
2172                                                 continue;
2173
2174                                         ale_real w = ((ii ? (1 - blx[0]) : blx[0]) * (jj ? (1 - blx[1]) : blx[1]));
2175                                         depth_weight += w;
2176                                         depth_val += w * d;
2177                                 }
2178                         }
2179
2180                         ale_real depth = depth_val / depth_weight;
2181
2182                         /*
2183                          * Handle encounter thresholds
2184                          */
2185
2186                         if (weights->pix(i, j)[0] < encounter_threshold) {
2187                                 im3->pix(i, j) = d2::pixel::zero() / d2::pixel::zero();
2188                         } else {
2189                                 im3->pix(i, j) = d2::pixel(1, 1, 1) * depth;
2190                         }
2191                 }
2192
2193                 delete weights;
2194                 delete im1;
2195                 delete im2;
2196
2197                 return im3;
2198         }
2199
2200         static const d2::image *depth(unsigned int n) {
2201
2202                 assert (n < d2::image_rw::count());
2203
2204                 pt _pt = align::projective(n);
2205
2206                 return depth(_pt, n);
2207         }
2208
2209
2210         /*
2211          * This function always performs exclusion.
2212          */
2213
2214         static space::node *most_visible_pointwise(d2::pixel *weight, space::iterate si, pt _pt, d2::point p) {
2215
2216                 space::node *result = NULL;
2217
2218                 while (!si.done()) {
2219                         space::traverse st = si.get();
2220
2221                         /*
2222                          * Prune certain regions known to be uninteresting.
2223                          */
2224
2225                         if (excluded(st) || !_pt.check_inclusion_scaled(st, p)) {
2226                                 si.cleave();
2227                                 continue;
2228                         }
2229
2230                         /*
2231                          * XXX: This could be more efficient, perhaps.
2232                          */
2233
2234                         if (spatial_info_map.count(st.get_node()) == 0) {
2235                                 si.next();
2236                                 continue;
2237                         }
2238
2239                         spatial_info sn = spatial_info_map[st.get_node()];
2240
2241                         /*
2242                          * Get information on the subspace.
2243                          */
2244
2245                         ale_real occupancy = sn.get_occupancy();
2246
2247                         /*
2248                          * Preserve current weight in order to check for
2249                          * modification by higher-resolution subspaces.
2250                          */
2251
2252                         d2::pixel old_weight = *weight;
2253
2254                         /*
2255                          * Check for higher resolution subspaces, and
2256                          * update the space iterator.
2257                          */
2258
2259                         if (st.get_node()->positive
2260                          || st.get_node()->negative) {
2261
2262                                 /*
2263                                  * Cleave space for the higher-resolution pass,
2264                                  * skipping the current space, since we will
2265                                  * process that afterward.
2266                                  */
2267
2268                                 space::iterate cleaved_space = si.cleave();
2269
2270                                 cleaved_space.next();
2271
2272                                 space::node *r = most_visible_pointwise(weight, cleaved_space, _pt, p);
2273
2274                                 if (old_weight[1] != (*weight)[1])
2275                                         result = r;
2276
2277                         } else {
2278                                 si.next();
2279                         }
2280
2281
2282                         /*
2283                          * Check for higher-resolution updates.
2284                          */
2285
2286                         if (old_weight != *weight)
2287                                 continue;
2288
2289                         /*
2290                          * Determine the probability of encounter.
2291                          */
2292
2293                         ale_pos encounter = (1 - (*weight)[0]) * occupancy;
2294
2295                         /*
2296                          * (*weight)[0] stores the cumulative weight; (*weight)[1] stores the maximum.
2297                          */
2298
2299                         if (encounter > (*weight)[1]) {
2300                                 result = st.get_node();
2301                                 (*weight)[1] = encounter;
2302                         }
2303
2304                         (*weight)[0] += encounter;
2305                 }
2306
2307                 return result;
2308         }
2309
2310         /*
2311          * This function performs exclusion iff SCALED is true.
2312          */
2313         static void  most_visible_generic(std::vector<space::node *> &results, d2::image *weights,
2314                         space::iterate si, pt _pt, int scaled) {
2315
2316                 assert (results.size() == weights->height() * weights->width());
2317
2318                 while (!si.done()) {
2319                         space::traverse st = si.get();
2320
2321                         if (scaled && excluded(st)) {
2322                                 si.cleave();
2323                                 continue;
2324                         }
2325
2326                         /*
2327                          * XXX: This could be more efficient, perhaps.
2328                          */
2329
2330                         if (spatial_info_map.count(st.get_node()) == 0) {
2331                                 si.next();
2332                                 continue;
2333                         }
2334
2335                         spatial_info sn = spatial_info_map[st.get_node()];
2336
2337                         /*
2338                          * Get information on the subspace.
2339                          */
2340
2341                         ale_real occupancy = sn.get_occupancy();
2342
2343                         /*
2344                          * Determine the view-local bounding box for the
2345                          * subspace.
2346                          */
2347
2348                         point bb[2];
2349
2350                         _pt.get_view_local_bb_scaled(st, bb);
2351
2352                         point min = bb[0];
2353                         point max = bb[1];
2354
2355                         /*
2356                          * Data structure to check modification of weights by
2357                          * higher-resolution subspaces.
2358                          */
2359
2360                         std::queue<d2::pixel> weight_queue;
2361
2362                         /*
2363                          * Check for higher resolution subspaces, and
2364                          * update the space iterator.
2365                          */
2366
2367                         if (st.get_node()->positive
2368                          || st.get_node()->negative) {
2369
2370                                 /*
2371                                  * Store information about current weights,
2372                                  * so we will know which areas have been
2373                                  * covered by higher-resolution subspaces.
2374                                  */
2375
2376                                 for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
2377                                 for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++)
2378                                         weight_queue.push(weights->get_pixel(i, j));
2379
2380                                 /*
2381                                  * Cleave space for the higher-resolution pass,
2382                                  * skipping the current space, since we will
2383                                  * process that afterward.
2384                                  */
2385
2386                                 space::iterate cleaved_space = si.cleave();
2387
2388                                 cleaved_space.next();
2389
2390                                 most_visible_generic(results, weights, cleaved_space, _pt, scaled);
2391
2392                         } else {
2393                                 si.next();
2394                         }
2395
2396
2397                         /*
2398                          * Iterate over pixels in the bounding box, finding
2399                          * pixels that intersect the subspace.  XXX: assume
2400                          * for now that all pixels in the bounding box
2401                          * intersect the subspace.
2402                          */
2403
2404                         for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
2405                         for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++) {
2406
2407                                 /*
2408                                  * Check for higher-resolution updates.
2409                                  */
2410
2411                                 if (weight_queue.size()) {
2412                                         if (weight_queue.front() != weights->get_pixel(i, j)) {
2413                                                 weight_queue.pop();
2414                                                 continue;
2415                                         }
2416                                         weight_queue.pop();
2417                                 }
2418
2419                                 /*
2420                                  * Determine the probability of encounter.
2421                                  */
2422
2423                                 ale_pos encounter = (1 - weights->get_pixel(i, j)[0]) * occupancy;
2424
2425                                 /*
2426                                  * weights[0] stores the cumulative weight; weights[1] stores the maximum.
2427                                  */
2428
2429                                 if (encounter > weights->get_pixel(i, j)[1]
2430                                  || results[i * weights->width() + j] == NULL) {
2431                                         results[i * weights->width() + j] = st.get_node();
2432                                         weights->chan(i, j, 1) = encounter;
2433                                 }
2434
2435                                 weights->chan(i, j, 0) += encounter;
2436                         }
2437                 }
2438         }
2439
2440         static std::vector<space::node *> most_visible_scaled(pt _pt) {
2441                 d2::image *weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2442                                 (int) floor(_pt.scaled_width()), 3);
2443                 std::vector<space::node *> results;
2444
2445                 results.resize(weights->height() * weights->width(), 0);
2446
2447                 most_visible_generic(results, weights, space::iterate(_pt.origin()), _pt, 1);
2448
2449                 return results;
2450         }
2451
2452         static std::vector<space::node *> most_visible_unscaled(pt _pt) {
2453                 d2::image *weights = new d2::image_ale_real((int) floor(_pt.unscaled_height()),
2454                                 (int) floor(_pt.unscaled_width()), 3);
2455                 std::vector<space::node *> results;
2456
2457                 results.resize(weights->height() * weights->width(), 0);
2458
2459                 most_visible_generic(results, weights, space::iterate(_pt.origin()), _pt, 0);
2460
2461                 return results;
2462         }
2463
2464         static const int visibility_search(const std::vector<space::node *> &fmv, space::node *mv) {
2465
2466                 if (mv == NULL)
2467                         return 0;
2468
2469                 if (std::binary_search(fmv.begin(), fmv.end(), mv))
2470                         return 1;
2471
2472                 return (visibility_search(fmv, mv->positive)
2473                      || visibility_search(fmv, mv->negative));
2474
2475         }
2476
2477         /*
2478          * Class to generate focal sample views.
2479          */
2480
2481         class view_generator {
2482
2483                 /*
2484                  * Original projective transformation.
2485                  */
2486
2487                 pt original_pt;
2488
2489                 /*
2490                  * Data type for shared view data.
2491                  */
2492
2493                 class shared_view {
2494                         pt _pt;
2495                         std::vector<space::node *> mv;
2496                         d2::image *color;
2497                         d2::image *color_weights;
2498                         const d2::image *_depth;
2499                         d2::image *median_depth;
2500                         d2::image *median_diff;
2501
2502                 public:
2503                         shared_view(pt _pt) {
2504                                 this->_pt = _pt;
2505                                 color = NULL;
2506                                 color_weights = NULL;
2507                                 _depth = NULL;
2508                                 median_depth = NULL;
2509                                 median_diff = NULL;
2510                         }
2511
2512                         shared_view(const shared_view &copy_origin) {
2513                                 _pt = copy_origin._pt;
2514                                 mv = copy_origin.mv;
2515                                 color = NULL;
2516                                 color_weights = NULL;
2517                                 _depth = NULL;
2518                                 median_depth = NULL;
2519                                 median_diff = NULL;
2520                         }
2521
2522                         ~shared_view() {
2523                                 delete color;
2524                                 delete _depth;
2525                                 delete color_weights;
2526                                 delete median_diff;
2527                                 delete median_depth;
2528                         }
2529
2530                         void get_view_recurse(d2::image *data, d2::image *weights, int type) {
2531                                 /*
2532                                  * Iterate through subspaces.
2533                                  */
2534
2535                                 space::iterate si(_pt.origin());
2536
2537                                 ui::get()->d3_render_status(0, 0, -1, -1, -1, -1, 0);
2538
2539                                 view_recurse(type, data, weights, si, _pt);
2540                         }
2541
2542                         void init_color() {
2543                                 color = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2544                                                                 (int) floor(_pt.scaled_width()), 3);
2545
2546                                 color_weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2547                                                                 (int) floor(_pt.scaled_width()), 3);
2548
2549                                 get_view_recurse(color, color_weights, 0);
2550                         }
2551
2552                         void init_depth() {
2553                                 _depth = depth(_pt, -1);
2554                         }
2555
2556                         void init_medians() {
2557                                 if (!_depth)
2558                                         init_depth();
2559
2560                                 assert(_depth);
2561
2562                                 median_diff = _depth->fcdiff_median((int) floor(diff_median_radius));
2563                                 median_depth = _depth->medians((int) floor(depth_median_radius));
2564
2565                                 assert(median_diff);
2566                                 assert(median_depth);
2567                         }
2568
2569                 public:
2570                         pt get_pt() {
2571                                 return _pt;
2572                         }
2573
2574                         space::node *get_most_visible(unsigned int i, unsigned int j) {
2575                                 unsigned int height = (int) floor(_pt.scaled_height());
2576                                 unsigned int width  = (int) floor(_pt.scaled_width());
2577
2578                                 if (i < 0 || i >= height
2579                                  || j < 0 || j >= width) {
2580                                         return NULL;
2581                                 }
2582
2583                                 if (mv.size() == 0) {
2584                                         mv = most_visible_scaled(_pt);
2585                                 }
2586
2587                                 assert (mv.size() > i * width + j);
2588
2589                                 return mv[i * width + j];
2590                         }
2591
2592                         space::node *get_most_visible(d2::point p) {
2593                                 unsigned int i = (unsigned int) round (p[0]);
2594                                 unsigned int j = (unsigned int) round (p[1]);
2595
2596                                 return get_most_visible(i, j);
2597                         }
2598
2599                         d2::pixel get_color(unsigned int i, unsigned int j) {
2600                                 if (color == NULL) {
2601                                         init_color();
2602                                 }
2603
2604                                 assert (color != NULL);
2605
2606                                 return color->get_pixel(i, j);
2607                         }
2608
2609                         d2::pixel get_depth(unsigned int i, unsigned int j) {
2610                                 if (_depth == NULL) {
2611                                         init_depth();
2612                                 }
2613
2614                                 assert (_depth != NULL);
2615
2616                                 return _depth->get_pixel(i, j);
2617                         }
2618
2619                         void get_median_depth_and_diff(d2::pixel *t, d2::pixel *f, unsigned int i, unsigned int j) {
2620                                 if (median_depth == NULL && median_diff == NULL)
2621                                         init_medians();
2622
2623                                 assert (median_depth && median_diff);
2624
2625                                 if (i < 0 || i >= median_depth->height()
2626                                  || j < 0 || j >= median_depth->width()) {
2627                                         *t = d2::pixel::undefined();
2628                                         *f = d2::pixel::undefined();
2629                                 } else {
2630                                         *t = median_depth->get_pixel(i, j);
2631                                         *f = median_diff->get_pixel(i, j);
2632                                 }
2633                         }
2634
2635                         void get_color_and_weight(d2::pixel *c, d2::pixel *w, d2::point p) {
2636                                 if (color == NULL) {
2637                                         init_color();
2638                                 }
2639
2640                                 assert (color != NULL);
2641
2642                                 if (!color->in_bounds(p)) {
2643                                         *c = d2::pixel::undefined();
2644                                         *w = d2::pixel::undefined();
2645                                 } else {
2646                                         *c = color->get_bl(p);
2647                                         *w = color_weights->get_bl(p);
2648                                 }
2649                         }
2650
2651                         d2::pixel get_depth(d2::point p) {
2652                                 if (_depth == NULL) {
2653                                         init_depth();
2654                                 }
2655
2656                                 assert (_depth != NULL);
2657
2658                                 if (!_depth->in_bounds(p)) {
2659                                         return d2::pixel::undefined();
2660                                 }
2661
2662                                 return _depth->get_bl(p);
2663                         }
2664
2665                         void get_median_depth_and_diff(d2::pixel *t, d2::pixel *f, d2::point p) {
2666                                 if (median_diff == NULL && median_depth == NULL)
2667                                         init_medians();
2668
2669                                 assert (median_diff != NULL && median_depth != NULL);
2670
2671                                 if (!median_diff->in_bounds(p)) {
2672                                         *t = d2::pixel::undefined();
2673                                         *f = d2::pixel::undefined();
2674                                 } else {
2675                                         *t = median_depth->get_bl(p);
2676                                         *f = median_diff->get_bl(p);
2677                                 }
2678                         }
2679
2680                 };
2681
2682                 /*
2683                  * Shared view array, indexed by aperture diameter and view number.
2684                  */
2685
2686                 std::map<ale_pos, std::vector<shared_view> > aperture_to_shared_views_map;
2687
2688                 /*
2689                  * Method to generate a new stochastic focal view.
2690                  */
2691
2692                 pt get_new_view(ale_pos aperture) {
2693
2694                         ale_pos ofx = aperture;
2695                         ale_pos ofy = aperture;
2696
2697                         while (ofx * ofx + ofy * ofy > aperture * aperture / 4) {
2698                                 ofx = (rand() * aperture) / RAND_MAX - aperture / 2;
2699                                 ofy = (rand() * aperture) / RAND_MAX - aperture / 2;
2700                         }
2701
2702                         /*
2703                          * Generate a new view from the given offset.
2704                          */
2705
2706                         point new_view = original_pt.cw(point(ofx, ofy, 0));
2707                         pt _pt_new = original_pt;
2708                         for (int d = 0; d < 3; d++)
2709                                 _pt_new.e().set_translation(d, -new_view[d]);
2710
2711                         return _pt_new;
2712                 }
2713
2714         public:
2715
2716                 /*
2717                  * Result type.
2718                  */
2719
2720                 class view {
2721                         shared_view *sv;
2722                         pt _pt;
2723
2724                 public:
2725
2726                         view(shared_view *sv, pt _pt = pt()) {
2727                                 this->sv = sv;
2728                                 if (sv) {
2729                                         this->_pt = sv->get_pt();
2730                                 } else {
2731                                         this->_pt = _pt;
2732                                 }
2733                         }
2734
2735                         pt get_pt() {
2736                                 return _pt;
2737                         }
2738
2739                         space::node *get_most_visible(unsigned int i, unsigned int j) {
2740                                 assert (sv);
2741                                 return sv->get_most_visible(i, j);
2742                         }
2743
2744                         space::node *get_most_visible(d2::point p) {
2745                                 if (sv) {
2746                                         return sv->get_most_visible(p);
2747                                 }
2748
2749                                 d2::pixel weight(0, 0, 0);
2750
2751                                 return most_visible_pointwise(&weight, space::iterate(_pt.origin()), _pt, p);
2752
2753                         }
2754
2755                         d2::pixel get_color(unsigned int i, unsigned int j) {
2756                                 return sv->get_color(i, j);
2757                         }
2758
2759                         void get_color_and_weight(d2::pixel *color, d2::pixel *weight, d2::point p) {
2760                                 if (sv) {
2761                                         sv->get_color_and_weight(color, weight, p);
2762                                         return;
2763                                 }
2764
2765                                 /*
2766                                  * Determine weight and color for the given point.
2767                                  */
2768
2769                                 d2::image *im_point = new d2::image_ale_real(1, 1, 3);
2770                                 d2::image *wt_point = new d2::image_ale_real(1, 1, 3);
2771
2772                                 view_recurse(0, im_point, wt_point, space::iterate(_pt.origin()), _pt, 1, p, p);
2773
2774                                 *color = im_point->pix(0, 0);
2775                                 *weight = wt_point->pix(0, 0);
2776
2777                                 delete im_point;
2778                                 delete wt_point;
2779
2780                                 return;
2781                         }
2782
2783                         d2::pixel get_depth(unsigned int i, unsigned int j) {
2784                                 assert(sv);
2785                                 return sv->get_depth(i, j);
2786                         }
2787
2788                         void get_median_depth_and_diff(d2::pixel *depth, d2::pixel *diff, unsigned int i, unsigned int j) {
2789                                 assert(sv);
2790                                 sv->get_median_depth_and_diff(depth, diff, i, j);
2791                         }
2792
2793                         void get_median_depth_and_diff(d2::pixel *_depth, d2::pixel *_diff, d2::point p) {
2794                                 if (sv) {
2795                                         sv->get_median_depth_and_diff(_depth, _diff, p);
2796                                         return;
2797                                 }
2798
2799                                 /*
2800                                  * Generate a local depth image of required radius.
2801                                  */
2802
2803                                 ale_pos radius = 1;
2804
2805                                 if (diff_median_radius + 1 > radius)
2806                                         radius = diff_median_radius + 1;
2807                                 if (depth_median_radius > radius)
2808                                         radius = depth_median_radius;
2809
2810                                 d2::point pl = p - d2::point(radius, radius);
2811                                 d2::point ph = p + d2::point(radius, radius);
2812                                 const d2::image *local_depth = depth(_pt, -1, 1, pl, ph);
2813
2814                                 /*
2815                                  * Find depth and diff at this point, check for
2816                                  * undefined values, and generate projections
2817                                  * of the image corners on the estimated normal
2818                                  * surface.
2819                                  */
2820
2821                                 d2::image *median_diffs = local_depth->fcdiff_median((int) floor(diff_median_radius));
2822                                 d2::image *median_depths = local_depth->medians((int) floor(depth_median_radius));
2823
2824                                 *_depth = median_depths->pix((int) radius, (int) radius);
2825                                 *_diff = median_diffs->pix((int) radius, (int) radius);
2826
2827                                 delete median_diffs;
2828                                 delete median_depths;
2829                                 delete local_depth;
2830                         }
2831                 };
2832
2833                 view get_view(ale_pos aperture, unsigned index, unsigned int randomization) {
2834                         if (randomization == 0) {
2835
2836                                 while (aperture_to_shared_views_map[aperture].size() <= index) {
2837                                         aperture_to_shared_views_map[aperture].push_back(shared_view(get_new_view(aperture)));
2838                                 }
2839
2840                                 return view(&(aperture_to_shared_views_map[aperture][index]));
2841                         }
2842
2843                         return view(NULL, get_new_view(aperture));
2844                 }
2845
2846                 view_generator(pt original_pt) {
2847                         this->original_pt = original_pt;
2848                 }
2849         };
2850
2851         /*
2852          * Unfiltered function
2853          */
2854         static const d2::image *view_nofilter_focus(pt _pt, int n) {
2855
2856                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
2857
2858                 if (n >= 0) {
2859                         assert((int) floor(d2::align::of(n).scaled_height())
2860                              == (int) floor(_pt.scaled_height()));
2861                         assert((int) floor(d2::align::of(n).scaled_width())
2862                              == (int) floor(_pt.scaled_width()));
2863                 }
2864
2865                 const d2::image *depths = depth(_pt, n);
2866
2867                 d2::image *im = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2868                                                (int) floor(_pt.scaled_width()), 3);
2869
2870                 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
2871
2872                 view_generator vg(_pt);
2873
2874                 for (unsigned int i = 0; i < im->height(); i++)
2875                 for (unsigned int j = 0; j < im->width();  j++) {
2876
2877                         focus::result _focus = focus::get(depths, i, j);
2878
2879                         if (!finite(_focus.focal_distance))
2880                                 continue;
2881
2882                         /*
2883                          * Data structures for calculating focal statistics.
2884                          */
2885
2886                         d2::pixel color, weight;
2887                         d2::image_weighted_median *iwm = NULL;
2888
2889                         if (_focus.statistic == 1) {
2890                                 iwm = new d2::image_weighted_median(1, 1, 3, _focus.sample_count);
2891                         }
2892
2893                         /*
2894                          * Iterate over views for this focus region.
2895                          */
2896
2897                         for (unsigned int v = 0; v < _focus.sample_count; v++) {
2898
2899                                 view_generator::view vw = vg.get_view(_focus.aperture, v, _focus.randomization);
2900
2901                                 ui::get()->d3_render_status(0, 1, -1, v, i, j, -1);
2902
2903
2904                                 /*
2905                                  * Map the focused point to the new view.
2906                                  */
2907
2908                                 point p = vw.get_pt().wp_scaled(_pt.pw_scaled(point(i, j, _focus.focal_distance)));
2909
2910                                 /*
2911                                  * Determine weight and color for the given point.
2912                                  */
2913
2914                                 d2::pixel view_weight, view_color;
2915
2916                                 vw.get_color_and_weight(&view_color, &view_weight, p.xy());
2917
2918                                 if (!color.finite() || !weight.finite())
2919                                         continue;
2920
2921                                 if (_focus.statistic == 0) {
2922                                         color += view_color;
2923                                         weight += view_weight;
2924                                 } else if (_focus.statistic == 1) {
2925                                         iwm->accumulate(0, 0, v, view_color, view_weight);
2926                                 } else
2927                                         assert(0);
2928                         }
2929
2930                         if (_focus.statistic == 1) {
2931                                 weight = iwm->get_weights()->get_pixel(0, 0);
2932                                 color = iwm->get_pixel(0, 0);
2933                                 delete iwm;
2934                         }
2935
2936                         if (weight.min_norm() < encounter_threshold) {
2937                                 im->pix(i, j) = d2::pixel::zero() / d2::pixel::zero();
2938                         } else if (normalize_weights)
2939                                 im->pix(i, j) = color / weight;
2940                         else
2941                                 im->pix(i, j) = color;
2942                 }
2943
2944                 delete depths;
2945
2946                 return im;
2947         }
2948
2949         /*
2950          * Unfiltered function
2951          */
2952         static const d2::image *view_nofilter(pt _pt, int n) {
2953
2954                 if (!focus::is_trivial())
2955                         return view_nofilter_focus(_pt, n);
2956
2957                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
2958
2959                 if (n >= 0) {
2960                         assert((int) floor(d2::align::of(n).scaled_height())
2961                              == (int) floor(_pt.scaled_height()));
2962                         assert((int) floor(d2::align::of(n).scaled_width())
2963                              == (int) floor(_pt.scaled_width()));
2964                 }
2965
2966                 const d2::image *depths = depth(_pt, n);
2967
2968                 d2::image *im = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2969                                                (int) floor(_pt.scaled_width()), 3);
2970
2971                 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
2972
2973                 /*
2974                  * Use adaptive subspace data.
2975                  */
2976
2977                 d2::image *weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2978                                                 (int) floor(_pt.scaled_width()), 3);
2979
2980                 /*
2981                  * Iterate through subspaces.
2982                  */
2983
2984                 space::iterate si(_pt.origin());
2985
2986                 ui::get()->d3_render_status(0, 0, -1, -1, -1, -1, 0);
2987
2988                 view_recurse(0, im, weights, si, _pt);
2989
2990                 for (unsigned int i = 0; i < im->height(); i++)
2991                 for (unsigned int j = 0; j < im->width();  j++) {
2992                         if (weights->pix(i, j).min_norm() < encounter_threshold
2993                          || (d3px_count > 0 && isnan(depths->pix(i, j)[0]))) {
2994                                 im->pix(i, j) = d2::pixel::zero() / d2::pixel::zero();
2995                                 weights->pix(i, j) = d2::pixel::zero();
2996                         } else if (normalize_weights)
2997                                 im->pix(i, j) /= weights->pix(i, j);
2998                 }
2999
3000                 delete weights;
3001
3002                 delete depths;
3003
3004                 return im;
3005         }
3006
3007         /*
3008          * Filtered function.
3009          */
3010         static const d2::image *view_filter_focus(pt _pt, int n) {
3011
3012                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
3013
3014                 /*
3015                  * Get depth image for focus region determination.
3016                  */
3017
3018                 const d2::image *depths = depth(_pt, n);
3019
3020                 unsigned int height = (unsigned int) floor(_pt.scaled_height());
3021                 unsigned int width = (unsigned int) floor(_pt.scaled_width());
3022
3023                 /*
3024                  * Prepare input frame data.
3025                  */
3026
3027                 if (tc_multiplier == 0)
3028                         al->open_all();
3029
3030                 pt *_ptf = new pt[al->count()];
3031                 std::vector<space::node *> *fmv = new std::vector<space::node *>[al->count()];
3032
3033                 for (unsigned int f = 0; f < al->count(); f++) {
3034                         _ptf[f] = al->get(f)->get_t(0);
3035                         fmv[f] = most_visible_unscaled(_ptf[f]);
3036                         std::sort(fmv[f].begin(), fmv[f].end());
3037                 }
3038
3039                 if (tc_multiplier == 0)
3040                         al->close_all();
3041
3042                 /*
3043                  * Open all files for rendering.
3044                  */
3045
3046                 d2::image_rw::open_all();
3047
3048                 /*
3049                  * Prepare data structures for averaging views, as we render
3050                  * each view separately.  This is spacewise inefficient, but
3051                  * is easy to implement given the current operation of the
3052                  * renderers.
3053                  */
3054
3055                 d2::image_weighted_avg *iwa;
3056
3057                 if (d3::focus::uses_medians()) {
3058                         iwa = new d2::image_weighted_median(height, width, 3, focus::max_samples());
3059                 } else {
3060                         iwa = new d2::image_weighted_simple(height, width, 3, new d2::invariant(NULL));
3061                 }
3062
3063                 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
3064
3065                 /*
3066                  * Prepare view generator.
3067                  */
3068
3069                 view_generator vg(_pt);
3070
3071                 /*
3072                  * Render views separately.  This is spacewise inefficient,
3073                  * but is easy to implement given the current operation of the
3074                  * renderers.
3075                  */
3076
3077                 for (unsigned int v = 0; v < focus::max_samples(); v++) {
3078
3079                         /*
3080                          * Generate a new 2D renderer for filtering.
3081                          */
3082
3083                         d2::render::reset();
3084                         d2::render *renderer = d2::render_parse::get(d3chain_type);
3085
3086                         renderer->init_point_renderer(height, width, 3);
3087
3088                         /*
3089                          * Iterate over output points.
3090                          */
3091
3092                         for (unsigned int i = 0; i < height; i++)
3093                         for (unsigned int j = 0; j < width; j++) {
3094
3095                                 focus::result _focus = focus::get(depths, i, j);
3096
3097                                 if (v >= _focus.sample_count)
3098                                         continue;
3099
3100                                 if (!finite(_focus.focal_distance))
3101                                         continue;
3102
3103                                 view_generator::view vw = vg.get_view(_focus.aperture, v, _focus.randomization);
3104
3105                                 pt _pt_new = vw.get_pt();
3106
3107                                 point p = _pt_new.wp_scaled(_pt.pw_scaled(point(i, j, _focus.focal_distance)));
3108
3109                                 /*
3110                                  * Determine the most-visible subspace.
3111                                  */
3112
3113                                 space::node *mv = vw.get_most_visible(p.xy());
3114
3115                                 if (mv == NULL)
3116                                         continue;
3117
3118                                 /*
3119                                  * Get median depth and diff.
3120                                  */
3121
3122                                 d2::pixel depth, diff;
3123
3124                                 vw.get_median_depth_and_diff(&depth, &diff, p.xy());
3125
3126                                 if (!depth.finite() || !diff.finite())
3127                                         continue;
3128
3129                                 point local_points[3] = {
3130                                         point(p[0],     p[1],     depth[0]),
3131                                         point(p[0] + 1, p[1],     depth[0] + diff[0]),
3132                                         point(p[0],     p[1] + 1, depth[0] + diff[1])
3133                                 };
3134
3135                                 /*
3136                                  * Iterate over files.
3137                                  */
3138
3139                                 for (unsigned int f = 0; f < d2::image_rw::count(); f++) {
3140
3141                                         ui::get()->d3_render_status(1, 1, f, v, i, j, -1);
3142
3143                                         if (!visibility_search(fmv[f], mv))
3144                                                 continue;
3145
3146                                         /*
3147                                          * Determine transformation at (i, j).  First
3148                                          * determine transformation from the output to
3149                                          * the input, then invert this, as we need the
3150                                          * inverse transformation for filtering.
3151                                          */
3152
3153                                         d2::point remote_points[3] = {
3154                                                 _ptf[f].wp_unscaled(_pt_new.pw_scaled(point(local_points[0]))).xy(),
3155                                                 _ptf[f].wp_unscaled(_pt_new.pw_scaled(point(local_points[1]))).xy(),
3156                                                 _ptf[f].wp_unscaled(_pt_new.pw_scaled(point(local_points[2]))).xy()
3157                                         };
3158
3159                                         /*
3160                                          * Forward matrix for the linear component of the
3161                                          * transformation.
3162                                          */
3163
3164                                         d2::point forward_matrix[2] = {
3165                                                 remote_points[1] - remote_points[0],
3166                                                 remote_points[2] - remote_points[0]
3167                                         };
3168
3169                                         /*
3170                                          * Inverse matrix for the linear component of
3171                                          * the transformation.  Calculate using the
3172                                          * determinant D.
3173                                          */
3174
3175                                         ale_pos D = forward_matrix[0][0] * forward_matrix[1][1]
3176                                                   - forward_matrix[0][1] * forward_matrix[1][0];
3177
3178                                         if (D == 0)
3179                                                 continue;
3180
3181                                         d2::point inverse_matrix[2] = {
3182                                                 d2::point( forward_matrix[1][1] / D, -forward_matrix[1][0] / D),
3183                                                 d2::point(-forward_matrix[0][1] / D,  forward_matrix[0][0] / D)
3184                                         };
3185
3186                                         /*
3187                                          * Determine the projective transformation parameters for the
3188                                          * inverse transformation.
3189                                          */
3190
3191                                         const d2::image *imf = d2::image_rw::get_open(f);
3192
3193                                         d2::transformation inv_t = d2::transformation::gpt_identity(imf, 1);
3194
3195                                         d2::point local_bounds[4];
3196
3197                                         for (int n = 0; n < 4; n++) {
3198                                                 d2::point remote_bound = d2::point((n == 1 || n == 2) ? imf->height() : 0,
3199                                                                                    (n == 2 || n == 3) ? imf->width()  : 0)
3200                                                                        - remote_points[0];
3201
3202                                                 local_bounds[n] = d2::point(i, j)
3203                                                                 + d2::point(remote_bound[0] * inverse_matrix[0][0]
3204                                                                           + remote_bound[1] * inverse_matrix[1][0],
3205                                                                             remote_bound[0] * inverse_matrix[0][1]
3206                                                                           + remote_bound[1] * inverse_matrix[1][1]);
3207
3208                                         }
3209
3210                                         if (!local_bounds[0].finite()
3211                                          || !local_bounds[1].finite()
3212                                          || !local_bounds[2].finite()
3213                                          || !local_bounds[3].finite())
3214                                                 continue;
3215
3216                                         inv_t.gpt_set(local_bounds);
3217
3218                                         /*
3219                                          * Perform render step for the given frame,
3220                                          * transformation, and point.
3221                                          */
3222
3223                                         renderer->point_render(i, j, f, inv_t);
3224                                 }
3225                         }
3226
3227                         renderer->finish_point_rendering();
3228
3229                         const d2::image *im = renderer->get_image();
3230                         const d2::image *df = renderer->get_defined();
3231
3232                         for (unsigned int i = 0; i < height; i++)
3233                         for (unsigned int j = 0; j < width; j++) {
3234                                 if (df->get_pixel(i, j).finite()
3235                                  && df->get_pixel(i, j)[0] > 0)
3236                                         iwa->accumulate(i, j, v, im->get_pixel(i, j), d2::pixel(1, 1, 1));
3237                         }
3238                 }
3239
3240                 /*
3241                  * Close all files and return the result.
3242                  */
3243
3244                 d2::image_rw::close_all();
3245
3246                 return iwa;
3247         }
3248
3249         static const d2::image *view_filter(pt _pt, int n) {
3250
3251                 if (!focus::is_trivial())
3252                         return view_filter_focus(_pt, n);
3253
3254                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
3255
3256                 /*
3257                  * Generate a new 2D renderer for filtering.
3258                  */
3259
3260                 d2::render::reset();
3261                 d2::render *renderer = d2::render_parse::get(d3chain_type);
3262
3263                 /*
3264                  * Get depth image in order to estimate normals (and hence
3265                  * transformations).
3266                  */
3267
3268                 const d2::image *depths = depth(_pt, n);
3269
3270                 d2::image *median_diffs = depths->fcdiff_median((int) floor(diff_median_radius));
3271                 d2::image *median_depths = depths->medians((int) floor(depth_median_radius));
3272
3273                 unsigned int height = (unsigned int) floor(_pt.scaled_height());
3274                 unsigned int width = (unsigned int) floor(_pt.scaled_width());
3275
3276                 renderer->init_point_renderer(height, width, 3);
3277
3278                 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
3279
3280                 std::vector<space::node *> mv = most_visible_scaled(_pt);
3281
3282                 for (unsigned int f = 0; f < d2::image_rw::count(); f++) {
3283
3284                         if (tc_multiplier == 0)
3285                                 al->open(f);
3286
3287                         pt _ptf = al->get(f)->get_t(0);
3288
3289                         std::vector<space::node *> fmv = most_visible_unscaled(_ptf);
3290                         std::sort(fmv.begin(), fmv.end());
3291
3292                         for (unsigned int i = 0; i < height; i++)
3293                         for (unsigned int j = 0; j < width; j++) {
3294
3295                                 ui::get()->d3_render_status(1, 0, f, -1, i, j, -1);
3296
3297                                 /*
3298                                  * Check visibility.
3299                                  */
3300
3301                                 int n = i * width + j;
3302
3303                                 if (!visibility_search(fmv, mv[n]))
3304                                         continue;
3305
3306                                 /*
3307                                  * Find depth and diff at this point, check for
3308                                  * undefined values, and generate projections
3309                                  * of the image corners on the estimated normal
3310                                  * surface.
3311                                  */
3312
3313                                 d2::pixel depth = median_depths->pix(i, j);
3314                                 d2::pixel diff = median_diffs->pix(i, j);
3315                                 // d2::pixel diff = d2::pixel(0, 0, 0);
3316
3317                                 if (!depth.finite() || !diff.finite())
3318                                         continue;
3319
3320                                 point local_points[3] = {
3321                                         point(i,     j,     depth[0]),
3322                                         point(i + 1, j,     depth[0] + diff[0]),
3323                                         point(i    , j + 1, depth[0] + diff[1])
3324                                 };
3325
3326                                 /*
3327                                  * Determine transformation at (i, j).  First
3328                                  * determine transformation from the output to
3329                                  * the input, then invert this, as we need the
3330                                  * inverse transformation for filtering.
3331                                  */
3332
3333                                 d2::point remote_points[3] = {
3334                                         _ptf.wp_unscaled(_pt.pw_scaled(point(local_points[0]))).xy(),
3335                                         _ptf.wp_unscaled(_pt.pw_scaled(point(local_points[1]))).xy(),
3336                                         _ptf.wp_unscaled(_pt.pw_scaled(point(local_points[2]))).xy()
3337                                 };
3338
3339                                 /*
3340                                  * Forward matrix for the linear component of the
3341                                  * transformation.
3342                                  */
3343
3344                                 d2::point forward_matrix[2] = {
3345                                         remote_points[1] - remote_points[0],
3346                                         remote_points[2] - remote_points[0]
3347                                 };
3348
3349                                 /*
3350                                  * Inverse matrix for the linear component of
3351                                  * the transformation.  Calculate using the
3352                                  * determinant D.
3353                                  */
3354
3355                                 ale_pos D = forward_matrix[0][0] * forward_matrix[1][1]
3356                                           - forward_matrix[0][1] * forward_matrix[1][0];
3357
3358                                 if (D == 0)
3359                                         continue;
3360
3361                                 d2::point inverse_matrix[2] = {
3362                                         d2::point( forward_matrix[1][1] / D, -forward_matrix[1][0] / D),
3363                                         d2::point(-forward_matrix[0][1] / D,  forward_matrix[0][0] / D)
3364                                 };
3365
3366                                 /*
3367                                  * Determine the projective transformation parameters for the
3368                                  * inverse transformation.
3369                                  */
3370
3371                                 const d2::image *imf = d2::image_rw::open(f);
3372
3373                                 d2::transformation inv_t = d2::transformation::gpt_identity(imf, 1);
3374
3375                                 d2::point local_bounds[4];
3376
3377                                 for (int n = 0; n < 4; n++) {
3378                                         d2::point remote_bound = d2::point((n == 1 || n == 2) ? imf->height() : 0,
3379                                                                            (n == 2 || n == 3) ? imf->width()  : 0)
3380                                                                - remote_points[0];
3381
3382                                         local_bounds[n] = local_points[0].xy()
3383                                                         + d2::point(remote_bound[0] * inverse_matrix[0][0]
3384                                                                   + remote_bound[1] * inverse_matrix[1][0],
3385                                                                     remote_bound[0] * inverse_matrix[0][1]
3386                                                                   + remote_bound[1] * inverse_matrix[1][1]);
3387                                 }
3388
3389                                 inv_t.gpt_set(local_bounds);
3390
3391                                 d2::image_rw::close(f);
3392
3393                                 /*
3394                                  * Perform render step for the given frame,
3395                                  * transformation, and point.
3396                                  */
3397
3398                                 d2::image_rw::open(f);
3399                                 renderer->point_render(i, j, f, inv_t);
3400                                 d2::image_rw::close(f);
3401                         }
3402
3403                         if (tc_multiplier == 0)
3404                                 al->close(f);
3405                 }
3406
3407                 renderer->finish_point_rendering();
3408
3409                 return renderer->get_image();
3410         }
3411
3412         /*
3413          * Generic function.
3414          */
3415         static const d2::image *view(pt _pt, int n = -1) {
3416
3417                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
3418
3419                 if (use_filter) {
3420                         return view_filter(_pt, n);
3421                 } else {
3422                         return view_nofilter(_pt, n);
3423                 }
3424         }
3425
3426         static void tcem(double _tcem) {
3427                 tc_multiplier = _tcem;
3428         }
3429
3430         static void oui(unsigned int _oui) {
3431                 ou_iterations = _oui;
3432         }
3433
3434         static void pa(unsigned int _pa) {
3435                 pairwise_ambiguity = _pa;
3436         }
3437
3438         static void pc(const char *_pc) {
3439                 pairwise_comparisons = _pc;
3440         }
3441
3442         static void d3px(int _d3px_count, double *_d3px_parameters) {
3443                 d3px_count = _d3px_count;
3444                 d3px_parameters = _d3px_parameters;
3445         }
3446
3447         static void fx(double _fx) {
3448                 falloff_exponent = _fx;
3449         }
3450
3451         static void nw() {
3452                 normalize_weights = 1;
3453         }
3454
3455         static void no_nw() {
3456                 normalize_weights = 0;
3457         }
3458
3459         static void nofilter() {
3460                 use_filter = 0;
3461         }
3462
3463         static void filter() {
3464                 use_filter = 1;
3465         }
3466
3467         static void set_filter_type(const char *type) {
3468                 d3chain_type = type;
3469         }
3470
3471         static void set_subspace_traverse() {
3472                 subspace_traverse = 1;
3473         }
3474
3475         static int excluded(point p) {
3476                 for (int n = 0; n < d3px_count; n++) {
3477                         double *region = d3px_parameters + (6 * n);
3478                         if (p[0] >= region[0]
3479                          && p[0] <= region[1]
3480                          && p[1] >= region[2]
3481                          && p[1] <= region[3]
3482                          && p[2] >= region[4]
3483                          && p[2] <= region[5])
3484                                 return 1;
3485                 }
3486
3487                 return 0;
3488         }
3489
3490         /*
3491          * This function returns true if a space is completely excluded.
3492          */
3493         static int excluded(const space::traverse &st) {
3494                 for (int n = 0; n < d3px_count; n++) {
3495                         double *region = d3px_parameters + (6 * n);
3496                         if (st.get_min()[0] >= region[0]
3497                          && st.get_max()[0] <= region[1]
3498                          && st.get_min()[1] >= region[2]
3499                          && st.get_max()[1] <= region[3]
3500                          && st.get_min()[2] >= region[4]
3501                          && st.get_max()[2] <= region[5])
3502                                 return 1;
3503                 }
3504
3505                 return 0;
3506         }
3507
3508         static const d2::image *view(unsigned int n) {
3509
3510                 assert (n < d2::image_rw::count());
3511
3512                 pt _pt = align::projective(n);
3513
3514                 return view(_pt, n);
3515         }
3516
3517         typedef struct {point iw; point ip, is;} analytic;
3518         typedef std::multimap<ale_real,analytic> score_map;
3519         typedef std::pair<ale_real,analytic> score_map_element;
3520
3521         /*
3522          * Make pt list.
3523          */
3524         static std::vector<pt> make_pt_list(const char *d_out[], const char *v_out[],
3525                         std::map<const char *, pt> *d3_depth_pt,
3526                         std::map<const char *, pt> *d3_output_pt) {
3527
3528                 std::vector<pt> result;
3529
3530                 for (unsigned int n = 0; n < d2::image_rw::count(); n++) {
3531                         if (d_out[n] || v_out[n]) {
3532                                 result.push_back(align::projective(n));
3533                         }
3534                 }
3535
3536                 for (std::map<const char *, pt>::iterator i = d3_depth_pt->begin(); i != d3_depth_pt->end(); i++) {
3537                         result.push_back(i->second);
3538                 }
3539
3540                 for (std::map<const char *, pt>::iterator i = d3_output_pt->begin(); i != d3_output_pt->end(); i++) {
3541                         result.push_back(i->second);
3542                 }
3543
3544                 return result;
3545         }
3546
3547         /*
3548          * Get a trilinear coordinate for an anisotropic candidate cell.
3549          */
3550         static ale_pos get_trilinear_coordinate(point min, point max, pt _pt) {
3551
3552                 d2::point local_min, local_max;
3553
3554                 local_min = _pt.wp_unscaled(min).xy();
3555                 local_max = _pt.wp_unscaled(min).xy();
3556
3557                 point cell[2] = {min, max};
3558
3559                 /*
3560                  * Determine the view-local extrema in 2 dimensions.
3561                  */
3562
3563                 for (int r = 1; r < 8; r++) {
3564                         point local = _pt.wp_unscaled(point(cell[r>>2][0], cell[(r>>1)%2][1], cell[r%2][2]));
3565
3566                         for (int d = 0; d < 2; d++) {
3567                                 if (local[d] < local_min[d])
3568                                         local_min[d] = local[d];
3569                                 if (local[d] > local_max[d])
3570                                         local_max[d] = local[d];
3571                                 if (isnan(local[d]))
3572                                         return local[d];
3573                         }
3574                 }
3575
3576                 ale_pos diameter = (local_max - local_min).norm();
3577
3578                 return log(diameter / sqrt(2)) / log(2);
3579         }
3580
3581         /*
3582          * Check whether a cell is visible from a given viewpoint.  This function
3583          * is guaranteed to return 1 when a cell is visible, but it is not guaranteed
3584          * to return 0 when a cell is invisible.
3585          */
3586         static int pt_might_be_visible(const pt &viewpoint, point min, point max) {
3587
3588                 int doc = (rand() % 100000) ? 0 : 1;
3589
3590                 if (doc)
3591                         fprintf(stderr, "checking visibility:\n");
3592
3593                 point cell[2] = {min, max};
3594
3595                 /*
3596                  * Cycle through all vertices of the cell to check certain
3597                  * properties.
3598                  */
3599                 int pos[3] = {0, 0, 0};
3600                 int neg[3] = {0, 0, 0};
3601                 for (int i = 0; i < 2; i++)
3602                 for (int j = 0; j < 2; j++)
3603                 for (int k = 0; k < 2; k++) {
3604                         point p = viewpoint.wp_unscaled(point(cell[i][0], cell[j][1], cell[k][2]));
3605
3606                         if (p[2] < 0 && viewpoint.unscaled_in_bounds(p))
3607                                 return 1;
3608
3609                         if (isnan(p[0])
3610                          || isnan(p[1])
3611                          || isnan(p[2]))
3612                                 return 1;
3613
3614                         if (p[2] > 0)
3615                                 for (int d = 0; d < 2; d++)
3616                                         p[d] *= -1;
3617
3618                         if (doc)
3619                                 fprintf(stderr, "\t[%f %f %f] --> [%f %f %f]\n",
3620                                                 cell[i][0], cell[j][1], cell[k][2],
3621                                                 p[0], p[1], p[2]);
3622
3623                         for (int d = 0; d < 3; d++)
3624                                 if (p[d] >= 0)
3625                                         pos[d] = 1;
3626
3627                         if (p[0] <= viewpoint.unscaled_height() - 1)
3628                                 neg[0] = 1;
3629
3630                         if (p[1] <= viewpoint.unscaled_width() - 1)
3631                                 neg[1] = 1;
3632
3633                         if (p[2] <= 0)
3634                                 neg[2] = 1;
3635                 }
3636
3637                 if (!neg[2])
3638                         return 0;
3639
3640                 if (!pos[0]
3641                  || !neg[0]
3642                  || !pos[1]
3643                  || !neg[1])
3644                         return 0;
3645
3646                 return 1;
3647         }
3648
3649         /*
3650          * Check whether a cell is output-visible.
3651          */
3652         static int output_might_be_visible(const std::vector<pt> &pt_outputs, point min, point max) {
3653                 for (unsigned int n = 0; n < pt_outputs.size(); n++)
3654                         if (pt_might_be_visible(pt_outputs[n], min, max))
3655                                 return 1;
3656                 return 0;
3657         }
3658
3659         /*
3660          * Check whether a cell is input-visible.
3661          */
3662         static int input_might_be_visible(unsigned int f, point min, point max) {
3663                 return pt_might_be_visible(align::projective(f), min, max);
3664         }
3665
3666         /*
3667          * Return true if a cell fails an output resolution bound.
3668          */
3669         static int fails_output_resolution_bound(point min, point max, const std::vector<pt> &pt_outputs) {
3670                 for (unsigned int n = 0; n < pt_outputs.size(); n++) {
3671
3672                         point p = pt_outputs[n].centroid(min, max);
3673
3674                         if (!p.defined())
3675                                 continue;
3676
3677                         if (get_trilinear_coordinate(min, max, pt_outputs[n]) < output_decimation_preferred)
3678                                 return 1;
3679                 }
3680
3681                 return 0;
3682         }
3683
3684         /*
3685          * Check lower-bound resolution constraints
3686          */
3687         static int exceeds_resolution_lower_bounds(unsigned int f1, unsigned int f2,
3688                         point min, point max, const std::vector<pt> &pt_outputs) {
3689
3690                 pt _pt = al->get(f1)->get_t(0);
3691
3692                 if (get_trilinear_coordinate(min, max, _pt) < input_decimation_lower)
3693                         return 1;
3694
3695                 if (fails_output_resolution_bound(min, max, pt_outputs))
3696                         return 0;
3697
3698                 if (get_trilinear_coordinate(min, max, _pt) < primary_decimation_upper)
3699                         return 1;
3700
3701                 return 0;
3702         }
3703
3704         /*
3705          * Try the candidate nearest to the specified cell.
3706          */
3707         static void try_nearest_candidate(unsigned int f1, unsigned int f2, candidates *c, point min, point max) {
3708                 point centroid = (max + min) / 2;
3709                 pt _pt[2] = { al->get(f1)->get_t(0), al->get(f2)->get_t(0) };
3710                 point p[2];
3711
3712                 // fprintf(stderr, "[tnc n=%f %f %f x=%f %f %f]\n", min[0], min[1], min[2], max[0], max[1], max[2]);
3713
3714                 /*
3715                  * Reject clipping plane violations.
3716                  */
3717
3718                 if (centroid[2] > front_clip
3719                  || centroid[2] < rear_clip)
3720                         return;
3721
3722                 /*
3723                  * Calculate projections.
3724                  */
3725
3726                 for (int n = 0; n < 2; n++) {
3727
3728                         p[n] = _pt[n].wp_unscaled(centroid);
3729
3730                         if (!_pt[n].unscaled_in_bounds(p[n]))
3731                                 return;
3732
3733                         // fprintf(stderr, ":");
3734
3735                         if (p[n][2] >= 0)
3736                                 return;
3737                 }
3738
3739
3740                 int tc = (int) round(get_trilinear_coordinate(min, max, _pt[0]));
3741                 int stc = (int) round(get_trilinear_coordinate(min, max, _pt[1]));
3742
3743                 while (tc < input_decimation_lower || stc < input_decimation_lower) {
3744                         tc++;
3745                         stc++;
3746                 }
3747
3748                 if (tc > primary_decimation_upper)
3749                         return;
3750
3751                 /*
3752                  * Calculate score from color match.  Assume for now
3753                  * that the transformation can be approximated locally
3754                  * with a translation.
3755                  */
3756
3757                 ale_pos score = 0;
3758                 ale_pos divisor = 0;
3759                 ale_pos l1_multiplier = 0.125;
3760                 lod_image *if1 = al->get(f1);
3761                 lod_image *if2 = al->get(f2);
3762
3763                 if (if1->in_bounds(p[0].xy())
3764                  && if2->in_bounds(p[1].xy())) {
3765                         divisor += 1 - l1_multiplier;
3766                         score += (1 - l1_multiplier)
3767                                * (if1->get_tl(p[0].xy(), tc) - if2->get_tl(p[1].xy(), stc)).normsq();
3768                 }
3769
3770                 for (int iii = -1; iii <= 1; iii++)
3771                 for (int jjj = -1; jjj <= 1; jjj++) {
3772                         d2::point t(iii, jjj);
3773
3774                         if (!if1->in_bounds(p[0].xy() + t)
3775                          || !if2->in_bounds(p[1].xy() + t))
3776                                 continue;
3777
3778                         divisor += l1_multiplier;
3779                         score   += l1_multiplier
3780                                  * (if1->get_tl(p[0].xy() + t, tc) - if2->get_tl(p[1].xy() + t, tc)).normsq();
3781
3782                 }
3783
3784                 /*
3785                  * Include third-camera contributions in the score.
3786                  */
3787
3788                 if (tc_multiplier != 0)
3789                 for (unsigned int n = 0; n < d2::image_rw::count(); n++) {
3790                         if (n == f1 || n == f2)
3791                                 continue;
3792
3793                         lod_image *ifn = al->get(n);
3794                         pt _ptn = ifn->get_t(0);
3795                         point pn = _ptn.wp_unscaled(centroid);
3796
3797                         if (!_ptn.unscaled_in_bounds(pn))
3798                                 continue;
3799
3800                         if (pn[2] >= 0)
3801                                 continue;
3802
3803                         ale_pos ttc = get_trilinear_coordinate(min, max, _ptn);
3804
3805                         divisor += tc_multiplier;
3806                         score   += tc_multiplier
3807                                  * (if1->get_tl(p[0].xy(), tc) - ifn->get_tl(pn.xy(), ttc)).normsq();
3808                 }
3809
3810                 c->add_candidate(p[0], tc, score / divisor);
3811         }
3812
3813         /*
3814          * Check for cells that are completely clipped.
3815          */
3816         static int completely_clipped(point min, point max) {
3817                 return (min[2] > front_clip
3818                      || max[2] < rear_clip);
3819         }
3820
3821         /*
3822          * Update extremum variables for cell points mapped to a particular view.
3823          */
3824         static void update_extrema(point min, point max, pt _pt, int *extreme_dim, ale_pos *extreme_ratio) {
3825
3826                 point local_min, local_max;
3827
3828                 local_min = _pt.wp_unscaled(min);
3829                 local_max = _pt.wp_unscaled(min);
3830
3831                 point cell[2] = {min, max};
3832
3833                 int near_vertex = 0;
3834
3835                 /*
3836                  * Determine the view-local extrema in all dimensions, and
3837                  * determine the vertex of closest z coordinate.
3838                  */
3839
3840                 for (int r = 1; r < 8; r++) {
3841                         point local = _pt.wp_unscaled(point(cell[r>>2][0], cell[(r>>1)%2][1], cell[r%2][2]));
3842
3843                         for (int d = 0; d < 3; d++) {
3844                                 if (local[d] < local_min[d])
3845                                         local_min[d] = local[d];
3846                                 if (local[d] > local_max[d])
3847                                         local_max[d] = local[d];
3848                         }
3849
3850                         if (local[2] == local_max[2])
3851                                 near_vertex = r;
3852                 }
3853
3854                 ale_pos diameter = (local_max.xy() - local_min.xy()).norm();
3855
3856                 /*
3857                  * Update extrema as necessary for each dimension.
3858                  */
3859
3860                 for (int d = 0; d < 3; d++) {
3861
3862                         int r = near_vertex;
3863
3864                         int p1[3] = {r>>2, (r>>1)%2, r%2};
3865                         int p2[3] = {r>>2, (r>>1)%2, r%2};
3866
3867                         p2[d] = 1 - p2[d];
3868
3869                         ale_pos local_distance = (_pt.wp_unscaled(point(cell[p1[0]][0], cell[p1[1]][1], cell[p1[2]][2])).xy()
3870                                                 - _pt.wp_unscaled(point(cell[p2[0]][0], cell[p2[1]][1], cell[p2[2]][2])).xy()).norm();
3871
3872                         if (local_distance / diameter > *extreme_ratio) {
3873                                 *extreme_ratio = local_distance / diameter;
3874                                 *extreme_dim = d;
3875                         }
3876                 }
3877         }
3878
3879         /*
3880          * Get the next split dimension.
3881          */
3882         static int get_next_split(int f1, int f2, point min, point max, const std::vector<pt> &pt_outputs) {
3883                 for (int d = 0; d < 3; d++)
3884                         if (isinf(min[d]) || isinf(max[d]))
3885                                 return space::traverse::get_next_split(min, max);
3886
3887                 int extreme_dim = 0;
3888                 ale_pos extreme_ratio = 0;
3889
3890                 update_extrema(min, max, al->get(f1)->get_t(0), &extreme_dim, &extreme_ratio);
3891                 update_extrema(min, max, al->get(f2)->get_t(0), &extreme_dim, &extreme_ratio);
3892
3893                 for (unsigned int n = 0; n < pt_outputs.size(); n++) {
3894                         update_extrema(min, max, pt_outputs[n], &extreme_dim, &extreme_ratio);
3895                 }
3896
3897                 return extreme_dim;
3898         }
3899
3900         /*
3901          * Find candidates for subspace creation.
3902          */
3903         static void find_candidates(unsigned int f1, unsigned int f2, candidates *c, point min, point max,
3904                         const std::vector<pt> &pt_outputs, int depth = 0) {
3905
3906                 int print = 0;
3907
3908                 if (min[0] < 20.0001 && max[0] > 20.0001
3909                  && min[1] < 20.0001 && max[1] > 20.0001
3910                  && min[2] < 0.0001 && max[2] > 0.0001)
3911                         print = 1;
3912
3913                 if (print) {
3914                         for (int i = depth; i > 0; i--) {
3915                                 fprintf(stderr, "+");
3916                         }
3917                         fprintf(stderr, "[fc n=%f %f %f x=%f %f %f]\n",
3918                                         min[0], min[1], min[2], max[0], max[1], max[2]);
3919                 }
3920
3921                 if (completely_clipped(min, max)) {
3922                         if (print)
3923                                 fprintf(stderr, "c");
3924                         return;
3925                 }
3926
3927                 if (!input_might_be_visible(f1, min, max)
3928                  || !input_might_be_visible(f2, min, max)) {
3929                         if (print)
3930                                 fprintf(stderr, "v");
3931                         return;
3932                 }
3933
3934                 if (output_clip && !output_might_be_visible(pt_outputs, min, max)) {
3935                         if (print)
3936                                 fprintf(stderr, "o");
3937                         return;
3938                 }
3939
3940                 if (exceeds_resolution_lower_bounds(f1, f2, min, max, pt_outputs)) {
3941                         if (!(rand() % 100000))
3942                                 fprintf(stderr, "([%f %f %f], [%f %f %f]) at %d\n",
3943                                                 min[0], min[1], min[2],
3944                                                 max[0], max[1], max[2],
3945                                                 __LINE__);
3946
3947                         if (print)
3948                                 fprintf(stderr, "t");
3949
3950                         try_nearest_candidate(f1, f2, c, min, max);
3951                         return;
3952                 }
3953
3954                 point new_cells[2][2];
3955
3956                 if (!space::traverse::get_next_cells(get_next_split(f1, f2, min, max, pt_outputs), min, max, new_cells)) {
3957                         if (print)
3958                                 fprintf(stderr, "n");
3959                         return;
3960                 }
3961
3962                 if (print) {
3963                         fprintf(stderr, "nc[0][0]=%f %f %f nc[0][1]=%f %f %f nc[1][0]=%f %f %f nc[1][1]=%f %f %f\n",
3964                                         new_cells[0][0][0],
3965                                         new_cells[0][0][1],
3966                                         new_cells[0][0][2],
3967                                         new_cells[0][1][0],
3968                                         new_cells[0][1][1],
3969                                         new_cells[0][1][2],
3970                                         new_cells[1][0][0],
3971                                         new_cells[1][0][1],
3972                                         new_cells[1][0][2],
3973                                         new_cells[1][1][0],
3974                                         new_cells[1][1][1],
3975                                         new_cells[1][1][2]);
3976                 }
3977
3978                 find_candidates(f1, f2, c, new_cells[0][0], new_cells[0][1], pt_outputs, depth + 1);
3979                 find_candidates(f1, f2, c, new_cells[1][0], new_cells[1][1], pt_outputs, depth + 1);
3980         }
3981
3982         /*
3983          * Generate a map from scores to 3D points for various depths at point (i, j) in f1, at
3984          * lowest resolution.
3985          */
3986         static score_map p2f_score_map(unsigned int f1, unsigned int f2, unsigned int i, unsigned int j) {
3987
3988                 score_map result;
3989
3990                 pt _pt1 = al->get(f1)->get_t(primary_decimation_upper);
3991                 pt _pt2 = al->get(f2)->get_t(primary_decimation_upper);
3992
3993                 const d2::image *if1 = al->get(f1)->get_image(primary_decimation_upper);
3994                 const d2::image *if2 = al->get(f2)->get_image(primary_decimation_upper);
3995                 ale_pos pdu_scale = pow(2, primary_decimation_upper);
3996
3997                 /*
3998                  * Get the pixel color in the primary frame
3999                  */
4000
4001                 // d2::pixel color_primary = if1->get_pixel(i, j);
4002
4003                 /*
4004                  * Map two depths to the secondary frame.
4005                  */
4006
4007                 point p1 = _pt2.wp_unscaled(_pt1.pw_unscaled(point(i, j, 1000)));
4008                 point p2 = _pt2.wp_unscaled(_pt1.pw_unscaled(point(i, j, -1000)));
4009
4010 //              fprintf(stderr, "%d->%d (%d, %d) point pair: (%d, %d, %d -> %f, %f), (%d, %d, %d -> %f, %f)\n",
4011 //                              f1, f2, i, j, i, j, 1000, p1[0], p1[1], i, j, -1000, p2[0], p2[1]);
4012 //              _pt1.debug_output();
4013 //              _pt2.debug_output();
4014
4015
4016                 /*
4017                  * For cases where the mapped points define a
4018                  * line and where points on the line fall
4019                  * within the defined area of the frame,
4020                  * determine the starting point for inspection.
4021                  * In other cases, continue to the next pixel.
4022                  */
4023
4024                 ale_pos diff_i = p2[0] - p1[0];
4025                 ale_pos diff_j = p2[1] - p1[1];
4026                 ale_pos slope = diff_j / diff_i;
4027
4028                 if (isnan(slope)) {
4029                         assert(0);
4030                         fprintf(stderr, "%d->%d (%d, %d) has undefined slope\n",
4031                                         f1, f2, i, j);
4032                         return result;
4033                 }
4034
4035                 /*
4036                  * Make absurdly large/small slopes either infinity, negative infinity, or zero.
4037                  */
4038
4039                 if (fabs(slope) > if2->width() * 100) {
4040                         double zero = 0;
4041                         double one  = 1;
4042                         double inf  = one / zero;
4043                         slope = inf;
4044                 } else if (slope < 1 / (double) if2->height() / 100
4045                         && slope > -1/ (double) if2->height() / 100) {
4046                         slope = 0;
4047                 }
4048
4049                 // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4050
4051                 ale_pos top_intersect = p1[1] - p1[0] * slope;
4052                 ale_pos lef_intersect = p1[0] - p1[1] / slope;
4053                 ale_pos rig_intersect = p1[0] - (p1[1] - if2->width() + 2) / slope;
4054                 ale_pos sp_i, sp_j;
4055
4056                 // fprintf(stderr, "slope == %f\n", slope);
4057
4058
4059                 if (slope == 0) {
4060                         // fprintf(stderr, "case 0\n");
4061                         sp_i = lef_intersect;
4062                         sp_j = 0;
4063                 } else if (finite(slope) && top_intersect >= 0 && top_intersect < if2->width() - 1) {
4064                         // fprintf(stderr, "case 1\n");
4065                         sp_i = 0;
4066                         sp_j = top_intersect;
4067                 } else if (slope > 0 && lef_intersect >= 0 && lef_intersect <= if2->height() - 1) {
4068                         // fprintf(stderr, "case 2\n");
4069                         sp_i = lef_intersect;
4070                         sp_j = 0;
4071                 } else if (slope < 0 && rig_intersect >= 0 && rig_intersect <= if2->height() - 1) {
4072                         // fprintf(stderr, "case 3\n");
4073                         sp_i = rig_intersect;
4074                         sp_j = if2->width() - 2;
4075                 } else {
4076                         // fprintf(stderr, "case 4\n");
4077                         // fprintf(stderr, "%d->%d (%d, %d) does not intersect the defined area\n",
4078                         //              f1, f2, i, j);
4079                         return result;
4080                 }
4081
4082
4083                 // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4084
4085                 /*
4086                  * Determine increment values for examining
4087                  * point, ensuring that incr_i is always
4088                  * positive.
4089                  */
4090
4091                 ale_pos incr_i, incr_j;
4092
4093                 if (fabs(diff_i) > fabs(diff_j)) {
4094                         incr_i = 1;
4095                         incr_j = slope;
4096                 } else if (slope > 0) {
4097                         incr_i = 1 / slope;
4098                         incr_j = 1;
4099                 } else {
4100                         incr_i = -1 / slope;
4101                         incr_j = -1;
4102                 }
4103
4104                 // fprintf(stderr, "%d->%d (%d, %d) increments are (%f, %f)\n",
4105                 //              f1, f2, i, j, incr_i, incr_j);
4106
4107                 /*
4108                  * Examine regions near the projected line.
4109                  */
4110
4111                 for (ale_pos ii = sp_i, jj = sp_j;
4112                         ii <= if2->height() - 1 && jj <= if2->width() - 1 && ii >= 0 && jj >= 0;
4113                         ii += incr_i, jj += incr_j) {
4114
4115                         // fprintf(stderr, "%d->%d (%d, %d) checking (%f, %f)\n",
4116                         //              f1, f2, i, j, ii, jj);
4117
4118 #if 0
4119                         /*
4120                          * Check for higher, lower, and nearby points.
4121                          *
4122                          *      Red   = 2^0
4123                          *      Green = 2^1
4124                          *      Blue  = 2^2
4125                          */
4126
4127                         int higher = 0, lower = 0, nearby = 0;
4128
4129                         for (int iii = 0; iii < 2; iii++)
4130                         for (int jjj = 0; jjj < 2; jjj++) {
4131                                 d2::pixel p = if2->get_pixel((int) floor(ii) + iii, (int) floor(jj) + jjj);
4132
4133                                 for (int k = 0; k < 3; k++) {
4134                                         int bitmask = (int) pow(2, k);
4135
4136                                         if (p[k] > color_primary[k])
4137                                                 higher |= bitmask;
4138                                         if (p[k] < color_primary[k])
4139                                                 lower  |= bitmask;
4140                                         if (fabs(p[k] - color_primary[k]) < nearness)
4141                                                 nearby |= bitmask;
4142                                 }
4143                         }
4144
4145                         /*
4146                          * If this is not a region of interest,
4147                          * then continue.
4148                          */
4149
4150
4151                         fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4152
4153                         // if (((higher & lower) | nearby) != 0x7)
4154                         //      continue;
4155 #endif
4156                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4157
4158                         // fprintf(stderr, "%d->%d (%d, %d) accepted (%f, %f)\n",
4159                         //              f1, f2, i, j, ii, jj);
4160
4161                         /*
4162                          * Create an orthonormal basis to
4163                          * determine line intersection.
4164                          */
4165
4166                         point bp0 = _pt1.pw_unscaled(point(i, j, 0));
4167                         point bp1 = _pt1.pw_unscaled(point(i, j, 10));
4168                         point bp2 = _pt2.pw_unscaled(point(ii, jj, 0));
4169
4170                         point foo = _pt1.wp_unscaled(bp0);
4171                         // fprintf(stderr, "(%d, %d, 0) transformed to world and back is: (%f, %f, %f)\n",
4172                         //              i, j, foo[0], foo[1], foo[2]);
4173
4174                         foo = _pt1.wp_unscaled(bp1);
4175                         // fprintf(stderr, "(%d, %d, 10)  transformed to world and back is: (%f, %f, %f)\n",
4176                         //              i, j, foo[0], foo[1], foo[2]);
4177
4178                         point b0  = (bp1 - bp0).normalize();
4179                         point b1n = bp2 - bp0;
4180                         point b1  = (b1n - b1n.dproduct(b0) * b0).normalize();
4181                         point b2  = point(0, 0, 0).xproduct(b0, b1).normalize(); // Should already have norm=1
4182
4183
4184                         foo = _pt1.wp_unscaled(bp0 + 30 * b0);
4185
4186                         /*
4187                          * Select a fourth point to define a second line.
4188                          */
4189
4190                         point p3  = _pt2.pw_unscaled(point(ii, jj, 10));
4191
4192                         /*
4193                          * Representation in the new basis.
4194                          */
4195
4196                         d2::point nbp0 = d2::point(bp0.dproduct(b0), bp0.dproduct(b1));
4197                         // d2::point nbp1 = d2::point(bp1.dproduct(b0), bp1.dproduct(b1));
4198                         d2::point nbp2 = d2::point(bp2.dproduct(b0), bp2.dproduct(b1));
4199                         d2::point np3  = d2::point( p3.dproduct(b0),  p3.dproduct(b1));
4200
4201                         /*
4202                          * Determine intersection of line
4203                          * (nbp0, nbp1), which is parallel to
4204                          * b0, with line (nbp2, np3).
4205                          */
4206
4207                         /*
4208                          * XXX: a stronger check would be
4209                          * better here, e.g., involving the
4210                          * ratio (np3[0] - nbp2[0]) / (np3[1] -
4211                          * nbp2[1]).  Also, acceptance of these
4212                          * cases is probably better than
4213                          * rejection.
4214                          */
4215
4216
4217                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4218
4219                         if (np3[1] - nbp2[1] == 0)
4220                                 continue;
4221
4222
4223                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4224
4225                         d2::point intersection = d2::point(nbp2[0]
4226                                         + (nbp0[1] - nbp2[1]) * (np3[0] - nbp2[0]) / (np3[1] - nbp2[1]),
4227                                         nbp0[1]);
4228
4229                         ale_pos b2_offset = b2.dproduct(bp0);
4230
4231                         /*
4232                          * Map the intersection back to the world
4233                          * basis.
4234                          */
4235
4236                         point iw = intersection[0] * b0 + intersection[1] * b1 + b2_offset * b2;
4237
4238                         /*
4239                          * Reject intersection points behind a
4240                          * camera.
4241                          */
4242
4243                         point icp = _pt1.wc(iw);
4244                         point ics = _pt2.wc(iw);
4245
4246
4247                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4248
4249                         if (icp[2] >= 0 || ics[2] >= 0)
4250                                 continue;
4251
4252
4253                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4254
4255                         /*
4256                          * Reject clipping plane violations.
4257                          */
4258
4259
4260                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4261
4262                         if (iw[2] > front_clip
4263                          || iw[2] < rear_clip)
4264                                 continue;
4265
4266
4267                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4268
4269                         /*
4270                          * Score the point.
4271                          */
4272
4273                         point ip = _pt1.wp_unscaled(iw);
4274
4275                         point is = _pt2.wp_unscaled(iw);
4276
4277                         analytic _a = { iw, ip, is };
4278
4279                         /*
4280                          * Calculate score from color match.  Assume for now
4281                          * that the transformation can be approximated locally
4282                          * with a translation.
4283                          */
4284
4285                         ale_pos score = 0;
4286                         ale_pos divisor = 0;
4287                         ale_pos l1_multiplier = 0.125;
4288
4289                         if (if1->in_bounds(ip.xy())
4290                          && if2->in_bounds(is.xy())
4291                          && !d2::render::is_excluded_f(ip.xy() * pdu_scale, f1)
4292                          && !d2::render::is_excluded_f(is.xy() * pdu_scale, f2)) {
4293                                 divisor += 1 - l1_multiplier;
4294                                 score += (1 - l1_multiplier)
4295                                        * (if1->get_bl(ip.xy()) - if2->get_bl(is.xy())).normsq();
4296                         }
4297
4298
4299                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4300
4301                         for (int iii = -1; iii <= 1; iii++)
4302                         for (int jjj = -1; jjj <= 1; jjj++) {
4303                                 d2::point t(iii, jjj);
4304
4305                                 if (!if1->in_bounds(ip.xy() + t)
4306                                  || !if2->in_bounds(is.xy() + t)
4307                                  || d2::render::is_excluded_f(ip.xy() * pdu_scale, f1)
4308                                  || d2::render::is_excluded_f(is.xy() * pdu_scale, f2))
4309                                         continue;
4310
4311                                 divisor += l1_multiplier;
4312                                 score   += l1_multiplier
4313                                          * (if1->get_bl(ip.xy() + t) - if2->get_bl(is.xy() + t)).normsq();
4314
4315                         }
4316
4317                         /*
4318                          * Include third-camera contributions in the score.
4319                          */
4320
4321                         if (tc_multiplier != 0)
4322                         for (unsigned int f = 0; f < d2::image_rw::count(); f++) {
4323                                 if (f == f1 || f == f2)
4324                                         continue;
4325
4326                                 const d2::image *if3 = al->get(f)->get_image(primary_decimation_upper);
4327                                 pt _pt3 = al->get(f)->get_t(primary_decimation_upper);
4328
4329                                 point p = _pt3.wp_unscaled(iw);
4330
4331                                 if (!if3->in_bounds(p.xy())
4332                                  || !if1->in_bounds(ip.xy())
4333                                  || d2::render::is_excluded_f(p.xy() * pdu_scale, f)
4334                                  || d2::render::is_excluded_f(ip.xy() * pdu_scale, f1))
4335                                         continue;
4336
4337                                 divisor += tc_multiplier;
4338                                 score   += tc_multiplier
4339                                          * (if1->get_bl(ip.xy()) - if3->get_bl(p.xy())).normsq();
4340                         }
4341
4342
4343
4344
4345                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4346
4347                         /*
4348                          * Reject points with undefined score.
4349                          */
4350
4351
4352                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4353
4354                         if (!finite(score / divisor))
4355                                 continue;
4356
4357
4358                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4359
4360 #if 0
4361                         /*
4362                          * XXX: reject points not on the z=-27.882252 plane.
4363                          */
4364
4365
4366                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4367
4368                         if (_a.ip[2] > -27 || _a.ip[2] < -28)
4369                                 continue;
4370 #endif
4371
4372
4373                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4374
4375                         /*
4376                          * Add the point to the score map.
4377                          */
4378
4379 //                      d2::pixel c_ip = if1->in_bounds(ip.xy()) ? if1->get_bl(ip.xy())
4380 //                                                               : d2::pixel();
4381 //                      d2::pixel c_is = if2->in_bounds(is.xy()) ? if2->get_bl(is.xy())
4382 //                                                               : d2::pixel();
4383
4384 //                      fprintf(stderr, "Candidate subspace: f1=%u f2=%u i=%u j=%u ii=%f jj=%f"
4385 //                                      "cp=[%f %f %f] cs=[%f %f %f]\n",
4386 //                                      f1, f2, i, j, ii, jj, c_ip[0], c_ip[1], c_ip[2],
4387 //                                      c_is[0], c_is[1], c_is[2]);
4388
4389                         result.insert(score_map_element(score / divisor, _a));
4390                 }
4391
4392 //              fprintf(stderr, "Iterating through the score map:\n");
4393 //
4394 //              for (score_map::iterator smi = result.begin(); smi != result.end(); smi++) {
4395 //                      fprintf(stderr, "%f ", smi->first);
4396 //              }
4397 //
4398 //              fprintf(stderr, "\n");
4399
4400                 return result;
4401         }
4402
4403
4404         /*
4405          * Attempt to refine space around a point, to high and low resolutions
4406          * resulting in two resolutions in total.
4407          */
4408
4409         static space::traverse refine_space(point iw, ale_pos target_dim, int use_filler) {
4410
4411                 space::traverse st = space::traverse::root();
4412
4413                 if (!st.includes(iw)) {
4414                         assert(0);
4415                         return st;
4416                 }
4417
4418                 int lr_done = !use_filler;
4419
4420                 /*
4421                  * Loop until all resolutions of interest have been generated.
4422                  */
4423
4424                 for(;;) {
4425
4426                         point p[2] = { st.get_min(), st.get_max() };
4427
4428                         ale_pos dim_max = 0;
4429
4430                         for (int d = 0; d < 3; d++) {
4431                                 ale_pos d_value = fabs(p[0][d] - p[1][d]);
4432                                 if (d_value > dim_max)
4433                                         dim_max = d_value;
4434                         }
4435
4436                         /*
4437                          * Generate any new desired spatial registers.
4438                          */
4439
4440                         for (int f = 0; f < 2; f++) {
4441
4442                                 /*
4443                                  * Low resolution
4444                                  */
4445
4446                                 if (dim_max < 2 * target_dim
4447                                  && lr_done == 0) {
4448                                         if (spatial_info_map.find(st.get_node()) == spatial_info_map.end()) {
4449                                                 spatial_info_map[st.get_node()];
4450                                                 ui::get()->d3_increment_spaces();
4451                                         }
4452                                         lr_done = 1;
4453                                 }
4454
4455                                 /*
4456                                  * High resolution.
4457                                  */
4458
4459                                 if (dim_max < target_dim) {
4460                                         if (spatial_info_map.find(st.get_node()) == spatial_info_map.end()) {
4461                                                 spatial_info_map[st.get_node()];
4462                                                 ui::get()->d3_increment_spaces();
4463                                         }
4464                                         return st;
4465                                 }
4466                         }
4467
4468                         /*
4469                          * Check precision before analyzing space further.
4470                          */
4471
4472                         if (st.precision_wall()) {
4473                                 fprintf(stderr, "\n\n*** Error: reached subspace precision wall ***\n\n");
4474                                 assert(0);
4475                                 return st;
4476                         }
4477
4478                         if (st.positive().includes(iw)) {
4479                                 st = st.positive();
4480                                 total_tsteps++;
4481                         } else if (st.negative().includes(iw)) {
4482                                 st = st.negative();
4483                                 total_tsteps++;
4484                         } else {
4485                                 fprintf(stderr, "failed iw = (%f, %f, %f)\n",
4486                                                 iw[0], iw[1], iw[2]);
4487                                 assert(0);
4488                         }
4489                 }
4490         }
4491
4492         /*
4493          * Calculate target dimension
4494          */
4495
4496         static ale_pos calc_target_dim(point iw, pt _pt, const char *d_out[], const char *v_out[],
4497                         std::map<const char *, pt> *d3_depth_pt,
4498                         std::map<const char *, pt> *d3_output_pt) {
4499
4500                 ale_pos result = _pt.distance_1d(iw, primary_decimation_upper);
4501
4502                 for (unsigned int n = 0; n < d2::image_rw::count(); n++) {
4503                         if (d_out[n] && align::projective(n).distance_1d(iw, 0) < result)
4504                                 result = align::projective(n).distance_1d(iw, 0);
4505                         if (v_out[n] && align::projective(n).distance_1d(iw, 0) < result)
4506                                 result = align::projective(n).distance_1d(iw, 0);
4507                 }
4508
4509                 for (std::map<const char *, pt>::iterator i = d3_output_pt->begin(); i != d3_output_pt->end(); i++) {
4510                         if (i->second.distance_1d(iw, 0) < result)
4511                                 result = i->second.distance_1d(iw, 0);
4512                 }
4513
4514                 for (std::map<const char *, pt>::iterator i = d3_depth_pt->begin(); i != d3_depth_pt->end(); i++) {
4515                         if (i->second.distance_1d(iw, 0) < result)
4516                                 result = i->second.distance_1d(iw, 0);
4517                 }
4518
4519                 assert (result > 0);
4520
4521                 return result;
4522         }
4523
4524         /*
4525          * Calculate level of detail for a given viewpoint.
4526          */
4527
4528         static int calc_lod(ale_pos depth1, pt _pt, ale_pos target_dim) {
4529                 return (int) round(_pt.trilinear_coordinate(depth1, target_dim * sqrt(2)));
4530         }
4531
4532         /*
4533          * Calculate depth range for a given pair of viewpoints.
4534          */
4535
4536         static ale_pos calc_depth_range(point iw, pt _pt1, pt _pt2) {
4537
4538                 point ip = _pt1.wp_unscaled(iw);
4539
4540                 ale_pos reference_change = fabs(ip[2] / 1000);
4541
4542                 point iw1 = _pt1.pw_scaled(ip + point(0, 0, reference_change));
4543                 point iw2 = _pt1.pw_scaled(ip - point(0, 0, reference_change));
4544
4545                 point is = _pt2.wc(iw);
4546                 point is1 = _pt2.wc(iw1);
4547                 point is2 = _pt2.wc(iw2);
4548
4549                 assert(is[2] < 0);
4550
4551                 ale_pos d1 = (is1.xy() - is.xy()).norm();
4552                 ale_pos d2 = (is2.xy() - is.xy()).norm();
4553
4554                 if (is1[2] < 0 && is2[2] < 0) {
4555
4556                         if (d1 > d2)
4557                                 return reference_change / d1;
4558                         else
4559                                 return reference_change / d2;
4560                 }
4561
4562                 if (is1[2] < 0)
4563                         return reference_change / d1;
4564
4565                 if (is2[2] < 0)
4566                         return reference_change / d2;
4567
4568                 return 0;
4569         }
4570
4571         /*
4572          * Calculate a refined point for a given set of parameters.
4573          */
4574
4575         static point get_refined_point(pt _pt1, pt _pt2, int i, int j,
4576                         int f1, int f2, int lod1, int lod2, ale_pos depth,
4577                         ale_pos depth_range) {
4578
4579                 d2::pixel comparison_color = al->get(f1)->get_image(lod1)->get_pixel(i, j);
4580
4581                 ale_pos best = -1;
4582                 ale_pos best_depth = depth;
4583
4584                 for (ale_pos d = depth - depth_range; d < depth + depth_range; d += depth_range / 10) {
4585
4586                         if (!(d < 0))
4587                                 continue;
4588
4589                         point iw = _pt1.pw_unscaled(point(i, j, d));
4590                         point is = _pt2.wp_unscaled(iw);
4591
4592                         if (!(is[2] < 0))
4593                                 continue;
4594
4595                         if (!al->get(f2)->get_image(lod2)->in_bounds(is.xy()))
4596                                 continue;
4597
4598                         ale_pos error = (comparison_color - al->get(f2)->get_image(lod2)->get_bl(is.xy())).norm();
4599
4600                         if (error < best || best == -1) {
4601                                 best = error;
4602                                 best_depth = d;
4603                         }
4604                 }
4605
4606                 return _pt1.pw_unscaled(point(i, j, best_depth));
4607         }
4608
4609         /*
4610          * Analyze space in a manner dependent on the score map.
4611          */
4612
4613         static void analyze_space_from_map(const char *d_out[], const char *v_out[],
4614                                        std::map<const char *, pt> *d3_depth_pt,
4615                                        std::map<const char *, pt> *d3_output_pt,
4616                                        unsigned int f1, unsigned int f2,
4617                                        unsigned int i, unsigned int j, score_map _sm, int use_filler) {
4618
4619                 int accumulated_ambiguity = 0;
4620                 int max_acc_amb = pairwise_ambiguity;
4621
4622                 pt _pt1 = al->get(f1)->get_t(0);
4623                 pt _pt2 = al->get(f2)->get_t(0);
4624
4625                 if (_pt1.scale_2d() != 1)
4626                         use_filler = 1;
4627
4628                 for(score_map::iterator smi = _sm.begin(); smi != _sm.end(); smi++) {
4629
4630                         point iw = smi->second.iw;
4631
4632                         if (accumulated_ambiguity++ >= max_acc_amb)
4633                                 break;
4634
4635                         total_ambiguity++;
4636
4637                         ale_pos depth1 = _pt1.wc(iw)[2];
4638                         ale_pos depth2 = _pt2.wc(iw)[2];
4639
4640                         ale_pos target_dim = calc_target_dim(iw, _pt1, d_out, v_out, d3_depth_pt, d3_output_pt);
4641
4642                         assert(target_dim > 0);
4643
4644                         int lod1 = calc_lod(depth1, _pt1, target_dim);
4645                         int lod2 = calc_lod(depth2, _pt2, target_dim);
4646
4647                         while (lod1 < input_decimation_lower
4648                             || lod2 < input_decimation_lower) {
4649                                 target_dim *= 2;
4650                                 lod1 = calc_lod(depth1, _pt1, target_dim);
4651                                 lod2 = calc_lod(depth2, _pt2, target_dim);
4652                         }
4653
4654
4655                         if (lod1 >= (int) al->get(f1)->count()
4656                          || lod2 >= (int) al->get(f2)->count())
4657                                 continue;
4658
4659                         int multiplier = (unsigned int) floor(pow(2, primary_decimation_upper - lod1));
4660
4661                         ale_pos depth_range = calc_depth_range(iw, _pt1, _pt2);
4662
4663                         pt _pt1_lod = al->get(f1)->get_t(lod1);
4664                         pt _pt2_lod = al->get(f2)->get_t(lod2);
4665
4666                         int im = i * multiplier;
4667                         int jm = j * multiplier;
4668
4669                         for (int ii = 0; ii < multiplier; ii += 1)
4670                         for (int jj = 0; jj < multiplier; jj += 1) {
4671
4672                                 point refined_point = get_refined_point(_pt1_lod, _pt2_lod, im + ii, jm + jj,
4673                                                 f1, f2, lod1, lod2, depth1, depth_range);
4674
4675                                 /*
4676                                  * Re-evaluate target dimension.
4677                                  */
4678
4679                                 ale_pos target_dim_ =
4680                                         calc_target_dim(refined_point, _pt1, d_out, v_out, d3_depth_pt, d3_output_pt);
4681
4682                                 ale_pos depth1_ = _pt1.wc(refined_point)[2];
4683                                 ale_pos depth2_ = _pt2.wc(refined_point)[2];
4684
4685                                 int lod1_ = calc_lod(depth1_, _pt1, target_dim_);
4686                                 int lod2_ = calc_lod(depth2_, _pt2, target_dim_);
4687
4688                                 while (lod1_ < input_decimation_lower
4689                                     || lod2_ < input_decimation_lower) {
4690                                         target_dim_ *= 2;
4691                                         lod1_ = calc_lod(depth1_, _pt1, target_dim_);
4692                                         lod2_ = calc_lod(depth2_, _pt2, target_dim_);
4693                                 }
4694
4695                                 /*
4696                                  * Attempt to refine space around the intersection point.
4697                                  */
4698
4699                                 space::traverse st =
4700                                         refine_space(refined_point, target_dim_, use_filler || _pt1.scale_2d() != 1);
4701
4702                                 assert(resolution_ok(al->get(f1)->get_t(0), al->get(f1)->get_t(0).trilinear_coordinate(st)));
4703                                 assert(resolution_ok(al->get(f2)->get_t(0), al->get(f2)->get_t(0).trilinear_coordinate(st)));
4704                         }
4705
4706                 }
4707         }
4708
4709
4710         /*
4711          * Initialize space and identify regions of interest for the adaptive
4712          * subspace model.
4713          */
4714         static void make_space(const char *d_out[], const char *v_out[],
4715                         std::map<const char *, pt> *d3_depth_pt,
4716                         std::map<const char *, pt> *d3_output_pt) {
4717
4718                 ui::get()->d3_total_spaces(0);
4719
4720                 /*
4721                  * Variable indicating whether low-resolution filler space
4722                  * is desired to avoid aliased gaps in surfaces.
4723                  */
4724
4725                 int use_filler = d3_depth_pt->size() != 0
4726                                    || d3_output_pt->size() != 0
4727                                    || output_decimation_preferred > 0
4728                                    || input_decimation_lower > 0
4729                                    || !focus::is_trivial()
4730                                    || !strcmp(pairwise_comparisons, "all");
4731
4732                 std::vector<pt> pt_outputs = make_pt_list(d_out, v_out, d3_depth_pt, d3_output_pt);
4733
4734                 /*
4735                  * Initialize root space.
4736                  */
4737
4738                 space::init_root();
4739
4740                 /*
4741                  * Special handling for experimental option 'subspace_traverse'.
4742                  */
4743
4744                 if (subspace_traverse) {
4745                         /*
4746                          * Subdivide space to resolve intensity matches between pairs
4747                          * of frames.
4748                          */
4749
4750                         for (unsigned int f1 = 0; f1 < d2::image_rw::count(); f1++) {
4751
4752                                 if (d3_depth_pt->size() == 0
4753                                  && d3_output_pt->size() == 0
4754                                  && d_out[f1] == NULL
4755                                  && v_out[f1] == NULL)
4756                                         continue;
4757
4758                                 if (tc_multiplier == 0)
4759                                         al->open(f1);
4760
4761                                 for (unsigned int f2 = 0; f2 < d2::image_rw::count(); f2++) {
4762
4763                                         if (f1 == f2)
4764                                                 continue;
4765
4766                                         if (tc_multiplier == 0)
4767                                                 al->open(f2);
4768
4769                                         candidates *c = new candidates(f1);
4770
4771                                         find_candidates(f1, f2, c, point::neginf(), point::posinf(), pt_outputs);
4772
4773
4774
4775                                         c->generate_subspaces();
4776
4777                                         if (tc_multiplier == 0)
4778                                                 al->close(f2);
4779                                 }
4780
4781                                 if (tc_multiplier == 0)
4782                                         al->close(f1);
4783                         }
4784
4785                         return;
4786                 }
4787
4788                 /*
4789                  * Subdivide space to resolve intensity matches between pairs
4790                  * of frames.
4791                  */
4792
4793                 for (unsigned int f1 = 0; f1 < d2::image_rw::count(); f1++)
4794                 for (unsigned int f2 = 0; f2 < d2::image_rw::count(); f2++) {
4795                         if (f1 == f2)
4796                                 continue;
4797
4798                         if (!d_out[f1] && !v_out[f1] && !d3_depth_pt->size()
4799                          && !d3_output_pt->size() && strcmp(pairwise_comparisons, "all"))
4800                                 continue;
4801
4802                         if (tc_multiplier == 0) {
4803                                 al->open(f1);
4804                                 al->open(f2);
4805                         }
4806
4807                         /*
4808                          * Iterate over all points in the primary frame.
4809                          */
4810
4811                         ale_pos pdu_scale = pow(2, primary_decimation_upper);
4812
4813                         for (unsigned int i = 0; i < al->get(f1)->get_image(primary_decimation_upper)->height(); i++)
4814                         for (unsigned int j = 0; j < al->get(f1)->get_image(primary_decimation_upper)->width();  j++) {
4815
4816                                 if (d2::render::is_excluded_f(d2::point(i, j) * pdu_scale, f1))
4817                                         continue;
4818
4819                                 ui::get()->d3_subdivision_status(f1, f2, i, j);
4820
4821                                 total_pixels++;
4822
4823                                 /*
4824                                  * Generate a map from scores to 3D points for
4825                                  * various depths in f1.
4826                                  */
4827
4828                                 score_map _sm = p2f_score_map(f1, f2, i, j);
4829
4830                                 /*
4831                                  * Analyze space in a manner dependent on the score map.
4832                                  */
4833
4834                                 analyze_space_from_map(d_out, v_out, d3_depth_pt, d3_output_pt,
4835                                                 f1, f2, i, j, _sm, use_filler);
4836
4837                         }
4838
4839                         /*
4840                          * This ordering may encourage image f1 to be cached.
4841                          */
4842
4843                         if (tc_multiplier == 0) {
4844                                 al->close(f2);
4845                                 al->close(f1);
4846                         }
4847                 }
4848         }
4849
4850
4851         /*
4852          * Update spatial information structures.
4853          *
4854          * XXX: the name of this function is horribly misleading.  There isn't
4855          * even a 'search depth' any longer, since there is no longer any
4856          * bounded DFS occurring.
4857          */
4858         static void reduce_cost_to_search_depth(d2::exposure *exp_out, int inc_bit) {
4859
4860                 /*
4861                  * Subspace model
4862                  */
4863
4864                 ui::get()->set_steps(ou_iterations);
4865
4866                 for (unsigned int i = 0; i < ou_iterations; i++) {
4867                         ui::get()->set_steps_completed(i);
4868                         spatial_info_update();
4869                 }
4870
4871         }
4872
4873 #if 0
4874         /*
4875          * Describe a scene to a renderer
4876          */
4877         static void describe(render *r) {
4878         }
4879 #endif
4880 };
4881
4882 #endif