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

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


//
// $Log: AuxDDTempl.c,v $
// Revision 1.7.1.2  1995/01/13  19:36:30  thoth
// CoreImage has been re-merged with Image.  All special stuff should be
// friends of FBI now.
//
// Revision 1.7.1.1  1995/01/09  17:48:14  thoth
// Boxy pointsets constructor replaced with function.
//
// Revision 1.7  1994/12/28  22:13:17  ljr
// Bogus return types fixed.
//
// Revision 1.6  1994/09/16  14:57:11  thoth
// more DOS-inspired renaming.
//
// Revision 1.5  1994/08/22  15:20:29  thoth
// DOS-inspired name rework.
//
// Revision 1.4  1994/07/25  17:26:39  thoth
// Name sanitization
//
// Revision 1.3  1994/04/15  13:15:11  thoth
// Flesh out the rest of the generic convolutions.
//
// Revision 1.2  1994/03/30  13:49:19  thoth
// generic convolutions will now compile.
//
// Revision 1.1  1994/03/22  16:08:28  thoth
// Initial revision
//

#include "AuxDDTempl.h"
#include "FBI.h"

#include "IntPoint.h"
#include "IntPS.h"
#include "PSIter.h"
#include "ImageIter.h"
#include "InvDT.h"

//
// Generic backward convolution
//

template <class R, class S, class T>
IA_Image<IA_Point<int>,R>
backw_generic_inv_core(IA_Point<int> src_infimum,
		    IA_Point<int> src_width,
		    const T *src_data, // length is prod(src_width)
		    const IA_Image<IA_Point<int>,S> &templ,
		    IA_Set<IA_Point<int> > dest_ps,
		    R (*circle)(T,S),
		    R (*gamma)(R,R),
		    R gamma_zero,
		    const S*)
{
    const int	dimen = src_width.dim();

    IA_Set<IA_Point<int> >	templ_ps = templ.domain();
    int		templ_sz = templ_ps.card();
    S	*const templ_data = new S[templ_sz];
    int	*const templ_offsets = new int[templ_sz];

    {
	S	*d_scan = templ_data;
	int	*o_scan = templ_offsets;
	IA_IPIter<IA_Point<int>,S>	iter(templ);
	IA_Point<int>	ip;
	while (iter(ip, *d_scan)) {
	    *o_scan = ip[0];
	    for (unsigned i=1; i<dimen; i++) {
		*o_scan *= src_width[i];
		*o_scan += ip[i];
	    }
	    d_scan++;
	    o_scan++;
	}
    }

    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    R *const	dest_data = new R[dest_ps.card()];
    R *	valp = dest_data;
    while (iter(ip)) {
	int	offset= ip[0] - src_infimum[0];
	unsigned i;
	for (i=1; i<dimen; i++) {
	    offset *= src_width[i];
	    offset += ip[i] - src_infimum[i];
	}
	const T *const base = src_data + offset;

	*valp = gamma_zero;
	for (i=0; i<templ_sz; i++) {
	    *valp = gamma( circle(base[templ_offsets[i]],
				  templ_data[i]), *valp);
	}

	valp++;
    }
    delete[] templ_offsets;
    delete[] templ_data;

    return IA_Image<IA_Point<int>,R>(dest_ps, dest_data, dest_ps.card(), 1);
}

template <class R, class S, class T>
IA_Image<IA_Point<int>, R>
backw_generic_inv(const IA_Image<IA_Point<int>, T> &img,
		  const IA_Image<IA_Point<int>, S> &invtempl,
		  const IA_Set<IA_Point<int> > &dest_ps,
		  R (*circle)(T,S),
		  R (*gamma)(R,R),
		  R gamma_zero,
		  T ring_zero,
		  const S*)
{
    IA_Point<int>	inf_ = dest_ps.inf()+invtempl.domain().inf();
    IA_Point<int>	sup_ = dest_ps.sup()+invtempl.domain().sup();
    IA_Set<IA_Point<int> >	src_ps = IA_boxy_pset(inf_,sup_);

    T	*const src_data = new T[src_ps.card()];

    zero_extend(IA_Image<IA_Point<int>,T>(img), IA_Set<IA_Point<int> >(src_ps),
		(T*)src_data, T(ring_zero));
    IA_Image<IA_Point<int>,R>	rval =
	backw_generic_inv_core(inf_, sup_-inf_ + 1, src_data,
			       invtempl, dest_ps,
			       circle, gamma, gamma_zero, (S*)0);

    delete[] src_data;

    return rval;
}

