// Emacs: -*- C++ -*-

//
//	Copyright 1993, Center for Computer Vision and Visualization,
//	University of Florida.  All rights reserved.
//


//
// $Log: CoreImage.c,v $
// Revision 1.10.1.2  1995/01/13  19:37:42  thoth
// CoreImage has been re-merged with Image.  All special stuff should be
// friends of FBI now.
//
// Revision 1.10.1.1  1994/12/28  15:47:07  thoth
// Conversion to classes for throw.
// restrict and extend are now ROFs in AuxImage.c.
//
// Revision 1.10  1994/09/16  14:57:11  thoth
// more DOS-inspired renaming.
//
// Revision 1.9  1994/03/14  15:48:33  thoth
// ValueSet has been replaced by container class Set.
//
// Revision 1.8  1994/02/26  16:02:08  thoth
// operator() of the container now relies upon the contained classes
// to check if p is contained in the domain.  This means that users
// could possibly write an erroneous operator().
//
// Revision 1.7  1994/01/31  15:40:44  thoth
// chi_contains and compositions have been moved into AuxImage.c.
//
// Revision 1.6  1994/01/07  15:10:00  thoth
// Image class is now CoreImage and named image types are
// Image<P,T>.
//
// Revision 1.5  1993/12/29  16:55:40  thoth
// Improved error handling.
// Pixel array constructor.
// pointarray valuearray constructor.
// Added FunctionI that accepts functions with simple T type instead of
// const T&.
// new chi_contains function.
//
// Revision 1.4  1993/11/17  18:12:59  thoth
// extensivep is now extensive
// translated is now translate.
//
// Revision 1.3  1993/10/05  19:33:08  thoth
// Move to new Closure naming convention.
//
// Revision 1.2  1993/09/27  15:50:21  thoth
// Some diagnostics for abort()s that need to be ia_throws.
// xlated is now translated.
//
// Revision 1.1  1993/09/15  12:53:48  thoth
// Initial revision
//

#include <iostream.h>
#include <string.h>

#include "Point.h"
#include "Set.h"
#include "SetStructure.h"
// if included here, DIIters screw up
//#include "ClosurePS.h"

#include "Image.h"

//#include "LazyDI.h"

#include "ImageIter.h"
// #include "IPIter.h"

#include "ErrorI.h"
#include "ConstI.h"
#include "VectorI.h"
#include "ClosureI.h"
#include "FuncI.h"

#include "ClosurePS.h"

#include "image_errors.h"


template <class P,class T>
void IA_Image<P,T>::set_and_reference_image(IA_BaseImage<P,T> *i)
{
    i->incr_ref();
    this->bip = i;
}

template < class P,class T>
void IA_Image<P,T>::disassociate_image()
//
// Adjusts the reference count of the BaseDI associated with *this.
// If this is the only reference to the BaseDI, it is deleted.
//
{
    if (this->bip->decr_ref() <= 0) {
	delete this->bip;
    }
    this->bip = 0;
}


//
//
//


template < class P,class T>
IA_Image<P,T>::~IA_Image()
// clean up the storage associated with this image
{
    this->disassociate_image();
}


template < class P,class T>
IA_Image<P,T>::IA_Image()
{
    set_and_reference_image(new IA_ErrorI<P,T>);
}


template < class P,class T>
IA_Image<P,T>::IA_Image(const IA_Image<P,T>& im)
//
// Copy ctor  increments the reference count of the BaseDI associated
// with the image being copied.
//
// im		UNmodified
//
{
    set_and_reference_image(im.bip);
}


template < class P,class T>
IA_Image<P,T>::IA_Image(const IA_Set<P>& p, const T& value)
//
// Constant valued image ctor
//
// value	UNmodified
// p		UNmodified
//
{
    set_and_reference_image(new IA_ConstI<P,T>(p, value));
}


template < class P,class T>
IA_Image<P,T>::IA_Image(const IA_Set<P>& p, T *values,
			unsigned sz, int giveaway)
