VoxelCoding.cpp

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

#include<ImLib3D/VoxelCoding.hpp>
#include<ImLib3D/Arithmetic.hpp>
#include<ImLib3D/Display.hpp>
#include<ImLib3D/ConvenienceProcessors.hpp>

void
IP3D::ShortestPathExtraction(const Mask3D &shape,const Vect3Di &start,const Vect3Di &end,vector<Vect3Di> &path,
                       Mask3D *pres,StructureElement *psel,
                       const LabelImage3D *pssCode0,
                       ShortestPathExtractionHelpers::StopCBType stopCB)
{
    // optional ssCode image
    const LabelImage3D *pssCode=pssCode0;
    if(!pssCode)
    {
        LabelImage3D *tmp=new LabelImage3D();
        IP3D::SSCode(shape,end,1,2,3,*tmp,&start);
        pssCode=tmp;
    }
    const LabelImage3D &ssCode=*pssCode;
       
    path.push_back(start);
    vector<Vect3Di> sel;
    if(!psel)
    {
        sel.push_back(Vect3Di(-1,0,0));
        sel.push_back(Vect3Di( 1,0,0));
        sel.push_back(Vect3Di(0,-1,0));
        sel.push_back(Vect3Di(0, 1,0));
        sel.push_back(Vect3Di(0,0,-1));
        sel.push_back(Vect3Di(0,0, 1));
    }
    else
    {
        for(uint i=0;i<psel->size();i++){sel.push_back(psel->GetPoint(i));}
    }

    for(uint it=0;;it++)
    {
        Vect3Di p=path.back();
        int smin=-1;
        int imin=-1;
        for(uint i=0;i<sel.size();i++)
        {
            Vect3Di npos=p+sel[i];
            if(ssCode.IsInside(npos) && ssCode(npos)>=0 && (imin<0 || ssCode(npos)<smin))
            {
                smin=ssCode(npos);
                imin=i;
            }
        }
        p=p+sel[imin];
        path.push_back(p);
        if(stopCB.first && stopCB.first(stopCB.second,p)){p.Show("CShortestPathExtraction:stoping at:","\n");break;}
        if(smin==0){break;}
    }

    if(!pssCode0){delete pssCode;}
}
void
IP3D::VoxelCode(const Mask3D &shape,const Mask3D &seeds,int m0,int m1,int m2,LabelImage3D& res,const Vect3Di *stopPoint)
{
    int metric[3]={m0,m1,m2};
    Mask3D s=seeds;
    s.Normalize();
    s.Invert();
    LabelImage3D codes;
    codes=s;
    int inf=(codes.Width()+codes.Height()+codes.Depth())*10;
    codes*=inf;
    

    // Build f,e,v neighborhoods
    StructureElement fNeighbors(MORPHO_Cross7);
    StructureElement planar19(MORPHO_Planar19);
    StructureElement cubic27(MORPHO_Cubic27);
    vector<Vect3Di> ePoints;
    for(uint i=0;i<planar19.size();i++)
    {
        if(fNeighbors.FindPoint(planar19[i])<0)
        {
            ePoints.push_back(planar19[i]);
        }
    }
    vector<Vect3Di> vPoints;
    for(uint i=0;i<cubic27.size();i++)
    {
        if(planar19.FindPoint(cubic27[i])<0)
        {
            vPoints.push_back(cubic27[i]);
        }
    }
    StructureElement eNeighbors(ePoints);
    StructureElement vNeighbors(vPoints);

    fNeighbors.RemoveCenter();
//      eNeighbors.RemoveCenter();
//      vNeighbors.RemoveCenter();

    StructureElement *neighbors[3];
    neighbors[0]=&fNeighbors;
    neighbors[1]=&eNeighbors;
    neighbors[2]=&vNeighbors;

    // keeping track of voxels that changed
    vector<Vect3Di> changedPts;// pts that changed on last iteration
    changedPts.reserve(shape.GetNVoxels());
    // initialize with seed points
    for(LabelImage3D::iteratorXYZ p=codes.begin();p!=codes.end();++p)
    {
        if(*p!=inf){changedPts.push_back(p.Pos());}
    }
    // image containing iteration at which a point last changed
    // this is used to avoid adding several times the same point to the changed list 
    // on the same iteration
    LabelImage3D changedIm(shape.Size());
    changedIm.Fill(-1);

    // iterate 
    for(int it=0;;it++)
    {
        bool changed=false;
        vector<Vect3Di> newChangedPts;
        newChangedPts.reserve(shape.GetNVoxels());
        LabelImage3D newCodes=codes;
        printf("it:%3d changedPts:%4d\n",it,changedPts.size());
        // for each point that changed on the last iter, propagate codes
        for(uint pt=0;pt<changedPts.size();pt++)
        {
            const Vect3Di &pos=changedPts[pt];
            // for f,e,v neighborhoods
            for(int n=0;n<3;n++)
            {
                // replace neighbors  code by new code if necessary
                for(uint i=0;i<neighbors[n]->size();i++)
                {
                    Vect3Di npos=pos+neighbors[n]->GetPoint(i);
                    int nval=codes(pos)+metric[n];
                    if(shape.IsInside(npos) && shape(npos) && nval<newCodes(npos))
                    {
                        newCodes(npos)=nval;
                        if(changedIm(npos)!=it){changedIm(npos)=it;newChangedPts.push_back(npos);}
                        changed=true;
                    }
                }
            }
        }
        codes=newCodes;
        if(0)
        {
            Mask3D changedDisp;
            IP3D::SimpleThresholds2(changedIm,it,it,changedDisp);
            ImLib3DDisplay::Show(changedDisp);
            ImLib3DDisplay::Show(codes);
        }
//          codes.WriteToFile(SPrintf("codes_%02d.im3D",it));
        if(stopPoint && changedIm(*stopPoint)!=-1){break;}
        if(!changed){break;}
        changedPts=newChangedPts;
    }
    // change inf to -1
    for(LabelImage3D::iteratorXYZ p=codes.begin();p!=codes.end();++p)
    {
        if(*p==inf){*p=-1;}
    }
    res=codes;
}

