// Copyright 2003, 2004, 2005, 2006 David Hilvert <dhilvert@auricle.dyndns.org>,
//                                                <dhilvert@ugcs.caltech.edu>

/*  This file is part of the Anti-Lamenessing Engine.

    The Anti-Lamenessing Engine is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    The Anti-Lamenessing Engine is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with the Anti-Lamenessing Engine; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

/*
 * d3/scene.h: Representation of a 3D scene.
 */

#ifndef __scene_h__
#define __scene_h__

#include "point.h"

/*
 * View angle multiplier.  
 *
 * Setting this to a value larger than one can be useful for debugging.
 */

#define VIEW_ANGLE_MULTIPLIER 1

class scene {

	/*
	 * Clipping planes
	 */
	static ale_pos front_clip;
	static ale_pos rear_clip;

	/*
	 * Decimation exponents for geometry
	 */
	static int primary_decimation_upper;
	static int input_decimation_lower;
	static int output_decimation_preferred;

	/*
	 * Output clipping
	 */
	static int output_clip;

	/*
	 * Model files
	 */
	static const char *load_model_name;
	static const char *save_model_name;

	/*
	 * Occupancy attenuation
	 */

	static double occ_att;

	/*
	 * Normalization of output by weight
	 */

	static int normalize_weights;

	/*
	 * Filtering data
	 */

	static int use_filter;
	static const char *d3chain_type;

	/*
	 * Falloff exponent
	 */

	static ale_real falloff_exponent;

	/*
	 * Third-camera error multiplier
	 */
	static double tc_multiplier;

	/*
	 * Occupancy update iterations
	 */
	static unsigned int ou_iterations;

	/*
	 * Pairwise ambiguity
	 */
	static unsigned int pairwise_ambiguity;

	/*
	 * Pairwise comparisons
	 */
	static const char *pairwise_comparisons;

	/*
	 * 3D Post-exclusion
	 */
	static int d3px_count;
	static double *d3px_parameters;

	/*
	 * Nearness threshold
	 */
	static const ale_real nearness;

	/*
	 * Encounter threshold for defined pixels.
	 */
	static double encounter_threshold;

	/*
	 * Median calculation radii.
	 */
	static double depth_median_radius;
	static double diff_median_radius;

	/*
	 * Flag for subspace traversal.
	 */
	static int subspace_traverse;

	/*
	 * Structure to hold input frame information at levels of 
	 * detail between full detail and full decimation.
	 */
	class lod_image {
		unsigned int f;
		unsigned int entries;
		std::vector<const d2::image *> im;
		std::vector<pt> transformation;

	public:
		/*
		 * Constructor
		 */
		lod_image(unsigned int _f) {

			pt _pt;
			
			f = _f;
			im.push_back(d2::image_rw::copy(f, "3D reference image"));
			assert(im.back());
			entries = 1;
			_pt = d3::align::projective(f);
			_pt.scale(1 / _pt.scale_2d());
			transformation.push_back(_pt);

			while (im.back()->height() > 4
			    && im.back()->width() > 4) {

				im.push_back(im.back()->scale_by_half("3D, reduced LOD"));
				assert(im.back());

				_pt.scale(1 / _pt.scale_2d() / pow((ale_pos) 2, entries));
				transformation.push_back(_pt);

				entries += 1;
			}
		}

		/*
		 * Get the number of scales
		 */
		unsigned int count() {
			return entries;
		}

		/*
		 * Get an image.
		 */
		const d2::image *get_image(unsigned int i) {
			assert(i < entries);
			return im[i];
		}

		int in_bounds(d2::point p) {
			return im[0]->in_bounds(p);
		}

		/*
		 * Get a 'trilinear' color.  We currently don't do interpolation
		 * between levels of detail; hence, it's discontinuous in tl_coord.
		 */
		d2::pixel get_tl(d2::point p, ale_pos tl_coord) {

			assert(in_bounds(p));

			tl_coord = round(tl_coord);

			if (tl_coord >= entries)
				tl_coord = entries;
			if (tl_coord < 0)
				tl_coord = 0;

			p = p / (ale_pos) pow(2, tl_coord);

			unsigned int itlc = (unsigned int) tl_coord;

			if (p[0] > im[itlc]->height() - 1)
				p[0] = im[itlc]->height() - 1;
			if (p[1] > im[itlc]->width() - 1)
				p[1] = im[itlc]->width() - 1;

			assert(p[0] >= 0);
			assert(p[1] >= 0);

			return im[itlc]->get_bl(p);
		}

		d2::pixel get_max_diff(d2::point p, ale_pos tl_coord) {
			assert(in_bounds(p));

			tl_coord = round(tl_coord);

			if (tl_coord >= entries)
				tl_coord = entries;
			if (tl_coord < 0)
				tl_coord = 0;

			p = p / (ale_pos) pow(2, tl_coord);

			unsigned int itlc = (unsigned int) tl_coord;

			if (p[0] > im[itlc]->height() - 1)
				p[0] = im[itlc]->height() - 1;
			if (p[1] > im[itlc]->width() - 1)
				p[1] = im[itlc]->width() - 1;

			assert(p[0] >= 0);
			assert(p[1] >= 0);

			return im[itlc]->get_max_diff(p);
		}

		/*
		 * Get the transformation
		 */
		pt get_t(unsigned int i) {
			assert(i >= 0);
			assert(i < entries);
			return transformation[i];
		}

		/*
		 * Get the camera origin in world coordinates
		 */
		point origin() {
			return transformation[0].origin();
		}

		/*
		 * Destructor
		 */
		~lod_image() {
			for (unsigned int i = 0; i < entries; i++) {
				delete im[i];
			}
		}
	};

	/*
	 * Structure to hold weight information for reference images.
	 */
	class ref_weights {
		unsigned int f;
		std::vector<d2::image *> weights;
		pt transformation;
		int tc_low;
		int tc_high;
		int initialized;

		void set_image(d2::image *im, ale_real value) {
			assert(im);
			for (unsigned int i = 0; i < im->height(); i++)
			for (unsigned int j = 0; j < im->width(); j++)
				im->set_pixel(i, j, d2::pixel(value, value, value));
		}

		d2::image *make_image(ale_pos sf, ale_real init_value = 0) {
			d2::image *result = d2::new_image_ale_real(
					(unsigned int) ceil(transformation.unscaled_height() * sf),
					(unsigned int) ceil(transformation.unscaled_width() * sf), 3);
			assert(result);

			if (init_value != 0)
				set_image(result, init_value);

			return result;
		}

	public:

		/*
		 * Explicit weight subtree
		 */
		struct subtree {
			ale_real node_value;
			subtree *children[2][2];

			subtree(ale_real nv, subtree *a, subtree *b, subtree *c, subtree *d) {
				node_value = nv;
				children[0][0] = a;
				children[0][1] = b;
				children[1][0] = c;
				children[1][1] = d;
			}

			~subtree() {
				for (int i = 0; i < 2; i++)
				for (int j = 0; j < 2; j++)
					delete children[i][j];
			}
		};

		/*
		 * Constructor
		 */
		ref_weights(unsigned int _f) {
			f = _f;
			transformation = d3::align::projective(f);
			initialized = 0;
		}

		/*
		 * Check spatial bounds.
		 */
		int in_spatial_bounds(point p) {

			if (!p.defined())
				return 0;

			if (p[0] < 0)
				return 0;
			if (p[1] < 0)
				return 0;
			if (p[0] > transformation.unscaled_height() - 1)
				return 0;
			if (p[1] > transformation.unscaled_width() - 1)
				return 0;
			if (p[2] >= 0)
				return 0;

			return 1;
		}

		int in_spatial_bounds(const space::traverse &t) {
			point p = transformation.centroid(t);
			return in_spatial_bounds(p);
		}

		/*
		 * Increase resolution to the given level.
		 */
		void increase_resolution(int tc, unsigned int i, unsigned int j) {
			d2::image *im = weights[tc - tc_low];
			assert(im);
			assert(i <= im->height() - 1);
			assert(j <= im->width() - 1);

			/*
			 * Check for the cases known to have no lower level of detail.
			 */

			if (im->get_chan(i, j, 0) == -1)
				return;

			if (tc == tc_high)
				return;

			increase_resolution(tc + 1, i / 2, j / 2);

			/*
			 * Load the lower-level image structure.
			 */

			d2::image *iim = weights[tc + 1 - tc_low];

			assert(iim);
			assert(i / 2 <= iim->height() - 1);
			assert(j / 2 <= iim->width() - 1);

			/*
			 * Check for the case where no lower level of detail is
			 * available.
			 */

			if (iim->get_chan(i / 2, j / 2, 0) == -1)
				return;

			/*
			 * Spread out the lower level of detail among (uninitialized)
			 * peer values.
			 */

			for (unsigned int ii = (i / 2) * 2; ii < (i / 2) * 2 + 1; ii++)
			for (unsigned int jj = (j / 2) * 2; jj < (j / 2) * 2 + 1; jj++) {
				assert(ii <= im->height() - 1);
				assert(jj <= im->width() - 1);
				assert(im->get_chan(ii, jj, 0) == 0);

				im->set_pixel(ii, jj, iim->get_pixel(i / 2, j / 2));
			}

			/*
			 * Set the lower level of detail to point here.
			 */

			iim->set_chan(i / 2, j / 2, 0, -1);
		}

		/*
		 * Add weights to positive higher-resolution pixels where
		 * found when their current values match the given subtree
		 * values; set negative pixels to zero and return 0 if no
		 * positive higher-resolution pixels are found.  
		 */
		int add_partial(int tc, unsigned int i, unsigned int j, ale_real weight, subtree *st) {
			d2::image *im = weights[tc - tc_low];
			assert(im);

			if (i == im->height() - 1
			 || j == im->width() - 1) {
				return 1;
			}

			assert(i <= im->height() - 1);
			assert(j <= im->width() - 1);

			/*
			 * Check for positive values.
			 */

			if (im->get_chan(i, j, 0) > 0) {
				if (st && st->node_value == im->get_pixel(i, j)[0])
					im->set_chan(i, j, 0, (ale_real) im->get_chan(i, j, 0) 
					                    + weight * (1 - im->get_pixel(i, j)[0]));
				return 1;
			}

			/*
			 * Handle the case where there are no higher levels of detail.
			 */

			if (tc == tc_low) {
				if (im->get_chan(i, j, 0) != 0) {
					fprintf(stderr, "failing assertion: im[%d]->pix(%d, %d)[0] == %g\n", tc, i, j, 
							(double) im->get_chan(i, j, 0));
				}
				assert(im->get_chan(i, j, 0) == 0);
				return 0;
			}

			/*
			 * Handle the case where higher levels of detail are available.
			 */

			int success[2][2];

			for (int ii = 0; ii < 2; ii++)
			for (int jj = 0; jj < 2; jj++)
				success[ii][jj] = add_partial(tc - 1, i * 2 + ii, j * 2 + jj, weight, 
						st ? st->children[ii][jj] : NULL);

			if (!success[0][0]
			 && !success[0][1]
			 && !success[1][0]
			 && !success[1][1]) {
				im->set_chan(i, j, 0, 0);
				return 0;
			}

			d2::image *iim = weights[tc - 1 - tc_low];
			assert(iim);

			for (int ii = 0; ii < 2; ii++)
			for (int jj = 0; jj < 2; jj++)
				if (success[ii][jj] == 0) {
					assert(i * 2 + ii < iim->height());
					assert(j * 2 + jj < iim->width());

					iim->set_chan(i * 2 + ii, j * 2 + jj, 0, weight);
				}

			im->set_chan(i, j, 0, -1);

			return 1;
		}

		/*
		 * Add weight.
		 */
		void add_weight(int tc, unsigned int i, unsigned int j, ale_real weight, subtree *st) {

			assert (weight >= 0);

			d2::image *im = weights[tc - tc_low];
			assert(im);

//			fprintf(stderr, "[aw tc=%d i=%d j=%d imax=%d jmax=%d]\n",
//					tc, i, j, im->height(), im->width());
			
			assert(i <= im->height() - 1);
			assert(j <= im->width() - 1);

			/*
			 * Increase resolution, if necessary
			 */

			increase_resolution(tc, i, j);

			/*
			 * Attempt to add the weight at levels of detail
			 * where weight is defined.
			 */

			if (add_partial(tc, i, j, weight, st))
				return;

			/*
			 * If no weights are defined at any level of detail,
			 * then set the weight here.
			 */

			im->set_chan(i, j, 0, weight);
		}

		void add_weight(int tc, d2::point p, ale_real weight, subtree *st) {

			assert (weight >= 0);

			p *= pow(2, -tc);

			unsigned int i = (unsigned int) floor(p[0]);
			unsigned int j = (unsigned int) floor(p[1]);

			add_weight(tc, i, j, weight, st);
		}

		void add_weight(const space::traverse &t, ale_real weight, subtree *st) {

			assert (weight >= 0);

			if (weight == 0)
				return;
			
			ale_pos tc = transformation.trilinear_coordinate(t);
			point p = transformation.centroid(t);
			assert(in_spatial_bounds(p));

			tc = round(tc);

			/*
			 * Establish a reasonable (?) upper bound on resolution.
			 */

			if (tc < input_decimation_lower) {
				weight /= pow(4, (input_decimation_lower - tc));
				tc = input_decimation_lower;
			}

			/*
			 * Initialize, if necessary.
			 */

			if (!initialized) {
				tc_low = tc_high = (int) tc;

				ale_pos sf = pow(2, -tc);

				weights.push_back(make_image(sf));

				initialized = 1;
			}

			/*
			 * Check resolution bounds
			 */

			assert (tc_low <= tc_high);

			/*
			 * Generate additional levels of detail, if necessary.
			 */

			while (tc < tc_low) {
				tc_low--;

				ale_pos sf = pow(2, -tc_low);

				weights.insert(weights.begin(), make_image(sf));
			}

			while (tc > tc_high) {
				tc_high++;

				ale_pos sf = pow(2, -tc_high);

				weights.push_back(make_image(sf, -1));
			}

			add_weight((int) tc, p.xy(), weight, st);
		}

		/*
		 * Get weight
		 */
		ale_real get_weight(int tc, unsigned int i, unsigned int j) {

//			fprintf(stderr, "[gw tc=%d i=%u j=%u tclow=%d tchigh=%d]\n", 
//					tc, i, j, tc_low, tc_high);

			if (tc < tc_low || !initialized)
				return 0;

			if (tc > tc_high) {
				return (get_weight(tc - 1, i * 2 + 0, j * 2 + 0)
				      + get_weight(tc - 1, i * 2 + 1, j * 2 + 0)
				      + get_weight(tc - 1, i * 2 + 1, j * 2 + 1)
				      + get_weight(tc - 1, i * 2 + 0, j * 2 + 1)) / 4;
			}

			assert(weights.size() > (unsigned int) (tc - tc_low));

			d2::image *im = weights[tc - tc_low];
			assert(im);

			if (i == im->height())
				return 1;
			if (j == im->width())
				return 1;

			assert(i < im->height());
			assert(j < im->width());

			if (im->get_chan(i, j, 0) == -1) {
				return (get_weight(tc - 1, i * 2 + 0, j * 2 + 0)
				      + get_weight(tc - 1, i * 2 + 1, j * 2 + 0)
				      + get_weight(tc - 1, i * 2 + 1, j * 2 + 1)
				      + get_weight(tc - 1, i * 2 + 0, j * 2 + 1)) / 4;
			}

			if (im->get_chan(i, j, 0) == 0) {
				if (tc == tc_high)
					return 0;
				if (weights[tc - tc_low + 1]->get_chan(i / 2, j / 2, 0) == -1)
					return 0;
				return get_weight(tc + 1, i / 2, j / 2);
			}

			return im->get_chan(i, j, 0);
		}

		ale_real get_weight(int tc, d2::point p) {

			p *= pow(2, -tc);

			unsigned int i = (unsigned int) floor(p[0]);
			unsigned int j = (unsigned int) floor(p[1]);

			return get_weight(tc, i, j);
		}

		ale_real get_weight(const space::traverse &t) {
			ale_pos tc = transformation.trilinear_coordinate(t);
			point p = transformation.centroid(t);
			assert(in_spatial_bounds(p));

			if (!initialized)
				return 0;

			tc = round(tc);

			if (tc < tc_low) {
				tc = tc_low;
			}

			return get_weight((int) tc, p.xy());
		}

