Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members
ImLib3D

MorphologicalOperators.cpp

Go to the documentation of this file.
00001 /* ImLib3D
00002  * Copyright (c) 2001, ULP-IPB Strasbourg.
00003  *
00004  * This program is free software; you can redistribute it and/or modify
00005  * it under the terms of the GNU General Public License as published by
00006  * the Free Software Foundation; either version 2 of the License, or (at
00007  * your option) any later version.
00008  * 
00009  * This program is distributed in the hope that it will be useful, but
00010  * WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * General Public License for more details.
00013  * 
00014  * You should have received a copy of the GNU General Public License
00015  * along with this program; if not, write to the Free Software
00016  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00017  */
00019 
00020 
00021 #include<ImLib3D/Image3D.hpp>
00022 #include<ImLib3D/Arithmetic.hpp>
00023 #include<ImLib3D/MorphologicalOperators.hpp>
00024 #include<ImLib3D/Threshold.hpp>
00025 #include<ImLib3D/ConvenienceProcessors.hpp>
00026 
00027 
00028 
00029 void
00030 StructureElement::CreateSpherePointList(int radius,vector<Vect3Di>& pl)
00031 {
00032     int i,j,k;
00033     for (i=-radius;i<=radius;i++)
00034         for (j=-radius;j<=radius;j++)
00035             for (k=-radius;k<=radius;k++)
00036             {
00037                 Vect3Di vect(i,j,k);
00038                 if (vect.Norm() <= radius)
00039                     pl.push_back(vect);
00040             }
00041 }
00042 
00043 void StructureElement::ComputeCoeffs()
00044 {
00045     coeffs.clear();
00046     coeffs.assign(maskPoints.size(),0);
00047     for(uint i=0;i<maskPoints.size();i++)
00048     {
00049         coeffs[i]=GetPoint(i).Norm();
00050     }
00051 }
00052 
00053 void 
00054 StructureElement::RemoveCenter()
00055 {
00056     (*mask)(*mask->properties.iCenter)=0;
00057     int center=FindPoint(Vect3Di(0,0,0));
00058     if(center<0){ThrowError("StructureElement::RemoveCenter cant find center??");}
00059     maskPoints.erase(maskPoints.begin()+center);
00060 }
00061 void 
00062 StructureElement::Show(const string &name) const
00063 {
00064     mprintf("StructureElement %s\n",name.c_str());
00065     uint i;
00066     mprintf("maskPoints: ");
00067     for(i=0;i<maskPoints.size();i++)
00068     {
00069         mprintf("%3d %3d %3d : ",maskPoints[i].x,maskPoints[i].y,maskPoints[i].z);
00070     }
00071     mprintf("\n");
00072     mprintf("coeffs: ");
00073     for(i=0;i<coeffs.size();i++)
00074     {
00075         mprintf("%f ",coeffs[i]);
00076     }
00077     mprintf("\n");
00078     mprintf("Mask: ");
00079     mask->properties.iCenter->Show("Mask center:","\n");
00080     Mask3D &m=*mask;
00081     for(int y=0;y<m.Height();y++)
00082     {
00083         for(int z=0;z<m.Depth();z++)
00084         {
00085             for(int x=0;x<m.Width();x++)
00086             {
00087                 mprintf("%c",(m(x,y,z) ? 'X': '.'));
00088             }
00089             mprintf("  ");
00090         }
00091         mprintf("\n");      
00092     }
00093 
00094 }
00095 
00096 
00097 StructureElement::StructureElement(const StructureElement &other,StructureElementType direction):
00098     maskType(other.maskType),
00099     mask(NULL)
00100 {
00101     
00102     switch(direction)
00103     {
00104     case MORPHO_ALLDIRECTIONS:
00105     {
00106         maskPoints.insert(maskPoints.begin(),other.maskPoints.begin(),other.maskPoints.end());
00107         ComputeCoeffs();
00108         CreateMask();
00109 //          Show("MORPHO_ALLDIRECTIONS");
00110     }
00111     break;
00112     case MORPHO_FORWARD:
00113     {
00114         StructureElement full(other,MORPHO_ALLDIRECTIONS);
00115         mask=new Mask3D(full.Mask());
00116         coeffs.clear();
00117         bool done=false;
00118         Vect3Di C=*(mask->properties.iCenter);
00119         for(Mask3D::iteratorXYZ p=(mask->begin());p!=mask->end();++p)
00120         {
00121             if(done){*p=0;continue;}
00122             if(p.Pos()==C){*p=0;done=true;continue;}
00123             if(*p){coeffs.push_back(full.GetCoeff(p.Pos()-C));}
00124         }
00125         CreateMaskPoints();
00126 //          Show("MORPHO_FORWARD");
00127     }
00128     break;
00129     case MORPHO_BACKWARD:
00130     {
00131         StructureElement full(other,MORPHO_ALLDIRECTIONS);
00132         mask=new Mask3D(full.Mask());
00133         Vect3Di C=*(mask->properties.iCenter);
00134         coeffs.clear();
00135         bool done=false;
00136 //           p=(mask->begin());p+=mask->GetNVoxels()-1;
00137         Mask3D::iteratorXYZ p=mask->end();
00138         do
00139         {
00140             p--;
00141             if(done){*p=0;continue;}
00142             if(p.Pos()==C){*p=0;done=true;continue;}
00143             if(*p){coeffs.insert(coeffs.begin(),full.GetCoeff(p.Pos()-C));}
00144         }
00145         while(p!=mask->begin());
00146         CreateMaskPoints();
00147 //          Show("MORPHO_BACKWARD");
00148     }   
00149     break;
00150     default:
00151         ThrowError("StructureElement::StructureElement unknown direction:%d",direction);
00152     }
00153     
00154 }
00155 StructureElement::StructureElement(const string &s):
00156     mask(NULL)
00157 {
00158     StructureElementType typ=MORPHO_CUSTOM;
00159     if(s=="CUSTOM")typ=             MORPHO_CUSTOM;                
00160     if(s=="FORWARD")typ=            MORPHO_FORWARD;               
00161     if(s=="BACKWARD")typ=           MORPHO_BACKWARD;              
00162     if(s=="ALLDIRECTIONS")typ=      MORPHO_ALLDIRECTIONS;         
00163     if(s=="Cubic27")typ=            MORPHO_Cubic27;               
00164     if(s=="Cross7")typ=             MORPHO_Cross7;                
00165     if(s=="Planar19")typ=           MORPHO_Planar19;              
00166     if(s=="Corner9")typ=            MORPHO_Corner9;              
00167     if(s=="Pts81")typ=              MORPHO_Pts81;                 
00168     if(s=="Pts14")typ=              MORPHO_Pts14;                 
00169     if(s=="Pts4")typ=               MORPHO_Pts4;                  
00170     if(s=="Pts10")typ=              MORPHO_Pts10;                 
00171     if(s=="Pts5")typ=               MORPHO_Pts5;                  
00172     if(s=="Pts41")typ=              MORPHO_Pts41;                    
00173     if(s=="XYCross5")typ=           MORPHO_XYCross5;                    
00174     if(s=="XZCross5")typ=           MORPHO_XZCross5;                    
00175     if(s=="YZCross5")typ=           MORPHO_YZCross5;                    
00176     if(s=="SphereR3")typ=           MORPHO_SphereR3;                    
00177     if(s=="SphereR4")typ=           MORPHO_SphereR4;                    
00178     if(typ==MORPHO_CUSTOM){ThrowError("StructureElement::StructureElement: bad StructureElementType:%s",s);}
00179     maskType=typ;
00180     CreateMaskPoints(maskType);
00181     CreateMask();       
00182 }
00183 void
00184 StructureElement::ComputeMaskPointsSize(Vect3Di &minv,Vect3Di &maxv) const
00185 {
00186     // now compute maximum and minimum excursion of mask
00187     minv=maskPoints[0];
00188     maxv=maskPoints[0];
00189     for(uint i=1;i<maskPoints.size();i++)
00190     {
00191         minv=min(minv,maskPoints[i]);
00192         maxv=max(maxv,maskPoints[i]);
00193     }
00194 }
00195 void
00196 StructureElement::CreateMask()
00197 {
00198     Vect3Di minv,maxv;
00199     ComputeMaskPointsSize(minv,maxv);
00200     Vect3Di size=maxv-minv+Vect3Di(1,1,1);
00201 //      size.Show("StructureElement::CreateMask size:","\n");
00202 //      minv.Show("minv",":");
00203 //      maxv.Show("maxv","\n");
00204 
00205     mask=new Mask3D(Size3D(size));
00206     mask->Fill(0);
00207     for(uint i=0;i<maskPoints.size();i++)
00208     {
00209         (*mask)(maskPoints[i]-minv)=1;
00210 //          (maskPoints[i]-minv).Show("maskpos:","\n");
00211     }
00212     mask->properties.iCenter=new Vect3Di(-minv);
00213 }
00214 
00215 void
00216 StructureElement::CreateMaskPoints()
00217 {
00218     for(Mask3D::iteratorXYZ p=mask->begin();p!=mask->end();++p)
00219     {
00220         if(*p){maskPoints.push_back(p.Pos()-(*mask->properties.iCenter));}
00221     }
00222 }
00223 /*******************************************
00224 ** --  masque_3d() ----------------------------
00225 **
00227 ********************************************/
00228 void
00229 StructureElement::CreateMaskPoints(StructureElementType maskType)
00230 {    
00231     switch (maskType) 
00232     {     
00233     case MORPHO_Cubic27 :                   
00234         /* 26 elements dans le voisinage */
00235         maskPoints.push_back(Vect3Di( 0, 0, 0));
00236         maskPoints.push_back(Vect3Di(-1,-1,-1));
00237         maskPoints.push_back(Vect3Di( 0,-1,-1));
00238         maskPoints.push_back(Vect3Di( 1,-1,-1));
00239         maskPoints.push_back(Vect3Di(-1, 0,-1));
00240         maskPoints.push_back(Vect3Di( 0, 0,-1));
00241         maskPoints.push_back(Vect3Di( 1, 0,-1));
00242         maskPoints.push_back(Vect3Di(-1, 1,-1));
00243         maskPoints.push_back(Vect3Di( 0, 1,-1));
00244         maskPoints.push_back(Vect3Di( 1, 1,-1));
00245         maskPoints.push_back(Vect3Di(-1,-1, 0));
00246         maskPoints.push_back(Vect3Di( 0,-1, 0));
00247         maskPoints.push_back(Vect3Di( 1,-1, 0));
00248         maskPoints.push_back(Vect3Di(-1, 0, 0));
00249         maskPoints.push_back(Vect3Di( 1, 0, 0));
00250         maskPoints.push_back(Vect3Di(-1, 1, 0));
00251         maskPoints.push_back(Vect3Di( 0, 1, 0));
00252         maskPoints.push_back(Vect3Di( 1, 1, 0));
00253         maskPoints.push_back(Vect3Di(-1,-1, 1));
00254         maskPoints.push_back(Vect3Di( 0,-1, 1));
00255         maskPoints.push_back(Vect3Di( 1,-1, 1));
00256         maskPoints.push_back(Vect3Di(-1, 0, 1));
00257         maskPoints.push_back(Vect3Di( 0, 0, 1));
00258         maskPoints.push_back(Vect3Di( 1, 0, 1));
00259         maskPoints.push_back(Vect3Di(-1, 1, 1));
00260         maskPoints.push_back(Vect3Di( 0, 1, 1));
00261         maskPoints.push_back(Vect3Di( 1, 1, 1));
00262         break;                       
00263         
00264     case MORPHO_Cross7 :                /* cas du voisinage croix */
00265         /* 6 elements dans le voisinage */
00266         maskPoints.push_back(Vect3Di( 0, 0, 0));
00267         maskPoints.push_back(Vect3Di( 0,-1, 0));
00268         maskPoints.push_back(Vect3Di(-1, 0, 0));
00269         maskPoints.push_back(Vect3Di( 0, 0,-1));
00270         maskPoints.push_back(Vect3Di( 0, 0, 1));
00271         maskPoints.push_back(Vect3Di( 1, 0, 0));
00272         maskPoints.push_back(Vect3Di( 0, 1, 0));        
00273         break;
00274     
00275     case MORPHO_Planar19 :                  /* cas du voisinage plans */
00276         /* 18 elements dans le voisinage */
00277         maskPoints.push_back(Vect3Di( 0, 0, 0));
00278         maskPoints.push_back(Vect3Di(-1,-1, 0));
00279         maskPoints.push_back(Vect3Di( 0,-1,-1));
00280         maskPoints.push_back(Vect3Di( 0,-1, 0));
00281         maskPoints.push_back(Vect3Di( 0,-1, 1));
00282         maskPoints.push_back(Vect3Di( 1,-1, 0));
00283         maskPoints.push_back(Vect3Di(-1, 0,-1));
00284         maskPoints.push_back(Vect3Di(-1, 0, 0));
00285         maskPoints.push_back(Vect3Di(-1, 0, 1));
00286         maskPoints.push_back(Vect3Di( 0, 0,-1));
00287         maskPoints.push_back(Vect3Di( 0, 0, 1));
00288         maskPoints.push_back(Vect3Di( 1, 0,-1));
00289         maskPoints.push_back(Vect3Di( 1, 0, 0));
00290         maskPoints.push_back(Vect3Di( 1, 0, 1));
00291         maskPoints.push_back(Vect3Di(-1, 1, 0));
00292         maskPoints.push_back(Vect3Di( 0, 1,-1));
00293         maskPoints.push_back(Vect3Di( 0, 1, 0));
00294         maskPoints.push_back(Vect3Di( 0, 1, 1));
00295         maskPoints.push_back(Vect3Di( 1, 1, 0));
00296         
00297         break;
00298         
00299     case MORPHO_Corner9 :               /* cas du voisinage coin */
00300         /* 8 elements dans le voisinage */
00301         maskPoints.push_back(Vect3Di( 0, 0, 0));
00302         maskPoints.push_back(Vect3Di(-1,-1,-1));
00303         maskPoints.push_back(Vect3Di(-1,-1, 1));
00304         maskPoints.push_back(Vect3Di( 1,-1,-1));
00305         maskPoints.push_back(Vect3Di( 1,-1, 1));
00306         maskPoints.push_back(Vect3Di(-1, 1,-1));
00307         maskPoints.push_back(Vect3Di(-1, 1, 1));
00308         maskPoints.push_back(Vect3Di( 1, 1,-1));
00309         maskPoints.push_back(Vect3Di( 1, 1, 1));
00310         
00311         break;
00312             
00313     case MORPHO_Pts81 :
00314         maskPoints.push_back(Vect3Di( 0, 0, 0));
00315         maskPoints.push_back(Vect3Di(-1,-2,-1));
00316         maskPoints.push_back(Vect3Di(-1,-2, 0));
00317         maskPoints.push_back(Vect3Di(-1,-2, 1));
00318         maskPoints.push_back(Vect3Di( 0,-2,-1));
00319         maskPoints.push_back(Vect3Di( 0,-2, 0));
00320         maskPoints.push_back(Vect3Di( 0,-2, 1));
00321         maskPoints.push_back(Vect3Di( 1,-2,-1));
00322         maskPoints.push_back(Vect3Di( 1,-2, 0));
00323         maskPoints.push_back(Vect3Di( 1,-2, 1));
00324         maskPoints.push_back(Vect3Di(-2,-1,-1));
00325         maskPoints.push_back(Vect3Di(-2,-1, 0));
00326         maskPoints.push_back(Vect3Di(-2,-1, 1));
00327         maskPoints.push_back(Vect3Di(-1,-1,-2));
00328         maskPoints.push_back(Vect3Di(-1,-1,-1));
00329         maskPoints.push_back(Vect3Di(-1,-1, 0));
00330         maskPoints.push_back(Vect3Di(-1,-1, 1));
00331         maskPoints.push_back(Vect3Di(-1,-1, 2));
00332         maskPoints.push_back(Vect3Di( 0,-1,-2));
00333         maskPoints.push_back(Vect3Di( 0,-1,-1));
00334         maskPoints.push_back(Vect3Di( 0,-1, 0));
00335         maskPoints.push_back(Vect3Di( 0,-1, 1));
00336         maskPoints.push_back(Vect3Di( 0,-1, 2));
00337         maskPoints.push_back(Vect3Di( 1,-1,-2));
00338         maskPoints.push_back(Vect3Di( 1,-1,-1));
00339         maskPoints.push_back(Vect3Di( 1,-1, 0));
00340         maskPoints.push_back(Vect3Di( 1,-1, 1));
00341         maskPoints.push_back(Vect3Di( 1,-1, 2));
00342         maskPoints.push_back(Vect3Di( 2,-1,-1));
00343         maskPoints.push_back(Vect3Di( 2,-1, 0));
00344         maskPoints.push_back(Vect3Di( 2,-1, 1));
00345         maskPoints.push_back(Vect3Di(-2, 0,-1));
00346         maskPoints.push_back(Vect3Di(-2, 0, 0));
00347         maskPoints.push_back(Vect3Di(-2, 0, 1));
00348         maskPoints.push_back(Vect3Di(-1, 0,-2));
00349         maskPoints.push_back(Vect3Di(-1, 0,-1));
00350         maskPoints.push_back(Vect3Di(-1, 0, 0));
00351         maskPoints.push_back(Vect3Di(-1, 0, 1));
00352         maskPoints.push_back(Vect3Di(-1, 0, 2));
00353         maskPoints.push_back(Vect3Di( 0, 0,-2));
00354         maskPoints.push_back(Vect3Di( 0, 0,-1));
00355         maskPoints.push_back(Vect3Di( 0, 0, 1));
00356         maskPoints.push_back(Vect3Di( 0, 0, 2));
00357         maskPoints.push_back(Vect3Di( 1, 0,-2));
00358         maskPoints.push_back(Vect3Di( 1, 0,-1));
00359         maskPoints.push_back(Vect3Di( 1, 0, 0));
00360         maskPoints.push_back(Vect3Di( 1, 0, 1));
00361         maskPoints.push_back(Vect3Di( 1, 0, 2));
00362         maskPoints.push_back(Vect3Di( 2, 0,-1));
00363         maskPoints.push_back(Vect3Di( 2, 0, 0));
00364         maskPoints.push_back(Vect3Di( 2, 0, 1));
00365         maskPoints.push_back(Vect3Di(-2, 1,-1));
00366         maskPoints.push_back(Vect3Di(-2, 1, 0));
00367         maskPoints.push_back(Vect3Di(-2, 1, 1));
00368         maskPoints.push_back(Vect3Di(-1, 1,-2));
00369         maskPoints.push_back(Vect3Di(-1, 1,-1));
00370         maskPoints.push_back(Vect3Di(-1, 1, 0));
00371         maskPoints.push_back(Vect3Di(-1, 1, 1));
00372         maskPoints.push_back(Vect3Di(-1, 1, 2));
00373         maskPoints.push_back(Vect3Di( 0, 1,-2));
00374         maskPoints.push_back(Vect3Di( 0, 1,-1));
00375         maskPoints.push_back(Vect3Di( 0, 1, 0));
00376         maskPoints.push_back(Vect3Di( 0, 1, 1));
00377         maskPoints.push_back(Vect3Di( 0, 1, 2));
00378         maskPoints.push_back(Vect3Di( 1, 1,-2));
00379         maskPoints.push_back(Vect3Di( 1, 1,-1));
00380         maskPoints.push_back(Vect3Di( 1, 1, 0));
00381         maskPoints.push_back(Vect3Di( 1, 1, 1));
00382         maskPoints.push_back(Vect3Di( 1, 1, 2));
00383         maskPoints.push_back(Vect3Di( 2, 1,-1));
00384         maskPoints.push_back(Vect3Di( 2, 1, 0));
00385         maskPoints.push_back(Vect3Di( 2, 1, 1));
00386         maskPoints.push_back(Vect3Di(-1, 2,-1));
00387         maskPoints.push_back(Vect3Di(-1, 2, 0));
00388         maskPoints.push_back(Vect3Di(-1, 2, 1));
00389         maskPoints.push_back(Vect3Di( 0, 2,-1));
00390         maskPoints.push_back(Vect3Di( 0, 2, 0));
00391         maskPoints.push_back(Vect3Di( 0, 2, 1));
00392         maskPoints.push_back(Vect3Di( 1, 2,-1));
00393         maskPoints.push_back(Vect3Di( 1, 2, 0));
00394         maskPoints.push_back(Vect3Di( 1, 2, 1));
00395         break;
00396             
00397     case MORPHO_Pts14 :
00398         maskPoints.push_back(Vect3Di( 0, 0, 0));
00399         maskPoints.push_back(Vect3Di(-1,-1,-1));
00400         maskPoints.push_back(Vect3Di(-1,-1, 0));
00401         maskPoints.push_back(Vect3Di(-1,-1, 1));
00402         maskPoints.push_back(Vect3Di( 0,-1,-1));
00403         maskPoints.push_back(Vect3Di( 0,-1, 0));
00404         maskPoints.push_back(Vect3Di( 0,-1, 1));
00405         maskPoints.push_back(Vect3Di( 1,-1,-1));
00406         maskPoints.push_back(Vect3Di( 1,-1, 0));
00407         maskPoints.push_back(Vect3Di( 1,-1, 1));
00408         maskPoints.push_back(Vect3Di(-1, 0,-1));
00409         maskPoints.push_back(Vect3Di(-1, 0, 0));
00410         maskPoints.push_back(Vect3Di(-1, 0, 1));
00411         maskPoints.push_back(Vect3Di( 0, 0,-1));
00412         break;
00413           
00414     case MORPHO_Pts4 :              
00415         maskPoints.push_back(Vect3Di( 0, 0, 0));
00416         maskPoints.push_back(Vect3Di( 0,-1, 0));
00417         maskPoints.push_back(Vect3Di(-1, 0, 0));
00418         maskPoints.push_back(Vect3Di( 0, 0,-1));
00419         break;
00420         
00421     case MORPHO_Pts10 :                 
00422         maskPoints.push_back(Vect3Di( 0, 0, 0));
00423         maskPoints.push_back(Vect3Di(-1,-1, 0));
00424         maskPoints.push_back(Vect3Di( 0,-1,-1));
00425         maskPoints.push_back(Vect3Di( 0,-1, 0));
00426         maskPoints.push_back(Vect3Di( 0,-1, 1));
00427         maskPoints.push_back(Vect3Di( 1,-1, 0));
00428         maskPoints.push_back(Vect3Di(-1, 0,-1));
00429         maskPoints.push_back(Vect3Di(-1, 0, 0));
00430         maskPoints.push_back(Vect3Di(-1, 0, 1));
00431         maskPoints.push_back(Vect3Di( 0, 0,-1));    
00432         break;
00433         
00434     case MORPHO_Pts5 :              
00435         maskPoints.push_back(Vect3Di( 0, 0, 0));
00436         maskPoints.push_back(Vect3Di(-1,-1,-1));
00437         maskPoints.push_back(Vect3Di(-1,-1, 1));
00438         maskPoints.push_back(Vect3Di( 1,-1,-1));
00439         maskPoints.push_back(Vect3Di( 1,-1, 1));
00440         break;
00441         
00442     case MORPHO_Pts41 :
00443         maskPoints.push_back(Vect3Di( 0, 0, 0));
00444         maskPoints.push_back(Vect3Di(-1,-2,-1));
00445         maskPoints.push_back(Vect3Di(-1,-2, 0));
00446         maskPoints.push_back(Vect3Di(-1,-2, 1));
00447         maskPoints.push_back(Vect3Di( 0,-2,-1));
00448         maskPoints.push_back(Vect3Di( 0,-2, 0));
00449         maskPoints.push_back(Vect3Di( 0,-2, 1));
00450         maskPoints.push_back(Vect3Di( 1,-2,-1));
00451         maskPoints.push_back(Vect3Di( 1,-2, 0));
00452         maskPoints.push_back(Vect3Di( 1,-2, 1));
00453         maskPoints.push_back(Vect3Di(-2,-1,-1));
00454         maskPoints.push_back(Vect3Di(-2,-1, 0));
00455         maskPoints.push_back(Vect3Di(-2,-1, 1));
00456         maskPoints.push_back(Vect3Di(-1,-1,-2));
00457         maskPoints.push_back(Vect3Di(-1,-1,-1));
00458         maskPoints.push_back(Vect3Di(-1,-1, 0));
00459         maskPoints.push_back(Vect3Di(-1,-1, 1));
00460         maskPoints.push_back(Vect3Di(-1,-1, 2));
00461         maskPoints.push_back(Vect3Di( 0,-1,-2));
00462         maskPoints.push_back(Vect3Di( 0,-1,-1));
00463         maskPoints.push_back(Vect3Di( 0,-1, 0));
00464         maskPoints.push_back(Vect3Di( 0,-1, 1));
00465         maskPoints.push_back(Vect3Di( 0,-1, 2));
00466         maskPoints.push_back(Vect3Di( 1,-1,-2));
00467         maskPoints.push_back(Vect3Di( 1,-1,-1));
00468         maskPoints.push_back(Vect3Di( 1,-1, 0));
00469         maskPoints.push_back(Vect3Di( 1,-1, 1));
00470         maskPoints.push_back(Vect3Di( 1,-1, 2));
00471         maskPoints.push_back(Vect3Di( 2,-1,-1));
00472         maskPoints.push_back(Vect3Di( 2,-1, 0));
00473         maskPoints.push_back(Vect3Di( 2,-1, 1));
00474         maskPoints.push_back(Vect3Di(-2, 0,-1));
00475         maskPoints.push_back(Vect3Di(-2, 0, 0));
00476         maskPoints.push_back(Vect3Di(-2, 0, 1));
00477         maskPoints.push_back(Vect3Di(-1, 0,-2));
00478         maskPoints.push_back(Vect3Di(-1, 0,-1));
00479         maskPoints.push_back(Vect3Di(-1, 0, 0));
00480         maskPoints.push_back(Vect3Di(-1, 0, 1));
00481         maskPoints.push_back(Vect3Di(-1, 0, 2));
00482         maskPoints.push_back(Vect3Di( 0, 0,-2));
00483         maskPoints.push_back(Vect3Di( 0, 0,-1));
00484         break;
00485     case MORPHO_XYCross5 :              
00486         maskPoints.push_back(Vect3Di( 0, 0, 0));
00487         maskPoints.push_back(Vect3Di( 1, 0, 0));
00488         maskPoints.push_back(Vect3Di(-1, 0, 0));
00489         maskPoints.push_back(Vect3Di( 0, 1, 0));
00490         maskPoints.push_back(Vect3Di( 0,-1, 0));
00491         break;
00492     case MORPHO_XZCross5 :              
00493         maskPoints.push_back(Vect3Di( 0, 0, 0));
00494         maskPoints.push_back(Vect3Di( 1, 0, 0));
00495         maskPoints.push_back(Vect3Di(-1, 0, 0));
00496         maskPoints.push_back(Vect3Di( 0, 0, 1));
00497         maskPoints.push_back(Vect3Di( 0, 0,-1));
00498         break;
00499     case MORPHO_YZCross5 :              
00500         maskPoints.push_back(Vect3Di( 0, 0, 0));
00501         maskPoints.push_back(Vect3Di( 0, 1, 0));
00502         maskPoints.push_back(Vect3Di( 0,-1, 0));
00503         maskPoints.push_back(Vect3Di( 0, 0, 1));
00504         maskPoints.push_back(Vect3Di( 0, 0,-1));
00505         break;
00506     case MORPHO_SphereR3 :              
00507         CreateSpherePointList(3,maskPoints);
00508         break;
00509     case MORPHO_SphereR4 :              
00510         CreateSpherePointList(4,maskPoints);
00511         break;
00512                 
00513     default:
00514         ThrowError("Unsupported StructureElementType:%d",maskType);
00515         break;
00516     }
00517 
00518 }
00519 
00520 
00521 
00522 #include<ImLib3D/Display.hpp>
00523 #include<ImLib3D/TestPatterns.hpp>
00524 
00525 void
00526 IP3D::DistanceTransform(const Mask3D& src, const StructureElement &mask, Image3Df& res,bool borderExtend)
00527 {
00528     res.Resize(src.Size());
00529     // initialize destiniation image
00530     // foreground=1 background=inf
00531     Mask3D::const_iteratorFast sp;
00532     Image3Df::iteratorFast rp;
00533     float BIGVALUE=1000000000.0;
00534     for(rp=res.begin(),sp=src.begin();rp!=res.end();++rp,sp++)
00535     {
00536         if(!*sp){*rp=0;}
00537         else{*rp=BIGVALUE;}
00538     }
00539     
00540     Vect3Di vmin(0,0,0);
00541     Vect3Di vmax=src.SizeV()-Vect3Di(1,1,1);
00542 
00543     // Forward pass
00544     {
00545         StructureElement forwardMask(mask,MORPHO_FORWARD);
00546         uint npts=forwardMask.size();
00547         ImageNeighbors<Image3Df> neighbors(res,forwardMask);
00548         Image3Df::iteratorXYZ rp;
00549         ImageProgress tracker("DistanceTransform-forward",rp);
00550         for(rp=res.begin();rp!=res.end();++rp)
00551         {
00552             if(!*rp){continue;}
00553             const Vect3Di &rpos=rp.Pos();
00554             if(neighbors.IsSafe(rpos))
00555             {
00556                 neighbors.MoveTo(rpos);
00557                 for(uint i=0;i<npts;i++)
00558                 {
00559                     *rp=min(*rp,neighbors[i]+forwardMask.GetCoeff(i));
00560                 }
00561             }
00562             else
00563             {
00564                 for(uint i=0;i<npts;i++)
00565                 {
00566                     Vect3Di pos=rp.Pos()+forwardMask[i];
00567                     if(!res.IsInside(pos))
00568                     {
00569                         if(borderExtend)
00570                         {
00571                             byte v=src(min(vmax,max(vmin,pos)));
00572                             if(!v){*rp=0;}
00573                             continue;
00574                         }
00575                         else
00576                         {*rp=0;break;}
00577                     }
00578                     else
00579                     {*rp=min(*rp,res(rp.Pos()+forwardMask[i])+forwardMask.GetCoeff(i));}
00580                 }
00581             }
00582         }   
00583     }
00584 
00585 //  res.WriteToFile("distanceTransformForwardpoass.im3D");
00586 
00587     // BackWard pass
00588     {
00589         StructureElement backwardMask(mask,MORPHO_BACKWARD);
00590         uint npts=backwardMask.size();
00591         ImageNeighbors<Image3Df> neighbors(res,backwardMask);
00592         Image3Df::iteratorXYZ rp=res.end();
00593         ImageProgress tracker("DistanceTransform-backward",rp);
00594         do
00595         {
00596             rp--;
00597             if(!*rp){continue;}
00598             const Vect3Di &rpos=rp.Pos();
00599             if(neighbors.IsSafe(rpos))
00600             {
00601                 neighbors.MoveBackwardsTo(rpos);
00602                 for(uint i=0;i<npts;i++)
00603                 {
00604                     *rp=min(*rp,neighbors[i]+backwardMask.GetCoeff(i));
00605                 }
00606             }
00607             else
00608             {
00609                 // is not inside safezone
00610                 for(uint i=0;i<backwardMask.size();i++)
00611                 {
00612                     Vect3Di pos=rp.Pos()+backwardMask[i];
00613                     if(!res.IsInside(pos))
00614                     {
00615                         // neighbor outside image
00616                         if(borderExtend)
00617                         {
00618                             byte v=src(min(vmax,max(vmin,pos)));
00619                             if(!v){*rp=0;}
00620                             continue;
00621                         }
00622                         else
00623                         {*rp=0;break;}
00624                     }
00625                     else// neighbor inside image 
00626                     {*rp=min(*rp,res(rp.Pos()+backwardMask[i])+backwardMask.GetCoeff(i));}                  
00627                 }
00628             }
00629         }
00630         while(rp!=res.begin());
00631     }
00632 }
00633 
00634 #ifdef NOTDEF
00635 void
00636 IP3D::ConnectedComponentLabelling(const Mask3D& src, const StructureElement &mask0, LabelImage3D& res)
00637 {
00638     // ALGO: For each pixel 
00639 
00640     vector<int> indx2region;
00641     indx2region.reserve(1000);
00642     // default indx=0 -> region=0 = background
00643     indx2region.push_back(0);
00644     int nbRegions=1;
00645 
00646     res=src;
00647 //      for(LabelImage3D::iteratorFast fp=res.begin();fp!=res.end();++fp){if(*fp){*fp=1;}}
00648 
00649     StructureElement mask(mask0,MORPHO_FORWARD);
00650     LabelImage3D::iteratorXYZ rp;
00651     vector<int> neigborIndexes;
00652     neigborIndexes.reserve(mask.size());
00653 
00654     for(rp=res.begin();rp!=res.end();++rp)
00655     {
00656         int talk=0;
00657         if(rp.Pos()==Vect3Di(48,5,76)){talk=1;}
00658 
00659         if(*rp)
00660         {// we found a non empty pixel
00661             // let's see if it has any neigbors
00662             neigborIndexes.clear();
00663             int currIndex;
00664             int prevIndex=-1;
00665             for(int i=0;i<mask.size();i++)
00666             {
00667                 if((currIndex=res.SafeValue(rp.Pos()+mask[i])))
00668                 {
00669                     // found a neigbor!
00670 //                      mprintf("found neigbor index:%2d ",currIndex);
00671 //                      (rp.Pos()+mask[i]).Show("pos:","\n");
00672                     if(prevIndex<0 || currIndex!=prevIndex)
00673                     {
00674                         neigborIndexes.push_back(currIndex);
00675                     }
00676                 }
00677             }
00678             if(neigborIndexes.size())
00679             {// found neigbors!
00680                 if(talk)
00681                 {
00682                     for(uint i=0;i<neigborIndexes.size();i++)
00683                     {
00684                         mprintf("found neigborIndexes:%d:%d\n",i,neigborIndexes[i]);
00685                     }
00686                 }
00687                 // set this pixel to the same value as one of those neigbors
00688                 int anIndexOfThisRegion=neigborIndexes[0];
00689                 *rp=anIndexOfThisRegion;
00690                 // make all neigbor indexes point to the same region
00691                 int thisRegion=indx2region[anIndexOfThisRegion];
00692                 if(talk)
00693                 {
00694                     mprintf("anIndexOfThisRegion:%d\n",anIndexOfThisRegion);
00695                     mprintf("thisRegion:%d\n",thisRegion);
00696                 }
00697                 for(uint i=1;i<neigborIndexes.size();i++)
00698                 {
00699                     if(talk)
00700                     {
00701                         mprintf("setting neigborIndexes:%d:%d: previous region:%d\n",i,neigborIndexes[i],indx2region[neigborIndexes[i]]);
00702                     }
00703                     int oldNeighborRegion=indx2region[neigborIndexes[i]];
00704                     if(oldNeighborRegion!=thisRegion)
00705                     {// merge 2 regions
00706                         for(uint j=0;j<indx2region.size();j++)
00707                         {
00708                             if(indx2region[j]==oldNeighborRegion)
00709                             {
00710                                 indx2region[j]=thisRegion;
00711                             }
00712                         }
00713                     }
00714                 }
00715             }
00716             else
00717             {// sorry, no neighbors 
00718                 // this lonesome pixel is a new region and index
00719                 *rp=indx2region.size();
00720                 int thisRegion=nbRegions++;
00721                 indx2region.push_back(thisRegion);
00722             }
00723         }else
00724         {
00725             *rp=0;// background index
00726         }
00727     }
00728 //      mprintf("CConnectedComponentLabelling:: region image\n");
00729 //      LabelImage3D regions=res;
00730 //      for(LabelImage3D::iteratorFast fp=regions.begin();fp!=regions.end();++fp){*fp=indx2region[*fp];}    
00731 //      ImLib3DDisplay::Show(regions);
00732 
00733     // ***** now compute new unique regions
00734     uint i;
00735     map<int,int> old2newregion;
00736     int nbNewRegions=0;
00737     // add  background region
00738     old2newregion[0]=nbNewRegions++;
00739     for(i=0;i<indx2region.size();i++)
00740     {
00741         int old=indx2region[i];
00742         if(old2newregion.find(old)==old2newregion.end())
00743         {
00744             old2newregion[old]=nbNewRegions++;
00745         }
00746     }
00747     // swithch to new regions
00748     for(i=0;i<indx2region.size();i++)
00749     {
00750         indx2region[i]=old2newregion[indx2region[i]];
00751     }   
00752 
00753 
00754     // ***** now lets switch from indexes to regions
00755     for(rp=res.begin();rp!=res.end();++rp)
00756     {
00757         *rp=indx2region[*rp];
00758     }
00759 }
00760 #endif// NOTDEF
00761 
00762 void
00763 IP3D::FillHoles(const Mask3D& src, const StructureElement &mask, Mask3D& res0)
00764 {
00765     ImageProcessorArgs<Mask3D> args(ImArgs(src),res0);
00766     args.SameSize();
00767     Mask3D &res=args.GetDest();
00768     
00769     // compute comonents connected to exterior
00770     Mask3D big(Size3D(src.SizeV()+Vect3Di(2,2,2)));
00771     big.Fill(0);
00772     IP3D::InsertImage(big,src,Vect3Di(1,1,1),big);
00773     big.Normalize();
00774 
00775     Mask3D inv(big.Size());
00776     inv.Fill(1);
00777     inv-=big;
00778     
00779 //      mprintf("CFillHoles:: before connected components: \n");
00780 //      ImLib3DDisplay::Show(inv);
00781 
00782     LabelImage3D lab;
00783     IP3D::ConnectedComponentLabelling(inv,mask,lab);
00784 //      mprintf("CFillHoles:: connected components: exterior=%d\n",lab(0,0,0));
00785 //      ImLib3DDisplay::Show(lab);
00786 
00787     // set to 1 connected connected to exterior
00788 //      int exterior=lab(0,0,0);
00789     // find which values are exterior (usefull for strange struct el)
00790     map<int,int> isExterior;
00791     for(LabelImage3D::iteratorXYZ ep=lab.begin();ep!=lab.end();++ep)
00792     {
00793         if(!src.IsInside(ep.Pos()-Vect3Di(1,1,1))){isExterior[*ep]=1;}
00794     }
00795     // now set res to 1 for non exterior values of lab
00796     RectZone3Di zone(Vect3Di(1,1,1),res.SizeV());
00797     LabelImage3D::iteratorZone llp(zone);
00798     Mask3D::iteratorFast rp;
00799     for(llp=lab.begin(),rp=res.begin();rp!=res.end();++rp,llp++)
00800     {
00801         if(isExterior.find(*llp)==isExterior.end()){*rp=1;}
00802         else{*rp=0;}
00803     }       
00804     mprintf("IP3D::FillHoles src:%x resimage:%x data:%x\n",(uint)&src,(uint)&res,(uint)res.GetData());
00805 }
00806 
00807 void
00808 IP3D::LargestConnectedComponent(const Mask3D &src,Mask3D &res0)
00809 {
00810     ImageProcessorArgs<Mask3D> args(ImArgs(src),res0);
00811     args.SameSize();
00812     Mask3D &res=args.GetDest();
00813 
00814     LabelImage3D imLabel;
00815     IP3D::ConnectedComponentLabelling(src,StructureElement(MORPHO_Cross7),imLabel);
00816     ImageValues<int> regions;
00817     IP3D::MakeValueList(imLabel,regions);
00818     // find bigest region
00819     int maxsize=-1;
00820     int maxval=-1;
00821     for(ImageValues<int>::MapType::iterator p=regions.values.begin();p!=regions.values.end();++p)
00822     {
00823         if((*p).second.size>maxsize && imLabel((*p).second.point)!=0){maxsize=(*p).second.size;maxval=(*p).first;}
00824     }
00825     printf("maxsize:%d maxval:%d\n",maxsize,maxval);
00826     IP3D::SimpleThresholds2(imLabel,maxval,maxval,res); 
00827 }
00828 
00829
Generated on Fri Jun 17 13:36:05 2005 for ImLib3D by  doxygen 1.4.2