template <class R, class S, class T>
IA_Image<IA_Point<int>, R>
generic_product (const IA_Image<IA_Point<int>, T> &img,
		 const IA_DDTemplate< IA_Image<IA_Point<int>, S> > &templ,
		 const IA_Set<IA_Point<int> > &dest_ps,
		 R (*circle)(T,S),
		 R (*gamma)(R,R),
		 R gamma_zero,
		 T ring_zero,
		 const S*)
{

    if (templ.is_a(IA_InvariantDT<IA_Image<IA_Point<int>, S> >::s_type())) {
	return backw_generic_inv
	    (img, ((IA_InvariantDT<IA_Image<IA_Point<int>, S> >*)IA_FBI<IA_Point<int>, IA_Point<int>, S, S, S>::extract_baseptr(templ))->value,
	     dest_ps, circle, gamma, gamma_zero, ring_zero,
	     (S*)0);
	//    } else if (templ.type() == IA_Lazy_REDUCT_DT<IA_Image<IA_Point<int>, S> >::s_type()) {
    } else {
	R *const	dest_data = new R[dest_ps.card()];
	R	*valp = dest_data;
	IA_PSIter<IA_Point<int> >	dest_iter(dest_ps);
	IA_Point<int>	base_ip;
	while (dest_iter(base_ip)) {
	    IA_Image<IA_Point<int>,S>	tv = templ(base_ip);

	    IA_IPIter<IA_Point<int>,S>	templ_iter(tv);
	    IA_Point<int>	templ_ip;
	    S	templ_val;

	    *valp = gamma_zero;

	    while ( templ_iter(templ_ip, templ_val) ) {
		const IA_Point<int>	ip = templ_ip; //+base_ip;
		if (! img.domain().contains(ip))
		    continue;
		*valp = gamma((*circle)(img(ip), templ_val), *valp);
	    }
	    valp++;
	}
	return IA_Image<IA_Point<int>,R>(dest_ps, dest_data,
				       dest_ps.card(), 1);
    }

}

