Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

itkBoxUtilities.h

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Program:   Insight Segmentation & Registration Toolkit
00004   Module:    $RCSfile: itkBoxUtilities.h,v $
00005   Language:  C++
00006   Date:      $Date: 2009-07-12 10:52:54 $
00007   Version:   $Revision: 1.8 $
00008 
00009   Copyright (c) Insight Software Consortium. All rights reserved.
00010   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
00011 
00012      This software is distributed WITHOUT ANY WARRANTY; without even 
00013      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
00014      PURPOSE.  See the above copyright notices for more information.
00015 
00016 =========================================================================*/
00017 
00018 #ifndef __itkBoxUtilities_h
00019 #define __itkBoxUtilities_h
00020 #include "itkProgressReporter.h"
00021 #include "itkShapedNeighborhoodIterator.h"
00022 #include "itkImageRegionConstIterator.h"
00023 #include "itkImageRegionIterator.h"
00024 #include "itkConstantBoundaryCondition.h"
00025 #include "itkImageRegionIteratorWithIndex.h"
00026 #include "itkImageRegionIterator.h"
00027 #include "itkImageRegionConstIterator.h"
00028 #include "itkOffset.h"
00029 #include "itkNumericTraits.h"
00030 #include "itkNeighborhoodAlgorithm.h"
00031 #include "itkShapedNeighborhoodIterator.h"
00032 #include "itkZeroFluxNeumannBoundaryCondition.h"
00033 
00034 namespace itk {
00035 
00036 template< class TIterator >
00037 TIterator*
00038 setConnectivityEarlyBox( TIterator* it, bool fullyConnected=false )
00039 {
00040   // activate the "previous" neighbours
00041   typename TIterator::OffsetType offset;
00042   it->ClearActiveList();
00043   if( !fullyConnected) 
00044     {
00045     // only activate the neighbors that are face connected
00046     // to the current pixel. do not include the center pixel
00047     offset.Fill( 0 );
00048     for( unsigned int d=0; d < TIterator::Dimension; ++d )
00049       {
00050       offset[d] = -1;
00051       it->ActivateOffset( offset );
00052       offset[d] = 0;
00053       }
00054     }
00055   else
00056     {
00057     // activate all neighbors that are face+edge+vertex
00058     // connected to the current pixel. do not include the center pixel
00059     unsigned int centerIndex = it->GetCenterNeighborhoodIndex();
00060     for( unsigned int d=0; d < centerIndex; d++ )
00061       {
00062       offset = it->GetOffset( d );
00063       // check for positives in any dimension
00064       bool keep = true;
00065       for (unsigned i = 0; i < TIterator::Dimension; i++)
00066         {
00067         if (offset[i] > 0)
00068           {
00069           keep = false;
00070           break;
00071           }
00072         }
00073       if (keep)
00074         {
00075         it->ActivateOffset( offset );
00076         }
00077       }
00078     offset.Fill(0);
00079     it->DeactivateOffset( offset );
00080     }
00081   return it;
00082 }
00083 
00084 template <class TInputImage, class TOutputImage>
00085 void
00086 BoxAccumulateFunction(const TInputImage * inputImage, 
00087                       const TOutputImage * outputImage, 
00088                       typename TInputImage::RegionType inputRegion,
00089                       typename TOutputImage::RegionType outputRegion,
00090                       ProgressReporter &progress)
00091 {
00092   // typedefs
00093   typedef TInputImage                              InputImageType;
00094   typedef typename TInputImage::RegionType         RegionType;
00095   typedef typename TInputImage::SizeType           SizeType;
00096   typedef typename TInputImage::IndexType          IndexType;
00097   typedef typename TInputImage::PixelType          PixelType;
00098   typedef typename TInputImage::OffsetType         OffsetType;
00099   typedef TOutputImage                             OutputImageType;
00100   typedef typename TOutputImage::PixelType         OutputPixelType;
00101   typedef typename TInputImage::PixelType          InputPixelType;
00102   
00103   typedef ImageRegionConstIterator<TInputImage>    InputIterator;
00104   typedef ImageRegionIterator<TOutputImage>        OutputIterator;
00105 
00106   typedef ShapedNeighborhoodIterator<TOutputImage> NOutputIterator;
00107   InputIterator inIt( inputImage, inputRegion);
00108   typename TInputImage::SizeType kernelRadius;
00109   kernelRadius.Fill(1);
00110 
00111   NOutputIterator noutIt( kernelRadius, outputImage, outputRegion);
00112   // this iterator is fully connected
00113   setConnectivityEarlyBox( &noutIt, true );
00114 
00115   ConstantBoundaryCondition<OutputImageType> oBC;
00116   oBC.SetConstant( NumericTraits< OutputPixelType >::Zero );
00117   noutIt.OverrideBoundaryCondition(&oBC);
00118   // This uses several iterators. An alternative and probably better
00119   // approach would be to copy the input to the output and convolve
00120   // with the following weights (in 2D)
00121   //   -(dim - 1)  1
00122   //       1       1
00123   // The result of each convolution needs to get written back to the
00124   // image being convolved so that the accumulation propogates
00125   // This should be implementable with neighborhood operators.
00126 
00127   std::vector<int> Weights;
00128   typename NOutputIterator::ConstIterator sIt;
00129   for( typename NOutputIterator::IndexListType::const_iterator idxIt = noutIt.GetActiveIndexList().begin();
00130        idxIt != noutIt.GetActiveIndexList().end();
00131        idxIt++ )
00132     {
00133     OffsetType offset = noutIt.GetOffset(*idxIt);
00134     int w = -1;
00135     for (unsigned k = 0; k < InputImageType::ImageDimension; k++)
00136       {
00137       if (offset[k] != 0)
00138         w *= offset[k];
00139       }
00140 //     std::cout << offset << "  " << w << std::endl;
00141     Weights.push_back(w);
00142     }
00143 
00144   for (inIt.GoToBegin(), noutIt.GoToBegin(); !noutIt.IsAtEnd(); ++inIt, ++noutIt )
00145     {
00146     OutputPixelType Sum = 0;
00147     int k;
00148     for (k = 0, sIt = noutIt.Begin(); !sIt.IsAtEnd();++sIt, ++k)
00149       {
00150       Sum += sIt.Get() * Weights[k];
00151       }
00152     noutIt.SetCenterPixel(Sum + inIt.Get());
00153     progress.CompletedPixel();
00154     }
00155 }
00156 
00157 // a function to generate corners of arbitary dimension box
00158 template <class ImType>
00159 std::vector<typename ImType::OffsetType> CornerOffsets(const ImType * im)
00160 {
00161   typedef ShapedNeighborhoodIterator<ImType> NIterator;
00162   typename ImType::SizeType unitradius;
00163   unitradius.Fill(1);
00164   NIterator N1(unitradius, im, im->GetRequestedRegion());
00165   unsigned int centerIndex = N1.GetCenterNeighborhoodIndex();
00166   typename NIterator::OffsetType offset;
00167   std::vector<ITK_TYPENAME ImType::OffsetType> result;
00168   for( unsigned int d=0; d < centerIndex*2 + 1; d++ )
00169     {
00170     offset = N1.GetOffset( d );
00171     // check whether this is a corner - corners have no zeros
00172     bool corner = true;
00173     for (unsigned k = 0; k < ImType::ImageDimension; k++)
00174       {
00175       if (offset[k] == 0)
00176         {
00177         corner = false;
00178         break;
00179         }
00180       }
00181     if (corner)
00182       result.push_back(offset);
00183     }
00184   return(result);
00185 }
00186 
00187 template <class TInputImage, class TOutputImage>
00188 void
00189 BoxMeanCalculatorFunction(const TInputImage * accImage, 
00190                           TOutputImage * outputImage, 
00191                           typename TInputImage::RegionType inputRegion,
00192                           typename TOutputImage::RegionType outputRegion,
00193                           typename TInputImage::SizeType Radius,
00194                           ProgressReporter &progress)
00195 {
00196   // typedefs
00197   typedef TInputImage                                         InputImageType;
00198   typedef typename TInputImage::RegionType                    RegionType;
00199   typedef typename TInputImage::SizeType                      SizeType;
00200   typedef typename TInputImage::IndexType                     IndexType;
00201   typedef typename TInputImage::PixelType                     PixelType;
00202   typedef typename TInputImage::OffsetType                    OffsetType;
00203   typedef typename TInputImage::OffsetValueType               OffsetValueType;
00204   typedef TOutputImage                                        OutputImageType;
00205   typedef typename TOutputImage::PixelType                    OutputPixelType;
00206   typedef typename TInputImage::PixelType                     InputPixelType;
00207    // use the face generator for speed
00208   typedef NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>
00209                                                               FaceCalculatorType;
00210   typedef typename FaceCalculatorType::FaceListType           FaceListType;
00211   typedef typename FaceCalculatorType::FaceListType::iterator FaceListTypeIt;
00212   FaceCalculatorType faceCalculator;
00213 
00214   FaceListType faceList;
00215   FaceListTypeIt fit;
00216   ZeroFluxNeumannBoundaryCondition<TInputImage> nbc;
00217 
00218   // this process is actually slightly asymmetric because we need to
00219   // subtract rectangles that are next to our kernel, not overlapping it
00220   SizeType kernelSize;
00221   SizeType internalRadius;
00222   SizeType RegionLimit;
00223 
00224   IndexType RegionStart = inputRegion.GetIndex();
00225   for( int i=0; i<TInputImage::ImageDimension; i++ )
00226     {
00227     kernelSize[i] = Radius[i] * 2 + 1;
00228     internalRadius[i] = Radius[i] + 1;
00229     RegionLimit[i] = inputRegion.GetSize()[i] + RegionStart[i] - 1;
00230     }
00231 
00232   typedef typename NumericTraits<OutputPixelType>::RealType AccPixType;
00233   // get a set of offsets to corners for a unit hypercube in this image
00234   std::vector<OffsetType> UnitCorners = CornerOffsets<TInputImage>(accImage);
00235   std::vector<OffsetType> RealCorners;
00236   std::vector<AccPixType> Weights;
00237   // now compute the weights
00238   for (unsigned k = 0; k < UnitCorners.size(); k++)
00239     {
00240     int prod = 1;
00241     OffsetType ThisCorner;
00242     for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00243       {
00244       prod *= UnitCorners[k][i];
00245       if (UnitCorners[k][i] > 0)
00246         {
00247         ThisCorner[i] = Radius[i];
00248         }
00249       else
00250         {
00251         ThisCorner[i] = -(Radius[i]+1);
00252         }
00253       }
00254     Weights.push_back((AccPixType)prod);
00255     RealCorners.push_back(ThisCorner);
00256     }
00257 
00258 
00259   faceList = faceCalculator(accImage, outputRegion, internalRadius);
00260   // start with the body region
00261   for (fit = faceList.begin(); fit != faceList.end(); ++fit)
00262     {
00263     if (fit == faceList.begin())
00264       {
00265       // this is the body region. This is meant to be an optimized
00266       // version that doesn't use neigborhood regions
00267       // compute the various offsets
00268       AccPixType pixelscount = 1;
00269       for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00270         {
00271         pixelscount *= (AccPixType)(2*Radius[i] + 1);
00272         }
00273       
00274       typedef ImageRegionIterator<OutputImageType>     OutputIteratorType;
00275       typedef ImageRegionConstIterator<InputImageType> InputIteratorType;
00276 
00277       typedef std::vector<InputIteratorType> CornerItVecType;
00278       CornerItVecType CornerItVec;
00279       // set up the iterators for each corner
00280       for (unsigned k = 0; k < RealCorners.size(); k++)
00281         {
00282         typename InputImageType::RegionType tReg=(*fit);
00283         tReg.SetIndex(tReg.GetIndex() + RealCorners[k]);
00284         InputIteratorType tempIt(accImage, tReg);
00285         tempIt.GoToBegin();
00286         CornerItVec.push_back(tempIt);
00287         }
00288       // set up the output iterator
00289       OutputIteratorType oIt(outputImage, *fit);
00290       // now do the work
00291       for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00292         {
00293         AccPixType Sum = 0;
00294         // check each corner
00295         for (unsigned k = 0; k < CornerItVec.size(); k++)
00296           {
00297           Sum += Weights[k] * CornerItVec[k].Get();
00298           // increment each corner iterator
00299           ++(CornerItVec[k]);
00300           }
00301         oIt.Set(static_cast<OutputPixelType>(Sum/pixelscount));
00302         progress.CompletedPixel();
00303         }
00304       }
00305     else
00306       {
00307       // now we need to deal with the border regions
00308       typedef ImageRegionIteratorWithIndex<OutputImageType> OutputIteratorType;
00309       OutputIteratorType oIt(outputImage, *fit);
00310       // now do the work
00311       for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00312         {
00313         // figure out the number of pixels in the box by creating an
00314         // equivalent region and cropping - this could probably be
00315         // included in the loop below.
00316         RegionType currentKernelRegion;
00317         currentKernelRegion.SetSize( kernelSize );
00318         // compute the region's index
00319         IndexType kernelRegionIdx = oIt.GetIndex();
00320         IndexType CentIndex = kernelRegionIdx;
00321         for( int i=0; i<TInputImage::ImageDimension; i++ )
00322           {
00323           kernelRegionIdx[i] -= Radius[i];
00324           }
00325         currentKernelRegion.SetIndex( kernelRegionIdx );
00326         currentKernelRegion.Crop( inputRegion );
00327         long edgepixelscount = currentKernelRegion.GetNumberOfPixels();
00328         AccPixType Sum = 0;
00329         // rules are : for each corner,
00330         //               for each dimension
00331         //                  if dimension offset is positive -> this is
00332         //                  a leading edge. Crop if outside the input
00333         //                  region 
00334         //                  if dimension offset is negative -> this is
00335         //                  a trailing edge. Ignore if it is outside
00336         //                  image region
00337         for (unsigned k = 0; k < RealCorners.size(); k++)
00338           {
00339           IndexType ThisCorner = CentIndex + RealCorners[k];
00340           bool IncludeCorner = true;
00341           for (unsigned j = 0; j < TInputImage::ImageDimension; j++)
00342             {
00343             if (UnitCorners[k][j] > 0)
00344               {
00345               // leading edge - crop it
00346               if( ThisCorner[j] > static_cast< OffsetValueType>( RegionLimit[j] ) )
00347                 {
00348                 ThisCorner[j] = static_cast< OffsetValueType>( RegionLimit[j] );
00349                 }
00350               }
00351             else
00352               {
00353               // trailing edge - check bounds
00354               if (ThisCorner[j] < RegionStart[j])
00355                 {
00356                 IncludeCorner = false;
00357                 break;
00358                 }
00359               }
00360             }
00361           if (IncludeCorner)
00362             {
00363             Sum += accImage->GetPixel(ThisCorner) * Weights[k];
00364             }
00365           }
00366 
00367         oIt.Set(static_cast<OutputPixelType>(Sum/(AccPixType)edgepixelscount));
00368         progress.CompletedPixel();
00369         }
00370       }
00371     }
00372 }
00373 
00374 
00375 template <class TInputImage, class TOutputImage>
00376 void
00377 BoxSigmaCalculatorFunction(const TInputImage * accImage, 
00378                           TOutputImage * outputImage, 
00379                           typename TInputImage::RegionType inputRegion,
00380                           typename TOutputImage::RegionType outputRegion,
00381                           typename TInputImage::SizeType Radius,
00382                           ProgressReporter &progress)
00383 {
00384   // typedefs
00385   typedef TInputImage                                         InputImageType;
00386   typedef typename TInputImage::RegionType                    RegionType;
00387   typedef typename TInputImage::SizeType                      SizeType;
00388   typedef typename TInputImage::IndexType                     IndexType;
00389   typedef typename TInputImage::PixelType                     PixelType;
00390   typedef typename TInputImage::OffsetType                    OffsetType;
00391   typedef typename TInputImage::OffsetValueType               OffsetValueType;
00392   typedef TOutputImage                                        OutputImageType;
00393   typedef typename TOutputImage::PixelType                    OutputPixelType;
00394   typedef typename TInputImage::PixelType                     InputPixelType;
00395   typedef typename InputPixelType::ValueType                  ValueType;
00396    // use the face generator for speed
00397   typedef typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>
00398                                                               FaceCalculatorType;
00399   typedef typename FaceCalculatorType::FaceListType           FaceListType;
00400   typedef typename FaceCalculatorType::FaceListType::iterator FaceListTypeIt;
00401   FaceCalculatorType faceCalculator;
00402 
00403   FaceListType faceList;
00404   FaceListTypeIt fit;
00405   ZeroFluxNeumannBoundaryCondition<TInputImage> nbc;
00406 
00407   // this process is actually slightly asymmetric because we need to
00408   // subtract rectangles that are next to our kernel, not overlapping it
00409   SizeType kernelSize;
00410   SizeType internalRadius;
00411   SizeType RegionLimit;
00412   IndexType RegionStart = inputRegion.GetIndex();
00413   for( int i=0; i<TInputImage::ImageDimension; i++ )
00414     {
00415     kernelSize[i] = Radius[i] * 2 + 1;
00416     internalRadius[i] = Radius[i] + 1;
00417     RegionLimit[i] = inputRegion.GetSize()[i] + RegionStart[i] - 1;
00418     }
00419 
00420   typedef typename NumericTraits<OutputPixelType>::RealType AccPixType;
00421   // get a set of offsets to corners for a unit hypercube in this image
00422   std::vector<OffsetType> UnitCorners = CornerOffsets<TInputImage>(accImage);
00423   std::vector<OffsetType> RealCorners;
00424   std::vector<AccPixType> Weights;
00425   // now compute the weights
00426   for (unsigned k = 0; k < UnitCorners.size(); k++)
00427     {
00428     int prod = 1;
00429     OffsetType ThisCorner;
00430     for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00431       {
00432       prod *= UnitCorners[k][i];
00433       if (UnitCorners[k][i] > 0)
00434         {
00435         ThisCorner[i] = Radius[i];
00436         }
00437       else
00438         {
00439         ThisCorner[i] = -(Radius[i]+1);
00440         }
00441       }
00442     Weights.push_back((AccPixType)prod);
00443     RealCorners.push_back(ThisCorner);
00444     }
00445 
00446 
00447   faceList = faceCalculator(accImage, outputRegion, internalRadius);
00448   // start with the body region
00449   for (fit = faceList.begin(); fit != faceList.end(); ++fit)
00450     {
00451     if (fit == faceList.begin())
00452       {
00453       // this is the body region. This is meant to be an optimized
00454       // version that doesn't use neigborhood regions
00455       // compute the various offsets
00456       AccPixType pixelscount = 1;
00457       for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00458         {
00459         pixelscount *= (AccPixType)(2*Radius[i] + 1);
00460         }
00461       
00462       typedef ImageRegionIterator<OutputImageType>     OutputIteratorType;
00463       typedef ImageRegionConstIterator<InputImageType> InputIteratorType;
00464 
00465       typedef std::vector<InputIteratorType> CornerItVecType;
00466       CornerItVecType CornerItVec;
00467       // set up the iterators for each corner
00468       for (unsigned k = 0; k < RealCorners.size(); k++)
00469         {
00470         typename InputImageType::RegionType tReg=(*fit);
00471         tReg.SetIndex(tReg.GetIndex() + RealCorners[k]);
00472         InputIteratorType tempIt(accImage, tReg);
00473         tempIt.GoToBegin();
00474         CornerItVec.push_back(tempIt);
00475         }
00476       // set up the output iterator
00477       OutputIteratorType oIt(outputImage, *fit);
00478       // now do the work
00479       for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00480         {
00481         AccPixType Sum = 0;
00482         AccPixType SquareSum = 0;
00483         // check each corner
00484         for (unsigned k = 0; k < CornerItVec.size(); k++)
00485           {
00486           const InputPixelType & i = CornerItVec[k].Get();
00487           Sum += Weights[k] * i[0];
00488           SquareSum += Weights[k] * i[1];
00489           // increment each corner iterator
00490           ++(CornerItVec[k]);
00491           }
00492 
00493         oIt.Set(static_cast<OutputPixelType>( vcl_sqrt( ( SquareSum - Sum*Sum/pixelscount ) / ( pixelscount -1 ) ) ) );
00494         progress.CompletedPixel();
00495         }
00496       }
00497     else
00498       {
00499       // now we need to deal with the border regions
00500       typedef ImageRegionIteratorWithIndex<OutputImageType> OutputIteratorType;
00501       OutputIteratorType oIt(outputImage, *fit);
00502       // now do the work
00503       for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00504         {
00505         // figure out the number of pixels in the box by creating an
00506         // equivalent region and cropping - this could probably be
00507         // included in the loop below.
00508         RegionType currentKernelRegion;
00509         currentKernelRegion.SetSize( kernelSize );
00510         // compute the region's index
00511         IndexType kernelRegionIdx = oIt.GetIndex();
00512         IndexType CentIndex = kernelRegionIdx;
00513         for( int i=0; i<TInputImage::ImageDimension; i++ )
00514           {
00515           kernelRegionIdx[i] -= Radius[i];
00516           }
00517         currentKernelRegion.SetIndex( kernelRegionIdx );
00518         currentKernelRegion.Crop( inputRegion );
00519         long edgepixelscount = currentKernelRegion.GetNumberOfPixels();
00520         AccPixType Sum = 0;
00521         AccPixType SquareSum = 0;
00522         // rules are : for each corner,
00523         //               for each dimension
00524         //                  if dimension offset is positive -> this is
00525         //                  a leading edge. Crop if outside the input
00526         //                  region 
00527         //                  if dimension offset is negative -> this is
00528         //                  a trailing edge. Ignore if it is outside
00529         //                  image region
00530         for (unsigned k = 0; k < RealCorners.size(); k++)
00531           {
00532           IndexType ThisCorner = CentIndex + RealCorners[k];
00533           bool IncludeCorner = true;
00534           for (unsigned j = 0; j < TInputImage::ImageDimension; j++)
00535             {
00536             if (UnitCorners[k][j] > 0)
00537               {
00538               // leading edge - crop it
00539               if( ThisCorner[j] > static_cast< OffsetValueType >( RegionLimit[j]) )
00540                 {
00541                 ThisCorner[j] = static_cast< OffsetValueType>( RegionLimit[j] );
00542                 }
00543               }
00544             else
00545               {
00546               // trailing edge - check bounds
00547               if (ThisCorner[j] < RegionStart[j])
00548                 {
00549                 IncludeCorner = false;
00550                 break;
00551                 }
00552               }
00553             }
00554           if (IncludeCorner)
00555             {
00556             const InputPixelType & i = accImage->GetPixel(ThisCorner);
00557             Sum += Weights[k] * i[0];
00558             SquareSum += Weights[k] * i[1];
00559             }
00560           }
00561 
00562         oIt.Set(static_cast<OutputPixelType>( vcl_sqrt( ( SquareSum - Sum*Sum/edgepixelscount ) / ( edgepixelscount -1 ) ) ) );
00563         progress.CompletedPixel();
00564         }
00565       }
00566     }
00567 }
00568 
00569 
00570 template <class TInputImage, class TOutputImage>
00571 void
00572 BoxSquareAccumulateFunction(const TInputImage * inputImage, 
00573                       TOutputImage * outputImage, 
00574                       typename TInputImage::RegionType inputRegion,
00575                       typename TOutputImage::RegionType outputRegion,
00576                       ProgressReporter &progress)
00577 {
00578   // typedefs
00579   typedef TInputImage                              InputImageType;
00580   typedef typename TInputImage::RegionType         RegionType;
00581   typedef typename TInputImage::SizeType           SizeType;
00582   typedef typename TInputImage::IndexType          IndexType;
00583   typedef typename TInputImage::PixelType          PixelType;
00584   typedef typename TInputImage::OffsetType         OffsetType;
00585   typedef TOutputImage                             OutputImageType;
00586   typedef typename TOutputImage::PixelType         OutputPixelType;
00587   typedef typename OutputPixelType::ValueType      ValueType;
00588   typedef typename TInputImage::PixelType          InputPixelType;
00589   
00590   typedef ImageRegionConstIterator<TInputImage>    InputIterator;
00591   typedef ImageRegionIterator<TOutputImage>        OutputIterator;
00592 
00593   typedef ShapedNeighborhoodIterator<TOutputImage> NOutputIterator;
00594   InputIterator inIt( inputImage, inputRegion);
00595   typename TInputImage::SizeType kernelRadius;
00596   kernelRadius.Fill(1);
00597 
00598   NOutputIterator noutIt( kernelRadius, outputImage, outputRegion);
00599   // this iterator is fully connected
00600   setConnectivityEarlyBox( &noutIt, true );
00601 
00602   ConstantBoundaryCondition<OutputImageType> oBC;
00603   oBC.SetConstant( NumericTraits< OutputPixelType >::Zero );
00604   noutIt.OverrideBoundaryCondition(&oBC);
00605   // This uses several iterators. An alternative and probably better
00606   // approach would be to copy the input to the output and convolve
00607   // with the following weights (in 2D)
00608   //   -(dim - 1)  1
00609   //       1       1
00610   // The result of each convolution needs to get written back to the
00611   // image being convolved so that the accumulation propogates
00612   // This should be implementable with neighborhood operators.
00613 
00614   std::vector<int> Weights;
00615   typename NOutputIterator::ConstIterator sIt;
00616   for( typename NOutputIterator::IndexListType::const_iterator idxIt = noutIt.GetActiveIndexList().begin();
00617        idxIt != noutIt.GetActiveIndexList().end();
00618        idxIt++ )
00619     {
00620     OffsetType offset = noutIt.GetOffset(*idxIt);
00621     int w = -1;
00622     for (unsigned k = 0; k < InputImageType::ImageDimension; k++)
00623       {
00624       if (offset[k] != 0)
00625         w *= offset[k];
00626       }
00627     Weights.push_back(w);
00628     }
00629 
00630   for (inIt.GoToBegin(), noutIt.GoToBegin(); !noutIt.IsAtEnd(); ++inIt, ++noutIt )
00631     {
00632     ValueType Sum = 0;
00633     ValueType SquareSum = 0;
00634     int k;
00635     for (k = 0, sIt = noutIt.Begin(); !sIt.IsAtEnd();++sIt, ++k)
00636       {
00637       const OutputPixelType & v = sIt.Get();
00638       Sum += v[0] * Weights[k];
00639       SquareSum += v[1] * Weights[k];
00640       }
00641     OutputPixelType o;
00642     const InputPixelType & i = inIt.Get();
00643     o[0] = Sum + i;
00644     o[1] = SquareSum + i*i;
00645     noutIt.SetCenterPixel( o );
00646     progress.CompletedPixel();
00647     }
00648 }
00649 
00650 
00651 } //namespace itk
00652 
00653 #endif
00654 

Generated at Tue Sep 15 02:26:18 2009 for ITK by doxygen 1.5.8 written by Dimitri van Heesch, © 1997-2000