PathGraph::Path::Path(PathGraph::Node *_first,PathGraph::Node *_last,const vector<Vect3Di> &_points):
    first(_first),last(_last),points(_points)
{
    first->point=points[0];
    last ->point=points.back();
    first->paths.push_back(this);
    last ->paths.push_back(this);
}
void
SkeletonExtraction::MedPathFind::ExtractPath(SkeletonExtraction *se)
{
    path.clear();
    Vect3Di start=cpath[0]->MedialPoint(cpath[0]->center);
    IP3D::ShortestPathExtraction(Mask3D(),start,Vect3Di(),path,NULL,NULL,&se->ssCode,
                                           IP3D::ShortestPathExtractionHelpers::StopCBType(&SkeletonExtraction::SpeStopCB,se));
    if(se->clustersIma(path.back())->stop>=0)
    {
        path.back().Show("ExtractPath: found stop at:","\n");
        id2=se->clustersIma(path.back())->stop;
    }
    else{   path.back().Show("ExtractPath: no stop at:","\n");}
    if(cpath.size()!=1){ThrowError("MedPathFind::ExtractPath Should only be start cluster in cpath");}
    for(uint i=1;i<path.size();i++)
    {
        cpath.push_back(se->clustersIma(path[i]));
    }
}
    
void
SkeletonExtraction::ExtractAllMedPaths(vector<MedPathFind> &medPaths)
{
    //init
    medPaths.clear();
    for(uint i=0;i<clusters.size();i++){clusters[i]->stop=-1;}
    // sorted LM clusters 
    for(uint i=0;i<clusters.size();i++)
    {
        if(clusters[i]->isLM){medPaths.push_back(MedPathFind(clusters[i]));}
    }
    sort(medPaths.begin(),medPaths.end());
    reverse(medPaths.begin(),medPaths.end());

    for(uint mi=0;mi<medPaths.size();mi++)
    {
//          MedPathFind &mpf=medPaths[i];
        printf("medpath %d:startcluster :%x val:%d\n",mi,(uint)medPaths[mi].cpath[0],medPaths[mi].cpath[0]->value);
        medPaths[mi].id=mi;
        medPaths[mi].ExtractPath(this);
        medPaths[mi].path[0].Show("first point:","\n");

        for(uint pti=0;pti<medPaths[mi].path.size();pti++)
        {
            Vect3Di pt=medPaths[mi].path[pti];
            CodeCluster *c=clustersIma(pt);
            Vect3Di mp=c->MedialPoint(pt);
            printf("medpath % 3d : cluster:%x pti:% 3d  ",mi,(uint)c,pti);
            pt.Show("pt:","     ");
            mp.Show("mp:","\n");

            // create new medpaths if this is a branching cluster
            vector<CodeCluster *> s;
            c->GetSuccessors(s);
            for(uint j=0;j<s.size();j++)
            {
                if(pti+1<medPaths[mi].path.size() && 
                   s[j]!=clustersIma(medPaths[mi].path[pti+1]))
                {
                    printf("branches %d\n",medPaths.size());
                    medPaths.push_back(MedPathFind(s[j]));
                    medPaths.back().id1=mi;
                    medPaths.back().parentCluster=c;
                }
            }
            medPaths[mi].path[pti]=mp;
            c->stop=mi;
        }
    }
}
void 
SkeletonExtraction::PathGraphFromMedPaths(vector<MedPathFind> &medPaths,PathGraph &pathGraph)
{

    ssCode.WriteToFile("ssCode.im3D");
    // display medPaths
    {
        Image3Df disp(clustersIma.Size());disp.Fill(-1);
        for(uint i=0;i<medPaths.size();i++)
        {
            MedPathFind &mpf=medPaths[i];
            for(uint j=0;j<mpf.path.size();j++)
            {
                disp(mpf.path[j])=i+j*.5/mpf.path.size();
            }
        }
        disp.WriteToFile("medpaths.im3D");
    }
    // add first point to branching medpaths
    for(uint i=0;i<medPaths.size();i++)
    {
        MedPathFind &mpf=medPaths[i];
        if(mpf.id1>=0)
        {
            mpf.cpath.insert(mpf.cpath.begin(),mpf.parentCluster);
            mpf. path.insert(mpf. path.begin(),mpf.parentCluster->MedialPoint(mpf.path[0]));
        }
    }

    // create empty branches for each point in each medpath
    for(uint i=0;i<medPaths.size();i++)
    {
        MedPathFind &mpf=medPaths[i];
        mpf.branches.assign(mpf.path.size(),vector<int>());
    }
    // mark branching points
    for(uint i=0;i<medPaths.size();i++)
    {
        MedPathFind &mpf=medPaths[i];
        if(mpf.id1>=0)
        {
            MedPathFind &other=medPaths[mpf.id1];
            other.branches[other.FindCluster(mpf.cpath[0])].push_back(i);
        }
        if(mpf.id2>=0)
        {
            MedPathFind &other=medPaths[mpf.id2];
            other.branches[other.FindCluster(mpf.cpath.back())].push_back(i);
        }
    }
    // display branching
    {
        for(uint i=0;i<medPaths.size();i++)
        {
            MedPathFind &mpf=medPaths[i];
            printf("medpath %3d : id1:%3d id2:%3d\n",i,mpf.id1,mpf.id2);
            for(uint j=0;j<mpf.path.size();j++)
            {
                printf("cluster %x, ",(uint)mpf.cpath[j]);
                mpf.path[j].Show("pt:"," ");
                if(mpf.branches[j].size()>0)
                {
                    printf("branches (%d):",mpf.branches[j].size());
                    for(uint k=0;k<mpf.branches[j].size();k++){printf("%d ",mpf.branches[j][k]);}
                }
                printf("\n");
            }
        }
    }
    // create nodes
    for(uint i=0;i<medPaths.size();i++)
    {
        MedPathFind &mpf=medPaths[i];
        mpf.branchNodes.assign(mpf.path.size(),NULL);
        if(mpf.id1<0){mpf.branchNodes[0]=pathGraph.AddNode();}
        for(uint j=1;j<mpf.path.size()-1;j++)
        {
            if(mpf.branches[j].size()>0)
            {
                mpf.branchNodes[j]=pathGraph.AddNode();
            }
        }
        if(mpf.id2<0){mpf.branchNodes[mpf.path.size()-1]=pathGraph.AddNode();}
    }
    // merge nodes, until all nodes set
    bool changed=false;
    while(!changed)
    {
        printf("merging\n");
        for(uint i=0;i<medPaths.size();i++)
        {
            MedPathFind &mpf=medPaths[i];
            if(mpf.id1>=0 && !mpf.branchNodes[0])
            {
                MedPathFind &other=medPaths[mpf.id1];
                int opos=other.FindCluster(mpf.cpath[0]);
                mpf.branchNodes[0]=other.branchNodes[opos];
                changed=true;
            }
            int l=mpf.path.size()-1;
            if(mpf.id2>=0 && !mpf.branchNodes[l])
            {
                MedPathFind &other=medPaths[mpf.id2];
                int opos=other.FindCluster(mpf.cpath[l]);
                mpf.branchNodes[l]=other.branchNodes[opos];
                changed=true;
            }
        }
    }
    // create paths
    for(uint i=0;i<medPaths.size();i++)
    {
        MedPathFind &mpf=medPaths[i];
        PathGraph::Node *first=mpf.branchNodes[0];
        int j0=0;
        for(uint j=1;j<mpf.path.size();j++)
        {
            PathGraph::Node *node=mpf.branchNodes[j];
            if(node)
            {
                pathGraph.AddPath(first,node,vector<Vect3Di>(mpf.path.begin()+j0,mpf.path.begin()+j+1));
                j0=j;
                first=node;
            }
        }
    }
    // reunite node points
    for(uint i=0;i<pathGraph.paths.size();i++)
    {
        pathGraph.paths[i]->points[0]    =pathGraph.paths[i]->first->point;
        pathGraph.paths[i]->points.back()=pathGraph.paths[i]->last ->point;
    }
    // display paths
    {
        Image3Df disp(clustersIma.Size());disp.Fill(-1);
        for(uint i=0;i<pathGraph.paths.size();i++)
        {
            vector<Vect3Di> &pts=pathGraph.paths[i]->points;
            for(uint j=0;j<pts.size();j++)
            {
                disp(pts[j])=i+j*.5/pts.size();
            }
        }
        disp.WriteToFile("paths.im3D");
    }

}

