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

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


//
// $Log: Nbh.c,v $
// Revision 1.9.1.3  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.9.1.2  1995/01/09  18:08:04  thoth
// Boxy pointset constructor replaced by function.
//
// Revision 1.9.1.1  1994/12/28  18:22:11  thoth
// New Closure and Function Nbh constructors.
//
// Revision 1.9  1994/12/22  17:37:10  fjsoria
// removed static data memeber template for double nbh.
//
// Revision 1.8  1994/12/22  16:27:20  fjsoria
// static data member template problem fixed for g++
//
// Revision 1.7  1994/09/16  14:57:11  thoth
// more DOS-inspired renaming.
//
// Revision 1.6  1994/08/22  15:20:29  thoth
// DOS-inspired name rework.
//
// Revision 1.5  1994/07/25  17:26:39  thoth
// Name sanitization
//
// Revision 1.4  1994/05/08  19:40:25  thoth
// New automatic neighborhood decomposition capability.
//
// Revision 1.3  1994/04/15  13:17:10  thoth
// Finish up the list of generic reductions.
//
// Revision 1.2  1994/03/30  13:56:46  thoth
// generic neighborhood reductions are now available.
//
// Revision 1.1  1994/01/31  16:43:38  thoth
// Initial revision
//

#include	<values.h>

#include "Image.h"
#include "ImageIter.h"
#include "VectorI.h"
#include "PSIter.h"

#include "Nbh.h"
#include "InvNbh.h"
#include "ClosureNbh.h"
#include "FuncNbh.h"

#include "image_errors.h"

template <class P, class Q>
void IA_Neighborhood<P,Q>::set_and_reference(IA_BaseNbh<P,Q> *p)
{
  p->incr_ref();
  bnbh = p;
}

template <class P, class Q>
void IA_Neighborhood<P,Q>::disassociate()
{
  if (bnbh->decr_ref() <= 0)
    delete bnbh;
  bnbh = 0;
}

#ifndef IA_NO_AUTOINST_OF_STATIC_DATA
template <class P, class Q>
int IA_Neighborhood<P,Q>::auto_decomposition=1;
#else
int IA_Neighborhood<IA_Point<int>,IA_Point<int> >::auto_decomposition=1;
// double nieghborhoods not implemented yet
//int IA_Neighborhood<IA_Point<double>,IA_Point<double> >::auto_decomposition=1;
#endif


template <class P, class Q>
IA_Neighborhood<P,Q>::IA_Neighborhood(const IA_Neighborhood<P,Q> &arg)
{
    set_and_reference(arg.bnbh);
}


template <class P, class Q>
IA_Neighborhood<P,Q>::IA_Neighborhood(unsigned dim, const IA_Set<P> &ps)
{
    set_and_reference
      (new IA_InvariantNbh<P,Q>(IA_Set<Q>::WhiteHole(dim), ps));
}


template <class P, class Q>
IA_Neighborhood<P,Q>::IA_Neighborhood(const IA_ClosureNbh<P,Q> &cl)
{
    set_and_reference (cl.clone_self());
}

template <class P, class Q>
IA_Neighborhood<P,Q>::IA_Neighborhood(const IA_Set<Q> &s, IA_Set<P>(*f)(Q))
{
    set_and_reference (new IA_FunctionNbh<P,Q>(s, f));
}

template <class P, class Q>
IA_Neighborhood<P,Q>::IA_Neighborhood(const IA_Set<Q> &s, IA_Set<P>(*f)(const Q&))
{
    set_and_reference (new IA_FunctionRefNbh<P,Q>(s, f));
}

template <class P, class Q>
IA_Neighborhood<P,Q>::~IA_Neighborhood()
{
  disassociate();
}

template <class P, class Q>
IA_Neighborhood<P,Q>& IA_Neighborhood<P,Q>::operator=
	(const IA_Neighborhood<P,Q> &arg)
{
  IA_BaseNbh<P,Q>	*temp = arg.bnbh;
  temp->incr_ref();

  disassociate();

  bnbh = temp;

  return *this;
}

#include "NbhOps.c"

//
// Generic backward neighborhood reduction
//