//
// Generic backward convolution with big Gamma
//

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>,Q>
backw_Generic_inv_core(IA_Point<int> src_infimum,
		       IA_Point<int> src_width,
		       const T *src_data, // length is prod(src_width)
		       const IA_Image<IA_Point<int>,S> &templ,
		       IA_Set<IA_Point<int> > dest_ps,
		       R (*circle)(T,S),
		       Q (*Gamma)(R*,unsigned),
		       const Q*,
		       const R*,
		       const S*)
{
    const int	dimen = src_width.dim();

    IA_Set<IA_Point<int> >	templ_ps = templ.domain();
    int		templ_sz = templ_ps.card();
    S	*const templ_data = new S[templ_sz];
    int	*const templ_offsets = new int[templ_sz];

    {
	S	*d_scan = templ_data;
	int	*o_scan = templ_offsets;
	IA_IPIter<IA_Point<int>,S>	iter(templ);
	IA_Point<int>	ip;
	while (iter(ip, *d_scan)) {
	    *o_scan = ip[0];
	    for (unsigned i=1; i<dimen; i++) {
		*o_scan *= src_width[i];
		*o_scan += ip[i];
	    }
	    d_scan++;
	    o_scan++;
	}
    }

    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    Q *const	dest_data = new R[dest_ps.card()];
    Q *	valp = dest_data;
    R	*a = new R[templ_sz];

    while (iter(ip)) {
	int	offset= ip[0] - src_infimum[0];
	unsigned i;
	for (i=1; i<dimen; i++) {
	    offset *= src_width[i];
	    offset += ip[i] - src_infimum[i];
	}
	const T *const base = src_data + offset;

	for (i=0; i<templ_sz; i++) {
	    a[i] = circle(base[templ_offsets[i]], templ_data[i]);
	}

	*(valp++) = Gamma(a, templ_sz);
    }
    delete[] a;

    delete[] templ_offsets;
    delete[] templ_data;

    return IA_Image<IA_Point<int>,Q>(dest_ps, dest_data, dest_ps.card(), 1);
}

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>, Q>
backw_Generic_inv(const IA_Image<IA_Point<int>, T> &img,
		  const IA_Image<IA_Point<int>, S> &invtempl,
		  const IA_Set<IA_Point<int> > &dest_ps,
		  R (*circle)(T,S),
		  Q (*Gamma)(R*,unsigned),
		  T ring_zero,
		  const Q*,
		  const R*,
		  const S*)
{
    IA_Point<int>	inf_ = dest_ps.inf()+invtempl.domain().inf();
    IA_Point<int>	sup_ = dest_ps.sup()+invtempl.domain().sup();
    IA_Set<IA_Point<int> >	src_ps = IA_boxy_pset(inf_,sup_);

    T	*const src_data = new T[src_ps.card()];

    zero_extend(IA_Image<IA_Point<int>,T>(img), IA_Set<IA_Point<int> >(src_ps),
		(T*)src_data, T(ring_zero));
    IA_Image<IA_Point<int>,Q>	rval =
	backw_Generic_inv_core(inf_, sup_-inf_ + 1, src_data,
			       invtempl, dest_ps,
			       circle, Gamma, (Q*)0, (R*)0, (S*)0);

    delete[] src_data;

    return rval;
}

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>, Q>
generic_product (const IA_Image<IA_Point<int>, T> &img,
		 const IA_DDTemplate< IA_Image<IA_Point<int>, S> > &templ,
		 const IA_Set<IA_Point<int> > &dest_ps,
		 R (*circle)(T,S),
		 Q (*Gamma)(R*,unsigned),
		 T ring_zero,
		 const Q*,
		 const R*,
		 const S*)
{

    if (templ.is_a(IA_InvariantDT<IA_Image<IA_Point<int>, S> >::s_type())) {
	return backw_Generic_inv
	    (img, (IA_Image<IA_Point<int>,S>&)((IA_InvariantDT<IA_Image<IA_Point<int>, S> >*)IA_FBI<IA_Point<int>, IA_Point<int>, S, S, S>::extract_baseptr(templ))->value,
	     dest_ps, (R(*)(T,S))circle,
	     (Q(*)(R*,unsigned))Gamma, T(ring_zero),
	     (const Q*)0, (const R*)0, (const S*)0);
	//    } else if (templ.type() == IA_Lazy_REDUCT_DT<IA_Image<IA_Point<int>, S> >::s_type()) {
    } else {
	Q *const	dest_data = new Q[dest_ps.card()];
	Q	*valp = dest_data;
	IA_PSIter<IA_Point<int> >	dest_iter(dest_ps);
	IA_Point<int>	base_ip;
	while (dest_iter(base_ip)) {
	    IA_Image<IA_Point<int>,S>	tv = templ(base_ip);

	    IA_IPIter<IA_Point<int>,S>	templ_iter(tv);
	    IA_Point<int>	templ_ip;
	    S	templ_val;

	    R	*a = new R[tv.card()];
	    int	idx=0;

	    while ( templ_iter(templ_ip, templ_val) ) {
		const IA_Point<int>	ip = templ_ip; //+base_ip;
		if (! img.domain().contains(ip))
		    continue;
		a[idx++] = (*circle)(img(ip), templ_val);
	    }
	    *(valp++) = Gamma(a, idx);
	    delete[] a;		// deleting every time is inefficient, rewrite
	}
	return IA_Image<IA_Point<int>,Q>(dest_ps, dest_data,
				       dest_ps.card(), 1);
    }

}

