Interpolator3D.hpp

Go to the documentation of this file.

/* 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.
 */

#ifndef _Interpolator3D_hpp
#define _Interpolator3D_hpp
#include<gsl/gsl_sf.h>


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 );
    }

    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)
    {;}

};

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) (1-fabs(x-(x0+i)))*(1-fabs(y-(y0+j)))*(1-fabs(z-(z0+k)))*ima.SafeValue(x0+i,y0+j,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;
    }
    
};


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= Zero<Im3DValue>();
        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+=coef*(*p);
            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

---