		/*
		 * Check whether a subtree is simple.
		 */
		int is_simple(subtree *s) {
			assert (s);

			if (s->node_value == -1
			 && s->children[0][0] == NULL
			 && s->children[0][1] == NULL
			 && s->children[1][0] == NULL
			 && s->children[1][1] == NULL)
				return 1;

			return 0;
		}

		/*
		 * Get a weight subtree.
		 */
		subtree *get_subtree(int tc, unsigned int i, unsigned int j) {

			/*
			 * tc > tc_high is handled recursively.
			 */

			if (tc > tc_high) {
				subtree *result = new subtree(-1, 
						get_subtree(tc - 1, i * 2 + 0, j * 2 + 0),
						get_subtree(tc - 1, i * 2 + 0, j * 2 + 1),
						get_subtree(tc - 1, i * 2 + 1, j * 2 + 0),
						get_subtree(tc - 1, i * 2 + 1, j * 2 + 1));

				if (is_simple(result)) {
					delete result;
					return NULL;
				}

				return result;
			}

			assert(tc >= tc_low);
			assert(weights[tc - tc_low]);

			d2::image *im = weights[tc - tc_low];

			/*
			 * Rectangular images will, in general, have
			 * out-of-bounds tree sections.  Handle this case.
			 */

			if (i >= im->height())
				return NULL;
			if (j >= im->width())
				return NULL;

			/*
			 * -1 weights are handled recursively
			 */

			if (im->get_chan(i, j, 0) == -1) {
				subtree *result = new subtree(-1, 
						get_subtree(tc - 1, i * 2 + 0, j * 2 + 0),
						get_subtree(tc - 1, i * 2 + 0, j * 2 + 1),
						get_subtree(tc - 1, i * 2 + 1, j * 2 + 0),
						get_subtree(tc - 1, i * 2 + 1, j * 2 + 1));

				if (is_simple(result)) {
					im->set_chan(i, j, 0, 0);
					delete result;
					return NULL;
				}

				return result;
			}

			/*
			 * Zero weights have NULL subtrees.
			 */

			if (im->get_chan(i, j, 0) == 0)
				return NULL;

			/*
			 * Handle the remaining case.
			 */

			return new subtree(im->get_chan(i, j, 0), NULL, NULL, NULL, NULL);
		}

		subtree *get_subtree(int tc, d2::point p) {
			p *= pow(2, -tc);

			unsigned int i = (unsigned int) floor(p[0]);
			unsigned int j = (unsigned int) floor(p[1]);

			return get_subtree(tc, i, j);
		}

		subtree *get_subtree(const space::traverse &t) {
			ale_pos tc = transformation.trilinear_coordinate(t);
			point p = transformation.centroid(t);
			assert(in_spatial_bounds(p));

			if (!initialized)
				return NULL;

			if (tc < input_decimation_lower)
				tc = input_decimation_lower;

			tc = round(tc);

			if (tc < tc_low)
				return NULL;

			return get_subtree((int) tc, p.xy());
		}

		/*
		 * Destructor
		 */
		~ref_weights() {
			for (unsigned int i = 0; i < weights.size(); i++) {
				delete weights[i];
			}
		}
	};

	/*
	 * Resolution check.
	 */
	static int resolution_ok(pt transformation, ale_pos tc) {

		if (pow(2, tc) > transformation.unscaled_height()
		 || pow(2, tc) > transformation.unscaled_width())
			return 0;

		if (tc < input_decimation_lower - 1.5)
			return 0;

		return 1;
	}

	/*
	 * Structure to hold input frame information at all levels of detail.
	 */
	class lod_images {

		/*
		 * All images.
		 */

		std::vector<lod_image *> images;

	public:

		lod_images() {
			images.resize(d2::image_rw::count(), NULL);
		}

		unsigned int count() {
			return d2::image_rw::count();
		}

		void open(unsigned int f) {
			assert (images[f] == NULL);

			if (images[f] == NULL)
				images[f] = new lod_image(f);
		}

		void open_all() {
			for (unsigned int f = 0; f < d2::image_rw::count(); f++)
				open(f);
		}

		lod_image *get(unsigned int f) {
			assert (images[f] != NULL);
			return images[f];
		}

		void close(unsigned int f) {
			assert (images[f] != NULL);
			delete images[f];
			images[f] = NULL;
		}

		void close_all() {
			for (unsigned int f = 0; f < d2::image_rw::count(); f++)
				close(f);
		}

		~lod_images() {
			close_all();
		}
	};

	/*
	 * All levels-of-detail
	 */

	static struct lod_images *al;

	/*
	 * Data structure for storing best encountered subspace candidates.
	 */
	class candidates {
		std::vector<std::vector<std::pair<ale_pos, ale_real> > > levels;
		int image_index;
		unsigned int height;
		unsigned int width;

		/*
		 * Point p is in world coordinates.
		 */
		void generate_subspace(point iw, ale_pos diagonal) {

//			fprintf(stderr, "[gs iw=%f %f %f d=%f]\n", 
//					iw[0], iw[1], iw[2], diagonal);

			space::traverse st = space::traverse::root();

			if (!st.includes(iw)) {
				assert(0);
				return;
			}

			int highres = 0;
			int lowres = 0;

			/*
			 * Loop until resolutions of interest have been generated.
			 */
			
			for(;;) {

				ale_pos current_diagonal = (st.get_max() - st.get_min()).norm();

				assert(!isnan(current_diagonal));

				/*
				 * Generate any new desired spatial registers.
				 */

				/*
				 * Inputs
				 */

				for (int f = 0; f < 2; f++) {

					/*
					 * Low resolution
					 */

					if (current_diagonal < 2 * diagonal
					 && lowres == 0) {
						if (spatial_info_map.find(st.get_node()) == spatial_info_map.end()) {
							spatial_info_map[st.get_node()];
							ui::get()->d3_increment_spaces();
						}
						lowres = 1;
					}

					/*
					 * High resolution.
					 */

					if (current_diagonal < diagonal
					 && highres == 0) {
						if (spatial_info_map.find(st.get_node()) == spatial_info_map.end()) {
							spatial_info_map[st.get_node()];
							ui::get()->d3_increment_spaces();
						}
						highres = 1;
					}
				}

				/*
				 * Check for completion
				 */

				if (highres && lowres)
					return;

				/*
				 * Check precision before analyzing space further.
				 */

				if (st.precision_wall()) {
					fprintf(stderr, "\n\n*** Error: reached subspace precision wall ***\n\n");
					assert(0);
					return;
				}

				if (st.positive().includes(iw)) {
					st = st.positive();
					total_tsteps++;
				} else if (st.negative().includes(iw)) {
					st = st.negative();
					total_tsteps++;
				} else {
					fprintf(stderr, "failed iw = (%f, %f, %f)\n", 
							(double) iw[0], (double) iw[1], (double) iw[2]);
					assert(0);
				}
			}
		}

	public:
		candidates(int f) {

			image_index = f;
			height = (unsigned int) al->get(f)->get_t(0).unscaled_height();
			width = (unsigned int) al->get(f)->get_t(0).unscaled_width();

			/*
			 * Is this necessary?
			 */

			levels.resize(primary_decimation_upper - input_decimation_lower + 1);
			for (int l = input_decimation_lower; l <= primary_decimation_upper; l++) {
				levels[l - input_decimation_lower].resize((unsigned int) (floor(height / pow(2, l))
							       * floor(width / pow(2, l))
							       * pairwise_ambiguity),
						 std::pair<ale_pos, ale_real>(0, 0));
			}
		}

		/*
		 * Point p is expected to be in local projective coordinates.
		 */

		void add_candidate(point p, int tc, ale_pos score) {
			assert(tc <= primary_decimation_upper);
			assert(tc >= input_decimation_lower);
			assert(p[2] < 0);
			assert(score >= 0);

			int i = (unsigned int) floor(p[0] / (ale_pos) pow(2, tc));
			int j = (unsigned int) floor(p[1] / (ale_pos) pow(2, tc));

			int swidth  = (int) floor(width / pow(2, tc));

			assert(j < swidth);
			assert(i < (int) floor(height / pow(2, tc)));

			for (unsigned int k = 0; k < pairwise_ambiguity; k++) {
				std::pair<ale_pos, ale_real> *pk =
					&(levels[tc - input_decimation_lower][i * swidth * pairwise_ambiguity + j * pairwise_ambiguity + k]);

				if (pk->first != 0 && score >= (ale_pos) pk->second)
					continue;

				if (i == 1 && j == 1 && tc == 4)
					fprintf(stderr, "[ac p2=%f score=%f]\n", (double) p[2], (double) score);

				ale_pos tp = pk->first;
				ale_real tr = pk->second;

				pk->first = p[2];
				pk->second = score;

				p[2] = tp;
				score = tr;

				if (p[2] == 0)
					break;
			}
		}

		/*
		 * Generate subspaces for candidates.
		 */

		void generate_subspaces() {

			fprintf(stderr, "+");
			for (int l = input_decimation_lower; l <= primary_decimation_upper; l++) {
				unsigned int sheight = (unsigned int) floor(height / pow(2, l));
				unsigned int swidth  = (unsigned int) floor(width  / pow(2, l));

				for (unsigned int i = 0; i < sheight; i++)
				for (unsigned int j = 0; j < swidth; j++)
				for (unsigned int k = 0; k < pairwise_ambiguity; k++) {
					std::pair<ale_pos, ale_real> *pk =
						&(levels[l - input_decimation_lower]
							[i * swidth * pairwise_ambiguity + j * pairwise_ambiguity + k]);

					if (pk->first == 0) {
						fprintf(stderr, "o");
						continue;
					} else {
						fprintf(stderr, "|");
					}

					ale_pos si = i * pow(2, l) + ((l > 0) ? pow(2, l - 1) : 0);
					ale_pos sj = j * pow(2, l) + ((l > 0) ? pow(2, l - 1) : 0);

//					fprintf(stderr, "[gss l=%d i=%d j=%d d=%g]\n", l, i, j, pk->first);

					point p = al->get(image_index)->get_t(0).pw_unscaled(point(si, sj, pk->first));

					generate_subspace(p, 
							al->get(image_index)->get_t(0).diagonal_distance_3d(pk->first, l));
				}
			}
		}
	};

	/*
	 * List for calculating weighted median.
	 */
	class wml {
		ale_real *data;
		unsigned int size;
		unsigned int used;

		ale_real &_w(unsigned int i) {
			assert(i <= used);
			return data[i * 2];
		}

		ale_real &_d(unsigned int i) {
			assert(i <= used);
			return data[i * 2 + 1];
		}

		void increase_capacity() {

			if (size > 0)
				size *= 2;
			else
				size = 1;

			data = (ale_real *) realloc(data, sizeof(ale_real) * 2 * (size * 1));

			assert(data);
			assert (size > used);

			if (!data) {
				fprintf(stderr, "Unable to allocate %d bytes of memory\n",
						sizeof(ale_real) * 2 * (size * 1));
				exit(1);
			} 
		}

		void insert_weight(unsigned int i, ale_real v, ale_real w) {
			assert(used < size);
			assert(used >= i);
			for (unsigned int j = used; j > i; j--) {
				_w(j) = _w(j - 1);
				_d(j) = _d(j - 1);
			}

			_w(i) = w;
			_d(i) = v;

			used++;
		}

	public:

		unsigned int get_size() {
			return size;
		}

		unsigned int get_used() {
			return used;
		}

		void print_info() {
			fprintf(stderr, "[st %p sz %u el", this, size);
			for (unsigned int i = 0; i < used; i++)
				fprintf(stderr, " (%f, %f)", (double) _d(i), (double) _w(i));
			fprintf(stderr, "]\n");
		}

		void clear() {
			used = 0;
		}

		void insert_weight(ale_real v, ale_real w) {
			for (unsigned int i = 0; i < used; i++) {
				if (_d(i) == v) {
					_w(i) += w;
					return;
				}
				if (_d(i) > v) {
					if (used == size)
						increase_capacity();
					insert_weight(i, v, w);
					return;
				}
			}
			if (used == size)
				increase_capacity();
			insert_weight(used, v, w);
		}

		/*
		 * Finds the median at half-weight, or between half-weight
		 * and zero-weight, depending on the attenuation value.
		 */

		ale_real find_median(double attenuation = 0) {

			assert(attenuation >= 0);
			assert(attenuation <= 1);

			ale_real zero1 = 0;
			ale_real zero2 = 0;
			ale_real undefined = zero1 / zero2;

			ale_accum weight_sum = 0;

			for (unsigned int i = 0; i < used; i++)
				weight_sum += _w(i);

//			if (weight_sum == 0)
//				return undefined;

			if (used == 0 || used == 1)
				return undefined;

			if (weight_sum == 0) {
				ale_accum data_sum = 0;
				for (unsigned int i = 0; i < used; i++)
					data_sum += _d(i);
				return data_sum / (ale_accum) used;
			}
					

			ale_accum midpoint = weight_sum * (ale_accum) (0.5 - 0.5 * attenuation);

			ale_accum weight_sum_2 = 0;

			for (unsigned int i = 0; i < used && weight_sum_2 < midpoint; i++) {
				weight_sum_2 += _w(i);

				if (weight_sum_2 > midpoint) {
					return _d(i);
				} else if (weight_sum_2 == midpoint) {
					assert (i + 1 < used);
					return (_d(i) + _d(i + 1)) / 2;
				} 
			}

			return undefined;
		}

		wml(int initial_size = 0) {

//			if (initial_size == 0) {
//				initial_size = (int) (d2::image_rw::count() * 1.5);
//			}

			size = initial_size;
			used = 0;

			if (size > 0) {
				data = (ale_real *) malloc(size * sizeof(ale_real) * 2);
				assert(data);
			} else {
				data = NULL;
			}
		}

		/*
		 * copy constructor.  This is required to avoid undesired frees.
		 */

		wml(const wml &w) {
			size = w.size;
			used = w.used;
			data = (ale_real *) malloc(size * sizeof(ale_real) * 2);
			assert(data);

			memcpy(data, w.data, size * sizeof(ale_real) * 2);
		}

		~wml() {
			free(data);
		}
	};

	/*
	 * Class for information regarding spatial regions of interest.
	 *
	 * This class is configured for convenience in cases where sampling is
	 * performed using an approximation of the fine:box:1,triangle:2 chain.
	 * In this case, the *_1 variables would store the fine data and the
	 * *_2 variables would store the coarse data.  Other uses are also
	 * possible.
	 */

	class spatial_info {
		/*
		 * Map channel value --> weight.
		 */
		wml color_weights_1[3];
		wml color_weights_2[3];

		/*
		 * Current color.
		 */
		d2::pixel color;

		/*
		 * Map occupancy value --> weight.
		 */
		wml occupancy_weights_1;
		wml occupancy_weights_2;

		/*
		 * Current occupancy value.
		 */
		ale_real occupancy;

		/*
		 * pocc/socc density
		 */

		unsigned int pocc_density;
		unsigned int socc_density;

		/*
		 * Insert a weight into a list.
		 */
		void insert_weight(wml *m, ale_real v, ale_real w) {
			m->insert_weight(v, w);
		}

		/*
		 * Find the median of a weighted list.  Uses NaN for undefined.
		 */
		ale_real find_median(wml *m, double attenuation = 0) {
			return m->find_median(attenuation);
		}

	public:
		/*
		 * Constructor.
		 */
		spatial_info() {
			color = d2::pixel::zero();
			occupancy = 0;
			pocc_density = 0;
			socc_density = 0;
		}

		/*
		 * Accumulate color; primary data set.
		 */
		void accumulate_color_1(int f, d2::pixel color, d2::pixel weight) {
			for (int k = 0; k < 3; k++)
				insert_weight(&color_weights_1[k], color[k], weight[k]);
		}

		/*
		 * Accumulate color; secondary data set.
		 */
		void accumulate_color_2(d2::pixel color, d2::pixel weight) {
			for (int k = 0; k < 3; k++)
				insert_weight(&color_weights_2[k], color[k], weight[k]);
		}

		/*
		 * Accumulate occupancy; primary data set.
		 */
		void accumulate_occupancy_1(int f, ale_real occupancy, ale_real weight) {
			insert_weight(&occupancy_weights_1, occupancy, weight);
		}