//
// Vector image constructor
//	Given values (a vector with elements of type T), IntPointSet p,
//	and sz (the number of Points in p), this creates a new image which
//	associates the elements in values with the points in p in the
//	iteration order of points p.  If the flag giveaway is true, we now
//	own the storage for values (and we give it to the VectorDI).
//
// p		UNmodified
{
    this->bip=0;	// in case we throw

    if (!p.extensive()){
	ia_throw(Pointset_RequireExtensive_Exception(__FILE__, __LINE__));
	set_and_reference_image(new IA_ErrorI<P,T>);
	return;
    }

    if (sz != p.card()) {
	ia_throw(Image_VectorConstructorSizeMismatch_Exception(__FILE__, __LINE__));
	set_and_reference_image(new IA_ErrorI<P,T>);
	return;
    }
    set_and_reference_image(new IA_VectorI<P,T>(p, values, sz, giveaway));
}


template < class P,class T>
IA_Image<P,T>::IA_Image(const IA_Set<P>& p,
			const T *values, unsigned sz)
//
// Vector image constructor
//	Given values (a vector with elements of type T), IntPointSet p,
//	and sz (the number of Points in p), this creates a new image which
//	associates the elements in values with the points in p in the
//	iteration order of points p.
//
// p		UNmodified
{
    this->bip=0;	// in case we throw

    if (!p.extensive()){
	ia_throw(Pointset_RequireExtensive_Exception(__FILE__, __LINE__));
	set_and_reference_image(new IA_ErrorI<P,T>);
	return;
    }

    if (sz != p.card()) {
	cerr << "image vector constructor, pointset size and array size mismatch"<<endl;
	ia_throw(Image_VectorConstructorSizeMismatch_Exception(__FILE__, __LINE__));
	set_and_reference_image(new IA_ErrorI<P,T>);
	return;
    }

    // We can cast away constness because we know the VectorDI
    // constructor won\'t hose us.
    set_and_reference_image(new IA_VectorI<P,T>(p, (T*)values, sz, 0));
}


//
//
//

template < class P,class T>
IA_VectorI<P,T> *
IA_Image<P,T>::sortedPixel_to_VectorI(const IA_Pixel<P,T> *pixels,
				      unsigned count)
{
    P * const	points = new P[count];
    T * const	vals = new T[count];
    for (int i=0; i<count; i++) {
	points[i] = pixels[i].point;
	vals[i] = pixels[i].value;
    }


    IA_Set<P>	ps(points[0].dim(), points, count);

    IA_VectorI<P,T> *	rval = new IA_VectorI<P,T>(ps, vals, count, 1);

    delete[] points;
    /* we've given the values away */

    return rval;
}

//
// constructor from Pixel arrays to images
//
template < class P,class T>
IA_Image<P,T>::IA_Image(const IA_Pixel<P,T> *pixels, unsigned count)
{
    // gotta check, dimensions, and what if they don't give us any pixels?
    // clean up later.

    IA_Pixel<P,T> * const	copy = new IA_Pixel<P,T>[count];
    for (int i=0; i<count; i++)
	copy[i] = pixels[i];

    qsort((char*)copy, count, sizeof(*copy), IA_Pixel<P,T>::vcompare);

    set_and_reference_image(sortedPixel_to_VectorI(copy, count));

    delete[] copy;
}


//
// constructor from an array of points with a corresponding array of values
//

template < class P,class T>
IA_Image<P,T>::IA_Image(const P *points, unsigned pcount,
			const T *values, unsigned vcount)
{
    if (pcount!=vcount)  {
	ia_throw(Image_VectorConstructorSizeMismatch_Exception(__FILE__, __LINE__));
	set_and_reference_image(new IA_ErrorI<P,T>);
	return;
    }

    if (pcount<1) {
	ia_throw(Image_PointsetTooSmall_Exception(__FILE__, __LINE__));
	set_and_reference_image(new IA_ErrorI<P,T>);
	return;
    }

    int	dim=points[0].dim();

    for (int i=1; i<pcount; i++)
	if (points[i].dim() != dim)
	    break;

    if (i<pcount) {
	ia_throw(Pointset_DimensionMismatch_Exception(__FILE__, __LINE__));
	set_and_reference_image(new IA_ErrorI<P,T>);
	return;
    }

    for (i=0; i+1<pcount; i++)
	if (pointcmp(points[i], points[i+1])>=0)
	    break;		// if it\'s not strictly increasing

    if (i+1<pcount) {
	// the point array is not sorted
	IA_Pixel<P,T> * const	pixels = new IA_Pixel<P,T>[pcount];
	for (int i=0; i<pcount; i++) {
	    pixels[i].point = points[i];
	    pixels[i].value = values[i];
	}

	qsort((char*)pixels, pcount, sizeof(*pixels), IA_Pixel<P,T>::vcompare);

	set_and_reference_image(sortedPixel_to_VectorI(pixels, pcount));

	delete[] pixels;
    } else {
	// the point array IS sorted
	IA_Set<P>	ps(dim, points, pcount);
	set_and_reference_image(new IA_VectorI<P,T>(ps, (T*)values, vcount,
						    0));
    }
}