//
// Generic backward convolution with big Gamma
//

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>,Q>
backw_Generic_inv_core3(IA_Point<int> src_infimum,
		       IA_Point<int> src_width,
		       const T *src_data, // length is prod(src_width)
		       IA_Set<IA_Point<int> > src_ps,
		       const IA_Image<IA_Point<int>,S> &templ,
		       IA_Set<IA_Point<int> > dest_ps,
		       R (*circle)(T,S),
		       Q (*Gamma)(R*,unsigned),
		       const Q*,
		       const R*,
		       const S*)
{
    const int	dimen = src_width.dim();

    IA_Set<IA_Point<int> >	templ_ps = templ.domain();
    int		templ_sz = templ_ps.card();
    S	*const templ_data = new S[templ_sz];
    int	*const templ_offsets = new int[templ_sz];
    IA_Point<int>	*const templ_points = new IA_Point<int>[templ_sz];

    {
	IA_IPIter<IA_Point<int>,S>	iter(templ);

	for (int i=0; iter(templ_points[i], templ_data[i]); i++) {
	    templ_offsets[i] = templ_points[i][0];
	    for (unsigned j=1; j<dimen; j++) {
		templ_offsets[i] *= src_width[j];
		templ_offsets[i] += templ_points[i][j];
	    }
	}
    }


    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    Q *const	dest_data = new R[dest_ps.card()];
    Q *	valp = dest_data;
    R	*a = new R[templ_sz];

    while (iter(ip)) {
	int	offset= ip[0] - src_infimum[0];
	unsigned i;
	for (i=1; i<dimen; i++) {
	    offset *= src_width[i];
	    offset += ip[i] - src_infimum[i];
	}
	const T *const base = src_data + offset;

	unsigned count=0;
	for (i=0; i<templ_sz; i++) {
	    if (src_ps.contains(ip + templ_points[i]))
		a[count++] = circle(base[templ_offsets[i]], templ_data[i]);
	}

	*(valp++) = Gamma(a, count);
    }
    delete[] a;

    delete[] templ_points;
    delete[] templ_offsets;
    delete[] templ_data;

    return IA_Image<IA_Point<int>,Q>(dest_ps, dest_data, dest_ps.card(), 1);
}

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>, Q>
backw_Generic_inv3(const IA_Image<IA_Point<int>, T> &img,
		  const IA_Image<IA_Point<int>, S> &invtempl,
		  const IA_Set<IA_Point<int> > &dest_ps,
		  R (*circle)(T,S),
		  Q (*Gamma)(R*,unsigned),
		  const Q*,
		  const R*,
		  const S*)
{
    IA_Point<int>	inf_ = dest_ps.inf()+invtempl.domain().inf();
    IA_Point<int>	sup_ = dest_ps.sup()+invtempl.domain().sup();
    IA_Set<IA_Point<int> >	src_ps = IA_boxy_pset(inf_,sup_);

    T	*const src_data = new T[src_ps.card()];

    T	ring_zero;
    zero_extend(IA_Image<IA_Point<int>,T>(img), IA_Set<IA_Point<int> >(src_ps),
		(T*)src_data, ring_zero);
    IA_Image<IA_Point<int>,Q>	rval =
	backw_Generic_inv_core3
	    (inf_, sup_-inf_ + 1, src_data, img.domain(),
	     invtempl, dest_ps,
	     circle, Gamma, (Q*)0, (R*)0, (S*)0);

    delete[] src_data;

    return rval;
}

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>, Q>
generic_product (const IA_Image<IA_Point<int>, T> &img,
		 const IA_DDTemplate< IA_Image<IA_Point<int>, S> > &templ,
		 const IA_Set<IA_Point<int> > &dest_ps,
		 R (*circle)(T,S),
		 Q (*Gamma)(R*,unsigned),
		 const Q*,
		 const R*,
		 const S*)
{

    if (templ.is_a(IA_InvariantDT<IA_Image<IA_Point<int>, S> >::s_type())) {
	return backw_Generic_inv3
	    (img, (IA_Image<IA_Point<int>,S>&)((IA_InvariantDT<IA_Image<IA_Point<int>, S> >*)IA_FBI<IA_Point<int>, IA_Point<int>, S, S, S>::extract_baseptr(templ))->value,
	     dest_ps, (R(*)(T,S))circle,
	     (Q(*)(R*,unsigned))Gamma,
	     (const Q*)0, (const R*)0, (const S*)0);
	//    } else if (templ.type() == IA_Lazy_REDUCT_DT<IA_Image<IA_Point<int>, S> >::s_type()) {
    } else {
	Q *const	dest_data = new Q[dest_ps.card()];
	Q	*valp = dest_data;
	IA_PSIter<IA_Point<int> >	dest_iter(dest_ps);
	IA_Point<int>	base_ip;
	while (dest_iter(base_ip)) {
	    IA_Image<IA_Point<int>,S>	tv = templ(base_ip);

	    IA_IPIter<IA_Point<int>,S>	templ_iter(tv);
	    IA_Point<int>	templ_ip;
	    S	templ_val;

	    R	*a = new R[tv.card()];
	    int	idx=0;

	    while ( templ_iter(templ_ip, templ_val) ) {
		const IA_Point<int>	ip = templ_ip; //+base_ip;
		if (! img.domain().contains(ip))
		    continue;
		a[idx++] = (*circle)(img(ip), templ_val);
	    }
	    *(valp++) = Gamma(a, idx);
	    delete[] a;		// deleting every time is inefficient, rewrite
	}
	return IA_Image<IA_Point<int>,Q>(dest_ps, dest_data,
				       dest_ps.card(), 1);
    }

}