void 
SkeletonExtraction::ExtractPaths(PathGraph &pathGraph)
{
    vector<MedPathFind> medPaths;
    ExtractAllMedPaths(medPaths);
    PathGraphFromMedPaths(medPaths,pathGraph);
    // fill in between points too far from each other
    for(uint i=0;i<pathGraph.paths.size();i++)
    {
        vector<Vect3Di> &pts=pathGraph.paths[i]->points;
        for(uint j=1;j<pts.size();j++)
        {
            Vect3Di dd=pts[j-1]-pts[j];
            if(abs(dd.x)>1 || abs(dd.y)>1 || abs(dd.z)>1)
            {
                vector<Vect3Di> inspts;
                ExtractLinkingMedPath(shape,pts[j-1],pts[j],inspts);
                if(inspts.size()>2)
                {
                    pts.insert(pts.begin()+j,inspts.begin()+1,inspts.end()-1);
                    j+=inspts.size()-2;
                }
            }
        }
    }
}

void 
SkeletonExtraction::CreateClusters()
{
    clusters.clear();
    clustersIma.Resize(ssCode.Size());
    clustersIma.Fill(NULL);

    LabelImage3D labels;
    int bg=-1;
    IP3D::ConnectedComponentLabelling(ssCode,StructureElement(MORPHO_Cubic27),labels,&bg);
//          ImLib3DDisplay::Show(labels);
    // for each unique label, create a new cluster  
    map<int,CodeCluster *> label2cluster;// used for keeping track of unique labels
    for(LabelImage3D::iteratorXYZ p=labels.begin();p!=labels.end();++p)
    {
        if(ssCode(p.Pos())<0){continue;}
        int label=*p;
        CodeCluster *cluster;
        // if we haven't seen this label yet, create a new cluster
        if(label2cluster.find(label)==label2cluster.end())
        {
            printf("new cluster:label %d value:%d\n",label,ssCode(p.Pos()));
            cluster=new CodeCluster(ssCode(p.Pos()));
            label2cluster[label]=cluster;
            clusters.push_back(cluster);
        }
        else{cluster=label2cluster[label];}
                
        clustersIma(p.Pos())=cluster;
        cluster->points.push_back(p.Pos());
    }

}
void 
SkeletonExtraction::InitClusters()
{
    for(uint i=0;i<clusters.size();i++)
    {
        clusters[i]->Init();
    }
    for(uint i=0;i<clusters.size();i++)
    {
        clusters[i]->FindNeighbors(clustersIma);
    }
    for(uint i=0;i<clusters.size();i++)
    {
        clusters[i]->FindType();
    }
    for(uint i=0;i<clusters.size();i++)
    {
        clusters[i]->FindMedialPoints(bsCode);
    }
}
SkeletonExtraction::SkeletonExtraction(const Mask3D &_shape,const Vect3Di *psourcePoint):shape(_shape)
{
    Vect3Di sourcePoint(-1,0,0);
    if(psourcePoint){sourcePoint=*psourcePoint;}
    else
    {
        // automatically find seed point
        for(Mask3D::const_iteratorXYZ p=shape.begin();p!=shape.end();++p){if(*p){sourcePoint=p.Pos();break;}}
        sourcePoint.Show("SkeletonExtraction: p0:","\n");
        if(sourcePoint.x<0){ThrowError("SkeletonExtraction::SkeletonExtraction shape is empty");}
        LabelImage3D ss0;
        IP3D::SSCode(shape,sourcePoint,1,2,3,ss0);
        Vect3Di dud;
        int dud0,dud1;
        IP3D::FindMinMax(ss0,dud0,dud1,&dud,&sourcePoint);
        sourcePoint.Show("SkeletonExtraction: auto sourcePoint:","\n");
    }
    IP3D::SSCode(shape,sourcePoint,1,2,3,ssCode);
    IP3D::BSCode(shape,1,2,3,bsCode);
    CreateClusters();
    InitClusters();
}
void 
SkeletonExtraction::ExtractAll(const Mask3D &shape,PathGraph &pathGraph)
{
    LabelImage3D components;
    IP3D::ConnectedComponentLabelling(shape,StructureElement(MORPHO_Cross7),components);
    int cmin,cmax;
    IP3D::FindMinMax(components,cmin,cmax);

    for(int i=1;i<=cmax;i++)
    {
        printf("SkeletonExtraction::ExtractAll indep component %d out of %d\n",i,cmax);
        Mask3D simpleShape;
        IP3D::SimpleThresholds2(components,i,i,simpleShape);
        SkeletonExtraction skel(simpleShape);
        skel.ExtractPaths(pathGraph);
    }
}
RectZone3Di
PathGraph::BoundingBox() const
{
    Vect3Di p0=paths[0]->points[0],p1=p0;
    for(uint i=0;i<paths.size();i++)
    {
        for(uint j=0;j<paths[i]->points.size();j++)
        {
            p0=min(p0,paths[i]->points[j]);
            p1=max(p1,paths[i]->points[j]);
        }
    }
    return RectZone3Di(p0,p1);
}
void
SkeletonExtraction::DrawPathGraph(const PathGraph &pathGraph,Mask3D &res)
{
    Vect3Di p0(0,0,0);
    if(res.Size()==Size3D(0,0,0))
    {
        RectZone3Di box=pathGraph.BoundingBox();
        res.Resize(Size3D(box.Width(),box.Height(),box.Depth()));
        p0=box.GetP0();
        res.properties.iCenter=new Vect3Di(p0);
    }
    for(uint i=0;i<pathGraph.paths.size();i++)
    {
        for(uint j=0;j<pathGraph.paths[i]->points.size();j++)
        {
            if(!res.IsInside(pathGraph.paths[i]->points[j]-p0))
            {
                ThrowError("SkeletonExtraction::DrawPathGraph point outside of image");
            }
            res(pathGraph.paths[i]->points[j]-p0)=1;
        }
    }
}
void 
SkeletonExtraction::ExtractLinkingMedPath(const Mask3D &shape,const Vect3Di &start,const Vect3Di &end,vector<Vect3Di> &res)
{
    if(start==end){res.push_back(start);res.push_back(end);return;}
    if(0)
    {// draw a line between both points
        int jmx=0,dmx=abs(end(0)-start(0));
        for(int j=1;j<3;j++)
        {
            int dd=abs(end(j)-start(j));
            if(dd>dmx){dmx=dd;jmx=j;}
        }
        int j=jmx;
        Vect3Di p0=start,p1=end;
        if(p1(j)<p0(j)){p1=start;p0=end;}
        for(int i=p0(j);i<=p1(j);i++)
        {
            Vect3Df p=(Vect3Df)p0+((i-p0(j))/(float)(p1(j)-p0(j)))*((Vect3Df)(p1-p0));
            res.push_back(floor(p+Vect3Df(.5,.5,.5)));
        }
    }
//      SkeletonExtraction skel(shape,&end);
//      vector<Vect3Di> path;
    res.clear();
    IP3D::ShortestPathExtraction(shape,start,end,res);
//      skel.MedPathFromPath(path,res);
}

---