template <class T>
IA_Image<IA_Point<int>, T> backw_nbh_inv_core1
(IA_Point<int> src_infimum,
 IA_Point<int> src_width,
 const T* src_data,
 const IA_Set<IA_Point<int> > &nbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 T (*gamma)(T, T),
 T gamma_zero)
{
    const int	dimen = src_width.dim();

    int		nbh_sz = nbh.card();
    int	*const nbh_offsets = new int[nbh_sz];

    {
	int	*o_scan = nbh_offsets;
	IA_PSIter<IA_Point<int> >	iter(nbh);
	IA_Point<int>	ip;
	while (iter(ip)) {
	    *o_scan = ip[0];
	    for (unsigned i=1; i<dimen; i++) {
		*o_scan *= src_width[i];
		*o_scan += ip[i];
	    }
	    o_scan++;
	}
    }

    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    T *const	dest_data = new T[dest_ps.card()];
    T *	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<nbh_sz; i++) {
	    *valp = gamma(base[nbh_offsets[i]], *valp);
	}
	valp++;
    }
    delete[] nbh_offsets;

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


template <class T>
IA_Image<IA_Point<int>, T> backw_nbh_inv1
(const IA_Image<IA_Point<int>, T> &img,
 const IA_Set<IA_Point<int> > &invnbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 T (*gamma)(T, T),
 T gamma_zero)
{
    IA_Point<int>	inf_ = dest_ps.inf()+invnbh.inf();
    IA_Point<int>	sup_ = dest_ps.sup()+invnbh.sup();
    IA_Set<IA_Point<int> >	src_ps = IA_boxy_pset(inf_,sup_);

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

    zero_extend(img, src_ps, (T*)src_data, gamma_zero);
    IA_Image<IA_Point<int>,T>	rval =
	backw_nbh_inv_core1(inf_, sup_-inf_ + 1, src_data, invnbh, dest_ps,
			    gamma, gamma_zero);

    delete[] src_data;

    return rval;
}

template <class T>
IA_Image<IA_Point<int>,T> neighborhood_reduction
(const IA_Image<IA_Point<int>, T> &img,
 const IA_Neighborhood<IA_Point<int>,IA_Point<int> > &nbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 T (*gamma)(T, T),
 T gamma_zero)
{
    if (nbh.type() == IA_InvariantNbh<IA_Point<int>,IA_Point<int> >::s_type()) {
	return backw_nbh_inv1
	    (img,
	     ((IA_InvariantNbh<IA_Point<int>,IA_Point<int> >*) IA_FBI<IA_Point<int>,IA_Point<int>, T, T, T>::extract_baseptr(nbh))->value,
	     dest_ps, gamma, gamma_zero);
    } else {
	T *const	dest_data = new T[dest_ps.card()];
	T	*valp = dest_data;
	IA_PSIter<IA_Point<int> >	dest_iter(dest_ps);
	IA_Point<int>	base_ip;
	while (dest_iter(base_ip)) {
	    IA_Set<IA_Point<int> >	ps = nbh(base_ip);

	    IA_PSIter<IA_Point<int> >	nbh_iter(ps);
	    IA_Point<int>	nbh_ip;

	    *valp = gamma_zero;

	    while ( nbh_iter(nbh_ip) ) {
		if (! img.domain().contains(nbh_ip))
		    continue;
		*valp = gamma(img(nbh_ip), *valp);
	    }
	    valp++;
	}
	return IA_Image<IA_Point<int>,T>(dest_ps, dest_data,
					   dest_ps.card(), 1);
    }
}

//
// Generic backward neighborhood reduction
//

template <class S, class T>
IA_Image<IA_Point<int>, S> backw_nbh_inv_core
(IA_Point<int> src_infimum,
 IA_Point<int> src_width,
 const T* src_data,
 const IA_Set<IA_Point<int> > &nbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 S (*Gamma)(T*, unsigned count),
 const S*)
{
    const int	dimen = src_width.dim();

    int		nbh_sz = nbh.card();
    int	*const nbh_offsets = new int[nbh_sz];

    {
	int	*o_scan = nbh_offsets;
	IA_PSIter<IA_Point<int> >	iter(nbh);
	IA_Point<int>	ip;
	while (iter(ip)) {
	    *o_scan = ip[0];
	    for (unsigned i=1; i<dimen; i++) {
		*o_scan *= src_width[i];
		*o_scan += ip[i];
	    }
	    o_scan++;
	}
    }

    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    S *const	dest_data = new S[dest_ps.card()];
    S *	valp = dest_data;
    T	*a = new T[nbh_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<nbh_sz; i++) {
	    a[i] = base[nbh_offsets[i]];
	}
	*(valp++) = Gamma(a, nbh_sz);
    }
    delete[] a;
    delete[] nbh_offsets;

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