//
// Generic forward convolution
//

template <class R, class S, class T>
IA_Image<IA_Point<int>,R>
forw_generic_inv_core(IA_Point<int> src_infimum,
		    IA_Point<int> src_width,
		    const T *src_data, // length is prod(src_width)
		    const IA_Image<IA_Point<int>,S> &templ,
		    IA_Set<IA_Point<int> > dest_ps,
		    R (*circle)(T,S),
		    R (*gamma)(R,R),
		    R gamma_zero,
		    const S*)
{
    const int	dimen = src_width.dim();

    IA_Set<IA_Point<int> >	templ_ps = templ.domain();
    int		templ_sz = templ_ps.card();
    S	*const templ_data = new S[templ_sz];
    int	*const templ_offsets = new int[templ_sz];

    {
	S	*d_scan = templ_data;
	int	*o_scan = templ_offsets;
	IA_IPIter<IA_Point<int>,S>	iter(templ);
	IA_Point<int>	ip;
	while (iter(ip, *d_scan)) {
	    *o_scan = -ip[0];
	    for (unsigned i=1; i<dimen; i++) {
		*o_scan *= src_width[i];
		*o_scan -= ip[i];
	    }
	    d_scan++;
	    o_scan++;
	}
    }

    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    R *const	dest_data = new R[dest_ps.card()];
    R *	valp = dest_data;
    while (iter(ip)) {
	int	offset= ip[0] - src_infimum[0];
	unsigned i;
	for (i=1; i<dimen; i++) {
	    offset *= src_width[i];
	    offset += ip[i] - src_infimum[i];
	}
	const T *const base = src_data + offset;

	*valp = gamma_zero;
	for (i=0; i<templ_sz; i++) {
	    *valp = gamma( circle(base[templ_offsets[i]],
				  templ_data[i]), *valp);
	}

	valp++;
    }
    delete[] templ_offsets;
    delete[] templ_data;

    return IA_Image<IA_Point<int>,R>(dest_ps, dest_data, dest_ps.card(), 1);
}

