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

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


static char DiscreteImagercsid[] = "$Id: DiscreteImage.c,v 1.14 93/06/10 17:23:59 thoth Exp $";

// $Log:	DiscreteImage.c,v $
// Revision 1.14  93/06/10  17:23:59  thoth
// Workaround for 3.0.1 bug.
// 
// Revision 1.13  93/05/27  11:28:22  thoth
// Copyright Notices
// 
// Revision 1.12  93/05/26  16:59:56  thoth
// Read and write of IAA image files.
// 
// Revision 1.11  93/05/18  21:39:00  thoth
// BaseDI range method returns ValueSet again.
// New Point-Point composition closure object.
// 
// Revision 1.10  93/04/29  11:18:38  thoth
// Support for point-point-function and image composition.
// 
// Revision 1.9  93/04/08  13:10:49  thoth
// removed cruft.
// restriction is now a virtual method of the BaseDI.  We don't
// need that goofy case statement any more.
// translation is also a virtual method.
// 
// Revision 1.8  93/03/10  13:32:36  thoth
// range method no longer returns valueset due to
// difficulties with template implementation.
// FunctionDI constructor was not truly generic.
// 
// Revision 1.7  93/02/23  14:10:14  thoth
// set_and_reference_image method encodes common initialization semantics.
// User-defined OtherDIs are now supported.
// 
// Revision 1.6  93/02/20  12:55:06  thoth
// Function images are no longer supported.
// ErrorDIs are created dynamically as needed.
// This is conceptually convenient.
// The pointset should go first in an image constructor.
// Vector constructors now can give away the value vector.
// Image restriction and extension should be fast now.  I can't
// be certain they work though...
// 
// Revision 1.5  93/02/17  11:00:18  thoth
// We now have the ability to give away the value vector.
// 
// Revision 1.4  93/02/08  12:57:57  thoth
// virtual destructor is now outlined.
// 
// Revision 1.3  93/01/19  14:56:03  jnw
// Modified to support IA_ naming convention
// 
// Revision 1.2  92/11/12  16:59:09  jnw
// Added restriction, extension, iteration, etc.
// 
// Revision 1.1  92/10/05  11:42:42  jnw
// Initial revision
// 

#include <stream.h>
#include	<netinet/in.h>	// for ntohl
#include <string.h>

#include "SetStructure.h"
// if included here, DIIters screw up
//#include "PredicatePS.h"

#include "DiscreteImage.h"

//#include "LazyDI.h"
#include "DIVIter.h"
#include "DIPIter.h"

#include "ErrorDI.h"
#include "ConstDI.h"
#include "VectorDI.h"
#include "OtherDI.h"
#include "FunctionRefDI.h"

#include "PredicatePS.h"

// I wish this worked
//template <class T> IA_ErrorDI<T> IA_DiscreteImage<T>::errorDI;


template <class T>
void IA_DiscreteImage<T>::set_and_reference_image(IA_BaseDI<T> *i)
{
    i->incr_ref();
    this->bdip = i;
}

template <class T>
void IA_DiscreteImage<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->bdip->decr_ref() <= 0) {
	delete this->bdip;
    }
    this->bdip = 0;
}

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

template <class T>
IA_DiscreteImage<T>::IA_DiscreteImage()
{
    set_and_reference_image(new IA_ErrorDI<T>);
}

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

template <class T>
IA_DiscreteImage<T>::IA_DiscreteImage(const IA_IntPointSet& p, const T& value)
//
// Constant valued image ctor
//
// value	UNmodified
// p		UNmodified
//
{
    set_and_reference_image(new IA_ConstDI<T>(p, value));
}

template <class T>
IA_DiscreteImage<T>::IA_DiscreteImage(const IA_IntPointSet& 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->bdip=0;	// in case we throw

    if (!p.extensivep()){
	ia_throw( IA::IMAGE_STRUCTURE_NOT_SUPPORTED, __FILE__, __LINE__);
	return;
    }

    if (sz != p.card()) {
	abort(); // throw
	return;
    }
    set_and_reference_image(new IA_VectorDI<T>(p, values, sz, giveaway));
}