template < class P,class T>
IA_Image<P,T>::IA_Image(const P *points, unsigned pcount,
			T *values, unsigned vcount, int giveaway)
{
    if (pcount!=vcount)  {
	ia_throw(Image_VectorConstructorSizeMismatch_Exception(__FILE__, __LINE__));
	set_and_reference_image(new IA_ErrorI<P,T>);
	return;
    }

    if (pcount<1) {
	IA::not_yet_implemented(__FILE__, __LINE__);
	abort();		// should ia_throw
	set_and_reference_image(new IA_ErrorI<P,T>);
    }

    int	dim=points[0].dim();

    for (int i=1; i<pcount; i++)
	if (points[i].dim() != dim)
	    break;

    if (i<pcount) {
	ia_throw(Pointset_DimensionMismatch_Exception(__FILE__, __LINE__));
	set_and_reference_image(new IA_ErrorI<P,T>);
	return;
    }

    for (i=0; i+1<pcount; i++)
	if (pointcmp(points[i], points[i+1])>=0)
	    break;		// if it\'s not strictly increasing

    if (i+1<pcount) {
	// the point array is not sorted
	IA_Pixel<P,T> * const	pixels = new IA_Pixel<P,T>[pcount];
	for (int i=0; i<pcount; i++) {
	    pixels[i].point = points[i];
	    pixels[i].value = values[i];
	}

	qsort((char*)pixels, pcount, sizeof(*pixels), IA_Pixel<P,T>::vcompare);

	set_and_reference_image(sortedPixel_to_VectorI(pixels, pcount));

	delete[] pixels;
	if (giveaway)
	    delete[] values;	// well, we couldn\'t use them as-is
    } else {
	// the point array IS sorted
	IA_Set<P>	ps(dim, points, pcount);
	set_and_reference_image(new IA_VectorI<P,T>(ps, values, vcount,
						    giveaway));
    }
}

template < class P,class T>
IA_Image<P,T>::IA_Image(const IA_ClosureI<P,T>&i)
// User-defined-image constructor.  They provide us with an object,
// we clone a copy onto the heap and point to it.
{
    set_and_reference_image(i.clone_self());
}

template < class P,class T>
IA_Image<P,T>::IA_Image(const IA_Set<P> &ps,
			T (*f)(const P&))
// User-defined-image constructor.  They provide us with a function
// and we wrap an implementation around it.
{
    set_and_reference_image(new IA_FunctionRefI<P,T>(ps,f));
}

template < class P,class T>
IA_Image<P,T>::IA_Image(const IA_Set<P> &ps,
			T (*f)( P))
// User-defined-image constructor.  They provide us with a function
// and we wrap an implementation around it.
{
    set_and_reference_image(new IA_FunctionI<P,T>(ps,f));
}

template < class P,class T>
IA_Image<P,T>&
IA_Image<P,T>::operator=(const IA_Image<P,T>& im)
//
// Image assignment operator
//
// im		UNmodified (however, the reference count of the image
//		       	referenced by im will be modified)
//
{
    IA_BaseImage<P,T> *temp = im.bip;
    temp->incr_ref();

    this->disassociate_image();

    this->bip = temp;

    return *this;
}

template < class P,class T>
IA_Image<P,T>& IA_Image<P,T>::operator =(const T& value)
//
// Image assignment by constant value
//
// value	UNmodified
//
{
    IA_Set<P> pset;
    pset = this->domain();
    this->disassociate_image();
    set_and_reference_image(new IA_ConstI<P,T>(pset, value));
    return *this;
}

