/* ImLib3D
 * Copyright (c) 2001, ULP-IPB Strasbourg.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
/*! \file Interpolator3D.hpp */
//! Interpolator3D
#ifndef _Interpolator3D_hpp
#define _Interpolator3D_hpp
#include<gsl/gsl_sf.h>

//! Virtual base class for interpolators
/*!
  Interpolators are associated to images using SetInterpolator of the Image3Dlinear classes.
 */
template<class Im3DValue>
class Interpolator3D
{
	typedef Image3Dlinear<Im3DValue> ImageType;
// support
public:
	static Interpolator3D *Create(const string &name);
	RectZone3Df support;
	RectZone3Di GetRectZone(float x, float y, float z)
	{
		RectZone3Df s1(support);
		s1.Translate(x,y,z);
		return RectZone3Di( s1, 0 );
	}
	//! Certains interpolators (ie. Spline) need to do precomputations on original image
	/*! Call this with originale image to do those precomputations,
	    Call it with null, to clean up.
	 */
	virtual void Reset(const ImageType *srcImage){;}
	virtual Im3DValue Value(const ImageType &ima,const float x,const float y,const float z)=0;
	virtual ~Interpolator3D(){;}
	Interpolator3D():support(0,0,0,0,0,0)
	{;}

};

//! Nearest neighbour interpolator
template<class Im3DValue>
class NNInterpolator : public Interpolator3D<Im3DValue>
{
	typedef Image3Dlinear<Im3DValue> ImageType;
public:
	virtual Im3DValue Value(const ImageType &ima,const float x,const float y,const float z)
	{
		return ima.SafeValue((int)(x+.5),(int)(y+.5),(int)(z+.5));
	}
	
};

template<class Im3DValue>
class TriLinearInterpolator : public Interpolator3D<Im3DValue>
{
	typedef Image3Dlinear<Im3DValue> ImageType;
public:
	virtual Im3DValue Value(const Image3Dlinear<Im3DValue> &ima,const float x,const float y,const float z)
	{
//  		float dx=x-floor(x);
//  		float dy=y-floor(y);
//  		float dz=z-floor(z);
		int x0=(int)x;
		int y0=(int)y;
		int z0=(int)z;
#define _linearincr(i,j,k) ima.SafeValue(x0+i,y0+j,z0+k)*(1-fabs(x-(x0+i)))*(1-fabs(y-(y0+j)))*(1-fabs(z-(z0+k)));
		Im3DValue res=_linearincr(0,0,0);
		res+=_linearincr(1,0,0);
		res+=_linearincr(0,1,0);
		res+=_linearincr(1,1,0);
		res+=_linearincr(0,0,1);
		res+=_linearincr(1,0,1);
		res+=_linearincr(0,1,1);
		res+=_linearincr(1,1,1);
		return res;
#undef _linearincr
//  		int tx = (int)(x+.5); int ty = (int)(y+.5); int tz = (int)(z+.5);
//  		float distx = fabs(x - tx);
//  		float disty = fabs(y - ty);
//  		float distz = fabs(z - tz);

//  		Im3DValue  vala = ima.SafeValue(tx  , ty  , tz  );
//  		Im3DValue  valb = ima.SafeValue(tx  , ty  , tz+1);
//  		Im3DValue  valc = ima.SafeValue(tx  , ty+1, tz  );
//  		Im3DValue  vald = ima.SafeValue(tx  , ty+1, tz+1);
//  		Im3DValue  vale = ima.SafeValue(tx+1, ty  , tz);
//  		Im3DValue  valf = ima.SafeValue(tx+1, ty  , tz+1);
//  		Im3DValue  valg = ima.SafeValue(tx+1, ty+1, tz  );
//  		Im3DValue  valh = ima.SafeValue(tx+1, ty+1, tz+1);
		
//  		Im3DValue val = 
//  			(vale - vala)*distx + 
//  			(valc - vala)*disty + 
//  			(valb - vala)*distz +
//  			(valg - vale - valc + vala)*distx*disty +
//  			(vald - valc - valb + vala)*disty*distz +
//  			(valf - vale - valb + vala)*distx*distz +
//  			(valh + vale + valc + valb - valg - valf - vald - vala)*distx*disty*distz + 
//  			vala;
//		return val;
	}
	
};


//! Truncated sinc interpolator. You should use spline unless you need huge precision 
template<class Im3DValue>
class TruncatedSinCInterpolator : public Interpolator3D<Im3DValue>
{
	typedef Image3Dlinear<Im3DValue> ImageType;
public:
	float radius;
	virtual Im3DValue Value(const Image3Dlinear<Im3DValue> &ima,const float x,const float y,const float z)
	{
		RectZone3Di zone((int) ceil(x-radius),(int) ceil(y-radius),(int) ceil(z-radius),
						 (int)floor(x+radius),(int)floor(y+radius),(int)floor(z+radius) );
		zone=zone.Intersection(ima.GetZone());
		Im3DValue res=ima.GetZero();
		if(!zone.IsOk()) {return res;}
		typename Image3Dlinear<Im3DValue>::const_iteratorZone p(zone);
		double totcoef=0;
		vector<vector<double> > dirCoefs;
		dirCoefs.assign(3,vector<double>());
		Vect3Di size=zone.SizeV();
		Vect3Di p0=zone.GetP0();
		Vect3Df pos(x,y,z);
		for(int i=0;i<3;i++)
		{
			dirCoefs[i].assign(size(i),0);
			for(int j=0;j<size(i);j++)
			{
				dirCoefs[i][j]=gsl_sf_sinc(pos(i)-(j+p0(i)));
			}
		}
		float r02=radius*radius;
		for(p=ima.begin();p!=ima.end();++p)
		{
  			float r2=(x-p.x)*(x-p.x)+(y-p.y)*(y-p.y)+(z-p.z)*(z-p.z);
			if(r2>r02){continue;}
			double coef=
				dirCoefs[0][p.x-p0.x]*
				dirCoefs[1][p.y-p0.y]*
				dirCoefs[2][p.z-p0.z] ;
			res+=*p*coef;
			totcoef+=coef;
		}
		// renorm (trunc sinc is not a true interpolating fct!!)
		res/=totcoef;
		return res;
	}
	TruncatedSinCInterpolator(float _radius=5):radius(_radius){;}
};


template<class Im3DValue>
inline Im3DValue
Image3Dlinear<Im3DValue>::Value( float x, float y, float z ) const 
{
	if(!interpolator){ThrowError("Image3Dlinear<Im3DValue>::Value(float,float,float) interpolator==NULL image:%x",this);}
	return interpolator->Value(*this,x,y,z);
}

// to avoid code bloat include this manually!
//#include<ImLib3D/Interpolator3D.hxx>
#endif // _Interpolator3D_hpp