template <class R, class S, class T>
IA_Image<IA_Point<int>, R>
forw_generic_inv(const IA_Image<IA_Point<int>, T> &img,
		  const IA_Image<IA_Point<int>, S> &invtempl,
		  const IA_Set<IA_Point<int> > &dest_ps,
		  R (*circle)(T,S),
		  R (*gamma)(R,R),
		  R gamma_zero,
		  T ring_zero,
		  const S*)
{
    IA_Point<int>	inf_ = dest_ps.inf()-invtempl.domain().sup();
    IA_Point<int>	sup_ = dest_ps.sup()-invtempl.domain().inf();
    IA_Set<IA_Point<int> >	src_ps = IA_boxy_pset(inf_,sup_);

    T	*const src_data = new T[src_ps.card()];

    zero_extend(IA_Image<IA_Point<int>,T>(img), IA_Set<IA_Point<int> >(src_ps),
		(T*)src_data, T(ring_zero));
    IA_Image<IA_Point<int>,R>	rval =
	forw_generic_inv_core(inf_, sup_-inf_ + 1, src_data,
			       invtempl, dest_ps,
			       circle, gamma, gamma_zero, (S*)0);

    delete[] src_data;

    return rval;
}

template <class R, class S, class T>
IA_Image<IA_Point<int>, R>
generic_product (const IA_DDTemplate< IA_Image<IA_Point<int>, S> > &templ,
		 const IA_Image<IA_Point<int>, T> &img,
		 const IA_Set<IA_Point<int> > &dest_ps,
		 R (*circle)(T,S),
		 R (*gamma)(R,R),
		 R gamma_zero,
		 T ring_zero,
		 const S*)
{

    if (templ.is_a(IA_InvariantDT<IA_Image<IA_Point<int>, S> >::s_type())) {
	return forw_generic_inv
	    (img, ((IA_InvariantDT<IA_Image<IA_Point<int>, S> >*)IA_FBI<IA_Point<int>, IA_Point<int>, S, S, S>::extract_baseptr(templ))->value,
	     dest_ps, circle, gamma, gamma_zero, ring_zero,
	     (S*)0);
    } else {
	unsigned	size = dest_ps.card();
	R *const	dest_data = new R[size];

	for (int i=0; i<size; i++) {
	    dest_data[i] = gamma_zero;
	}

	IA_PSIter<IA_Point<int> >	src_iter(img.domain());
	IA_Point<int>	base_ip;
	while (src_iter(base_ip)) {
	    IA_Image<IA_Point<int>,S>	tv = templ(base_ip);

	    IA_IPIter<IA_Point<int>,S>	templ_iter(tv);
	    IA_Point<int>	templ_ip;
	    S	templ_val;
	    T	ival = img(base_ip);

	    while ( templ_iter(templ_ip, templ_val) ) {
		const IA_Point<int>	ip = templ_ip; //+base_ip;
		if (! dest_ps.contains(templ_ip))
		    continue;
		unsigned	offset = dest_ps.index(templ_ip);
		R	*valp = dest_data + offset;
		*valp = gamma(circle(templ_val, ival),*valp);
	    }
	}
	return IA_Image<IA_Point<int>,R>(dest_ps, dest_data,
				       size, 1);
    }

}