template <class S, class T>
IA_Image<IA_Point<int>, S> backw_nbh_inv
(const IA_Image<IA_Point<int>, T> &img,
 const IA_Set<IA_Point<int> > &invnbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 S (*Gamma)(T*, unsigned count),
 T gamma_zero,
 const S*)
{
    IA_Point<int>	inf_ = dest_ps.inf()+invnbh.inf();
    IA_Point<int>	sup_ = dest_ps.sup()+invnbh.sup();
    IA_Set<IA_Point<int> >	src_ps = IA_boxy_pset(inf_,sup_);

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

    zero_extend(img, src_ps, (T*)src_data, gamma_zero);
    IA_Image<IA_Point<int>,S>	rval =
	backw_nbh_inv_core(inf_, sup_-inf_ + 1, src_data, invnbh, dest_ps,
			   Gamma, (S*)0);

    delete[] src_data;

    return rval;
}

template <class S, class T>
IA_Image<IA_Point<int>,S> neighborhood_reduction
(const IA_Image<IA_Point<int>, T> &img,
 const IA_Neighborhood<IA_Point<int>,IA_Point<int> > &nbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 T gamma_zero,
 S (*Gamma)(T*, unsigned count),
 const S*)
{
    if (nbh.type() == IA_InvariantNbh<IA_Point<int>,IA_Point<int> >::s_type()) {
	return backw_nbh_inv
	    (img,
	     ((IA_InvariantNbh<IA_Point<int>,IA_Point<int> >*) IA_FBI<IA_Point<int>,IA_Point<int>, T, S, T>::extract_baseptr(nbh))->value,
	     dest_ps, Gamma, gamma_zero, (const S*)0);
    } else {
	S *const	dest_data = new S[dest_ps.card()];
	S	*valp = dest_data;
	IA_PSIter<IA_Point<int> >	dest_iter(dest_ps);
	IA_Point<int>	base_ip;
	while (dest_iter(base_ip)) {
	    IA_Set<IA_Point<int> >	ps = nbh(base_ip);

	    IA_PSIter<IA_Point<int> >	nbh_iter(ps);
	    IA_Point<int>	nbh_ip;

	    T	*const a = new T[ps.card()];
	    unsigned	count=0;

	    while ( nbh_iter(nbh_ip) ) {
		if (! img.domain().contains(nbh_ip))
		    continue;
		a[count++] = img(nbh_ip);
	    }
	    *(valp++) = Gamma(a, count);
	    delete[] a;
	}
	return IA_Image<IA_Point<int>,S>(dest_ps, dest_data,
					   dest_ps.card(), 1);
    }
}

//
// Generic backward neighborhood reduction
//

template <class S, class T>
IA_Image<IA_Point<int>, S> backw_nbh_inv_core3
(IA_Point<int> src_infimum,
 IA_Point<int> src_width,
 const T* src_data,
 const IA_Set<IA_Point<int> > &src_ps,
 const IA_Set<IA_Point<int> > &nbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 S (*Gamma)(T*, unsigned count),
 const S*)
{
    const int	dimen = src_width.dim();

    int		nbh_sz = nbh.card();
    int	*const nbh_offsets = new int[nbh_sz];
    IA_Point<int>	*nbh_points = new IA_Point<int>[nbh_sz];

    {
	int	j;
	IA_PSIter<IA_Point<int> >	iter(nbh);
	for (j=0; iter(nbh_points[j]); j++) {
	    nbh_offsets[j] = nbh_points[j][0];
	    for (unsigned i=1; i<dimen; i++) {
		nbh_offsets[j] *= src_width[i];
		nbh_offsets[j] += nbh_points[j][i];
	    }
	}
    }

    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    S *const	dest_data = new S[dest_ps.card()];
    S *	valp = dest_data;
    T	*a = new T[nbh_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<nbh_sz; i++) {
	    if (src_ps.contains(ip+nbh_points[i]))
		a[count++] = base[nbh_offsets[i]];
	}
	*(valp++) = Gamma(a, count);
    }
    delete[] a;
    delete[] nbh_points;
    delete[] nbh_offsets;

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