template < class P,class T>
T& IA_Image<P,T>::operator [](const P& p)
//
// Return a reference to a pixel value.  The user could potentially
// modify the value through the reference.  Since images are reference
// counted we use a copy-on-(potential)write scheme.
//
{
#if 0
    if (!this->domain().contains(p)) {
	ia_throw(IA::IMAGE_POINT_OUTSIDE_DOMAIN, __FILE__, __LINE__);
    }
#endif
    if (this->bip->get_ref() > 1 || !this->bip->extensive()) {

	// If reference count > 1, we need to clone the base image to avoid
	//	multiple writers.
	// If bip is not extensive, we need to clone the base image to have
	//	an address we can return.

	IA_BaseImage<P,T>	*temp = bip->clone_self_extensively();

	this->disassociate_image();	// We don\'t need this anymore.

	set_and_reference_image(temp);
    }
    return (*this->bip)[p];
}

#if 0

template <class P,class T>
IA_Image<IA_Point<int>,T> IA_Image<P,T>::restrict(const IA_Set<IA_Point<int> >& ips)
    const
{
    IA_BaseImage<IA_Point<int>,T>	*rval = bip->restricted_to(ips);
    return rval;
}

template <class P,class T>
IA_Image<P,T> IA_Image<P,T>::restrict(const IA_Set<IA_Point<double> >& ips)
    const
{
    IA_BaseImage<P,T>	*rval = bip->restricted_to(ips);
    if (rval)
	return rval;
    else
	return *this;
}

template < class P,class T>
IA_Image<P,T> IA_Image<P,T>::restrict(const IA_Set<T>& vs)
    const
{
    if (!this->domain().extensive()) {
	ia_throw( Pointset_RequireExtensive_Exception(__FILE__, __LINE__));
	return IA_Image<P,T>();
    }
    
    const	unsigned sz = this->domain().card();
    P *point_vec = new P[sz];
    T *value_vec = new T[sz];
    
    IA_IPIter<P,T> iter(*this);
    P point;
    T value;
    int new_size=0;
    
    // Copy all points in *this\'s domain that have values in the argument
    // Set.
    while(iter(point,value)){
	if (vs.contains(value)){
	    value_vec[new_size] = value;
	    point_vec[new_size] = point;
	    new_size++;
	}
    }
    
    IA_Image<P,T> rval = IA_Image<P,T>
	(IA_Set<P>(this->domain().dim(), point_vec, new_size),
	 value_vec, new_size, 1);
    
    delete[] point_vec;
    
    return rval;
}

#endif

//
// The following functions implement image extension
//

template <class P, class T>
void IA_Image<P,T>::vec_vec_2_vec_structcopy
(const T **lhs_data, const T **rhs_data, T **dest, const IA_SetStructure &ss)
//
// Copy values from *lhs_data and *rhs_data into *dest according to
// SetStructure ss.  If there is a value from both, the **lhs_data
// is preferred.
//
// *lhs_data	MODIFIED
// *rhs_data	MODIFIED
// *dest	MODIFIED
// ss		not modified
{
    for(int i = 0; i < ss.nintervals(); i++) {
	IA_ss_interval ssi = ss.retrieve_interval(i);
	if (ssi.substructure == IA_SetStructure::BOTH) {
	    for (int j=0; j<ssi.count; j++) {
		*((*dest)++) = *((*lhs_data)++);
		(*rhs_data)++;
	    }
	} else if (ssi.substructure == IA_SetStructure::FIRST_ONLY) {
	    for (int j=0; j<ssi.count; j++) {
		*((*dest)++) = *((*lhs_data)++);
	    }
	} else if (ssi.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (int j=0; j<ssi.count; j++) {
		*((*dest)++) = *((*rhs_data)++);
	    }
	} else {
	    for (int j=0; j<ssi.count; j++) {
		vec_vec_2_vec_structcopy(lhs_data, rhs_data, dest,
					 ssi.substructure);
	    }
	}
    }
}

template <class P, class T>
void IA_Image<P,T>::vec_iter_2_vec_structcopy
(
 const T	**lhs_data,
#ifdef FUNKY_ERROR
 void *rhs_v,
#else
 IA_IVIter<P,T> *rhs_iter,
#endif
 T **dest, const IA_SetStructure &ss)
