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