VoxelCoding.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 _VoxelCoding_hpp
#define _VoxelCoding_hpp
#include<ImLib3D/MorphologicalOperators.hpp>
#include<ImLib3D/Threshold.hpp>


namespace IP3D
{
        void VoxelCode(const Mask3D &shape,const Mask3D &seeds,int m0,int m1,int m2,LabelImage3D& res,
                      const Vect3Di *stopPoint=NULL);
};

namespace IP3D
{
    inline void SSCode(const Mask3D &shape,const Vect3Di &sourcePoint,int m0,int m1,int m2,LabelImage3D& res,
                      const Vect3Di *stopPoint=NULL)
    {
        Mask3D seeds(shape.Size());
        seeds.Fill(0);
        seeds(sourcePoint)=1;
        VoxelCode(shape,seeds,m0,m1,m2,res,stopPoint);
    }
};
namespace IP3D
{
    inline void BSCode(const Mask3D &shape,int m0,int m1,int m2,LabelImage3D& res)
    {
        Mask3D seeds=shape;
        Mask3D t;
        Erosion(shape,MORPHO_Cross7,t);
        seeds-=t;
        VoxelCode(shape,seeds,m0,m1,m2,res);
    }
};
class CodeCluster;

namespace IP3D
{
namespace ShortestPathExtractionHelpers
{
    typedef bool (*StopFctType)(void *dat,const Vect3Di &p);
    typedef pair<StopFctType,void*> StopCBType;
}

    void ShortestPathExtraction(const Mask3D &shape,const Vect3Di &start,const Vect3Di &end,vector<Vect3Di> &path,
                                Mask3D *pres=NULL,StructureElement *psel=NULL,
                                const LabelImage3D *pssCode0=NULL,
                                ShortestPathExtractionHelpers::StopCBType stopCB=
                                ShortestPathExtractionHelpers::StopCBType(NULL,NULL));


};


class CodeCluster
{
    friend class SkeletonExtraction;
public:
    int value;
    Vect3Df center;
    vector<CodeCluster *> neighbors;
    vector<Vect3Di> points;
    bool isLM;
    bool isMerging;
    bool isDividing;
    bool isBranching;

    int stop;
private:
    CodeCluster(int _value):
        value(_value),stop(0)
    {
    }
    void Init()
    {
        center=points[0];
        for(uint i=1;i<points.size();i++){center+=points[i];}
        center/=points.size();
    }
    void FindType()
    {
        isLM=true;
        for(uint i=0;i<neighbors.size();i++)
        {
            if(neighbors[i]->value>value){isLM=false;break;}
        }
        int ct=0;
        for(uint i=0;i<neighbors.size();i++)
        {
            if(neighbors[i]->value==value-1){ct++;}
        }
        if(ct>=2){isMerging=true;}else{isMerging=false;}
        ct=0;
        for(uint i=0;i<neighbors.size();i++)
        {
            if(neighbors[i]->value==value+1){ct++;}
        }
        if(ct>=2){isDividing=true;}else{isDividing=false;}
        if(isDividing || isMerging){isBranching=true;}else{isBranching=false;}
    }
public:
    void GetSuccessors(vector<CodeCluster *> &successors)
    {
        for(uint i=0;i<neighbors.size();i++)
        {
            if(neighbors[i]->value<value){successors.push_back(neighbors[i]);}
        }
    } 
    vector<Vect3Di> medialPoints;
    void FindMedialPoints(const LabelImage3D &bsCode)
    {
        int mx=bsCode(points[0]);
        for(uint i=1;i<points.size();i++){mx=max(mx,bsCode(points[i]));}
        for(uint i=0;i<points.size();i++){if(bsCode(points[i])==mx){medialPoints.push_back(points[i]);}}
    }
    Vect3Di MedialPoint(const Vect3Df &ref=Vect3Di(0,0,0))
    {
        int imin=0;
        float d2=ref.Dist2(medialPoints[0]);
        for(uint i=1;i<medialPoints.size();i++)
        {
            float d=ref.Dist2(medialPoints[i]);
            if(d<d2){d2=d;imin=i;}
        }
        return medialPoints[imin];
    }
    void FindNeighbors(Container3D<CodeCluster *> &clustersIma)
    {
        // for each point in this cluster, check all if each of its neighbors is in another cluster
        // and add it if it's not there
        StructureElement sel(MORPHO_Cubic27);
        sel.RemoveCenter();
        for(uint pt=0;pt<points.size();pt++)
        {
            for(uint i=0;i<sel.size();i++)
            {
                Vect3Di npos=points[pt]+sel[i];
                if(!clustersIma.IsInside(npos)){continue;}
                CodeCluster *neighbor=clustersIma(npos);
                if(neighbor && neighbor->value!=value && abs(neighbor->value-value)<2)
                {
                    if(find(neighbors.begin(),neighbors.end(),neighbor)==neighbors.end())
                    {
                        neighbors.push_back(neighbor);
                    }
                }
            }
                
        }
//          printf("neighbors for cluster %x %03d\n",this,value);
//          for(uint i=0;i<neighbors.size();i++)
//          {
//              printf("cluster %x %03d:neighbor %2d: %x %03d\n",this,value,i,neighbors[i],neighbors[i]->value);
//          }
    }
    
};
class PathGraph
{
public:
    class Node;
    class Path
    {
        friend class PathGraph;
    public:
        Node *first,*last;
        vector<Vect3Di> points;
    private:
        Path(Node *_first,Node *_last,const vector<Vect3Di> &_points);
    };
    class Node
    {
        friend class PathGraph;
    public:
        Vect3Di point;
        vector<Path *> paths;
    private:
        Node(){;}
    };
    vector<Node *> nodes;
    vector<Path *> paths;
    Node *AddNode(){Node *node=new Node();nodes.push_back(node);return node;}
    RectZone3Di BoundingBox() const;
    Path *AddPath(Node *_first,Node *_last,const vector<Vect3Di> &_points)
    {
        Path *path=new Path(_first,_last,_points);paths.push_back(path);return path;
    }
};

