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                          * Reject out-of-bounds spaces.
1589                          */
1590                         if (!weights->in_spatial_bounds(st)) {
1591                                 si.next();
1592                                 continue;
1593                         }
1594
1595                         /*
1596                          * Skip spaces with no color information.
1597                          *
1598                          * XXX: This could be more efficient, perhaps.
1599                          */
1600
1601                         if (spatial_info_map.count(st.get_node()) == 0) {
1602                                 si.next();
1603                                 continue;
1604                         }
1605
1606                         ui::get()->d3_increment_space_num();
1607
1608
1609                         /*
1610                          * Get in-bounds centroid, if one exists.
1611                          */
1612
1613                         point p = al->get(f)->get_t(0).centroid(st);
1614
1615                         if (!p.defined()) {
1616                                 si.next();
1617                                 continue;
1618                         }
1619
1620                         /*
1621                          * Get information on the subspace.
1622                          */
1623
1624                         spatial_info *sn = &spatial_info_map[st.get_node()];
1625                         d2::pixel color = sn->get_color();
1626                         ale_real occupancy = sn->get_occupancy();
1627
1628                         /*
1629                          * Store current weight so we can later check for
1630                          * modification by higher-resolution subspaces.
1631                          */
1632
1633                         ref_weights::subtree *tree = weights->get_subtree(st);
1634
1635                         /*
1636                          * Check for higher resolution subspaces, and
1637                          * update the space iterator.
1638                          */
1639
1640                         if (st.get_node()->positive
1641                          || st.get_node()->negative) {
1642
1643                                 /*
1644                                  * Cleave space for the higher-resolution pass,
1645                                  * skipping the current space, since we will
1646                                  * process that later.
1647                                  */
1648
1649                                 space::iterate cleaved_space = si.cleave();
1650
1651                                 cleaved_space.next();
1652
1653                                 subspace_info_update(cleaved_space, f, weights);
1654
1655                         } else {
1656                                 si.next();
1657                         }
1658
1659                         /*
1660                          * Add new data on the subspace and update weights.
1661                          */
1662
1663                         ale_pos tc = al->get(f)->get_t(0).trilinear_coordinate(st);
1664                         d2::pixel pcolor = al->get(f)->get_tl(p.xy(), tc);
1665                         d2::pixel colordiff = (color - pcolor) * (ale_real) 256;
1666
1667                         if (falloff_exponent != 0) {
1668                                 d2::pixel max_diff = al->get(f)->get_max_diff(p.xy(), tc) * (ale_real) 256;
1669
1670                                 for (int k = 0; k < 3; k++)
1671                                 if (max_diff[k] > 1)
1672                                         colordiff[k] /= pow(max_diff[k], falloff_exponent);
1673                         }
1674
1675                         /*
1676                          * Determine the probability of encounter.
1677                          */
1678
1679                         d2::pixel encounter = d2::pixel(1, 1, 1) * (1 - weights->get_weight(st));
1680
1681                         /*
1682                          * Update weights
1683                          */
1684
1685                         weights->add_weight(st, occupancy, tree);
1686
1687                         /*
1688                          * Delete the subtree, if necessary.
1689                          */
1690
1691                         delete tree;
1692
1693                         /*
1694                          * Check for cases in which the subspace should not be
1695                          * updated.
1696                          */
1697
1698                         if (!resolution_ok(al->get(f)->get_t(0), tc))
1699                                 continue;
1700
1701                         if (d2::render::is_excluded_f(p.xy(), f))
1702                                 continue;
1703
1704                         /*
1705                          * Update subspace.
1706                          */
1707
1708                         sn->accumulate_color_1(f, pcolor, encounter);
1709                         d2::pixel channel_occ = pexp(-colordiff * colordiff);
1710
1711                         ale_accum occ = channel_occ[0];
1712
1713                         for (int k = 1; k < 3; k++)
1714                                 if (channel_occ[k] < occ)
1715                                         occ = channel_occ[k];
1716
1717                         sn->accumulate_occupancy_1(f, occ, encounter[0]);
1718
1719                 }
1720         }
1721
1722         /*
1723          * Run a single iteration of the spatial_info update cycle.
1724          */
1725         static void spatial_info_update() {
1726                 /*
1727                  * Iterate through each frame.
1728                  */
1729                 for (unsigned int f = 0; f < d2::image_rw::count(); f++) {
1730
1731                         ui::get()->d3_occupancy_status(f);
1732
1733                         /*
1734                          * Open the frame and transformation.
1735                          */
1736
1737                         if (tc_multiplier == 0)
1738                                 al->open(f);
1739
1740                         /*
1741                          * Allocate weights data structure for storing encounter
1742                          * probabilities.
1743                          */
1744
1745                         ref_weights *weights = new ref_weights(f);
1746
1747                         /*
1748                          * Call subspace_info_update for the root space.
1749                          */
1750
1751                         subspace_info_update(space::iterate(al->get(f)->origin()), f, weights);
1752
1753                         /*
1754                          * Free weights.
1755                          */
1756
1757                         delete weights;
1758
1759                         /*
1760                          * Close the frame and transformation.
1761                          */
1762
1763                         if (tc_multiplier == 0)
1764                                 al->close(f);
1765                 }
1766
1767                 /*
1768                  * Update all spatial_info structures.
1769                  */
1770                 for (spatial_info_map_t::iterator i = spatial_info_map.begin(); i != spatial_info_map.end(); i++) {
1771                         i->second.update_color();
1772                         i->second.update_occupancy();
1773
1774 //                      d2::pixel color = i->second.get_color();
1775
1776 //                      fprintf(stderr, "space p=%p updated to c=[%f %f %f] o=%f\n",
1777 //                                      i->first, color[0], color[1], color[2],
1778 //                                      i->second.get_occupancy());
1779                 }
1780         }
1781
1782         /*
1783          * Support function for view() and depth().  This function
1784          * always performs exclusion.
1785          */
1786
1787         static const void view_recurse(int type, d2::image *im, d2::image *weights, space::iterate si, pt _pt,
1788                         int prune = 0, d2::point pl = d2::point(0, 0), d2::point ph = d2::point(0, 0)) {
1789                 while (!si.done()) {
1790                         space::traverse st = si.get();
1791
1792                         /*
1793                          * Remove excluded regions.
1794                          */
1795
1796                         if (excluded(st)) {
1797                                 si.cleave();
1798                                 continue;
1799                         }
1800
1801                         /*
1802                          * Prune.
1803                          */
1804
1805                         if (prune && !_pt.check_inclusion_scaled(st, pl, ph)) {
1806                                 si.cleave();
1807                                 continue;
1808                         }
1809
1810                         /*
1811                          * XXX: This could be more efficient, perhaps.
1812                          */
1813
1814                         if (spatial_info_map.count(st.get_node()) == 0) {
1815                                 si.next();
1816                                 continue;
1817                         }
1818
1819                         ui::get()->d3_increment_space_num();
1820
1821                         spatial_info sn = spatial_info_map[st.get_node()];
1822
1823                         /*
1824                          * Get information on the subspace.
1825                          */
1826
1827                         d2::pixel color = sn.get_color();
1828                         // d2::pixel color = d2::pixel(1, 1, 1) * (double) (((unsigned int) (st.get_node()) / sizeof(space)) % 65535);
1829                         ale_real occupancy = sn.get_occupancy();
1830
1831                         /*
1832                          * Determine the view-local bounding box for the
1833                          * subspace.
1834                          */
1835
1836                         point bb[2];
1837
1838                         _pt.get_view_local_bb_scaled(st, bb);
1839
1840                         point min = bb[0];
1841                         point max = bb[1];
1842
1843                         if (prune) {
1844                                 if (min[0] > ph[0]
1845                                  || min[1] > ph[1]
1846                                  || max[0] < pl[0]
1847                                  || max[1] < pl[1]) {
1848                                         si.next();
1849                                         continue;
1850                                 }
1851
1852                                 if (min[0] < pl[0])
1853                                         min[0] = pl[0];
1854                                 if (min[1] < pl[1])
1855                                         min[1] = pl[1];
1856                                 if (max[0] > ph[0])
1857                                         max[0] = ph[0];
1858                                 if (max[1] > ph[1])
1859                                         max[1] = ph[1];
1860
1861                                 min[0] -= pl[0];
1862                                 min[1] -= pl[1];
1863                                 max[0] -= pl[0];
1864                                 max[1] -= pl[1];
1865                         }
1866
1867                         /*
1868                          * Data structure to check modification of weights by
1869                          * higher-resolution subspaces.
1870                          */
1871
1872                         std::queue<d2::pixel> weight_queue;
1873
1874                         /*
1875                          * Check for higher resolution subspaces, and
1876                          * update the space iterator.
1877                          */
1878
1879                         if (st.get_node()->positive
1880                          || st.get_node()->negative) {
1881
1882                                 /*
1883                                  * Store information about current weights,
1884                                  * so we will know which areas have been
1885                                  * covered by higher-resolution subspaces.
1886                                  */
1887
1888                                 for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
1889                                 for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++)
1890                                         weight_queue.push(weights->get_pixel(i, j));
1891
1892                                 /*
1893                                  * Cleave space for the higher-resolution pass,
1894                                  * skipping the current space, since we will
1895                                  * process that afterward.
1896                                  */
1897
1898                                 space::iterate cleaved_space = si.cleave();
1899
1900                                 cleaved_space.next();
1901
1902                                 view_recurse(type, im, weights, cleaved_space, _pt, prune, pl, ph);
1903
1904                         } else {
1905                                 si.next();
1906                         }
1907
1908
1909                         /*
1910                          * Iterate over pixels in the bounding box, finding
1911                          * pixels that intersect the subspace.  XXX: assume
1912                          * for now that all pixels in the bounding box
1913                          * intersect the subspace.
1914                          */
1915
1916                         for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
1917                         for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++) {
1918
1919                                 /*
1920                                  * Check for higher-resolution updates.
1921                                  */
1922
1923                                 if (weight_queue.size()) {
1924                                         if (weight_queue.front() != weights->get_pixel(i, j)) {
1925                                                 weight_queue.pop();
1926                                                 continue;
1927                                         }
1928                                         weight_queue.pop();
1929                                 }
1930
1931                                 /*
1932                                  * Determine the probability of encounter.
1933                                  */
1934
1935                                 d2::pixel encounter = (d2::pixel(1, 1, 1)
1936                                                      - weights->get_pixel(i, j))
1937                                                     * occupancy;
1938
1939                                 /*
1940                                  * Update images.
1941                                  */
1942
1943                                 if (type == 0) {
1944
1945                                         /*
1946                                          * Color view
1947                                          */
1948
1949                                         weights->pix(i, j) += encounter;
1950                                         im->pix(i, j)      += encounter * color;
1951
1952                                 } else if (type == 1) {
1953
1954                                         /*
1955                                          * Weighted (transparent) depth display
1956                                          */
1957
1958                                         ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
1959                                         weights->pix(i, j) += encounter;
1960                                         im->pix(i, j)      += encounter * depth_value;
1961
1962                                 } else if (type == 2) {
1963
1964                                         /*
1965                                          * Ambiguity (ambivalence) measure.
1966                                          */
1967
1968                                         weights->pix(i, j) = d2::pixel(1, 1, 1);
1969                                         im->pix(i, j) += 0.1 * d2::pixel(1, 1, 1);
1970
1971                                 } else if (type == 3) {
1972
1973                                         /*
1974                                          * Closeness measure.
1975                                          */
1976
1977                                         ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
1978                                         if (weights->pix(i, j)[0] == 0) {
1979                                                 weights->pix(i, j) = d2::pixel(1, 1, 1);
1980                                                 im->pix(i, j) = d2::pixel(1, 1, 1) * depth_value;
1981                                         } else if (im->pix(i, j)[2] < depth_value) {
1982                                                 im->pix(i, j) = d2::pixel(1, 1, 1) * depth_value;
1983                                         } else {
1984                                                 continue;
1985                                         }
1986
1987                                 } else if (type == 4) {
1988
1989                                         /*
1990                                          * Weighted (transparent) contribution display
1991                                          */
1992
1993                                         ale_pos contribution_value = sn.get_pocc_density() /* + sn.get_socc_density() */;
1994                                         weights->pix(i, j) += encounter;
1995                                         im->pix(i, j)      += encounter * contribution_value;
1996
1997                                         assert (finite(encounter[0]));
1998                                         assert (finite(contribution_value));
1999
2000                                 } else if (type == 5) {
2001
2002                                         /*
2003                                          * Weighted (transparent) occupancy display
2004                                          */
2005
2006                                         ale_pos contribution_value = occupancy;
2007                                         weights->pix(i, j) += encounter;
2008                                         im->pix(i, j)      += encounter * contribution_value;
2009
2010                                 } else if (type == 6) {
2011
2012                                         /*
2013                                          * (Depth, xres, yres) triple
2014                                          */
2015
2016                                         ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
2017                                         weights->pix(i, j)[0] += encounter[0];
2018                                         if (weights->pix(i, j)[1] < encounter[0]) {
2019                                                 weights->pix(i, j)[1] = encounter[0];
2020                                                 im->pix(i, j)[0] = weights->pix(i, j)[1] * depth_value;
2021                                                 im->pix(i, j)[1] = max[0] - min[0];
2022                                                 im->pix(i, j)[2] = max[1] - min[1];
2023                                         }
2024
2025                                 } else if (type == 7) {
2026
2027                                         /*
2028                                          * (xoff, yoff, 0) triple
2029                                          */
2030
2031                                         weights->pix(i, j)[0] += encounter[0];
2032                                         if (weights->pix(i, j)[1] < encounter[0]) {
2033                                                 weights->pix(i, j)[1] = encounter[0];
2034                                                 im->pix(i, j)[0] = i - min[0];
2035                                                 im->pix(i, j)[1] = j - min[1];
2036                                                 im->pix(i, j)[2] = 0;
2037                                         }
2038
2039                                 } else if (type == 8) {
2040
2041                                         /*
2042                                          * Value = 1 for any intersected space.
2043                                          */
2044
2045                                         weights->pix(i, j) = d2::pixel(1, 1, 1);
2046                                         im->pix(i, j) = d2::pixel(1, 1, 1);
2047
2048                                 } else if (type == 9) {
2049
2050                                         /*
2051                                          * Number of contributions for the nearest space.
2052                                          */
2053
2054                                         if (weights->pix(i, j)[0] == 1)
2055                                                 continue;
2056
2057                                         weights->pix(i, j) = d2::pixel(1, 1, 1);
2058                                         im->pix(i, j) = d2::pixel(1, 1, 1) * (sn.get_pocc_density() * 0.1);
2059
2060                                 } else
2061                                         assert(0);
2062                         }
2063                 }
2064         }
2065
2066         /*
2067          * Generate an depth image from a specified view.
2068          */
2069         static const d2::image *depth(pt _pt, int n = -1, int prune = 0,
2070                         d2::point pl = d2::point(0, 0), d2::point ph = d2::point(0, 0)) {
2071                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
2072
2073                 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
2074
2075                 if (n >= 0) {
2076                         assert((int) floor(d2::align::of(n).scaled_height())
2077                              == (int) floor(_pt.scaled_height()));
2078                         assert((int) floor(d2::align::of(n).scaled_width())
2079                              == (int) floor(_pt.scaled_width()));
2080                 }
2081
2082                 d2::image *im1, *im2, *im3, *weights;;
2083
2084                 if (prune) {
2085
2086                         im1 = new d2::image_ale_real((int) floor(ph[0] - pl[0]) + 1,
2087                                         (int) floor(ph[1] - pl[1]) + 1, 3);
2088
2089                         im2 = new d2::image_ale_real((int) floor(ph[0] - pl[0]) + 1,
2090                                         (int) floor(ph[1] - pl[1]) + 1, 3);
2091
2092                         im3 = new d2::image_ale_real((int) floor(ph[0] - pl[0]) + 1,
2093                                         (int) floor(ph[1] - pl[1]) + 1, 3);
2094
2095                         weights = new d2::image_ale_real((int) floor(ph[0] - pl[0]) + 1,
2096                                         (int) floor(ph[1] - pl[1]) + 1, 3);
2097
2098                 } else {
2099
2100                         im1 = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2101                                        (int) floor(_pt.scaled_width()), 3);
2102
2103                         im2 = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2104                                        (int) floor(_pt.scaled_width()), 3);
2105
2106                         im3 = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2107                                        (int) floor(_pt.scaled_width()), 3);
2108
2109                         weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2110                                                         (int) floor(_pt.scaled_width()), 3);
2111                 }
2112
2113                 /*
2114                  * Iterate through subspaces.
2115                  */
2116
2117                 space::iterate si(_pt.origin());
2118
2119                 view_recurse(6, im1, weights, si, _pt, prune, pl, ph);
2120
2121                 delete weights;
2122
2123                 if (prune) {
2124                         weights = new d2::image_ale_real((int) floor(ph[0] - pl[0]) + 1,
2125                                         (int) floor(ph[1] - pl[1]) + 1, 3);
2126                 } else {
2127                         weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2128                                                         (int) floor(_pt.scaled_width()), 3);
2129                 }
2130
2131 #if 1
2132                 view_recurse(7, im2, weights, si, _pt, prune, pl, ph);
2133 #else
2134                 view_recurse(8, im2, weights, si, _pt, prune, pl, ph);
2135                 return im2;
2136 #endif
2137
2138                 /*
2139                  * Normalize depths by weights
2140                  */
2141
2142                 if (normalize_weights)
2143                 for (unsigned int i = 0; i < im1->height(); i++)
2144                 for (unsigned int j = 0; j < im1->width();  j++)
2145                         im1->pix(i, j)[0] /= weights->pix(i, j)[1];
2146
2147
2148                 for (unsigned int i = 0; i < im1->height(); i++)
2149                 for (unsigned int j = 0; j < im1->width();  j++) {
2150
2151                         /*
2152                          * Handle interpolation.
2153                          */
2154
2155                         d2::point x;
2156                         d2::point blx;
2157                         d2::point res(im1->pix(i, j)[1], im1->pix(i, j)[2]);
2158
2159                         for (int d = 0; d < 2; d++) {
2160
2161                                 if (im2->pix(i, j)[d] < res[d] / 2)
2162                                         x[d] = (ale_pos) (d?j:i) - res[d] / 2 - im2->pix(i, j)[d];
2163                                 else
2164                                         x[d] = (ale_pos) (d?j:i) + res[d] / 2 - im2->pix(i, j)[d];
2165
2166                                 blx[d] = 1 - ((d?j:i) - x[d]) / res[d];
2167                         }
2168
2169                         ale_real depth_val = 0;
2170                         ale_real depth_weight = 0;
2171
2172                         for (int ii = 0; ii < 2; ii++)
2173                         for (int jj = 0; jj < 2; jj++) {
2174                                 d2::point p = x + d2::point(ii, jj) * res;
2175                                 if (im1->in_bounds(p)) {
2176
2177                                         ale_real d = im1->get_bl(p)[0];
2178
2179                                         if (isnan(d))
2180                                                 continue;
2181
2182                                         ale_real w = ((ii ? (1 - blx[0]) : blx[0]) * (jj ? (1 - blx[1]) : blx[1]));
2183                                         depth_weight += w;
2184                                         depth_val += w * d;
2185                                 }
2186                         }
2187
2188                         ale_real depth = depth_val / depth_weight;
2189
2190                         /*
2191                          * Handle encounter thresholds
2192                          */
2193
2194                         if (weights->pix(i, j)[0] < encounter_threshold) {
2195                                 im3->pix(i, j) = d2::pixel::zero() / d2::pixel::zero();
2196                         } else {
2197                                 im3->pix(i, j) = d2::pixel(1, 1, 1) * depth;
2198                         }
2199                 }
2200
2201                 delete weights;
2202                 delete im1;
2203                 delete im2;
2204
2205                 return im3;
2206         }
2207
2208         static const d2::image *depth(unsigned int n) {
2209
2210                 assert (n < d2::image_rw::count());
2211
2212                 pt _pt = align::projective(n);
2213
2214                 return depth(_pt, n);
2215         }
2216
2217
2218         /*
2219          * This function always performs exclusion.
2220          */
2221
2222         static space::node *most_visible_pointwise(d2::pixel *weight, space::iterate si, pt _pt, d2::point p) {
2223
2224                 space::node *result = NULL;
2225
2226                 while (!si.done()) {
2227                         space::traverse st = si.get();
2228
2229                         /*
2230                          * Prune certain regions known to be uninteresting.
2231                          */
2232
2233                         if (excluded(st) || !_pt.check_inclusion_scaled(st, p)) {
2234                                 si.cleave();
2235                                 continue;
2236                         }
2237
2238                         /*
2239                          * XXX: This could be more efficient, perhaps.
2240                          */
2241
2242                         if (spatial_info_map.count(st.get_node()) == 0) {
2243                                 si.next();
2244                                 continue;
2245                         }
2246
2247                         spatial_info sn = spatial_info_map[st.get_node()];
2248
2249                         /*
2250                          * Get information on the subspace.
2251                          */
2252
2253                         ale_real occupancy = sn.get_occupancy();
2254
2255                         /*
2256                          * Preserve current weight in order to check for
2257                          * modification by higher-resolution subspaces.
2258                          */
2259
2260                         d2::pixel old_weight = *weight;
2261
2262                         /*
2263                          * Check for higher resolution subspaces, and
2264                          * update the space iterator.
2265                          */
2266
2267                         if (st.get_node()->positive
2268                          || st.get_node()->negative) {
2269
2270                                 /*
2271                                  * Cleave space for the higher-resolution pass,
2272                                  * skipping the current space, since we will
2273                                  * process that afterward.
2274                                  */
2275
2276                                 space::iterate cleaved_space = si.cleave();
2277
2278                                 cleaved_space.next();
2279
2280                                 space::node *r = most_visible_pointwise(weight, cleaved_space, _pt, p);
2281
2282                                 if (old_weight[1] != (*weight)[1])
2283                                         result = r;
2284
2285                         } else {
2286                                 si.next();
2287                         }
2288
2289
2290                         /*
2291                          * Check for higher-resolution updates.
2292                          */
2293
2294                         if (old_weight != *weight)
2295                                 continue;
2296
2297                         /*
2298                          * Determine the probability of encounter.
2299                          */
2300
2301                         ale_pos encounter = (1 - (*weight)[0]) * occupancy;
2302
2303                         /*
2304                          * (*weight)[0] stores the cumulative weight; (*weight)[1] stores the maximum.
2305                          */
2306
2307                         if (encounter > (*weight)[1]) {
2308                                 result = st.get_node();
2309                                 (*weight)[1] = encounter;
2310                         }
2311
2312                         (*weight)[0] += encounter;
2313                 }
2314
2315                 return result;
2316         }
2317
2318         /*
2319          * This function performs exclusion iff SCALED is true.
2320          */
2321         static void  most_visible_generic(std::vector<space::node *> &results, d2::image *weights,
2322                         space::iterate si, pt _pt, int scaled) {
2323
2324                 assert (results.size() == weights->height() * weights->width());
2325
2326                 while (!si.done()) {
2327                         space::traverse st = si.get();
2328
2329                         if (scaled && excluded(st)) {
2330                                 si.cleave();
2331                                 continue;
2332                         }
2333
2334                         /*
2335                          * XXX: This could be more efficient, perhaps.
2336                          */
2337
2338                         if (spatial_info_map.count(st.get_node()) == 0) {
2339                                 si.next();
2340                                 continue;
2341                         }
2342
2343                         spatial_info sn = spatial_info_map[st.get_node()];
2344
2345                         /*
2346                          * Get information on the subspace.
2347                          */
2348
2349                         ale_real occupancy = sn.get_occupancy();
2350
2351                         /*
2352                          * Determine the view-local bounding box for the
2353                          * subspace.
2354                          */
2355
2356                         point bb[2];
2357
2358                         _pt.get_view_local_bb_scaled(st, bb);
2359
2360                         point min = bb[0];
2361                         point max = bb[1];
2362
2363                         /*
2364                          * Data structure to check modification of weights by
2365                          * higher-resolution subspaces.
2366                          */
2367
2368                         std::queue<d2::pixel> weight_queue;
2369
2370                         /*
2371                          * Check for higher resolution subspaces, and
2372                          * update the space iterator.
2373                          */
2374
2375                         if (st.get_node()->positive
2376                          || st.get_node()->negative) {
2377
2378                                 /*
2379                                  * Store information about current weights,
2380                                  * so we will know which areas have been
2381                                  * covered by higher-resolution subspaces.
2382                                  */
2383
2384                                 for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
2385                                 for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++)
2386                                         weight_queue.push(weights->get_pixel(i, j));
2387
2388                                 /*
2389                                  * Cleave space for the higher-resolution pass,
2390                                  * skipping the current space, since we will
2391                                  * process that afterward.
2392                                  */
2393
2394                                 space::iterate cleaved_space = si.cleave();
2395
2396                                 cleaved_space.next();
2397
2398                                 most_visible_generic(results, weights, cleaved_space, _pt, scaled);
2399
2400                         } else {
2401                                 si.next();
2402                         }
2403
2404
2405                         /*
2406                          * Iterate over pixels in the bounding box, finding
2407                          * pixels that intersect the subspace.  XXX: assume
2408                          * for now that all pixels in the bounding box
2409                          * intersect the subspace.
2410                          */
2411
2412                         for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
2413                         for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++) {
2414
2415                                 /*
2416                                  * Check for higher-resolution updates.
2417                                  */
2418
2419                                 if (weight_queue.size()) {
2420                                         if (weight_queue.front() != weights->get_pixel(i, j)) {
2421                                                 weight_queue.pop();
2422                                                 continue;
2423                                         }
2424                                         weight_queue.pop();
2425                                 }
2426
2427                                 /*
2428                                  * Determine the probability of encounter.
2429                                  */
2430
2431                                 ale_pos encounter = (1 - weights->get_pixel(i, j)[0]) * occupancy;
2432
2433                                 /*
2434                                  * weights[0] stores the cumulative weight; weights[1] stores the maximum.
2435                                  */
2436
2437                                 if (encounter > weights->get_pixel(i, j)[1]
2438                                  || results[i * weights->width() + j] == NULL) {
2439                                         results[i * weights->width() + j] = st.get_node();
2440                                         weights->chan(i, j, 1) = encounter;
2441                                 }
2442
2443                                 weights->chan(i, j, 0) += encounter;
2444                         }
2445                 }
2446         }
2447
2448         static std::vector<space::node *> most_visible_scaled(pt _pt) {
2449                 d2::image *weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2450                                 (int) floor(_pt.scaled_width()), 3);
2451                 std::vector<space::node *> results;
2452
2453                 results.resize(weights->height() * weights->width(), 0);
2454
2455                 most_visible_generic(results, weights, space::iterate(_pt.origin()), _pt, 1);
2456
2457                 return results;
2458         }
2459
2460         static std::vector<space::node *> most_visible_unscaled(pt _pt) {
2461                 d2::image *weights = new d2::image_ale_real((int) floor(_pt.unscaled_height()),
2462                                 (int) floor(_pt.unscaled_width()), 3);
2463                 std::vector<space::node *> results;
2464
2465                 results.resize(weights->height() * weights->width(), 0);
2466
2467                 most_visible_generic(results, weights, space::iterate(_pt.origin()), _pt, 0);
2468
2469                 return results;
2470         }
2471
2472         static const int visibility_search(const std::vector<space::node *> &fmv, space::node *mv) {
2473
2474                 if (mv == NULL)
2475                         return 0;
2476
2477                 if (std::binary_search(fmv.begin(), fmv.end(), mv))
2478                         return 1;
2479
2480                 return (visibility_search(fmv, mv->positive)
2481                      || visibility_search(fmv, mv->negative));
2482
2483         }
2484
2485         /*
2486          * Class to generate focal sample views.
2487          */
2488
2489         class view_generator {
2490
2491                 /*
2492                  * Original projective transformation.
2493                  */
2494
2495                 pt original_pt;
2496
2497                 /*
2498                  * Data type for shared view data.
2499                  */
2500
2501                 class shared_view {
2502                         pt _pt;
2503                         std::vector<space::node *> mv;
2504                         d2::image *color;
2505                         d2::image *color_weights;
2506                         const d2::image *_depth;
2507                         d2::image *median_depth;
2508                         d2::image *median_diff;
2509
2510                 public:
2511                         shared_view(pt _pt) {
2512                                 this->_pt = _pt;
2513                                 color = NULL;
2514                                 color_weights = NULL;
2515                                 _depth = NULL;
2516                                 median_depth = NULL;
2517                                 median_diff = NULL;
2518                         }
2519
2520                         shared_view(const shared_view &copy_origin) {
2521                                 _pt = copy_origin._pt;
2522                                 mv = copy_origin.mv;
2523                                 color = NULL;
2524                                 color_weights = NULL;
2525                                 _depth = NULL;
2526                                 median_depth = NULL;
2527                                 median_diff = NULL;
2528                         }
2529
2530                         ~shared_view() {
2531                                 delete color;
2532                                 delete _depth;
2533                                 delete color_weights;
2534                                 delete median_diff;
2535                                 delete median_depth;
2536                         }
2537
2538                         void get_view_recurse(d2::image *data, d2::image *weights, int type) {
2539                                 /*
2540                                  * Iterate through subspaces.
2541                                  */
2542
2543                                 space::iterate si(_pt.origin());
2544
2545                                 ui::get()->d3_render_status(0, 0, -1, -1, -1, -1, 0);
2546
2547                                 view_recurse(type, data, weights, si, _pt);
2548                         }
2549
2550                         void init_color() {
2551                                 color = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2552                                                                 (int) floor(_pt.scaled_width()), 3);
2553
2554                                 color_weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2555                                                                 (int) floor(_pt.scaled_width()), 3);
2556
2557                                 get_view_recurse(color, color_weights, 0);
2558                         }
2559
2560                         void init_depth() {
2561                                 _depth = depth(_pt, -1);
2562                         }
2563
2564                         void init_medians() {
2565                                 if (!_depth)
2566                                         init_depth();
2567
2568                                 assert(_depth);
2569
2570                                 median_diff = _depth->fcdiff_median((int) floor(diff_median_radius));
2571                                 median_depth = _depth->medians((int) floor(depth_median_radius));
2572
2573                                 assert(median_diff);
2574                                 assert(median_depth);
2575                         }
2576
2577                 public:
2578                         pt get_pt() {
2579                                 return _pt;
2580                         }
2581
2582                         space::node *get_most_visible(unsigned int i, unsigned int j) {
2583                                 unsigned int height = (int) floor(_pt.scaled_height());
2584                                 unsigned int width  = (int) floor(_pt.scaled_width());
2585
2586                                 if (i < 0 || i >= height
2587                                  || j < 0 || j >= width) {
2588                                         return NULL;
2589                                 }
2590
2591                                 if (mv.size() == 0) {
2592                                         mv = most_visible_scaled(_pt);
2593                                 }
2594
2595                                 assert (mv.size() > i * width + j);
2596
2597                                 return mv[i * width + j];
2598                         }
2599
2600                         space::node *get_most_visible(d2::point p) {
2601                                 unsigned int i = (unsigned int) round (p[0]);
2602                                 unsigned int j = (unsigned int) round (p[1]);
2603
2604                                 return get_most_visible(i, j);
2605                         }
2606
2607                         d2::pixel get_color(unsigned int i, unsigned int j) {
2608                                 if (color == NULL) {
2609                                         init_color();
2610                                 }
2611
2612                                 assert (color != NULL);
2613
2614                                 return color->get_pixel(i, j);
2615                         }
2616
2617                         d2::pixel get_depth(unsigned int i, unsigned int j) {
2618                                 if (_depth == NULL) {
2619                                         init_depth();
2620                                 }
2621
2622                                 assert (_depth != NULL);
2623
2624                                 return _depth->get_pixel(i, j);
2625                         }
2626
2627                         void get_median_depth_and_diff(d2::pixel *t, d2::pixel *f, unsigned int i, unsigned int j) {
2628                                 if (median_depth == NULL && median_diff == NULL)
2629                                         init_medians();
2630
2631                                 assert (median_depth && median_diff);
2632
2633                                 if (i < 0 || i >= median_depth->height()
2634                                  || j < 0 || j >= median_depth->width()) {
2635                                         *t = d2::pixel::undefined();
2636                                         *f = d2::pixel::undefined();
2637                                 } else {
2638                                         *t = median_depth->get_pixel(i, j);
2639                                         *f = median_diff->get_pixel(i, j);
2640                                 }
2641                         }
2642
2643                         void get_color_and_weight(d2::pixel *c, d2::pixel *w, d2::point p) {
2644                                 if (color == NULL) {
2645                                         init_color();
2646                                 }
2647
2648                                 assert (color != NULL);
2649
2650                                 if (!color->in_bounds(p)) {
2651                                         *c = d2::pixel::undefined();
2652                                         *w = d2::pixel::undefined();
2653                                 } else {
2654                                         *c = color->get_bl(p);
2655                                         *w = color_weights->get_bl(p);
2656                                 }
2657                         }
2658
2659                         d2::pixel get_depth(d2::point p) {
2660                                 if (_depth == NULL) {
2661                                         init_depth();
2662                                 }
2663
2664                                 assert (_depth != NULL);
2665
2666                                 if (!_depth->in_bounds(p)) {
2667                                         return d2::pixel::undefined();
2668                                 }
2669
2670                                 return _depth->get_bl(p);
2671                         }
2672
2673                         void get_median_depth_and_diff(d2::pixel *t, d2::pixel *f, d2::point p) {
2674                                 if (median_diff == NULL && median_depth == NULL)
2675                                         init_medians();
2676
2677                                 assert (median_diff != NULL && median_depth != NULL);
2678
2679                                 if (!median_diff->in_bounds(p)) {
2680                                         *t = d2::pixel::undefined();
2681                                         *f = d2::pixel::undefined();
2682                                 } else {
2683                                         *t = median_depth->get_bl(p);
2684                                         *f = median_diff->get_bl(p);
2685                                 }
2686                         }
2687
2688                 };
2689
2690                 /*
2691                  * Shared view array, indexed by aperture diameter and view number.
2692                  */
2693
2694                 std::map<ale_pos, std::vector<shared_view> > aperture_to_shared_views_map;
2695
2696                 /*
2697                  * Method to generate a new stochastic focal view.
2698                  */
2699
2700                 pt get_new_view(ale_pos aperture) {
2701
2702                         ale_pos ofx = aperture;
2703                         ale_pos ofy = aperture;
2704
2705                         while (ofx * ofx + ofy * ofy > aperture * aperture / 4) {
2706                                 ofx = (rand() * aperture) / RAND_MAX - aperture / 2;
2707                                 ofy = (rand() * aperture) / RAND_MAX - aperture / 2;
2708                         }
2709
2710                         /*
2711                          * Generate a new view from the given offset.
2712                          */
2713
2714                         point new_view = original_pt.cw(point(ofx, ofy, 0));
2715                         pt _pt_new = original_pt;
2716                         for (int d = 0; d < 3; d++)
2717                                 _pt_new.e().set_translation(d, -new_view[d]);
2718
2719                         return _pt_new;
2720                 }
2721
2722         public:
2723
2724                 /*
2725                  * Result type.
2726                  */
2727
2728                 class view {
2729                         shared_view *sv;
2730                         pt _pt;
2731
2732                 public:
2733
2734                         view(shared_view *sv, pt _pt = pt()) {
2735                                 this->sv = sv;
2736                                 if (sv) {
2737                                         this->_pt = sv->get_pt();
2738                                 } else {
2739                                         this->_pt = _pt;
2740                                 }
2741                         }
2742
2743                         pt get_pt() {
2744                                 return _pt;
2745                         }
2746
2747                         space::node *get_most_visible(unsigned int i, unsigned int j) {
2748                                 assert (sv);
2749                                 return sv->get_most_visible(i, j);
2750                         }
2751
2752                         space::node *get_most_visible(d2::point p) {
2753                                 if (sv) {
2754                                         return sv->get_most_visible(p);
2755                                 }
2756
2757                                 d2::pixel weight(0, 0, 0);
2758
2759                                 return most_visible_pointwise(&weight, space::iterate(_pt.origin()), _pt, p);
2760
2761                         }
2762
2763                         d2::pixel get_color(unsigned int i, unsigned int j) {
2764                                 return sv->get_color(i, j);
2765                         }
2766
2767                         void get_color_and_weight(d2::pixel *color, d2::pixel *weight, d2::point p) {
2768                                 if (sv) {
2769                                         sv->get_color_and_weight(color, weight, p);
2770                                         return;
2771                                 }
2772
2773                                 /*
2774                                  * Determine weight and color for the given point.
2775                                  */
2776
2777                                 d2::image *im_point = new d2::image_ale_real(1, 1, 3);
2778                                 d2::image *wt_point = new d2::image_ale_real(1, 1, 3);
2779
2780                                 view_recurse(0, im_point, wt_point, space::iterate(_pt.origin()), _pt, 1, p, p);
2781
2782                                 *color = im_point->pix(0, 0);
2783                                 *weight = wt_point->pix(0, 0);
2784
2785                                 delete im_point;
2786                                 delete wt_point;
2787
2788                                 return;
2789                         }
2790
2791                         d2::pixel get_depth(unsigned int i, unsigned int j) {
2792                                 assert(sv);
2793                                 return sv->get_depth(i, j);
2794                         }
2795
2796                         void get_median_depth_and_diff(d2::pixel *depth, d2::pixel *diff, unsigned int i, unsigned int j) {
2797                                 assert(sv);
2798                                 sv->get_median_depth_and_diff(depth, diff, i, j);
2799                         }
2800
2801                         void get_median_depth_and_diff(d2::pixel *_depth, d2::pixel *_diff, d2::point p) {
2802                                 if (sv) {
2803                                         sv->get_median_depth_and_diff(_depth, _diff, p);
2804                                         return;
2805                                 }
2806
2807                                 /*
2808                                  * Generate a local depth image of required radius.
2809                                  */
2810
2811                                 ale_pos radius = 1;
2812
2813                                 if (diff_median_radius + 1 > radius)
2814                                         radius = diff_median_radius + 1;
2815                                 if (depth_median_radius > radius)
2816                                         radius = depth_median_radius;
2817
2818                                 d2::point pl = p - d2::point(radius, radius);
2819                                 d2::point ph = p + d2::point(radius, radius);
2820                                 const d2::image *local_depth = depth(_pt, -1, 1, pl, ph);
2821
2822                                 /*
2823                                  * Find depth and diff at this point, check for
2824                                  * undefined values, and generate projections
2825                                  * of the image corners on the estimated normal
2826                                  * surface.
2827                                  */
2828
2829                                 d2::image *median_diffs = local_depth->fcdiff_median((int) floor(diff_median_radius));
2830                                 d2::image *median_depths = local_depth->medians((int) floor(depth_median_radius));
2831
2832                                 *_depth = median_depths->pix((int) radius, (int) radius);
2833                                 *_diff = median_diffs->pix((int) radius, (int) radius);
2834
2835                                 delete median_diffs;
2836                                 delete median_depths;
2837                                 delete local_depth;
2838                         }
2839                 };
2840
2841                 view get_view(ale_pos aperture, unsigned index, unsigned int randomization) {
2842                         if (randomization == 0) {
2843
2844                                 while (aperture_to_shared_views_map[aperture].size() <= index) {
2845                                         aperture_to_shared_views_map[aperture].push_back(shared_view(get_new_view(aperture)));
2846                                 }
2847
2848                                 return view(&(aperture_to_shared_views_map[aperture][index]));
2849                         }
2850
2851                         return view(NULL, get_new_view(aperture));
2852                 }
2853
2854                 view_generator(pt original_pt) {
2855                         this->original_pt = original_pt;
2856                 }
2857         };
2858
2859         /*
2860          * Unfiltered function
2861          */
2862         static const d2::image *view_nofilter_focus(pt _pt, int n) {
2863
2864                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
2865
2866                 if (n >= 0) {
2867                         assert((int) floor(d2::align::of(n).scaled_height())
2868                              == (int) floor(_pt.scaled_height()));
2869                         assert((int) floor(d2::align::of(n).scaled_width())
2870                              == (int) floor(_pt.scaled_width()));
2871                 }
2872
2873                 const d2::image *depths = depth(_pt, n);
2874
2875                 d2::image *im = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2876                                                (int) floor(_pt.scaled_width()), 3);
2877
2878                 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
2879
2880                 view_generator vg(_pt);
2881
2882                 for (unsigned int i = 0; i < im->height(); i++)
2883                 for (unsigned int j = 0; j < im->width();  j++) {
2884
2885                         focus::result _focus = focus::get(depths, i, j);
2886
2887                         if (!finite(_focus.focal_distance))
2888                                 continue;
2889
2890                         /*
2891                          * Data structures for calculating focal statistics.
2892                          */
2893
2894                         d2::pixel color, weight;
2895                         d2::image_weighted_median *iwm = NULL;
2896
2897                         if (_focus.statistic == 1) {
2898                                 iwm = new d2::image_weighted_median(1, 1, 3, _focus.sample_count);
2899                         }
2900
2901                         /*
2902                          * Iterate over views for this focus region.
2903                          */
2904
2905                         for (unsigned int v = 0; v < _focus.sample_count; v++) {
2906
2907                                 view_generator::view vw = vg.get_view(_focus.aperture, v, _focus.randomization);
2908
2909                                 ui::get()->d3_render_status(0, 1, -1, v, i, j, -1);
2910
2911
2912                                 /*
2913                                  * Map the focused point to the new view.
2914                                  */
2915
2916                                 point p = vw.get_pt().wp_scaled(_pt.pw_scaled(point(i, j, _focus.focal_distance)));
2917
2918                                 /*
2919                                  * Determine weight and color for the given point.
2920                                  */
2921
2922                                 d2::pixel view_weight, view_color;
2923
2924                                 vw.get_color_and_weight(&view_color, &view_weight, p.xy());
2925
2926                                 if (!color.finite() || !weight.finite())
2927                                         continue;
2928
2929                                 if (_focus.statistic == 0) {
2930                                         color += view_color;
2931                                         weight += view_weight;
2932                                 } else if (_focus.statistic == 1) {
2933                                         iwm->accumulate(0, 0, v, view_color, view_weight);
2934                                 } else
2935                                         assert(0);
2936                         }
2937
2938                         if (_focus.statistic == 1) {
2939                                 weight = iwm->get_weights()->get_pixel(0, 0);
2940                                 color = iwm->get_pixel(0, 0);
2941                                 delete iwm;
2942                         }
2943
2944                         if (weight.min_norm() < encounter_threshold) {
2945                                 im->pix(i, j) = d2::pixel::zero() / d2::pixel::zero();
2946                         } else if (normalize_weights)
2947                                 im->pix(i, j) = color / weight;
2948                         else
2949                                 im->pix(i, j) = color;
2950                 }
2951
2952                 delete depths;
2953
2954                 return im;
2955         }
2956
2957         /*
2958          * Unfiltered function
2959          */
2960         static const d2::image *view_nofilter(pt _pt, int n) {
2961
2962                 if (!focus::is_trivial())
2963                         return view_nofilter_focus(_pt, n);
2964
2965                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
2966
2967                 if (n >= 0) {
2968                         assert((int) floor(d2::align::of(n).scaled_height())
2969                              == (int) floor(_pt.scaled_height()));
2970                         assert((int) floor(d2::align::of(n).scaled_width())
2971                              == (int) floor(_pt.scaled_width()));
2972                 }
2973
2974                 const d2::image *depths = depth(_pt, n);
2975
2976                 d2::image *im = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2977                                                (int) floor(_pt.scaled_width()), 3);
2978
2979                 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
2980
2981                 /*
2982                  * Use adaptive subspace data.
2983                  */
2984
2985                 d2::image *weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2986                                                 (int) floor(_pt.scaled_width()), 3);
2987
2988                 /*
2989                  * Iterate through subspaces.
2990                  */
2991
2992                 space::iterate si(_pt.origin());
2993
2994                 ui::get()->d3_render_status(0, 0, -1, -1, -1, -1, 0);
2995
2996                 view_recurse(0, im, weights, si, _pt);
2997
2998                 for (unsigned int i = 0; i < im->height(); i++)
2999                 for (unsigned int j = 0; j < im->width();  j++) {
3000                         if (weights->pix(i, j).min_norm() < encounter_threshold
3001                          || (d3px_count > 0 && isnan(depths->pix(i, j)[0]))) {
3002                                 im->pix(i, j) = d2::pixel::zero() / d2::pixel::zero();
3003                                 weights->pix(i, j) = d2::pixel::zero();
3004                         } else if (normalize_weights)
3005                                 im->pix(i, j) /= weights->pix(i, j);
3006                 }
3007
3008                 delete weights;
3009
3010                 delete depths;
3011
3012                 return im;
3013         }
3014
3015         /*
3016          * Filtered function.
3017          */
3018         static const d2::image *view_filter_focus(pt _pt, int n) {
3019
3020                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
3021
3022                 /*
3023                  * Get depth image for focus region determination.
3024                  */
3025
3026                 const d2::image *depths = depth(_pt, n);
3027
3028                 unsigned int height = (unsigned int) floor(_pt.scaled_height());
3029                 unsigned int width = (unsigned int) floor(_pt.scaled_width());
3030
3031                 /*
3032                  * Prepare input frame data.
3033                  */
3034
3035                 if (tc_multiplier == 0)
3036                         al->open_all();
3037
3038                 pt *_ptf = new pt[al->count()];
3039                 std::vector<space::node *> *fmv = new std::vector<space::node *>[al->count()];
3040
3041                 for (unsigned int f = 0; f < al->count(); f++) {
3042                         _ptf[f] = al->get(f)->get_t(0);
3043                         fmv[f] = most_visible_unscaled(_ptf[f]);
3044                         std::sort(fmv[f].begin(), fmv[f].end());
3045                 }
3046
3047                 if (tc_multiplier == 0)
3048                         al->close_all();
3049
3050                 /*
3051                  * Open all files for rendering.
3052                  */
3053
3054                 d2::image_rw::open_all();
3055
3056                 /*
3057                  * Prepare data structures for averaging views, as we render
3058                  * each view separately.  This is spacewise inefficient, but
3059                  * is easy to implement given the current operation of the
3060                  * renderers.
3061                  */
3062
3063                 d2::image_weighted_avg *iwa;
3064
3065                 if (d3::focus::uses_medians()) {
3066                         iwa = new d2::image_weighted_median(height, width, 3, focus::max_samples());
3067                 } else {
3068                         iwa = new d2::image_weighted_simple(height, width, 3, new d2::invariant(NULL));
3069                 }
3070
3071                 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
3072
3073                 /*
3074                  * Prepare view generator.
3075                  */
3076
3077                 view_generator vg(_pt);
3078
3079                 /*
3080                  * Render views separately.  This is spacewise inefficient,
3081                  * but is easy to implement given the current operation of the
3082                  * renderers.
3083                  */
3084
3085                 for (unsigned int v = 0; v < focus::max_samples(); v++) {
3086
3087                         /*
3088                          * Generate a new 2D renderer for filtering.
3089                          */
3090
3091                         d2::render::reset();
3092                         d2::render *renderer = d2::render_parse::get(d3chain_type);
3093
3094                         renderer->init_point_renderer(height, width, 3);
3095
3096                         /*
3097                          * Iterate over output points.
3098                          */
3099
3100                         for (unsigned int i = 0; i < height; i++)
3101                         for (unsigned int j = 0; j < width; j++) {
3102
3103                                 focus::result _focus = focus::get(depths, i, j);
3104
3105                                 if (v >= _focus.sample_count)
3106                                         continue;
3107
3108                                 if (!finite(_focus.focal_distance))
3109                                         continue;
3110
3111                                 view_generator::view vw = vg.get_view(_focus.aperture, v, _focus.randomization);
3112
3113                                 pt _pt_new = vw.get_pt();
3114
3115                                 point p = _pt_new.wp_scaled(_pt.pw_scaled(point(i, j, _focus.focal_distance)));
3116
3117                                 /*
3118                                  * Determine the most-visible subspace.
3119                                  */
3120
3121                                 space::node *mv = vw.get_most_visible(p.xy());
3122
3123                                 if (mv == NULL)
3124                                         continue;
3125
3126                                 /*
3127                                  * Get median depth and diff.
3128                                  */
3129
3130                                 d2::pixel depth, diff;
3131
3132                                 vw.get_median_depth_and_diff(&depth, &diff, p.xy());
3133
3134                                 if (!depth.finite() || !diff.finite())
3135                                         continue;
3136
3137                                 point local_points[3] = {
3138                                         point(p[0],     p[1],     depth[0]),
3139                                         point(p[0] + 1, p[1],     depth[0] + diff[0]),
3140                                         point(p[0],     p[1] + 1, depth[0] + diff[1])
3141                                 };
3142
3143                                 /*
3144                                  * Iterate over files.
3145                                  */
3146
3147                                 for (unsigned int f = 0; f < d2::image_rw::count(); f++) {
3148
3149                                         ui::get()->d3_render_status(1, 1, f, v, i, j, -1);
3150
3151                                         if (!visibility_search(fmv[f], mv))
3152                                                 continue;
3153
3154                                         /*
3155                                          * Determine transformation at (i, j).  First
3156                                          * determine transformation from the output to
3157                                          * the input, then invert this, as we need the
3158                                          * inverse transformation for filtering.
3159                                          */
3160
3161                                         d2::point remote_points[3] = {
3162                                                 _ptf[f].wp_unscaled(_pt_new.pw_scaled(point(local_points[0]))).xy(),
3163                                                 _ptf[f].wp_unscaled(_pt_new.pw_scaled(point(local_points[1]))).xy(),
3164                                                 _ptf[f].wp_unscaled(_pt_new.pw_scaled(point(local_points[2]))).xy()
3165                                         };
3166
3167                                         /*
3168                                          * Forward matrix for the linear component of the
3169                                          * transformation.
3170                                          */
3171
3172                                         d2::point forward_matrix[2] = {
3173                                                 remote_points[1] - remote_points[0],
3174                                                 remote_points[2] - remote_points[0]
3175                                         };
3176
3177                                         /*
3178                                          * Inverse matrix for the linear component of
3179                                          * the transformation.  Calculate using the
3180                                          * determinant D.
3181                                          */
3182
3183                                         ale_pos D = forward_matrix[0][0] * forward_matrix[1][1]
3184                                                   - forward_matrix[0][1] * forward_matrix[1][0];
3185
3186                                         if (D == 0)
3187                                                 continue;
3188
3189                                         d2::point inverse_matrix[2] = {
3190                                                 d2::point( forward_matrix[1][1] / D, -forward_matrix[1][0] / D),
3191                                                 d2::point(-forward_matrix[0][1] / D,  forward_matrix[0][0] / D)
3192                                         };
3193
3194                                         /*
3195                                          * Determine the projective transformation parameters for the
3196                                          * inverse transformation.
3197                                          */
3198
3199                                         const d2::image *imf = d2::image_rw::get_open(f);
3200
3201                                         d2::transformation inv_t = d2::transformation::gpt_identity(imf, 1);
3202
3203                                         d2::point local_bounds[4];
3204
3205                                         for (int n = 0; n < 4; n++) {
3206                                                 d2::point remote_bound = d2::point((n == 1 || n == 2) ? imf->height() : 0,
3207                                                                                    (n == 2 || n == 3) ? imf->width()  : 0)
3208                                                                        - remote_points[0];
3209
3210                                                 local_bounds[n] = d2::point(i, j)
3211                                                                 + d2::point(remote_bound[0] * inverse_matrix[0][0]
3212                                                                           + remote_bound[1] * inverse_matrix[1][0],
3213                                                                             remote_bound[0] * inverse_matrix[0][1]
3214                                                                           + remote_bound[1] * inverse_matrix[1][1]);
3215
3216                                         }
3217
3218                                         if (!local_bounds[0].finite()
3219                                          || !local_bounds[1].finite()
3220                                          || !local_bounds[2].finite()
3221                                          || !local_bounds[3].finite())
3222                                                 continue;
3223
3224                                         inv_t.gpt_set(local_bounds);
3225
3226                                         /*
3227                                          * Perform render step for the given frame,
3228                                          * transformation, and point.
3229                                          */
3230
3231                                         renderer->point_render(i, j, f, inv_t);
3232                                 }
3233                         }
3234
3235                         renderer->finish_point_rendering();
3236
3237                         const d2::image *im = renderer->get_image();
3238                         const d2::image *df = renderer->get_defined();
3239
3240                         for (unsigned int i = 0; i < height; i++)
3241                         for (unsigned int j = 0; j < width; j++) {
3242                                 if (df->get_pixel(i, j).finite()
3243                                  && df->get_pixel(i, j)[0] > 0)
3244                                         iwa->accumulate(i, j, v, im->get_pixel(i, j), d2::pixel(1, 1, 1));
3245                         }
3246                 }
3247
3248                 /*
3249                  * Close all files and return the result.
3250                  */
3251
3252                 d2::image_rw::close_all();
3253
3254                 return iwa;
3255         }
3256
3257         static const d2::image *view_filter(pt _pt, int n) {
3258
3259                 if (!focus::is_trivial())
3260                         return view_filter_focus(_pt, n);
3261
3262                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
3263
3264                 /*
3265                  * Generate a new 2D renderer for filtering.
3266                  */
3267
3268                 d2::render::reset();
3269                 d2::render *renderer = d2::render_parse::get(d3chain_type);
3270
3271                 /*
3272                  * Get depth image in order to estimate normals (and hence
3273                  * transformations).
3274                  */
3275
3276                 const d2::image *depths = depth(_pt, n);
3277
3278                 d2::image *median_diffs = depths->fcdiff_median((int) floor(diff_median_radius));
3279                 d2::image *median_depths = depths->medians((int) floor(depth_median_radius));
3280
3281                 unsigned int height = (unsigned int) floor(_pt.scaled_height());
3282                 unsigned int width = (unsigned int) floor(_pt.scaled_width());
3283
3284                 renderer->init_point_renderer(height, width, 3);
3285
3286                 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
3287
3288                 std::vector<space::node *> mv = most_visible_scaled(_pt);
3289
3290                 for (unsigned int f = 0; f < d2::image_rw::count(); f++) {
3291
3292                         if (tc_multiplier == 0)
3293                                 al->open(f);
3294
3295                         pt _ptf = al->get(f)->get_t(0);
3296
3297                         std::vector<space::node *> fmv = most_visible_unscaled(_ptf);
3298                         std::sort(fmv.begin(), fmv.end());
3299
3300                         for (unsigned int i = 0; i < height; i++)
3301                         for (unsigned int j = 0; j < width; j++) {
3302
3303                                 ui::get()->d3_render_status(1, 0, f, -1, i, j, -1);
3304
3305                                 /*
3306                                  * Check visibility.
3307                                  */
3308
3309                                 int n = i * width + j;
3310
3311                                 if (!visibility_search(fmv, mv[n]))
3312                                         continue;
3313
3314                                 /*
3315                                  * Find depth and diff at this point, check for
3316                                  * undefined values, and generate projections
3317                                  * of the image corners on the estimated normal
3318                                  * surface.
3319                                  */
3320
3321                                 d2::pixel depth = median_depths->pix(i, j);
3322                                 d2::pixel diff = median_diffs->pix(i, j);
3323                                 // d2::pixel diff = d2::pixel(0, 0, 0);
3324
3325                                 if (!depth.finite() || !diff.finite())
3326                                         continue;
3327
3328                                 point local_points[3] = {
3329                                         point(i,     j,     depth[0]),
3330                                         point(i + 1, j,     depth[0] + diff[0]),
3331                                         point(i    , j + 1, depth[0] + diff[1])
3332                                 };
3333
3334                                 /*
3335                                  * Determine transformation at (i, j).  First
3336                                  * determine transformation from the output to
3337                                  * the input, then invert this, as we need the
3338                                  * inverse transformation for filtering.
3339                                  */
3340
3341                                 d2::point remote_points[3] = {
3342                                         _ptf.wp_unscaled(_pt.pw_scaled(point(local_points[0]))).xy(),
3343                                         _ptf.wp_unscaled(_pt.pw_scaled(point(local_points[1]))).xy(),
3344                                         _ptf.wp_unscaled(_pt.pw_scaled(point(local_points[2]))).xy()
3345                                 };
3346
3347                                 /*
3348                                  * Forward matrix for the linear component of the
3349                                  * transformation.
3350                                  */
3351
3352                                 d2::point forward_matrix[2] = {
3353                                         remote_points[1] - remote_points[0],
3354                                         remote_points[2] - remote_points[0]
3355                                 };
3356
3357                                 /*
3358                                  * Inverse matrix for the linear component of
3359                                  * the transformation.  Calculate using the
3360                                  * determinant D.
3361                                  */
3362
3363                                 ale_pos D = forward_matrix[0][0] * forward_matrix[1][1]
3364                                           - forward_matrix[0][1] * forward_matrix[1][0];
3365
3366                                 if (D == 0)
3367                                         continue;
3368
3369                                 d2::point inverse_matrix[2] = {
3370                                         d2::point( forward_matrix[1][1] / D, -forward_matrix[1][0] / D),
3371                                         d2::point(-forward_matrix[0][1] / D,  forward_matrix[0][0] / D)
3372                                 };
3373
3374                                 /*
3375                                  * Determine the projective transformation parameters for the
3376                                  * inverse transformation.
3377                                  */
3378
3379                                 const d2::image *imf = d2::image_rw::open(f);
3380
3381                                 d2::transformation inv_t = d2::transformation::gpt_identity(imf, 1);
3382
3383                                 d2::point local_bounds[4];
3384
3385                                 for (int n = 0; n < 4; n++) {
3386                                         d2::point remote_bound = d2::point((n == 1 || n == 2) ? imf->height() : 0,
3387                                                                            (n == 2 || n == 3) ? imf->width()  : 0)
3388                                                                - remote_points[0];
3389
3390                                         local_bounds[n] = local_points[0].xy()
3391                                                         + d2::point(remote_bound[0] * inverse_matrix[0][0]
3392                                                                   + remote_bound[1] * inverse_matrix[1][0],
3393                                                                     remote_bound[0] * inverse_matrix[0][1]
3394                                                                   + remote_bound[1] * inverse_matrix[1][1]);
3395                                 }
3396
3397                                 inv_t.gpt_set(local_bounds);
3398
3399                                 d2::image_rw::close(f);
3400
3401                                 /*
3402                                  * Perform render step for the given frame,
3403                                  * transformation, and point.
3404                                  */
3405
3406                                 d2::image_rw::open(f);
3407                                 renderer->point_render(i, j, f, inv_t);
3408                                 d2::image_rw::close(f);
3409                         }
3410
3411                         if (tc_multiplier == 0)
3412                                 al->close(f);
3413                 }
3414
3415                 renderer->finish_point_rendering();
3416
3417                 return renderer->get_image();
3418         }
3419
3420         /*
3421          * Generic function.
3422          */
3423         static const d2::image *view(pt _pt, int n = -1) {
3424
3425                 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
3426
3427                 if (use_filter) {
3428                         return view_filter(_pt, n);
3429                 } else {
3430                         return view_nofilter(_pt, n);
3431                 }
3432         }
3433
3434         static void tcem(double _tcem) {
3435                 tc_multiplier = _tcem;
3436         }
3437
3438         static void oui(unsigned int _oui) {
3439                 ou_iterations = _oui;
3440         }
3441
3442         static void pa(unsigned int _pa) {
3443                 pairwise_ambiguity = _pa;
3444         }
3445
3446         static void pc(const char *_pc) {
3447                 pairwise_comparisons = _pc;
3448         }
3449
3450         static void d3px(int _d3px_count, double *_d3px_parameters) {
3451                 d3px_count = _d3px_count;
3452                 d3px_parameters = _d3px_parameters;
3453         }
3454
3455         static void fx(double _fx) {
3456                 falloff_exponent = _fx;
3457         }
3458
3459         static void nw() {
3460                 normalize_weights = 1;
3461         }
3462
3463         static void no_nw() {
3464                 normalize_weights = 0;
3465         }
3466
3467         static void nofilter() {
3468                 use_filter = 0;
3469         }
3470
3471         static void filter() {
3472                 use_filter = 1;
3473         }
3474
3475         static void set_filter_type(const char *type) {
3476                 d3chain_type = type;
3477         }
3478
3479         static void set_subspace_traverse() {
3480                 subspace_traverse = 1;
3481         }
3482
3483         static int excluded(point p) {
3484                 for (int n = 0; n < d3px_count; n++) {
3485                         double *region = d3px_parameters + (6 * n);
3486                         if (p[0] >= region[0]
3487                          && p[0] <= region[1]
3488                          && p[1] >= region[2]
3489                          && p[1] <= region[3]
3490                          && p[2] >= region[4]
3491                          && p[2] <= region[5])
3492                                 return 1;
3493                 }
3494
3495                 return 0;
3496         }
3497
3498         /*
3499          * This function returns true if a space is completely excluded.
3500          */
3501         static int excluded(const space::traverse &st) {
3502                 for (int n = 0; n < d3px_count; n++) {
3503                         double *region = d3px_parameters + (6 * n);
3504                         if (st.get_min()[0] >= region[0]
3505                          && st.get_max()[0] <= region[1]
3506                          && st.get_min()[1] >= region[2]
3507                          && st.get_max()[1] <= region[3]
3508                          && st.get_min()[2] >= region[4]
3509                          && st.get_max()[2] <= region[5])
3510                                 return 1;
3511                 }
3512
3513                 return 0;
3514         }
3515
3516         static const d2::image *view(unsigned int n) {
3517
3518                 assert (n < d2::image_rw::count());
3519
3520                 pt _pt = align::projective(n);
3521
3522                 return view(_pt, n);
3523         }
3524
3525         typedef struct {point iw; point ip, is;} analytic;
3526         typedef std::multimap<ale_real,analytic> score_map;
3527         typedef std::pair<ale_real,analytic> score_map_element;
3528
3529         /*
3530          * Make pt list.
3531          */
3532         static std::vector<pt> make_pt_list(const char *d_out[], const char *v_out[],
3533                         std::map<const char *, pt> *d3_depth_pt,
3534                         std::map<const char *, pt> *d3_output_pt) {
3535
3536                 std::vector<pt> result;
3537
3538                 for (unsigned int n = 0; n < d2::image_rw::count(); n++) {
3539                         if (d_out[n] || v_out[n]) {
3540                                 result.push_back(align::projective(n));
3541                         }
3542                 }
3543
3544                 for (std::map<const char *, pt>::iterator i = d3_depth_pt->begin(); i != d3_depth_pt->end(); i++) {
3545                         result.push_back(i->second);
3546                 }
3547
3548                 for (std::map<const char *, pt>::iterator i = d3_output_pt->begin(); i != d3_output_pt->end(); i++) {
3549                         result.push_back(i->second);
3550                 }
3551
3552                 return result;
3553         }
3554
3555         /*
3556          * Get a trilinear coordinate for an anisotropic candidate cell.
3557          */
3558         static ale_pos get_trilinear_coordinate(point min, point max, pt _pt) {
3559
3560                 d2::point local_min, local_max;
3561
3562                 local_min = _pt.wp_unscaled(min).xy();
3563                 local_max = _pt.wp_unscaled(min).xy();
3564
3565                 point cell[2] = {min, max};
3566
3567                 /*
3568                  * Determine the view-local extrema in 2 dimensions.
3569                  */
3570
3571                 for (int r = 1; r < 8; r++) {
3572                         point local = _pt.wp_unscaled(point(cell[r>>2][0], cell[(r>>1)%2][1], cell[r%2][2]));
3573
3574                         for (int d = 0; d < 2; d++) {
3575                                 if (local[d] < local_min[d])
3576                                         local_min[d] = local[d];
3577                                 if (local[d] > local_max[d])
3578                                         local_max[d] = local[d];
3579                                 if (isnan(local[d]))
3580                                         return local[d];
3581                         }
3582                 }
3583
3584                 ale_pos diameter = (local_max - local_min).norm();
3585
3586                 return log(diameter / sqrt(2)) / log(2);
3587         }
3588
3589         /*
3590          * Check whether a cell is visible from a given viewpoint.  This function
3591          * is guaranteed to return 1 when a cell is visible, but it is not guaranteed
3592          * to return 0 when a cell is invisible.
3593          */
3594         static int pt_might_be_visible(const pt &viewpoint, point min, point max) {
3595
3596                 int doc = (rand() % 100000) ? 0 : 1;
3597
3598                 if (doc)
3599                         fprintf(stderr, "checking visibility:\n");
3600
3601                 point cell[2] = {min, max};
3602
3603                 /*
3604                  * Cycle through all vertices of the cell to check certain
3605                  * properties.
3606                  */
3607                 int pos[3] = {0, 0, 0};
3608                 int neg[3] = {0, 0, 0};
3609                 for (int i = 0; i < 2; i++)
3610                 for (int j = 0; j < 2; j++)
3611                 for (int k = 0; k < 2; k++) {
3612                         point p = viewpoint.wp_unscaled(point(cell[i][0], cell[j][1], cell[k][2]));
3613
3614                         if (p[2] < 0 && viewpoint.unscaled_in_bounds(p))
3615                                 return 1;
3616
3617                         if (isnan(p[0])
3618                          || isnan(p[1])
3619                          || isnan(p[2]))
3620                                 return 1;
3621
3622                         if (p[2] > 0)
3623                                 for (int d = 0; d < 2; d++)
3624                                         p[d] *= -1;
3625
3626                         if (doc)
3627                                 fprintf(stderr, "\t[%f %f %f] --> [%f %f %f]\n",
3628                                                 cell[i][0], cell[j][1], cell[k][2],
3629                                                 p[0], p[1], p[2]);
3630
3631                         for (int d = 0; d < 3; d++)
3632                                 if (p[d] >= 0)
3633                                         pos[d] = 1;
3634
3635                         if (p[0] <= viewpoint.unscaled_height() - 1)
3636                                 neg[0] = 1;
3637
3638                         if (p[1] <= viewpoint.unscaled_width() - 1)
3639                                 neg[1] = 1;
3640
3641                         if (p[2] <= 0)
3642                                 neg[2] = 1;
3643                 }
3644
3645                 if (!neg[2])
3646                         return 0;
3647
3648                 if (!pos[0]
3649                  || !neg[0]
3650                  || !pos[1]
3651                  || !neg[1])
3652                         return 0;
3653
3654                 return 1;
3655         }
3656
3657         /*
3658          * Check whether a cell is output-visible.
3659          */
3660         static int output_might_be_visible(const std::vector<pt> &pt_outputs, point min, point max) {
3661                 for (unsigned int n = 0; n < pt_outputs.size(); n++)
3662                         if (pt_might_be_visible(pt_outputs[n], min, max))
3663                                 return 1;
3664                 return 0;
3665         }
3666
3667         /*
3668          * Check whether a cell is input-visible.
3669          */
3670         static int input_might_be_visible(unsigned int f, point min, point max) {
3671                 return pt_might_be_visible(align::projective(f), min, max);
3672         }
3673
3674         /*
3675          * Return true if a cell fails an output resolution bound.
3676          */
3677         static int fails_output_resolution_bound(point min, point max, const std::vector<pt> &pt_outputs) {
3678                 for (unsigned int n = 0; n < pt_outputs.size(); n++) {
3679
3680                         point p = pt_outputs[n].centroid(min, max);
3681
3682                         if (!p.defined())
3683                                 continue;
3684
3685                         if (get_trilinear_coordinate(min, max, pt_outputs[n]) < output_decimation_preferred)
3686                                 return 1;
3687                 }
3688
3689                 return 0;
3690         }
3691
3692         /*
3693          * Check lower-bound resolution constraints
3694          */
3695         static int exceeds_resolution_lower_bounds(unsigned int f1, unsigned int f2,
3696                         point min, point max, const std::vector<pt> &pt_outputs) {
3697
3698                 pt _pt = al->get(f1)->get_t(0);
3699
3700                 if (get_trilinear_coordinate(min, max, _pt) < input_decimation_lower)
3701                         return 1;
3702
3703                 if (fails_output_resolution_bound(min, max, pt_outputs))
3704                         return 0;
3705
3706                 if (get_trilinear_coordinate(min, max, _pt) < primary_decimation_upper)
3707                         return 1;
3708
3709                 return 0;
3710         }
3711
3712         /*
3713          * Try the candidate nearest to the specified cell.
3714          */
3715         static void try_nearest_candidate(unsigned int f1, unsigned int f2, candidates *c, point min, point max) {
3716                 point centroid = (max + min) / 2;
3717                 pt _pt[2] = { al->get(f1)->get_t(0), al->get(f2)->get_t(0) };
3718                 point p[2];
3719
3720                 // fprintf(stderr, "[tnc n=%f %f %f x=%f %f %f]\n", min[0], min[1], min[2], max[0], max[1], max[2]);
3721
3722                 /*
3723                  * Reject clipping plane violations.
3724                  */
3725
3726                 if (centroid[2] > front_clip
3727                  || centroid[2] < rear_clip)
3728                         return;
3729
3730                 /*
3731                  * Calculate projections.
3732                  */
3733
3734                 for (int n = 0; n < 2; n++) {
3735
3736                         p[n] = _pt[n].wp_unscaled(centroid);
3737
3738                         if (!_pt[n].unscaled_in_bounds(p[n]))
3739                                 return;
3740
3741                         // fprintf(stderr, ":");
3742
3743                         if (p[n][2] >= 0)
3744                                 return;
3745                 }
3746
3747
3748                 int tc = (int) round(get_trilinear_coordinate(min, max, _pt[0]));
3749                 int stc = (int) round(get_trilinear_coordinate(min, max, _pt[1]));
3750
3751                 while (tc < input_decimation_lower || stc < input_decimation_lower) {
3752                         tc++;
3753                         stc++;
3754                 }
3755
3756                 if (tc > primary_decimation_upper)
3757                         return;
3758
3759                 /*
3760                  * Calculate score from color match.  Assume for now
3761                  * that the transformation can be approximated locally
3762                  * with a translation.
3763                  */
3764
3765                 ale_pos score = 0;
3766                 ale_pos divisor = 0;
3767                 ale_pos l1_multiplier = 0.125;
3768                 lod_image *if1 = al->get(f1);
3769                 lod_image *if2 = al->get(f2);
3770
3771                 if (if1->in_bounds(p[0].xy())
3772                  && if2->in_bounds(p[1].xy())) {
3773                         divisor += 1 - l1_multiplier;
3774                         score += (1 - l1_multiplier)
3775                                * (if1->get_tl(p[0].xy(), tc) - if2->get_tl(p[1].xy(), stc)).normsq();
3776                 }
3777
3778                 for (int iii = -1; iii <= 1; iii++)
3779                 for (int jjj = -1; jjj <= 1; jjj++) {
3780                         d2::point t(iii, jjj);
3781
3782                         if (!if1->in_bounds(p[0].xy() + t)
3783                          || !if2->in_bounds(p[1].xy() + t))
3784                                 continue;
3785
3786                         divisor += l1_multiplier;
3787                         score   += l1_multiplier
3788                                  * (if1->get_tl(p[0].xy() + t, tc) - if2->get_tl(p[1].xy() + t, tc)).normsq();
3789
3790                 }
3791
3792                 /*
3793                  * Include third-camera contributions in the score.
3794                  */
3795
3796                 if (tc_multiplier != 0)
3797                 for (unsigned int n = 0; n < d2::image_rw::count(); n++) {
3798                         if (n == f1 || n == f2)
3799                                 continue;
3800
3801                         lod_image *ifn = al->get(n);
3802                         pt _ptn = ifn->get_t(0);
3803                         point pn = _ptn.wp_unscaled(centroid);
3804
3805                         if (!_ptn.unscaled_in_bounds(pn))
3806                                 continue;
3807
3808                         if (pn[2] >= 0)
3809                                 continue;
3810
3811                         ale_pos ttc = get_trilinear_coordinate(min, max, _ptn);
3812
3813                         divisor += tc_multiplier;
3814                         score   += tc_multiplier
3815                                  * (if1->get_tl(p[0].xy(), tc) - ifn->get_tl(pn.xy(), ttc)).normsq();
3816                 }
3817
3818                 c->add_candidate(p[0], tc, score / divisor);
3819         }
3820
3821         /*
3822          * Check for cells that are completely clipped.
3823          */
3824         static int completely_clipped(point min, point max) {
3825                 return (min[2] > front_clip
3826                      || max[2] < rear_clip);
3827         }
3828
3829         /*
3830          * Update extremum variables for cell points mapped to a particular view.
3831          */
3832         static void update_extrema(point min, point max, pt _pt, int *extreme_dim, ale_pos *extreme_ratio) {
3833
3834                 point local_min, local_max;
3835
3836                 local_min = _pt.wp_unscaled(min);
3837                 local_max = _pt.wp_unscaled(min);
3838
3839                 point cell[2] = {min, max};
3840
3841                 int near_vertex = 0;
3842
3843                 /*
3844                  * Determine the view-local extrema in all dimensions, and
3845                  * determine the vertex of closest z coordinate.
3846                  */
3847
3848                 for (int r = 1; r < 8; r++) {
3849                         point local = _pt.wp_unscaled(point(cell[r>>2][0], cell[(r>>1)%2][1], cell[r%2][2]));
3850
3851                         for (int d = 0; d < 3; d++) {
3852                                 if (local[d] < local_min[d])
3853                                         local_min[d] = local[d];
3854                                 if (local[d] > local_max[d])
3855                                         local_max[d] = local[d];
3856                         }
3857
3858                         if (local[2] == local_max[2])
3859                                 near_vertex = r;
3860                 }
3861
3862                 ale_pos diameter = (local_max.xy() - local_min.xy()).norm();
3863
3864                 /*
3865                  * Update extrema as necessary for each dimension.
3866                  */
3867
3868                 for (int d = 0; d < 3; d++) {
3869
3870                         int r = near_vertex;
3871
3872                         int p1[3] = {r>>2, (r>>1)%2, r%2};
3873                         int p2[3] = {r>>2, (r>>1)%2, r%2};
3874
3875                         p2[d] = 1 - p2[d];
3876
3877                         ale_pos local_distance = (_pt.wp_unscaled(point(cell[p1[0]][0], cell[p1[1]][1], cell[p1[2]][2])).xy()
3878                                                 - _pt.wp_unscaled(point(cell[p2[0]][0], cell[p2[1]][1], cell[p2[2]][2])).xy()).norm();
3879
3880                         if (local_distance / diameter > *extreme_ratio) {
3881                                 *extreme_ratio = local_distance / diameter;
3882                                 *extreme_dim = d;
3883                         }
3884                 }
3885         }
3886
3887         /*
3888          * Get the next split dimension.
3889          */
3890         static int get_next_split(int f1, int f2, point min, point max, const std::vector<pt> &pt_outputs) {
3891                 for (int d = 0; d < 3; d++)
3892                         if (isinf(min[d]) || isinf(max[d]))
3893                                 return space::traverse::get_next_split(min, max);
3894
3895                 int extreme_dim = 0;
3896                 ale_pos extreme_ratio = 0;
3897
3898                 update_extrema(min, max, al->get(f1)->get_t(0), &extreme_dim, &extreme_ratio);
3899                 update_extrema(min, max, al->get(f2)->get_t(0), &extreme_dim, &extreme_ratio);
3900
3901                 for (unsigned int n = 0; n < pt_outputs.size(); n++) {
3902                         update_extrema(min, max, pt_outputs[n], &extreme_dim, &extreme_ratio);
3903                 }
3904
3905                 return extreme_dim;
3906         }
3907
3908         /*
3909          * Find candidates for subspace creation.
3910          */
3911         static void find_candidates(unsigned int f1, unsigned int f2, candidates *c, point min, point max,
3912                         const std::vector<pt> &pt_outputs, int depth = 0) {
3913
3914                 int print = 0;
3915
3916                 if (min[0] < 20.0001 && max[0] > 20.0001
3917                  && min[1] < 20.0001 && max[1] > 20.0001
3918                  && min[2] < 0.0001 && max[2] > 0.0001)
3919                         print = 1;
3920
3921                 if (print) {
3922                         for (int i = depth; i > 0; i--) {
3923                                 fprintf(stderr, "+");
3924                         }
3925                         fprintf(stderr, "[fc n=%f %f %f x=%f %f %f]\n",
3926                                         min[0], min[1], min[2], max[0], max[1], max[2]);
3927                 }
3928
3929                 if (completely_clipped(min, max)) {
3930                         if (print)
3931                                 fprintf(stderr, "c");
3932                         return;
3933                 }
3934
3935                 if (!input_might_be_visible(f1, min, max)
3936                  || !input_might_be_visible(f2, min, max)) {
3937                         if (print)
3938                                 fprintf(stderr, "v");
3939                         return;
3940                 }
3941
3942                 if (output_clip && !output_might_be_visible(pt_outputs, min, max)) {
3943                         if (print)
3944                                 fprintf(stderr, "o");
3945                         return;
3946                 }
3947
3948                 if (exceeds_resolution_lower_bounds(f1, f2, min, max, pt_outputs)) {
3949                         if (!(rand() % 100000))
3950                                 fprintf(stderr, "([%f %f %f], [%f %f %f]) at %d\n",
3951                                                 min[0], min[1], min[2],
3952                                                 max[0], max[1], max[2],
3953                                                 __LINE__);
3954
3955                         if (print)
3956                                 fprintf(stderr, "t");
3957
3958                         try_nearest_candidate(f1, f2, c, min, max);
3959                         return;
3960                 }
3961
3962                 point new_cells[2][2];
3963
3964                 if (!space::traverse::get_next_cells(get_next_split(f1, f2, min, max, pt_outputs), min, max, new_cells)) {
3965                         if (print)
3966                                 fprintf(stderr, "n");
3967                         return;
3968                 }
3969
3970                 if (print) {
3971                         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",
3972                                         new_cells[0][0][0],
3973                                         new_cells[0][0][1],
3974                                         new_cells[0][0][2],
3975                                         new_cells[0][1][0],
3976                                         new_cells[0][1][1],
3977                                         new_cells[0][1][2],
3978                                         new_cells[1][0][0],
3979                                         new_cells[1][0][1],
3980                                         new_cells[1][0][2],
3981                                         new_cells[1][1][0],
3982                                         new_cells[1][1][1],
3983                                         new_cells[1][1][2]);
3984                 }
3985
3986                 find_candidates(f1, f2, c, new_cells[0][0], new_cells[0][1], pt_outputs, depth + 1);
3987                 find_candidates(f1, f2, c, new_cells[1][0], new_cells[1][1], pt_outputs, depth + 1);
3988         }
3989
3990         /*
3991          * Generate a map from scores to 3D points for various depths at point (i, j) in f1, at
3992          * lowest resolution.
3993          */
3994         static score_map p2f_score_map(unsigned int f1, unsigned int f2, unsigned int i, unsigned int j) {
3995
3996                 score_map result;
3997
3998                 pt _pt1 = al->get(f1)->get_t(primary_decimation_upper);
3999                 pt _pt2 = al->get(f2)->get_t(primary_decimation_upper);
4000
4001                 const d2::image *if1 = al->get(f1)->get_image(primary_decimation_upper);
4002                 const d2::image *if2 = al->get(f2)->get_image(primary_decimation_upper);
4003                 ale_pos pdu_scale = pow(2, primary_decimation_upper);
4004
4005                 /*
4006                  * Get the pixel color in the primary frame
4007                  */
4008
4009                 // d2::pixel color_primary = if1->get_pixel(i, j);
4010
4011                 /*
4012                  * Map two depths to the secondary frame.
4013                  */
4014
4015                 point p1 = _pt2.wp_unscaled(_pt1.pw_unscaled(point(i, j, 1000)));
4016                 point p2 = _pt2.wp_unscaled(_pt1.pw_unscaled(point(i, j, -1000)));
4017
4018 //              fprintf(stderr, "%d->%d (%d, %d) point pair: (%d, %d, %d -> %f, %f), (%d, %d, %d -> %f, %f)\n",
4019 //                              f1, f2, i, j, i, j, 1000, p1[0], p1[1], i, j, -1000, p2[0], p2[1]);
4020 //              _pt1.debug_output();
4021 //              _pt2.debug_output();
4022
4023
4024                 /*
4025                  * For cases where the mapped points define a
4026                  * line and where points on the line fall
4027                  * within the defined area of the frame,
4028                  * determine the starting point for inspection.
4029                  * In other cases, continue to the next pixel.
4030                  */
4031
4032                 ale_pos diff_i = p2[0] - p1[0];
4033                 ale_pos diff_j = p2[1] - p1[1];
4034                 ale_pos slope = diff_j / diff_i;
4035
4036                 if (isnan(slope)) {
4037                         assert(0);
4038                         fprintf(stderr, "%d->%d (%d, %d) has undefined slope\n",
4039                                         f1, f2, i, j);
4040                         return result;
4041                 }
4042
4043                 /*
4044                  * Make absurdly large/small slopes either infinity, negative infinity, or zero.
4045                  */
4046
4047                 if (fabs(slope) > if2->width() * 100) {
4048                         double zero = 0;
4049                         double one  = 1;
4050                         double inf  = one / zero;
4051                         slope = inf;
4052                 } else if (slope < 1 / (double) if2->height() / 100
4053                         && slope > -1/ (double) if2->height() / 100) {
4054                         slope = 0;
4055                 }
4056
4057                 // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4058
4059                 ale_pos top_intersect = p1[1] - p1[0] * slope;
4060                 ale_pos lef_intersect = p1[0] - p1[1] / slope;
4061                 ale_pos rig_intersect = p1[0] - (p1[1] - if2->width() + 2) / slope;
4062                 ale_pos sp_i, sp_j;
4063
4064                 // fprintf(stderr, "slope == %f\n", slope);
4065
4066
4067                 if (slope == 0) {
4068                         // fprintf(stderr, "case 0\n");
4069                         sp_i = lef_intersect;
4070                         sp_j = 0;
4071                 } else if (finite(slope) && top_intersect >= 0 && top_intersect < if2->width() - 1) {
4072                         // fprintf(stderr, "case 1\n");
4073                         sp_i = 0;
4074                         sp_j = top_intersect;
4075                 } else if (slope > 0 && lef_intersect >= 0 && lef_intersect <= if2->height() - 1) {
4076                         // fprintf(stderr, "case 2\n");
4077                         sp_i = lef_intersect;
4078                         sp_j = 0;
4079                 } else if (slope < 0 && rig_intersect >= 0 && rig_intersect <= if2->height() - 1) {
4080                         // fprintf(stderr, "case 3\n");
4081                         sp_i = rig_intersect;
4082                         sp_j = if2->width() - 2;
4083                 } else {
4084                         // fprintf(stderr, "case 4\n");
4085                         // fprintf(stderr, "%d->%d (%d, %d) does not intersect the defined area\n",
4086                         //              f1, f2, i, j);
4087                         return result;
4088                 }
4089
4090
4091                 // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4092
4093                 /*
4094                  * Determine increment values for examining
4095                  * point, ensuring that incr_i is always
4096                  * positive.
4097                  */
4098
4099                 ale_pos incr_i, incr_j;
4100
4101                 if (fabs(diff_i) > fabs(diff_j)) {
4102                         incr_i = 1;
4103                         incr_j = slope;
4104                 } else if (slope > 0) {
4105                         incr_i = 1 / slope;
4106                         incr_j = 1;
4107                 } else {
4108                         incr_i = -1 / slope;
4109                         incr_j = -1;
4110                 }
4111
4112                 // fprintf(stderr, "%d->%d (%d, %d) increments are (%f, %f)\n",
4113                 //              f1, f2, i, j, incr_i, incr_j);
4114
4115                 /*
4116                  * Examine regions near the projected line.
4117                  */
4118
4119                 for (ale_pos ii = sp_i, jj = sp_j;
4120                         ii <= if2->height() - 1 && jj <= if2->width() - 1 && ii >= 0 && jj >= 0;
4121                         ii += incr_i, jj += incr_j) {
4122
4123                         // fprintf(stderr, "%d->%d (%d, %d) checking (%f, %f)\n",
4124                         //              f1, f2, i, j, ii, jj);
4125
4126 #if 0
4127                         /*
4128                          * Check for higher, lower, and nearby points.
4129                          *
4130                          *      Red   = 2^0
4131                          *      Green = 2^1
4132                          *      Blue  = 2^2
4133                          */
4134
4135                         int higher = 0, lower = 0, nearby = 0;
4136
4137                         for (int iii = 0; iii < 2; iii++)
4138                         for (int jjj = 0; jjj < 2; jjj++) {
4139                                 d2::pixel p = if2->get_pixel((int) floor(ii) + iii, (int) floor(jj) + jjj);
4140
4141                                 for (int k = 0; k < 3; k++) {
4142                                         int bitmask = (int) pow(2, k);
4143
4144                                         if (p[k] > color_primary[k])
4145                                                 higher |= bitmask;
4146                                         if (p[k] < color_primary[k])
4147                                                 lower  |= bitmask;
4148                                         if (fabs(p[k] - color_primary[k]) < nearness)
4149                                                 nearby |= bitmask;
4150                                 }
4151                         }
4152
4153                         /*
4154                          * If this is not a region of interest,
4155                          * then continue.
4156                          */
4157
4158
4159                         fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4160
4161                         // if (((higher & lower) | nearby) != 0x7)
4162                         //      continue;
4163 #endif
4164                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4165
4166                         // fprintf(stderr, "%d->%d (%d, %d) accepted (%f, %f)\n",
4167                         //              f1, f2, i, j, ii, jj);
4168
4169                         /*
4170                          * Create an orthonormal basis to
4171                          * determine line intersection.
4172                          */
4173
4174                         point bp0 = _pt1.pw_unscaled(point(i, j, 0));
4175                         point bp1 = _pt1.pw_unscaled(point(i, j, 10));
4176                         point bp2 = _pt2.pw_unscaled(point(ii, jj, 0));
4177
4178                         point foo = _pt1.wp_unscaled(bp0);
4179                         // fprintf(stderr, "(%d, %d, 0) transformed to world and back is: (%f, %f, %f)\n",
4180                         //              i, j, foo[0], foo[1], foo[2]);
4181
4182                         foo = _pt1.wp_unscaled(bp1);
4183                         // fprintf(stderr, "(%d, %d, 10)  transformed to world and back is: (%f, %f, %f)\n",
4184                         //              i, j, foo[0], foo[1], foo[2]);
4185
4186                         point b0  = (bp1 - bp0).normalize();
4187                         point b1n = bp2 - bp0;
4188                         point b1  = (b1n - b1n.dproduct(b0) * b0).normalize();
4189                         point b2  = point(0, 0, 0).xproduct(b0, b1).normalize(); // Should already have norm=1
4190
4191
4192                         foo = _pt1.wp_unscaled(bp0 + 30 * b0);
4193
4194                         /*
4195                          * Select a fourth point to define a second line.
4196                          */
4197
4198                         point p3  = _pt2.pw_unscaled(point(ii, jj, 10));
4199
4200                         /*
4201                          * Representation in the new basis.
4202                          */
4203
4204                         d2::point nbp0 = d2::point(bp0.dproduct(b0), bp0.dproduct(b1));
4205                         // d2::point nbp1 = d2::point(bp1.dproduct(b0), bp1.dproduct(b1));
4206                         d2::point nbp2 = d2::point(bp2.dproduct(b0), bp2.dproduct(b1));
4207                         d2::point np3  = d2::point( p3.dproduct(b0),  p3.dproduct(b1));
4208
4209                         /*
4210                          * Determine intersection of line
4211                          * (nbp0, nbp1), which is parallel to
4212                          * b0, with line (nbp2, np3).
4213                          */
4214
4215                         /*
4216                          * XXX: a stronger check would be
4217                          * better here, e.g., involving the
4218                          * ratio (np3[0] - nbp2[0]) / (np3[1] -
4219                          * nbp2[1]).  Also, acceptance of these
4220                          * cases is probably better than
4221                          * rejection.
4222                          */
4223
4224
4225                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4226
4227                         if (np3[1] - nbp2[1] == 0)
4228                                 continue;
4229
4230
4231                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4232
4233                         d2::point intersection = d2::point(nbp2[0]
4234                                         + (nbp0[1] - nbp2[1]) * (np3[0] - nbp2[0]) / (np3[1] - nbp2[1]),
4235                                         nbp0[1]);
4236
4237                         ale_pos b2_offset = b2.dproduct(bp0);
4238
4239                         /*
4240                          * Map the intersection back to the world
4241                          * basis.
4242                          */
4243
4244                         point iw = intersection[0] * b0 + intersection[1] * b1 + b2_offset * b2;
4245
4246                         /*
4247                          * Reject intersection points behind a
4248                          * camera.
4249                          */
4250
4251                         point icp = _pt1.wc(iw);
4252                         point ics = _pt2.wc(iw);
4253
4254
4255                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4256
4257                         if (icp[2] >= 0 || ics[2] >= 0)
4258                                 continue;
4259
4260
4261                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4262
4263                         /*
4264                          * Reject clipping plane violations.
4265                          */
4266
4267
4268                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4269
4270                         if (iw[2] > front_clip
4271                          || iw[2] < rear_clip)
4272                                 continue;
4273
4274
4275                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4276
4277                         /*
4278                          * Score the point.
4279                          */
4280
4281                         point ip = _pt1.wp_unscaled(iw);
4282
4283                         point is = _pt2.wp_unscaled(iw);
4284
4285                         analytic _a = { iw, ip, is };
4286
4287                         /*
4288                          * Calculate score from color match.  Assume for now
4289                          * that the transformation can be approximated locally
4290                          * with a translation.
4291                          */
4292
4293                         ale_pos score = 0;
4294                         ale_pos divisor = 0;
4295                         ale_pos l1_multiplier = 0.125;
4296
4297                         if (if1->in_bounds(ip.xy())
4298                          && if2->in_bounds(is.xy())
4299                          && !d2::render::is_excluded_f(ip.xy() * pdu_scale, f1)
4300                          && !d2::render::is_excluded_f(is.xy() * pdu_scale, f2)) {
4301                                 divisor += 1 - l1_multiplier;
4302                                 score += (1 - l1_multiplier)
4303                                        * (if1->get_bl(ip.xy()) - if2->get_bl(is.xy())).normsq();
4304                         }
4305
4306
4307                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4308
4309                         for (int iii = -1; iii <= 1; iii++)
4310                         for (int jjj = -1; jjj <= 1; jjj++) {
4311                                 d2::point t(iii, jjj);
4312
4313                                 if (!if1->in_bounds(ip.xy() + t)
4314                                  || !if2->in_bounds(is.xy() + t)
4315                                  || d2::render::is_excluded_f(ip.xy() * pdu_scale, f1)
4316                                  || d2::render::is_excluded_f(is.xy() * pdu_scale, f2))
4317                                         continue;
4318
4319                                 divisor += l1_multiplier;
4320                                 score   += l1_multiplier
4321                                          * (if1->get_bl(ip.xy() + t) - if2->get_bl(is.xy() + t)).normsq();
4322
4323                         }
4324
4325                         /*
4326                          * Include third-camera contributions in the score.
4327                          */
4328
4329                         if (tc_multiplier != 0)
4330                         for (unsigned int f = 0; f < d2::image_rw::count(); f++) {
4331                                 if (f == f1 || f == f2)
4332                                         continue;
4333
4334                                 const d2::image *if3 = al->get(f)->get_image(primary_decimation_upper);
4335                                 pt _pt3 = al->get(f)->get_t(primary_decimation_upper);
4336
4337                                 point p = _pt3.wp_unscaled(iw);
4338
4339                                 if (!if3->in_bounds(p.xy())
4340                                  || !if1->in_bounds(ip.xy())
4341                                  || d2::render::is_excluded_f(p.xy() * pdu_scale, f)
4342                                  || d2::render::is_excluded_f(ip.xy() * pdu_scale, f1))
4343                                         continue;
4344
4345                                 divisor += tc_multiplier;
4346                                 score   += tc_multiplier
4347                                          * (if1->get_bl(ip.xy()) - if3->get_bl(p.xy())).normsq();
4348                         }
4349
4350
4351
4352
4353                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4354
4355                         /*
4356                          * Reject points with undefined score.
4357                          */
4358
4359
4360                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4361
4362                         if (!finite(score / divisor))
4363                                 continue;
4364
4365
4366                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4367
4368 #if 0
4369                         /*
4370                          * XXX: reject points not on the z=-27.882252 plane.
4371                          */
4372
4373
4374                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4375
4376                         if (_a.ip[2] > -27 || _a.ip[2] < -28)
4377                                 continue;
4378 #endif
4379
4380
4381                         // fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);
4382
4383                         /*
4384                          * Add the point to the score map.
4385                          */
4386
4387 //                      d2::pixel c_ip = if1->in_bounds(ip.xy()) ? if1->get_bl(ip.xy())
4388 //                                                               : d2::pixel();
4389 //                      d2::pixel c_is = if2->in_bounds(is.xy()) ? if2->get_bl(is.xy())
4390 //                                                               : d2::pixel();
4391
4392 //                      fprintf(stderr, "Candidate subspace: f1=%u f2=%u i=%u j=%u ii=%f jj=%f"
4393 //                                      "cp=[%f %f %f] cs=[%f %f %f]\n",
4394 //                                      f1, f2, i, j, ii, jj, c_ip[0], c_ip[1], c_ip[2],
4395 //                                      c_is[0], c_is[1], c_is[2]);
4396
4397                         result.insert(score_map_element(score / divisor, _a));
4398                 }
4399
4400 //              fprintf(stderr, "Iterating through the score map:\n");
4401 //
4402 //              for (score_map::iterator smi = result.begin(); smi != result.end(); smi++) {
4403 //                      fprintf(stderr, "%f ", smi->first);
4404 //              }
4405 //
4406 //              fprintf(stderr, "\n");
4407
4408                 return result;
4409         }
4410
4411
4412         /*
4413          * Attempt to refine space around a point, to high and low resolutions
4414          * resulting in two resolutions in total.
4415          */
4416
4417         static space::traverse refine_space(point iw, ale_pos target_dim, int use_filler) {
4418
4419                 space::traverse st = space::traverse::root();
4420
4421                 if (!st.includes(iw)) {
4422                         assert(0);
4423                         return st;
4424                 }
4425
4426                 int lr_done = !use_filler;
4427
4428                 /*
4429                  * Loop until all resolutions of interest have been generated.
4430                  */
4431
4432                 for(;;) {
4433
4434                         point p[2] = { st.get_min(), st.get_max() };
4435
4436                         ale_pos dim_max = 0;
4437
4438                         for (int d = 0; d < 3; d++) {
4439                                 ale_pos d_value = fabs(p[0][d] - p[1][d]);
4440                                 if (d_value > dim_max)
4441                                         dim_max = d_value;
4442                         }
4443
4444                         /*
4445                          * Generate any new desired spatial registers.
4446                          */
4447
4448                         for (int f = 0; f < 2; f++) {
4449
4450                                 /*
4451                                  * Low resolution
4452                                  */
4453
4454                                 if (dim_max < 2 * target_dim
4455                                  && lr_done == 0) {
4456                                         if (spatial_info_map.find(st.get_node()) == spatial_info_map.end()) {
4457                                                 spatial_info_map[st.get_node()];
4458                                                 ui::get()->d3_increment_spaces();
4459                                         }
4460                                         lr_done = 1;
4461                                 }
4462
4463                                 /*
4464                                  * High resolution.
4465                                  */
4466
4467                                 if (dim_max < target_dim) {
4468                                         if (spatial_info_map.find(st.get_node()) == spatial_info_map.end()) {
4469                                                 spatial_info_map[st.get_node()];
4470                                                 ui::get()->d3_increment_spaces();
4471                                         }
4472                                         return st;
4473                                 }
4474                         }
4475
4476                         /*
4477                          * Check precision before analyzing space further.
4478                          */
4479
4480                         if (st.precision_wall()) {
4481                                 fprintf(stderr, "\n\n*** Error: reached subspace precision wall ***\n\n");
4482                                 assert(0);
4483                                 return st;
4484                         }
4485
4486                         if (st.positive().includes(iw)) {
4487                                 st = st.positive();
4488                                 total_tsteps++;
4489                         } else if (st.negative().includes(iw)) {
4490                                 st = st.negative();
4491                                 total_tsteps++;
4492                         } else {
4493                                 fprintf(stderr, "failed iw = (%f, %f, %f)\n",
4494                                                 iw[0], iw[1], iw[2]);
4495                                 assert(0);
4496                         }
4497                 }
4498         }
4499
4500         /*
4501          * Calculate target dimension
4502          */
4503
4504         static ale_pos calc_target_dim(point iw, pt _pt, const char *d_out[], const char *v_out[],
4505                         std::map<const char *, pt> *d3_depth_pt,
4506                         std::map<const char *, pt> *d3_output_pt) {
4507
4508                 ale_pos result = _pt.distance_1d(iw, primary_decimation_upper);
4509
4510                 for (unsigned int n = 0; n < d2::image_rw::count(); n++) {
4511                         if (d_out[n] && align::projective(n).distance_1d(iw, 0) < result)
4512                                 result = align::projective(n).distance_1d(iw, 0);
4513                         if (v_out[n] && align::projective(n).distance_1d(iw, 0) < result)
4514                                 result = align::projective(n).distance_1d(iw, 0);
4515                 }
4516
4517                 for (std::map<const char *, pt>::iterator i = d3_output_pt->begin(); i != d3_output_pt->end(); i++) {
4518                         if (i->second.distance_1d(iw, 0) < result)
4519                                 result = i->second.distance_1d(iw, 0);
4520                 }
4521
4522                 for (std::map<const char *, pt>::iterator i = d3_depth_pt->begin(); i != d3_depth_pt->end(); i++) {
4523                         if (i->second.distance_1d(iw, 0) < result)
4524                                 result = i->second.distance_1d(iw, 0);
4525                 }
4526
4527                 assert (result > 0);
4528
4529                 return result;
4530         }
4531
4532         /*
4533          * Calculate level of detail for a given viewpoint.
4534          */
4535
4536         static int calc_lod(ale_pos depth1, pt _pt, ale_pos target_dim) {
4537                 return (int) round(_pt.trilinear_coordinate(depth1, target_dim * sqrt(2)));
4538         }
4539
4540         /*
4541          * Calculate depth range for a given pair of viewpoints.
4542          */
4543
4544         static ale_pos calc_depth_range(point iw, pt _pt1, pt _pt2) {
4545
4546                 point ip = _pt1.wp_unscaled(iw);
4547
4548                 ale_pos reference_change = fabs(ip[2] / 1000);
4549
4550                 point iw1 = _pt1.pw_scaled(ip + point(0, 0, reference_change));
4551                 point iw2 = _pt1.pw_scaled(ip - point(0, 0, reference_change));
4552
4553                 point is = _pt2.wc(iw);
4554                 point is1 = _pt2.wc(iw1);
4555                 point is2 = _pt2.wc(iw2);
4556
4557                 assert(is[2] < 0);
4558
4559                 ale_pos d1 = (is1.xy() - is.xy()).norm();
4560                 ale_pos d2 = (is2.xy() - is.xy()).norm();
4561
4562                 // assert (reference_change > 0);
4563                 // assert (d1 > 0 || d2 > 0);
4564
4565                 if (is1[2] < 0 && is2[2] < 0) {
4566
4567                         if (d1 > d2)
4568                                 return reference_change / d1;
4569                         else
4570                                 return reference_change / d2;
4571                 }
4572
4573                 if (is1[2] < 0)
4574                         return reference_change / d1;
4575
4576                 if (is2[2] < 0)
4577                         return reference_change / d2;
4578
4579                 return 0;
4580         }
4581
4582         /*
4583          * Calculate a refined point for a given set of parameters.
4584          */
4585
4586         static point get_refined_point(pt _pt1, pt _pt2, int i, int j,
4587                         int f1, int f2, int lod1, int lod2, ale_pos depth,
4588                         ale_pos depth_range) {
4589
4590                 d2::pixel comparison_color = al->get(f1)->get_image(lod1)->get_pixel(i, j);
4591
4592                 ale_pos best = -1;
4593                 ale_pos best_depth = depth;
4594
4595                 assert (depth_range > 0);
4596
4597                 if (fabs(depth_range) < fabs(depth / 10000))
4598                         return _pt1.pw_unscaled(point(i, j, depth));
4599
4600                 for (ale_pos d = depth - depth_range; d < depth + depth_range; d += depth_range / 10) {
4601
4602                         if (!(d < 0))
4603                                 continue;
4604
4605                         point iw = _pt1.pw_unscaled(point(i, j, d));
4606                         point is = _pt2.wp_unscaled(iw);
4607
4608                         if (!(is[2] < 0))
4609                                 continue;
4610
4611                         if (!al->get(f2)->get_image(lod2)->in_bounds(is.xy()))
4612                                 continue;
4613
4614                         ale_pos error = (comparison_color - al->get(f2)->get_image(lod2)->get_bl(is.xy())).norm();
4615
4616                         if (error < best || best == -1) {
4617                                 best = error;
4618                                 best_depth = d;
4619                         }
4620                 }
4621
4622                 return _pt1.pw_unscaled(point(i, j, best_depth));
4623         }
4624
4625         /*
4626          * Analyze space in a manner dependent on the score map.
4627          */
4628
4629         static void analyze_space_from_map(const char *d_out[], const char *v_out[],
4630                                        std::map<const char *, pt> *d3_depth_pt,
4631                                        std::map<const char *, pt> *d3_output_pt,
4632                                        unsigned int f1, unsigned int f2,
4633                                        unsigned int i, unsigned int j, score_map _sm, int use_filler) {
4634
4635                 int accumulated_ambiguity = 0;
4636                 int max_acc_amb = pairwise_ambiguity;
4637
4638                 pt _pt1 = al->get(f1)->get_t(0);
4639                 pt _pt2 = al->get(f2)->get_t(0);
4640
4641                 if (_pt1.scale_2d() != 1)
4642                         use_filler = 1;
4643
4644                 for(score_map::iterator smi = _sm.begin(); smi != _sm.end(); smi++) {
4645
4646                         point iw = smi->second.iw;
4647
4648                         if (accumulated_ambiguity++ >= max_acc_amb)
4649                                 break;
4650
4651                         total_ambiguity++;
4652
4653                         ale_pos depth1 = _pt1.wc(iw)[2];
4654                         ale_pos depth2 = _pt2.wc(iw)[2];
4655
4656                         ale_pos target_dim = calc_target_dim(iw, _pt1, d_out, v_out, d3_depth_pt, d3_output_pt);
4657
4658                         assert(target_dim > 0);
4659
4660                         int lod1 = calc_lod(depth1, _pt1, target_dim);
4661                         int lod2 = calc_lod(depth2, _pt2, target_dim);
4662
4663                         while (lod1 < input_decimation_lower
4664                             || lod2 < input_decimation_lower) {
4665                                 target_dim *= 2;
4666                                 lod1 = calc_lod(depth1, _pt1, target_dim);
4667                                 lod2 = calc_lod(depth2, _pt2, target_dim);
4668                         }
4669
4670
4671                         if (lod1 >= (int) al->get(f1)->count()
4672                          || lod2 >= (int) al->get(f2)->count())
4673                                 continue;
4674
4675                         int multiplier = (unsigned int) floor(pow(2, primary_decimation_upper - lod1));
4676
4677                         ale_pos depth_range = calc_depth_range(iw, _pt1, _pt2);
4678
4679                         assert (depth_range > 0);
4680
4681                         pt _pt1_lod = al->get(f1)->get_t(lod1);
4682                         pt _pt2_lod = al->get(f2)->get_t(lod2);
4683
4684                         int im = i * multiplier;
4685                         int jm = j * multiplier;
4686
4687                         for (int ii = 0; ii < multiplier; ii += 1)
4688                         for (int jj = 0; jj < multiplier; jj += 1) {
4689
4690                                 point refined_point = get_refined_point(_pt1_lod, _pt2_lod, im + ii, jm + jj,
4691                                                 f1, f2, lod1, lod2, depth1, depth_range);
4692
4693                                 /*
4694                                  * Re-evaluate target dimension.
4695                                  */
4696
4697                                 ale_pos target_dim_ =
4698                                         calc_target_dim(refined_point, _pt1, d_out, v_out, d3_depth_pt, d3_output_pt);
4699
4700                                 ale_pos depth1_ = _pt1.wc(refined_point)[2];
4701                                 ale_pos depth2_ = _pt2.wc(refined_point)[2];
4702
4703                                 int lod1_ = calc_lod(depth1_, _pt1, target_dim_);
4704                                 int lod2_ = calc_lod(depth2_, _pt2, target_dim_);
4705
4706                                 while (lod1_ < input_decimation_lower
4707                                     || lod2_ < input_decimation_lower) {
4708                                         target_dim_ *= 2;
4709                                         lod1_ = calc_lod(depth1_, _pt1, target_dim_);
4710                                         lod2_ = calc_lod(depth2_, _pt2, target_dim_);
4711                                 }
4712
4713                                 /*
4714                                  * Attempt to refine space around the intersection point.
4715                                  */
4716
4717                                 space::traverse st =
4718                                         refine_space(refined_point, target_dim_, use_filler || _pt1.scale_2d() != 1);
4719
4720 //                              if (!resolution_ok(al->get(f1)->get_t(0), al->get(f1)->get_t(0).trilinear_coordinate(st))) {
4721 //                                      pt transformation = al->get(f1)->get_t(0);
4722 //                                      ale_pos tc = al->get(f1)->get_t(0).trilinear_coordinate(st);
4723 //
4724 //                                      fprintf(stderr, "Resolution not ok.\n");
4725 //                                      fprintf(stderr, "pow(2, tc)=%f\n", pow(2, tc));
4726 //                                      fprintf(stderr, "transformation.unscaled_height()=%f\n",
4727 //                                                      transformation.unscaled_height());
4728 //                                      fprintf(stderr, "transformation.unscaled_width()=%f\n",
4729 //                                                      transformation.unscaled_width());
4730 //                                      fprintf(stderr, "tc=%f", tc);
4731 //                                      fprintf(stderr, "input_decimation_lower - 1.5 = %f\n",
4732 //                                                      input_decimation_lower - 1.5);
4733 //
4734 //                              }
4735 //
4736 //                              assert(resolution_ok(al->get(f1)->get_t(0), al->get(f1)->get_t(0).trilinear_coordinate(st)));
4737 //                              assert(resolution_ok(al->get(f2)->get_t(0), al->get(f2)->get_t(0).trilinear_coordinate(st)));
4738                         }
4739
4740                 }
4741         }
4742
4743
4744         /*
4745          * Initialize space and identify regions of interest for the adaptive
4746          * subspace model.
4747          */
4748         static void make_space(const char *d_out[], const char *v_out[],
4749                         std::map<const char *, pt> *d3_depth_pt,
4750                         std::map<const char *, pt> *d3_output_pt) {
4751
4752                 ui::get()->d3_total_spaces(0);
4753
4754                 /*
4755                  * Variable indicating whether low-resolution filler space
4756                  * is desired to avoid aliased gaps in surfaces.
4757                  */
4758
4759                 int use_filler = d3_depth_pt->size() != 0
4760                                    || d3_output_pt->size() != 0
4761                                    || output_decimation_preferred > 0
4762                                    || input_decimation_lower > 0
4763                                    || !focus::is_trivial()
4764                                    || !strcmp(pairwise_comparisons, "all");
4765
4766                 std::vector<pt> pt_outputs = make_pt_list(d_out, v_out, d3_depth_pt, d3_output_pt);
4767
4768                 /*
4769                  * Initialize root space.
4770                  */
4771
4772                 space::init_root();
4773
4774                 /*
4775                  * Special handling for experimental option 'subspace_traverse'.
4776                  */
4777
4778                 if (subspace_traverse) {
4779                         /*
4780                          * Subdivide space to resolve intensity matches between pairs
4781                          * of frames.
4782                          */
4783
4784                         for (unsigned int f1 = 0; f1 < d2::image_rw::count(); f1++) {
4785
4786                                 if (d3_depth_pt->size() == 0
4787                                  && d3_output_pt->size() == 0
4788                                  && d_out[f1] == NULL
4789                                  && v_out[f1] == NULL)
4790                                         continue;
4791
4792                                 if (tc_multiplier == 0)
4793                                         al->open(f1);
4794
4795                                 for (unsigned int f2 = 0; f2 < d2::image_rw::count(); f2++) {
4796
4797                                         if (f1 == f2)
4798                                                 continue;
4799
4800                                         if (tc_multiplier == 0)
4801                                                 al->open(f2);
4802
4803                                         candidates *c = new candidates(f1);
4804
4805                                         find_candidates(f1, f2, c, point::neginf(), point::posinf(), pt_outputs);
4806
4807
4808
4809                                         c->generate_subspaces();
4810
4811                                         if (tc_multiplier == 0)
4812                                                 al->close(f2);
4813                                 }
4814
4815                                 if (tc_multiplier == 0)
4816                                         al->close(f1);
4817                         }
4818
4819                         return;
4820                 }
4821
4822                 /*
4823                  * Subdivide space to resolve intensity matches between pairs
4824                  * of frames.
4825                  */
4826
4827                 for (unsigned int f1 = 0; f1 < d2::image_rw::count(); f1++)
4828                 for (unsigned int f2 = 0; f2 < d2::image_rw::count(); f2++) {
4829                         if (f1 == f2)
4830                                 continue;
4831
4832                         if (!d_out[f1] && !v_out[f1] && !d3_depth_pt->size()
4833                          && !d3_output_pt->size() && strcmp(pairwise_comparisons, "all"))
4834                                 continue;
4835
4836                         if (tc_multiplier == 0) {
4837                                 al->open(f1);
4838                                 al->open(f2);
4839                         }
4840
4841                         /*
4842                          * Iterate over all points in the primary frame.
4843                          */
4844
4845                         ale_pos pdu_scale = pow(2, primary_decimation_upper);
4846
4847                         for (unsigned int i = 0; i < al->get(f1)->get_image(primary_decimation_upper)->height(); i++)
4848                         for (unsigned int j = 0; j < al->get(f1)->get_image(primary_decimation_upper)->width();  j++) {
4849
4850                                 if (d2::render::is_excluded_f(d2::point(i, j) * pdu_scale, f1))
4851                                         continue;
4852
4853                                 ui::get()->d3_subdivision_status(f1, f2, i, j);
4854
4855                                 total_pixels++;
4856
4857                                 /*
4858                                  * Generate a map from scores to 3D points for
4859                                  * various depths in f1.
4860                                  */
4861
4862                                 score_map _sm = p2f_score_map(f1, f2, i, j);
4863
4864                                 /*
4865                                  * Analyze space in a manner dependent on the score map.
4866                                  */
4867
4868                                 analyze_space_from_map(d_out, v_out, d3_depth_pt, d3_output_pt,
4869                                                 f1, f2, i, j, _sm, use_filler);
4870
4871                         }
4872
4873                         /*
4874                          * This ordering may encourage image f1 to be cached.
4875                          */
4876
4877                         if (tc_multiplier == 0) {
4878                                 al->close(f2);
4879                                 al->close(f1);
4880                         }
4881                 }
4882         }
4883
4884
4885         /*
4886          * Update spatial information structures.
4887          *
4888          * XXX: the name of this function is horribly misleading.  There isn't
4889          * even a 'search depth' any longer, since there is no longer any
4890          * bounded DFS occurring.
4891          */
4892         static void reduce_cost_to_search_depth(d2::exposure *exp_out, int inc_bit) {
4893
4894                 /*
4895                  * Subspace model
4896                  */
4897
4898                 ui::get()->set_steps(ou_iterations);
4899
4900                 for (unsigned int i = 0; i < ou_iterations; i++) {
4901                         ui::get()->set_steps_completed(i);
4902                         spatial_info_update();
4903                 }
4904
4905         }
4906
4907 #if 0
4908         /*
4909          * Describe a scene to a renderer
4910          */
4911         static void describe(render *r) {
4912         }
4913 #endif
4914 };
4915
4916 #endif