//
// Generic forward convolution with big Gamma
//

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>,Q>
forw_Generic_inv_core(IA_Point<int> src_infimum,
		       IA_Point<int> src_width,
		       const T *src_data, // length is prod(src_width)
		       const IA_Image<IA_Point<int>,S> &templ,
		       IA_Set<IA_Point<int> > dest_ps,
		       R (*circle)(T,S),
		       Q (*Gamma)(R*,unsigned),
		       const Q*,
		       const R*,
		       const S*)
{
    const int	dimen = src_width.dim();

    IA_Set<IA_Point<int> >	templ_ps = templ.domain();
    int		templ_sz = templ_ps.card();
    S	*const templ_data = new S[templ_sz];
    int	*const templ_offsets = new int[templ_sz];

    {
	S	*d_scan = templ_data;
	int	*o_scan = templ_offsets;
	IA_IPIter<IA_Point<int>,S>	iter(templ);
	IA_Point<int>	ip;
	while (iter(ip, *d_scan)) {
	    *o_scan = -ip[0];
	    for (unsigned i=1; i<dimen; i++) {
		*o_scan *= src_width[i];
		*o_scan -= ip[i];
	    }
	    d_scan++;
	    o_scan++;
	}
    }

    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    Q *const	dest_data = new R[dest_ps.card()];
    Q *	valp = dest_data;
    R	*a = new R[templ_sz];

    while (iter(ip)) {
	int	offset= ip[0] - src_infimum[0];
	unsigned i;
	for (i=1; i<dimen; i++) {
	    offset *= src_width[i];
	    offset += ip[i] - src_infimum[i];
	}
	const T *const base = src_data + offset;

	for (i=0; i<templ_sz; i++) {
	    a[i] = circle(base[templ_offsets[i]], templ_data[i]);
	}

	*(valp++) = Gamma(a, templ_sz);
    }
    delete[] a;

    delete[] templ_offsets;
    delete[] templ_data;

    return IA_Image<IA_Point<int>,Q>(dest_ps, dest_data, dest_ps.card(), 1);
}

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>, Q>
forw_Generic_inv(const IA_Image<IA_Point<int>, T> &img,
		  const IA_Image<IA_Point<int>, S> &invtempl,
		  const IA_Set<IA_Point<int> > &dest_ps,
		  R (*circle)(T,S),
		  Q (*Gamma)(R*,unsigned),
		  T ring_zero,
		  const Q*,
		  const R*,
		  const S*)
{
    IA_Point<int>	inf_ = dest_ps.inf()-invtempl.domain().sup();
    IA_Point<int>	sup_ = dest_ps.sup()-invtempl.domain().inf();
    IA_Set<IA_Point<int> >	src_ps = IA_boxy_pset(inf_,sup_);

    T	*const src_data = new T[src_ps.card()];

    zero_extend(IA_Image<IA_Point<int>,T>(img), IA_Set<IA_Point<int> >(src_ps),
		(T*)src_data, T(ring_zero));
    IA_Image<IA_Point<int>,Q>	rval =
	forw_Generic_inv_core(inf_, sup_-inf_ + 1, src_data,
			       invtempl, dest_ps,
			       circle, Gamma, (Q*)0, (R*)0, (S*)0);

    delete[] src_data;

    return rval;
}

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>, Q>
generic_product (const IA_DDTemplate< IA_Image<IA_Point<int>, S> > &templ,
		 const IA_Image<IA_Point<int>, T> &img,
		 const IA_Set<IA_Point<int> > &dest_ps,
		 R (*circle)(T,S),
		 Q (*Gamma)(R*,unsigned),
		 T ring_zero,
		 const Q*,
		 const R*,
		 const S*)
{

    if (templ.is_a(IA_InvariantDT<IA_Image<IA_Point<int>, S> >::s_type())) {
	return forw_Generic_inv
	    (img, (IA_Image<IA_Point<int>,S>&)((IA_InvariantDT<IA_Image<IA_Point<int>, S> >*)IA_FBI<IA_Point<int>, IA_Point<int>, S, S, S>::extract_baseptr(templ))->value,
	     dest_ps, (R(*)(T,S))circle,
	     (Q(*)(R*,unsigned))Gamma, T(ring_zero),
	     (const Q*)0, (const R*)0, (const S*)0);
	//    } else if (templ.type() == IA_Lazy_REDUCT_DT<IA_Image<IA_Point<int>, S> >::s_type()) {
    } else {
	IA::not_yet_implemented(__FILE__, __LINE__);
    }

}