template <class T>
IA_DiscreteImage<T>::IA_DiscreteImage(const IA_IntPointSet& 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->bdip=0;	// in case we throw

    if (!p.extensivep()){
	ia_throw( IA::IMAGE_STRUCTURE_NOT_SUPPORTED
	      , __FILE__, __LINE__);
	return;
    }

    if (sz != p.card()) {
	abort(); // throw
	return;
    }

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

template <class T>
IA_DiscreteImage<T>::IA_DiscreteImage(const IA_OtherDI<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 T>
IA_DiscreteImage<T>::IA_DiscreteImage(const IA_IntPointSet &ps,
				      T (*f)(const IA_IntPoint&))
// User-defined-image constructor.  They provide us with a function
// and we wrap an implementation around it.
{
    set_and_reference_image(new IA_FunctionRefDI<T>(ps,f));
}

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

    this->disassociate_image();

    this->bdip = temp;

    return *this;
}

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

template <class T>
T& IA_DiscreteImage<T>::operator [](const IA_IntPoint& 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 (!this->domain().contains(p)) {
	abort();		// throw
    }
    if (this->bdip->get_ref() > 1 || !this->bdip->extensivep()) {

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

	IA_BaseDI<T>	*temp = bdip->clone_self_extensively();

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

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


template <class T>
IA_DiscreteImage<T> IA_DiscreteImage<T>::restrict(const IA_IntPointSet& ips)
    const
{
    IA_BaseDI<T>	*rval = bdip->restricted_to(ips);
    if (rval)
	return rval;
    else
	return *this;
}

template <class T>
IA_DiscreteImage<T> IA_DiscreteImage<T>::restrict(const IA_ValueSet<T>& vs)
    const
{
    if (!this->domain().extensivep()) {
	ia_throw( IA::PSET_REQUIRE_EXTENSIVE, __FILE__, __LINE__);
    }
    
    const	unsigned sz = this->domain().card();
    IA_IntPoint *point_vec = new IA_IntPoint[sz];
    T *value_vec = new T[sz];
    
    IA_DIPIter<T> iter(*this);
    IA_IntPoint point;
    T value;
    int new_size=0;
    
    // Copy all points in *this\'s domain that have values in the argument
    // ValueSet.
    while(iter(point,value)){
	if (vs.contains(value)){
	    value_vec[new_size] = value;
	    point_vec[new_size] = point;
	    new_size++;
	}
    }
    
    IA_DiscreteImage<T> rval = IA_DiscreteImage<T>
	(IA_IntPointSet(this->domain().dim(), point_vec, new_size),
	 value_vec, new_size, 1);
    
    delete[] point_vec;
    
    return rval;
}

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

template <class T>
void IA_DiscreteImage<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 T>
void IA_DiscreteImage<T>::vec_iter_2_vec_structcopy
(
 const T	**lhs_data,
#ifdef FUNKY_ERROR
 void *rhs_v,
#else
 IA_DIVIter<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_DIVIter<T>	*rhs_iter = (IA_DIVIter<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))
		    abort(); // internal error
		*((*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))
		    abort(); // internal error
		*((*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 T>
void IA_DiscreteImage<T>::iter_vec_2_vec_structcopy
(
#ifdef FUNKY_ERROR
 void *lhs_v,
#else
 IA_DIVIter<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_DIVIter<T>	*lhs_iter = (IA_DIVIter<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))
		    abort(); // internal error
		*((*dest)++) = temp;
	    }
	} else if (ssi.substructure == IA_SetStructure::FIRST_ONLY) {
	    for (int j=0; j<ssi.count; j++) {
		if (!(*lhs_iter)(temp))
		    abort(); // internal error
		*((*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 T>
void IA_DiscreteImage<T>::iter_iter_2_vec_structcopy
(
#ifdef FUNKY_ERROR
 void *lhs_v, void *rhs_v,
#else
 IA_DIVIter<T> *lhs_iter, IA_DIVIter<T> *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_DIVIter<T>	*lhs_iter = (IA_DIVIter<T> *)lhs_v;
    IA_DIVIter<T>	*rhs_iter = (IA_DIVIter<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))
		    abort(); // internal error
		if (!(*lhs_iter)(temp)) // throw away the rhs value
		    abort(); // internal error
		*((*dest)++) = temp;
	    }
	} else if (ssi.substructure == IA_SetStructure::FIRST_ONLY) {
	    for (int j=0; j<ssi.count; j++) {
		if (!(*lhs_iter)(temp))
		    abort(); // internal error
		*((*dest)++) = temp;
	    }
	} else if (ssi.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (int j=0; j<ssi.count; j++) {
		if (!(*rhs_iter)(temp))
		    abort(); // internal error
		*((*dest)++) = temp;
	    }
	} else {
	    for (int j=0; j<ssi.count; j++) {
		iter_iter_2_vec_structcopy(lhs_iter, rhs_iter, dest,
					   ssi.substructure);
	    }
	}
    }
}


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

    IA_SetStructure ss;
    IA_IntPointSet 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_VectorDI<T>::s_type()) {
	const IA_VectorDI<T> *lhs = (IA_VectorDI<T>*)this->bdip;
	const T	*lhs_data = lhs->vec;

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

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

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

//
//
//

template <class T>
IA_DiscreteImage<T> IA_DiscreteImage<T>::xlated(const IA_IntPoint &p) const
{
    return bdip->xlated_by(p);
}
//
//
//

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

template <class T>
T IA_DiscreteImage<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_DIVIter<T> iter(*this);

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

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

    if (this->card()>0) {
	IA_DIVIter<T> iter(*this);

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


template <class T>
ostream& operator <<(ostream& o, const IA_DiscreteImage<T> &im)
{
    return o << *(im.bdip);
}

//
//
//

inline int slurpint(istream &is) {
    int	k;
    is.read((char*)&k, 4);
    return ntohl(k);
}

template <class T>
IA_DiscreteImage<T> IA_DiscreteImage<T>::read_IAA(istream &header,
						      istream &values)
{
    // start parsing the header

    if ( ! ( header.get() == 0x61 &&
	     header.get() == 0x69 &&
	     header.get() == 0x41 &&
	     header.get() == 0x49) ) {
	// format error;
	ia_throw(IA::IAAIO_FORMAT_ERROR, __FILE__, __LINE__);
	return IA_DiscreteImage<T>();
    }

    // this is all ugly.  The IAA image file spec is not specced.
    //  We should probably talk about network byte order and specify the
    // sizes of the fields.

    int	k;
    k = slurpint(header);
    switch (k) {
    case 1:
	break;
    default:
	// wrong version of the file format
	ia_throw(IA::IAAIO_FORMAT_ERROR, __FILE__, __LINE__);
	return IA_DiscreteImage<T>();
    }

    k = slurpint(header);	// bits per value
    if (k!=8*sizeof(T)) {
	// wrong sort.
	ia_throw(IA::TYPE_MISMATCH, __FILE__, __LINE__);
	return IA_DiscreteImage<T>();
    }

    header.seekg(5*4, ios::cur); // skip reserved bytes

    const int nheadelem = slurpint(header);
    const unsigned valsoffset = slurpint(header);
    const int dim = slurpint(header);
    const int code = slurpint(header);
    IA_IntPointSet	domain;
    switch (code) {
    case 0:
	domain = IA_BlackHole(dim);
	// how do we handle the rest of the header?
	return IA_DiscreteImage<T>(domain, (T)0);
	break;
    case 1:
	{
	    int	*inf = new int[dim];
	    int	*sup = new int[dim];
	    for (k=0; k<dim; k++)
		inf[k] = slurpint(header);
	    for (k=0; k<dim; k++)
		sup[k] = slurpint(header);
	    domain = IA_IntPointSet(IA_IntPoint(inf, dim, 1),
				    IA_IntPoint(sup, dim, 1));
	}
	break;
    default:
	// pure or impure interval?
	// abort(); // throw
	ia_throw(IA::IAAIO_FORMAT_ERROR, __FILE__, __LINE__);
	return IA_DiscreteImage<T>();
    }

    const int card = slurpint(header);

    if (card != domain.card()) {
	// file consistency error
	// abort(); throw
	ia_throw(IA::IAAIO_FORMAT_ERROR, __FILE__, __LINE__);
	return IA_DiscreteImage<T>();
    }

    if (!header) {
	// error somewhere in the reading process
	ia_throw(IA::IAAIO_FORMAT_ERROR, __FILE__, __LINE__);
	return IA_DiscreteImage<T>();
    }

    T	*valvec = new T[card];
    values.read((char*)valvec, card * sizeof(T));

    return IA_DiscreteImage<T>(domain, valvec, card, 1);
}

template <class T>
IA_DiscreteImage<T> IA_DiscreteImage<T>::read_IAA(const char *fname1,
						      const char *fname2)
{
    ifstream header(fname1), values(fname2);
    return read_IAA(header, values);
}

template <class T>
IA_DiscreteImage<T> IA_DiscreteImage<T>::read_IAA(const char *fname)
{
    // create some istreams
    char	*temp = new char[strlen(fname)+5];
    strcpy(temp, fname);
    strcat(temp, ".iaa1");
    ifstream	header(temp);
    strcpy(temp, fname);
    strcat(temp, ".iaa2");
    ifstream	values(temp);
    delete [] temp;

    return read_IAA(header, values);
}

//
//
//


inline void spewint(ostream &is, int i) {
    int	k = htonl(i);
    is.write((char*)&k, 4);
}

template <class T>
void IA_DiscreteImage<T>::write_IAA(ostream &header, ostream &values) const
{
    header.put(0x61);		// magic number
    header.put(0x69);
    header.put(0x41);
    header.put(0x49);

    spewint(header, 1);		// version
    spewint(header, 8*sizeof(T)); // bits per value

    spewint(header, 0);		// reserved
    spewint(header, 0);
    spewint(header, 0);
    spewint(header, 0);
    spewint(header, 0);

    spewint(header, 18);	// how do we compute this?

    spewint(header, 1);		// data offset

    const int dim = domain().dim();
    spewint(header, dim);
    if (domain().empty()) {
	// this stuff is probably wrong
	spewint(header, 0);
	spewint(header, 0); // cardinality
	// we\'re done
    } else if (domain().boxy()) {
	spewint(header, 1);
	int	k;
	IA_IntPoint	p;

	p=domain().inf();
	for (k=0; k<dim; k++)
	    spewint(header, p[k]);

	p=domain().sup();
	for (k=0; k<dim; k++)
	    spewint(header, p[k]);

	const int card=domain().card();
	spewint(header, card);

	values.write((const char*)value_array().data_start(), card*sizeof(T));

    } else {
	ia_throw(IA::IAAIO_FORMAT_ERROR, __FILE__, __LINE__);
	abort(); // unhandled domain type
    }

}

template <class T>
void IA_DiscreteImage<T>::write_IAA(const char *fname1, const char *fname2)
    const
{
    ofstream header(fname1), values(fname2);
    write_IAA(header, values);
}

template <class T>
void IA_DiscreteImage<T>::write_IAA(const char *fname) const
{
    // create some ostreams
    char	*temp = new char[strlen(fname)+6];
    strcpy(temp, fname);
    strcat(temp, ".iaa1");
    ofstream	header(temp);
    strcpy(temp, fname);
    strcat(temp, ".iaa2");
    ofstream	values(temp);
    delete [] temp;

    write_IAA(header, values);
}


//
//
//

class IA_PPCompositionIPS: public IA_PSPred<IA_IntPoint> {
    IA_IntPoint	(*f)(const IA_IntPoint&);
    IA_IntPointSet	base;
  public:
    IA_PPCompositionIPS(IA_IntPoint (*f_)(const IA_IntPoint&),
			const IA_IntPointSet &base_)
    :IA_PSPred<IA_IntPoint>(base_.dim()), f(f_), base(base_) {
    }
    IA_PSPred<IA_IntPoint> *clone_self() const {
	return new IA_PPCompositionIPS(*this);
    }
    int contains(const IA_IntPoint &ip) const {
	// Uglier than Rosanne bug
	IA_IntPoint (*f2)(const IA_IntPoint&) = f;
	return base.contains(f2(ip));
    }
};

template <class T>
class IA_PPCompositionDI: public IA_OtherDI<T> {
    IA_IntPoint	(*f)(const IA_IntPoint&);
    IA_DiscreteImage<T>	im;
  public:
    IA_PPCompositionDI( IA_IntPoint (*f_)(const IA_IntPoint&),
		       const IA_DiscreteImage<T> &im_)
    :IA_OtherDI<T>(IA_PPCompositionIPS(f_, im_.domain())), f(f_), im(im_) {
    }
    T operator()(const IA_IntPoint&p) const {
	IA_IntPoint (*f2)(const IA_IntPoint &) = f;
	return im(f2(p));
    }
    IA_OtherDI<T>* clone_self() const {
	return new IA_PPCompositionDI<T>(*this);
    }
#if 0
    T *range(unsigned *n) const {
	return im.range(n);
    }
#endif
};

template <class T>
IA_DiscreteImage<T> compose(const IA_DiscreteImage<T> &i,
			    IA_IntPoint (*f)(const IA_IntPoint&))
{
    return IA_PPCompositionDI<T>(f, i);
}

//
//
//


class IA_PPCompClosureIPS: public IA_PSPred<IA_IntPoint> {
    IA_PPCompClosure<IA_IntPoint, IA_IntPoint> *f;
    IA_IntPointSet	base;
  public:
    IA_PPCompClosureIPS( const IA_PPCompClosure<IA_IntPoint, IA_IntPoint> &f_,
			const IA_IntPointSet &base_)
    :IA_PSPred<IA_IntPoint>(base_.dim()), f(f_.clone_self()), base(base_) {
    }
    IA_PPCompClosureIPS(const IA_PPCompClosureIPS &other)
    :IA_PSPred<IA_IntPoint>(other.base.dim()), f(other.f->clone_self()),
     base(other.base) {
    }
    ~IA_PPCompClosureIPS() { delete f; }

    IA_PSPred<IA_IntPoint> *clone_self() const {
	return new IA_PPCompClosureIPS(*this);
    }
    int contains(const IA_IntPoint &ip) const {
	return base.contains((*f)(ip));
    }
};

template <class T>
class IA_PPCompClosureDI: public IA_OtherDI<T> {
    IA_PPCompClosure<IA_IntPoint, IA_IntPoint> *f;
    IA_DiscreteImage<T>	im;
  public:
    IA_PPCompClosureDI( const IA_PPCompClosure<IA_IntPoint, IA_IntPoint> &f_,
		       const IA_DiscreteImage<T> &im_)
    :IA_OtherDI<T>(IA_PPCompClosureIPS(f_, im_.domain())),
     f(f_.clone_self()), im(im_) {
    }
    IA_PPCompClosureDI(const IA_PPCompClosureDI &other) 
    :IA_OtherDI<T>(IA_PPCompClosureIPS(*other.f, other.im.domain())),
     f(other.f->clone_self()), im(other.im) { }
    ~IA_PPCompClosureDI() { delete f; }

    T operator()(const IA_IntPoint&p) const {
	return im((*f)(p));
    }
    IA_OtherDI<T>* clone_self() const {
	return new IA_PPCompClosureDI<T>(*this);
    }
#if 0
    T *range(unsigned *n) const {
	return im.range(n);
    }
#endif
};


template <class T>
IA_DiscreteImage<T> compose
(const IA_DiscreteImage<T>& im,
 const IA_PPCompClosure<IA_IntPoint, IA_IntPoint> &f)
{
    return IA_PPCompClosureDI<T>(f, im);
}

//
//
//

template <class T>
IA_DiscreteImage<T> compose(T (*f)(T), const IA_DiscreteImage<T> &i)
{
    IA_DIVIter<T>	iter(i);
    const int size = i.card();
    T * const vals = new T[size];

    {
	T	*scan=vals;
	T	temp;
	while (iter(temp))
	    *(scan++) = f(temp);
    }

    return IA_DiscreteImage<T>(i.domain(), vals, size, 1);
}

template <class S,class T>
IA_DiscreteImage<S> compose(S (*f)(T), const IA_DiscreteImage<T> &i, const S*)
{
    IA_DIVIter<T>	iter(i);
    const int size = i.card();
    S * const vals = new S[size];

    {
	S	*scan=vals;
	T	temp;
	while (iter(temp))
	    *(scan++) = f(temp);
    }

    return IA_DiscreteImage<S>(i.domain(), vals, size, 1);
}

template <class R, class S, class T>
IA_DiscreteImage<R> compose(R (*f)(S), const IA_DiscreteImage<T> &i, const R*, const S*)
{
    IA_DIVIter<T>	iter(i);
    const int size = i.card();
    R * const vals = new R[size];

    {
	R	*scan=vals;
	T	temp;
	while (iter(temp))
	    *(scan++) = f(temp);
    }

    return IA_DiscreteImage<R>(i.domain(), vals, size, 1);
}

template <class Q, class R, class S, class T>
IA_DiscreteImage<Q> compose(R (*f)(S), const IA_DiscreteImage<T> &i, const Q*, const R*, const S*)
{
    IA_DIVIter<T>	iter(i);
    const int size = i.card();
    Q * const vals = new Q[size];

    {
	Q	*scan=vals;
	T	temp;
	while (iter(temp))
	    *(scan++) = f(temp);
    }

    return IA_DiscreteImage<Q>(i.domain(), vals, size, 1);
}
