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

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


//
// $Log: AuxImage.c,v $
// Revision 1.8.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.8.1.1  1994/12/28  15:38:25  thoth
// Image extension, and restriction are now friends.
//
// Revision 1.8  1994/09/16  14:57:59  thoth
// Compositions now live in their own include files.
//
// Revision 1.7  1994/08/22  15:20:29  thoth
// DOS-inspired name rework.
//
// Revision 1.6  1994/07/25  17:24:13  thoth
// Name sanitization
//
// Revision 1.5  1994/03/30  13:50:00  thoth
// zero_extension should be a template function.
//
// Revision 1.4  1994/03/14  15:48:33  thoth
// ValueSet has been replaced by container class Set.
//
// Revision 1.3  1994/02/26  16:01:09  thoth
// New 3-argument spatial composition functions.  Allow destination
// pointsets to be specified.
//
// Revision 1.2  1994/02/01  17:01:37  thoth
// Fix some screwy instantiation problems.
//
// Revision 1.1  1994/01/31  16:43:38  thoth
// Initial revision
//

#include "Point.h"

#include "Set.h"
#include "IntPS.h"
#include "DblPS.h"
#include "PSIter.h"
#include "ClosurePS.h"

#include "VectorI.h"
#include "ConstI.h"

#include "Image.h"
#include "ClosureI.h"
#include "ImageIter.h"
#include "image_errors.h"

template <class T>
void zero_extend_vec_scan(const IA_SetStructure &ss,
				 const T **src,
				 T **dest,
				 T zero)
{
    for (unsigned i=0; i<ss.nintervals(); i++) {
	IA_ss_interval	temp(ss.retrieve_interval(i));
	if (temp.substructure == IA_SetStructure::FIRST_ONLY) {
	    (*src) += temp.count;
	} else if (temp.substructure == IA_SetStructure::BOTH) {
	    for (unsigned j=0; j<temp.count; j++) {
		*((*dest)++) = *((*src)++);
	    }
	} else if (temp.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (unsigned j=0; j<temp.count; j++)
		*((*dest)++) = zero;
	} else {
	    for (unsigned j=0; j<temp.count; j++)
		zero_extend_vec_scan(temp.substructure, src, dest, zero);
	}
    }
}

template <class T>
void zero_extend_iter_scan(const IA_SetStructure &ss,
				  IA_IVIter<IA_Point<int>,T> *srciter,
				  T **dest,
				  T zero)
{
    for (unsigned i=0; i<ss.nintervals(); i++) {
	IA_ss_interval	temp(ss.retrieve_interval(i));
	if (temp.substructure == IA_SetStructure::FIRST_ONLY) {
	    T	blah;
	    for (unsigned j=0; j<temp.count; j++) 
		(*srciter)(blah);
	} else if (temp.substructure == IA_SetStructure::BOTH) {
	    T	blah;
	    for (unsigned j=0; j<temp.count; j++) {
		(*srciter)(blah);
		*((*dest)++) = blah;
	    }
	} else if (temp.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (unsigned j=0; j<temp.count; j++)
		*((*dest)++) = zero;
	} else {
	    for (unsigned j=0; j<temp.count; j++)
		zero_extend_iter_scan(temp.substructure, srciter, dest, zero);
	}
    }
}

template <class T>
void zero_extend(const IA_Image<IA_Point<int>,T> &srcimg,
		 const IA_Set<IA_Point<int> > &dest_ps, T *dest_data,
		 T zero)
{
    IA_SetStructure	ss;
    intersect_with_dualstruct(srcimg.domain(), dest_ps, &ss);

    if (srcimg.type() == IA_VectorI<IA_Point<int>,T>::s_type()) {
	T	*srcdata = ((IA_VectorI<IA_Point<int>,T>*)srcimg.bip)->vec;
	zero_extend_vec_scan(ss, &srcdata, &dest_data, zero);
    } else {
	IA_IVIter<IA_Point<int>,T>	iter(srcimg);
	zero_extend_iter_scan(ss, &iter, &dest_data, zero);
    }
}