template <class S, class T>
IA_Image<IA_Point<int>, S> backw_nbh_inv3
(const IA_Image<IA_Point<int>, T> &img,
 const IA_Set<IA_Point<int> > &invnbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 S (*Gamma)(T*, unsigned count),
 const S*)
{
    IA_Point<int>	inf_ = dest_ps.inf()+invnbh.inf();
    IA_Point<int>	sup_ = dest_ps.sup()+invnbh.sup();
    IA_Set<IA_Point<int> >	src_ps = IA_boxy_pset(inf_,sup_);

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

    zero_extend(img, src_ps, (T*)src_data, gamma_zero);
    IA_Image<IA_Point<int>,S>	rval =
	backw_nbh_inv_core3(inf_, sup_-inf_ + 1, src_data, img.domain(),
			    invnbh, dest_ps,
			    Gamma, (S*)0);

    delete[] src_data;

    return rval;
}

template <class S, class T>
IA_Image<IA_Point<int>,S> neighborhood_reduction
(const IA_Image<IA_Point<int>, T> &img,
 const IA_Neighborhood<IA_Point<int>,IA_Point<int> > &nbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 S (*Gamma)(T*, unsigned count),
 const S*)
{
    if (nbh.type() == IA_InvariantNbh<IA_Point<int>,IA_Point<int> >::s_type()) {
	return backw_nbh_inv3
	    (img,
	     ((IA_InvariantNbh<IA_Point<int>,IA_Point<int> >*) IA_FBI<IA_Point<int>,IA_Point<int>, T, S, T>::extract_baseptr(nbh))->value,
	     dest_ps, Gamma, (const S*)0);
    } else {
	S *const	dest_data = new S[dest_ps.card()];
	S	*valp = dest_data;
	IA_PSIter<IA_Point<int> >	dest_iter(dest_ps);
	IA_Point<int>	base_ip;
	while (dest_iter(base_ip)) {
	    IA_Set<IA_Point<int> >	ps = nbh(base_ip);

	    IA_PSIter<IA_Point<int> >	nbh_iter(ps);
	    IA_Point<int>	nbh_ip;

	    T	*const a = new T[ps.card()];
	    unsigned	count=0;

	    while ( nbh_iter(nbh_ip) ) {
		if (! img.domain().contains(nbh_ip))
		    continue;
		a[count++] = img(nbh_ip);
	    }
	    *(valp++) = Gamma(a, count);
	    delete[] a;
	}
	return IA_Image<IA_Point<int>,S>(dest_ps, dest_data,
					   dest_ps.card(), 1);
    }
}

//
// Generic forward neighborhood reduction
//

template <class T>
IA_Image<IA_Point<int>, T> forw_nbh_inv_core1
(IA_Point<int> src_infimum,
 IA_Point<int> src_width,
 const T* src_data,
 const IA_Set<IA_Point<int> > &nbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 T (*gamma)(T, T),
 T gamma_zero)
{
    const int	dimen = src_width.dim();

    int		nbh_sz = nbh.card();
    int	*const nbh_offsets = new int[nbh_sz];

    {
	int	*o_scan = nbh_offsets;
	IA_PSIter<IA_Point<int> >	iter(nbh);
	IA_Point<int>	ip;
	while (iter(ip)) {
	    *o_scan = -ip[0];
	    for (unsigned i=1; i<dimen; i++) {
		*o_scan *= src_width[i];
		*o_scan -= ip[i];
	    }
	    o_scan++;
	}
    }

    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    T *const	dest_data = new T[dest_ps.card()];
    T *	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<nbh_sz; i++) {
	    *valp = gamma(base[nbh_offsets[i]], *valp);
	}
	valp++;
    }
    delete[] nbh_offsets;

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