		/*
		 * Accumulate occupancy; secondary data set.
		 */
		void accumulate_occupancy_2(ale_real occupancy, ale_real weight) {
			insert_weight(&occupancy_weights_2, occupancy, weight);
			
			if (occupancy == 0 || occupancy_weights_2.get_size() < 96)
				return;

			// fprintf(stderr, "%p updated socc with: %f %f\n", this, occupancy, weight);
			// occupancy_weights_2.print_info();
		}

		/*
		 * Update color (and clear accumulation structures).
		 */
		void update_color() {
			for (int d = 0; d < 3; d++) {
				ale_real c = find_median(&color_weights_1[d]);
				if (isnan(c))
					c = find_median(&color_weights_2[d]);
				if (isnan(c))
					c = 0;

				color[d] = c;

				color_weights_1[d].clear();
				color_weights_2[d].clear();
			}
		}

		/*
		 * Update occupancy (and clear accumulation structures).
		 */
		void update_occupancy() {
			ale_real o = find_median(&occupancy_weights_1, occ_att);
			if (isnan(o))
				o = find_median(&occupancy_weights_2, occ_att);
			if (isnan(o))
				o = 0;

			occupancy = o;

			pocc_density = occupancy_weights_1.get_used();
			socc_density = occupancy_weights_2.get_used();

			occupancy_weights_1.clear();
			occupancy_weights_2.clear();

		}

		/*
		 * Get current color.
		 */
		d2::pixel get_color() {
			return color;
		}

		/*
		 * Get current occupancy.
		 */
		ale_real get_occupancy() {
			assert (finite(occupancy));
			return occupancy;
		}

		/*
		 * Get primary color density.
		 */

		unsigned int get_pocc_density() {
			return pocc_density;
		}

		unsigned int get_socc_density() {
			return socc_density;
		}
	};

	/*
	 * Map spatial regions of interest to spatial info structures.  XXX:
	 * This may get very poor cache behavior in comparison with, say, an
	 * array.  Unfortunately, there is no immediately obvious array
	 * representation.  If some kind of array representation were adopted,
	 * it would probably cluster regions of similar depth from the
	 * perspective of the typical camera.  In particular, for a
	 * stereoscopic view, depth ordering for two random points tends to be
	 * similar between cameras, I think.  Unfortunately, it is never
	 * identical for all points (unless cameras are co-located).  One
	 * possible approach would be to order based on, say, camera 0's idea
	 * of depth.
	 */

#if !defined(HASH_MAP_GNU) && !defined(HASH_MAP_STD)
	typedef std::map<struct space::node *, spatial_info> spatial_info_map_t;
#elif defined(HASH_MAP_GNU)
        struct node_hash
        {
                size_t operator()(struct space::node *n) const
                {
                        return __gnu_cxx::hash<long>()((long) n);
                }
        };
	typedef __gnu_cxx::hash_map<struct space::node *, spatial_info, node_hash > spatial_info_map_t;
#elif defined(HASH_MAP_STD)
	typedef std::hash_map<struct space::node *, spatial_info> spatial_info_map_t;
#endif

	static spatial_info_map_t spatial_info_map;

public:

	/*
	 * Debugging variables.
	 */

	static unsigned long total_ambiguity;
	static unsigned long total_pixels;
	static unsigned long total_divisions;
	static unsigned long total_tsteps;

	/*
	 * Member functions
	 */

	static void et(double et_parameter) {
		encounter_threshold = et_parameter;
	}

	static void dmr(double dmr_parameter) {
		depth_median_radius = dmr_parameter;
	}

	static void fmr(double fmr_parameter) {
		diff_median_radius = fmr_parameter;
	}

	static void load_model(const char *name) {
		load_model_name = name;
	}

	static void save_model(const char *name) {
		save_model_name = name;
	}

	static void fc(ale_pos fc) {
		front_clip = fc;
	}

	static void di_upper(ale_pos _dgi) {
		primary_decimation_upper = (int) round(_dgi);
	}

	static void do_try(ale_pos _dgo) {
		output_decimation_preferred = (int) round(_dgo);
	}

	static void di_lower(ale_pos _idiv) {
		input_decimation_lower = (int) round(_idiv);
	}

	static void oc() {
		output_clip = 1;
	}

	static void no_oc() {
		output_clip = 0;
	}

	static void rc(ale_pos rc) {
		rear_clip = rc;
	}

	/*
	 * Initialize 3D scene from 2D scene, using 2D and 3D alignment
	 * information.
	 */
	static void init_from_d2() {

		/*
		 * Rear clip value of 0 is converted to infinity.
		 */

		if (rear_clip == 0) {
			ale_pos one = +1;
			ale_pos zero = +0;

			rear_clip = one / zero;
			assert(isinf(rear_clip) && rear_clip > 0);
		}

		/*
		 * Scale and translate clipping plane depths.
		 */

		ale_pos cp_scalar = d3::align::projective(0).wc(point(0, 0, 0))[2];

		front_clip = front_clip * cp_scalar - cp_scalar;
		rear_clip = rear_clip * cp_scalar - cp_scalar;

		/*
		 * Allocate image structures.
		 */

		al = new lod_images;

		if (tc_multiplier != 0) {
			al->open_all();
		}
	}

	/*
	 * Perform spatial_info updating on a given subspace, for given
	 * parameters.
	 */
	static void subspace_info_update(space::iterate si, int f, ref_weights *weights) {

		while(!si.done()) {

			space::traverse st = si.get();

			/*
			 * Reject out-of-bounds spaces.
			 */
			if (!weights->in_spatial_bounds(st)) {
				si.next();
				continue;
			}
				
			/*
			 * Skip spaces with no color information.
			 *
			 * XXX: This could be more efficient, perhaps.
			 */

			if (spatial_info_map.count(st.get_node()) == 0) {
				si.next();
				continue;
			}

			ui::get()->d3_increment_space_num();


			/*
			 * Get in-bounds centroid, if one exists.
			 */

			point p = al->get(f)->get_t(0).centroid(st);

			if (!p.defined()) {
				si.next();
				continue;
			}
				
			/*
			 * Get information on the subspace.
			 */

			spatial_info *sn = &spatial_info_map[st.get_node()];
			d2::pixel color = sn->get_color();
			ale_real occupancy = sn->get_occupancy();

			/*
			 * Store current weight so we can later check for
			 * modification by higher-resolution subspaces.
			 */

			ref_weights::subtree *tree = weights->get_subtree(st);

			/*
			 * Check for higher resolution subspaces, and
			 * update the space iterator.
			 */

			if (st.get_node()->positive
			 || st.get_node()->negative) {

				/*
				 * Cleave space for the higher-resolution pass,
				 * skipping the current space, since we will
				 * process that later.
				 */

				space::iterate cleaved_space = si.cleave();

				cleaved_space.next();

				subspace_info_update(cleaved_space, f, weights);

			} else {
				si.next();
			}

			/*
			 * Add new data on the subspace and update weights.
			 */

			ale_pos tc = al->get(f)->get_t(0).trilinear_coordinate(st);
			d2::pixel pcolor = al->get(f)->get_tl(p.xy(), tc);
			d2::pixel colordiff = (color - pcolor) * (ale_real) 256;

			if (falloff_exponent != 0) {
				d2::pixel max_diff = al->get(f)->get_max_diff(p.xy(), tc) * (ale_real) 256;

				for (int k = 0; k < 3; k++)
				if (max_diff[k] > 1)
					colordiff[k] /= pow(max_diff[k], falloff_exponent);
			}

			/*
			 * Determine the probability of encounter.
			 */

			d2::pixel encounter = d2::pixel(1, 1, 1) * (1 - weights->get_weight(st));

			/*
			 * Update weights
			 */

			weights->add_weight(st, occupancy, tree);

			/*
			 * Delete the subtree, if necessary.
			 */

			delete tree;

			/*
			 * Check for cases in which the subspace should not be
			 * updated.
			 */

			if (!resolution_ok(al->get(f)->get_t(0), tc))
				continue;

			if (d2::render::is_excluded_f(p.xy(), f))
				continue;

			/*
			 * Update subspace.
			 */

			sn->accumulate_color_1(f, pcolor, encounter);
			d2::pixel channel_occ = pexp(-colordiff * colordiff);

			ale_real occ = channel_occ[0];

			for (int k = 1; k < 3; k++)
				if (channel_occ[k] < occ)
					occ = channel_occ[k];

			sn->accumulate_occupancy_1(f, occ, encounter[0]);

		}
	}

	/*
	 * Run a single iteration of the spatial_info update cycle.
	 */
	static void spatial_info_update() {
		/*
		 * Iterate through each frame.
		 */
		for (unsigned int f = 0; f < d2::image_rw::count(); f++) {

			ui::get()->d3_occupancy_status(f);

			/*
			 * Open the frame and transformation.
			 */

			if (tc_multiplier == 0)
				al->open(f);

			/*
			 * Allocate weights data structure for storing encounter
			 * probabilities.
			 */

			ref_weights *weights = new ref_weights(f);

			/*
			 * Call subspace_info_update for the root space.
			 */

			subspace_info_update(space::iterate(al->get(f)->origin()), f, weights);

			/*
			 * Free weights.
			 */

			delete weights;

			/*
			 * Close the frame and transformation.
			 */

			if (tc_multiplier == 0)
				al->close(f);
		}

		/*
		 * Update all spatial_info structures.
		 */
		for (spatial_info_map_t::iterator i = spatial_info_map.begin(); i != spatial_info_map.end(); i++) {
			i->second.update_color();
			i->second.update_occupancy();

//			d2::pixel color = i->second.get_color();

//			fprintf(stderr, "space p=%p updated to c=[%f %f %f] o=%f\n",
//					i->first, color[0], color[1], color[2], 
//					i->second.get_occupancy());
		}
	}

	/*
	 * Support function for view() and depth().  This function
	 * always performs exclusion.
	 */

	static const void view_recurse(int type, d2::image *im, d2::image *weights, space::iterate si, pt _pt,
			int prune = 0, d2::point pl = d2::point(0, 0), d2::point ph = d2::point(0, 0)) {
		while (!si.done()) {
			space::traverse st = si.get();

			/*
			 * Remove excluded regions.
			 */

			if (excluded(st)) {
				si.cleave();
				continue;
			}

			/*
			 * Prune.
			 */

			if (prune && !_pt.check_inclusion_scaled(st, pl, ph)) {
				si.cleave();
				continue;
			}

			/*
			 * XXX: This could be more efficient, perhaps.
			 */

			if (spatial_info_map.count(st.get_node()) == 0) {
				si.next();
				continue;
			}

			ui::get()->d3_increment_space_num();

			spatial_info sn = spatial_info_map[st.get_node()];

			/*
			 * Get information on the subspace.
			 */

			d2::pixel color = sn.get_color();
			// d2::pixel color = d2::pixel(1, 1, 1) * (double) (((unsigned int) (st.get_node()) / sizeof(space)) % 65535);
			ale_real occupancy = sn.get_occupancy();

			/*
			 * Determine the view-local bounding box for the
			 * subspace.
			 */

			point bb[2];

			_pt.get_view_local_bb_scaled(st, bb);

			point min = bb[0];
			point max = bb[1];

			if (prune) {
				if (min[0] > ph[0]
				 || min[1] > ph[1]
				 || max[0] < pl[0]
				 || max[1] < pl[1]) {
					si.next();
					continue;
				}

				if (min[0] < pl[0])
					min[0] = pl[0];
				if (min[1] < pl[1])
					min[1] = pl[1];
				if (max[0] > ph[0])
					max[0] = ph[0];
				if (max[1] > ph[1])
					max[1] = ph[1];

				min[0] -= pl[0];
				min[1] -= pl[1];
				max[0] -= pl[0];
				max[1] -= pl[1];
			}

			/*
			 * Data structure to check modification of weights by
			 * higher-resolution subspaces.
			 */

			std::queue<d2::pixel> weight_queue;

			/*
			 * Check for higher resolution subspaces, and
			 * update the space iterator.
			 */

			if (st.get_node()->positive
			 || st.get_node()->negative) {

				/*
				 * Store information about current weights,
				 * so we will know which areas have been
				 * covered by higher-resolution subspaces.
				 */

				for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
				for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++)
					weight_queue.push(weights->get_pixel(i, j));
				
				/*
				 * Cleave space for the higher-resolution pass,
				 * skipping the current space, since we will
				 * process that afterward.
				 */

				space::iterate cleaved_space = si.cleave();

				cleaved_space.next();

				view_recurse(type, im, weights, cleaved_space, _pt, prune, pl, ph);

			} else {
				si.next();
			}
				

			/*
			 * Iterate over pixels in the bounding box, finding
			 * pixels that intersect the subspace.  XXX: assume
			 * for now that all pixels in the bounding box
			 * intersect the subspace.
			 */

			for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
			for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++) {

				/*
				 * Check for higher-resolution updates.
				 */

				if (weight_queue.size()) {
					if (weight_queue.front() != weights->get_pixel(i, j)) {
						weight_queue.pop();
						continue;
					}
					weight_queue.pop();
				}

				/*
				 * Determine the probability of encounter.
				 */

				d2::pixel encounter = (d2::pixel(1, 1, 1) 
						     - weights->get_pixel(i, j)) 
					            * occupancy;

				/*
				 * Update images.
				 */

				if (type == 0) {

					/*
					 * Color view
					 */

					weights->set_pixel(i, j, (d2::pixel) weights->get_pixel(i, j) 
					                       + encounter);
					im->set_pixel(i, j, (d2::pixel) im->get_pixel(i, j)
					                  + encounter * color);

				} else if (type == 1) {

					/*
					 * Weighted (transparent) depth display
					 */

					ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
					weights->set_pixel(i, j, (d2::pixel) weights->get_pixel(i, j)
					                       + encounter);
					im->set_pixel(i, j, (d2::pixel) im->get_pixel(i, j)
					                  + encounter * (ale_real) depth_value);

				} else if (type == 2) {

					/*
					 * Ambiguity (ambivalence) measure.
					 */

					weights->set_pixel(i, j, d2::pixel(1, 1, 1));
					im->set_pixel(i, j, (d2::pixel) im->get_pixel(i, j)
					                  + 0.1 * d2::pixel(1, 1, 1));

				} else if (type == 3) {

					/*
					 * Closeness measure.
					 */

					ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
					if (weights->get_chan(i, j, 0) == 0) {
						weights->set_pixel(i, j, d2::pixel(1, 1, 1));
						im->set_pixel(i, j, d2::pixel(1, 1, 1)
						                  * (ale_real) depth_value);
					} else if (im->get_chan(i, j, 2) < (ale_sreal) depth_value) {
						im->set_pixel(i, j, d2::pixel(1, 1, 1) 
						              * (ale_real) depth_value);
					} else {
						continue;
					}

				} else if (type == 4) {

					/*
					 * Weighted (transparent) contribution display
					 */

					ale_pos contribution_value = sn.get_pocc_density() /* + sn.get_socc_density() */;
					weights->set_pixel(i, j, (d2::pixel) weights->get_pixel(i, j)
					                       + encounter);
					im->set_pixel(i, j, (d2::pixel) im->get_pixel(i, j) 
					                  + encounter * (ale_real) contribution_value);

					assert (finite(encounter[0]));
					assert (finite(contribution_value));

				} else if (type == 5) {

					/*
					 * Weighted (transparent) occupancy display
					 */

					ale_real contribution_value = occupancy;
					weights->set_pixel(i, j, (d2::pixel) weights->get_pixel(i, j)
					                       + encounter);
					im->set_pixel(i, j, (d2::pixel) im->get_pixel(i, j)
					                  + encounter * contribution_value);

				} else if (type == 6) {
					
					/*
					 * (Depth, xres, yres) triple
					 */

					ale_pos depth_value = _pt.wp_scaled(st.get_min())[2];
					weights->set_chan(i, j, 0, weights->get_chan(i, j, 0)
					                         + encounter[0]);
					if (weights->get_pixel(i, j)[1] < encounter[0]) {
						weights->set_chan(i, j, 1, encounter[0]);
						im->set_pixel(i, j, d2::pixel(
							weights->get_pixel(i, j)[1] * (ale_real) depth_value,
							ale_pos_to_real(max[0] - min[0]),
							ale_pos_to_real(max[1] - min[1])));
					}

				} else if (type == 7) {
					
					/*
					 * (xoff, yoff, 0) triple
					 */

					weights->set_chan(i, j, 0,
						weights->get_chan(i, j, 0) + encounter[0]);
					if (weights->get_chan(i, j, 1) < (ale_sreal) encounter[0]) {
						weights->set_chan(i, j, 1, encounter[0]);
						im->set_pixel(i, j, d2::pixel(
							ale_pos_to_real(i - min[0]),
							ale_pos_to_real(j - min[1]),
							0));
					}

				} else if (type == 8) {

					/*
					 * Value = 1 for any intersected space.
					 */

					weights->set_pixel(i, j, d2::pixel(1, 1, 1));
					im->set_pixel(i, j, d2::pixel(1, 1, 1));

				} else if (type == 9) {

					/*
					 * Number of contributions for the nearest space.
					 */

					if (weights->get_chan(i, j, 0) == 1)
						continue;

					weights->set_pixel(i, j, d2::pixel(1, 1, 1));
					im->set_pixel(i, j, d2::pixel(1, 1, 1) * (sn.get_pocc_density() * 0.1));

				} else 
					assert(0);
			}
		}
	}

	/*
	 * Generate an depth image from a specified view.
	 */
	static const d2::image *depth(pt _pt, int n = -1, int prune = 0, 
			d2::point pl = d2::point(0, 0), d2::point ph = d2::point(0, 0)) {
		assert ((unsigned int) n < d2::image_rw::count() || n < 0);

		_pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);

		if (n >= 0) {
			assert((int) floor(d2::align::of(n).scaled_height())
			     == (int) floor(_pt.scaled_height()));
			assert((int) floor(d2::align::of(n).scaled_width())
			     == (int) floor(_pt.scaled_width()));
		}

		d2::image *im1, *im2, *im3, *weights;;

		if (prune) {

			im1 = d2::new_image_ale_real((int) floor(ph[0] - pl[0]) + 1,
					(int) floor(ph[1] - pl[1]) + 1, 3);

			im2 = d2::new_image_ale_real((int) floor(ph[0] - pl[0]) + 1,
					(int) floor(ph[1] - pl[1]) + 1, 3);

			im3 = d2::new_image_ale_real((int) floor(ph[0] - pl[0]) + 1,
					(int) floor(ph[1] - pl[1]) + 1, 3);

			weights = d2::new_image_ale_real((int) floor(ph[0] - pl[0]) + 1,
					(int) floor(ph[1] - pl[1]) + 1, 3);

		} else {

			im1 = d2::new_image_ale_real((int) floor(_pt.scaled_height()),
			     	       (int) floor(_pt.scaled_width()), 3);

			im2 = d2::new_image_ale_real((int) floor(_pt.scaled_height()),
			     	       (int) floor(_pt.scaled_width()), 3);

			im3 = d2::new_image_ale_real((int) floor(_pt.scaled_height()),
			     	       (int) floor(_pt.scaled_width()), 3);

			weights = d2::new_image_ale_real((int) floor(_pt.scaled_height()),
							(int) floor(_pt.scaled_width()), 3);
		}

		/*
		 * Iterate through subspaces.
		 */

		space::iterate si(_pt.origin());

		view_recurse(6, im1, weights, si, _pt, prune, pl, ph);

		delete weights;

		if (prune) {
			weights = d2::new_image_ale_real((int) floor(ph[0] - pl[0]) + 1,
					(int) floor(ph[1] - pl[1]) + 1, 3);
		} else {
			weights = d2::new_image_ale_real((int) floor(_pt.scaled_height()),
							(int) floor(_pt.scaled_width()), 3);
		}

