MorphologicalOperators.hxx

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 _MorphologicalOperators_hxx
#define _MorphologicalOperators_hxx

#include<ImLib3D/Image3D.hpp>
#include<ImLib3D/MorphologicalOperators.hpp>
#include<algorithm>
#include<ImLib3D/ImageProgress.hpp>

template <class ImageType>
void
IP3D::MedianFilter(const ImageType& src, const StructureElement &mask, ImageType& res0)
{
    ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
    args.SameSize();
    ImageType &res=args.GetDest();

    // create vector for containing neighborValues
    vector<typename ImageType::value_type> neighborValues;
    neighborValues.assign(mask.size(),src(0,0,0));
    int imedian=mask.size()/2;

    typename ImageType::iteratorXYZ rp;
    ImageProgress tracker("MedianFilter",rp);
    for(rp=res.begin();rp!=res.end();rp++)
    {
        // find neighborValues
        neighborValues[0]=src.SafeValue(rp.Pos()+mask.GetPoint(0));
        for(uint i=1;i<mask.size();i++)
        {
            neighborValues[i]=src.SafeValue(rp.Pos()+mask[i]);
        }
        // now find median of neighborValues
        sort(neighborValues.begin(),neighborValues.end());
        *rp=neighborValues[imedian];
    }   

    //args.Finished();
}


template <class ImageType>
void
IP3D::SharpeningFilter(const ImageType& src, const StructureElement &mask, ImageType& res0)
{
    typedef typename ImageType::value_type Im3DValue;
    ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
    args.SameSize();
    ImageType &res=args.GetDest();

    // create vector for containing neighborValues
    vector<Im3DValue> neighborValues;
    neighborValues.assign(mask.size(),src(0,0,0));

    typename ImageType::iteratorXYZ rp;
    ImageProgress tracker("SharpeningFilter",rp);
    for(rp=res.begin();rp!=res.end();rp++)
    {
        *rp=src(rp.Pos());
        // find neighborValues
        neighborValues[0]=src.SafeValue(rp.Pos()+mask.GetPoint(0));
        for(uint i=1;i<mask.size();i++)
        {
            neighborValues[i]=src.SafeValue(rp.Pos()+mask[i]);
        }
        Im3DValue min = neighborValues[0];
        Im3DValue max = neighborValues[0];
        int minindex=0;
        int maxindex=0; 
        // find min and max neighbor-value
        for(uint i=1;i<mask.size();i++)
        {
            if (neighborValues[i] > max)
            {max = neighborValues[i];maxindex=i;}
            if (neighborValues[i] < min)
            {min = neighborValues[i];minindex=i;}
        }
        // now assign current to min or max, whichever value is closest:
        if ( ((*rp)-min) < (max-(*rp)))
            *rp = min;
        else if ( ((*rp)-min) > (max-(*rp)))
            *rp = max;
        // or leave unchanged when distance the same
    }   
    //args.Finished();
}


template <class ImageType>
void
IP3D::Erosion(const ImageType& src, const StructureElement &mask, ImageType& res0)
{
    ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
    args.SameSize();
    ImageType &res=args.GetDest();

    ImageNeighbors<const ImageType> neighbors(src,mask);
    uint npts=mask.size();
    typename ImageType::iteratorXYZ rp;
    ImageProgress tracker("Erosion",rp);
    for(rp=res.begin();rp!=res.end();rp++)
    {
        const Vect3Di &pos=rp.Pos();
        if(neighbors.IsSafe(pos))
        {
            neighbors.MoveTo(pos);
            *rp=neighbors[0];
            for(uint i=1;i<npts;i++)
            {
                *rp=min(*rp,neighbors[i]);
            }
        }
        else
        {
            *rp=src.SafeValue(rp.Pos()+mask.GetPoint(0));
            for(uint i=1;i<npts;i++)
            {
                *rp=min(*rp,src.SafeValue(rp.Pos()+mask[i]));
            }
        }
    }   
}

template <class ImageType>
void
IP3D::Dilation(const ImageType& src, const StructureElement &mask, ImageType& res0)
{
    ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
    args.SameSize();
    ImageType &res=args.GetDest(0);
    typename ImageType::iteratorXYZ rp;
    ImageProgress tracker("Dilation",rp);
    for(rp=res.begin();rp!=res.end();rp++)
    {
        *rp=src.SafeValue(rp.Pos()+mask.GetPoint(0));
        for(uint i=1;i<mask.size();i++)
        {
            *rp=max(*rp,src.SafeValue(rp.Pos()+mask[i]));
        }
    }   
}


template <class ImageType>
void
IP3D::MakeValueList(const ImageType& src, ImageValues<typename ImageType::value_type>& res)
{
    typename ImageType::const_iteratorXYZ sp;
    for(sp=src.begin();sp!=src.end();sp++)
    {
        res.values[*sp].point=sp.Pos();
        res.values[*sp].size++;
    }   
}