template <class T>
IA_Image<IA_Point<int>, T> forw_nbh_inv1
(const IA_Image<IA_Point<int>, T> &img,
 const IA_Set<IA_Point<int> > &invnbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 T (*gamma)(T, T),
 T gamma_zero)
{
    IA_Point<int>	inf_ = dest_ps.inf()-invnbh.sup();
    IA_Point<int>	sup_ = dest_ps.sup()-invnbh.inf();
    IA_Set<IA_Point<int> >	src_ps = IA_boxy_pset(inf_,sup_);

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

    zero_extend(img, src_ps, (T*)src_data, gamma_zero);
    IA_Image<IA_Point<int>,T>	rval =
	forw_nbh_inv_core1(inf_, sup_-inf_ + 1, src_data, invnbh, dest_ps,
			    gamma, gamma_zero);

    delete[] src_data;

    return rval;
}

template <class T>
IA_Image<IA_Point<int>,T> neighborhood_reduction
(const IA_Neighborhood<IA_Point<int>,IA_Point<int> > &nbh,
 const IA_Image<IA_Point<int>, T> &img,
 const IA_Set<IA_Point<int> > &dest_ps,
 T (*gamma)(T, T),
 T gamma_zero)
{
    if (nbh.type() == IA_InvariantNbh<IA_Point<int>,IA_Point<int> >::s_type()) {
	return forw_nbh_inv1
	    (img,
	     ((IA_InvariantNbh<IA_Point<int>,IA_Point<int> >*) IA_FBI<IA_Point<int>,IA_Point<int>, T, T, T>::extract_baseptr(nbh))->value,
	     dest_ps, gamma, gamma_zero);
    } else {
	unsigned	size = dest_ps.card();
	T *const	dest_data = new T[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_Set<IA_Point<int> >	ps = nbh(base_ip);

	    IA_PSIter<IA_Point<int> >	nbh_iter(ps);
	    IA_Point<int>	nbh_ip;
	    T	ival = img(base_ip);

	    while ( nbh_iter(nbh_ip) ) {
		if (! dest_ps.contains(nbh_ip))
		    continue;
		unsigned offset = dest_ps.index(nbh_ip);
		dest_data[offset] = gamma(ival, dest_data[offset]);
	    }
	}
	return IA_Image<IA_Point<int>,T>(dest_ps, dest_data,
				       dest_ps.card(), 1);
    }
}

//
// Generic forward neighborhood reduction
//

template <class S, class T>
IA_Image<IA_Point<int>, S> forw_nbh_inv_core
(IA_Point<int> src_infimum,
 IA_Point<int> src_width,
 const T* src_data,
 const IA_Set<IA_Point<int> > &nbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 S (*Gamma)(T*, unsigned count),
 const S*)
{
    const int	dimen = src_width.dim();

    int		nbh_sz = nbh.card();
    int	*const nbh_offsets = new int[nbh_sz];

    {
	int	*o_scan = nbh_offsets;
	IA_PSIter<IA_Point<int> >	iter(nbh);
	IA_Point<int>	ip;
	while (iter(ip)) {
	    *o_scan = -ip[0];
	    for (unsigned i=1; i<dimen; i++) {
		*o_scan *= src_width[i];
		*o_scan -= ip[i];
	    }
	    o_scan++;
	}
    }

    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    S *const	dest_data = new S[dest_ps.card()];
    S *	valp = dest_data;
    T	*a = new T[nbh_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<nbh_sz; i++) {
	    a[i] = base[nbh_offsets[i]];
	}
	*(valp++) = Gamma(a, nbh_sz);
    }
    delete[] a;
    delete[] nbh_offsets;

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


template <class S, class T>
IA_Image<IA_Point<int>, S> forw_nbh_inv
(const IA_Image<IA_Point<int>, T> &img,
 const IA_Set<IA_Point<int> > &invnbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 S (*Gamma)(T*, unsigned count),
 T gamma_zero,
 const S*)
{
    IA_Point<int>	inf_ = dest_ps.inf()-invnbh.sup();
    IA_Point<int>	sup_ = dest_ps.sup()-invnbh.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), src_ps, (T*)src_data, gamma_zero);
    IA_Image<IA_Point<int>,S>	rval =
	forw_nbh_inv_core(inf_, sup_-inf_ + 1, src_data, invnbh, dest_ps,
			   Gamma, (S*)0);

    delete[] src_data;

    return rval;
}

