X.Zhuang/zxhproj/zxh_morphology_8h_source.html

00001
00002 /*=========================================================================
00003
00004   Program:   ZXH Registration Software
00005   Author:    Xiahai Zhuang
00006   Module:    $RCSfle: zxhMorphology.h    $
00007   Language:  C++
00008   Date:      $Date: From  2007-07 $
00009   Version:   $Revision: 2.0 $
00010
00011 =========================================================================*/
00012
00013 #ifndef zxhMorphology_h
00014 #define zxhMorphology_h
00015
00016 #include <string>
00017 #include <iostream>
00018 #include <fstream>
00019 #include <stdio.h>
00020 #include <stdlib.h>
00021 #include "zxhImageData.h"
00022 #include "zxhImageModelingLinear.h"
00023 #include "zxhImageGipl.h"
00024
00031
00032 template <class PixelType=PixelTypeDefault>
00033 class zxhMorphologyT
00034 {
00035 public:
00037     zxhMorphologyT():m_BackGround(ZXH_Background),m_ForeGround(ZXH_Foreground),m_pImage(0){};
00039     zxhMorphologyT(PixelType background,PixelType foreground):m_BackGround(background),m_ForeGround(foreground),m_pImage(0){};
00041     virtual ~zxhMorphologyT(){};
00043     zxhMorphologyT<PixelType> * Clone(zxhMorphologyT<PixelType>*& pRet)
00044     {
00045         if(pRet==0) pRet=new zxhMorphologyT<PixelType>(m_BackGround,m_ForeGround);
00046         else pRet->SetBinaryValue(m_BackGround,m_ForeGround);
00047         pRet->SetImage(m_pImage);
00048         return pRet;
00049     };
00051     void    SetBinaryValue(PixelType background,PixelType foreground)
00052     {m_BackGround   = background;m_ForeGround   = foreground;}
00054     void    SetImage(zxhImageDataT<PixelType> * pImage)     {m_pImage=pImage;};
00056     zxhImageDataT<PixelType> * GetImage()           {return m_pImage;};
00057
00060     bool    Dilation(int steps)
00061     {Dilation3D(m_pImage,float(steps),m_BackGround,m_ForeGround);};
00065     static bool Dilation3D(zxhImageDataT<PixelType> *pImage, float steps, PixelType background );
00067     static bool DilateStep3D(zxhImageDataT<PixelType> *pImage, int steps,PixelType background, PixelType foreground=ZXH_Foreground );
00068     static bool ErodeStep3D( zxhImageDataT<PixelType> *pImage, int step, PixelType background, PixelType foreground=ZXH_Foreground );
00070     static bool Erosion3D(zxhImageDataT<PixelType> *pImage, float steps,PixelType background ) ;
00072     static bool Open3D(zxhImageDataT<PixelType> *pImage, float steps,PixelType background ) ;
00074     static bool Close3D(zxhImageDataT<PixelType> *pImage, float steps,PixelType background ) ;
00076     static bool MaskImageOpAnd3D( zxhImageDataT<PixelType>* pDes, zxhImageDataT<PixelType>*pSrc ) ;
00078     static bool EuclideanTransformBk3D(zxhImageDataT<PixelType>*pImage,zxhImageDataT<float>*pOut,PixelType background);
00079
00081     static bool Threshold3d( zxhImageDataT<PixelType>* pDes, PixelType RangeFrom, PixelType RangeTo, PixelType foreground, PixelType background );
00083     //static bool SqEuclideanTransformBk3D(zxhImageDataT<PixelType>*pImage,zxhImageDataT<float>*pOut,PixelType background);
00084
00085 protected:
00086     PixelType       m_BackGround;
00087     PixelType       m_ForeGround;
00088     zxhImageDataT<PixelType>    *m_pImage;
00089 };
00090
00091 typedef zxhMorphologyT<PixelTypeDefault> zxhMorphology;
00092
00093 // implementation -----------------------------------------------------------------------
00094 template <class PixelType>
00095 bool zxhMorphologyT<PixelType>::Threshold3d( zxhImageDataT<PixelType>* pDes, PixelType RangeFrom, PixelType RangeTo, PixelType foreground, PixelType background)
00096 {
00097     long int lsize = pDes->GetNumberOfPixels() ;
00098     for(long int ip=0; ip<lsize; ++ip )
00099         if( pDes->GetImageData( ip )>=RangeFrom && pDes->GetImageData(ip)<=RangeTo )
00100             pDes->SetImageData( ip, foreground ) ;
00101         else pDes->SetImageData( ip, background ) ;
00102     return true;
00103 }
00105 template <class PixelType>
00106 bool zxhMorphologyT<PixelType>::EuclideanTransformBk3D(zxhImageDataT<PixelType>*pImage,zxhImageDataT<float>*pOut,PixelType background)
00107 {
00108     if(pImage->GetDimension()!=3) return false;
00109     zxhImageDataT<bool>* p1 = new zxhImageDataT<bool>();//background == false
00110     zxhImageDataT<bool>* p2 = 0;
00111     float tsp[] = {1,1,1,1} ; pImage->GetImageSpacing( tsp[0], tsp[1], tsp[2], tsp[3] ) ;
00112     zxhImageInfo imageinfo ;  pImage->GetImageInfo( &imageinfo ) ;
00113     p1->NewImage( pImage->GetDimension(), pImage->GetImageSize(), tsp, &imageinfo ) ;
00114     if(pOut==0)
00115     {   pOut->NewImage(p1->GetDimension(),p1->GetImageSize(), tsp, &imageinfo ); }
00116     if(pOut->GetDimension()!=3)
00117         return false;
00118     const int *size = p1->GetImageSize();
00119     float maxdistance   = size[0]*tsp[0] + size[1]*tsp[1] + size[2]*tsp[2] ;
00120     bool bFinish    = true;
00121     int *pBoundary1 = new int [ (size[0]*size[1] + size[1]*size[2] + size[0]*size[2]) * (3*6*2) ] ;// currently max subregions 8
00122     int *pBoundary2 = new int [ (size[0]*size[1] + size[1]*size[2] + size[0]*size[2]) * (3*6*2) ] ;
00123     int nBoundary1  = 0 ;
00124     int nBoundary2  = 0;
00125     for( int x=0 ; x<size[0] ; ++x )
00126     for( int y=0 ; y<size[1] ; ++y )
00127     for( int z=0 ; z<size[2] ; ++z )
00128     {
00129         if( pImage->GetPixelGreyscale( x,y,z ) == background )
00130             pOut->SetPixelByGreyscale( x, y, z, 0, maxdistance ) ;
00131         else
00132             pOut->SetPixelByGreyscale( x, y, z, 0, 0 ) ;
00133     }
00134     for( int x=1 ; x<size[0]-1 ; ++x )
00135     for( int y=1 ; y<size[1]-1 ; ++y )
00136     for( int z=1 ; z<size[2]-1 ; ++z )
00137     {
00138         if( pImage->GetPixelGreyscale( x,y,z ) == background )
00139         {
00140             p1->SetPixelByGreyscale( x, y, z, 0, false ) ;
00141         }
00142         else //not background, find the boundary
00143         {
00144             p1->SetPixelByGreyscale( x, y, z, 0, true ) ;
00145             // 6 neighbour points
00146             if( pImage->GetPixelGreyscale( x-1, y, z ) == background || pImage->GetPixelGreyscale( x+1, y, z ) == background ||
00147                 pImage->GetPixelGreyscale( x, y-1, z ) == background || pImage->GetPixelGreyscale( x, y+1, z ) == background ||
00148                 pImage->GetPixelGreyscale( x, y, z-1 ) == background || pImage->GetPixelGreyscale( x, y, z+1 ) == background )
00149             {
00150                 pBoundary1[ nBoundary1 * 3 ]    = x ;
00151                 pBoundary1[ nBoundary1 * 3 + 1] = y ;
00152                 pBoundary1[ nBoundary1 * 3 + 2] = z ;
00153                 ++nBoundary1 ;
00154             }
00155         }
00156     }
00157     for( int x = 0 ; x < size[0] ; ++x )
00158     for( int y = 0 ; y < size[1] ; ++y )
00159     {
00160         p1->SetPixelByGreyscale( x, y, 0, 0, pImage->GetPixelGreyscale(x,y,0) != background );
00161         p1->SetPixelByGreyscale( x, y, size[2]-1, 0, pImage->GetPixelGreyscale(x,y,size[2]-1) != background );
00162     }
00163     for( int y = 0 ; y < size[1] ; ++y )
00164     for( int z = 0 ; z < size[2] ; ++z )
00165     {
00166         p1->SetPixelByGreyscale( 0, y, z, 0, pImage->GetPixelGreyscale(0,y,z) != background );
00167         p1->SetPixelByGreyscale( size[0]-1, y, z, 0, pImage->GetPixelGreyscale(size[0]-1,y,z) != background );
00168     }
00169     for( int x = 0 ; x < size[0] ; ++x )
00170     for( int z = 0 ; z < size[2] ; ++z )
00171     {
00172         p1->SetPixelByGreyscale( x, 0, z, 0, pImage->GetPixelGreyscale( x, 0, z ) != background );
00173         p1->SetPixelByGreyscale( x, size[1]-1, z, 0, pImage->GetPixelGreyscale( x, size[1]-1, z ) != background );
00174     }
00175     p2  = p1->Clone( p1 );
00176
00177     float currdis = 0;
00178     int currx,curry,currz;
00179     int index[6][3] = {-1,0,0, 1,0,0, 0,-1,0, 0,1,0, 0,0,-1, 0,0,1};
00180     float spacingdis[6] = { tsp[0],tsp[0], tsp[1],tsp[1], tsp[2],tsp[2] };
00181     while( nBoundary1 > 0 )
00182     {
00183         int x, y, z;
00184         nBoundary2  = 0;
00185         for( int iB = 0 ; iB < nBoundary1 ; ++iB )
00186         {
00187             currx = pBoundary1[ 3*iB ];
00188             curry = pBoundary1[ 3*iB + 1 ];
00189             currz = pBoundary1[ 3*iB + 2 ];
00190             currdis = pOut->GetPixelGreyscale( currx, curry, currz );
00191             for( int cuindex = 0 ; cuindex < 6 ; ++cuindex )
00192             {
00193                 x   = currx + index[cuindex][0] ;
00194                 y   = curry + index[cuindex][1] ;
00195                 z   = currz + index[cuindex][2] ;
00196                 x = x < 0 ? 0 : x ;
00197                 x = x > size[0]-1 ? size[0]-1 : x ;
00198                 y = y < 0 ? 0 : y ;
00199                 y = y > size[1]-1 ? size[1]-1 : y ;
00200                 z = z < 0 ? 0 : z ;
00201                 z = z > size[2]-1 ? size[2]-1 : z ;
00202                 if( p1->GetPixelGreyscale( x, y, z )== false )
00203                 {
00204                     if( spacingdis[ cuindex ] + currdis < pOut->GetPixelGreyscale( x, y, z ) )
00205                         pOut->SetPixelByGreyscale( x, y, z, 0, spacingdis[ cuindex ] + currdis);
00206                     if( p2->GetPixelGreyscale( x, y, z )== false )
00207                     {
00208                         pBoundary2[ nBoundary2 * 3 ]    = x ;
00209                         pBoundary2[ nBoundary2 * 3 + 1] = y ;
00210                         pBoundary2[ nBoundary2 * 3 + 2] = z ;
00211                         ++nBoundary2 ;
00212                         p2->SetPixelByGreyscale( currx+index[cuindex][0], curry+index[cuindex][1], currz+index[cuindex][2], 0, true );
00213                     }
00214                 }
00215             }
00216         }
00217         nBoundary1 = nBoundary2 ;
00218         int * ptemp = pBoundary1 ;
00219         pBoundary1  = pBoundary2 ;
00220         pBoundary2  = ptemp ;
00221         for( int iB = 0 ; iB < nBoundary1 ; ++iB )
00222             p1->SetPixelByGreyscale( pBoundary1[ 3*iB ], pBoundary1[ 3*iB +1 ], pBoundary1[ 3*iB +2 ] , 0, true ) ;
00223     }
00224     delete p1, p2, pBoundary1, pBoundary2 ;
00225     return true ;
00226 }
00227
00228
00230 template <class PixelType>
00231 bool zxhMorphologyT<PixelType>::Dilation3D(zxhImageDataT<PixelType> *pImage, float step, PixelType background )
00232 {
00233     if( pImage==0 || pImage->GetDimension()<2 || (step<1&&step>-1) ) return false;
00234
00235     if( step<=-1 ) return  zxhMorphologyT<PixelType>::Erosion3D( pImage, -step, background ) ;
00236
00237     zxhImageDataT<PixelType> *pBk = 0 ;
00238     pBk = pImage->Clone( pBk ) ;
00239     PixelType value ;
00240     int *kernal_index[3], iLength= 1, nLength=0 ;
00241     if( pImage->GetDimension()==2 )
00242         iLength = int( (step*2+1)*(step*2+1) ) ;
00243     else iLength= int( (step*2+1)*(step*2+1)*(step*2+1) ) ;
00244     kernal_index[0] = new int [iLength];
00245     kernal_index[1] = new int [iLength];
00246     kernal_index[2] = new int [iLength];
00247     for( int i=0; i<iLength; ++i )
00248         kernal_index[0][i]=kernal_index[1][i]=kernal_index[2][i]=0 ;
00249     int centre = int(step) ;
00250
00251     if( pImage->GetDimension()==2 ) //2D
00252     {
00253         nLength = 0 ;
00254         for( int ix=-centre; ix<=centre; ++ix )
00255         for( int iy=-centre; iy<=centre; ++iy )
00256         {
00257             if( (ix)*(ix)+(iy)*(iy) <=step*step )
00258             {
00259                 kernal_index[0][nLength] = ix ;
00260                 kernal_index[1][nLength] = iy ;
00261                 nLength++;
00262             }
00263         }
00264         for( int iy=1; iy<=pImage->GetImageSize()[1]-1-1 ; ++iy )
00265         for( int ix=1; ix<=pImage->GetImageSize()[0]-1-1 ; ++ix )
00266         {
00267             value = pBk->GetPixelGreyscale(ix,iy,0,0) ;
00268             if( value != background )
00269             {
00270                 if( pBk->GetPixelGreyscale( ix+1,iy ) == background
00271                     || pBk->GetPixelGreyscale( ix-1,iy ) == background
00272                     || pBk->GetPixelGreyscale( ix,iy+1 ) == background
00273                     || pBk->GetPixelGreyscale( ix,iy-1 ) == background )
00274                 {
00275                     for( int ip=0; ip<nLength; ++ip )
00276                         if( pImage->InsideImage( ix+kernal_index[0][ip], iy+kernal_index[1][ip], 0,0 ) )
00277                             pImage->SetPixelByGreyscale( ix+kernal_index[0][ip], iy+kernal_index[1][ip], 0, 0, value );
00278                 }
00279             }
00280         }
00281     }
00282     else
00283     {
00284         nLength = 0 ;
00285         for( int iz=-centre; iz<=centre; ++iz )
00286         for( int ix=-centre; ix<=centre; ++ix )
00287         for( int iy=-centre; iy<=centre; ++iy )
00288         {
00289             if( (ix)*(ix)+(iy)*(iy)+iz*iz <=step*step )
00290             {
00291                 kernal_index[0][nLength] = ix ;
00292                 kernal_index[1][nLength] = iy ;
00293                 kernal_index[2][nLength] = iz ;
00294                 nLength++;
00295             }
00296         }
00297         for( int it=0; it<=pImage->GetImageSize()[3]-1; ++it )
00298         for( int iz=1; iz<=pImage->GetImageSize()[2]-1-1 ; ++iz )
00299         for( int iy=1; iy<=pImage->GetImageSize()[1]-1-1 ; ++iy )
00300         for( int ix=1; ix<=pImage->GetImageSize()[0]-1-1 ; ++ix )
00301         {
00302             value = pBk->GetPixelGreyscale(ix,iy,iz,it) ;
00303             if( value != background )
00304             {
00305                 if( pBk->GetPixelGreyscale( ix+1,iy,iz,it ) == background
00306                     || pBk->GetPixelGreyscale( ix-1,iy,iz,it ) == background
00307                     || pBk->GetPixelGreyscale( ix,iy+1,iz,it ) == background
00308                     || pBk->GetPixelGreyscale( ix,iy-1,iz,it ) == background
00309                     || pBk->GetPixelGreyscale( ix,iy,iz+1,it ) == background
00310                     || pBk->GetPixelGreyscale( ix,iy,iz-1,it ) == background )
00311                 {
00312                     for( int ip=0; ip<nLength; ++ip )
00313                         if( pImage->InsideImage( ix+kernal_index[0][ip], iy+kernal_index[1][ip], iz+kernal_index[2][ip], it ) )
00314                             pImage->SetPixelByGreyscale( ix+kernal_index[0][ip], iy+kernal_index[1][ip], iz+kernal_index[2][ip], it, value );
00315                 }
00316             }
00317         }
00318     }
00319     delete pBk ;
00320     delete []kernal_index[0];
00321     delete []kernal_index[1];
00322     delete []kernal_index[2];
00323     return true;
00324 }
00326 template <class PixelType>
00327 bool zxhMorphologyT<PixelType>::DilateStep3D(zxhImageDataT<PixelType> *pImage,
00328                     int step, PixelType background, PixelType foreground )
00329 {
00330     if( pImage==0 || pImage->GetDimension()<2 || (step<1&&step>-1) ) return false;
00331
00332     if( step<=-1 ) zxhMorphologyT<PixelType>::ErodeStep3D( pImage, -step, background, foreground ) ;
00333
00334     zxhImageDataT<PixelType> * pBK=0;
00335     for(int istep=0; istep<step; ++istep)
00336     {
00337         pBK = pImage->Clone( pBK );
00338         if( pImage->GetDimension()==2 ) //2D
00339         {
00340             int iz=0, it=0 ;
00341             for( int iy=1; iy<=pImage->GetImageSize()[1]-1-1 ; ++iy )
00342             for( int ix=1; ix<=pImage->GetImageSize()[0]-1-1 ; ++ix )
00343             {
00344                 if( pBK->GetPixelGreyscale(ix,iy,iz,it) != background )
00345                 {
00346                     pImage->SetPixelByGreyscale( ix+1,iy,iz,it, foreground ) ;
00347                     pImage->SetPixelByGreyscale( ix-1,iy,iz,it, foreground ) ;
00348                     pImage->SetPixelByGreyscale( ix,iy+1,iz,it, foreground ) ;
00349                     pImage->SetPixelByGreyscale( ix,iy-1,iz,it, foreground ) ;
00350                 }
00351             }
00352         }
00353         else
00354         {
00355             for( int it=0; it<=pImage->GetImageSize()[3]-1; ++it )
00356             for( int iz=1; iz<=pImage->GetImageSize()[2]-1-1 ; ++iz )
00357             for( int iy=1; iy<=pImage->GetImageSize()[1]-1-1 ; ++iy )
00358             for( int ix=1; ix<=pImage->GetImageSize()[0]-1-1 ; ++ix )
00359             {
00360                 if( pBK->GetPixelGreyscale(ix,iy,iz,it) != background )
00361                 {
00362                     pImage->SetPixelByGreyscale( ix+1,iy,iz,it, foreground ) ;
00363                     pImage->SetPixelByGreyscale( ix-1,iy,iz,it, foreground ) ;
00364                     pImage->SetPixelByGreyscale( ix,iy+1,iz,it, foreground ) ;
00365                     pImage->SetPixelByGreyscale( ix,iy-1,iz,it, foreground ) ;
00366                     pImage->SetPixelByGreyscale( ix,iy,iz+1,it, foreground ) ;
00367                     pImage->SetPixelByGreyscale( ix,iy,iz-1,it, foreground ) ;
00368
00369                 }
00370             }
00371         }
00372     }
00373     delete pBK;
00374     return true;
00375 }
00377 template <class PixelType>
00378 bool zxhMorphologyT<PixelType>::ErodeStep3D(zxhImageDataT<PixelType> *pImage,
00379                     int step, PixelType background, PixelType foreground )
00380 {
00381     if( pImage==0 || pImage->GetDimension()<2 || (step<1&&step>-1) ) return false;
00382
00383     if( step<=-1 ) zxhMorphologyT<PixelType>::DilateStep3D( pImage, -step, background, foreground ) ;
00384
00385     zxhImageDataT<PixelType> * pBK=0;
00386     for(int istep=0; istep<step; ++istep)
00387     {
00388         pBK = pImage->Clone( pBK );
00389         if( pImage->GetDimension()==2 ) //2D
00390         {
00391             for( int iy=0; iy<=pImage->GetImageSize()[1]-1; iy+=(pImage->GetImageSize()[1]-1) )
00392             for( int ix=0; ix<=pImage->GetImageSize()[0]-1; ix+=1 )
00393                 pImage->SetPixelByGreyscale( ix, iy, 0, 0, background );
00394
00395             for( int ix=0; ix<=pImage->GetImageSize()[0]-1; ix+= (pImage->GetImageSize()[0]-1))
00396             for( int iy=0; iy<=pImage->GetImageSize()[1]-1; iy+=1 )
00397                 pImage->SetPixelByGreyscale( ix, iy, 0, 0, background );
00398
00399             int iz=0, it=0 ;
00400             for( int iy=1; iy<=pImage->GetImageSize()[1]-1-1 ; ++iy )
00401             for( int ix=1; ix<=pImage->GetImageSize()[0]-1-1 ; ++ix )
00402             {
00403                 if( pBK->GetPixelGreyscale(ix,iy,iz,it) == background )
00404                 {
00405                     pImage->SetPixelByGreyscale( ix+1,iy,iz,it, background ) ;
00406                     pImage->SetPixelByGreyscale( ix-1,iy,iz,it, background ) ;
00407                     pImage->SetPixelByGreyscale( ix,iy+1,iz,it, background ) ;
00408                     pImage->SetPixelByGreyscale( ix,iy-1,iz,it, background ) ;
00409                 }
00410             }
00411         }
00412         else
00413         {
00414             for( int it=0; it<=pImage->GetImageSize()[3]-1; ++it ) // erode roi boundary
00415             {
00416                 for( int iz=0; iz<=pImage->GetImageSize()[2]-1; iz+= (pImage->GetImageSize()[2]-1))
00417                 for( int iy=0; iy<=pImage->GetImageSize()[1]-1; iy+=1 )
00418                 for( int ix=0; ix<=pImage->GetImageSize()[0]-1; ix+=1 )
00419                     pImage->SetPixelByGreyscale( ix, iy, iz, it, background );
00420
00421                 for( int iy=0; iy<=pImage->GetImageSize()[1]-1; iy+= (pImage->GetImageSize()[1]-1) )
00422                 for( int ix=0; ix<=pImage->GetImageSize()[0]-1; ix+=1 )
00423                 for( int iz=0; iz<=pImage->GetImageSize()[2]-1; iz+=1 )
00424                     pImage->SetPixelByGreyscale( ix, iy, iz, it, background );
00425
00426                 for( int ix=0; ix<=pImage->GetImageSize()[0]-1; ix+= (pImage->GetImageSize()[0]-1))
00427                 for( int iy=0; iy<=pImage->GetImageSize()[1]-1; iy+=1 )
00428                 for( int iz=0; iz<=pImage->GetImageSize()[2]-1; iz+=1 )
00429                     pImage->SetPixelByGreyscale( ix, iy, iz, it, background );
00430             }
00431
00432             for( int it=0; it<=pImage->GetImageSize()[3]-1; ++it )
00433             for( int iz=1; iz<=pImage->GetImageSize()[2]-1-1 ; ++iz )
00434             for( int iy=1; iy<=pImage->GetImageSize()[1]-1-1 ; ++iy )
00435             for( int ix=1; ix<=pImage->GetImageSize()[0]-1-1 ; ++ix )
00436             {
00437                 if( pBK->GetPixelGreyscale(ix,iy,iz,it) == background )
00438                 {
00439                     pImage->SetPixelByGreyscale( ix+1,iy,iz,it, background ) ;
00440                     pImage->SetPixelByGreyscale( ix-1,iy,iz,it, background ) ;
00441                     pImage->SetPixelByGreyscale( ix,iy+1,iz,it, background ) ;
00442                     pImage->SetPixelByGreyscale( ix,iy-1,iz,it, background ) ;
00443                     pImage->SetPixelByGreyscale( ix,iy,iz+1,it, background ) ;
00444                     pImage->SetPixelByGreyscale( ix,iy,iz-1,it, background ) ;
00445
00446                 }
00447             }
00448         }
00449     }
00450     delete pBK;
00451     return true;
00452 }
00453
00454
00455 template <class PixelType>
00456 bool zxhMorphologyT<PixelType>::Erosion3D( zxhImageDataT<PixelType> *pImage, float step,PixelType background )
00457 {
00458     if( pImage==0 || pImage->GetDimension()<2 || (step<1&&step>-1) ) return false;
00459
00460     if( step<=-1 ) zxhMorphologyT<PixelType>::Dilation3D( pImage, -step, background ) ;
00461
00462     int *kernal_index[3], iLength= 1, nLength=0 ;
00463     if( pImage->GetDimension()==2 )
00464         iLength = int( (step*2+1)*(step*2+1) ) ;
00465     else iLength= int( (step*2+1)*(step*2+1)*(step*2+1) ) ;
00466     kernal_index[0] = new int [iLength];
00467     kernal_index[1] = new int [iLength];
00468     kernal_index[2] = new int [iLength];
00469     for( int i=0; i<iLength; ++i )
00470         kernal_index[0][i]=kernal_index[1][i]=kernal_index[2][i]=0 ;
00471     int centre = int(step-1) ;
00472     zxhImageDataT<PixelType> *pBk = 0 ;
00473     pBk = pImage->Clone( pBk ) ;
00474     PixelType value ;
00475
00476     if( pImage->GetDimension()==2 )
00477     {
00478         nLength = 0 ;
00479         int it=0, iz=0 ;
00480         for( int iy=-centre; iy<=centre; ++iy )
00481         for( int ix=-centre; ix<=centre; ++ix )
00482             if( ix*ix+iy*iy <= (step-1)*(step-1) )
00483             {
00484                 kernal_index[0][nLength] = ix ;
00485                 kernal_index[1][nLength] = iy ;
00486                 ++nLength ;
00487             }
00488         {
00489             for( int iy=0; iy<=pImage->GetImageSize()[1]-1; iy+=(pImage->GetImageSize()[1]-1) )
00490             for( int ix=0; ix<=pImage->GetImageSize()[0]-1; ix+=1 )
00491                 pImage->SetPixelByGreyscale( ix, iy, iz, it, background );
00492
00493             for( int ix=0; ix<=pImage->GetImageSize()[0]-1; ix+= (pImage->GetImageSize()[0]-1))
00494             for( int iy=0; iy<=pImage->GetImageSize()[1]-1; iy+=1 )
00495                 pImage->SetPixelByGreyscale( ix, iy, iz, it, background );
00496         }
00497         for( int iy=1; iy<=pImage->GetImageSize()[1]-1-1 ; ++iy )
00498         for( int ix=1; ix<=pImage->GetImageSize()[0]-1-1 ; ++ix )
00499         {
00500             value = pBk->GetPixelGreyscale(ix,iy,iz,it) ;
00501             if( value != background )
00502             {
00503                 if( pBk->GetPixelGreyscale( ix+1,iy,iz,it ) == background
00504                     || pBk->GetPixelGreyscale( ix-1,iy,iz,it ) == background
00505                     || pBk->GetPixelGreyscale( ix,iy+1,iz,it ) == background
00506                     || pBk->GetPixelGreyscale( ix,iy-1,iz,it ) == background )
00507                 {
00508                     for( int ip=0; ip<nLength; ++ip )
00509                         if( pImage->InsideImage( ix+kernal_index[0][ip], iy+kernal_index[1][ip], iz, it ) )
00510                             pImage->SetPixelByGreyscale( ix+kernal_index[0][ip], iy+kernal_index[1][ip], iz, it, background );
00511                 }
00512             }
00513         }
00514     }
00515     else
00516     {
00517         nLength = 0 ;
00518         for( int iz=-centre; iz<=centre; ++iz )
00519         for( int iy=-centre; iy<=centre; ++iy )
00520         for( int ix=-centre; ix<=centre; ++ix )
00521             if( ix*ix+iy*iy+iz*iz <= (step-1)*(step-1) )
00522             {
00523                 kernal_index[0][nLength] = ix ;
00524                 kernal_index[1][nLength] = iy ;
00525                 kernal_index[2][nLength] = iz ;
00526                 ++nLength ;
00527             }
00528         for( int it=0; it<=pImage->GetImageSize()[3]-1; ++it ) // erode roi boundary
00529         {
00530             for( int iz=0; iz<=pImage->GetImageSize()[2]-1; iz+= (pImage->GetImageSize()[2]-1))
00531             for( int iy=0; iy<=pImage->GetImageSize()[1]-1; iy+=1 )
00532             for( int ix=0; ix<=pImage->GetImageSize()[0]-1; ix+=1 )
00533                 pImage->SetPixelByGreyscale( ix, iy, iz, it, background );
00534
00535             for( int iy=0; iy<=pImage->GetImageSize()[1]-1; iy+= (pImage->GetImageSize()[1]-1) )
00536             for( int ix=0; ix<=pImage->GetImageSize()[0]-1; ix+=1 )
00537             for( int iz=0; iz<=pImage->GetImageSize()[2]-1; iz+=1 )
00538                 pImage->SetPixelByGreyscale( ix, iy, iz, it, background );
00539
00540             for( int ix=0; ix<=pImage->GetImageSize()[0]-1; ix+= (pImage->GetImageSize()[0]-1))
00541             for( int iy=0; iy<=pImage->GetImageSize()[1]-1; iy+=1 )
00542             for( int iz=0; iz<=pImage->GetImageSize()[2]-1; iz+=1 )
00543                 pImage->SetPixelByGreyscale( ix, iy, iz, it, background );
00544         }
00545
00546         for( int it=0; it<=pImage->GetImageSize()[3]-1; ++it )
00547         for( int iz=1; iz<=pImage->GetImageSize()[2]-1-1 ; ++iz )
00548         for( int iy=1; iy<=pImage->GetImageSize()[1]-1-1 ; ++iy )
00549         for( int ix=1; ix<=pImage->GetImageSize()[0]-1-1 ; ++ix )
00550         {
00551             value = pBk->GetPixelGreyscale(ix,iy,iz,it) ;
00552             if( value != background )
00553             {
00554                 if( pBk->GetPixelGreyscale( ix+1,iy,iz,it ) == background
00555                     || pBk->GetPixelGreyscale( ix-1,iy,iz,it ) == background
00556                     || pBk->GetPixelGreyscale( ix,iy+1,iz,it ) == background
00557                     || pBk->GetPixelGreyscale( ix,iy-1,iz,it ) == background
00558                     || pBk->GetPixelGreyscale( ix,iy,iz+1,it ) == background
00559                     || pBk->GetPixelGreyscale( ix,iy,iz-1,it ) == background )
00560                 {
00561                     for( int ip=0; ip<nLength; ++ip )
00562                         if( pImage->InsideImage( ix+kernal_index[0][ip], iy+kernal_index[1][ip], iz+kernal_index[2][ip], it ) )
00563                             pImage->SetPixelByGreyscale( ix+kernal_index[0][ip], iy+kernal_index[1][ip], iz+kernal_index[2][ip], it, background );
00564                 }
00565             }
00566         }
00567     }
00568     delete []kernal_index[0] ;
00569     delete []kernal_index[1] ;
00570     delete []kernal_index[2] ;
00571     delete pBk;
00572     return true;
00573 }
00575 template <class PixelType>
00576 bool zxhMorphologyT<PixelType>::Close3D( zxhImageDataT<PixelType> *pImage, float steps,PixelType background )
00577 {
00578     if( pImage==0 || pImage->GetDimension()<2) return false;
00579     zxhImageDataT<PixelType> *pBk = 0 ;
00580     pBk = pImage->Clone( pBk ) ;
00581     {
00582         Dilation3D( pImage, steps, background ) ;
00583         Erosion3D( pImage, steps, background ) ;
00584     }
00585     return true;
00586 }
00588 template <class PixelType>
00589 bool zxhMorphologyT<PixelType>::Open3D( zxhImageDataT<PixelType> *pImage, float steps,PixelType background )
00590 {
00591     if( pImage==0 || pImage->GetDimension()<2) return false;
00592     zxhImageDataT<PixelType> *pBk = 0 ;
00593     pBk = pImage->Clone( pBk ) ;
00594     {
00595         Erosion3D( pImage, steps, background ) ;
00596         Dilation3D( pImage, steps, background ) ;
00597     }
00598     return true;
00599 }
00600
00602 template <class PixelType>
00603 bool zxhMorphologyT<PixelType>::MaskImageOpAnd3D( zxhImageDataT<PixelType>* pDes, zxhImageDataT<PixelType>*pSrc )
00604 {
00605     if( pSrc==0 || pSrc->GetImageData()==0 || pDes==0 || pDes->GetImageData()==0 ) return false ;
00606     bool same = pDes->ImageNvoxelDimensionSame( pSrc ) ;
00607     if( same )
00608     {
00609         long int lsize = pDes->GetNumberOfPixels() ;
00610         for( long int ip=0; ip<lsize; ip++ )
00611             if( pSrc->GetImageData( ip ) == ZXH_Background )
00612                 pDes->SetImageData( ip, ZXH_Background ) ;
00613     }
00614     else
00615     {
00616         zxhImageModelingLinearT<PixelType,float> model ;
00617         model.SetImage( pSrc );
00618         for( int it=0; it<=pDes->GetImageSize()[3]-1 ; ++it )
00619         for( int iz=0; iz<=pDes->GetImageSize()[2]-1 ; ++iz )
00620         for( int iy=0; iy<=pDes->GetImageSize()[1]-1 ; ++iy )
00621         for( int ix=0; ix<=pDes->GetImageSize()[0]-1 ; ++ix )
00622         {
00623             float World[] = {ix,iy,iz,it} ;
00624             pDes->ImageToWorld( World ) ;
00625             if( model.GetPixelFloatValueWithCheckByWorld( World[0], World[1], World[2], World[3] ) == ZXH_Background )
00626                 pDes->SetPixelByGreyscale( ix,iy,iz,it, ZXH_Background ) ;
00627         }
00628     }
00629     return true;
00630 }
00631 #endif //zxhMorphology_h
00632
00633
00634