//
// Copy values from *lhs_data and *rhs_iter into *dest according to
// SetStructure ss.  If there is a value from both, the lhs data
// is preferred.
//
// *lhs_data	MODIFIED
// *rhs_iter	MODIFIED
// *dest	MODIFIED
// ss		not modified
{
#ifdef FUNKY_ERROR
    IA_IVIter<P,T>	*rhs_iter = (IA_IVIter<P,T> *)rhs_v;
#endif
    T	temp;
    for(int i = 0; i < ss.nintervals(); i++) {
	IA_ss_interval ssi = ss.retrieve_interval(i);
	if (ssi.substructure == IA_SetStructure::BOTH) {
	    for (int j=0; j<ssi.count; j++) {
		if (!(*rhs_iter)(temp))
		    IA::internal_error(__FILE__, __LINE__);
		*((*dest)++) = *((*lhs_data)++);
	    }
	} else if (ssi.substructure == IA_SetStructure::FIRST_ONLY) {
	    for (int j=0; j<ssi.count; j++) {
		*((*dest)++) = *((*lhs_data)++);
	    }
	} else if (ssi.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (int j=0; j<ssi.count; j++) {
		if (!(*rhs_iter)(temp))
		    IA::internal_error(__FILE__, __LINE__);
		*((*dest)++) = temp;
	    }
	} else {
	    for (int j=0; j<ssi.count; j++) {
		vec_iter_2_vec_structcopy(lhs_data, rhs_iter, dest,
					   ssi.substructure);
	    }
	}
    }
}

template <class P, class T>
void IA_Image<P,T>::iter_vec_2_vec_structcopy
(
#ifdef FUNKY_ERROR
 void *lhs_v,
#else
 IA_IVIter<P,T> *lhs_iter,
#endif
 const T**rhs_data,
 T **dest, const IA_SetStructure &ss)
//
// Copy values from *lhs_iter and *rhs_data into *dest according to
// SetStructure ss.  If there is a value from both, the lhs data
// is preferred.
//
// *lhs_iter	MODIFIED
// *rhs_data	MODIFIED
// *dest	MODIFIED
// ss		not modified
{
#ifdef FUNKY_ERROR
    IA_IVIter<P,T>	*lhs_iter = (IA_IVIter<P,T> *)lhs_v;
#endif
    T	temp;
    for(int i = 0; i < ss.nintervals(); i++) {
	IA_ss_interval ssi = ss.retrieve_interval(i);
	if (ssi.substructure == IA_SetStructure::BOTH) {
	    for (int j=0; j<ssi.count; j++) {
		(*rhs_data)++;
		if (!(*lhs_iter)(temp))
		    IA::internal_error(__FILE__,__LINE__);
		*((*dest)++) = temp;
	    }
	} else if (ssi.substructure == IA_SetStructure::FIRST_ONLY) {
	    for (int j=0; j<ssi.count; j++) {
		if (!(*lhs_iter)(temp))
		    IA::internal_error(__FILE__,__LINE__);
		*((*dest)++) = temp;
	    }
	} else if (ssi.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (int j=0; j<ssi.count; j++) {
		*((*dest)++) = *((*rhs_data)++);
	    }
	} else {
	    for (int j=0; j<ssi.count; j++) {
		iter_vec_2_vec_structcopy(lhs_iter, rhs_data, dest,
					   ssi.substructure);
	    }
	}
    }
}

template <class P, class T>
void IA_Image<P,T>::iter_iter_2_vec_structcopy
(
#ifdef FUNKY_ERROR
 void *lhs_v, void *rhs_v,
#else
 IA_IVIter< IA_IVIter<T,T> *lhs_iter> *rhs_iter,
#endif
 T **dest, const IA_SetStructure &ss)
