?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Return weight of 1 for just-out-of-bounds locations.
[Ale.git] / d3 / scene.h
blob3ae26a517da282d9c8a99e6aedfc86e20c1a0f46
1 // Copyright 2003, 2004, 2005 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 * Falloff exponent
78 */
79
80 static double falloff_exponent;
81
82 /*
83 * Third-camera error multiplier
84 */
85 static double tc_multiplier;
86
87 /*
88 * Occupancy update iterations
89 */
90 static unsigned int ou_iterations;
91
92 /*
93 * Pairwise ambiguity
94 */
95 static unsigned int pairwise_ambiguity;
96
97 /*
98 * Pairwise comparisons
99 */
100 static const char *pairwise_comparisons;
101
102 /*
103 * 3D Post-exclusion
104 */
105 static int d3px_count;
106 static double *d3px_parameters;
107
108 /*
109 * Nearness threshold
110 */
111 static const ale_real nearness;
112
113 /*
114 * Encounter threshold for defined pixels.
115 */
116 static double encounter_threshold;
117
118 /*
119 * Structure to hold input frame information at levels of
120 * detail between full detail and full decimation.
121 */
122 class lod_image {
123 unsigned int f;
124 unsigned int entries;
125 std::vector<const d2::image *> im;
126 std::vector<pt> transformation;
127
128 public:
129 /*
130 * Constructor
131 */
132 lod_image(unsigned int _f) {
133
134 pt _pt;
135
136 f = _f;
137 im.push_back(d2::image_rw::copy(f, "3D reference image"));
138 assert(im.back());
139 entries = 1;
140 _pt = d3::align::projective(f);
141 _pt.scale(1 / _pt.scale_2d());
142 transformation.push_back(_pt);
143
144 while (im.back()->height() > 4
145 && im.back()->width() > 4) {
146
147 im.push_back(im.back()->scale_by_half("3D, reduced LOD"));
148 assert(im.back());
149
150 _pt.scale(1 / _pt.scale_2d() / pow(2, entries));
151 transformation.push_back(_pt);
152
153 entries += 1;
154 }
155 }
156
157 /*
158 * Get the number of scales
159 */
160 unsigned int count() {
161 return entries;
162 }
163
164 /*
165 * Get an image.
166 */
167 const d2::image *get_image(unsigned int i) {
168 assert(i >= 0);
169 assert(i < entries);
170 return im[i];
171 }
172
173 int in_bounds(d2::point p) {
174 return im[0]->in_bounds(p);
175 }
176
177 /*
178 * Get a 'trilinear' color. We currently don't do interpolation
179 * between levels of detail; hence, it's discontinuous in tl_coord.
180 */
181 d2::pixel get_tl(d2::point p, ale_pos tl_coord) {
182
183 assert(in_bounds(p));
184
185 tl_coord = round(tl_coord);
186
187 if (tl_coord >= entries)
188 tl_coord = entries;
189 if (tl_coord < 0)
190 tl_coord = 0;
191
192 p = p / (ale_pos) pow(2, tl_coord);
193
194 unsigned int itlc = (unsigned int) tl_coord;
195
196 if (p[0] > im[itlc]->height() - 1)
197 p[0] = im[itlc]->height() - 1;
198 if (p[1] > im[itlc]->width() - 1)
199 p[1] = im[itlc]->width() - 1;
200
201 assert(p[0] >= 0);
202 assert(p[1] >= 0);
203
204 return im[itlc]->get_bl(p);
205 }
206
207 d2::pixel get_max_diff(d2::point p, ale_pos tl_coord) {
208 assert(in_bounds(p));
209
210 tl_coord = round(tl_coord);
211
212 if (tl_coord >= entries)
213 tl_coord = entries;
214 if (tl_coord < 0)
215 tl_coord = 0;
216
217 p = p / (ale_pos) pow(2, tl_coord);
218
219 unsigned int itlc = (unsigned int) tl_coord;
220
221 if (p[0] > im[itlc]->height() - 1)
222 p[0] = im[itlc]->height() - 1;
223 if (p[1] > im[itlc]->width() - 1)
224 p[1] = im[itlc]->width() - 1;
225
226 assert(p[0] >= 0);
227 assert(p[1] >= 0);
228
229 return im[itlc]->get_max_diff(p);
230 }
231
232 /*
233 * Get the transformation
234 */
235 pt get_t(unsigned int i) {
236 assert(i >= 0);
237 assert(i < entries);
238 return transformation[i];
239 }
240
241 /*
242 * Get the camera origin in world coordinates
243 */
244 point origin() {
245 return transformation[0].origin();
246 }
247
248 /*
249 * Destructor
250 */
251 ~lod_image() {
252 for (unsigned int i = 0; i < entries; i++) {
253 delete im[i];
254 }
255 }
256 };
257
258 /*
259 * Structure to hold weight information for reference images.
260 */
261 class ref_weights {
262 unsigned int f;
263 std::vector<d2::image *> weights;
264 pt transformation;
265 int tc_low;
266 int tc_high;
267 int initialized;
268
269 void set_image(d2::image *im, ale_real value) {
270 assert(im);
271 for (unsigned int i = 0; i < im->height(); i++)
272 for (unsigned int j = 0; j < im->width(); j++)
273 im->pix(i, j) = d2::pixel(value, value, value);
274 }
275
276 d2::image *make_image(ale_pos sf, ale_real init_value = 0) {
277 d2::image *result = new d2::image_ale_real(
278 (unsigned int) ceil(transformation.unscaled_height() * sf),
279 (unsigned int) ceil(transformation.unscaled_width() * sf), 3);
280 assert(result);
281
282 if (init_value != 0)
283 set_image(result, init_value);
284
285 fprintf(stderr, "make_image(%f)\n", sf);
286
287 if (initialized) {
288
289 fprintf(stderr, "Current image list: ");
290
291 for (int i = tc_low; i <= tc_high; i++) {
292 fprintf(stderr, "[tc=%d imax=%d jmax=%d] ",
293 i, weights[i - tc_low]->height(),
294 weights[i - tc_low]->width());
295 }
296
297 fprintf(stderr, "\n");
298
299 }
300
301 fprintf(stderr, "New image: [tc=%f imax=%d jmax=%d]\n",
302 sf, result->height(), result->width());
303
304
305 return result;
306 }
307
308 public:
309
310 /*
311 * Explicit weight subtree
312 */
313 struct subtree {
314 ale_real node_value;
315 subtree *children[2][2];
316
317 subtree(ale_real nv, subtree *a, subtree *b, subtree *c, subtree *d) {
318 node_value = nv;
319 children[0][0] = a;
320 children[0][1] = b;
321 children[1][0] = c;
322 children[1][1] = d;
323 }
324
325 ~subtree() {
326 for (int i = 0; i < 2; i++)
327 for (int j = 0; j < 2; j++)
328 delete children[i][j];
329 }
330 };
331
332 /*
333 * Constructor
334 */
335 ref_weights(unsigned int _f) {
336 f = _f;
337 transformation = d3::align::projective(f);
338 initialized = 0;
339 }
340
341 /*
342 * Check spatial bounds.
343 */
344 int in_spatial_bounds(point p) {
345
346 if (!p.defined())
347 return 0;
348
349 if (p[0] < 0)
350 return 0;
351 if (p[1] < 0)
352 return 0;
353 if (p[0] > transformation.unscaled_height() - 1)
354 return 0;
355 if (p[1] > transformation.unscaled_width() - 1)
356 return 0;
357 if (p[2] >= 0)
358 return 0;
359
360 return 1;
361 }
362
363 int in_spatial_bounds(const space::traverse &t) {
364 point p = transformation.centroid(t);
365 return in_spatial_bounds(p);
366 }
367
368 /*
369 * Increase resolution to the given level.
370 */
371 void increase_resolution(int tc, unsigned int i, unsigned int j) {
372 d2::image *im = weights[tc - tc_low];
373 assert(im);
374 assert(i <= im->height() - 1);
375 assert(j <= im->width() - 1);
376
377 /*
378 * Check for the cases known to have no lower level of detail.
379 */
380
381 if (im->pix(i, j)[0] == -1)
382 return;
383
384 if (tc == tc_high)
385 return;
386
387 increase_resolution(tc + 1, i / 2, j / 2);
388
389 /*
390 * Load the lower-level image structure.
391 */
392
393 d2::image *iim = weights[tc + 1 - tc_low];
394
395 assert(iim);
396 assert(i / 2 <= iim->height() - 1);
397 assert(j / 2 <= iim->width() - 1);
398
399 /*
400 * Check for the case where no lower level of detail is
401 * available.
402 */
403
404 if (iim->pix(i / 2, j / 2)[0] == -1)
405 return;
406
407 /*
408 * Spread out the lower level of detail among (uninitialized)
409 * peer values.
410 */
411
412 for (unsigned int ii = (i / 2) * 2; ii < (i / 2) * 2 + 1; ii++)
413 for (unsigned int jj = (j / 2) * 2; jj < (j / 2) * 2 + 1; jj++) {
414 assert(ii <= im->height() - 1);
415 assert(jj <= im->width() - 1);
416 assert(im->pix(ii, jj)[0] == 0);
417
418 im->pix(ii, jj) = iim->pix(i / 2, j / 2);
419 }
420
421 /*
422 * Set the lower level of detail to point here.
423 */
424
425 iim->pix(i / 2, j / 2)[0] = -1;
426 }
427
428 /*
429 * Add weights to positive higher-resolution pixels where
430 * found when their current values match the given subtree
431 * values; set negative pixels to zero and return 0 if no
432 * positive higher-resolution pixels are found.
433 */
434 int add_partial(int tc, unsigned int i, unsigned int j, ale_real weight, subtree *st) {
435 d2::image *im = weights[tc - tc_low];
436 assert(im);
437
438 fprintf(stderr, "[ap tc=%d i=%d j=%d imax=%d jmax=%d]\n",
439 tc, i, j, im->height(), im->width());
440
441 if (i == im->height() - 1
442 || j == im->width() - 1) {
443 return 1;
444 }
445
446 assert(i <= im->height() - 1);
447 assert(j <= im->width() - 1);
448
449 /*
450 * Check for positive values.
451 */
452
453 if (im->pix(i, j)[0] > 0) {
454 if (st && st->node_value == im->pix(i, j)[0])
455 im->pix(i, j)[0] += weight;
456 return 1;
457 }
458
459 /*
460 * Handle the case where there are no higher levels of detail.
461 */
462
463 if (tc == tc_low) {
464 assert(im->pix(i, j)[0] == 0);
465 return 0;
466 }
467
468 /*
469 * Handle the case where higher levels of detail are available.
470 */
471
472 int success[2][2];
473
474 for (int ii = 0; ii < 2; ii++)
475 for (int jj = 0; jj < 2; jj++)
476 success[ii][jj] = add_partial(tc - 1, i * 2 + ii, j * 2 + jj, weight,
477 st ? st->children[ii][jj] : NULL);
478
479 if (!success[0][0]
480 && !success[0][1]
481 && !success[1][0]
482 && !success[1][1]) {
483 im->pix(i, j)[0] = 0;
484 return 0;
485 }
486
487 d2::image *iim = weights[tc - 1 - tc_low];
488 assert(iim);
489
490 for (int ii = 0; ii < 2; ii++)
491 for (int jj = 0; jj < 2; jj++)
492 if (success[ii][jj] == 0) {
493 assert(i * 2 + ii < iim->height());
494 assert(j * 2 + jj < iim->width());
495
496 im->pix(i * 2 + ii, j * 2 + jj)[0] = weight;
497 }
498
499 im->pix(i, j)[0] = -1;
500
501 return 1;
502 }
503
504 /*
505 * Add weight.
506 */
507 void add_weight(int tc, unsigned int i, unsigned int j, ale_real weight, subtree *st) {
508
509 d2::image *im = weights[tc - tc_low];
510 assert(im);
511
512 fprintf(stderr, "[aw tc=%d i=%d j=%d imax=%d jmax=%d]\n",
513 tc, i, j, im->height(), im->width());
514
515 assert(i <= im->height() - 1);
516 assert(j <= im->width() - 1);
517
518 /*
519 * Increase resolution, if necessary
520 */
521
522 increase_resolution(tc, i, j);
523
524 /*
525 * Attempt to add the weight at levels of detail
526 * where weight is defined.
527 */
528
529 if (add_partial(tc, i, j, weight, st))
530 return;
531
532 /*
533 * If no weights are defined at any level of detail,
534 * then set the weight here.
535 */
536
537 im->pix(i, j)[0] = weight;
538 }
539
540 void add_weight(int tc, d2::point p, ale_real weight, subtree *st) {
541
542 p *= pow(2, -tc);
543
544 unsigned int i = (unsigned int) floor(p[0]);
545 unsigned int j = (unsigned int) floor(p[1]);
546
547 add_weight(tc, i, j, weight, st);
548 }
549
550 void add_weight(const space::traverse &t, ale_real weight, subtree *st) {
551
552 if (weight == 0)
553 return;
554
555 ale_pos tc = transformation.trilinear_coordinate(t);
556 point p = transformation.centroid(t);
557 assert(in_spatial_bounds(p));
558
559 tc = round(tc);
560
561 /*
562 * Establish a reasonable (?) upper bound on resolution.
563 */
564
565 if (tc < input_decimation_lower) {
566 weight /= pow(4, (input_decimation_lower - tc));
567 tc = input_decimation_lower;
568 }
569
570 /*
571 * Initialize, if necessary.
572 */
573
574 if (!initialized) {
575 tc_low = tc_high = (int) tc;
576
577 ale_real sf = pow(2, -tc);
578
579 weights.push_back(make_image(sf));
580
581 initialized = 1;
582 }
583
584 /*
585 * Check resolution bounds
586 */
587
588 assert (tc_low <= tc_high);
589
590 /*
591 * Generate additional levels of detail, if necessary.
592 */
593
594 while (tc < tc_low) {
595 tc_low--;
596
597 ale_real sf = pow(2, -tc_low);
598
599 weights.insert(weights.begin(), make_image(sf));
600 }
601
602 while (tc > tc_high) {
603 tc_high++;
604
605 ale_real sf = pow(2, -tc_high);
606
607 weights.push_back(make_image(sf, -1));
608 }
609
610 add_weight((int) tc, p.xy(), weight, st);
611 }
612
613 /*
614 * Get weight
615 */
616 ale_real get_weight(int tc, unsigned int i, unsigned int j) {
617
618 fprintf(stderr, "[gw tc=%d i=%u j=%u tclow=%d tchigh=%d]\n",
619 tc, i, j, tc_low, tc_high);
620
621 if (tc < tc_low || !initialized)
622 return 0;
623
624 if (tc > tc_high) {
625 return (get_weight(tc - 1, i * 2 + 0, j * 2 + 0)
626 + get_weight(tc - 1, i * 2 + 1, j * 2 + 0)
627 + get_weight(tc - 1, i * 2 + 1, j * 2 + 1)
628 + get_weight(tc - 1, i * 2 + 0, j * 2 + 1)) / 4;
629 }
630
631 assert(weights.size() > (unsigned int) (tc - tc_low));
632
633 d2::image *im = weights[tc - tc_low];
634 assert(im);
635 assert(i < im->height());
636 assert(j < im->width());
637
638 if (im->pix(i, j)[0] == -1) {
639 return (get_weight(tc - 1, i * 2 + 0, j * 2 + 0)
640 + get_weight(tc - 1, i * 2 + 1, j * 2 + 0)
641 + get_weight(tc - 1, i * 2 + 1, j * 2 + 1)
642 + get_weight(tc - 1, i * 2 + 0, j * 2 + 1)) / 4;
643 }
644
645 if (im->pix(i, j)[0] == 0) {
646 if (tc == tc_high)
647 return 0;
648 if (weights[tc - tc_low + 1]->pix(i / 2, j / 2)[0] == -1)
649 return 0;
650 return get_weight(tc + 1, i / 2, j / 2);
651 }
652
653 return im->pix(i, j)[0];
654 }
655
656 ale_real get_weight(int tc, d2::point p) {
657
658 p *= pow(2, -tc);
659
660 unsigned int i = (unsigned int) floor(p[0]);
661 unsigned int j = (unsigned int) floor(p[1]);
662
663 return get_weight(tc, i, j);
664 }
665
666 ale_real get_weight(const space::traverse &t) {
667 ale_pos tc = transformation.trilinear_coordinate(t);
668 point p = transformation.centroid(t);
669 assert(in_spatial_bounds(p));
670
671 if (!initialized)
672 return 0;
673
674 tc = round(tc);
675
676 if (tc < tc_low) {
677 tc = tc_low;
678 }
679
680 if (tc <= tc_high) {
681 return get_weight((int) tc, p.xy());
682 }
683
684 assert(tc > tc_high);
685
686 int multiplier = (int) pow(2, (tc - tc_high));
687 ale_real result = 0;
688
689 for (int i = 0; i < multiplier; i++)
690 for (int j = 0; j < multiplier; j++) {
691 result += get_weight(tc_high,
692 (unsigned int) floor(p[0]) * multiplier + i,
693 (unsigned int) floor(p[1]) * multiplier + j);
694 }
695
696 return result;
697 }
698
699 /*
700 * Check whether a subtree is simple.
701 */
702 int is_simple(subtree *s) {
703 assert (s);
704
705 if (s->node_value == -1
706 && s->children[0][0] == NULL
707 && s->children[0][1] == NULL
708 && s->children[1][0] == NULL
709 && s->children[1][1] == NULL)
710 return 1;
711
712 return 0;
713 }
714
715 /*
716 * Get a weight subtree.
717 */
718 subtree *get_subtree(int tc, unsigned int i, unsigned int j) {
719
720 /*
721 * tc > tc_high is handled recursively.
722 */
723
724 if (tc > tc_high) {
725 subtree *result = new subtree(-1,
726 get_subtree(tc - 1, i * 2 + 0, j * 2 + 0),
727 get_subtree(tc - 1, i * 2 + 0, j * 2 + 1),
728 get_subtree(tc - 1, i * 2 + 1, j * 2 + 0),
729 get_subtree(tc - 1, i * 2 + 1, j * 2 + 1));
730
731 if (is_simple(result)) {
732 delete result;
733 return NULL;
734 }
735
736 return result;
737 }
738
739 assert(tc >= tc_low);
740 assert(weights[tc - tc_low]);
741
742 d2::image *im = weights[tc - tc_low];
743
744 /*
745 * Rectangular images will, in general, have
746 * out-of-bounds tree sections. Handle this case.
747 */
748
749 if (i >= im->height())
750 return NULL;
751 if (j >= im->width())
752 return NULL;
753
754 /*
755 * -1 weights are handled recursively
756 */
757
758 if (im->pix(i, j)[0] == -1) {
759 subtree *result = new subtree(-1,
760 get_subtree(tc - 1, i * 2 + 0, j * 2 + 0),
761 get_subtree(tc - 1, i * 2 + 0, j * 2 + 1),
762 get_subtree(tc - 1, i * 2 + 1, j * 2 + 0),
763 get_subtree(tc - 1, i * 2 + 1, j * 2 + 1));
764
765 if (is_simple(result)) {
766 im->pix(i, j)[0] = 0;
767 delete result;
768 return NULL;
769 }
770
771 return result;
772 }
773
774 /*
775 * Zero weights have NULL subtrees.
776 */
777
778 if (im->pix(i, j)[0] == 0)
779 return NULL;
780
781 /*
782 * Handle the remaining case.
783 */
784
785 return new subtree(im->pix(i, j)[0], NULL, NULL, NULL, NULL);
786 }
787
788 subtree *get_subtree(int tc, d2::point p) {
789 p *= pow(2, -tc);
790
791 unsigned int i = (unsigned int) floor(p[0]);
792 unsigned int j = (unsigned int) floor(p[1]);
793
794 return get_subtree(tc, i, j);
795 }
796
797 subtree *get_subtree(const space::traverse &t) {
798 ale_pos tc = transformation.trilinear_coordinate(t);
799 point p = transformation.centroid(t);
800 assert(in_spatial_bounds(p));
801
802 if (!initialized)
803 return NULL;
804
805 if (tc < input_decimation_lower)
806 tc = input_decimation_lower;
807
808 tc = round(tc);
809
810 if (tc < tc_low)
811 return NULL;
812
813 return get_subtree((int) tc, p.xy());
814 }
815
816 /*
817 * Destructor
818 */
819 ~ref_weights() {
820 for (unsigned int i = 0; i < weights.size(); i++) {
821 delete weights[i];
822 }
823 }
824 };
825
826 /*
827 * Resolution check.
828 */
829 static int resolution_ok(pt transformation, ale_pos tc) {
830
831 if (pow(2, tc) > transformation.unscaled_height()
832 || pow(2, tc) > transformation.unscaled_width())
833 return 0;
834
835 if (tc < input_decimation_lower)
836 return 0;
837
838 return 1;
839 }
840
841 /*
842 * Structure to hold input frame information at all levels of detail.
843 */
844 class lod_images {
845
846 /*
847 * All images.
848 */
849
850 std::vector<lod_image *> images;
851
852 public:
853
854 lod_images() {
855 images.resize(d2::image_rw::count(), NULL);
856 }
857
858 void open(unsigned int f) {
859 assert (images[f] == NULL);
860
861 if (images[f] == NULL)
862 images[f] = new lod_image(f);
863 }
864
865 void open_all() {
866 for (unsigned int f = 0; f < d2::image_rw::count(); f++)
867 open(f);
868 }
869
870 lod_image *get(unsigned int f) {
871 assert (images[f] != NULL);
872 return images[f];
873 }
874
875 void close(unsigned int f) {
876 assert (images[f] != NULL);
877 delete images[f];
878 images[f] = NULL;
879 }
880
881 void close_all() {
882 for (unsigned int f = 0; f < d2::image_rw::count(); f++)
883 close(f);
884 }
885
886 ~lod_images() {
887 close_all();
888 }
889 };
890
891 /*
892 * All levels-of-detail
893 */
894
895 static struct lod_images *al;
896
897 /*
898 * Data structure for storing best encountered subspace candidates.
899 */
900 class candidates {
901 std::vector<std::vector<std::pair<ale_pos, ale_real> > > levels;
902 int image_index;
903 unsigned int height;
904 unsigned int width;
905
906 /*
907 * Point p is in world coordinates.
908 */
909 void generate_subspace(point iw, ale_pos diagonal) {
910
911 fprintf(stderr, "[gs iw=%f %f %f d=%f]\n",
912 iw[0], iw[1], iw[2], diagonal);
913
914 fprintf(stderr, "*");
915
916 space::traverse st = space::traverse::root();
917
918 if (!st.includes(iw)) {
919 assert(0);
920 return;
921 }
922
923 int highres = 0;
924 int lowres = 0;
925
926 /*
927 * Loop until resolutions of interest have been generated.
928 */
929
930 for(;;) {
931
932 ale_pos current_diagonal = (st.get_max() - st.get_min()).norm();
933
934 assert(!isnan(current_diagonal));
935
936 /*
937 * Generate any new desired spatial registers.
938 */
939
940 /*
941 * Inputs
942 */
943
944 for (int f = 0; f < 2; f++) {
945
946 /*
947 * Low resolution
948 */
949
950 if (current_diagonal < 2 * diagonal
951 && lowres == 0) {
952 spatial_info_map[st.get_node()];
953 lowres = 1;
954 }
955
956 /*
957 * High resolution.
958 */
959
960 if (current_diagonal < diagonal
961 && highres == 0) {
962 spatial_info_map[st.get_node()];
963 highres = 1;
964 }
965 }
966
967 /*
968 * Check for completion
969 */
970
971 if (highres && lowres)
972 return;
973
974 /*
975 * Check precision before analyzing space further.
976 */
977
978 if (st.precision_wall()) {
979 fprintf(stderr, "\n\n*** Error: reached subspace precision wall ***\n\n");
980 assert(0);
981 return;
982 }
983
984 if (st.positive().includes(iw)) {
985 st = st.positive();
986 total_tsteps++;
987 } else if (st.negative().includes(iw)) {
988 st = st.negative();
989 total_tsteps++;
990 } else {
991 fprintf(stderr, "failed iw = (%f, %f, %f)\n",
992 iw[0], iw[1], iw[2]);
993 assert(0);
994 }
995 }
996 }
997
998 public:
999 candidates(int f) {
1000
1001 image_index = f;
1002 height = (unsigned int) al->get(f)->get_t(0).unscaled_height();
1003 width = (unsigned int) al->get(f)->get_t(0).unscaled_width();
1004
1005 /*
1006 * Is this necessary?
1007 */
1008
1009 levels.resize(primary_decimation_upper - input_decimation_lower + 1);
1010 for (int l = input_decimation_lower; l <= primary_decimation_upper; l++) {
1011 levels[l - input_decimation_lower].resize((unsigned int) (floor(height / pow(2, l))
1012 * floor(width / pow(2, l))
1013 * pairwise_ambiguity),
1014 std::pair<ale_pos, ale_real>(0, 0));
1015 }
1016 }
1017
1018 /*
1019 * Point p is expected to be in local projective coordinates.
1020 */
1021
1022 void add_candidate(point p, int tc, ale_real score) {
1023 assert(tc <= primary_decimation_upper);
1024 assert(tc >= input_decimation_lower);
1025 assert(p[2] < 0);
1026 assert(score >= 0);
1027
1028 int i = (unsigned int) floor(p[0] / pow(2, tc));
1029 int j = (unsigned int) floor(p[1] / pow(2, tc));
1030
1031 int sheight = (int) floor(height / pow(2, tc));
1032 int swidth = (int) floor(width / pow(2, tc));
1033
1034 assert(i < sheight);
1035 assert(j < swidth);
1036
1037 for (unsigned int k = 0; k < pairwise_ambiguity; k++) {
1038 std::pair<ale_pos, ale_real> *pk =
1039 &(levels[tc - input_decimation_lower][i * swidth * pairwise_ambiguity + j * pairwise_ambiguity + k]);
1040
1041 if (pk->first != 0 && score >= pk->second)
1042 continue;
1043
1044 if (i == 1 && j == 1 && tc == 4)
1045 fprintf(stderr, "[ac p2=%f score=%f]\n", p[2], score);
1046
1047 ale_pos tp = pk->first;
1048 ale_real tr = pk->second;
1049
1050 pk->first = p[2];
1051 pk->second = score;
1052
1053 p[2] = tp;
1054 score = tr;
1055
1056 if (p[2] == 0)
1057 break;
1058 }
1059 }
1060
1061 /*
1062 * Generate subspaces for candidates.
1063 */
1064
1065 void generate_subspaces() {
1066
1067 fprintf(stderr, "+");
1068 for (int l = input_decimation_lower; l <= primary_decimation_upper; l++) {
1069 unsigned int sheight = (unsigned int) floor(height / pow(2, l));
1070 unsigned int swidth = (unsigned int) floor(width / pow(2, l));
1071
1072 for (unsigned int i = 0; i < sheight; i++)
1073 for (unsigned int j = 0; j < swidth; j++)
1074 for (unsigned int k = 0; k < pairwise_ambiguity; k++) {
1075 std::pair<ale_pos, ale_real> *pk =
1076 &(levels[l - input_decimation_lower]
1077 [i * swidth * pairwise_ambiguity + j * pairwise_ambiguity + k]);
1078
1079 if (pk->first == 0) {
1080 fprintf(stderr, "o");
1081 continue;
1082 } else {
1083 fprintf(stderr, "|");
1084 }
1085
1086 ale_pos si = i * pow(2, l) + ((l > 0) ? pow(2, l - 1) : 0);
1087 ale_pos sj = j * pow(2, l) + ((l > 0) ? pow(2, l - 1) : 0);
1088
1089 fprintf(stderr, "[gss l=%d i=%d j=%d d=%g]\n", l, i, j, pk->first);
1090
1091 point p = al->get(image_index)->get_t(0).pw_unscaled(point(si, sj, pk->first));
1092
1093 generate_subspace(p,
1094 al->get(image_index)->get_t(0).diagonal_distance_3d(pk->first, l));
1095 }
1096 }
1097 }
1098 };
1099
1100 /*
1101 * List for calculating weighted median.
1102 */
1103 class wml {
1104 ale_real *data;
1105 unsigned int size;
1106 unsigned int used;
1107
1108 ale_real &_w(unsigned int i) {
1109 assert(i <= used);
1110 return data[i * 2];
1111 }
1112
1113 ale_real &_d(unsigned int i) {
1114 assert(i <= used);
1115 return data[i * 2 + 1];
1116 }
1117
1118 void increase_capacity() {
1119
1120 if (size > 0)
1121 size *= 2;
1122 else
1123 size = 1;
1124
1125 data = (ale_real *) realloc(data, sizeof(ale_real) * 2 * (size * 1));
1126
1127 assert(data);
1128 assert (size > used);
1129
1130 if (!data) {
1131 fprintf(stderr, "Unable to allocate %d bytes of memory\n",
1132 sizeof(ale_real) * 2 * (size * 1));
1133 exit(1);
1134 }
1135 }
1136
1137 void insert_weight(unsigned int i, ale_real v, ale_real w) {
1138 assert(used < size);
1139 assert(used >= i);
1140 for (unsigned int j = used; j > i; j--) {
1141 _w(j) = _w(j - 1);
1142 _d(j) = _d(j - 1);
1143 }
1144
1145 _w(i) = w;
1146 _d(i) = v;
1147
1148 used++;
1149 }
1150
1151 public:
1152
1153 unsigned int get_size() {
1154 return size;
1155 }
1156
1157 unsigned int get_used() {
1158 return used;
1159 }
1160
1161 void print_info() {
1162 fprintf(stderr, "[st %p sz %u el", this, size);
1163 for (unsigned int i = 0; i < used; i++)
1164 fprintf(stderr, " (%f, %f)", _d(i), _w(i));
1165 fprintf(stderr, "]\n");
1166 }
1167
1168 void clear() {
1169 used = 0;
1170 }
1171
1172 void insert_weight(ale_real v, ale_real w) {
1173 for (unsigned int i = 0; i < used; i++) {
1174 if (_d(i) == v) {
1175 _w(i) += w;
1176 return;
1177 }
1178 if (_d(i) > v) {
1179 if (used == size)
1180 increase_capacity();
1181 insert_weight(i, v, w);
1182 return;
1183 }
1184 }
1185 if (used == size)
1186 increase_capacity();
1187 insert_weight(used, v, w);
1188 }
1189
1190 /*
1191 * Finds the median at half-weight, or between half-weight
1192 * and zero-weight, depending on the attenuation value.
1193 */
1194
1195 ale_real find_median(double attenuation = 0) {
1196
1197 assert(attenuation >= 0);
1198 assert(attenuation <= 1);
1199
1200 ale_real zero1 = 0;
1201 ale_real zero2 = 0;
1202 ale_real undefined = zero1 / zero2;
1203
1204 ale_accum weight_sum = 0;
1205
1206 for (unsigned int i = 0; i < used; i++)
1207 weight_sum += _w(i);
1208
1209 // if (weight_sum == 0)
1210 // return undefined;
1211
1212 if (used == 0 || used == 1)
1213 return undefined;
1214
1215 if (weight_sum == 0) {
1216 ale_accum data_sum = 0;
1217 for (unsigned int i = 0; i < used; i++)
1218 data_sum += _d(i);
1219 return data_sum / used;
1220 }
1221
1222
1223 ale_accum midpoint = weight_sum * (0.5 - 0.5 * attenuation);
1224
1225 ale_accum weight_sum_2 = 0;
1226
1227 for (unsigned int i = 0; i < used && weight_sum_2 < midpoint; i++) {
1228 weight_sum_2 += _w(i);
1229
1230 if (weight_sum_2 > midpoint) {
1231 return _d(i);
1232 } else if (weight_sum_2 == midpoint) {
1233 assert (i + 1 < used);
1234 return (_d(i) + _d(i + 1)) / 2;
1235 }
1236 }
1237
1238 return undefined;
1239 }
1240
1241 wml(int initial_size = 0) {
1242
1243 // if (initial_size == 0) {
1244 // initial_size = (int) (d2::image_rw::count() * 1.5);
1245 // }
1246
1247 size = initial_size;
1248 used = 0;
1249
1250 if (size > 0) {
1251 data = (ale_real *) malloc(size * sizeof(ale_real) * 2);
1252 assert(data);
1253 } else {
1254 data = NULL;
1255 }
1256 }
1257
1258 /*
1259 * copy constructor. This is required to avoid undesired frees.
1260 */
1261
1262 wml(const wml &w) {
1263 size = w.size;
1264 used = w.used;
1265 data = (ale_real *) malloc(size * sizeof(ale_real) * 2);
1266 assert(data);
1267
1268 memcpy(data, w.data, size * sizeof(ale_real) * 2);
1269 }
1270
1271 ~wml() {
1272 free(data);
1273 }
1274 };
1275
1276 /*
1277 * Class for information regarding spatial regions of interest.
1278 *
1279 * This class is configured for convenience in cases where sampling is
1280 * performed using an approximation of the fine:box:1,triangle:2 chain.
1281 * In this case, the *_1 variables would store the fine data and the
1282 * *_2 variables would store the coarse data. Other uses are also
1283 * possible.
1284 */
1285
1286 class spatial_info {
1287 /*
1288 * Map channel value --> weight.
1289 */
1290 wml color_weights_1[3];
1291 wml color_weights_2[3];
1292
1293 /*
1294 * Current color.
1295 */
1296 d2::pixel color;
1297
1298 /*
1299 * Map occupancy value --> weight.
1300 */
1301 wml occupancy_weights_1;
1302 wml occupancy_weights_2;
1303
1304 /*
1305 * Current occupancy value.
1306 */
1307 ale_real occupancy;
1308
1309 /*
1310 * pocc/socc density
1311 */
1312
1313 unsigned int pocc_density;
1314 unsigned int socc_density;
1315
1316 /*
1317 * Insert a weight into a list.
1318 */
1319 void insert_weight(wml *m, ale_real v, ale_real w) {
1320 m->insert_weight(v, w);
1321 }
1322
1323 /*
1324 * Find the median of a weighted list. Uses NaN for undefined.
1325 */
1326 ale_real find_median(wml *m, double attenuation = 0) {
1327 return m->find_median(attenuation);
1328 }
1329
1330 public:
1331 /*
1332 * Constructor.
1333 */
1334 spatial_info() {
1335 color = d2::pixel::zero();
1336 occupancy = 0;
1337 pocc_density = 0;
1338 socc_density = 0;
1339 }
1340
1341 /*
1342 * Accumulate color; primary data set.
1343 */
1344 void accumulate_color_1(int f, d2::pixel color, d2::pixel weight) {
1345 for (int k = 0; k < 3; k++)
1346 insert_weight(&color_weights_1[k], color[k], weight[k]);
1347 }
1348
1349 /*
1350 * Accumulate color; secondary data set.
1351 */
1352 void accumulate_color_2(d2::pixel color, d2::pixel weight) {
1353 for (int k = 0; k < 3; k++)
1354 insert_weight(&color_weights_2[k], color[k], weight[k]);
1355 }
1356
1357 /*
1358 * Accumulate occupancy; primary data set.
1359 */
1360 void accumulate_occupancy_1(int f, ale_real occupancy, ale_real weight) {
1361 insert_weight(&occupancy_weights_1, occupancy, weight);
1362 }
1363
1364 /*
1365 * Accumulate occupancy; secondary data set.
1366 */
1367 void accumulate_occupancy_2(ale_real occupancy, ale_real weight) {
1368 insert_weight(&occupancy_weights_2, occupancy, weight);
1369
1370 if (occupancy == 0 || occupancy_weights_2.get_size() < 96)
1371 return;
1372
1373 // fprintf(stderr, "%p updated socc with: %f %f\n", this, occupancy, weight);
1374 // occupancy_weights_2.print_info();
1375 }
1376
1377 /*
1378 * Update color (and clear accumulation structures).
1379 */
1380 void update_color() {
1381 for (int d = 0; d < 3; d++) {
1382 ale_real c = find_median(&color_weights_1[d]);
1383 if (isnan(c))
1384 c = find_median(&color_weights_2[d]);
1385 if (isnan(c))
1386 c = 0;
1387
1388 color[d] = c;
1389
1390 color_weights_1[d].clear();
1391 color_weights_2[d].clear();
1392 }
1393 }
1394
1395 /*
1396 * Update occupancy (and clear accumulation structures).
1397 */
1398 void update_occupancy() {
1399 ale_real o = find_median(&occupancy_weights_1, occ_att);
1400 if (isnan(o))
1401 o = find_median(&occupancy_weights_2, occ_att);
1402 if (isnan(o))
1403 o = 0;
1404
1405 occupancy = o;
1406
1407 pocc_density = occupancy_weights_1.get_used();
1408 socc_density = occupancy_weights_2.get_used();
1409
1410 occupancy_weights_1.clear();
1411 occupancy_weights_2.clear();
1412
1413 }
1414
1415 /*
1416 * Get current color.
1417 */
1418 d2::pixel get_color() {
1419 return color;
1420 }
1421
1422 /*
1423 * Get current occupancy.
1424 */
1425 ale_real get_occupancy() {
1426 return occupancy;
1427 }
1428
1429 /*
1430 * Get primary color density.
1431 */
1432
1433 unsigned int get_pocc_density() {
1434 return pocc_density;
1435 }
1436
1437 unsigned int get_socc_density() {
1438 return socc_density;
1439 }
1440 };
1441
1442 /*
1443 * Map spatial regions of interest to spatial info structures. XXX:
1444 * This may get very poor cache behavior in comparison with, say, an
1445 * array. Unfortunately, there is no immediately obvious array
1446 * representation. If some kind of array representation were adopted,
1447 * it would probably cluster regions of similar depth from the
1448 * perspective of the typical camera. In particular, for a
1449 * stereoscopic view, depth ordering for two random points tends to be
1450 * similar between cameras, I think. Unfortunately, it is never
1451 * identical for all points (unless cameras are co-located). One
1452 * possible approach would be to order based on, say, camera 0's idea
1453 * of depth.
1454 */
1455
1456 typedef std::map<struct space::node *, spatial_info> spatial_info_map_t;
1457
1458 static spatial_info_map_t spatial_info_map;
1459
1460 public:
1461
1462 /*
1463 * Debugging variables.
1464 */
1465
1466 static unsigned long total_ambiguity;
1467 static unsigned long total_pixels;
1468 static unsigned long total_divisions;
1469 static unsigned long total_tsteps;
1470
1471 /*
1472 * Member functions
1473 */
1474
1475 static void et(double et_parameter) {
1476 encounter_threshold = et_parameter;
1477 }
1478
1479 static void load_model(const char *name) {
1480 load_model_name = name;
1481 }
1482
1483 static void save_model(const char *name) {
1484 save_model_name = name;
1485 }
1486
1487 static void fc(ale_pos fc) {
1488 front_clip = fc;
1489 }
1490
1491 static void di_upper(ale_pos _dgi) {
1492 primary_decimation_upper = (int) round(_dgi);
1493 }
1494
1495 static void do_try(ale_pos _dgo) {
1496 output_decimation_preferred = (int) round(_dgo);
1497 }
1498
1499 static void di_lower(ale_pos _idiv) {
1500 input_decimation_lower = (int) round(_idiv);
1501 }
1502
1503 static void oc() {
1504 output_clip = 1;
1505 }
1506
1507 static void no_oc() {
1508 output_clip = 0;
1509 }
1510
1511 static void rc(ale_pos rc) {
1512 rear_clip = rc;
1513 }
1514
1515 /*
1516 * Initialize 3D scene from 2D scene, using 2D and 3D alignment
1517 * information.
1518 */
1519 static void init_from_d2() {
1520
1521 /*
1522 * Rear clip value of 0 is converted to infinity.
1523 */
1524
1525 if (rear_clip == 0) {
1526 ale_pos one = +1;
1527 ale_pos zero = +0;
1528
1529 rear_clip = one / zero;
1530 assert(isinf(rear_clip) == +1);
1531 }
1532
1533 /*
1534 * Scale and translate clipping plane depths.
1535 */
1536
1537 ale_pos cp_scalar = d3::align::projective(0).wc(point(0, 0, 0))[2];
1538
1539 front_clip = front_clip * cp_scalar - cp_scalar;
1540 rear_clip = rear_clip * cp_scalar - cp_scalar;
1541
1542 fprintf(stderr, "front_clip=%f rear_clip=%f\n", front_clip, rear_clip);
1543
1544 /*
1545 * Allocate image structures.
1546 */
1547
1548 al = new lod_images;
1549
1550 if (tc_multiplier != 0) {
1551 al->open_all();
1552 }
1553 }
1554
1555 /*
1556 * Perform spatial_info updating on a given subspace, for given
1557 * parameters.
1558 */
1559 static void subspace_info_update(space::iterate si, int f, ref_weights *weights) {
1560 while(!si.done()) {
1561
1562 space::traverse st = si.get();
1563
1564 /*
1565 * Skip spaces with no color information.
1566 *
1567 * XXX: This could be more efficient, perhaps.
1568 */
1569
1570 if (spatial_info_map.count(st.get_node()) == 0) {
1571 si.next();
1572 continue;
1573 }
1574
1575 /*
1576 * Get in-bounds centroid, if one exists.
1577 */
1578
1579 point p = al->get(f)->get_t(0).centroid(st);
1580
1581 if (!p.defined()) {
1582 si.next();
1583 continue;
1584 }
1585
1586 /*
1587 * Get information on the subspace.
1588 */
1589
1590 spatial_info *sn = &spatial_info_map[st.get_node()];
1591 d2::pixel color = sn->get_color();
1592 ale_real occupancy = sn->get_occupancy();
1593
1594 /*
1595 * Store current weight so we can later check for
1596 * modification by higher-resolution subspaces.
1597 */
1598
1599 ref_weights::subtree *tree = weights->get_subtree(st);
1600
1601 /*
1602 * Check for higher resolution subspaces, and
1603 * update the space iterator.
1604 */
1605
1606 if (st.get_node()->positive
1607 || st.get_node()->negative) {
1608
1609 /*
1610 * Cleave space for the higher-resolution pass,
1611 * skipping the current space, since we will
1612 * process that later.
1613 */
1614
1615 space::iterate cleaved_space = si.cleave();
1616
1617 cleaved_space.next();
1618
1619 subspace_info_update(cleaved_space, f, weights);
1620
1621 } else {
1622 si.next();
1623 }
1624
1625 /*
1626 * Add new data on the subspace and update weights.
1627 */
1628
1629 ale_pos tc = al->get(f)->get_t(0).trilinear_coordinate(st);
1630 d2::pixel pcolor = al->get(f)->get_tl(p.xy(), tc);
1631 d2::pixel colordiff = (color - pcolor) * (ale_real) 256;
1632
1633 if (falloff_exponent != 0) {
1634 d2::pixel max_diff = al->get(f)->get_max_diff(p.xy(), tc) * (ale_real) 256;
1635
1636 for (int k = 0; k < 3; k++)
1637 if (max_diff[k] > 1)
1638 colordiff[k] /= pow(max_diff[k], falloff_exponent);
1639 }
1640
1641 /*
1642 * Determine the probability of
1643 * encounter, divided by the occupancy.
1644 */
1645
1646 d2::pixel encounter = d2::pixel(1, 1, 1) * (1 - weights->get_weight(st));
1647
1648 /*
1649 * Update weights
1650 */
1651
1652 weights->add_weight(st, (encounter * occupancy)[0], tree);
1653
1654 /*
1655 * Delete the subtree, if necessary.
1656 */
1657
1658 delete tree;
1659
1660 /*
1661 * Check for cases in which the subspace should not be
1662 * updated.
1663 */
1664
1665 if (!resolution_ok(al->get(f)->get_t(0), tc))
1666 return;
1667
1668 /*
1669 * Update subspace.
1670 */
1671
1672 sn->accumulate_color_1(f, pcolor, encounter);
1673 d2::pixel channel_occ = pexp(-colordiff * colordiff);
1674
1675 ale_accum occ = channel_occ[0];
1676
1677 for (int k = 1; k < 3; k++)
1678 if (channel_occ[k] < occ)
1679 occ = channel_occ[k];
1680
1681 sn->accumulate_occupancy_1(f, occ, encounter[0]);
1682
1683 }
1684 }
1685
1686 /*
1687 * Run a single iteration of the spatial_info update cycle.
1688 */
1689 static void spatial_info_update() {
1690 /*
1691 * Iterate through each frame.
1692 */
1693 for (unsigned int f = 0; f < d2::image_rw::count(); f++) {
1694
1695 /*
1696 * Open the frame and transformation.
1697 */
1698
1699 if (tc_multiplier == 0)
1700 al->open(f);
1701
1702 /*
1703 * Allocate weights data structure for storing encounter
1704 * probabilities.
1705 */
1706
1707 ref_weights *weights = new ref_weights(f);
1708
1709 /*
1710 * Call subspace_info_update for the root space.
1711 */
1712
1713 subspace_info_update(space::iterate(al->get(f)->origin()), f, weights);
1714
1715 /*
1716 * Free weights.
1717 */
1718
1719 delete weights;
1720
1721 /*
1722 * Close the frame and transformation.
1723 */
1724
1725 if (tc_multiplier == 0)
1726 al->close(f);
1727 }
1728
1729 /*
1730 * Update all spatial_info structures.
1731 */
1732 for (spatial_info_map_t::iterator i = spatial_info_map.begin(); i != spatial_info_map.end(); i++) {
1733 i->second.update_color();
1734 i->second.update_occupancy();
1735
1736 // d2::pixel color = i->second.get_color();
1737
1738 // fprintf(stderr, "space p=%p updated to c=[%f %f %f] o=%f\n",
1739 // i->first, color[0], color[1], color[2],
1740 // i->second.get_occupancy());
1741 }
1742 }
1743
1744 /*
1745 * Support function for view() and depth().
1746 */
1747
1748 static const void view_recurse(int type, d2::image *im, d2::image *weights, space::iterate si, pt _pt) {
1749 while (!si.done()) {
1750 space::traverse st = si.get();
1751
1752 /*
1753 * XXX: This could be more efficient, perhaps.
1754 */
1755
1756 if (spatial_info_map.count(st.get_node()) == 0) {
1757 si.next();
1758 continue;
1759 }
1760
1761 spatial_info sn = spatial_info_map[st.get_node()];
1762
1763 /*
1764 * Get information on the subspace.
1765 */
1766
1767 d2::pixel color = sn.get_color();
1768 // d2::pixel color = d2::pixel(1, 1, 1) * (double) (((unsigned int) (st.get_node()) / sizeof(space)) % 65535);
1769 ale_real occupancy = sn.get_occupancy();
1770
1771 /*
1772 * Determine the view-local bounding box for the
1773 * subspace.
1774 */
1775
1776 point bb[2];
1777
1778 _pt.get_view_local_bb_scaled(st, bb);
1779
1780 point min = bb[0];
1781 point max = bb[1];
1782
1783 /*
1784 * Data structure to check modification of weights by
1785 * higher-resolution subspaces.
1786 */
1787
1788 std::queue<d2::pixel> weight_queue;
1789
1790 /*
1791 * Check for higher resolution subspaces, and
1792 * update the space iterator.
1793 */
1794
1795 if (st.get_node()->positive
1796 || st.get_node()->negative) {
1797
1798 /*
1799 * Store information about current weights,
1800 * so we will know which areas have been
1801 * covered by higher-resolution subspaces.
1802 */
1803
1804 for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
1805 for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++)
1806 weight_queue.push(weights->get_pixel(i, j));
1807
1808 /*
1809 * Cleave space for the higher-resolution pass,
1810 * skipping the current space, since we will
1811 * process that afterward.
1812 */
1813
1814 space::iterate cleaved_space = si.cleave();
1815
1816 cleaved_space.next();
1817
1818 view_recurse(type, im, weights, cleaved_space, _pt);
1819
1820 } else {
1821 si.next();
1822 }
1823
1824
1825 /*
1826 * Iterate over pixels in the bounding box, finding
1827 * pixels that intersect the subspace. XXX: assume
1828 * for now that all pixels in the bounding box
1829 * intersect the subspace.
1830 */
1831
1832 for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
1833 for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++) {
1834
1835 /*
1836 * Check for higher-resolution updates.
1837 */
1838
1839 if (weight_queue.size()) {
1840 if (weight_queue.front() != weights->get_pixel(i, j)) {
1841 weight_queue.pop();
1842 continue;
1843 }
1844 weight_queue.pop();
1845 }
1846
1847 /*
1848 * Determine the probability of encounter.
1849 */
1850
1851 d2::pixel encounter = (d2::pixel(1, 1, 1) - weights->get_pixel(i, j)) * occupancy;
1852
1853 /*
1854 * Update images.
1855 */
1856
1857 if (type == 0) {
1858
1859 /*
1860 * Color view
1861 */
1862
1863 weights->pix(i, j) += encounter;
1864 im->pix(i, j) += encounter * color;
1865
1866 } else if (type == 1) {
1867
1868 /*
1869 * Weighted (transparent) depth display
1870 */
1871
1872 ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
1873 weights->pix(i, j) += encounter;
1874 im->pix(i, j) += encounter * depth_value;
1875
1876 } else if (type == 2) {
1877
1878 /*
1879 * Ambiguity (ambivalence) measure.
1880 */
1881
1882 weights->pix(i, j) = d2::pixel(1, 1, 1);
1883 im->pix(i, j) += 0.1 * d2::pixel(1, 1, 1);
1884
1885 } else if (type == 3) {
1886
1887 /*
1888 * Closeness measure.
1889 */
1890
1891 ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
1892 if (weights->pix(i, j)[0] == 0) {
1893 weights->pix(i, j) = d2::pixel(1, 1, 1);
1894 im->pix(i, j) = d2::pixel(1, 1, 1) * depth_value;
1895 } else if (im->pix(i, j)[2] < depth_value) {
1896 im->pix(i, j) = d2::pixel(1, 1, 1) * depth_value;
1897 } else {
1898 continue;
1899 }
1900
1901 } else if (type == 4) {
1902
1903 /*
1904 * Weighted (transparent) contribution display
1905 */
1906
1907 ale_pos contribution_value = sn.get_pocc_density() + sn.get_socc_density();
1908 weights->pix(i, j) += encounter;
1909 im->pix(i, j) += encounter * contribution_value;
1910
1911 } else if (type == 5) {
1912
1913 /*
1914 * Weighted (transparent) occupancy display
1915 */
1916
1917 ale_pos contribution_value = occupancy;
1918 weights->pix(i, j) += encounter;
1919 im->pix(i, j) += encounter * contribution_value;
1920
1921 } else if (type == 6) {
1922
1923 /*
1924 * (Depth, xres, yres) triple
1925 */
1926
1927 ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
1928 weights->pix(i, j)[0] += encounter[0];
1929 if (weights->pix(i, j)[1] < encounter[0]) {
1930 weights->pix(i, j)[1] = encounter[0];
1931 im->pix(i, j)[0] = weights->pix(i, j)[1] * depth_value;
1932 im->pix(i, j)[1] = max[0] - min[0];
1933 im->pix(i, j)[2] = max[1] - min[1];
1934 }
1935
1936 } else if (type == 7) {
1937
1938 /*
1939 * (xoff, yoff, 0) triple
1940 */
1941
1942 weights->pix(i, j)[0] += encounter[0];
1943 if (weights->pix(i, j)[1] < encounter[0]) {
1944 weights->pix(i, j)[1] = encounter[0];
1945 im->pix(i, j)[0] = i - min[0];
1946 im->pix(i, j)[1] = j - min[1];
1947 im->pix(i, j)[2] = 0;
1948 }
1949
1950 } else
1951 assert(0);
1952 }
1953 }
1954 }
1955
1956 /*
1957 * Generate an depth image from a specified view.
1958 */
1959 static const d2::image *depth(pt _pt, int n = -1) {
1960 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
1961
1962 if (n >= 0) {
1963 assert((int) floor(d2::align::of(n).scaled_height())
1964 == (int) floor(_pt.scaled_height()));
1965 assert((int) floor(d2::align::of(n).scaled_width())
1966 == (int) floor(_pt.scaled_width()));
1967 }
1968
1969 d2::image *im1 = new d2::image_ale_real((int) floor(_pt.scaled_height()),
1970 (int) floor(_pt.scaled_width()), 3);
1971
1972 d2::image *im2 = new d2::image_ale_real((int) floor(_pt.scaled_height()),
1973 (int) floor(_pt.scaled_width()), 3);
1974
1975 d2::image *im3 = new d2::image_ale_real((int) floor(_pt.scaled_height()),
1976 (int) floor(_pt.scaled_width()), 3);
1977
1978 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
1979
1980 /*
1981 * Use adaptive subspace data.
1982 */
1983
1984 d2::image *weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
1985 (int) floor(_pt.scaled_width()), 3);
1986
1987 /*
1988 * Iterate through subspaces.
1989 */
1990
1991 space::iterate si(_pt.origin());
1992
1993 view_recurse(6, im1, weights, si, _pt);
1994
1995 delete weights;
1996 weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
1997 (int) floor(_pt.scaled_width()), 3);
1998
1999 view_recurse(7, im2, weights, si, _pt);
2000
2001 /*
2002 * Normalize depths by weights
2003 */
2004
2005 if (normalize_weights)
2006 for (unsigned int i = 0; i < im1->height(); i++)
2007 for (unsigned int j = 0; j < im1->width(); j++)
2008 im1->pix(i, j)[0] /= weights->pix(i, j)[1];
2009
2010
2011 for (unsigned int i = 0; i < im1->height(); i++)
2012 for (unsigned int j = 0; j < im1->width(); j++) {
2013
2014 /*
2015 * Handle interpolation.
2016 */
2017
2018 d2::point x;
2019 d2::point blx;
2020 d2::point res(im1->pix(i, j)[1], im1->pix(i, j)[2]);
2021
2022 for (int d = 0; d < 2; d++) {
2023
2024 if (im2->pix(i, j)[d] < res[d] / 2)
2025 x[d] = (ale_pos) (d?j:i) - res[d] / 2 - im2->pix(i, j)[d];
2026 else
2027 x[d] = (ale_pos) (d?j:i) + res[d] / 2 - im2->pix(i, j)[d];
2028
2029 blx[d] = 1 - ((d?j:i) - x[d]) / res[d];
2030 }
2031
2032 ale_real depth_val = 0;
2033 ale_real depth_weight = 0;
2034
2035 for (int ii = 0; ii < 2; ii++)
2036 for (int jj = 0; jj < 2; jj++) {
2037 d2::point p = x + d2::point(ii, jj) * res;
2038 if (im1->in_bounds(p)) {
2039
2040 ale_real d = im1->get_bl(p)[0];
2041
2042 if (isnan(d))
2043 continue;
2044
2045 ale_real w = ((ii ? (1 - blx[0]) : blx[0]) * (jj ? (1 - blx[1]) : blx[1]));
2046 depth_weight += w;
2047 depth_val += w * d;
2048 }
2049 }
2050
2051 ale_real depth = depth_val / depth_weight;
2052
2053 /*
2054 * Handle exclusions and encounter thresholds
2055 */
2056
2057 point w = _pt.pw_scaled(point(i, j, depth));
2058
2059 if (weights->pix(i, j)[0] < encounter_threshold || excluded(w)) {
2060 im3->pix(i, j) = d2::pixel::zero() / d2::pixel::zero();
2061 } else {
2062 im3->pix(i, j) = d2::pixel(1, 1, 1) * depth;
2063 }
2064 }
2065
2066 delete weights;
2067 delete im1;
2068 delete im2;
2069
2070 return im3;
2071 }
2072
2073 static const d2::image *depth(unsigned int n) {
2074
2075 assert (n < d2::image_rw::count());
2076
2077 pt _pt = align::projective(n);
2078
2079 return depth(_pt, n);
2080 }
2081
2082 /*
2083 * Generate an image from a specified view.
2084 */
2085 static const d2::image *view(pt _pt, int n = -1) {
2086 assert ((unsigned int) n < d2::image_rw::count() || n < 0);
2087
2088 if (n >= 0) {
2089 assert((int) floor(d2::align::of(n).scaled_height())
2090 == (int) floor(_pt.scaled_height()));
2091 assert((int) floor(d2::align::of(n).scaled_width())
2092 == (int) floor(_pt.scaled_width()));
2093 }
2094
2095 const d2::image *depths = NULL;
2096
2097 if (d3px_count > 0)
2098 depths = depth(_pt, n);
2099
2100 d2::image *im = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2101 (int) floor(_pt.scaled_width()), 3);
2102
2103 _pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);
2104
2105 /*
2106 * Use adaptive subspace data.
2107 */
2108
2109 d2::image *weights = new d2::image_ale_real((int) floor(_pt.scaled_height()),
2110 (int) floor(_pt.scaled_width()), 3);
2111
2112 /*
2113 * Iterate through subspaces.
2114 */
2115
2116 space::iterate si(_pt.origin());
2117
2118 view_recurse(0, im, weights, si, _pt);
2119
2120 for (unsigned int i = 0; i < im->height(); i++)
2121 for (unsigned int j = 0; j < im->width(); j++) {
2122 if (weights->pix(i, j).min_norm() < encounter_threshold
2123 || (d3px_count > 0 && isnan(depths->pix(i, j)[0]))) {
2124 im->pix(i, j) = d2::pixel::zero() / d2::pixel::zero();
2125 weights->pix(i, j) = d2::pixel::zero();
2126 } else if (normalize_weights)
2127 im->pix(i, j) /= weights->pix(i, j);
2128 }
2129
2130 if (d3px_count > 0)
2131 delete depths;
2132
2133 delete weights;
2134
2135 return im;
2136 }
2137
2138 static void tcem(double _tcem) {
2139 tc_multiplier = _tcem;
2140 }
2141
2142 static void oui(unsigned int _oui) {
2143 ou_iterations = _oui;
2144 }
2145
2146 static void pa(unsigned int _pa) {
2147 pairwise_ambiguity = _pa;
2148 }
2149
2150 static void pc(const char *_pc) {
2151 pairwise_comparisons = _pc;
2152 }
2153
2154 static void d3px(int _d3px_count, double *_d3px_parameters) {
2155 d3px_count = _d3px_count;
2156 d3px_parameters = _d3px_parameters;
2157 }
2158
2159 static void fx(double _fx) {
2160 falloff_exponent = _fx;
2161 }
2162
2163 static void nw() {
2164 normalize_weights = 1;
2165 }
2166
2167 static void no_nw() {
2168 normalize_weights = 0;
2169 }
2170
2171 static int excluded(point p) {
2172 for (int n = 0; n < d3px_count; n++) {
2173 double *region = d3px_parameters + (6 * n);
2174 if (p[0] >= region[0]
2175 && p[0] <= region[1]
2176 && p[1] >= region[2]
2177 && p[1] <= region[3]
2178 && p[2] >= region[4]
2179 && p[2] <= region[5])
2180 return 1;
2181 }
2182
2183 return 0;
2184 }
2185
2186 static const d2::image *view(unsigned int n) {
2187
2188 assert (n < d2::image_rw::count());
2189
2190 pt _pt = align::projective(n);
2191
2192 return view(_pt, n);
2193 }
2194
2195 /*
2196 * Add specified control points to the model
2197 */
2198 static void add_control_points() {
2199 }
2200
2201 typedef struct {point iw; point ip, is;} analytic;
2202 typedef std::multimap<ale_real,analytic> score_map;
2203 typedef std::pair<ale_real,analytic> score_map_element;
2204
2205 /*
2206 * Make pt list.
2207 */
2208 static std::vector<pt> make_pt_list(const char *d_out[], const char *v_out[],
2209 std::map<const char *, pt> *d3_depth_pt,
2210 std::map<const char *, pt> *d3_output_pt) {
2211
2212 std::vector<pt> result;
2213
2214 for (unsigned int n = 0; n < d2::image_rw::count(); n++) {
2215 if (d_out[n] || v_out[n]) {
2216 result.push_back(align::projective(n));
2217 }
2218 }
2219
2220 for (std::map<const char *, pt>::iterator i = d3_depth_pt->begin(); i != d3_depth_pt->end(); i++) {
2221 result.push_back(i->second);
2222 }
2223
2224 for (std::map<const char *, pt>::iterator i = d3_output_pt->begin(); i != d3_output_pt->end(); i++) {
2225 result.push_back(i->second);
2226 }
2227
2228 return result;
2229 }
2230
2231 /*
2232 * Get a trilinear coordinate for an anisotropic candidate cell.
2233 */
2234 static ale_pos get_trilinear_coordinate(point min, point max, pt _pt) {
2235
2236 d2::point local_min, local_max;
2237
2238 local_min = _pt.wp_unscaled(min).xy();
2239 local_max = _pt.wp_unscaled(min).xy();
2240
2241 point cell[2] = {min, max};
2242
2243 /*
2244 * Determine the view-local extrema in 2 dimensions.
2245 */
2246
2247 for (int r = 1; r < 8; r++) {
2248 point local = _pt.wp_unscaled(point(cell[r>>2][0], cell[(r>>1)%2][1], cell[r%2][2]));
2249
2250 for (int d = 0; d < 2; d++) {
2251 if (local[d] < local_min[d])
2252 local_min[d] = local[d];
2253 if (local[d] > local_max[d])
2254 local_max[d] = local[d];
2255 if (isnan(local[d]))
2256 return local[d];
2257 }
2258 }
2259
2260 ale_pos diameter = (local_max - local_min).norm();
2261
2262 return log(diameter / sqrt(2)) / log(2);
2263 }
2264
2265 /*
2266 * Check whether a cell is visible from a given viewpoint. This function
2267 * is guaranteed to return 1 when a cell is visible, but it is not guaranteed
2268 * to return 0 when a cell is invisible.
2269 */
2270 static int pt_might_be_visible(const pt &viewpoint, point min, point max) {
2271
2272 int doc = (rand() % 100000) ? 0 : 1;
2273
2274 if (doc)
2275 fprintf(stderr, "checking visibility:\n");
2276
2277 point cell[2] = {min, max};
2278
2279 /*
2280 * Cycle through all vertices of the cell to check certain
2281 * properties.
2282 */
2283 int pos[3] = {0, 0, 0};
2284 int neg[3] = {0, 0, 0};
2285 for (int i = 0; i < 2; i++)
2286 for (int j = 0; j < 2; j++)
2287 for (int k = 0; k < 2; k++) {
2288 point p = viewpoint.wp_unscaled(point(cell[i][0], cell[j][1], cell[k][2]));
2289
2290 if (p[2] < 0 && viewpoint.unscaled_in_bounds(p))
2291 return 1;
2292
2293 if (isnan(p[0])
2294 || isnan(p[1])
2295 || isnan(p[2]))
2296 return 1;
2297
2298 if (p[2] > 0)
2299 for (int d = 0; d < 2; d++)
2300 p[d] *= -1;
2301
2302 if (doc)
2303 fprintf(stderr, "\t[%f %f %f] --> [%f %f %f]\n",
2304 cell[i][0], cell[j][1], cell[k][2],
2305 p[0], p[1], p[2]);
2306
2307 for (int d = 0; d < 3; d++)
2308 if (p[d] >= 0)
2309 pos[d] = 1;
2310
2311 if (p[0] <= viewpoint.unscaled_height() - 1)
2312 neg[0] = 1;
2313
2314 if (p[1] <= viewpoint.unscaled_width() - 1)
2315 neg[1] = 1;
2316
2317 if (p[2] <= 0)
2318 neg[2] = 1;
2319 }
2320
2321 if (!neg[2])
2322 return 0;
2323
2324 if (!pos[0]
2325 || !neg[0]
2326 || !pos[1]
2327 || !neg[1])
2328 return 0;
2329
2330 return 1;
2331 }
2332
2333 /*
2334 * Check whether a cell is output-visible.
2335 */
2336 static int output_might_be_visible(const std::vector<pt> &pt_outputs, point min, point max) {
2337 for (unsigned int n = 0; n < pt_outputs.size(); n++)
2338 if (pt_might_be_visible(pt_outputs[n], min, max))
2339 return 1;
2340 return 0;
2341 }
2342
2343 /*
2344 * Check whether a cell is input-visible.
2345 */
2346 static int input_might_be_visible(unsigned int f, point min, point max) {
2347 return pt_might_be_visible(align::projective(f), min, max);
2348 }
2349
2350 /*
2351 * Return true if a cell fails an output resolution bound.
2352 */
2353 static int fails_output_resolution_bound(point min, point max, const std::vector<pt> &pt_outputs) {
2354 for (unsigned int n = 0; n < pt_outputs.size(); n++) {
2355
2356 point p = pt_outputs[n].centroid(min, max);
2357
2358 if (!p.defined())
2359 continue;
2360
2361 if (get_trilinear_coordinate(min, max, pt_outputs[n]) < output_decimation_preferred)
2362 return 1;
2363 }
2364
2365 return 0;
2366 }
2367
2368 /*
2369 * Check lower-bound resolution constraints
2370 */
2371 static int exceeds_resolution_lower_bounds(unsigned int f1, unsigned int f2,
2372 point min, point max, const std::vector<pt> &pt_outputs) {
2373
2374 pt _pt = al->get(f1)->get_t(0);
2375 point p = _pt.centroid(min, max);
2376
2377 if (get_trilinear_coordinate(min, max, _pt) < input_decimation_lower)
2378 return 1;
2379
2380 if (fails_output_resolution_bound(min, max, pt_outputs))
2381 return 0;
2382
2383 if (get_trilinear_coordinate(min, max, _pt) < primary_decimation_upper)
2384 return 1;
2385
2386 return 0;
2387 }
2388
2389 /*
2390 * Try the candidate nearest to the specified cell.
2391 */
2392 static void try_nearest_candidate(unsigned int f1, unsigned int f2, candidates *c, point min, point max) {
2393 point centroid = (max + min) / 2;
2394 pt _pt[2] = { al->get(f1)->get_t(0), al->get(f2)->get_t(0) };
2395 point p[2];
2396
2397 // fprintf(stderr, "[tnc n=%f %f %f x=%f %f %f]\n", min[0], min[1], min[2], max[0], max[1], max[2]);
2398
2399 /*
2400 * Reject clipping plane violations.
2401 */
2402
2403 if (centroid[2] > front_clip
2404 || centroid[2] < rear_clip)
2405 return;
2406
2407 /*
2408 * Calculate projections.
2409 */
2410
2411 for (int n = 0; n < 2; n++) {
2412
2413 p[n] = _pt[n].wp_unscaled(centroid);
2414
2415 if (!_pt[n].unscaled_in_bounds(p[n]))
2416 return;
2417
2418 // fprintf(stderr, ":");
2419
2420 if (p[n][2] >= 0)
2421 return;
2422 }
2423
2424
2425 int tc = (int) round(get_trilinear_coordinate(min, max, _pt[0]));
2426 int stc = (int) round(get_trilinear_coordinate(min, max, _pt[1]));
2427
2428 while (tc < input_decimation_lower || stc < input_decimation_lower) {
2429 tc++;
2430 stc++;
2431 }
2432
2433 if (tc > primary_decimation_upper)
2434 return;
2435
2436 /*
2437 * Calculate score from color match. Assume for now
2438 * that the transformation can be approximated locally
2439 * with a translation.
2440 */
2441
2442 ale_pos score = 0;
2443 ale_pos divisor = 0;
2444 ale_pos l1_multiplier = 0.125;
2445 lod_image *if1 = al->get(f1);
2446 lod_image *if2 = al->get(f2);
2447
2448 if (if1->in_bounds(p[0].xy())
2449 && if2->in_bounds(p[1].xy())) {
2450 divisor += 1 - l1_multiplier;
2451 score += (1 - l1_multiplier)
2452 * (if1->get_tl(p[0].xy(), tc) - if2->get_tl(p[1].xy(), stc)).normsq();
2453 }
2454
2455 for (int iii = -1; iii <= 1; iii++)
2456 for (int jjj = -1; jjj <= 1; jjj++) {
2457 d2::point t(iii, jjj);
2458
2459 if (!if1->in_bounds(p[0].xy() + t)
2460 || !if2->in_bounds(p[1].xy() + t))
2461 continue;
2462
2463 divisor += l1_multiplier;
2464 score += l1_multiplier
2465 * (if1->get_tl(p[0].xy() + t, tc) - if2->get_tl(p[1].xy() + t, tc)).normsq();
2466
2467 }
2468
2469 /*
2470 * Include third-camera contributions in the score.
2471 */
2472
2473 if (tc_multiplier != 0)
2474 for (unsigned int n = 0; n < d2::image_rw::count(); n++) {
2475 if (n == f1 || n == f2)
2476 continue;
2477
2478 lod_image *ifn = al->get(n);
2479 pt _ptn = ifn->get_t(0);
2480 point pn = _ptn.wp_unscaled(centroid);
2481
2482 if (!_ptn.unscaled_in_bounds(pn))
2483 continue;
2484
2485 if (pn[2] >= 0)
2486 continue;
2487
2488 ale_pos ttc = get_trilinear_coordinate(min, max, _ptn);
2489
2490 divisor += tc_multiplier;
2491 score += tc_multiplier
2492 * (if1->get_tl(p[0].xy(), tc) - ifn->get_tl(pn.xy(), ttc)).normsq();
2493 }
2494
2495 c->add_candidate(p[0], tc, score / divisor);
2496 }
2497
2498 /*
2499 * Check for cells that are completely clipped.
2500 */
2501 static int completely_clipped(point min, point max) {
2502 return (min[2] > front_clip
2503 || max[2] < rear_clip);
2504 }
2505
2506 /*
2507 * Update extremum variables for cell points mapped to a particular view.
2508 */
2509 static void update_extrema(point min, point max, pt _pt, int *extreme_dim, ale_pos *extreme_ratio) {
2510
2511 point local_min, local_max;
2512
2513 local_min = _pt.wp_unscaled(min);
2514 local_max = _pt.wp_unscaled(min);
2515
2516 point cell[2] = {min, max};
2517
2518 int near_vertex = 0;
2519
2520 /*
2521 * Determine the view-local extrema in all dimensions, and
2522 * determine the vertex of closest z coordinate.
2523 */
2524
2525 for (int r = 1; r < 8; r++) {
2526 point local = _pt.wp_unscaled(point(cell[r>>2][0], cell[(r>>1)%2][1], cell[r%2][2]));
2527
2528 for (int d = 0; d < 3; d++) {
2529 if (local[d] < local_min[d])
2530 local_min[d] = local[d];
2531 if (local[d] > local_max[d])
2532 local_max[d] = local[d];
2533 }
2534
2535 if (local[2] == local_max[2])
2536 near_vertex = r;
2537 }
2538
2539 ale_pos diameter = (local_max.xy() - local_min.xy()).norm();
2540
2541 /*
2542 * Update extrema as necessary for each dimension.
2543 */
2544
2545 for (int d = 0; d < 3; d++) {
2546
2547 int r = near_vertex;
2548
2549 int p1[3] = {r>>2, (r>>1)%2, r%2};
2550 int p2[3] = {r>>2, (r>>1)%2, r%2};
2551
2552 p2[d] = 1 - p2[d];
2553
2554 ale_pos local_distance = (_pt.wp_unscaled(point(cell[p1[0]][0], cell[p1[1]][1], cell[p1[2]][2])).xy()
2555 - _pt.wp_unscaled(point(cell[p2[0]][0], cell[p2[1]][1], cell[p2[2]][2])).xy()).norm();
2556
2557 if (local_distance / diameter > *extreme_ratio) {
2558 *extreme_ratio = local_distance / diameter;
2559 *extreme_dim = d;
2560 }
2561 }
2562 }
2563
2564 /*
2565 * Get the next split dimension.
2566 */
2567 static int get_next_split(int f1, int f2, point min, point max, const std::vector<pt> &pt_outputs) {
2568 for (int d = 0; d < 3; d++)
2569 if (isinf(min[d]) || isinf(max[d]))
2570 return space::traverse::get_next_split(min, max);
2571
2572 int extreme_dim = 0;
2573 ale_pos extreme_ratio = 0;
2574
2575 update_extrema(min, max, al->get(f1)->get_t(0), &extreme_dim, &extreme_ratio);
2576 update_extrema(min, max, al->get(f2)->get_t(0), &extreme_dim, &extreme_ratio);
2577
2578 for (unsigned int n = 0; n < pt_outputs.size(); n++) {
2579 update_extrema(min, max, pt_outputs[n], &extreme_dim, &extreme_ratio);
2580 }
2581
2582 return extreme_dim;
2583 }
2584
2585 /*
2586 * Find candidates for subspace creation.
2587 */
2588 static void find_candidates(unsigned int f1, unsigned int f2, candidates *c, point min, point max,
2589 const std::vector<pt> &pt_outputs, int depth = 0) {
2590
2591 int print = 0;
2592
2593 if (min[0] < 20.0001 && max[0] > 20.0001
2594 && min[1] < 20.0001 && max[1] > 20.0001
2595 && min[2] < 0.0001 && max[2] > 0.0001)
2596 print = 1;
2597
2598 if (print) {
2599 for (int i = depth; i > 0; i--) {
2600 fprintf(stderr, "+");
2601 }
2602 fprintf(stderr, "[fc n=%f %f %f x=%f %f %f]\n",
2603 min[0], min[1], min[2], max[0], max[1], max[2]);
2604 }
2605
2606 if (completely_clipped(min, max)) {
2607 if (print)
2608 fprintf(stderr, "c");
2609 return;
2610 }
2611
2612 if (!input_might_be_visible(f1, min, max)
2613 || !input_might_be_visible(f2, min, max)) {
2614 if (print)
2615 fprintf(stderr, "v");
2616 return;
2617 }
2618
2619 if (output_clip && !output_might_be_visible(pt_outputs, min, max)) {
2620 if (print)
2621 fprintf(stderr, "o");
2622 return;
2623 }
2624
2625 if (exceeds_resolution_lower_bounds(f1, f2, min, max, pt_outputs)) {
2626 if (!(rand() % 100000))
2627 fprintf(stderr, "([%f %f %f], [%f %f %f]) at %d\n",
2628 min[0], min[1], min[2],
2629 max[0], max[1], max[2],
2630 __LINE__);
2631
2632 if (print)
2633 fprintf(stderr, "t");
2634
2635 try_nearest_candidate(f1, f2, c, min, max);
2636 return;
2637 }
2638
2639 point new_cells[2][2];
2640
2641 if (!space::traverse::get_next_cells(get_next_split(f1, f2, min, max, pt_outputs), min, max, new_cells)) {
2642 if (print)
2643 fprintf(stderr, "n");
2644 return;
2645 }
2646
2647 if (print) {
2648 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",
2649 new_cells[0][0][0],
2650 new_cells[0][0][1],
2651 new_cells[0][0][2],
2652 new_cells[0][1][0],
2653 new_cells[0][1][1],
2654 new_cells[0][1][2],
2655 new_cells[1][0][0],
2656 new_cells[1][0][1],
2657 new_cells[1][0][2],
2658 new_cells[1][1][0],
2659 new_cells[1][1][1],
2660 new_cells[1][1][2]);
2661 }
2662
2663 find_candidates(f1, f2, c, new_cells[0][0], new_cells[0][1], pt_outputs, depth + 1);
2664 find_candidates(f1, f2, c, new_cells[1][0], new_cells[1][1], pt_outputs, depth + 1);
2665 }
2666
2667 /*
2668 * Initialize space and identify regions of interest for the adaptive
2669 * subspace model.
2670 */
2671 static void make_space(const char *d_out[], const char *v_out[],
2672 std::map<const char *, pt> *d3_depth_pt,
2673 std::map<const char *, pt> *d3_output_pt) {
2674
2675 fprintf(stderr, "Subdividing 3D space");
2676
2677 std::vector<pt> pt_outputs = make_pt_list(d_out, v_out, d3_depth_pt, d3_output_pt);
2678
2679 /*
2680 * Initialize root space.
2681 */
2682
2683 space::init_root();
2684
2685 /*
2686 * Subdivide space to resolve intensity matches between pairs
2687 * of frames.
2688 */
2689
2690 for (unsigned int f1 = 0; f1 < d2::image_rw::count(); f1++) {
2691
2692 if (d3_depth_pt->size() == 0
2693 && d3_output_pt->size() == 0
2694 && d_out[f1] == NULL
2695 && v_out[f1] == NULL)
2696 continue;
2697
2698 if (tc_multiplier == 0)
2699 al->open(f1);
2700
2701 for (unsigned int f2 = 0; f2 < d2::image_rw::count(); f2++) {
2702
2703 if (f1 == f2)
2704 continue;
2705
2706 if (tc_multiplier == 0)
2707 al->open(f2);
2708
2709 candidates *c = new candidates(f1);
2710
2711 find_candidates(f1, f2, c, point::neginf(), point::posinf(), pt_outputs);
2712
2713
2714
2715 c->generate_subspaces();
2716
2717 if (tc_multiplier == 0)
2718 al->close(f2);
2719 }
2720
2721 if (tc_multiplier == 0)
2722 al->close(f1);
2723 }
2724
2725 fprintf(stderr, "Final spatial map size: %u\n", spatial_info_map.size());
2726
2727 fprintf(stderr, ".\n");
2728 }
2729
2730
2731 /*
2732 * Update spatial information structures.
2733 *
2734 * XXX: the name of this function is horribly misleading. There isn't
2735 * even a 'search depth' any longer, since there is no longer any
2736 * bounded DFS occurring.
2737 */
2738 static void reduce_cost_to_search_depth(d2::exposure *exp_out, int inc_bit) {
2739
2740 /*
2741 * Subspace model
2742 */
2743
2744 for (unsigned int i = 0; i < ou_iterations; i++)
2745 spatial_info_update();
2746
2747 fprintf(stderr, "Final spatial map size: %u\n", spatial_info_map.size());
2748 }
2749
2750 #if 0
2751 /*
2752 * Describe a scene to a renderer
2753 */
2754 static void describe(render *r) {
2755 }
2756 #endif
2757 };
2758
2759 #endif
synthetic capture engine and renderer