#if 1
		view_recurse(7, im2, weights, si, _pt, prune, pl, ph);
#else
		view_recurse(8, im2, weights, si, _pt, prune, pl, ph);
		return im2;
#endif

		/*
		 * Normalize depths by weights
		 */

		if (normalize_weights)
		for (unsigned int i = 0; i < im1->height(); i++)
		for (unsigned int j = 0; j < im1->width();  j++)
			im1->set_chan(i, j, 0, im1->get_chan(i, j, 0) / weights->get_chan(i, j, 1));

	
		for (unsigned int i = 0; i < im1->height(); i++)
		for (unsigned int j = 0; j < im1->width();  j++) {

			/*
			 * Handle interpolation.
			 */

			d2::point x;
			d2::point blx;
			d2::point res((double) im1->get_chan(i, j, 1), 
			              (double) im1->get_chan(i, j, 2));

			for (int d = 0; d < 2; d++) {

				if (im2->get_chan(i, j, d) < (ale_sreal) res[d] / 2)
					x[d] = (ale_pos) (d?j:i) - res[d] / 2 - (ale_pos) im2->get_chan(i, j, d);
				else
					x[d] = (ale_pos) (d?j:i) + res[d] / 2 - (ale_pos) im2->get_chan(i, j, d);

				blx[d] = 1 - ((d?j:i) - x[d]) / res[d];
			}

			ale_real depth_val = 0;
			ale_real depth_weight = 0;

			for (int ii = 0; ii < 2; ii++)
			for (int jj = 0; jj < 2; jj++) {
				d2::point p = x + d2::point(ii, jj) * res;
				if (im1->in_bounds(p)) {

					ale_real d = im1->get_bl(p)[0];

					if (isnan(d))
						continue;

					ale_real w = ale_pos_to_real((ii ? (1 - blx[0]) : blx[0]) * (jj ? (1 - blx[1]) : blx[1]));
					depth_weight += w;
					depth_val += w * d;
				}
			}

			ale_real depth = depth_val / depth_weight;

			/*
			 * Handle encounter thresholds
			 */

			if (weights->get_chan(i, j, 0) < encounter_threshold) {
				im3->set_pixel(i, j, d2::pixel::zero() / d2::pixel::zero());
			} else {
				im3->set_pixel(i, j, d2::pixel(1, 1, 1) * depth);
			}
		}

		delete weights;
		delete im1;
		delete im2;

		return im3;
	}

	static const d2::image *depth(unsigned int n) {

		assert (n < d2::image_rw::count());

		pt _pt = align::projective(n);

		return depth(_pt, n);
	}


	/*
	 * This function always performs exclusion.
	 */

	static space::node *most_visible_pointwise(d2::pixel *weight, space::iterate si, pt _pt, d2::point p) {

		space::node *result = NULL;

		while (!si.done()) {
			space::traverse st = si.get();

			/*
			 * Prune certain regions known to be uninteresting.
			 */

			if (excluded(st) || !_pt.check_inclusion_scaled(st, p)) {
				si.cleave();
				continue;
			}

			/*
			 * XXX: This could be more efficient, perhaps.
			 */

			if (spatial_info_map.count(st.get_node()) == 0) {
				si.next();
				continue;
			}

			spatial_info sn = spatial_info_map[st.get_node()];

			/*
			 * Get information on the subspace.
			 */

			ale_real occupancy = sn.get_occupancy();

			/*
			 * Preserve current weight in order to check for
			 * modification by higher-resolution subspaces.
			 */

			d2::pixel old_weight = *weight;

			/*
			 * Check for higher resolution subspaces, and
			 * update the space iterator.
			 */

			if (st.get_node()->positive
			 || st.get_node()->negative) {

				/*
				 * Cleave space for the higher-resolution pass,
				 * skipping the current space, since we will
				 * process that afterward.
				 */

				space::iterate cleaved_space = si.cleave();

				cleaved_space.next();

				space::node *r = most_visible_pointwise(weight, cleaved_space, _pt, p);

				if (old_weight[1] != (*weight)[1])
					result = r;

			} else {
				si.next();
			}
				

			/*
			 * Check for higher-resolution updates.
			 */

			if (old_weight != *weight)
				continue;

			/*
			 * Determine the probability of encounter.
			 */

			ale_real encounter = (1 - (*weight)[0]) * occupancy;

			/*
			 * (*weight)[0] stores the cumulative weight; (*weight)[1] stores the maximum.
			 */

			if (encounter > (*weight)[1]) {
				result = st.get_node();
				(*weight)[1] = encounter;
			}

			(*weight)[0] += encounter;
		}

		return result;
	}

	/*
	 * This function performs exclusion iff SCALED is true.
	 */
	static void  most_visible_generic(std::vector<space::node *> &results, d2::image *weights, 
			space::iterate si, pt _pt, int scaled) {

		assert (results.size() == weights->height() * weights->width());

		while (!si.done()) {
			space::traverse st = si.get();

			if (scaled && excluded(st)) {
				si.cleave();
				continue;
			}

			/*
			 * XXX: This could be more efficient, perhaps.
			 */

			if (spatial_info_map.count(st.get_node()) == 0) {
				si.next();
				continue;
			}

			spatial_info sn = spatial_info_map[st.get_node()];

			/*
			 * Get information on the subspace.
			 */

			ale_real occupancy = sn.get_occupancy();

			/*
			 * Determine the view-local bounding box for the
			 * subspace.
			 */

			point bb[2];

			_pt.get_view_local_bb_scaled(st, bb);

			point min = bb[0];
			point max = bb[1];

			/*
			 * Data structure to check modification of weights by
			 * higher-resolution subspaces.
			 */

			std::queue<d2::pixel> weight_queue;

			/*
			 * Check for higher resolution subspaces, and
			 * update the space iterator.
			 */

			if (st.get_node()->positive
			 || st.get_node()->negative) {

				/*
				 * Store information about current weights,
				 * so we will know which areas have been
				 * covered by higher-resolution subspaces.
				 */

				for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
				for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++)
					weight_queue.push(weights->get_pixel(i, j));
				
				/*
				 * Cleave space for the higher-resolution pass,
				 * skipping the current space, since we will
				 * process that afterward.
				 */

				space::iterate cleaved_space = si.cleave();

				cleaved_space.next();

				most_visible_generic(results, weights, cleaved_space, _pt, scaled);

			} else {
				si.next();
			}
				

			/*
			 * Iterate over pixels in the bounding box, finding
			 * pixels that intersect the subspace.  XXX: assume
			 * for now that all pixels in the bounding box
			 * intersect the subspace.
			 */

			for (int i = (int) ceil(min[0]); i <= (int) floor(max[0]); i++)
			for (int j = (int) ceil(min[1]); j <= (int) floor(max[1]); j++) {

				/*
				 * Check for higher-resolution updates.
				 */

				if (weight_queue.size()) {
					if (weight_queue.front() != weights->get_pixel(i, j)) {
						weight_queue.pop();
						continue;
					}
					weight_queue.pop();
				}

				/*
				 * Determine the probability of encounter.
				 */

				ale_real encounter = (1 - weights->get_pixel(i, j)[0]) * occupancy;

				/*
				 * weights[0] stores the cumulative weight; weights[1] stores the maximum.
				 */

				if (encounter > weights->get_pixel(i, j)[1]
				 || results[i * weights->width() + j] == NULL) {
					results[i * weights->width() + j] = st.get_node();
					weights->set_chan(i, j, 1, encounter);
				}

				weights->set_chan(i, j, 0, weights->get_chan(i, j, 0) + encounter);
			}
		}
	}

	static std::vector<space::node *> most_visible_scaled(pt _pt) {
		d2::image *weights = d2::new_image_ale_real((int) floor(_pt.scaled_height()), 
				(int) floor(_pt.scaled_width()), 3);
		std::vector<space::node *> results;

		results.resize(weights->height() * weights->width(), 0);
	
		most_visible_generic(results, weights, space::iterate(_pt.origin()), _pt, 1);
		
		return results;
	}

	static std::vector<space::node *> most_visible_unscaled(pt _pt) {
		d2::image *weights = d2::new_image_ale_real((int) floor(_pt.unscaled_height()), 
				(int) floor(_pt.unscaled_width()), 3);
		std::vector<space::node *> results;
		
		results.resize(weights->height() * weights->width(), 0);

		most_visible_generic(results, weights, space::iterate(_pt.origin()), _pt, 0);
		
		return results;
	}

	static const int visibility_search(const std::vector<space::node *> &fmv, space::node *mv) {

		if (mv == NULL)
			return 0;

		if (std::binary_search(fmv.begin(), fmv.end(), mv))
			return 1;

		return (visibility_search(fmv, mv->positive)
		     || visibility_search(fmv, mv->negative));

	}

	/*
	 * Class to generate focal sample views.
	 */

	class view_generator {

		/*
		 * Original projective transformation.
		 */

		pt original_pt;

		/*
		 * Data type for shared view data.
		 */

		class shared_view {
			pt _pt;
			std::vector<space::node *> mv;
			d2::image *color;
			d2::image *color_weights;
			const d2::image *_depth;
			d2::image *median_depth;
			d2::image *median_diff;

		public:
			shared_view(pt _pt) {
				this->_pt = _pt;
				color = NULL;
				color_weights = NULL;
				_depth = NULL;
				median_depth = NULL;
				median_diff = NULL;
			}

			shared_view(const shared_view &copy_origin) {
				_pt = copy_origin._pt;
				mv = copy_origin.mv;
				color = NULL;
				color_weights = NULL;
				_depth = NULL;
				median_depth = NULL;
				median_diff = NULL;
			}

			~shared_view() {
				delete color;
				delete _depth;
				delete color_weights;
				delete median_diff;
				delete median_depth;
			}

			void get_view_recurse(d2::image *data, d2::image *weights, int type) {
				/*
				 * Iterate through subspaces.
				 */

				space::iterate si(_pt.origin());

				ui::get()->d3_render_status(0, 0, -1, -1, -1, -1, 0);

				view_recurse(type, data, weights, si, _pt);
			}

			void init_color() {
				color = d2::new_image_ale_real((int) floor(_pt.scaled_height()),
								(int) floor(_pt.scaled_width()), 3);

				color_weights = d2::new_image_ale_real((int) floor(_pt.scaled_height()),
								(int) floor(_pt.scaled_width()), 3);

				get_view_recurse(color, color_weights, 0);
			}

			void init_depth() {
				_depth = depth(_pt, -1);
			}

			void init_medians() {
				if (!_depth)
					init_depth();

				assert(_depth);

				median_diff = _depth->fcdiff_median((int) floor(diff_median_radius));
				median_depth = _depth->medians((int) floor(depth_median_radius));

				assert(median_diff);
				assert(median_depth);
			}

		public:
			pt get_pt() {
				return _pt;
			}

			space::node *get_most_visible(unsigned int i, unsigned int j) {
				unsigned int height = (int) floor(_pt.scaled_height());
				unsigned int width  = (int) floor(_pt.scaled_width());

				if (i >= height
				 || j >= width) {
					return NULL;
				}

				if (mv.size() == 0) {
					mv = most_visible_scaled(_pt);
				}

				assert (mv.size() > i * width + j);

				return mv[i * width + j];
			}

			space::node *get_most_visible(d2::point p) {
				unsigned int i = (unsigned int) round (p[0]);
				unsigned int j = (unsigned int) round (p[1]);

				return get_most_visible(i, j);
			}

			d2::pixel get_color(unsigned int i, unsigned int j) {
				if (color == NULL) {
					init_color();
				}

				assert (color != NULL);

				return color->get_pixel(i, j);
			}

			d2::pixel get_depth(unsigned int i, unsigned int j) {
				if (_depth == NULL) {
					init_depth();
				}

				assert (_depth != NULL);

				return _depth->get_pixel(i, j);
			}

			void get_median_depth_and_diff(d2::pixel *t, d2::pixel *f, unsigned int i, unsigned int j) {
				if (median_depth == NULL && median_diff == NULL) 
					init_medians();

				assert (median_depth && median_diff);

				if (i >= median_depth->height()
				 || j >= median_depth->width()) {
					*t = d2::pixel::undefined();
					*f = d2::pixel::undefined();
				} else {
					*t = median_depth->get_pixel(i, j);
					*f = median_diff->get_pixel(i, j);
				}
			}

			void get_color_and_weight(d2::pixel *c, d2::pixel *w, d2::point p) {
				if (color == NULL) {
					init_color();
				}

				assert (color != NULL);

				if (!color->in_bounds(p)) {
					*c = d2::pixel::undefined();
					*w = d2::pixel::undefined();
				} else {
					*c = color->get_bl(p);
					*w = color_weights->get_bl(p);
				}
			}

			d2::pixel get_depth(d2::point p) {
				if (_depth == NULL) {
					init_depth();
				}

				assert (_depth != NULL);

				if (!_depth->in_bounds(p)) {
					return d2::pixel::undefined();
				}

				return _depth->get_bl(p);
			}

			void get_median_depth_and_diff(d2::pixel *t, d2::pixel *f, d2::point p) {
				if (median_diff == NULL && median_depth == NULL)
					init_medians();

				assert (median_diff != NULL && median_depth != NULL);

				if (!median_diff->in_bounds(p)) {
					*t = d2::pixel::undefined();
					*f = d2::pixel::undefined();
				} else {
					*t = median_depth->get_bl(p);
					*f = median_diff->get_bl(p);
				}
			}

		};

		/*
		 * Shared view array, indexed by aperture diameter and view number.
		 */

		std::map<ale_pos, std::vector<shared_view> > aperture_to_shared_views_map;

		/*
		 * Method to generate a new stochastic focal view. 
		 */

		pt get_new_view(ale_pos aperture) {

			ale_pos ofx = aperture;
			ale_pos ofy = aperture;

			while (ofx * ofx + ofy * ofy > aperture * aperture / 4) {
				ofx = (rand() * aperture) / RAND_MAX - aperture / 2;
				ofy = (rand() * aperture) / RAND_MAX - aperture / 2;
			}

			/*
			 * Generate a new view from the given offset.
			 */

			point new_view = original_pt.cw(point(ofx, ofy, 0));
			pt _pt_new = original_pt;
			for (int d = 0; d < 3; d++)
				_pt_new.e().set_translation(d, -new_view[d]);

			return _pt_new;
		}

	public:

		/*
		 * Result type.
		 */

		class view {
			shared_view *sv;
			pt _pt;

		public:

			view(shared_view *sv, pt _pt = pt()) {
				this->sv = sv;
				if (sv) {
					this->_pt = sv->get_pt();
				} else {
					this->_pt = _pt;
				}
			}

			pt get_pt() {
				return _pt;
			}

			space::node *get_most_visible(unsigned int i, unsigned int j) {
				assert (sv);
				return sv->get_most_visible(i, j);
			}

			space::node *get_most_visible(d2::point p) {
				if (sv) {
					return sv->get_most_visible(p);
				}

				d2::pixel weight(0, 0, 0);

				return most_visible_pointwise(&weight, space::iterate(_pt.origin()), _pt, p);

			}

			d2::pixel get_color(unsigned int i, unsigned int j) {
				return sv->get_color(i, j);
			}

			void get_color_and_weight(d2::pixel *color, d2::pixel *weight, d2::point p) {
				if (sv) {
					sv->get_color_and_weight(color, weight, p);
					return;
				}

				/*
				 * Determine weight and color for the given point.
				 */

				d2::image *im_point = d2::new_image_ale_real(1, 1, 3);
				d2::image *wt_point = d2::new_image_ale_real(1, 1, 3);

				view_recurse(0, im_point, wt_point, space::iterate(_pt.origin()), _pt, 1, p, p);

				*color = im_point->get_pixel(0, 0);
				*weight = wt_point->get_pixel(0, 0);

				delete im_point;
				delete wt_point;

				return;
			}

			d2::pixel get_depth(unsigned int i, unsigned int j) {
				assert(sv);
				return sv->get_depth(i, j);
			}

			void get_median_depth_and_diff(d2::pixel *depth, d2::pixel *diff, unsigned int i, unsigned int j) {
				assert(sv);
				sv->get_median_depth_and_diff(depth, diff, i, j);
			}

			void get_median_depth_and_diff(d2::pixel *_depth, d2::pixel *_diff, d2::point p) {
				if (sv) {
					sv->get_median_depth_and_diff(_depth, _diff, p);
					return;
				}

				/*
				 * Generate a local depth image of required radius.
				 */

				ale_pos radius = 1;

				if (diff_median_radius + 1 > radius)
					radius = diff_median_radius + 1;
				if (depth_median_radius > radius)
					radius = depth_median_radius;

				d2::point pl = p - d2::point(radius, radius);
				d2::point ph = p + d2::point(radius, radius);
				const d2::image *local_depth = depth(_pt, -1, 1, pl, ph);

				/*
				 * Find depth and diff at this point, check for
				 * undefined values, and generate projections
				 * of the image corners on the estimated normal
				 * surface.
				 */

				d2::image *median_diffs = local_depth->fcdiff_median((int) floor(diff_median_radius));
				d2::image *median_depths = local_depth->medians((int) floor(depth_median_radius));

				*_depth = median_depths->get_pixel((int) radius, (int) radius);
				*_diff = median_diffs->get_pixel((int) radius, (int) radius);

				delete median_diffs;
				delete median_depths;
				delete local_depth;
			}
		};

		view get_view(ale_pos aperture, unsigned index, unsigned int randomization) {
			if (randomization == 0) {

				while (aperture_to_shared_views_map[aperture].size() <= index) {
					aperture_to_shared_views_map[aperture].push_back(shared_view(get_new_view(aperture)));
				}

				return view(&(aperture_to_shared_views_map[aperture][index]));
			}

			return view(NULL, get_new_view(aperture));
		}

		view_generator(pt original_pt) {
			this->original_pt = original_pt;
		}
	};

	/*
	 * Unfiltered function
	 */
	static const d2::image *view_nofilter_focus(pt _pt, int n) {

		assert ((unsigned int) n < d2::image_rw::count() || n < 0);

		if (n >= 0) {
			assert((int) floor(d2::align::of(n).scaled_height())
			     == (int) floor(_pt.scaled_height()));
			assert((int) floor(d2::align::of(n).scaled_width())
			     == (int) floor(_pt.scaled_width()));
		}

		const d2::image *depths = depth(_pt, n);

		d2::image *im = d2::new_image_ale_real((int) floor(_pt.scaled_height()),
					       (int) floor(_pt.scaled_width()), 3);

		_pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);

		view_generator vg(_pt);

		for (unsigned int i = 0; i < im->height(); i++)
		for (unsigned int j = 0; j < im->width();  j++) {

			focus::result _focus = focus::get(depths, i, j);

			if (!finite(_focus.focal_distance))
				continue;

			/*
			 * Data structures for calculating focal statistics.
			 */
			
			d2::pixel color, weight;
			d2::image_weighted_median *iwm = NULL;

			if (_focus.statistic == 1) {
				iwm = new d2::image_weighted_median(1, 1, 3, _focus.sample_count);
			}

			/*
			 * Iterate over views for this focus region.
			 */

			for (unsigned int v = 0; v < _focus.sample_count; v++) {

				view_generator::view vw = vg.get_view(_focus.aperture, v, _focus.randomization);

				ui::get()->d3_render_status(0, 1, -1, v, i, j, -1);


				/*
				 * Map the focused point to the new view.
				 */

				point p = vw.get_pt().wp_scaled(_pt.pw_scaled(point(i, j, _focus.focal_distance)));

				/*
				 * Determine weight and color for the given point.
				 */

				d2::pixel view_weight, view_color;

				vw.get_color_and_weight(&view_color, &view_weight, p.xy());

				if (!color.finite() || !weight.finite())
					continue;

				if (_focus.statistic == 0) {
					color += view_color;
					weight += view_weight;
				} else if (_focus.statistic == 1) {
					iwm->accumulate(0, 0, v, view_color, view_weight);
				} else
					assert(0);
			}

			if (_focus.statistic == 1) {
				weight = iwm->get_weights()->get_pixel(0, 0);
				color = iwm->get_pixel(0, 0);
				delete iwm;
			}

			if (weight.min_norm() < encounter_threshold) {
				im->set_pixel(i, j, d2::pixel::zero() / d2::pixel::zero());
			} else if (normalize_weights)
				im->set_pixel(i, j, color / weight);
			else
				im->set_pixel(i, j, color);
		}

		delete depths;

		return im;
	}

	/*
	 * Unfiltered function
	 */
	static const d2::image *view_nofilter(pt _pt, int n) {

		if (!focus::is_trivial())
			return view_nofilter_focus(_pt, n);

		assert ((unsigned int) n < d2::image_rw::count() || n < 0);

		if (n >= 0) {
			assert((int) floor(d2::align::of(n).scaled_height())
			     == (int) floor(_pt.scaled_height()));
			assert((int) floor(d2::align::of(n).scaled_width())
			     == (int) floor(_pt.scaled_width()));
		}

		const d2::image *depths = depth(_pt, n);

		d2::image *im = d2::new_image_ale_real((int) floor(_pt.scaled_height()),
					       (int) floor(_pt.scaled_width()), 3);

		_pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);

		/*
		 * Use adaptive subspace data.
		 */

		d2::image *weights = d2::new_image_ale_real((int) floor(_pt.scaled_height()),
						(int) floor(_pt.scaled_width()), 3);

		/*
		 * Iterate through subspaces.
		 */

		space::iterate si(_pt.origin());

		ui::get()->d3_render_status(0, 0, -1, -1, -1, -1, 0);

		view_recurse(0, im, weights, si, _pt);

		for (unsigned int i = 0; i < im->height(); i++)
		for (unsigned int j = 0; j < im->width();  j++) {
			if (weights->get_pixel(i, j).min_norm() < encounter_threshold
			 || (d3px_count > 0 && isnan(depths->get_chan(i, j, 0)))) {
				im->set_pixel(i, j, d2::pixel::zero() / d2::pixel::zero());
				weights->set_pixel(i, j, d2::pixel::zero());
			} else if (normalize_weights)
				im->set_pixel(i, j, (d2::pixel) im->get_pixel(i, j) 
				                  / (d2::pixel) weights->get_pixel(i, j));
		}

		delete weights;

		delete depths;

		return im;
	}

	/*
	 * Filtered function.
	 */
	static const d2::image *view_filter_focus(pt _pt, int n) {

		assert ((unsigned int) n < d2::image_rw::count() || n < 0);

		/*
		 * Get depth image for focus region determination.
		 */

		const d2::image *depths = depth(_pt, n);

		unsigned int height = (unsigned int) floor(_pt.scaled_height());
		unsigned int width = (unsigned int) floor(_pt.scaled_width());

		/*
		 * Prepare input frame data.
		 */

		if (tc_multiplier == 0)
			al->open_all();

		pt *_ptf = new pt[al->count()];
		std::vector<space::node *> *fmv = new std::vector<space::node *>[al->count()];

		for (unsigned int f = 0; f < al->count(); f++) {
			_ptf[f] = al->get(f)->get_t(0);
			fmv[f] = most_visible_unscaled(_ptf[f]);
			std::sort(fmv[f].begin(), fmv[f].end());
		}

		if (tc_multiplier == 0)
			al->close_all();

		/*
		 * Open all files for rendering.
		 */

		d2::image_rw::open_all();

		/*
		 * Prepare data structures for averaging views, as we render
		 * each view separately.  This is spacewise inefficient, but
		 * is easy to implement given the current operation of the
		 * renderers.
		 */

		d2::image_weighted_avg *iwa;

		if (d3::focus::uses_medians()) {
			iwa = new d2::image_weighted_median(height, width, 3, focus::max_samples());
		} else {
			iwa = new d2::image_weighted_simple(height, width, 3, new d2::invariant(NULL));
		}

		_pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);

		/*
		 * Prepare view generator.
		 */

		view_generator vg(_pt);

		/*
		 * Render views separately.  This is spacewise inefficient,
		 * but is easy to implement given the current operation of the
		 * renderers.
		 */

		for (unsigned int v = 0; v < focus::max_samples(); v++) {

			/*
			 * Generate a new 2D renderer for filtering.
			 */

			d2::render::reset();
			d2::render *renderer = d2::render_parse::get(d3chain_type);

			renderer->init_point_renderer(height, width, 3);

			/*
			 * Iterate over output points.
			 */
			
			for (unsigned int i = 0; i < height; i++)
			for (unsigned int j = 0; j < width; j++) {

				focus::result _focus = focus::get(depths, i, j);

				if (v >= _focus.sample_count)
					continue;

				if (!finite(_focus.focal_distance))
					continue;

				view_generator::view vw = vg.get_view(_focus.aperture, v, _focus.randomization);

				pt _pt_new = vw.get_pt();

				point p = _pt_new.wp_scaled(_pt.pw_scaled(point(i, j, _focus.focal_distance)));

				/*
				 * Determine the most-visible subspace.
				 */

				space::node *mv = vw.get_most_visible(p.xy());

				if (mv == NULL)
					continue;

				/*
				 * Get median depth and diff.
				 */

				d2::pixel depth, diff;

				vw.get_median_depth_and_diff(&depth, &diff, p.xy());

				if (!depth.finite() || !diff.finite())
					continue;

				point local_points[3] = { 
					point(p[0],     p[1],     ale_real_to_pos(depth[0])),
					point(p[0] + 1, p[1],     ale_real_to_pos(depth[0] + diff[0])),
					point(p[0],     p[1] + 1, ale_real_to_pos(depth[0] + diff[1]))
				};

				/*
				 * Iterate over files.
				 */

				for (unsigned int f = 0; f < d2::image_rw::count(); f++) {

					ui::get()->d3_render_status(1, 1, f, v, i, j, -1);

					if (!visibility_search(fmv[f], mv))
						continue;

					/*
					 * Determine transformation at (i, j).  First
					 * determine transformation from the output to
					 * the input, then invert this, as we need the
					 * inverse transformation for filtering.
					 */

					d2::point remote_points[3] = {
						_ptf[f].wp_unscaled(_pt_new.pw_scaled(point(local_points[0]))).xy(),
						_ptf[f].wp_unscaled(_pt_new.pw_scaled(point(local_points[1]))).xy(),
						_ptf[f].wp_unscaled(_pt_new.pw_scaled(point(local_points[2]))).xy()
					};

					/*
					 * Forward matrix for the linear component of the 
					 * transformation.
					 */

					d2::point forward_matrix[2] = {
						remote_points[1] - remote_points[0],
						remote_points[2] - remote_points[0]
					};

					/*
					 * Inverse matrix for the linear component of
					 * the transformation.  Calculate using the
					 * determinant D.
					 */

					ale_pos D = forward_matrix[0][0] * forward_matrix[1][1]
						  - forward_matrix[0][1] * forward_matrix[1][0];

					if (D == 0)
						continue;

					d2::point inverse_matrix[2] = {
						d2::point( forward_matrix[1][1] / D, -forward_matrix[1][0] / D),
						d2::point(-forward_matrix[0][1] / D,  forward_matrix[0][0] / D)
					};

					/*
					 * Determine the projective transformation parameters for the
					 * inverse transformation.
					 */
					
					const d2::image *imf = d2::image_rw::get_open(f);

					d2::transformation inv_t = d2::transformation::gpt_identity(imf, 1);

					d2::point local_bounds[4];

					for (int n = 0; n < 4; n++) {
						d2::point remote_bound = d2::point((n == 1 || n == 2) ? imf->height() : 0,
										   (n == 2 || n == 3) ? imf->width()  : 0)
								       - remote_points[0];

						local_bounds[n] = d2::point(i, j)
								+ d2::point(remote_bound[0] * inverse_matrix[0][0]
									  + remote_bound[1] * inverse_matrix[1][0],
									    remote_bound[0] * inverse_matrix[0][1]
									  + remote_bound[1] * inverse_matrix[1][1]);

					}

					if (!local_bounds[0].finite()
					 || !local_bounds[1].finite()
					 || !local_bounds[2].finite()
					 || !local_bounds[3].finite())
						continue;

					inv_t.gpt_set(local_bounds);

					/*
					 * Perform render step for the given frame,
					 * transformation, and point.
					 */

					renderer->point_render(i, j, f, inv_t);
				}
			}

			renderer->finish_point_rendering();

			const d2::image *im = renderer->get_image();
			const d2::image *df = renderer->get_defined();

			for (unsigned int i = 0; i < height; i++)
			for (unsigned int j = 0; j < width; j++) {
				if (((d2::pixel) df->get_pixel(i, j)).finite()
				 && df->get_pixel(i, j)[0] > 0)
					iwa->accumulate(i, j, v, im->get_pixel(i, j), d2::pixel(1, 1, 1));
			}
		}

		/*
		 * Close all files and return the result.
		 */

		d2::image_rw::close_all();

		return iwa;
	}

	static const d2::image *view_filter(pt _pt, int n) {

		if (!focus::is_trivial())
			return view_filter_focus(_pt, n);

		assert ((unsigned int) n < d2::image_rw::count() || n < 0);

		/*
		 * Generate a new 2D renderer for filtering.
		 */

		d2::render::reset();
		d2::render *renderer = d2::render_parse::get(d3chain_type);

		/*
		 * Get depth image in order to estimate normals (and hence
		 * transformations).
		 */

		const d2::image *depths = depth(_pt, n);

		d2::image *median_diffs = depths->fcdiff_median((int) floor(diff_median_radius));
		d2::image *median_depths = depths->medians((int) floor(depth_median_radius));

		unsigned int height = (unsigned int) floor(_pt.scaled_height());
		unsigned int width = (unsigned int) floor(_pt.scaled_width());

		renderer->init_point_renderer(height, width, 3);

		_pt.view_angle(_pt.view_angle() * VIEW_ANGLE_MULTIPLIER);

		std::vector<space::node *> mv = most_visible_scaled(_pt);

		for (unsigned int f = 0; f < d2::image_rw::count(); f++) {

			if (tc_multiplier == 0)
				al->open(f);

			pt _ptf = al->get(f)->get_t(0);

			std::vector<space::node *> fmv = most_visible_unscaled(_ptf);
			std::sort(fmv.begin(), fmv.end());

			for (unsigned int i = 0; i < height; i++)
			for (unsigned int j = 0; j < width; j++) {

				ui::get()->d3_render_status(1, 0, f, -1, i, j, -1);

				/*
				 * Check visibility.
				 */

				int n = i * width + j;

				if (!visibility_search(fmv, mv[n]))
					continue;

				/*
				 * Find depth and diff at this point, check for
				 * undefined values, and generate projections
				 * of the image corners on the estimated normal
				 * surface.
				 */

				d2::pixel depth = median_depths->get_pixel(i, j);
				d2::pixel diff = median_diffs->get_pixel(i, j);
				// d2::pixel diff = d2::pixel(0, 0, 0);

				if (!depth.finite() || !diff.finite())
					continue;

				point local_points[3] = { 
					point(i,     j,     ale_real_to_pos(depth[0])),
				        point(i + 1, j,     ale_real_to_pos(depth[0] + diff[0])),
				        point(i    , j + 1, ale_real_to_pos(depth[0] + diff[1]))
				};

				/*
				 * Determine transformation at (i, j).  First
				 * determine transformation from the output to
				 * the input, then invert this, as we need the
				 * inverse transformation for filtering.
				 */

				d2::point remote_points[3] = {
					_ptf.wp_unscaled(_pt.pw_scaled(point(local_points[0]))).xy(),
					_ptf.wp_unscaled(_pt.pw_scaled(point(local_points[1]))).xy(),
					_ptf.wp_unscaled(_pt.pw_scaled(point(local_points[2]))).xy()
				};

				/*
				 * Forward matrix for the linear component of the 
				 * transformation.
				 */

				d2::point forward_matrix[2] = {
					remote_points[1] - remote_points[0],
					remote_points[2] - remote_points[0]
				};

				/*
				 * Inverse matrix for the linear component of
				 * the transformation.  Calculate using the
				 * determinant D.
				 */

				ale_pos D = forward_matrix[0][0] * forward_matrix[1][1]
					  - forward_matrix[0][1] * forward_matrix[1][0];

				if (D == 0)
					continue;

				d2::point inverse_matrix[2] = {
					d2::point( forward_matrix[1][1] / D, -forward_matrix[1][0] / D),
					d2::point(-forward_matrix[0][1] / D,  forward_matrix[0][0] / D)
				};

				/*
				 * Determine the projective transformation parameters for the
				 * inverse transformation.
				 */
				
				const d2::image *imf = d2::image_rw::open(f);

				d2::transformation inv_t = d2::transformation::gpt_identity(imf, 1);

				d2::point local_bounds[4];

				for (int n = 0; n < 4; n++) {
					d2::point remote_bound = d2::point((n == 1 || n == 2) ? imf->height() : 0,
							                   (n == 2 || n == 3) ? imf->width()  : 0)
							       - remote_points[0];

					local_bounds[n] = local_points[0].xy()
						        + d2::point(remote_bound[0] * inverse_matrix[0][0]
							          + remote_bound[1] * inverse_matrix[1][0],
								    remote_bound[0] * inverse_matrix[0][1]
								  + remote_bound[1] * inverse_matrix[1][1]);
				}

				inv_t.gpt_set(local_bounds);

				d2::image_rw::close(f);

				/*
				 * Perform render step for the given frame,
				 * transformation, and point.
				 */

				d2::image_rw::open(f);
				renderer->point_render(i, j, f, inv_t);
				d2::image_rw::close(f);
			}

			if (tc_multiplier == 0) 
				al->close(f);
		}

		renderer->finish_point_rendering();

		return renderer->get_image();
	}

	/*
	 * Generic function.
	 */
	static const d2::image *view(pt _pt, int n = -1) {

		assert ((unsigned int) n < d2::image_rw::count() || n < 0);

		if (use_filter) {
			return view_filter(_pt, n);
		} else {
			return view_nofilter(_pt, n);
		}
	}
		
	static void tcem(double _tcem) {
		tc_multiplier = _tcem;
	}

	static void oui(unsigned int _oui) {
		ou_iterations = _oui;
	}

	static void pa(unsigned int _pa) {
		pairwise_ambiguity = _pa;
	}

	static void pc(const char *_pc) {
		pairwise_comparisons = _pc;
	}

	static void d3px(int _d3px_count, double *_d3px_parameters) {
		d3px_count = _d3px_count;
		d3px_parameters = _d3px_parameters;
	}

	static void fx(double _fx) {
		falloff_exponent = _fx;
	}

	static void nw() {
		normalize_weights = 1;
	}

	static void no_nw() {
		normalize_weights = 0;
	}

	static void nofilter() {
		use_filter = 0;
	}

	static void filter() {
		use_filter = 1;
	}

	static void set_filter_type(const char *type) {
		d3chain_type = type;
	}

	static void set_subspace_traverse() {
		subspace_traverse = 1;
	}

	static int excluded(point p) {
		for (int n = 0; n < d3px_count; n++) {
			double *region = d3px_parameters + (6 * n);
			if (p[0] >= region[0]
			 && p[0] <= region[1]
			 && p[1] >= region[2]
			 && p[1] <= region[3]
			 && p[2] >= region[4]
			 && p[2] <= region[5])
				return 1;
		}

		return 0;
	}

	/*
	 * This function returns true if a space is completely excluded.
	 */
	static int excluded(const space::traverse &st) {
		for (int n = 0; n < d3px_count; n++) {
			double *region = d3px_parameters + (6 * n);
			if (st.get_min()[0] >= region[0]
			 && st.get_max()[0] <= region[1]
			 && st.get_min()[1] >= region[2]
			 && st.get_max()[1] <= region[3]
			 && st.get_min()[2] >= region[4]
			 && st.get_max()[2] <= region[5])
				return 1;
		}

		return 0;
	}

	static const d2::image *view(unsigned int n) {

		assert (n < d2::image_rw::count());

		pt _pt = align::projective(n);

		return view(_pt, n);
	}

	typedef struct {point iw; point ip, is;} analytic;
	typedef std::multimap<ale_real,analytic> score_map;
	typedef std::pair<ale_real,analytic> score_map_element;

	/*
	 * Make pt list.
	 */
	static std::vector<pt> make_pt_list(const char *d_out[], const char *v_out[],
			std::map<const char *, pt> *d3_depth_pt,
			std::map<const char *, pt> *d3_output_pt) {

		std::vector<pt> result;

		for (unsigned int n = 0; n < d2::image_rw::count(); n++) {
			if (d_out[n] || v_out[n]) {
				result.push_back(align::projective(n));
			}
		}

		for (std::map<const char *, pt>::iterator i = d3_depth_pt->begin(); i != d3_depth_pt->end(); i++) {
			result.push_back(i->second);
		}

		for (std::map<const char *, pt>::iterator i = d3_output_pt->begin(); i != d3_output_pt->end(); i++) {
			result.push_back(i->second);
		}

		return result;
	}

	/*
	 * Get a trilinear coordinate for an anisotropic candidate cell.
	 */
	static ale_pos get_trilinear_coordinate(point min, point max, pt _pt) {

		d2::point local_min, local_max;

		local_min = _pt.wp_unscaled(min).xy();
		local_max = _pt.wp_unscaled(min).xy();

		point cell[2] = {min, max};

		/*
		 * Determine the view-local extrema in 2 dimensions.
		 */

		for (int r = 1; r < 8; r++) {
			point local = _pt.wp_unscaled(point(cell[r>>2][0], cell[(r>>1)%2][1], cell[r%2][2]));
			
			for (int d = 0; d < 2; d++) {
				if (local[d] < local_min[d])
					local_min[d] = local[d];
				if (local[d] > local_max[d])
					local_max[d] = local[d];
				if (isnan(local[d])) 
					return local[d];
			}
		}

		ale_pos diameter = (local_max - local_min).norm();

		return log((double) diameter / sqrt(2)) / log(2);
	}

	/*
	 * Check whether a cell is visible from a given viewpoint.  This function
	 * is guaranteed to return 1 when a cell is visible, but it is not guaranteed
	 * to return 0 when a cell is invisible.
	 */
	static int pt_might_be_visible(const pt &viewpoint, point min, point max) {

		int doc = (rand() % 100000) ? 0 : 1;

		if (doc)
			fprintf(stderr, "checking visibility:\n");

		point cell[2] = {min, max};

		/*
		 * Cycle through all vertices of the cell to check certain
		 * properties.
		 */
		int pos[3] = {0, 0, 0};
		int neg[3] = {0, 0, 0};
		for (int i = 0; i < 2; i++)
		for (int j = 0; j < 2; j++)
		for (int k = 0; k < 2; k++) {
			point p = viewpoint.wp_unscaled(point(cell[i][0], cell[j][1], cell[k][2]));

			if (p[2] < 0 && viewpoint.unscaled_in_bounds(p))
				return 1;

			if (isnan(p[0])
			 || isnan(p[1])
			 || isnan(p[2]))
				return 1;

			if (p[2] > 0)
				for (int d = 0; d < 2; d++)
					p[d] *= -1;

			if (doc)
				fprintf(stderr, "\t[%f %f %f] --> [%f %f %f]\n", 
						(double) cell[i][0], (double) cell[j][1], (double) cell[k][2],
						(double) p[0], (double) p[1], (double) p[2]);

			for (int d = 0; d < 3; d++)
				if (p[d] >= 0)
					pos[d] = 1;

			if (p[0] <= viewpoint.unscaled_height() - 1)
				neg[0] = 1;

			if (p[1] <= viewpoint.unscaled_width() - 1)
				neg[1] = 1;

			if (p[2] <= 0)
				neg[2] = 1;
		}
		
		if (!neg[2])
			return 0;

		if (!pos[0]
		 || !neg[0]
		 || !pos[1]
		 || !neg[1])
			return 0;

		return 1;
	}

	/*
	 * Check whether a cell is output-visible.
	 */
	static int output_might_be_visible(const std::vector<pt> &pt_outputs, point min, point max) {
		for (unsigned int n = 0; n < pt_outputs.size(); n++)
			if (pt_might_be_visible(pt_outputs[n], min, max))
				return 1;
		return 0;
	}

	/*
	 * Check whether a cell is input-visible.
	 */
	static int input_might_be_visible(unsigned int f, point min, point max) {
		return pt_might_be_visible(align::projective(f), min, max);
	}

	/*
	 * Return true if a cell fails an output resolution bound.
	 */
	static int fails_output_resolution_bound(point min, point max, const std::vector<pt> &pt_outputs) {
		for (unsigned int n = 0; n < pt_outputs.size(); n++) {

			point p = pt_outputs[n].centroid(min, max);

			if (!p.defined())
				continue;

			if (get_trilinear_coordinate(min, max, pt_outputs[n]) < output_decimation_preferred)
				return 1;
		}
		
		return 0;
	}

	/*
	 * Check lower-bound resolution constraints 
	 */
	static int exceeds_resolution_lower_bounds(unsigned int f1, unsigned int f2,
			point min, point max, const std::vector<pt> &pt_outputs) {

		pt _pt = al->get(f1)->get_t(0);

		if (get_trilinear_coordinate(min, max, _pt) < input_decimation_lower)
			return 1;

		if (fails_output_resolution_bound(min, max, pt_outputs))
			return 0;

		if (get_trilinear_coordinate(min, max, _pt) < primary_decimation_upper)
			return 1;

		return 0;
	}

	/*
	 * Try the candidate nearest to the specified cell.
	 */
	static void try_nearest_candidate(unsigned int f1, unsigned int f2, candidates *c, point min, point max) {
		point centroid = (max + min) / 2;
		pt _pt[2] = { al->get(f1)->get_t(0), al->get(f2)->get_t(0) };
		point p[2];

		// fprintf(stderr, "[tnc n=%f %f %f x=%f %f %f]\n", min[0], min[1], min[2], max[0], max[1], max[2]);

		/*
		 * Reject clipping plane violations.
		 */

		if (centroid[2] > front_clip
		 || centroid[2] < rear_clip)
			return;

		/*
		 * Calculate projections.
		 */

		for (int n = 0; n < 2; n++) {

			p[n] = _pt[n].wp_unscaled(centroid);

			if (!_pt[n].unscaled_in_bounds(p[n]))
				return;

			// fprintf(stderr, ":");

			if (p[n][2] >= 0)
				return;
		}


		int tc = (int) round(get_trilinear_coordinate(min, max, _pt[0]));
		int stc = (int) round(get_trilinear_coordinate(min, max, _pt[1]));

		while (tc < input_decimation_lower || stc < input_decimation_lower) {
			tc++;
			stc++;
		}

		if (tc > primary_decimation_upper)
			return;

		/*
		 * Calculate score from color match.  Assume for now
		 * that the transformation can be approximated locally
		 * with a translation.
		 */

		ale_pos score = 0;
		ale_pos divisor = 0;
		ale_real l1_multiplier = 0.125;
		lod_image *if1 = al->get(f1);
		lod_image *if2 = al->get(f2);

		if (if1->in_bounds(p[0].xy())
		 && if2->in_bounds(p[1].xy())) {
			divisor += ale_real_to_pos(1 - l1_multiplier);
			score += ale_real_to_pos((1 - l1_multiplier)
			       * (if1->get_tl(p[0].xy(), tc) - if2->get_tl(p[1].xy(), stc)).normsq());
		}

		for (int iii = -1; iii <= 1; iii++)
		for (int jjj = -1; jjj <= 1; jjj++) {
			d2::point t(iii, jjj);

			if (!if1->in_bounds(p[0].xy() + t)
			 || !if2->in_bounds(p[1].xy() + t))
				continue;

			divisor += ale_real_to_pos(l1_multiplier);
			score   += ale_real_to_pos(l1_multiplier
				 * (if1->get_tl(p[0].xy() + t, tc) - if2->get_tl(p[1].xy() + t, tc)).normsq());
				 
		}

		/*
		 * Include third-camera contributions in the score.
		 */

		if (tc_multiplier != 0)
		for (unsigned int n = 0; n < d2::image_rw::count(); n++) {
			if (n == f1 || n == f2)
				continue;

			lod_image *ifn = al->get(n);
			pt _ptn = ifn->get_t(0);
			point pn = _ptn.wp_unscaled(centroid);

			if (!_ptn.unscaled_in_bounds(pn))
				continue;

			if (pn[2] >= 0)
				continue;

			ale_pos ttc = get_trilinear_coordinate(min, max, _ptn);

			divisor += tc_multiplier;
			score   += tc_multiplier
				 * (if1->get_tl(p[0].xy(), tc) - ifn->get_tl(pn.xy(), ttc)).normsq();
		}

		c->add_candidate(p[0], tc, score / divisor);
	}

	/*
	 * Check for cells that are completely clipped.
	 */
	static int completely_clipped(point min, point max) {
		return (min[2] > front_clip
		     || max[2] < rear_clip);
	}

	/*
	 * Update extremum variables for cell points mapped to a particular view.
	 */
	static void update_extrema(point min, point max, pt _pt, int *extreme_dim, ale_pos *extreme_ratio) {

		point local_min, local_max;

		local_min = _pt.wp_unscaled(min);
		local_max = _pt.wp_unscaled(min);

		point cell[2] = {min, max};

		int near_vertex = 0;

		/*
		 * Determine the view-local extrema in all dimensions, and
		 * determine the vertex of closest z coordinate.
		 */

		for (int r = 1; r < 8; r++) {
			point local = _pt.wp_unscaled(point(cell[r>>2][0], cell[(r>>1)%2][1], cell[r%2][2]));
			
			for (int d = 0; d < 3; d++) {
				if (local[d] < local_min[d])
					local_min[d] = local[d];
				if (local[d] > local_max[d])
					local_max[d] = local[d];
			}

			if (local[2] == local_max[2])
				near_vertex = r;
		}

		ale_pos diameter = (local_max.xy() - local_min.xy()).norm();

		/*
		 * Update extrema as necessary for each dimension.
		 */

		for (int d = 0; d < 3; d++) {

			int r = near_vertex;

			int p1[3] = {r>>2, (r>>1)%2, r%2};
			int p2[3] = {r>>2, (r>>1)%2, r%2};

			p2[d] = 1 - p2[d];

			ale_pos local_distance = (_pt.wp_unscaled(point(cell[p1[0]][0], cell[p1[1]][1], cell[p1[2]][2])).xy()
						- _pt.wp_unscaled(point(cell[p2[0]][0], cell[p2[1]][1], cell[p2[2]][2])).xy()).norm();

			if (local_distance / diameter > *extreme_ratio) {
				*extreme_ratio = local_distance / diameter;
				*extreme_dim = d;
			}
		}
	}

	/*
	 * Get the next split dimension.
	 */
	static int get_next_split(int f1, int f2, point min, point max, const std::vector<pt> &pt_outputs) {
		for (int d = 0; d < 3; d++)
			if (isinf(min[d]) || isinf(max[d]))
				return space::traverse::get_next_split(min, max);

		int extreme_dim = 0;
		ale_pos extreme_ratio = 0;

		update_extrema(min, max, al->get(f1)->get_t(0), &extreme_dim, &extreme_ratio);
		update_extrema(min, max, al->get(f2)->get_t(0), &extreme_dim, &extreme_ratio);

		for (unsigned int n = 0; n < pt_outputs.size(); n++) {
			update_extrema(min, max, pt_outputs[n], &extreme_dim, &extreme_ratio);
		}

		return extreme_dim;
	}

	/*
	 * Find candidates for subspace creation.
	 */
	static void find_candidates(unsigned int f1, unsigned int f2, candidates *c, point min, point max,
			const std::vector<pt> &pt_outputs, int depth = 0) {

		int print = 0;

		if (min[0] < 20.0001 && max[0] > 20.0001
		 && min[1] < 20.0001 && max[1] > 20.0001
		 && min[2] < 0.0001 && max[2] > 0.0001)
			print = 1;

		if (print) {
			for (int i = depth; i > 0; i--) {
				fprintf(stderr, "+");
			}
			fprintf(stderr, "[fc n=%f %f %f x=%f %f %f]\n",
					(double) min[0], (double) min[1], (double) min[2], (double) max[0], (double) max[1], (double) max[2]);
		}

		if (completely_clipped(min, max)) {
			if (print)
				fprintf(stderr, "c");
			return;
		}

		if (!input_might_be_visible(f1, min, max)
		 || !input_might_be_visible(f2, min, max)) {
			if (print)
				fprintf(stderr, "v");
			return;
		}

		if (output_clip && !output_might_be_visible(pt_outputs, min, max)) {
			if (print)
				fprintf(stderr, "o");
			return;
		}

		if (exceeds_resolution_lower_bounds(f1, f2, min, max, pt_outputs)) {
			if (!(rand() % 100000))
				fprintf(stderr, "([%f %f %f], [%f %f %f]) at %d\n", 
						(double) min[0], (double) min[1], (double) min[2],
						(double) max[0], (double) max[1], (double) max[2],
						__LINE__);

			if (print)
				fprintf(stderr, "t");

			try_nearest_candidate(f1, f2, c, min, max);
			return;
		}

		point new_cells[2][2];

		if (!space::traverse::get_next_cells(get_next_split(f1, f2, min, max, pt_outputs), min, max, new_cells)) {
			if (print)
				fprintf(stderr, "n");
			return;
		}

		if (print) {
			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",
					(double) new_cells[0][0][0],
					(double) new_cells[0][0][1],
					(double) new_cells[0][0][2],
					(double) new_cells[0][1][0],
					(double) new_cells[0][1][1],
					(double) new_cells[0][1][2],
					(double) new_cells[1][0][0],
					(double) new_cells[1][0][1],
					(double) new_cells[1][0][2],
					(double) new_cells[1][1][0],
					(double) new_cells[1][1][1],
					(double) new_cells[1][1][2]);
		}

		find_candidates(f1, f2, c, new_cells[0][0], new_cells[0][1], pt_outputs, depth + 1);
		find_candidates(f1, f2, c, new_cells[1][0], new_cells[1][1], pt_outputs, depth + 1);
	}

	/*
	 * Generate a map from scores to 3D points for various depths at point (i, j) in f1, at 
	 * lowest resolution.
	 */
	static score_map p2f_score_map(unsigned int f1, unsigned int f2, unsigned int i, unsigned int j) {

		score_map result;

		pt _pt1 = al->get(f1)->get_t(primary_decimation_upper);
		pt _pt2 = al->get(f2)->get_t(primary_decimation_upper);

		const d2::image *if1 = al->get(f1)->get_image(primary_decimation_upper);
		const d2::image *if2 = al->get(f2)->get_image(primary_decimation_upper);
		ale_pos pdu_scale = pow(2, primary_decimation_upper);

		/*
		 * Get the pixel color in the primary frame
		 */

		// d2::pixel color_primary = if1->get_pixel(i, j);

		/*
		 * Map two depths to the secondary frame.
		 */

		point p1 = _pt2.wp_unscaled(_pt1.pw_unscaled(point(i, j, 1000)));
		point p2 = _pt2.wp_unscaled(_pt1.pw_unscaled(point(i, j, -1000)));

//		fprintf(stderr, "%d->%d (%d, %d) point pair: (%d, %d, %d -> %f, %f), (%d, %d, %d -> %f, %f)\n",
//				f1, f2, i, j, i, j, 1000, p1[0], p1[1], i, j, -1000, p2[0], p2[1]);
//		_pt1.debug_output();
//		_pt2.debug_output();


		/*
		 * For cases where the mapped points define a
		 * line and where points on the line fall
		 * within the defined area of the frame,
		 * determine the starting point for inspection.
		 * In other cases, continue to the next pixel.
		 */

		ale_pos diff_i = p2[0] - p1[0];
		ale_pos diff_j = p2[1] - p1[1];
		ale_pos slope = diff_j / diff_i;

		if (isnan(slope)) {
			assert(0);
			fprintf(stderr, "%d->%d (%d, %d) has undefined slope\n",
					f1, f2, i, j);
			return result;
		}

		/*
		 * Make absurdly large/small slopes either infinity, negative infinity, or zero.
		 */

		if (fabs(slope) > if2->width() * 100) {
			double zero = 0;
			double one  = 1;
			double inf  = one / zero;
			slope = inf;
		} else if (slope < 1 / (double) if2->height() / 100
			&& slope > -1/ (double) if2->height() / 100) {
			slope = 0;
		}

		// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

		ale_pos top_intersect = p1[1] - p1[0] * slope;
		ale_pos lef_intersect = p1[0] - p1[1] / slope;
		ale_pos rig_intersect = p1[0] - (p1[1] - if2->width() + 2) / slope;
		ale_pos sp_i, sp_j;

		// fprintf(stderr, "slope == %f\n", slope);


		if (slope == 0) {
			// fprintf(stderr, "case 0\n");
			sp_i = lef_intersect;
			sp_j = 0;
		} else if (finite(slope) && top_intersect >= 0 && top_intersect < if2->width() - 1) {
			// fprintf(stderr, "case 1\n");
			sp_i = 0;
			sp_j = top_intersect;
		} else if (slope > 0 && lef_intersect >= 0 && lef_intersect <= if2->height() - 1) {
			// fprintf(stderr, "case 2\n");
			sp_i = lef_intersect;
			sp_j = 0;
		} else if (slope < 0 && rig_intersect >= 0 && rig_intersect <= if2->height() - 1) {
			// fprintf(stderr, "case 3\n");
			sp_i = rig_intersect;
			sp_j = if2->width() - 2;
		} else {
			// fprintf(stderr, "case 4\n");
			// fprintf(stderr, "%d->%d (%d, %d) does not intersect the defined area\n",
			//		f1, f2, i, j);
			return result;
		}


		// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

		/*
		 * Determine increment values for examining
		 * point, ensuring that incr_i is always
		 * positive.
		 */

		ale_pos incr_i, incr_j;

		if (fabs(diff_i) > fabs(diff_j)) {
			incr_i = 1;
			incr_j = slope;
		} else if (slope > 0) {
			incr_i = 1 / slope;
			incr_j = 1;
		} else {
			incr_i = -1 / slope;
			incr_j = -1;
		}
		
		// fprintf(stderr, "%d->%d (%d, %d) increments are (%f, %f)\n",
		//		f1, f2, i, j, incr_i, incr_j);

		/*
		 * Examine regions near the projected line.
		 */

		for (ale_pos ii = sp_i, jj = sp_j; 
			ii <= if2->height() - 1 && jj <= if2->width() - 1 && ii >= 0 && jj >= 0; 
			ii += incr_i, jj += incr_j) {

			// fprintf(stderr, "%d->%d (%d, %d) checking (%f, %f)\n", 
			//		f1, f2, i, j, ii, jj);

#if 0
			/*
			 * Check for higher, lower, and nearby points.
			 *
			 *	Red   = 2^0
			 *	Green = 2^1
			 *	Blue  = 2^2
			 */

			int higher = 0, lower = 0, nearby = 0;

			for (int iii = 0; iii < 2; iii++)
			for (int jjj = 0; jjj < 2; jjj++) {
				d2::pixel p = if2->get_pixel((int) floor(ii) + iii, (int) floor(jj) + jjj);

				for (int k = 0; k < 3; k++) {
					int bitmask = (int) pow(2, k);

					if (p[k] > color_primary[k])
						higher |= bitmask;
					if (p[k] < color_primary[k])
						lower  |= bitmask;
					if (fabs(p[k] - color_primary[k]) < nearness)
						nearby |= bitmask;
				}
			}

			/*
			 * If this is not a region of interest,
			 * then continue.
			 */


			fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			// if (((higher & lower) | nearby) != 0x7)
			//	continue;
#endif
			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			// fprintf(stderr, "%d->%d (%d, %d) accepted (%f, %f)\n", 
			//		f1, f2, i, j, ii, jj);

			/*
			 * Create an orthonormal basis to
			 * determine line intersection.
			 */

			point bp0 = _pt1.pw_unscaled(point(i, j, 0));
			point bp1 = _pt1.pw_unscaled(point(i, j, 10));
			point bp2 = _pt2.pw_unscaled(point(ii, jj, 0));

			point foo = _pt1.wp_unscaled(bp0);
			// fprintf(stderr, "(%d, %d, 0) transformed to world and back is: (%f, %f, %f)\n",
			//		i, j, foo[0], foo[1], foo[2]);

			foo = _pt1.wp_unscaled(bp1);
			// fprintf(stderr, "(%d, %d, 10)  transformed to world and back is: (%f, %f, %f)\n",
			//		i, j, foo[0], foo[1], foo[2]);

			point b0  = (bp1 - bp0).normalize();
			point b1n = bp2 - bp0;
			point b1  = (b1n - b1n.dproduct(b0) * b0).normalize();
			point b2  = point(0, 0, 0).xproduct(b0, b1).normalize(); // Should already have norm=1
			

			foo = _pt1.wp_unscaled(bp0 + 30 * b0);

			/*
			 * Select a fourth point to define a second line.
			 */

			point p3  = _pt2.pw_unscaled(point(ii, jj, 10));

			/*
			 * Representation in the new basis.
			 */

			d2::point nbp0 = d2::point(bp0.dproduct(b0), bp0.dproduct(b1));
			// d2::point nbp1 = d2::point(bp1.dproduct(b0), bp1.dproduct(b1));
			d2::point nbp2 = d2::point(bp2.dproduct(b0), bp2.dproduct(b1));
			d2::point np3  = d2::point( p3.dproduct(b0),  p3.dproduct(b1));

			/*
			 * Determine intersection of line
			 * (nbp0, nbp1), which is parallel to
			 * b0, with line (nbp2, np3).
			 */

			/*
			 * XXX: a stronger check would be
			 * better here, e.g., involving the
			 * ratio (np3[0] - nbp2[0]) / (np3[1] -
			 * nbp2[1]).  Also, acceptance of these
			 * cases is probably better than
			 * rejection.
			 */


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			if (np3[1] - nbp2[1] == 0)
				continue;


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			d2::point intersection = d2::point(nbp2[0] 
					+ (nbp0[1] - nbp2[1]) * (np3[0] - nbp2[0]) / (np3[1] - nbp2[1]),
					nbp0[1]);

			ale_pos b2_offset = b2.dproduct(bp0);

			/*
			 * Map the intersection back to the world
			 * basis.
			 */

			point iw = intersection[0] * b0 + intersection[1] * b1 + b2_offset * b2;

			/*
			 * Reject intersection points behind a
			 * camera.
			 */

			point icp = _pt1.wc(iw);
			point ics = _pt2.wc(iw);


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			if (icp[2] >= 0 || ics[2] >= 0)
				continue;


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			/*
			 * Reject clipping plane violations.
			 */


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			if (iw[2] > front_clip
			 || iw[2] < rear_clip)
				continue;


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			/*
			 * Score the point.
			 */

			point ip = _pt1.wp_unscaled(iw);

			point is = _pt2.wp_unscaled(iw);

			analytic _a = { iw, ip, is };

			/*
			 * Calculate score from color match.  Assume for now
			 * that the transformation can be approximated locally
			 * with a translation.
			 */

			ale_pos score = 0;
			ale_pos divisor = 0;
			ale_pos l1_multiplier = 0.125;

			if (if1->in_bounds(ip.xy())
			 && if2->in_bounds(is.xy())
			 && !d2::render::is_excluded_f(ip.xy() * pdu_scale, f1)
			 && !d2::render::is_excluded_f(is.xy() * pdu_scale, f2)) {
				divisor += 1 - l1_multiplier;
				score += (1 - l1_multiplier)
				       * (ale_pos) ((if1->get_bl(ip.xy()) - if2->get_bl(is.xy())).normsq());
			}


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			for (int iii = -1; iii <= 1; iii++)
			for (int jjj = -1; jjj <= 1; jjj++) {
				d2::point t(iii, jjj);

				if (!if1->in_bounds(ip.xy() + t)
				 || !if2->in_bounds(is.xy() + t)
				 || d2::render::is_excluded_f(ip.xy() * pdu_scale, f1)
				 || d2::render::is_excluded_f(is.xy() * pdu_scale, f2))
					continue;

				divisor += l1_multiplier;
				score   += l1_multiplier
					 * (ale_pos) ((if1->get_bl(ip.xy() + t) - if2->get_bl(is.xy() + t)).normsq());
					 
			}

			/*
			 * Include third-camera contributions in the score.
			 */

			if (tc_multiplier != 0)
			for (unsigned int f = 0; f < d2::image_rw::count(); f++) {
				if (f == f1 || f == f2)
					continue;

				const d2::image *if3 = al->get(f)->get_image(primary_decimation_upper);
				pt _pt3 = al->get(f)->get_t(primary_decimation_upper);

				point p = _pt3.wp_unscaled(iw);

				if (!if3->in_bounds(p.xy())
				 || !if1->in_bounds(ip.xy())
				 || d2::render::is_excluded_f(p.xy() * pdu_scale, f)
				 || d2::render::is_excluded_f(ip.xy() * pdu_scale, f1))
					continue;

				divisor += tc_multiplier;
				score   += tc_multiplier
					 * (if1->get_bl(ip.xy()) - if3->get_bl(p.xy())).normsq();
			}

			


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			/*
			 * Reject points with undefined score.
			 */


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			if (!finite(score / divisor))
				continue;


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

#if 0
			/*
			 * XXX: reject points not on the z=-27.882252 plane.
			 */


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			if (_a.ip[2] > -27 || _a.ip[2] < -28)
				continue;
#endif


			// fprintf(stderr, "score map (%u, %u) line %u\n", i, j, __LINE__);

			/*
			 * Add the point to the score map.
			 */

//			d2::pixel c_ip = if1->in_bounds(ip.xy()) ? if1->get_bl(ip.xy())
//								 : d2::pixel();
//			d2::pixel c_is = if2->in_bounds(is.xy()) ? if2->get_bl(is.xy())
//								 : d2::pixel();

//			fprintf(stderr, "Candidate subspace: f1=%u f2=%u i=%u j=%u ii=%f jj=%f"
//					"cp=[%f %f %f] cs=[%f %f %f]\n",
//					f1, f2, i, j, ii, jj, c_ip[0], c_ip[1], c_ip[2],
//					c_is[0], c_is[1], c_is[2]);

			result.insert(score_map_element(score / divisor, _a));
		}

//		fprintf(stderr, "Iterating through the score map:\n");
//
//		for (score_map::iterator smi = result.begin(); smi != result.end(); smi++) {
//			fprintf(stderr, "%f ", smi->first);
//		}
//
//		fprintf(stderr, "\n");

		return result;
	}


	/*
	 * Attempt to refine space around a point, to high and low resolutions
	 * resulting in two resolutions in total.
	 */

	static space::traverse refine_space(point iw, ale_pos target_dim, int use_filler) {

		space::traverse st = space::traverse::root();

		if (!st.includes(iw)) {
			assert(0);
			return st;
		}

		int lr_done = !use_filler;

		/*
		 * Loop until all resolutions of interest have been generated.
		 */
		
		for(;;) {

			point p[2] = { st.get_min(), st.get_max() };

			ale_pos dim_max = 0;

			for (int d = 0; d < 3; d++) {
				ale_pos d_value = fabs(p[0][d] - p[1][d]);
				if (d_value > dim_max)
					dim_max = d_value;
			}

			/*
			 * Generate any new desired spatial registers.
			 */

			for (int f = 0; f < 2; f++) {

				/*
				 * Low resolution
				 */

				if (dim_max < 2 * target_dim
				 && lr_done == 0) {
					if (spatial_info_map.find(st.get_node()) == spatial_info_map.end()) {
						spatial_info_map[st.get_node()];
						ui::get()->d3_increment_spaces();
					}
					lr_done = 1;
				}

				/*
				 * High resolution.
				 */

				if (dim_max < target_dim) {
					if (spatial_info_map.find(st.get_node()) == spatial_info_map.end()) {
						spatial_info_map[st.get_node()];
						ui::get()->d3_increment_spaces();
					}
					return st;
				}
			}

			/*
			 * Check precision before analyzing space further.
			 */

			if (st.precision_wall()) {
				fprintf(stderr, "\n\n*** Error: reached subspace precision wall ***\n\n");
				assert(0);
				return st;
			}

			if (st.positive().includes(iw)) {
				st = st.positive();
				total_tsteps++;
			} else if (st.negative().includes(iw)) {
				st = st.negative();
				total_tsteps++;
			} else {
				fprintf(stderr, "failed iw = (%f, %f, %f)\n", 
						(double) iw[0], (double) iw[1], (double) iw[2]);
				assert(0);
			}
		}
	}

	/*
	 * Calculate target dimension
	 */

	static ale_pos calc_target_dim(point iw, pt _pt, const char *d_out[], const char *v_out[], 
			std::map<const char *, pt> *d3_depth_pt, 
			std::map<const char *, pt> *d3_output_pt) {

		ale_pos result = _pt.distance_1d(iw, primary_decimation_upper);

		for (unsigned int n = 0; n < d2::image_rw::count(); n++) {
			if (d_out[n] && align::projective(n).distance_1d(iw, 0) < result)
				result = align::projective(n).distance_1d(iw, 0);
			if (v_out[n] && align::projective(n).distance_1d(iw, 0) < result)
				result = align::projective(n).distance_1d(iw, 0);
		}

		for (std::map<const char *, pt>::iterator i = d3_output_pt->begin(); i != d3_output_pt->end(); i++) {
			if (i->second.distance_1d(iw, 0) < result)
				result = i->second.distance_1d(iw, 0);
		}

		for (std::map<const char *, pt>::iterator i = d3_depth_pt->begin(); i != d3_depth_pt->end(); i++) {
			if (i->second.distance_1d(iw, 0) < result)
				result = i->second.distance_1d(iw, 0);
		}

		assert (result > 0);

		return result;
	}

	/*
	 * Calculate level of detail for a given viewpoint.
	 */

	static int calc_lod(ale_pos depth1, pt _pt, ale_pos target_dim) {
		return (int) round(_pt.trilinear_coordinate(depth1, target_dim * (ale_pos) sqrt(2)));
	}

	/*
	 * Calculate depth range for a given pair of viewpoints.
	 */

	static ale_pos calc_depth_range(point iw, pt _pt1, pt _pt2) {

		point ip = _pt1.wp_unscaled(iw);

		ale_pos reference_change = fabs(ip[2] / 1000);

		point iw1 = _pt1.pw_scaled(ip + point(0, 0, reference_change));
		point iw2 = _pt1.pw_scaled(ip - point(0, 0, reference_change));

		point is = _pt2.wc(iw);
		point is1 = _pt2.wc(iw1);
		point is2 = _pt2.wc(iw2);

		assert(is[2] < 0);

		ale_pos d1 = (is1.xy() - is.xy()).norm();
		ale_pos d2 = (is2.xy() - is.xy()).norm();

		// assert (reference_change > 0);
		// assert (d1 > 0 || d2 > 0);

		if (is1[2] < 0 && is2[2] < 0) {

			if (d1 > d2)
				return reference_change / d1;
			else
				return reference_change / d2;
		}

		if (is1[2] < 0)
			return reference_change / d1;

		if (is2[2] < 0)
			return reference_change / d2;

		return 0;
	}

	/*
	 * Calculate a refined point for a given set of parameters.
	 */

	static point get_refined_point(pt _pt1, pt _pt2, int i, int j, 
			int f1, int f2, int lod1, int lod2, ale_pos depth,
			ale_pos depth_range) {

		d2::pixel comparison_color = al->get(f1)->get_image(lod1)->get_pixel(i, j);

		ale_pos best = -1;
		ale_pos best_depth = depth;

		assert (depth_range > 0);

		if (fabs(depth_range) < fabs(depth / 10000))
			return _pt1.pw_unscaled(point(i, j, depth));

		for (ale_pos d = depth - depth_range; d < depth + depth_range; d += depth_range / 10) {

			if (!(d < 0))
				continue;
			
			point iw = _pt1.pw_unscaled(point(i, j, d));
			point is = _pt2.wp_unscaled(iw);

			if (!(is[2] < 0))
				continue;

			if (!al->get(f2)->get_image(lod2)->in_bounds(is.xy()))
				continue;

			ale_pos error = (comparison_color - al->get(f2)->get_image(lod2)->get_bl(is.xy())).norm();

			if (error < best || best == -1) {
				best = error;
				best_depth = d;
			}
		}

		return _pt1.pw_unscaled(point(i, j, best_depth));
	}

	/*
	 * Analyze space in a manner dependent on the score map.
	 */

	static void analyze_space_from_map(const char *d_out[], const char *v_out[],
				       std::map<const char *, pt> *d3_depth_pt,
				       std::map<const char *, pt> *d3_output_pt,
				       unsigned int f1, unsigned int f2, 
				       unsigned int i, unsigned int j, score_map _sm, int use_filler) {

		int accumulated_ambiguity = 0;
		int max_acc_amb = pairwise_ambiguity;

		pt _pt1 = al->get(f1)->get_t(0);
		pt _pt2 = al->get(f2)->get_t(0);

		if (_pt1.scale_2d() != 1)
			use_filler = 1;

		for(score_map::iterator smi = _sm.begin(); smi != _sm.end(); smi++) {

			point iw = smi->second.iw;

			if (accumulated_ambiguity++ >= max_acc_amb)
				break;

			total_ambiguity++;

			ale_pos depth1 = _pt1.wc(iw)[2];
			ale_pos depth2 = _pt2.wc(iw)[2];

			ale_pos target_dim = calc_target_dim(iw, _pt1, d_out, v_out, d3_depth_pt, d3_output_pt);

			assert(target_dim > 0);

			int lod1 = calc_lod(depth1, _pt1, target_dim);
			int lod2 = calc_lod(depth2, _pt2, target_dim);

			while (lod1 < input_decimation_lower
			    || lod2 < input_decimation_lower) {
				target_dim *= 2;
				lod1 = calc_lod(depth1, _pt1, target_dim);
				lod2 = calc_lod(depth2, _pt2, target_dim);
			}


			if (lod1 >= (int) al->get(f1)->count()
			 || lod2 >= (int) al->get(f2)->count())
				continue;

			int multiplier = (unsigned int) floor(pow(2, primary_decimation_upper - lod1));

			ale_pos depth_range = calc_depth_range(iw, _pt1, _pt2);

			assert (depth_range > 0);

			pt _pt1_lod = al->get(f1)->get_t(lod1);
			pt _pt2_lod = al->get(f2)->get_t(lod2);

			int im = i * multiplier;
			int jm = j * multiplier;

			for (int ii = 0; ii < multiplier; ii += 1)
			for (int jj = 0; jj < multiplier; jj += 1) {

				point refined_point = get_refined_point(_pt1_lod, _pt2_lod, im + ii, jm + jj, 
						f1, f2, lod1, lod2, depth1, depth_range);

				/*
				 * Re-evaluate target dimension.
				 */

				ale_pos target_dim_ = 
					calc_target_dim(refined_point, _pt1, d_out, v_out, d3_depth_pt, d3_output_pt);

				ale_pos depth1_ = _pt1.wc(refined_point)[2];
				ale_pos depth2_ = _pt2.wc(refined_point)[2];

				int lod1_ = calc_lod(depth1_, _pt1, target_dim_);
				int lod2_ = calc_lod(depth2_, _pt2, target_dim_);

				while (lod1_ < input_decimation_lower
				    || lod2_ < input_decimation_lower) {
					target_dim_ *= 2;
					lod1_ = calc_lod(depth1_, _pt1, target_dim_);
					lod2_ = calc_lod(depth2_, _pt2, target_dim_);
				}

				/*
				 * Attempt to refine space around the intersection point.
				 */

				space::traverse st = 
					refine_space(refined_point, target_dim_, use_filler || _pt1.scale_2d() != 1);

// 				if (!resolution_ok(al->get(f1)->get_t(0), al->get(f1)->get_t(0).trilinear_coordinate(st))) {
// 					pt transformation = al->get(f1)->get_t(0);
// 					ale_pos tc = al->get(f1)->get_t(0).trilinear_coordinate(st);
// 
// 					fprintf(stderr, "Resolution not ok.\n");  
// 					fprintf(stderr, "pow(2, tc)=%f\n", pow(2, tc));
// 					fprintf(stderr, "transformation.unscaled_height()=%f\n", 
// 							transformation.unscaled_height());
// 					fprintf(stderr, "transformation.unscaled_width()=%f\n", 
// 							transformation.unscaled_width());
// 					fprintf(stderr, "tc=%f", tc);
// 					fprintf(stderr, "input_decimation_lower - 1.5 = %f\n", 
// 							input_decimation_lower - 1.5);
// 
// 				}
// 
// 				assert(resolution_ok(al->get(f1)->get_t(0), al->get(f1)->get_t(0).trilinear_coordinate(st)));
// 				assert(resolution_ok(al->get(f2)->get_t(0), al->get(f2)->get_t(0).trilinear_coordinate(st)));
			}

		}
	}


	/*
	 * Initialize space and identify regions of interest for the adaptive
	 * subspace model.
	 */
	static void make_space(const char *d_out[], const char *v_out[],
			std::map<const char *, pt> *d3_depth_pt,
			std::map<const char *, pt> *d3_output_pt) {

		ui::get()->d3_total_spaces(0);

		/*
		 * Variable indicating whether low-resolution filler space
		 * is desired to avoid aliased gaps in surfaces.
		 */

		int use_filler = d3_depth_pt->size() != 0
				   || d3_output_pt->size() != 0
				   || output_decimation_preferred > 0
				   || input_decimation_lower > 0
				   || !focus::is_trivial()
				   || !strcmp(pairwise_comparisons, "all");

		std::vector<pt> pt_outputs = make_pt_list(d_out, v_out, d3_depth_pt, d3_output_pt);

		/*
		 * Initialize root space.
		 */

		space::init_root();

		/*
		 * Special handling for experimental option 'subspace_traverse'.
		 */

		if (subspace_traverse) {
			/*
			 * Subdivide space to resolve intensity matches between pairs
			 * of frames.
			 */

			for (unsigned int f1 = 0; f1 < d2::image_rw::count(); f1++) {

				if (d3_depth_pt->size() == 0
				 && d3_output_pt->size() == 0
				 && d_out[f1] == NULL
				 && v_out[f1] == NULL)
					continue;

				if (tc_multiplier == 0)
					al->open(f1);

				for (unsigned int f2 = 0; f2 < d2::image_rw::count(); f2++) {

					if (f1 == f2)
						continue;

					if (tc_multiplier == 0)
						al->open(f2);

					candidates *c = new candidates(f1);

					find_candidates(f1, f2, c, point::neginf(), point::posinf(), pt_outputs);



					c->generate_subspaces();

					if (tc_multiplier == 0)
						al->close(f2);
				}

				if (tc_multiplier == 0)
					al->close(f1);
			}

			return;
		}

		/*
		 * Subdivide space to resolve intensity matches between pairs
		 * of frames.
		 */

		for (unsigned int f1 = 0; f1 < d2::image_rw::count(); f1++)
		for (unsigned int f2 = 0; f2 < d2::image_rw::count(); f2++) {
			if (f1 == f2)
				continue;

			if (!d_out[f1] && !v_out[f1] && !d3_depth_pt->size()
			 && !d3_output_pt->size() && strcmp(pairwise_comparisons, "all"))
				continue;

			if (tc_multiplier == 0) {
				al->open(f1);
				al->open(f2);
			}

			/*
			 * Iterate over all points in the primary frame.
			 */

			ale_pos pdu_scale = pow(2, primary_decimation_upper);

			for (unsigned int i = 0; i < al->get(f1)->get_image(primary_decimation_upper)->height(); i++)
			for (unsigned int j = 0; j < al->get(f1)->get_image(primary_decimation_upper)->width();  j++) {

				if (d2::render::is_excluded_f(d2::point(i, j) * pdu_scale, f1))
					continue;

				ui::get()->d3_subdivision_status(f1, f2, i, j);

				total_pixels++;

				/*
				 * Generate a map from scores to 3D points for
				 * various depths in f1.
				 */

				score_map _sm = p2f_score_map(f1, f2, i, j);

				/*
				 * Analyze space in a manner dependent on the score map.
				 */

				analyze_space_from_map(d_out, v_out, d3_depth_pt, d3_output_pt, 
						f1, f2, i, j, _sm, use_filler);

			}

			/*
			 * This ordering may encourage image f1 to be cached.
			 */

			if (tc_multiplier == 0) {
				al->close(f2);
				al->close(f1);
			}
		}
	}


	/*
	 * Update spatial information structures.
	 *
	 * XXX: the name of this function is horribly misleading.  There isn't
	 * even a 'search depth' any longer, since there is no longer any
	 * bounded DFS occurring.
	 */
	static void reduce_cost_to_search_depth(d2::exposure *exp_out, int inc_bit) {

		/*
		 * Subspace model
		 */

		ui::get()->set_steps(ou_iterations);

		for (unsigned int i = 0; i < ou_iterations; i++) {
			ui::get()->set_steps_completed(i);
			spatial_info_update();
		}

	}

#if 0
	/*
	 * Describe a scene to a renderer
	 */
	static void describe(render *r) {
	}
#endif
};

#endif