// Copyright 2002, 2003, 2004 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
*/

/*
 * image.h: Abstract base class for the internal representations of images used
 * by ALE.
 */

#ifndef __image_h__
#define __image_h__

#include "point.h"
#include "pixel.h"
#include "exposure/exposure.h"

#define IMAGE_BAYER_NONE 0

/*
 * This constant indicates that some other default value should be filled in.
 */

#define IMAGE_BAYER_DEFAULT 0x8

/*
 * Do not change these values without inspecting
 * image_bayer_ale_real::r_*_offset().
 */
#define IMAGE_BAYER_RGBG 0x4 /* binary 100 */
#define IMAGE_BAYER_GBGR 0x5 /* binary 101 */
#define IMAGE_BAYER_GRGB 0x6 /* binary 110 */
#define IMAGE_BAYER_BGRG 0x7 /* binary 111 */

class image : protected exposure::listener {
protected:
	static double resident;
	unsigned int _dimx, _dimy, _depth;
	point _offset;
	const char *name;
	mutable exposure *_exp;
	unsigned int bayer;
private:
	/*
	 * Memoized function variables.  We may want to change these even when
	 * *this is constant.
	 */
	mutable int _apm_memo;
	mutable ale_real _apm;
	mutable int _accm_memo;
	mutable pixel _accm;
	mutable int _acm_memo;
	mutable pixel _acm;

	void avg_channel_clamped_magnitude_memo() const {
		unsigned int i, j, k;
		pixel_accum accumulator;
		pixel_accum divisor;

		if (_accm_memo)
			return;

		_accm_memo = 1;

		accumulator = pixel_accum(0, 0, 0);

		for (i = 0; i < _dimy;  i++)
		for (j = 0; j < _dimx;  j++) {
			pixel value = get_pixel(i, j);

			for (k = 0; k < _depth; k++)
				if (finite(value[k])) {
					if (value[k] > 1)
						value[k] = 1;
					if (value[k] < 0)
						value[k] = 0;
					accumulator[k] += value[k];
					divisor[k] += 1;
				}
		}

		accumulator /= divisor;

		_accm = accumulator;
	}

	void avg_channel_magnitude_memo() const {
		unsigned int i, j, k;
		pixel_accum accumulator;
		pixel_accum divisor;

		if (_acm_memo)
			return;

		_acm_memo = 1;

		accumulator = pixel_accum(0, 0, 0);

		for (i = 0; i < _dimy;  i++)
		for (j = 0; j < _dimx;  j++) {
			pixel value = get_pixel(i, j);

			for (k = 0; k < _depth; k++) 
				if (finite(value[k])) {
					accumulator[k] += value[k];
					divisor[k] += 1;
				}
		}

		accumulator /= divisor;

		_acm = accumulator;
	}

protected:
	void image_updated() {
		_apm_memo = 0;
		_acm_memo = 0;
		_accm_memo = 0;
	}

public:
	static void set_resident(double r) {
		resident = r;
	}

	static double get_resident() {
		return resident;
	}

	image (unsigned int dimy, unsigned int dimx, unsigned int depth, 
			const char *name = "anonymous", exposure *_exp = NULL,
			unsigned int bayer = IMAGE_BAYER_NONE) {

		assert (depth == 3);
		_depth = 3;

		_dimx = dimx;
		_dimy = dimy;
		_offset = point(0, 0);
		_apm_memo = 0;
		_acm_memo = 0;
		_accm_memo = 0;
		this->name = name;
		this->_exp = _exp;
		this->bayer = bayer;

		if (_exp != NULL)
			_exp->add_listener(this, name);
	}

	unsigned int get_bayer() const {
		return bayer;
	}

	virtual char get_channels(int i, int j) const {
		return 0x7;
	}

	virtual unsigned int bayer_color(unsigned int i, unsigned int j) const {
		assert(0);
	}

	double storage_size() const {
		if (bayer != IMAGE_BAYER_NONE)
			return _dimx * _dimy * sizeof(ale_real);

		return 3 * _dimx * _dimy * sizeof(ale_real);
	}

