#include "PSIter.h"

template <class T>
IA_BaseImage<IA_Point<int>,T> *restrictD_toD
  (
   const IA_BaseImage<IA_Point<int>,T> *im ,
   const IA_Set<IA_Point<int> > &ips
   )
{
    IA_Set<IA_Point<int> >	result_domain = im->domain() & ips;
    // intersect_with_structure(im.domain(), ips, &ss);
#if 0
    if (result_domain == im->domain())
	// clue to reuse this object
	return 0;
#endif

    const unsigned	sz = result_domain.card();
    T * const	dest = new T[sz];
    {
	T	*scan = dest;
	IA_PSIter<IA_Point<int> >	iter(result_domain);
	IA_Point<int>	ip;
	while (iter(ip))
	    *(scan++) = (*im)(ip);
    }

    return new IA_VectorI<IA_Point<int>,T>(result_domain, dest, sz, 1);
  
}

///

template <class T>
IA_BaseImage<IA_Point<int>,T> *restrictC_toD
  (
   const IA_BaseImage<IA_Point<double>,T> *im,
   const IA_Set<IA_Point<int> > &ips
   )
{
  IA_PSIter<IA_IntPoint>	iter(ips);
  int	sz = ips.card();
  int	len=0;
  IA_Point<int>	*aip = new IA_Point<int>[sz];
  T	*vals = new T[sz];
  int	rejects=0;

  while (iter(aip[len])) {
    if (!im->domain().contains(aip[len])) {
      rejects++;
      continue;
    }

    vals[len] = (*im)(aip[len]);
    len++;
  }

  IA_Set<IA_Point<int> >	newips(ips);
  if (rejects)
      newips = IA_Set<IA_Point<int> >(ips.dim(), aip, len);

  delete[] aip;

  return new IA_VectorI<IA_Point<int>,T>(newips, vals, len, 1);
}