//
//  The chi_contains function between an image and a valueset.
//

template <class T>
IA_Image<IA_Point<int>, IA_Bit>
vec_chi_contains(const IA_VectorI<IA_Point<int>, T> *lhs, const IA_Set<T> &rhs)
{
    const IA_Set<IA_Point<int> > domain = lhs->domain();
    const unsigned sz = domain.card();
    IA_Bit *const vec = new IA_Bit[sz];

    {
	IA_Bit *ds = vec;
	const T *ls = lhs->vec;
	for (unsigned i=0; i<sz; i++) {
	    *(ds++) = rhs.contains(*ls);
	    ls++;
	}
    }
    return IA_Image<IA_Point<int>, IA_Bit>(domain, vec, sz, 1);
}


// scalar-scalar operation is trivial


template <class T>
IA_Image<IA_Point<int>, IA_Bit>
iter_chi_contains(const IA_Image<IA_Point<int>, T> &lhs, const IA_Set<T> &rhs)
{
    const IA_Set<IA_Point<int> > domain = lhs.domain();
    const unsigned sz = domain.card();
    IA_Bit *const vec = new IA_Bit[sz];

    {
	IA_Bit *ds = vec;
	T ltmp;
	IA_IVIter<IA_Point<int>, T>	liter(lhs);
	for (unsigned i=0; i<sz; i++) {
	    liter(ltmp);
	    *(ds++) = rhs.contains(ltmp);
	}
    }
    return IA_Image<IA_Point<int>, IA_Bit>(domain, vec, sz, 1);
}


template <class T>
IA_Image<IA_Point<int>, IA_Bit>
chi_contains ( const IA_Image<IA_Point<int>, T> &lhs, const IA_Set<T> &rhs)
{
    if (lhs.type() == IA_VectorI<IA_Point<int>, T>::s_type()) {
	return vec_chi_contains((IA_VectorI<IA_Point<int>, T> *)lhs.bip,rhs);
    } else if (lhs.type() == IA_ConstI<IA_Point<int>, T>::s_type()) {
	return IA_Image<IA_Point<int>, IA_Bit> (lhs.domain(), rhs.contains(((IA_ConstI<IA_Point<int>, T> *)lhs.bip)->value));
    } else {
	return iter_chi_contains(lhs,rhs);
    }
}

//
//
//

// transpose of a boxy 2D image

template <class T>
IA_Image<IA_Point<int>, T>
vec_transpose(const IA_VectorI<IA_Point<int>, T> *img)
{
    const IA_Set<IA_Point<int> > domain = img->domain();
    const unsigned sz = domain.card();
    T *const vec = new T[sz];
    unsigned	width = domain.sup()[0] - domain.inf()[0] + 1;
    unsigned	height = domain.sup()[1] - domain.inf()[1] + 1;

    {
	const T *p = img->vec;
	for (unsigned x=0; x<width; x++) {
	    for (unsigned y=0; y<height; y++) {
		vec[y*width + x] = p[x*height + y];
	    }
	}
    }
    return IA_Image<IA_Point<int>, T>(transpose(domain), vec, sz, 1);
}

template <class T>
IA_Image<IA_Point<int>, T>
iter_transpose(const IA_Image<IA_Point<int>, T> &img)
{
    const IA_Set<IA_Point<int> > domain = img.domain();
    const unsigned sz = domain.card();
    T *const vec = new T[sz];
    unsigned	width = domain.sup()[0] - domain.inf()[0] + 1;
    unsigned	height = domain.sup()[1] - domain.inf()[1] + 1;

    {
	T ltmp;
	IA_IVIter<IA_Point<int>, T>	liter(img);
	for (unsigned x=0; x<width; x++) {
	    for (unsigned y=0; y<height; y++) {
		liter(ltmp);
		vec[y*width + x] = ltmp;
	    }
	}
    }
    return IA_Image<IA_Point<int>, T>(transpose(domain), vec, sz, 1);
}

