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: 2008-09-29 18:36:37 $
00007   Version:   $Revision: 1.7 $
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 TOutputImage                                        OutputImageType;
00204   typedef typename TOutputImage::PixelType                    OutputPixelType;
00205   typedef typename TInputImage::PixelType                     InputPixelType;
00206    // use the face generator for speed
00207   typedef NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>
00208                                                               FaceCalculatorType;
00209   typedef typename FaceCalculatorType::FaceListType           FaceListType;
00210   typedef typename FaceCalculatorType::FaceListType::iterator FaceListTypeIt;
00211   FaceCalculatorType faceCalculator;
00212 
00213   FaceListType faceList;
00214   FaceListTypeIt fit;
00215   ZeroFluxNeumannBoundaryCondition<TInputImage> nbc;
00216 
00217   // this process is actually slightly asymmetric because we need to
00218   // subtract rectangles that are next to our kernel, not overlapping it
00219   SizeType kernelSize, internalRadius, RegionLimit;
00220   IndexType RegionStart = inputRegion.GetIndex();
00221   for( int i=0; i<TInputImage::ImageDimension; i++ )
00222     {
00223     kernelSize[i] = Radius[i] * 2 + 1;
00224     internalRadius[i] = Radius[i] + 1;
00225     RegionLimit[i] = inputRegion.GetSize()[i] + RegionStart[i] - 1;
00226     }
00227 
00228   typedef typename NumericTraits<OutputPixelType>::RealType AccPixType;
00229   // get a set of offsets to corners for a unit hypercube in this image
00230   std::vector<OffsetType> UnitCorners = CornerOffsets<TInputImage>(accImage);
00231   std::vector<OffsetType> RealCorners;
00232   std::vector<AccPixType> Weights;
00233   // now compute the weights
00234   for (unsigned k = 0; k < UnitCorners.size(); k++)
00235     {
00236     int prod = 1;
00237     OffsetType ThisCorner;
00238     for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00239       {
00240       prod *= UnitCorners[k][i];
00241       if (UnitCorners[k][i] > 0)
00242         {
00243         ThisCorner[i] = Radius[i];
00244         }
00245       else
00246         {
00247         ThisCorner[i] = -(Radius[i]+1);
00248         }
00249       }
00250     Weights.push_back((AccPixType)prod);
00251     RealCorners.push_back(ThisCorner);
00252     }
00253 
00254 
00255   faceList = faceCalculator(accImage, outputRegion, internalRadius);
00256   // start with the body region
00257   for (fit = faceList.begin(); fit != faceList.end(); ++fit)
00258     {
00259     if (fit == faceList.begin())
00260       {
00261       // this is the body region. This is meant to be an optimized
00262       // version that doesn't use neigborhood regions
00263       // compute the various offsets
00264       AccPixType pixelscount = 1;
00265       for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00266         {
00267         pixelscount *= (AccPixType)(2*Radius[i] + 1);
00268         }
00269       
00270       typedef ImageRegionIterator<OutputImageType>     OutputIteratorType;
00271       typedef ImageRegionConstIterator<InputImageType> InputIteratorType;
00272 
00273       typedef std::vector<InputIteratorType> CornerItVecType;
00274       CornerItVecType CornerItVec;
00275       // set up the iterators for each corner
00276       for (unsigned k = 0; k < RealCorners.size(); k++)
00277         {
00278         typename InputImageType::RegionType tReg=(*fit);
00279         tReg.SetIndex(tReg.GetIndex() + RealCorners[k]);
00280         InputIteratorType tempIt(accImage, tReg);
00281         tempIt.GoToBegin();
00282         CornerItVec.push_back(tempIt);
00283         }
00284       // set up the output iterator
00285       OutputIteratorType oIt(outputImage, *fit);
00286       // now do the work
00287       for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00288         {
00289         AccPixType Sum = 0;
00290         // check each corner
00291         for (unsigned k = 0; k < CornerItVec.size(); k++)
00292           {
00293           Sum += Weights[k] * CornerItVec[k].Get();
00294           // increment each corner iterator
00295           ++(CornerItVec[k]);
00296           }
00297         oIt.Set(static_cast<OutputPixelType>(Sum/pixelscount));
00298         progress.CompletedPixel();
00299         }
00300       }
00301     else
00302       {
00303       // now we need to deal with the border regions
00304       typedef ImageRegionIteratorWithIndex<OutputImageType> OutputIteratorType;
00305       OutputIteratorType oIt(outputImage, *fit);
00306       // now do the work
00307       for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00308         {
00309         // figure out the number of pixels in the box by creating an
00310         // equivalent region and cropping - this could probably be
00311         // included in the loop below.
00312         RegionType currentKernelRegion;
00313         currentKernelRegion.SetSize( kernelSize );
00314         // compute the region's index
00315         IndexType kernelRegionIdx = oIt.GetIndex();
00316         IndexType CentIndex = kernelRegionIdx;
00317         for( int i=0; i<TInputImage::ImageDimension; i++ )
00318           {
00319           kernelRegionIdx[i] -= Radius[i];
00320           }
00321         currentKernelRegion.SetIndex( kernelRegionIdx );
00322         currentKernelRegion.Crop( inputRegion );
00323         long edgepixelscount = currentKernelRegion.GetNumberOfPixels();
00324         AccPixType Sum = 0;
00325         // rules are : for each corner,
00326         //               for each dimension
00327         //                  if dimension offset is positive -> this is
00328         //                  a leading edge. Crop if outside the input
00329         //                  region 
00330         //                  if dimension offset is negative -> this is
00331         //                  a trailing edge. Ignore if it is outside
00332         //                  image region
00333         for (unsigned k = 0; k < RealCorners.size(); k++)
00334           {
00335           IndexType ThisCorner = CentIndex + RealCorners[k];
00336           bool IncludeCorner = true;
00337           for (unsigned j = 0; j < TInputImage::ImageDimension; j++)
00338             {
00339             if (UnitCorners[k][j] > 0)
00340               {
00341               // leading edge - crop it
00342               ThisCorner[j] = vnl_math_min(ThisCorner[j], (long)RegionLimit[j]);
00343               }
00344             else
00345               {
00346               // trailing edge - check bounds
00347               if (ThisCorner[j] < RegionStart[j])
00348                 {
00349                 IncludeCorner = false;
00350                 break;
00351                 }
00352               }
00353             }
00354           if (IncludeCorner)
00355             {
00356             Sum += accImage->GetPixel(ThisCorner) * Weights[k];
00357             }
00358           }
00359 
00360         oIt.Set(static_cast<OutputPixelType>(Sum/(AccPixType)edgepixelscount));
00361         progress.CompletedPixel();
00362         }
00363       }
00364     }
00365 }
00366 
00367 
00368 template <class TInputImage, class TOutputImage>
00369 void
00370 BoxSigmaCalculatorFunction(const TInputImage * accImage, 
00371                           TOutputImage * outputImage, 
00372                           typename TInputImage::RegionType inputRegion,
00373                           typename TOutputImage::RegionType outputRegion,
00374                           typename TInputImage::SizeType Radius,
00375                           ProgressReporter &progress)
00376 {
00377   // typedefs
00378   typedef TInputImage                                         InputImageType;
00379   typedef typename TInputImage::RegionType                    RegionType;
00380   typedef typename TInputImage::SizeType                      SizeType;
00381   typedef typename TInputImage::IndexType                     IndexType;
00382   typedef typename TInputImage::PixelType                     PixelType;
00383   typedef typename TInputImage::OffsetType                    OffsetType;
00384   typedef TOutputImage                                        OutputImageType;
00385   typedef typename TOutputImage::PixelType                    OutputPixelType;
00386   typedef typename TInputImage::PixelType                     InputPixelType;
00387   typedef typename InputPixelType::ValueType                  ValueType;
00388    // use the face generator for speed
00389   typedef typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>
00390                                                               FaceCalculatorType;
00391   typedef typename FaceCalculatorType::FaceListType           FaceListType;
00392   typedef typename FaceCalculatorType::FaceListType::iterator FaceListTypeIt;
00393   FaceCalculatorType faceCalculator;
00394 
00395   FaceListType faceList;
00396   FaceListTypeIt fit;
00397   ZeroFluxNeumannBoundaryCondition<TInputImage> nbc;
00398 
00399   // this process is actually slightly asymmetric because we need to
00400   // subtract rectangles that are next to our kernel, not overlapping it
00401   SizeType kernelSize, internalRadius, RegionLimit;
00402   IndexType RegionStart = inputRegion.GetIndex();
00403   for( int i=0; i<TInputImage::ImageDimension; i++ )
00404     {
00405     kernelSize[i] = Radius[i] * 2 + 1;
00406     internalRadius[i] = Radius[i] + 1;
00407     RegionLimit[i] = inputRegion.GetSize()[i] + RegionStart[i] - 1;
00408     }
00409 
00410   typedef typename NumericTraits<OutputPixelType>::RealType AccPixType;
00411   // get a set of offsets to corners for a unit hypercube in this image
00412   std::vector<OffsetType> UnitCorners = CornerOffsets<TInputImage>(accImage);
00413   std::vector<OffsetType> RealCorners;
00414   std::vector<AccPixType> Weights;
00415   // now compute the weights
00416   for (unsigned k = 0; k < UnitCorners.size(); k++)
00417     {
00418     int prod = 1;
00419     OffsetType ThisCorner;
00420     for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00421       {
00422       prod *= UnitCorners[k][i];
00423       if (UnitCorners[k][i] > 0)
00424         {
00425         ThisCorner[i] = Radius[i];
00426         }
00427       else
00428         {
00429         ThisCorner[i] = -(Radius[i]+1);
00430         }
00431       }
00432     Weights.push_back((AccPixType)prod);
00433     RealCorners.push_back(ThisCorner);
00434     }
00435 
00436 
00437   faceList = faceCalculator(accImage, outputRegion, internalRadius);
00438   // start with the body region
00439   for (fit = faceList.begin(); fit != faceList.end(); ++fit)
00440     {
00441     if (fit == faceList.begin())
00442       {
00443       // this is the body region. This is meant to be an optimized
00444       // version that doesn't use neigborhood regions
00445       // compute the various offsets
00446       AccPixType pixelscount = 1;
00447       for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00448         {
00449         pixelscount *= (AccPixType)(2*Radius[i] + 1);
00450         }
00451       
00452       typedef ImageRegionIterator<OutputImageType>     OutputIteratorType;
00453       typedef ImageRegionConstIterator<InputImageType> InputIteratorType;
00454 
00455       typedef std::vector<InputIteratorType> CornerItVecType;
00456       CornerItVecType CornerItVec;
00457       // set up the iterators for each corner
00458       for (unsigned k = 0; k < RealCorners.size(); k++)
00459         {
00460         typename InputImageType::RegionType tReg=(*fit);
00461         tReg.SetIndex(tReg.GetIndex() + RealCorners[k]);
00462         InputIteratorType tempIt(accImage, tReg);
00463         tempIt.GoToBegin();
00464         CornerItVec.push_back(tempIt);
00465         }
00466       // set up the output iterator
00467       OutputIteratorType oIt(outputImage, *fit);
00468       // now do the work
00469       for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00470         {
00471         AccPixType Sum = 0;
00472         AccPixType SquareSum = 0;
00473         // check each corner
00474         for (unsigned k = 0; k < CornerItVec.size(); k++)
00475           {
00476           const InputPixelType & i = CornerItVec[k].Get();
00477           Sum += Weights[k] * i[0];
00478           SquareSum += Weights[k] * i[1];
00479           // increment each corner iterator
00480           ++(CornerItVec[k]);
00481           }
00482 
00483         oIt.Set(static_cast<OutputPixelType>( vcl_sqrt( ( SquareSum - Sum*Sum/pixelscount ) / ( pixelscount -1 ) ) ) );
00484         progress.CompletedPixel();
00485         }
00486       }
00487     else
00488       {
00489       // now we need to deal with the border regions
00490       typedef ImageRegionIteratorWithIndex<OutputImageType> OutputIteratorType;
00491       OutputIteratorType oIt(outputImage, *fit);
00492       // now do the work
00493       for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00494         {
00495         // figure out the number of pixels in the box by creating an
00496         // equivalent region and cropping - this could probably be
00497         // included in the loop below.
00498         RegionType currentKernelRegion;
00499         currentKernelRegion.SetSize( kernelSize );
00500         // compute the region's index
00501         IndexType kernelRegionIdx = oIt.GetIndex();
00502         IndexType CentIndex = kernelRegionIdx;
00503         for( int i=0; i<TInputImage::ImageDimension; i++ )
00504           {
00505           kernelRegionIdx[i] -= Radius[i];
00506           }
00507         currentKernelRegion.SetIndex( kernelRegionIdx );
00508         currentKernelRegion.Crop( inputRegion );
00509         long edgepixelscount = currentKernelRegion.GetNumberOfPixels();
00510         AccPixType Sum = 0;
00511         AccPixType SquareSum = 0;
00512         // rules are : for each corner,
00513         //               for each dimension
00514         //                  if dimension offset is positive -> this is
00515         //                  a leading edge. Crop if outside the input
00516         //                  region 
00517         //                  if dimension offset is negative -> this is
00518         //                  a trailing edge. Ignore if it is outside
00519         //                  image region
00520         for (unsigned k = 0; k < RealCorners.size(); k++)
00521           {
00522           IndexType ThisCorner = CentIndex + RealCorners[k];
00523           bool IncludeCorner = true;
00524           for (unsigned j = 0; j < TInputImage::ImageDimension; j++)
00525             {
00526             if (UnitCorners[k][j] > 0)
00527               {
00528               // leading edge - crop it
00529               ThisCorner[j] = vnl_math_min(ThisCorner[j], (long)RegionLimit[j]);
00530               }
00531             else
00532               {
00533               // trailing edge - check bounds
00534               if (ThisCorner[j] < RegionStart[j])
00535                 {
00536                 IncludeCorner = false;
00537                 break;
00538                 }
00539               }
00540             }
00541           if (IncludeCorner)
00542             {
00543             const InputPixelType & i = accImage->GetPixel(ThisCorner);
00544             Sum += Weights[k] * i[0];
00545             SquareSum += Weights[k] * i[1];
00546             }
00547           }
00548 
00549         oIt.Set(static_cast<OutputPixelType>( vcl_sqrt( ( SquareSum - Sum*Sum/edgepixelscount ) / ( edgepixelscount -1 ) ) ) );
00550         progress.CompletedPixel();
00551         }
00552       }
00553     }
00554 }
00555 
00556 
00557 template <class TInputImage, class TOutputImage>
00558 void
00559 BoxSquareAccumulateFunction(const TInputImage * inputImage, 
00560                       TOutputImage * outputImage, 
00561                       typename TInputImage::RegionType inputRegion,
00562                       typename TOutputImage::RegionType outputRegion,
00563                       ProgressReporter &progress)
00564 {
00565   // typedefs
00566   typedef TInputImage                              InputImageType;
00567   typedef typename TInputImage::RegionType         RegionType;
00568   typedef typename TInputImage::SizeType           SizeType;
00569   typedef typename TInputImage::IndexType          IndexType;
00570   typedef typename TInputImage::PixelType          PixelType;
00571   typedef typename TInputImage::OffsetType         OffsetType;
00572   typedef TOutputImage                             OutputImageType;
00573   typedef typename TOutputImage::PixelType         OutputPixelType;
00574   typedef typename OutputPixelType::ValueType      ValueType;
00575   typedef typename TInputImage::PixelType          InputPixelType;
00576   
00577   typedef ImageRegionConstIterator<TInputImage>    InputIterator;
00578   typedef ImageRegionIterator<TOutputImage>        OutputIterator;
00579 
00580   typedef ShapedNeighborhoodIterator<TOutputImage> NOutputIterator;
00581   InputIterator inIt( inputImage, inputRegion);
00582   typename TInputImage::SizeType kernelRadius;
00583   kernelRadius.Fill(1);
00584 
00585   NOutputIterator noutIt( kernelRadius, outputImage, outputRegion);
00586   // this iterator is fully connected
00587   setConnectivityEarlyBox( &noutIt, true );
00588 
00589   ConstantBoundaryCondition<OutputImageType> oBC;
00590   oBC.SetConstant( NumericTraits< OutputPixelType >::Zero );
00591   noutIt.OverrideBoundaryCondition(&oBC);
00592   // This uses several iterators. An alternative and probably better
00593   // approach would be to copy the input to the output and convolve
00594   // with the following weights (in 2D)
00595   //   -(dim - 1)  1
00596   //       1       1
00597   // The result of each convolution needs to get written back to the
00598   // image being convolved so that the accumulation propogates
00599   // This should be implementable with neighborhood operators.
00600 
00601   std::vector<int> Weights;
00602   typename NOutputIterator::ConstIterator sIt;
00603   for( typename NOutputIterator::IndexListType::const_iterator idxIt = noutIt.GetActiveIndexList().begin();
00604        idxIt != noutIt.GetActiveIndexList().end();
00605        idxIt++ )
00606     {
00607     OffsetType offset = noutIt.GetOffset(*idxIt);
00608     int w = -1;
00609     for (unsigned k = 0; k < InputImageType::ImageDimension; k++)
00610       {
00611       if (offset[k] != 0)
00612         w *= offset[k];
00613       }
00614     Weights.push_back(w);
00615     }
00616 
00617   for (inIt.GoToBegin(), noutIt.GoToBegin(); !noutIt.IsAtEnd(); ++inIt, ++noutIt )
00618     {
00619     ValueType Sum = 0;
00620     ValueType SquareSum = 0;
00621     int k;
00622     for (k = 0, sIt = noutIt.Begin(); !sIt.IsAtEnd();++sIt, ++k)
00623       {
00624       const OutputPixelType & v = sIt.Get();
00625       Sum += v[0] * Weights[k];
00626       SquareSum += v[1] * Weights[k];
00627       }
00628     OutputPixelType o;
00629     const InputPixelType & i = inIt.Get();
00630     o[0] = Sum + i;
00631     o[1] = SquareSum + i*i;
00632     noutIt.SetCenterPixel( o );
00633     progress.CompletedPixel();
00634     }
00635 }
00636 
00637 
00638 } //namespace itk
00639 
00640 #endif
00641 

Generated at Sat Feb 28 12:07:07 2009 for ITK by doxygen 1.5.6 written by Dimitri van Heesch, © 1997-2000