	exposure &exp() const {
		return *_exp;
	}

	point offset() const {
		return _offset;
	}

	void set_offset(int i, int j) {
		_offset[0] = i;
		_offset[1] = j;
	}

	void set_offset(point p) {
		_offset = p;
	}

	unsigned int width() const {
		return _dimx;
	}

	unsigned int height() const {
		return _dimy;
	}

	unsigned int depth() const {
		return _depth;
	}

	virtual void set_pixel(unsigned int y, unsigned int x, spixel p) = 0;
	
	virtual spixel get_pixel(unsigned int y, unsigned int x) const = 0;

	virtual spixel get_raw_pixel(unsigned int y, unsigned int x) const {
		return ((const image *)this)->get_pixel(y, x);
	}

	virtual void add_pixel(unsigned int y, unsigned int x, pixel p) {
		assert(0);
	}

	virtual void mul_pixel(unsigned int y, unsigned int x, pixel p) {
		assert(0);
	}

	virtual void div_pixel(unsigned int y, unsigned int x, pixel p) {
		assert(0);
	}

	virtual void add_chan(unsigned int y, unsigned int x, unsigned int k, ale_real c) {
		assert(0);
	}

	virtual void div_chan(unsigned int y, unsigned int x, unsigned int k, ale_real c) {
		assert(0);
	}

	virtual void set_chan(unsigned int y, unsigned int x, unsigned int k, ale_sreal c) = 0;

	virtual ale_sreal get_chan(unsigned int y, unsigned int x, unsigned int k) const = 0;

	ale_real maxval() const {
		ale_real result = get_pixel(0, 0)[0];

		for (unsigned int i = 0; i < _dimy; i++)
		for (unsigned int j = 0; j < _dimx; j++) {
			pixel p = get_pixel(i, j);

			for (unsigned int k = 0; k < _depth; k++)
				if (p[k] > result || !finite(result))
					result = p[k];
		}

		return result;
	}

	ale_real minval() const {
		ale_real result = get_pixel(0, 0)[0];

		for (unsigned int i = 0; i < _dimy; i++)
		for (unsigned int j = 0; j < _dimx; j++) {
			pixel p = get_pixel(i, j);

			for (unsigned int k = 0; k < _depth; k++)
				if (p[k] < result || !finite(result))
					result = p[k];
		}

		return result;
	}

	/*
	 * Get the maximum difference among adjacent pixels.
	 */

	pixel get_max_diff(unsigned int i, unsigned int j) const {
		assert(i <= _dimy - 1);
		assert(j <= _dimx - 1);

		pixel max = get_pixel(i, j), min = get_pixel(i, j);

		for (int ii = -1; ii <= 1; ii++)
		for (int jj = -1; jj <= 1; jj++) {
			int iii = i + ii;
			int jjj = j + jj;

			if (iii < 0)
				continue;
			if (jjj < 0)
				continue;
			if ((unsigned int) iii > _dimy - 1)
				continue;
			if ((unsigned int) jjj > _dimx - 1)
				continue;

			pixel p = get_pixel(iii, jjj);

			for (int d = 0; d < 3; d++) {
				if (p[d] > max[d])
					max[d] = p[d];
				if (p[d] < min[d])
					min[d] = p[d];
			}
		}

		return max - min;
	}

	pixel get_max_diff(point x) const {
		assert (x[0] >= 0);
		assert (x[1] >= 0);
		assert (x[0] <= _dimy - 1);
		assert (x[1] <= _dimx - 1);

		unsigned int i = (unsigned int) round(x[0]);
		unsigned int j = (unsigned int) round(x[1]);

		return get_max_diff(i, j);
	}

	int in_bounds(point x) const {
		if (x[0] < 0
		 || x[1] < 0
		 || x[0] > height() - 1
		 || x[1] > width() - 1)
			return 0;
		if (!x.defined())
			return 0;
		return 1;
	}