template <class ImageType>
void
IP3D::ConnectedComponentLabelling(const ImageType& src, const StructureElement &mask0, LabelImage3D& res,typename ImageType::value_type *background0)
{
    // ALGO: For each pixel 
    typedef typename ImageType::value_type Im3DValue;
    pair<bool,Im3DValue> backgroundInfo=DefaultBackground<ImageType>();
    bool hasBackground=backgroundInfo.first;
    Im3DValue background=backgroundInfo.second;
    vector<int> indx2region;
    indx2region.reserve(1000);
    // default indx=0 -> region=0 = background
    indx2region.push_back(0);
    int nbRegions=1;

    res.Resize(src.Size());
    res.Fill(0);
//      for(LabelImage3D::iteratorFast fp=res.begin();fp!=res.end();fp++){if(*fp){*fp=1;}}

    StructureElement mask(mask0,MORPHO_FORWARD);
    LabelImage3D::iteratorXYZ rp;
    vector<int> neigborIndexes;
    neigborIndexes.reserve(mask.size());

    ImageProgress tracker("ConnectedComponentLabelling",rp);
    for(rp=res.begin();rp!=res.end();rp++)
    {
        int talk=0;
        if(rp.Pos()==Vect3Di(48,5,76)){talk=1;}
        const Im3DValue &srcval=src(rp.Pos());
        if(hasBackground && background==srcval){*rp=0;continue;}
        // let's see if this voxel has any neigbors
        neigborIndexes.clear();
        for(uint i=0;i<mask.size();i++)
        {
//          if((currIndex=res.SafeValue(rp.Pos()+mask[i])))
            Vect3Di npos=rp.Pos()+mask[i];
            if(src.IsInside(npos) && src(npos)==srcval)
            {// found a neigbor!
//                      mprintf("found neigbor index:%2d ",currIndex);
//                      (rp.Pos()+mask[i]).Show("pos:","\n");
                int neighborIndex=res(npos);
                neigborIndexes.push_back(neighborIndex);
            }
        }
        if(neigborIndexes.size())
        {// found neigbors!
            if(talk)
            {
                for(uint i=0;i<neigborIndexes.size();i++)
                {
                    mprintf("found neigborIndexes:%d:%d\n",i,neigborIndexes[i]);
                }
            }
            // set this pixel to the same value as one of those neigbors
            int anIndexOfThisRegion=neigborIndexes[0];
            *rp=anIndexOfThisRegion;
            // make all neigbor indexes point to the same region
            int thisRegion=indx2region[anIndexOfThisRegion];
            if(talk)
            {
                mprintf("anIndexOfThisRegion:%d\n",anIndexOfThisRegion);
                mprintf("thisRegion:%d\n",thisRegion);
            }
            for(uint i=1;i<neigborIndexes.size();i++)
            {
                if(talk)
                {
                    mprintf("setting neigborIndexes:%d:%d: previous region:%d\n",i,neigborIndexes[i],indx2region[neigborIndexes[i]]);
                }
                int oldNeighborRegion=indx2region[neigborIndexes[i]];
                if(oldNeighborRegion!=thisRegion)
                {// merge 2 regions
                    for(uint j=0;j<indx2region.size();j++)
                    {
                        if(indx2region[j]==oldNeighborRegion)
                        {
                            indx2region[j]=thisRegion;
                        }
                    }
                }
            }
        }
        else
        {// sorry, no neighbors 
            // this lonesome pixel is a new region and index
            *rp=indx2region.size();
            int thisRegion=nbRegions++;
            indx2region.push_back(thisRegion);
        }
//      }else
//      {
//          *rp=0;// background index
//      }
    }
//      mprintf("CConnectedComponentLabelling:: region image\n");
//      LabelImage3D regions=res;
//      for(LabelImage3D::iteratorFast fp=regions.begin();fp!=regions.end();fp++){*fp=indx2region[*fp];}    
//      ImLib3DDisplay::Show(regions);

    // ***** now compute new unique regions
    uint i;
    map<int,int> old2newregion;
    int nbNewRegions=0;
    // add  background region
    old2newregion[0]=nbNewRegions++;
    for(i=0;i<indx2region.size();i++)
    {
        int old=indx2region[i];
        if(old2newregion.find(old)==old2newregion.end())
        {
            old2newregion[old]=nbNewRegions++;
        }
    }
    // switch to new regions
    for(i=0;i<indx2region.size();i++)
    {
        indx2region[i]=old2newregion[indx2region[i]];
    }   


    // ***** now let's switch from indexes to regions
    for(rp=res.begin();rp!=res.end();rp++)
    {
        *rp=indx2region[*rp];
    }

}
#endif // _MorphologicalOperators_hxx

---