//
// Generic forward convolution with big Gamma
//

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>,Q>
forw_Generic_inv_core3(IA_Point<int> src_infimum,
		       IA_Point<int> src_width,
		       const T *src_data, // length is prod(src_width)
		       IA_Set<IA_Point<int> > src_ps,
		       const IA_Image<IA_Point<int>,S> &templ,
		       IA_Set<IA_Point<int> > dest_ps,
		       R (*circle)(T,S),
		       Q (*Gamma)(R*,unsigned),
		       const Q*,
		       const R*,
		       const S*)
{
    const int	dimen = src_width.dim();

    IA_Set<IA_Point<int> >	templ_ps = templ.domain();
    int		templ_sz = templ_ps.card();
    S	*const templ_data = new S[templ_sz];
    int	*const templ_offsets = new int[templ_sz];
    IA_Point<int>	*const templ_points = new IA_Point<int>[templ_sz];

    {
	int	i;
	IA_IPIter<IA_Point<int>,S>	iter(templ);

	for (i=0; iter(templ_points[i], templ_data[i]); i++) {
	    templ_offsets[i] = -templ_points[i][0];
	    for (unsigned j=1; j<dimen; j++) {
		templ_offsets[i] *= src_width[j];
		templ_offsets[i] -= templ_points[i][j];
	    }
	}
    }

    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    Q *const	dest_data = new R[dest_ps.card()];
    Q *	valp = dest_data;
    R	*a = new R[templ_sz];

    while (iter(ip)) {
	int	offset= ip[0] - src_infimum[0];
	unsigned i;
	for (i=1; i<dimen; i++) {
	    offset *= src_width[i];
	    offset += ip[i] - src_infimum[i];
	}
	const T *const base = src_data + offset;

	unsigned count=0;
	for (i=0; i<templ_sz; i++) {
	    if (src_ps.contains(ip - templ_points[i]))
		a[count++] = circle(base[templ_offsets[i]], templ_data[i]);
	}

	*(valp++) = Gamma(a, count);
    }
    delete[] a;

    delete[] templ_points;
    delete[] templ_offsets;
    delete[] templ_data;

    return IA_Image<IA_Point<int>,Q>(dest_ps, dest_data, dest_ps.card(), 1);
}

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>, Q>
forw_Generic_inv3(const IA_Image<IA_Point<int>, T> &img,
		  const IA_Image<IA_Point<int>, S> &invtempl,
		  const IA_Set<IA_Point<int> > &dest_ps,
		  R (*circle)(T,S),
		  Q (*Gamma)(R*,unsigned),
		  const Q*,
		  const R*,
		  const S*)
{
    IA_Point<int>	inf_ = dest_ps.inf()-invtempl.domain().sup();
    IA_Point<int>	sup_ = dest_ps.sup()-invtempl.domain().inf();
    IA_Set<IA_Point<int> >	src_ps = IA_boxy_pset(inf_,sup_);

    T	*const src_data = new T[src_ps.card()];

    T	ring_zero;
    zero_extend(IA_Image<IA_Point<int>,T>(img), IA_Set<IA_Point<int> >(src_ps),
		(T*)src_data, ring_zero);
    IA_Image<IA_Point<int>,Q>	rval =
	forw_Generic_inv_core3(inf_, sup_-inf_ + 1, src_data,
			       img.domain(), invtempl, dest_ps,
			       circle, Gamma, (Q*)0, (R*)0, (S*)0);

    delete[] src_data;

    return rval;
}

template <class Q, class R, class S, class T>
IA_Image<IA_Point<int>, Q>
generic_product (const IA_DDTemplate< IA_Image<IA_Point<int>, S> > &templ,
		 const IA_Image<IA_Point<int>, T> &img,
		 const IA_Set<IA_Point<int> > &dest_ps,
		 R (*circle)(T,S),
		 Q (*Gamma)(R*,unsigned),
		 const Q*,
		 const R*,
		 const S*)
{

    if (templ.is_a(IA_InvariantDT<IA_Image<IA_Point<int>, S> >::s_type())) {
	return forw_Generic_inv3
	    (img, (IA_Image<IA_Point<int>,S>&)((IA_InvariantDT<IA_Image<IA_Point<int>, S> >*)IA_FBI<IA_Point<int>, IA_Point<int>, S, S, S>::extract_baseptr(templ))->value,
	     dest_ps, (R(*)(T,S))circle,
	     (Q(*)(R*,unsigned))Gamma,
	     (const Q*)0, (const R*)0, (const S*)0);
	//    } else if (templ.type() == IA_Lazy_REDUCT_DT<IA_Image<IA_Point<int>, S> >::s_type()) {
    } else {
	IA::not_yet_implemented(__FILE__, __LINE__);
    }

}