	/*
	 * Get a color value at a given position using bilinear interpolation between the
	 * four nearest pixels.
	 */
	pixel get_bl(point x, int defined = 0) const {
//		fprintf(stderr, "get_bl x=%f %f\n", (double) x[0], (double) x[1]);

		pixel result;

		assert (x[0] >= 0);
		assert (x[1] >= 0);
		assert (x[0] <= _dimy - 1);
		assert (x[1] <= _dimx - 1);

		int lx = (int) floor(x[1]);
		int hx = (int) floor(x[1]) + 1;
		int ly = (int) floor(x[0]);
		int hy = (int) floor(x[0]) + 1;

//		fprintf(stderr, "get_bl l=%d %d h=%d %d\n", ly, lx, hy, hx);

		pixel neighbor[4];
		ale_real factor[4];

		neighbor[0] = get_pixel(ly, lx);
		neighbor[1] = get_pixel(hy % _dimy, lx);
		neighbor[2] = get_pixel(hy % _dimy, hx % _dimx);
		neighbor[3] = get_pixel(ly, hx % _dimx);

//		for (int d = 0; d < 4; d++)
//			fprintf(stderr, "neighbor_%d=%f %f %f\n", d,
//					(double) neighbor[d][0],
//					(double) neighbor[d][1],
//					(double) neighbor[d][2]);

		factor[0] = (ale_real) (hx - x[1]) * (ale_real) (hy - x[0]);
		factor[1] = (ale_real) (hx - x[1]) * (ale_real) (x[0] - ly);
		factor[2] = (ale_real) (x[1] - lx) * (ale_real) (x[0] - ly);
		factor[3] = (ale_real) (x[1] - lx) * (ale_real) (hy - x[0]);

//		for (int d = 0; d < 4; d++)
//			fprintf(stderr, "factor_%d=%f\n", d,
//					(double) factor[d]);

		/*
		 * Use bilinear and/or geometric interpolation
		 */

		if (defined == 0) {
			result = pixel(0, 0, 0);

			for (int n = 0; n < 4; n++)
				result += factor[n] * neighbor[n];
		} else {
#if 0
			/*
			 * Calculating the geometric mean may be expensive on
			 * some platforms (e.g., those without floating-point
			 * support.
			 */

			result = pixel(1, 1, 1);

			for (int n = 0; n < 4; n++)
				result *= ppow(neighbor[n], factor[n]);
#else
			/*
			 * Taking the minimum value may be cheaper than
			 * calculating a geometric mean.
			 */

			result = neighbor[0];

			for (int n = 1; n < 4; n++) 
			for (int k = 0; k < 3; k++) 
				if (neighbor[n][k] < result[k])
					result[k] = neighbor[n][k];
#endif
		}

//		fprintf(stderr, "result=%f %f %f\n", 
//				(double) result[0],
//				(double) result[1],
//				(double) result[2]);

		return result;
	}

	pixel get_scaled_bl(point x, ale_pos f, int defined = 0) const {
		point scaled(
				x[0]/f <= height() - 1
			               ? (x[0]/f)
				       : (ale_pos) (height() - 1),
				x[1]/f <= width() - 1
				       ? (x[1]/f)
				       : (ale_pos) (width() - 1));

		return get_bl(scaled, defined);
	}

		
	/*
	 * Make a new image suitable for receiving scaled values.
	 */
	virtual image *scale_generator(int height, int width, int depth, const char *name) const = 0;

	/*
	 * Generate an image of medians within a given radius
	 */

	image *medians(int radius) const {

		assert (radius >= 0);
		
		image *is = scale_generator(height(), width(), depth(), "median");
		assert(is);

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

			std::vector<ale_real> p[3];

			for (int ii = -radius; ii <= radius; ii++)
			for (int jj = -radius; jj <= radius; jj++) {
				int iii = i + ii;
				int jjj = j + jj;

				if (in_bounds(point(iii, jjj)))
					for (int k = 0; k < 3; k++)
						if (finite(get_pixel(iii, jjj)[k]))
							p[k].push_back(get_pixel(iii, jjj)[k]);
			}

			is->set_pixel(i, j, d2::pixel::undefined());

			for (int k = 0; k < 3; k++) {
				std::sort(p[k].begin(), p[k].end());

				unsigned int pkc = p[k].size();

				if (pkc == 0)
					continue;

				if (pkc % 2 == 0) 
					is->set_chan(i, j, k, 
						(p[k][pkc / 2] + p[k][pkc / 2 - 1]) / 2);
				else
					is->set_chan(i, j, k, 
						p[k][pkc / 2]);
			}
		}