//
// Copy values from *lhs_iter and *rhs_iter into *dest according to
// SetStructure ss.  If there is a value from both, the lhs data
// is preferred.
//
// *lhs_iter	MODIFIED
// *rhs_iter	MODIFIED
// *dest	MODIFIED
// ss		not modified
{
#ifdef FUNKY_ERROR
    IA_IVIter<P,T>	*lhs_iter = (IA_IVIter<P,T> *)lhs_v;
    IA_IVIter<P,T>	*rhs_iter = (IA_IVIter<P,T> *)rhs_v;
#endif
    T	temp;
    for(int i = 0; i < ss.nintervals(); i++) {
	IA_ss_interval ssi = ss.retrieve_interval(i);
	if (ssi.substructure == IA_SetStructure::BOTH) {
	    for (int j=0; j<ssi.count; j++) {
		if (!(*rhs_iter)(temp))
		    IA::internal_error(__FILE__,__LINE__);
		if (!(*lhs_iter)(temp)) // throw away the rhs value
		    IA::internal_error(__FILE__,__LINE__);
		*((*dest)++) = temp;
	    }
	} else if (ssi.substructure == IA_SetStructure::FIRST_ONLY) {
	    for (int j=0; j<ssi.count; j++) {
		if (!(*lhs_iter)(temp))
		    IA::internal_error(__FILE__,__LINE__);
		*((*dest)++) = temp;
	    }
	} else if (ssi.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (int j=0; j<ssi.count; j++) {
		if (!(*rhs_iter)(temp))
		    IA::internal_error(__FILE__,__LINE__);
		*((*dest)++) = temp;
	    }
	} else {
	    for (int j=0; j<ssi.count; j++) {
		iter_iter_2_vec_structcopy(lhs_iter, rhs_iter, dest,
					   ssi.substructure);
	    }
	}
    }
}


#if 0 // implemented elsewhere

template <class P, class T>
IA_Image<P,T>
IA_Image<P,T>::extend(const IA_Image<P,T>& imr) const
//
// Returns the extension of image *this to image imr.
//
// imr		UNmodified
//
{
    if(!this->domain().extensive() || !imr.domain().extensive()) {
	// Need LazyExtensionImage
	ia_throw(Pointset_RequireExtensive_Exception(__FILE__, __LINE__));
	return IA_Image<P,T>();
    }

    IA_SetStructure ss;
    IA_Set<P> newps = union_with_structure(this->domain(),
					   imr.domain(), &ss);
    const unsigned	sz = newps.card();
    T	* const dest = new T[sz];
    T	*scan = dest;

    if (this->type()==IA_VectorI<P,T>::s_type()) {
	const IA_VectorI<P,T> *lhs = (IA_VectorI<P,T>*)this->bip;
	const T	*lhs_data = lhs->vec;

	if (imr.type()==IA_VectorI<P,T>::s_type()) {
	    const T	*rhs_data = ((IA_VectorI<P,T>*)imr.bip)->vec;
	    vec_vec_2_vec_structcopy(&lhs_data, &rhs_data, &scan, ss);
	} else {
	    IA_IVIter<P,T>	rhs_iter(imr);
	    vec_iter_2_vec_structcopy(&lhs_data, &rhs_iter, &scan, ss);
	}
    } else {
	IA_IVIter<P,T>	lhs_iter(*this);

	if (imr.type()==IA_VectorI<P,T>::s_type()) {
	    const T	*rhs_data = ((IA_VectorI<P,T>*)imr.bip)->vec;
	    iter_vec_2_vec_structcopy(&lhs_iter, &rhs_data, &scan, ss);
	} else {
	    IA_IVIter<P,T>	rhs_iter(imr);
	    iter_iter_2_vec_structcopy(&lhs_iter, &rhs_iter, &scan, ss);
	}
    }

    return IA_Image<P,T>(newps, dest, sz, 1);
}

#endif // NYI

//
//
//

#if 1
template < class P,class T>
IA_Image<P,T> IA_Image<P,T>::translate(const P &p) const
{
    return bip->xlated_by(p);
}
#endif

//
//
//

template < class P,class T>
IA_Set<T> IA_Image<P,T>::range() const
{
#if 1
    return bip->range();
#else
    unsigned	sz;
    T	*vec = bip->range(&sz);
    IA_Set<T>	rval(vec,sz);
    delete[] vec;
    return rval;
#endif
}

template < class P,class T>
T IA_Image<P,T>::reduce(T (*gamma)(const T&, const T&), T identity) const
//
// Reduces *this with function gamma.
//
{
    T result(identity), new_value;

    if (this->card()>0) {
	IA_IVIter<P,T> iter(*this);

	iter(result);
	while(iter(new_value)) {
	    result = gamma( result, new_value);
	}
    }
    return result;
}

template < class P,class T>
T IA_Image<P,T>::reduce(T (*gamma) (T, T), T identity) const
//
// Reduces *this with function gamma.
//
{
    T result(identity), new_value;

    if (this->card()>0) {
	IA_IVIter<P,T> iter(*this);

	iter(result);
	while(iter(new_value)) {
	    result = gamma( result, new_value);
	}
    }
    return result;
}

//
//
//

#include "ImgRestrictAssn.c"
