zxhproj v 2.2
zxhproj

zxhMultiCannyEdge.h

00001 
00002 /*=========================================================================
00003 
00004   Program:   ZXH Registration Software
00005   Author:    Xiahai Zhuang
00006   Module:    $RCSfle: zxhMultiCannyEdge.h    $
00007   Language:  C++
00008   Date:      $Date: From  2007-11 $
00009   Version:   $Revision: 2.0 $
00010 
00011 =========================================================================*/
00012 
00013 #ifndef ZXHMULTICANNYEDGE_H
00014 #define ZXHMULTICANNYEDGE_H
00015 
00016 #include "zxh.h"
00017 #include "zxhEdgeBase.h"
00018 #include "zxhEdgeCanny.h"
00019 #include "zxhImageData.h"
00020 #include "zxhImageModelingBase.h"
00021 #include "zxhImageModelingLinear.h"
00022 
00027 class zxhMultiCannyEdge :
00028     public zxhEdgeCanny
00029 {
00030 public:
00031     static const int EdgePoint = ZXH_Foreground  ;
00032     static const int FoundEdge = -ZXH_Foreground  ;
00033     static const int InactivateGrayRange = -ZXH_Foreground ;
00035     zxhMultiCannyEdge(void);
00037     virtual ~zxhMultiCannyEdge(void);
00038 
00040     virtual bool    Evaluate()
00041     {
00042         int i = SearchMultiCannyEdges( m_fSearchRange, m_fSearchStepLength, //unit mm
00043             m_pMaskImage, m_afGrayRange) ;
00044         SmoothBoundary( m_iMainTimes, m_fNeighbour, m_fAccuracy );
00045         return i>0;
00046     }
00048     virtual void    SetSurfacBaseImage( zxhImageData* p )   {   m_pSurfaceBase=p;   } ;
00050     virtual void    SetAtlasBaseImage( zxhImageData* p )    {   m_pAtlasBase=p;     } ;
00052     virtual void    SetBluredImage( zxhImageData* p )       {   m_pImageBlur=p;     } ;
00054     virtual void    SetNumberMultiPartialEdgeImages(int i)  {   m_nPartialEdges=i;  } ;
00056     virtual void    SetMultiPartialEdgeImages(zxhImageData**p)  {   m_pPartialEdges=p;  } ;
00058     virtual zxhImageDataT<float>*   GetDisplaceVectorFromSurface()      {   return m_aDVsurf;   } ;
00060     virtual zxhImageDataT<float>*   GetDisplaceVectorFromCorrect()      {   return m_aDVcor;    } ;
00062     virtual void    SetSearchingRange( float r )            {   m_fSearchRange=r;   };
00064     virtual void    SetSearchStepLength( float l )          {   m_fSearchStepLength = l ;   };
00066     virtual void    SetMaskImage( zxhImageData* p )         {   m_pMaskImage = p;   } ;
00068     virtual void    SetEdgePointGrayRange( float f, float t )   { m_afGrayRange[0] = f; m_afGrayRange[1] = t; }
00070     virtual void    SetMaintainTimes( int it )              {   m_iMainTimes = it;  };
00072     virtual void    SetMaintainNeighbour( float f )         {   m_fNeighbour = f;   };
00074     virtual void    SetMaintainAccuracy( float f )          {   m_fAccuracy = f;    };
00075 
00077     virtual long int    SearchMultiCannyEdges( float fSearchRange, float fSearchStepLength, //unit mm
00078                                     zxhImageData* pMask, float aGrayRange[2]);
00080     virtual long int    SmoothBoundary( int nMainTimes, float fNeighbour, float fAccuracy );
00081 
00083     virtual long int    RemovingJaggedShape( int nMainTimes, float fNeighbour, float fAccuracy ) ;
00085     virtual long int    FillingJaggedShape( int nMainTimes, float fNeighbour, float fAccuracy )
00086     { return -1 ; std::cerr<<"error: havenot implemented\n";}
00087 
00088 protected:
00089     // surface and atlas and partialedge images are same size and spacing
00091     zxhImageData*   m_pSurfaceBase ;
00093     zxhImageData*   m_pAtlasBase ;
00095     zxhImageData*   m_pImageBlur;
00097     int m_nPartialEdges ;
00099     zxhImageData**  m_pPartialEdges;
00101     zxhImageDataT<float>    m_aDVsurf[3];
00103     zxhImageDataT<float>    m_aDVcor[3];
00104 
00106     float       m_fSearchRange ;
00108     float       m_fSearchStepLength ;
00110     zxhImageData*       m_pMaskImage ;
00112     float       m_afGrayRange[2] ;
00114     int         m_iMainTimes ;
00116     float       m_fNeighbour ;
00118     float       m_fAccuracy ;
00119 };
00120 
00121 
00122 zxhMultiCannyEdge::zxhMultiCannyEdge(void):
00123 m_pSurfaceBase(0),
00124 m_pAtlasBase(0),
00125 m_pImageBlur(0),
00126 m_nPartialEdges(0),
00127 m_pPartialEdges(0),
00128 m_fSearchRange (5),
00129 m_fSearchStepLength (1),
00130 m_pMaskImage (0),
00131 m_iMainTimes (1),
00132 m_fNeighbour (5),
00133 m_fAccuracy (2)
00134 {
00135     m_afGrayRange[0] = InactivateGrayRange ;
00136     m_afGrayRange[1] = InactivateGrayRange ;
00137 }
00138 
00139 zxhMultiCannyEdge::~zxhMultiCannyEdge(void)
00140 {
00141 }
00142 
00143 long int zxhMultiCannyEdge::SearchMultiCannyEdges( float fSearchRange, float fSearchStepLength,
00144                                               zxhImageData* pMask, float aGrayRange[2])
00145 {
00146     if( m_pSurfaceBase == 0 )
00147     {
00148         std::cerr<<"error: input the surface base image \n" ;
00149         return -1 ;
00150     }
00151 
00152     int size[] = {1,1,1,1} ; this->m_pSurfaceBase->GetImageSize( size[0],size[1],size[2],size[3] ) ;
00153     float sp[] = {1,1,1,1} ; this->m_pSurfaceBase->GetImageSpacing( sp[0],sp[1],sp[2],sp[3] ) ;
00154     int dimension = this->m_pSurfaceBase->GetDimension() ;
00155 
00156     zxhImageInfo imageinfo ;
00157 
00158     zxhImageModelingLinear modelInput, modelBlur, modelAtlas;
00159     zxhImageModelingLinearT<float,float> modelCannyImage ;
00160     modelInput.SetImage( m_pImage ) ;
00161 
00162     this->m_pSurfaceBase->GetImageInfo( &imageinfo ) ;
00163 
00164     if( m_pImage->GetImageInfo()->SameDimSizeSpacingAs( &imageinfo ) )
00165     {
00166         modelInput.ResizeImage( size, sp ) ;
00167         if( glbVerboseOutput>0 ) std::cout<<"success: resize the input base image \n" ;
00168     }
00169 
00170     if( m_pImageBlur != 0 )
00171     {
00172         modelBlur.SetImage(m_pImageBlur);
00173         if( m_pImageBlur->GetImageInfo()->SameDimSizeSpacingAs( &imageinfo ) )
00174         {
00175             modelBlur.ResizeImage( size, sp ) ;
00176             if( glbVerboseOutput>0 ) std::cout<<"success: resize the blured base image \n" ;
00177         }
00178     }
00179     else modelBlur.SetImage( modelInput.GetImage() ) ;
00180 
00181     this->SetEdgeImageSpacing( sp ) ;
00182     if( this->EvaluateCannyImage() == false )
00183         return false ;
00184     modelCannyImage.SetImage( &this->m_CannyImage ) ;
00185     modelAtlas.SetImage( this->m_pAtlasBase ) ;
00186 
00187     this->m_pImage->GetImageInfo( &imageinfo ) ;
00188 
00189     this->m_EdgeImage.NewImage( dimension, size, sp, &imageinfo ) ;
00190     for( int id=0; id<3; ++id )
00191     {
00192         m_aDVsurf[id].NewImage(dimension, size, sp, &imageinfo ) ;
00193         m_aDVcor[id].NewImage(dimension, size, sp, &imageinfo ) ;
00194     }
00195     float NV[4] = {0,0,0,0};
00196     float maxstep = fSearchRange/sp[0], steplength = fSearchStepLength/sp[0] ;
00197     long int iPoints = 0 ;
00198     for( int it=0; it<size[3]; ++it ) // search for edge points
00199     for( int ix=0; ix<size[0]; ++ix )
00200     for( int iy=0; iy<size[1]; ++iy )
00201     for( int iz=0; iz<size[2]; ++iz )
00202     {
00203         if( m_pSurfaceBase->GetPixelGreyscale( ix,iy,iz,it ) != ZXH_Background ) // surface point
00204         {
00205             // search from aCannyEdge, to pCannyImage
00206             float adjustco[]={ float(ix), float(iy), float(iz), float(it) } ;
00207             //m_pSurfaceBase->AdjustToPhysical( adjustco ) ;
00208             //pimgMRIBlur->AdjustPhysical2Grid( adjustco ) ;
00209             modelBlur.GetPixelGradient( NV, adjustco[0], adjustco[1], adjustco[2], adjustco[3] );
00210             zxh::NormaliseVector( NV, dimension ) ;
00211             bool find = false ;
00212             float edge[4]= {0,0,0, float(it)};
00213             for( int iCannyEdge = 0; iCannyEdge<m_nPartialEdges && find==false; ++iCannyEdge )
00214             {
00215                 for( float fpos = 0 ; fpos< maxstep ; fpos+=steplength )
00216                 {
00217                     edge[0] = ix+NV[0]*fpos;
00218                     edge[1] = iy+NV[1]*fpos;
00219                     edge[2] = iz+NV[2]*fpos;
00220                     bool bmask = true;
00221                     if( pMask != 0 )
00222                         bmask = pMask->InsideImage(edge[0], edge[1], edge[2], edge[3])
00223                                 && pMask->GetPixelGreyscale(edge[0], edge[1], edge[2], edge[3]) != ZXH_Background ;
00224                     if( bmask
00225                         && m_pPartialEdges[iCannyEdge]->InsideImage( edge[0], edge[1], edge[2], edge[3] )
00226                         && m_pPartialEdges[iCannyEdge]->GetPixelGreyscale( edge[0], edge[1], edge[2], edge[3] ) != ZXH_Background )
00227                     {// find a canny edge candidate point
00228                         if( m_pAtlasBase!= 0 )
00229                         {
00230                             if( (modelAtlas.GetPixelFloatValueWithCheck( float(ix)+NV[0], float(iy)+NV[1], float(iz)+NV[2], float(it) )
00231                                   - modelAtlas.GetPixelFloatValueWithCheck( float(ix)-NV[0], float(iy)-NV[1], float(iz)-NV[2], float(it) ))
00232                                 *(modelInput.GetPixelFloatValueWithCheck( edge[0]+NV[0], edge[1]+NV[1], edge[2]+NV[2], float(it) )
00233                                   - modelInput.GetPixelFloatValueWithCheck( edge[0]-NV[0], edge[1]-NV[1], edge[2]-NV[2], float(it) ))
00234                                 > 0 ) // same surface point because similar gradient of intensity
00235                             { find = true ; }
00236                         }
00237                         else { find = true; }
00238                         if( find==true && aGrayRange[0]!=InactivateGrayRange && aGrayRange[1]!=InactivateGrayRange )
00239                         {
00240                             if( modelBlur.GetPixelFloatValueWithCheck(edge[0], edge[1], edge[2], edge[3]) < aGrayRange[0]
00241                                 || modelBlur.GetPixelFloatValueWithCheck(edge[0], edge[1], edge[2], edge[3]) > aGrayRange[1] )
00242                                 find = false ;
00243                         }
00244                         if( find == true ) break;
00245                     }
00246                     edge[0] = ix-NV[0]*fpos;
00247                     edge[1] = iy-NV[1]*fpos;
00248                     edge[2] = iz-NV[2]*fpos;
00249                     bmask = true;
00250                     if( pMask != 0 )
00251                         bmask = pMask->InsideImage(edge[0], edge[1], edge[2], edge[3])
00252                                 && pMask->GetPixelGreyscale(edge[0], edge[1], edge[2], edge[3]) != ZXH_Background ;
00253                     if( bmask
00254                         && m_pPartialEdges[iCannyEdge]->InsideImage( edge[0], edge[1], edge[2], edge[3] )
00255                         && m_pPartialEdges[iCannyEdge]->GetPixelGreyscale( edge[0], edge[1], edge[2], edge[3] ) != ZXH_Background )
00256                     {// find a canny edge candidate point
00257                         if( m_pAtlasBase!= 0 )
00258                         {
00259                             if( (modelAtlas.GetPixelFloatValueWithCheck( float(ix)+NV[0], float(iy)+NV[1], float(iz)+NV[2], float(it) )
00260                                   - modelAtlas.GetPixelFloatValueWithCheck( float(ix)-NV[0], float(iy)-NV[1], float(iz)-NV[2], float(it)))
00261                                 *(modelInput.GetPixelFloatValueWithCheck( edge[0]+NV[0], edge[1]+NV[1], edge[2]+NV[2], float(it) )
00262                                   - modelInput.GetPixelFloatValueWithCheck( edge[0]-NV[0], edge[1]-NV[1], edge[2]-NV[2], float(it) ))
00263                                 > 0 ) // same surface point because similar gradient of intensity
00264                             { find = true ; }
00265                         }
00266                         else { find = true;}
00267                         if( find==true && aGrayRange[0]!=InactivateGrayRange && aGrayRange[1]!=InactivateGrayRange )
00268                         {
00269                             if( modelBlur.GetPixelFloatValueWithCheck(edge[0], edge[1], edge[2], edge[3]) < aGrayRange[0]
00270                                 || modelBlur.GetPixelFloatValueWithCheck(edge[0], edge[1], edge[2], edge[3]) > aGrayRange[1] )
00271                                 find = false ;
00272                         }
00273                         if( find == true ) break;
00274                     }
00275                 }
00276             }
00277             if( find == false ) // search from the canny image for optimal edge
00278             {
00279                 for( float fpos = 0 ; fpos< maxstep ; fpos+=steplength )
00280                 {
00281                     if( modelCannyImage.GetPixelFloatValueWithCheck( float(ix)+NV[0]*fpos, float(iy)+NV[1]*fpos, float(iz)+NV[2]*fpos, float(it) )
00282                        *modelCannyImage.GetPixelFloatValueWithCheck( float(ix)+NV[0]*(fpos+steplength), float(iy)+NV[1]*(fpos+steplength), float(iz)+NV[2]*(fpos+steplength), float(it) )
00283                        < 0 )
00284                     {
00285                         edge[0] = ix+NV[0]*(fpos+steplength*0.5f);
00286                         edge[1] = iy+NV[1]*(fpos+steplength*0.5f);
00287                         edge[2] = iz+NV[2]*(fpos+steplength*0.5f);
00288                         bool bmask = true;
00289                         if( pMask != 0 )
00290                             bmask = pMask->InsideImage(edge[0], edge[1], edge[2], edge[3])
00291                                     && pMask->GetPixelGreyscale(edge[0], edge[1], edge[2], edge[3]) != ZXH_Background ;
00292                         if( bmask)
00293                         {
00294                             if( m_pAtlasBase!= 0 )
00295                             {
00296                                 if( (modelAtlas.GetPixelFloatValueWithCheck( float(ix)+NV[0], float(iy)+NV[1], float(iz)+NV[2], float(it) )
00297                                       - modelAtlas.GetPixelFloatValueWithCheck( float(ix)-NV[0], float(iy)-NV[1], float(iz)-NV[2], float(it) ))
00298                                     *(modelInput.GetPixelFloatValueWithCheck( edge[0]+NV[0], edge[1]+NV[1], edge[2]+NV[2], float(it) )
00299                                       - modelInput.GetPixelFloatValueWithCheck( edge[0]-NV[0], edge[1]-NV[1], edge[2]-NV[2], float(it) ))
00300                                     > 0 ) // same surface point because similar gradient of intensity
00301                                 { find = true ; ; }
00302                             }
00303                             else { find = true;  ;}
00304                             if( find==true && aGrayRange[0]!=InactivateGrayRange && aGrayRange[1]!=InactivateGrayRange  )
00305                             {
00306                                 if( modelBlur.GetPixelFloatValueWithCheck(edge[0], edge[1], edge[2], edge[3]) < aGrayRange[0]
00307                                     || modelBlur.GetPixelFloatValueWithCheck(edge[0], edge[1], edge[2], edge[3]) > aGrayRange[1] )
00308                                     find = false ;
00309                             }
00310                             if( find == true ) break;
00311                         }
00312                     }
00313                     if( modelCannyImage.GetPixelFloatValueWithCheck( float(ix)-NV[0]*fpos, float(iy)-NV[1]*fpos, float(iz)-NV[2]*fpos, float(it) )
00314                        *modelCannyImage.GetPixelFloatValueWithCheck( float(ix)-NV[0]*(fpos+steplength), float(iy)-NV[1]*(fpos+steplength), float(iz)-NV[2]*(fpos+steplength), float(it) )
00315                        < 0 )
00316                     {
00317                         edge[0] = ix-NV[0]*(fpos+steplength*0.5f);
00318                         edge[1] = iy-NV[1]*(fpos+steplength*0.5f);
00319                         edge[2] = iz-NV[2]*(fpos+steplength*0.5f);
00320                         bool bmask = true;
00321                         if( pMask != 0 )
00322                             bmask = pMask->InsideImage(edge[0], edge[1], edge[2], edge[3])
00323                                     && pMask->GetPixelGreyscale(edge[0], edge[1], edge[2], edge[3]) != ZXH_Background ;
00324                         if( bmask )
00325                         {
00326                             if( m_pAtlasBase!= 0 )
00327                             {
00328                                 if( (modelAtlas.GetPixelFloatValueWithCheck( float(ix)+NV[0], float(iy)+NV[1], float(iz)+NV[2], float(it) )
00329                                       - modelAtlas.GetPixelFloatValueWithCheck( float(ix)-NV[0], float(iy)-NV[1], float(iz)-NV[2], float(it) ))
00330                                     *(modelInput.GetPixelFloatValueWithCheck( edge[0]+NV[0], edge[1]+NV[1], edge[2]+NV[2], float(it) )
00331                                       - modelInput.GetPixelFloatValueWithCheck( edge[0]-NV[0], edge[1]-NV[1], edge[2]-NV[2], float(it) ))
00332                                     > 0 ) // same surface point because similar gradient of intensity
00333                                 { find = true ; ; }
00334                             }
00335                             else { find = true; }
00336                             if( find==true && aGrayRange[0]!=InactivateGrayRange && aGrayRange[1]!=InactivateGrayRange  )
00337                             {
00338                                 if( modelBlur.GetPixelFloatValueWithCheck(edge[0], edge[1], edge[2], edge[3]) < aGrayRange[0]
00339                                     || modelBlur.GetPixelFloatValueWithCheck(edge[0], edge[1], edge[2], edge[3]) > aGrayRange[1] )
00340                                     find = false ;
00341                             }
00342                             if( find == true ) break;
00343                         }
00344                     }
00345                 }
00346             }
00347             if( find == true && m_EdgeImage.InsideImage( edge[0], edge[1], edge[2], edge[3] ) ) // the edge
00348             {
00349                 m_EdgeImage.SetPixelByGreyscale( static_cast<int>(edge[0]), static_cast<int>(edge[1]), 
00350                                                  static_cast<int>(edge[2]), it, EdgePoint ) ;
00351                 //set the displacement
00352                 m_pSurfaceBase->SetPixelByGreyscale( ix,iy,iz,it, FoundEdge) ; // have found correct surface point
00353                 for( int id=0; id<3; ++id )
00354                 {
00355                     float dis = (edge[id]-adjustco[id])*sp[id] ;
00356                     m_aDVsurf[id].SetPixelByGreyscale( ix,iy,iz,it, dis ) ;
00357                     m_aDVcor[id].SetPixelByGreyscale( static_cast<int>(edge[0]), static_cast<int>(edge[1]), 
00358                                                       static_cast<int>(edge[2]), static_cast<int>(edge[3]), -dis ) ;
00359                 }
00360                 iPoints ++ ;
00361             }
00362         } // end of surface point check
00363     }
00364     return iPoints ;
00365 }
00366 long int zxhMultiCannyEdge::SmoothBoundary( int nMainTimes, float fNeighbour, float fAccuracy )
00367 {
00374     if( this->m_EdgeImage.GetImageData() == 0 )
00375         return -1 ;
00376     float sp[] = {1,1,1,1} ; m_pSurfaceBase->GetImageSpacing( sp[0],sp[1],sp[2],sp[3] ) ;
00377     if( sp[0] > 0 )
00378         fNeighbour /= sp[0] ;
00379 
00380     int size[] = {1,1,1,1} ;
00381     m_EdgeImage.GetImageSize( size[0],size[1],size[2],size[3] ) ;
00382 
00383     float NV[4] = {0,0,0,0} ;
00384     int iMainPoints = 0 ;
00385     for( int imt=0; imt< nMainTimes; ++imt )
00386     {
00387         for( int ix=0; ix<size[0]; ++ix )
00388         for( int iy=0; iy<size[1]; ++iy )
00389         for( int iz=0; iz<size[2]; ++iz )
00390         for( int it=0; it<size[3]; ++it ) // shape maintanance - Smooth boundary
00391         {
00392             if( m_pSurfaceBase->GetPixelGreyscale( ix,iy,iz,it ) != ZXH_Background )   // for found or not-found
00393             {
00394                 NV[0]=0 ; NV[1]=0 ; NV[2]=0, NV[3]= 0 ;
00395                 int x , y, z , nSurfPnt=0 ;
00396                 for( int nx=-static_cast<int>(fNeighbour); nx<= static_cast<int>(fNeighbour); ++nx )
00397                 for( int ny=-static_cast<int>(fNeighbour); ny<= static_cast<int>(fNeighbour); ++ny )
00398                 for( int nz=-static_cast<int>(fNeighbour); nz<= static_cast<int>(fNeighbour); ++nz ) 
00399                 {
00400                     x= ix+nx ; y= iy+ny ; z= iz+nz ;
00401                     if( m_pSurfaceBase->InsideImage( x, y, z, it )
00402                         && m_pSurfaceBase->GetPixelGreyscale( x, y, z, it ) == FoundEdge ) // find the edge
00403                     {
00404                         NV[0] += m_aDVsurf[0].GetPixelGreyscale( x,y,z,it ) ;
00405                         NV[1] += m_aDVsurf[1].GetPixelGreyscale( x,y,z,it ) ;
00406                         NV[2] += m_aDVsurf[2].GetPixelGreyscale( x,y,z,it ) ;
00407                         ++nSurfPnt ;
00408                     }
00409                 }
00410                 if( nSurfPnt !=0 )
00411                 {
00412                     NV[0] /= nSurfPnt ;
00413                     NV[1] /= nSurfPnt ;
00414                     NV[2] /= nSurfPnt ;
00415                 }
00416                 if(( (NV[0]-m_aDVsurf[0].GetPixelGreyscale( ix,iy,iz,it ))*(NV[0]-m_aDVsurf[0].GetPixelGreyscale( ix,iy,iz,it ))
00417                     +(NV[1]-m_aDVsurf[1].GetPixelGreyscale( ix,iy,iz,it ))*(NV[1]-m_aDVsurf[1].GetPixelGreyscale( ix,iy,iz,it ))
00418                     +(NV[2]-m_aDVsurf[2].GetPixelGreyscale( ix,iy,iz,it ))*(NV[2]-m_aDVsurf[2].GetPixelGreyscale( ix,iy,iz,it ))
00419                     ) > fAccuracy*fAccuracy )
00420                 { // take the neighbour average displacement instead
00421                     iMainPoints ++ ;
00422                     float edge[4] = { float(ix)+ m_aDVsurf[0].GetPixelGreyscale( ix,iy,iz,it )/ sp[0],
00423                                     float(iy)+ m_aDVsurf[1].GetPixelGreyscale( ix,iy,iz,it )/sp[1],
00424                                     float(iz)+ m_aDVsurf[2].GetPixelGreyscale( ix,iy,iz,it )/sp[2], float(it) } ;
00425                     m_EdgeImage.SetPixelByGreyscale( static_cast<int>(edge[0]), static_cast<int>(edge[1]), 
00426                                                      static_cast<int>(edge[2]), it, ZXH_Background ) ;
00427 
00428                     m_pSurfaceBase->SetPixelByGreyscale( ix,iy,iz,it, FoundEdge ) ; // set find the edge
00429                     m_EdgeImage.SetPixelByGreyscale( 
00430                                     ix+static_cast<int>(NV[0]/sp[0]), iy+static_cast<int>(NV[1]/sp[1]),
00431                                     iz+static_cast<int>(NV[2]/sp[2]), it, EdgePoint ) ;
00432                     for( int id=0; id<3; ++id )
00433                     {
00434                         m_aDVsurf[id].SetPixelByGreyscale( ix,iy,iz,it, NV[id] ) ;
00435                         m_aDVcor[id].SetPixelByGreyscale( 
00436                             ix+static_cast<int>(NV[0]/sp[0]), iy+static_cast<int>(NV[1]/sp[1]),
00437                             iz+static_cast<int>(NV[2]/sp[2]), it, -NV[id] ) ;
00438                     }
00439                 }
00440             }
00441         } // end of  shape maintanance - Smooth boundary
00442         //if( glbVerboseOutput>0 ) std::cout<<"success: finish smoothing the " << iMainPoints <<" far corrected surface points and no correct surface points by using neighbour average instead\n" ;
00443     }// end of repeat of shape maintanance
00444     return iMainPoints ;
00445 }
00446 long int    zxhMultiCannyEdge::RemovingJaggedShape( int nMainTimes, float fNeighbour, float fAccuracy )
00447 {
00448     if( this->m_EdgeImage.GetImageData() == 0 )
00449         return -1 ;
00450 
00451     float sp[] = {1,1,1,1} ; m_pSurfaceBase->GetImageSpacing( sp[0],sp[1],sp[2],sp[3] ) ;
00452 
00453     if( sp[0]>0 )
00454         fNeighbour /= sp[0] ;
00455     int size[] = {1,1,1,1};
00456     m_pSurfaceBase->GetImageSize( size[0],size[1],size[2],size[3] ) ;
00457     long int iMainPoints =  0 ;
00458     float NV[ImageDimensionMax] = {0,0,0,0} ;
00459     for( int imt=0; imt< nMainTimes; ++imt )
00460     {
00461         for( int ix=0; ix<size[0]; ++ix )
00462         for( int iy=0; iy<size[1]; ++iy )
00463         for( int iz=0; iz<size[2]; ++iz )
00464         for( int it=0; it<size[3]; ++it ) //removing the far corrected surface points
00465         {
00466             if( m_pSurfaceBase->GetPixelGreyscale( ix,iy,iz,it ) == FoundEdge )
00467             {
00468                 NV[0]=0 ; NV[1]=0 ; NV[2]=0, NV[3]= 0 ;
00469                 int x , y, z , nSurfPnt=0 ;
00470                 for( int nx=-static_cast<int>(fNeighbour); nx<= static_cast<int>(fNeighbour); ++nx )
00471                 for( int ny=-static_cast<int>(fNeighbour); ny<= static_cast<int>(fNeighbour); ++ny )
00472                 for( int nz=-static_cast<int>(fNeighbour); nz<= static_cast<int>(fNeighbour); ++nz )
00473                 {
00474                     x= ix+nx ; y= iy+ny ; z= iz+nz ;
00475                     if( m_pSurfaceBase->InsideImage( x, y, z, it )
00476                         && m_pSurfaceBase->GetPixelGreyscale( x, y, z, it ) == FoundEdge)
00477                     {
00478                         NV[0] += m_aDVsurf[0].GetPixelGreyscale( x,y,z,it ) ;
00479                         NV[1] += m_aDVsurf[1].GetPixelGreyscale( x,y,z,it ) ;
00480                         NV[2] += m_aDVsurf[2].GetPixelGreyscale( x,y,z,it ) ;
00481                         ++nSurfPnt ;
00482                     }
00483                 }
00484                 if( nSurfPnt !=0 )
00485                 {
00486                     NV[0] /= nSurfPnt ;
00487                     NV[1] /= nSurfPnt ;
00488                     NV[2] /= nSurfPnt ;
00489                 }
00490                 if(( (NV[0]-m_aDVsurf[0].GetPixelGreyscale( ix,iy,iz,it ))*(NV[0]-m_aDVsurf[0].GetPixelGreyscale( ix,iy,iz,it ))
00491                     +(NV[1]-m_aDVsurf[1].GetPixelGreyscale( ix,iy,iz,it ))*(NV[1]-m_aDVsurf[1].GetPixelGreyscale( ix,iy,iz,it ))
00492                     +(NV[2]-m_aDVsurf[2].GetPixelGreyscale( ix,iy,iz,it ))*(NV[2]-m_aDVsurf[2].GetPixelGreyscale( ix,iy,iz,it ))
00493                     ) > fAccuracy*fAccuracy )
00494                 {// remove these points as candidates
00495                     iMainPoints ++ ;
00496                     m_pSurfaceBase->SetPixelByGreyscale( ix,iy,iz,it, ZXH_Foreground ) ;
00497                     float edge[4] = { float(ix)+ m_aDVsurf[0].GetPixelGreyscale( ix,iy,iz,it )/sp[0],
00498                         float(iy)+ m_aDVsurf[1].GetPixelGreyscale( ix,iy,iz,it )/sp[1],
00499                         float(iz)+ m_aDVsurf[2].GetPixelGreyscale( ix,iy,iz,it )/sp[2], float(it) } ;
00500                     m_EdgeImage.SetPixelByGreyscale( static_cast<int>(edge[0]), static_cast<int>(edge[1]), 
00501                                                      static_cast<int>(edge[2]), it, ZXH_Background ) ;
00502                     for( int id=0; id<3; ++id )
00503                     {
00504                         m_aDVsurf[id].SetPixelByGreyscale( ix,iy,iz,it, 0 ) ;
00505                         m_aDVcor[id].SetPixelByGreyscale( static_cast<int>(edge[0]), static_cast<int>(edge[1]), 
00506                                                           static_cast<int>(edge[2]), it, 0 ) ;
00507                     }
00508                 }//end of > fAccuracy
00509             }//end of if is surface point
00510         } // end of removing the far corrected surface points
00511     } // end of repeat of nMainTimes
00512     return iMainPoints ;
00513 }
00514 
00515 #endif //ZXHMULTICANNYEDGE_H
00516 
 All Classes Namespaces Functions Variables Typedefs