		return is;
	}

	/*
	 * Generate an image of differences of the first channel.  The first
	 * coordinate differences are stored in the first channel, second in the
	 * second channel.
	 */

	image *fcdiffs() const {
		image *is = scale_generator(height(), width(), depth(), "diff");

		assert(is);

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

			if (i + 1 < height()
			 && i > 0
			 && !finite(get_chan(i, j, 0))) {

				is->set_chan(i, j, 0, (get_chan(i + 1, j, 0) 
				                     - get_chan(i - 1, j, 0)) / 2);

			} else if (i + 1 < height()
			        && i > 0
			        && finite(get_chan(i + 1, j, 0))
			        && finite(get_chan(i - 1, j, 0))) {

				is->set_chan(i, j, 0, ((get_chan(i, j, 0) - get_chan(i - 1, j, 0))
						     + (get_chan(i + 1, j, 0) - get_chan(i, j, 0))) / 2);

			} else if (i + 1 < height()
				&& finite(get_chan(i + 1, j, 0))) {

				is->set_chan(i, j, 0, get_chan(i + 1, j, 0) - get_chan(i, j, 0));

			} else if (i > 0
				&& finite(get_chan(i - 1, j, 0))) {
				 
				is->set_chan(i, j, 0, get_chan(i, j, 0) - get_chan(i - 1, j, 0));

			} else {
				is->set_chan(i, j, 0, 0);
			}

			if (j + 1 < width()
			 && j > 0
			 && !finite(get_chan(i, j, 0))) {

				is->set_chan(i, j, 1, (get_chan(i, j + 1, 0) - get_chan(i, j - 1, 0)) / 2);

			} else if (j + 1 < width()
			        && j > 0
			        && finite(get_chan(i, j + 1, 0))
			        && finite(get_chan(i, j - 1, 0))) {

				is->set_chan(i, j, 1, ((get_chan(i, j, 0) - get_chan(i, j - 1, 0))
						     + (get_chan(i, j + 1, 0) - get_chan(i, j, 0))) / 2);

			} else if (j + 1 < width() && finite(get_chan(i, j + 1, 0))) {

				is->set_chan(i, j, 1, get_chan(i, j + 1, 0) - get_chan(i, j, 0));

			} else if (j > 0 && finite(get_chan(i, j - 1, 0))) {
				 
				is->set_chan(i, j, 1, get_chan(i, j, 0) - get_chan(i, j - 1, 0));

			} else {
				is->set_chan(i, j, 1, 0);
			}
		}

		return is;
	}

	/*
	 * Generate an image of median (within a given radius) difference of the
	 * first channel.
	 */

	image *fcdiff_median(int radius) const {
		image *diff = fcdiffs();

		assert(diff);

		image *median = diff->medians(radius);

		assert(median);

		delete diff;

		return median;
	}

	/*
	 * Scale by half.  We use the following filter:
	 *
	 * 	1/16	1/8	1/16
	 * 	1/8	1/4	1/8
	 * 	1/16	1/8	1/16
	 *
	 * At the edges, these values are normalized so that the sum of the
	 * weights of contributing pixels is 1.
	 */
	class scale_by_half_threaded : public thread::decompose_domain {
		image *is;
		const image *iu;
	protected:
		void subdomain_algorithm(unsigned int thread,
				int i_min, int i_max, int j_min, int j_max) {

			ale_real _0625 = (ale_real) 0.0625;
			ale_real _125  = (ale_real) 0.125;
			ale_real _25   = (ale_real) 0.25;
			ale_real _0    = (ale_real) 0;

			unsigned int ui_min = (unsigned int) i_min;
			unsigned int ui_max = (unsigned int) i_max;
			unsigned int uj_min = (unsigned int) j_min;
			unsigned int uj_max = (unsigned int) j_max;

			for (unsigned int i = ui_min; i < ui_max; i++)
			for (unsigned int j = uj_min; j < uj_max; j++) {
				is->set_pixel(i, j, 

				      ( ( ((i > 0 && j > 0) 
					    ? iu->get_pixel(2 * i - 1, 2 * j - 1) * _0625
					    : pixel(0, 0, 0))
					+ ((i > 0)
					    ? iu->get_pixel(2 * i - 1, 2 * j) * _125
					    : pixel(0, 0, 0))
					+ ((i > 0 && j < is->width() - 1)
					    ? iu->get_pixel(2 * i - 1, 2 * j + 1) * _0625
					    : pixel(0, 0, 0))
					+ ((j > 0)
					    ? iu->get_pixel(2 * i, 2 * j - 1) * _125
					    : pixel(0, 0, 0))
					+ iu->get_pixel(2 * i, 2 * j) * _25
					+ ((j < is->width() - 1)
					    ? iu->get_pixel(2 * i, 2 * j + 1) * _125
					    : pixel(0, 0, 0))
					+ ((i < is->height() - 1 && j > 0)
					    ? iu->get_pixel(2 * i + 1, 2 * j - 1) * _0625
					    : pixel(0, 0, 0))
					+ ((i < is->height() - 1)
					    ? iu->get_pixel(2 * i + 1, 2 * j) * _125
					    : pixel(0, 0, 0))
					+ ((i < is->height() && j < is->width() - 1)
					    ? iu->get_pixel(2 * i + 1, 2 * j + 1) * _0625 
					    : pixel(0, 0, 0)))

				     /

					( ((i > 0 && j > 0) 
					    ? _0625
					    : _0)
					+ ((i > 0)
					    ? _125
					    : _0)
					+ ((i > 0 && j < is->width() - 1)
					    ? _0625
					    : _0)
					+ ((j > 0)
					    ? _125
					    : _0)
					+ _25
					+ ((j < is->width() - 1)
					    ? _125
					    : _0)
					+ ((i < is->height() - 1 && j > 0)
					    ? _0625
					    : _0)
					+ ((i < is->height() - 1)
					    ? _125
					    : _0)
					+ ((i < is->height() && j < is->width() - 1)
					    ? _0625
					    : _0) ) ) );
			}
		}

	public:
		scale_by_half_threaded(image *_is, const image *_iu) 
					: decompose_domain(0, _is->height(),
					                   0, _is->width()) {
			is = _is;	
			iu = _iu;
		}
	};

	image *scale_by_half(const char *name) const {
		ale_pos f = 0.5;

		image *is = scale_generator(
			(int) floor(height() * (double) f), 
			(int) floor(width()  * (double) f), depth(), name);

		assert(is);

		scale_by_half_threaded sbht(is, this);
		sbht.run();

		is->_offset = point(_offset[0] * f, _offset[1] * f);

		return is;
	}

	/*
	 * Scale by half.  This function uses externally-provided weights,
	 * multiplied by the following filter:
	 *
	 * 	1/16	1/8	1/16
	 * 	1/8	1/4	1/8
	 * 	1/16	1/8	1/16
	 *
	 * Values are normalized so that the sum of the weights of contributing
	 * pixels is 1.
	 */
	image *scale_by_half(const image *weights, const char *name) const {

		if (weights == NULL)
			return scale_by_half(name);

		ale_pos f = 0.5;

		image *is = scale_generator(
			(int) floor(height() * (double) f), 
			(int) floor(width()  * (double) f), depth(), name);

		assert(is);

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

			pixel value = pixel

			      ( ( ((i > 0 && j > 0) 
				    ? (pixel) get_pixel(2 * i - 1, 2 * j - 1) 
				    * (pixel) weights->get_pixel(2 * i - 1, 2 * j - 1) 
				    * (ale_real) 0.0625
				    : pixel(0, 0, 0))
				+ ((i > 0)
				    ? (pixel) get_pixel(2 * i - 1, 2 * j) 
				    * (pixel) weights->get_pixel(2 * i - 1, 2 * j)
				    * 0.125
				    : pixel(0, 0, 0))
				+ ((i > 0 && j < is->width() - 1)
				    ? (pixel) get_pixel(2 * i - 1, 2 * j + 1) 
				    * (pixel) weights->get_pixel(2 * i - 1, 2 * j + 1)
				    * 0.0625
				    : pixel(0, 0, 0))
				+ ((j > 0)
				    ? (pixel) get_pixel(2 * i, 2 * j - 1) 
				    * (pixel) weights->get_pixel(2 * i, 2 * j - 1)
				    * 0.125
				    : pixel(0, 0, 0))
				+ get_pixel(2 * i, 2 * j) 
				    * (pixel) weights->get_pixel(2 * i, 2 * j)
				    * 0.25
				+ ((j < is->width() - 1)
				    ? (pixel) get_pixel(2 * i, 2 * j + 1) 
				    * (pixel) weights->get_pixel(2 * i, 2 * j + 1)
				    * 0.125
				    : pixel(0, 0, 0))
				+ ((i < is->height() - 1 && j > 0)
				    ? (pixel) get_pixel(2 * i + 1, 2 * j - 1) 
				    * (pixel) weights->get_pixel(2 * i + 1, 2 * j - 1)
				    * 0.0625
				    : pixel(0, 0, 0))
				+ ((i < is->height() - 1)
				    ? (pixel) get_pixel(2 * i + 1, 2 * j) 
				    * (pixel) weights->get_pixel(2 * i + 1, 2 * j)
				    * 0.125
				    : pixel(0, 0, 0))
				+ ((i < is->height() && j < is->width() - 1)
				    ? (pixel) get_pixel(2 * i + 1, 2 * j + 1) 
				    * (pixel) weights->get_pixel(2 * i + 1, 2 * j + 1)
				    * 0.0625 
				    : pixel(0, 0, 0)))

			     /

				( ((i > 0 && j > 0) 
				    ? weights->get_pixel(2 * i - 1, 2 * j - 1)
				    * 0.0625
				    : pixel(0, 0, 0))
				+ ((i > 0)
				    ? weights->get_pixel(2 * i - 1, 2 * j)
				    * 0.125
				    : pixel(0, 0, 0))
				+ ((i > 0 && j < is->width() - 1)
				    ? weights->get_pixel(2 * i - 1, 2 * j + 1)
				    * 0.0625
				    : pixel(0, 0, 0))
				+ ((j > 0)
				    ? weights->get_pixel(2 * i, 2 * j - 1)
				    * 0.125
				    : pixel(0, 0, 0))
				+ weights->get_pixel(2 * i, 2 * j) 
				* 0.25
				+ ((j < is->width() - 1)
				    ? weights->get_pixel(2 * i, 2 * j + 1)
				    * 0.125
				    : pixel(0, 0, 0))
				+ ((i < is->height() - 1 && j > 0)
				    ? weights->get_pixel(2 * i + 1, 2 * j - 1)
				    * 0.0625
				    : pixel(0, 0, 0))
				+ ((i < is->height() - 1)
				    ? weights->get_pixel(2 * i + 1, 2 * j)
				    * 0.125
				    : pixel(0, 0, 0))
				+ ((i < is->height() && j < is->width() - 1)
				    ? weights->get_pixel(2 * i + 1, 2 * j + 1)
				    * 0.0625
				    : pixel(0, 0, 0)) ) );

			for (int k = 0; k < 3; k++)
				if (!finite(value[k]))
					value[k] = 0;

			is->set_pixel(i, j, value);
		}

		is->_offset = point(_offset[0] * f, _offset[1] * f);

		return is;
	}

	/*
	 * Scale an image definition array by 1/2.
	 *
	 * ALE considers an image definition array as a special kind of image
	 * weight array (typedefs of which should appear below the definition
	 * of this class).  ALE uses nonzero pixel values to mean 'defined' and
	 * zero values to mean 'undefined'.  Through this interpretation, the
	 * image weight array implementation that ALE uses allows image weight
	 * arrays to also serve as image definition arrays.
	 *
	 * Whereas scaling of image weight arrays is not generally obvious in
	 * either purpose or method, ALE requires that image definition arrays
	 * be scalable.  (Note that in the special case where weight is treated
	 * as certainty, using a geometric mean is probably correct.)
	 *
	 * We currently use a geometric mean to implement scaling of
	 * definition arrays.
	 */

	class defined_scale_by_half_threaded : public thread::decompose_domain {
		image *is;
		const image *iu;
	protected:
		void subdomain_algorithm(unsigned int thread,
				int i_min, int i_max, int j_min, int j_max) {

#if 0
			ale_real _0625 = (ale_real) 0.0625;
			ale_real _125  = (ale_real) 0.125;
			ale_real _25   = (ale_real) 0.25;
#endif

			int ui_min = (int) i_min;
			int ui_max = (int) i_max;
			int uj_min = (int) j_min;
			int uj_max = (int) j_max;

			for (int i = ui_min; i < ui_max; i++)
			for (int j = uj_min; j < uj_max; j++) {

#if 0

				/*
				 * Calculate a geometric mean; this approach
				 * may be expensive on some platforms (e.g.,
				 * those without floating-point support in
				 * hardware).
				 */

				pixel value = pixel

				      ( ( ((i > 0 && j > 0) 
					    ? ppow(iu->get_pixel(2 * i - 1, 2 * j - 1), _0625) 
					    : pixel(0, 0, 0))
					* ((i > 0)
					    ? ppow(iu->get_pixel(2 * i - 1, 2 * j), _125)
					    : pixel(0, 0, 0))
					* ((i > 0 && j < is->width() - 1)
					    ? ppow(iu->get_pixel(2 * i - 1, 2 * j + 1), _0625)
					    : pixel(0, 0, 0))
					* ((j > 0)
					    ? ppow(iu->get_pixel(2 * i, 2 * j - 1), _125)
					    : pixel(0, 0, 0))
					* ppow(iu->get_pixel(2 * i, 2 * j), _25) 
					* ((j < is->width() - 1)
					    ? ppow(iu->get_pixel(2 * i, 2 * j + 1), _125)
					    : pixel(0, 0, 0))
					* ((i < is->height() - 1 && j > 0)
					    ? ppow(iu->get_pixel(2 * i + 1, 2 * j - 1), _0625)
					    : pixel(0, 0, 0))
					* ((i < is->height() - 1)
					    ? ppow(iu->get_pixel(2 * i + 1, 2 * j), _125)
					    : pixel(0, 0, 0))
					* ((i < is->height() && j < is->width() - 1)
					    ? ppow(iu->get_pixel(2 * i + 1, 2 * j + 1), _0625)
					    : pixel(0, 0, 0))));
#else

				pixel value = iu->get_pixel(2 * i, 2 * j);

				for (int ii = 2 * i - 1; ii <= 2 * i + 1; ii++)
				for (int jj = 2 * j - 1; jj <= 2 * j + 1; jj++) {
					if (ii < 0
					 || jj < 0
					 || ii > (int) iu->height() - 1
					 || jj > (int) iu->height() - 1)
						continue;

					pixel value2 = iu->get_pixel(ii, jj);

					for (int k = 0; k < 3; k++)
						if (value2[k] < value[k]
						 || !finite(value2[k]))   /* propagate non-finites */
							value[k] = value2[k];
				}

#endif


				for (int k = 0; k < 3; k++)
					if (!finite(value[k]))
						value[k] = 0;

				is->set_pixel(i, j, value);
			}
		}

	public:
		defined_scale_by_half_threaded(image *_is, const image *_iu) 
					: decompose_domain(0, _is->height(),
					                   0, _is->width()) {
			is = _is;	
			iu = _iu;
		}
	};

	image *defined_scale_by_half(const char *name) const {
		ale_pos f = 0.5;

		image *is = scale_generator(
			(int) floor(height() * (double) f), 
			(int) floor(width()  * (double) f), depth(), name);

		assert(is);

		defined_scale_by_half_threaded dsbht(is, this);
		dsbht.run();

		is->_offset = point(_offset[0] * f, _offset[1] * f);

		return is;
	}

	/*
	 * Return an image scaled by some factor != 1.0, using bilinear
	 * interpolation.
	 */
	image *scale(ale_pos f, const char *name, int defined = 0) const {

		/*
		 * We probably don't want to scale images by a factor of 1.0,
		 * or by non-positive values.
		 */
		assert (f != 1.0 && f > 0);

		if (f > 1.0) {
			image *is = scale_generator(
				(int) floor(height() * (double) f), 
				(int) floor(width()  * (double) f), depth(), name);

			assert(is);

			unsigned int i, j, k;

			for (i = 0; i < is->height(); i++)
			for (j = 0; j < is->width(); j++) 
			for (k = 0; k < is->depth(); k++)
				is->set_pixel(i, j,
					get_scaled_bl(point(i, j), f, defined));

			is->_offset = point(_offset[0] * f, _offset[1] * f);

			return is;
		} else if (f == 0.5) {
			if (defined == 0)
				return scale_by_half(name);
			else
				return defined_scale_by_half(name);
		} else {
			image *is = scale(2*f, name, defined);
			image *result = is->scale(0.5, name, defined);
			delete is;
			return result;
		}

	}

	/*
	 * Extend the image area to the top, bottom, left, and right,
	 * initializing the new image areas with black pixels.  Negative values
	 * shrink the image.
	 */
	virtual image *_extend(int top, int bottom, int left, int right) = 0;

	static void extend(image **i, int top, int bottom, int left, int right) {
		image *is = (*i)->_extend(top, bottom, left, right);

		if (is != NULL) {
			delete (*i);
			*i = is;
		}
	}

	/*
	 * Clone 
	 */
	image *clone(const char *name) const {
		image *ic = scale_generator(
			height(), width(), depth(), name);

		assert(ic);

		for (unsigned int i = 0; i < height(); i++)
		for (unsigned int j = 0; j < width();  j++)
			ic->set_pixel(i, j, 
				get_pixel(i, j));


		ic->_offset = _offset;

		ic->_apm_memo = _apm_memo;
		ic->_acm_memo = _acm_memo;
		ic->_accm_memo = _accm_memo;
		ic->_apm = _apm;
		ic->_acm = _acm;
		ic->_accm = _accm;

		return ic;
	}

	/*
	 * Calculate the average (mean) clamped magnitude of a channel across
	 * all pixels in an image.  The magnitude is clamped to the range of 
	 * real inputs.
	 */
	ale_real avg_channel_clamped_magnitude(unsigned int k) const {

		/*
		 * This is a memoized function
		 */

		assert (k < _depth);

		avg_channel_clamped_magnitude_memo();
		return _accm[k];
	}

	pixel avg_channel_clamped_magnitude() const {
		avg_channel_clamped_magnitude_memo();
		return _accm;
	}

	/*
	 * Calculate the average (mean) magnitude of a channel across all
	 * pixels in an image.
	 */
	ale_real avg_channel_magnitude(unsigned int k) const {

		/*
		 * This is a memoized function
		 */

		assert (k < _depth);

		avg_channel_magnitude_memo();
		return _acm[k];
	}

	pixel avg_channel_magnitude() const {
		avg_channel_magnitude_memo();
		return _acm;
	}

	/*
	 * Calculate the average (mean) magnitude of a pixel (where magnitude
	 * is defined as the mean of the channel values).
	 */
	ale_real avg_pixel_magnitude() const {
		unsigned int i, j, k;

		ale_accum accumulator;
		ale_accum divisor = 0;

		if (_apm_memo)
			return _apm;

		_apm_memo = 1;
		accumulator = 0;

		for (i = 0; i < _dimy;  i++)
		for (j = 0; j < _dimx;  j++) {
			pixel value = get_pixel(i, j);

			for (k = 0; k < _depth; k++) 
				if (finite(value[k])) {
					accumulator += value[k];
					divisor += 1;
				}
		}

		accumulator /= divisor;

		_apm = accumulator;

		return _apm;
	}

	virtual ~image() {
	}
};

#endif