Filter.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 _Filter_hpp
#define _Filter_hpp

#include<iterator>
#include<ImLib3D/Image3Dlinear.hpp>
#include<algorithm>


// ****************** Filter ************************ 

// **************************************************
class Filter : public Image3Df
{
private:
//  Vect3Di center;
public:

    Filter(const Filter&  fi, const Vect3Di& centerpoint):Image3Df(fi)
    {properties.iCenter=new Vect3Di(centerpoint);}
    Filter(const Size3D & size, const Vect3Di& centerpoint):Image3Df(size)
    {properties.iCenter=new Vect3Di(centerpoint);}
    Filter(int _width, int _height, int _depth, const Vect3Di centerpoint):Image3Df(_width, _height, _depth)
    {properties.iCenter=new Vect3Di(centerpoint);}
    Filter(const string& fname):Image3Df()
    {Streamable::ReadFromFile(fname);}

    Vect3Di GetCenter() const {return *properties.iCenter;}
    void SetCenter(const Vect3Di& centerpoint){*properties.iCenter = centerpoint;}
};


class SeparableFilter
{
protected:
    void Resize()
    {
        xcomponent.assign(size.width,0);
        ycomponent.assign(size.height,0);
        zcomponent.assign(size.depth,0);
    }
public:
    SeparableFilter(Size3D  _size,Vect3Di _iCenter):size(_size),iCenter(_iCenter){Resize();}
    SeparableFilter(Vect3Di _size,Vect3Di _iCenter):size(_size),iCenter(_iCenter){Resize();}
    vector<float> xcomponent;
    vector<float> ycomponent;
    vector<float> zcomponent;
    Size3D size;
    Vect3Di iCenter;
};

template<class RealType>
void BinomialFilterMask(uint order, vector<RealType>& filtermask)
{
    if (order == 0)
    {
        filtermask.assign(1,1);
        return;
    }
    if (order == 1)
    {
        filtermask.assign(2,1);
        return;
    }
    vector<RealType> upperRow(order+1,0);
    vector<RealType> lowerRow(order+1,0);
    upperRow[0] = 1;
    for (uint i=1; i<order+1; i++)
    {
        lowerRow[0]=1;
        for (uint j=1; j<i; j++)
        {
            lowerRow[j]=upperRow[j-1]+upperRow[j];
        }
        lowerRow[i]=1;
        upperRow = lowerRow;
    }

    double sum=0;
    for (typename vector<RealType>::iterator p = lowerRow.begin(); p!= lowerRow.end(); ++p)
    {
        sum += *p;
    }
    if (sum==0)ThrowError("Runtime error in BinomialFilterMask (sum==0)");

    for (typename vector<RealType>::iterator p = lowerRow.begin(); p!= lowerRow.end(); ++p)
    {
        *p /= sum;
    }
    filtermask = lowerRow;
}

class BinomialFilter : public SeparableFilter
{
public:
    BinomialFilter(int order) : SeparableFilter(Size3D(order+1, order+1, order+1), Vect3Di(order/2, order/2, order/2))
    {
        vector<double> filtermask;
        BinomialFilterMask(order, filtermask);// calculate using double for higher precision(pb with high order filters)

        // prune at +- 4.5 sigma (sigma*sigme=order/4)
        double sigma = sqrt((double)order/4.0);
        int support = (int)(4.5*sigma + 1.0);
        int diff = filtermask.size() - (support*2+1);
        if (diff > 0) // (support*2+1) > filtermask.size()
        {
            vector<double> fullmask = filtermask;
            filtermask.resize(support*2+1);
            for (int i = 0; i < (2*support+1) ; ++i)
                filtermask[i] = fullmask[diff/2 + i];
            if (diff%2)
                filtermask.push_back(fullmask[diff/2+2*support+1]);
        }

        // assign components
        xcomponent.resize(filtermask.size());
        copy(filtermask.begin(),filtermask.end(),xcomponent.begin());
        ycomponent=zcomponent=xcomponent;
        size=Size3D(xcomponent.size(), ycomponent.size(), zcomponent.size());
        //      assert(filtermask.size()%2); // if not impair result will be deplaced one voxel towards the origin
        int c=(filtermask.size()-1)/2;
        iCenter = Vect3Di(c,c,c);
    }
};

inline ostream& operator<<(ostream& s, SeparableFilter& filt)
{
    s << "xVect: " << endl;
    copy(filt.xcomponent.begin(),filt.xcomponent.end(),ostream_iterator<float>(s," "));
    s << endl;
    s << "yVect: " << endl;
    copy(filt.ycomponent.begin(),filt.ycomponent.end(),ostream_iterator<float>(s," "));
    s << endl;
    s << "zVect: " << endl;
    copy(filt.zcomponent.begin(),filt.zcomponent.end(),ostream_iterator<float>(s," "));
    s << endl;
    s << "Center: " << filt.iCenter << endl;
    return s;
}


#endif //_Filter_hpp

---