class SkeletonExtraction
{
    static bool SpeStopCB(void *dat,const Vect3Di &p)
    {
        return ((SkeletonExtraction *)dat)->clustersIma(p)->stop!=-1;
    }
    class MedPathFind
    {
    public:
        CodeCluster *parentCluster;
        int id;
        int id1,id2;
        vector<Vect3Di> path;
        vector<CodeCluster *> cpath;
        vector<vector<int> > branches;
        vector<PathGraph::Node *> branchNodes;
        void ExtractPath(SkeletonExtraction *se);
        int FindCluster(CodeCluster *cluster)
        {
            int res=-1;
            for(uint i=0;i<cpath.size();i++){if(cpath[i]==cluster){res=i;break;}}
            return res;
        }
        int operator<(const MedPathFind &other) const {return cpath[0]->value<other.cpath[0]->value;}
        MedPathFind(CodeCluster *startCluster):
            id1(-1),id2(-1)
        {
            cpath.push_back(startCluster);
        }
    };
    friend class SkeletonExtraction::MedPathFind;
    const Mask3D &shape;
    LabelImage3D ssCode;
    LabelImage3D bsCode;
    vector<CodeCluster *> clusters;
    Container3D<CodeCluster *> clustersIma;
    void CreateClusters();
    void InitClusters();

    void PathGraphFromMedPaths(vector<MedPathFind> &medPaths,PathGraph &pathGraph);
    void ExtractAllMedPaths(vector<MedPathFind> &medPaths);
    void ExtractPaths(PathGraph &pathGraph);

    SkeletonExtraction(const Mask3D &_shape,const Vect3Di *psourcePoint=NULL);
    void MedPathFromPath(const vector<Vect3Di> &path,vector<Vect3Di> &medPath)
    {
        medPath.push_back(path[0]);
        for(uint i=1;i<path.size();i++)
        {
            medPath.push_back(clustersIma(path[i])->MedialPoint(medPath[i-1]));
        }
    }
public:
    static void ExtractAll(const Mask3D &shape,PathGraph &pathGraph);
    static void ExtractLinkingMedPath(const Mask3D &shape,const Vect3Di &start,const Vect3Di &end,vector<Vect3Di> &res);
    static void DrawPathGraph(const PathGraph &pathGraph,Mask3D &res);
};

namespace IP3D
{
    inline void Skeleton(const Mask3D &shape,Mask3D& res)
    {
        PathGraph pathGraph;
        SkeletonExtraction::ExtractAll(shape,pathGraph);
        res.Resize(shape.Size());
        res.Fill(0);
        SkeletonExtraction::DrawPathGraph(pathGraph,res);
    }
};

#endif// _VoxelCoding_hpp

---