template <class T> IA_Image<IA_Point<int>,T>
transpose(const IA_Image<IA_Point<int>, T> & img)
{
    if (!img.domain().boxy()
	|| img.domain().dim()!=2) {
	ia_throw(Image_BadTranspose_Exception(__FILE__, __LINE__));
	return IA_Image<IA_Point<int>,T>();
    }

    if (img.type() == IA_VectorI<IA_Point<int>, T>::s_type()) {
	return vec_transpose((IA_VectorI<IA_Point<int>, T> *)img.bip);

    } else if (img.type() == IA_ConstI<IA_Point<int>, T>::s_type()) {
	return IA_Image<IA_Point<int>, T>
	    (transpose(img.domain()),
	     ((IA_ConstI<IA_Point<int>, T> *)img.bip)->value);

    } else {
	return iter_transpose(img);
    }
}

//
//
//


#if 1 // NYI

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

template <class P, class T>
IA_Image<P,T>
extend(const IA_Image<P,T>& a, const IA_Image<P,T>& b)
//
// Returns the extension of image a to image b.
//
{
    if(!a.domain().extensive() || !b.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(a.domain(),
					   b.domain(), &ss);
    const unsigned	sz = newps.card();
    T	* const dest = new T[sz];
    T	*scan = dest;

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

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

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

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

#endif // NYI

//
//

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

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

template <class P, class T>
IA_Image<P, T>
value_restrict(const IA_Image<P,T> &i, const IA_Set<T>&vs)
{
    if (!i.domain().extensive()) {
	ia_throw(Pointset_RequireExtensive_Exception(__FILE__, __LINE__));
    }
    
    const	unsigned sz = i.domain().card();
    P *point_vec = new P[sz];
    T *value_vec = new T[sz];
    
    IA_IPIter<P,T> iter(i);
    P point;
    T value;
    int new_size=0;
    
    // Copy all points in i\'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>(i.domain().dim(), point_vec, new_size),
	 value_vec, new_size, 1);
    
    delete[] point_vec;
    
    return rval;
}

//
//
//

template <class Q, class P>
class IA_PPCompositionPS: public IA_ClosurePS<P> {
    // This class is used to define a pointset that contains
    // a point of type P if the pointset "base" contains f(p).
    // It is a Point-Point composition PointSet.
    Q	(*f)(const P&);
    IA_Set<Q>	base;
  public:
    IA_PPCompositionPS(Q (*f_)(const P&),
			const IA_Set<Q> &base_)
    :IA_ClosurePS<P>(base_.dim()), f(f_), base(base_) {
    }
    IA_ClosurePS<P> *clone_self() const {
	return new IA_PPCompositionPS(*this);
    }
    // given:
    // c = compose(ps, f);
    // then
    // c.contains(p)   iff ps.contains(f(p))
    int contains(const P &ip) const {
	// Uglier than Rosanne bug
	Q (*f2)(const P&) = f;
	return base.contains(f2(ip));
    }
};

template <class P, class Q, class T>
class IA_PPCompositionI: public IA_ClosureI<P,T> {
    //  This defines an image whose values are taken
    // from a source image after applying a function to the point.
    // It is a Point-Point composition image.
    Q	(*f)(const P&);
    IA_Image<Q,T>	im;
  public:
    IA_PPCompositionI( const IA_Image<Q,T> &im_,
		       Q (*f_)(const P&))
    :IA_ClosureI<P,T>(IA_PPCompositionPS<Q,P>(f_, im_.domain())),
     f(f_), im(im_) { }
    // given:
    // c = compose(i, f);
    // then
    // c(p) is defined by i(f(p))
    T operator()(const P&p) const {
	Q (*f2)(const P &) = f;
	return im(f2(p));
    }
    IA_ClosureI<P,T>* clone_self() const {
	return new IA_PPCompositionI<P,Q,T>(*this);
    }
#if 0
    T *range(unsigned *n) const {
	return im.range(n);
    }
#endif
};

template <class P, class T>
IA_Image<P,T> compose(const IA_Image<P,T> &i,
		      P (*f)(const P&))
{
    return IA_PPCompositionI<P,P,T>(i, f);
}

template <class P, class T>
IA_Image<P,T> compose(const IA_Image<P,T> &i,
			  P (*f)(const P&),
			  const IA_Set<P> &ps)
{
    IA_PSIter<P>	iter(ps);
    P	p;
    const unsigned	size=ps.card();
    T * const	vec = new T[size];

    {
	T	*dst = vec;
	while (iter(p)) {
	    *(dst++) = i(f(p));
	}
    }
    return IA_Image<P,T>(ps, vec, size, 1);
}

//
//
//

#include "Closure.h"

template <class Q, class P>
class IA_PPCompClosurePS: public IA_ClosurePS<P> {
    IA_Closure< P,Q> *f;
    IA_Set<Q>	base;
  public:
    IA_PPCompClosurePS( const IA_Closure< P,Q> &f_,
			const IA_Set<Q> &base_)
    :IA_ClosurePS<P>(base_.dim()), f(f_.clone_self()), base(base_) {
    }
    IA_PPCompClosurePS(const IA_PPCompClosurePS &other)
    :IA_ClosurePS<P>(other.base.dim()), f(other.f->clone_self()),
     base(other.base) {
	 // f has been cloned
    }
    ~IA_PPCompClosurePS() { delete f; }

    IA_ClosurePS<P> *clone_self() const {
	return new IA_PPCompClosurePS(*this);
    }
    int contains(const P &ip) const {
	return base.contains((*f)(ip));
    }
};

template <class P, class Q, class T>
class IA_PPCompClosureI: public IA_ClosureI<P,T> {
    IA_Closure<P,Q> *f;
    IA_Image<Q,T>	im;
  public:
    IA_PPCompClosureI( const IA_Closure< P,Q> &f_,
		       const IA_Image<Q,T> &im_)
    :IA_ClosureI<P,T>(IA_PPCompClosurePS<P,Q>(f_, im_.domain())),
     f(f_.clone_self()), im(im_) { }
    IA_PPCompClosureI(const IA_PPCompClosureI &other) 
    :IA_ClosureI<P,T>(IA_PPCompClosurePS<P,Q>(*other.f, other.im.domain())),
     f(other.f->clone_self()), im(other.im) { }
    ~IA_PPCompClosureI() { delete f; }

    T operator()(const P&p) const {
	return im((*f)(p));
    }
    IA_ClosureI<P,T>* clone_self() const {
	return new IA_PPCompClosureI<P,Q,T>(*this);
    }
#if 0
    T *range(unsigned *n) const {
	return im.range(n);
    }
#endif
};


template <class P, class T>
IA_Image<P,T> compose
(const IA_Image<P,T>& im,
 const IA_Closure<P,P> &f)
{
    return IA_PPCompClosureI<P,P,T>(f, im);
}

template <class P, class T>
IA_Image<P,T> compose(const IA_Image<P,T> &i,
			  const IA_Closure<P,P> &f,
			  const IA_Set<P> &ps)
{
    IA_PSIter<P>	iter(ps);
    P	p;
    const unsigned	size=ps.card();
    T * const	vec = new T[size];

    {
	T	*dst = vec;
	while (iter(p)) {
	    *(dst++) = i(f(p));
	}
    }
    return IA_Image<P,T>(ps, vec, size, 1);
}


//
//
//

template <class P, class T>
IA_Image<P,T> reflect(const IA_Image<P,T>&arg)
{
    return arg.bip->reflection();
}

//
//
//

#if 0
template < class P,class T>
ostream& operator <<(ostream& o, const IA_Image<P,T> &im)
{
    return o << *(im.bip);
}
#endif