template <class S, class T>
IA_Image<IA_Point<int>,S> neighborhood_reduction
(const IA_Neighborhood<IA_Point<int>,IA_Point<int> > &nbh,
 const IA_Image<IA_Point<int>, T> &img,
 const IA_Set<IA_Point<int> > &dest_ps,
 T gamma_zero,
 S (*Gamma)(T*, unsigned count),
 const S*)
{
    if (nbh.type() == IA_InvariantNbh<IA_Point<int>,IA_Point<int> >::s_type()) {
	return forw_nbh_inv
	    (img,
	     ((IA_InvariantNbh<IA_Point<int>,IA_Point<int> >*) IA_FBI<IA_Point<int>,IA_Point<int>, T, S, T>::extract_baseptr(nbh))->value,
	     dest_ps, Gamma, gamma_zero, (const S*)0);
    } else {
	IA::not_yet_implemented(__FILE__,__LINE__);
    }
}

//
// Generic forward neighborhood reduction
//

template <class S, class T>
IA_Image<IA_Point<int>, S> forw_nbh_inv_core3
(IA_Point<int> src_infimum,
 IA_Point<int> src_width,
 const T* src_data,
 const IA_Set<IA_Point<int> > &src_ps,
 const IA_Set<IA_Point<int> > &nbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 S (*Gamma)(T*, unsigned count),
 const S*)
{
    const int	dimen = src_width.dim();

    int		nbh_sz = nbh.card();
    int	*const nbh_offsets = new int[nbh_sz];
    IA_Point<int>	*nbh_points = new IA_Point<int>[nbh_sz];

    {
	int	j;
	IA_PSIter<IA_Point<int> >	iter(nbh);
	for (j=0; iter(nbh_points[j]); j++) {
	    nbh_offsets[j] = -nbh_points[j][0];
	    for (unsigned i=1; i<dimen; i++) {
		nbh_offsets[j] *= src_width[i];
		nbh_offsets[j] -= nbh_points[j][i];
	    }
	}
    }


    IA_PSIter<IA_Point<int> >	iter(dest_ps);
    IA_Point<int>	ip;
    S *const	dest_data = new S[dest_ps.card()];
    S *	valp = dest_data;
    T	*a = new T[nbh_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<nbh_sz; i++) {
	    if (src_ps.contains(ip - nbh_points[i]))
		a[count++] = base[nbh_offsets[i]];
	}
	*(valp++) = Gamma(a, count);
    }
    delete[] a;
    delete[] nbh_points;
    delete[] nbh_offsets;

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


template <class S, class T>
IA_Image<IA_Point<int>, S> forw_nbh_inv3
(const IA_Image<IA_Point<int>, T> &img,
 const IA_Set<IA_Point<int> > &invnbh,
 const IA_Set<IA_Point<int> > &dest_ps,
 S (*Gamma)(T*, unsigned count),
 const S*)
{
    IA_Point<int>	inf_ = dest_ps.inf()-invnbh.sup();
    IA_Point<int>	sup_ = dest_ps.sup()-invnbh.inf();
    IA_Set<IA_Point<int> >	src_ps = IA_boxy_pset(inf_,sup_);

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

    T	gamma_zero;
    zero_extend(img, src_ps, (T*)src_data, gamma_zero);
    IA_Image<IA_Point<int>,S>	rval =
	forw_nbh_inv_core3(inf_, sup_-inf_ + 1, src_data, img.domain(),
			   invnbh, dest_ps,
			   Gamma, (S*)0);

    delete[] src_data;

    return rval;
}

template <class S, class T>
IA_Image<IA_Point<int>,S> neighborhood_reduction
(const IA_Neighborhood<IA_Point<int>,IA_Point<int> > &nbh,
 const IA_Image<IA_Point<int>, T> &img,
 const IA_Set<IA_Point<int> > &dest_ps,
 S (*Gamma)(T*, unsigned count),
 const S*)
{
    if (nbh.type() == IA_InvariantNbh<IA_Point<int>,IA_Point<int> >::s_type()) {
	return forw_nbh_inv3
	    (img,
	     ((IA_InvariantNbh<IA_Point<int>,IA_Point<int> >*) IA_FBI<IA_Point<int>,IA_Point<int>, T, S, T>::extract_baseptr(nbh))->value,
	     dest_ps, Gamma, (const S*)0);
    } else {
	IA::not_yet_implemented(__FILE__,__LINE__);
    }
}
