ITK  4.3.0
Insight Segmentation and Registration Toolkit
itkBoxUtilities.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef __itkBoxUtilities_h
19 #define __itkBoxUtilities_h
20 #include "itkProgressReporter.h"
21 #include "itkShapedNeighborhoodIterator.h"
22 #include "itkImageRegionIterator.h"
23 #include "itkConstantBoundaryCondition.h"
24 #include "itkImageRegionIteratorWithIndex.h"
25 #include "itkImageRegionIterator.h"
26 #include "itkOffset.h"
27 #include "itkNeighborhoodAlgorithm.h"
28 #include "itkShapedNeighborhoodIterator.h"
29 #include "itkZeroFluxNeumannBoundaryCondition.h"
30 
31 /*
32  *
33  * This code was contributed in the Insight Journal paper:
34  * "Efficient implementation of kernel filtering"
35  * by Beare R., Lehmann G
36  * http://hdl.handle.net/1926/555
37  * http://www.insight-journal.org/browse/publication/160
38  *
39  */
40 
41 
42 namespace itk
43 {
44 template< class TIterator >
45 TIterator *
46 setConnectivityEarlyBox(TIterator *it, bool fullyConnected = false)
47 {
48  // activate the "previous" neighbours
49  typename TIterator::OffsetType offset;
50  it->ClearActiveList();
51  if ( !fullyConnected )
52  {
53  // only activate the neighbors that are face connected
54  // to the current pixel. do not include the center pixel
55  offset.Fill(0);
56  for ( unsigned int d = 0; d < TIterator::Dimension; ++d )
57  {
58  offset[d] = -1;
59  it->ActivateOffset(offset);
60  offset[d] = 0;
61  }
62  }
63  else
64  {
65  // activate all neighbors that are face+edge+vertex
66  // connected to the current pixel. do not include the center pixel
67  unsigned int centerIndex = it->GetCenterNeighborhoodIndex();
68  for ( unsigned int d = 0; d < centerIndex; d++ )
69  {
70  offset = it->GetOffset(d);
71  // check for positives in any dimension
72  bool keep = true;
73  for ( unsigned int i = 0; i < TIterator::Dimension; i++ )
74  {
75  if ( offset[i] > 0 )
76  {
77  keep = false;
78  break;
79  }
80  }
81  if ( keep )
82  {
83  it->ActivateOffset(offset);
84  }
85  }
86  offset.Fill(0);
87  it->DeactivateOffset(offset);
88  }
89  return it;
90 }
91 
92 template< class TInputImage, class TOutputImage >
93 void
94 BoxAccumulateFunction(const TInputImage *inputImage,
95  const TOutputImage *outputImage,
96  typename TInputImage::RegionType inputRegion,
97  typename TOutputImage::RegionType outputRegion,
98  ProgressReporter & progress)
99 {
100  // typedefs
101  typedef TInputImage InputImageType;
102  typedef typename TInputImage::RegionType RegionType;
103  typedef typename TInputImage::SizeType SizeType;
104  typedef typename TInputImage::IndexType IndexType;
105  typedef typename TInputImage::PixelType PixelType;
106  typedef typename TInputImage::OffsetType OffsetType;
107  typedef TOutputImage OutputImageType;
108  typedef typename TOutputImage::PixelType OutputPixelType;
109  typedef typename TInputImage::PixelType InputPixelType;
110 
111  typedef ImageRegionConstIterator< TInputImage > InputIterator;
112  typedef ImageRegionIterator< TOutputImage > OutputIterator;
113 
114  typedef ShapedNeighborhoodIterator< TOutputImage > NOutputIterator;
115  InputIterator inIt(inputImage, inputRegion);
116  typename TInputImage::SizeType kernelRadius;
117  kernelRadius.Fill(1);
118 
119  NOutputIterator noutIt(kernelRadius, outputImage, outputRegion);
120  // this iterator is fully connected
121  setConnectivityEarlyBox(&noutIt, true);
122 
123  ConstantBoundaryCondition< OutputImageType > oBC;
124  oBC.SetConstant(NumericTraits< OutputPixelType >::Zero);
125  noutIt.OverrideBoundaryCondition(&oBC);
126  // This uses several iterators. An alternative and probably better
127  // approach would be to copy the input to the output and convolve
128  // with the following weights (in 2D)
129  // -(dim - 1) 1
130  // 1 1
131  // The result of each convolution needs to get written back to the
132  // image being convolved so that the accumulation propagates
133  // This should be implementable with neighborhood operators.
134 
135  std::vector< int > Weights;
136  typename NOutputIterator::ConstIterator sIt;
137  for ( typename NOutputIterator::IndexListType::const_iterator idxIt = noutIt.GetActiveIndexList().begin();
138  idxIt != noutIt.GetActiveIndexList().end();
139  idxIt++ )
140  {
141  OffsetType offset = noutIt.GetOffset(*idxIt);
142  int w = -1;
143  for ( unsigned int k = 0; k < InputImageType::ImageDimension; k++ )
144  {
145  if ( offset[k] != 0 )
146  {
147  w *= offset[k];
148  }
149  }
150 // std::cout << offset << " " << w << std::endl;
151  Weights.push_back(w);
152  }
153 
154  for ( inIt.GoToBegin(), noutIt.GoToBegin(); !noutIt.IsAtEnd(); ++inIt, ++noutIt )
155  {
156  OutputPixelType Sum = 0;
157  int k;
158  for ( k = 0, sIt = noutIt.Begin(); !sIt.IsAtEnd(); ++sIt, ++k )
159  {
160  Sum += sIt.Get() * Weights[k];
161  }
162  noutIt.SetCenterPixel( Sum + inIt.Get() );
163  progress.CompletedPixel();
164  }
165 }
166 
167 // a function to generate corners of arbitrary dimension box
168 template< class ImType >
169 std::vector< typename ImType::OffsetType > CornerOffsets(const ImType *im)
170 {
171  typedef ShapedNeighborhoodIterator< ImType > NIterator;
172  typename ImType::SizeType unitradius;
173  unitradius.Fill(1);
174  NIterator N1( unitradius, im, im->GetRequestedRegion() );
175  unsigned int centerIndex = N1.GetCenterNeighborhoodIndex();
176  typename NIterator::OffsetType offset;
177  std::vector< typename ImType::OffsetType > result;
178  for ( unsigned int d = 0; d < centerIndex * 2 + 1; d++ )
179  {
180  offset = N1.GetOffset(d);
181  // check whether this is a corner - corners have no zeros
182  bool corner = true;
183  for ( unsigned int k = 0; k < ImType::ImageDimension; k++ )
184  {
185  if ( offset[k] == 0 )
186  {
187  corner = false;
188  break;
189  }
190  }
191  if ( corner )
192  {
193  result.push_back(offset);
194  }
195  }
196  return ( result );
197 }
198 
199 template< class TInputImage, class TOutputImage >
200 void
201 BoxMeanCalculatorFunction(const TInputImage *accImage,
202  TOutputImage *outputImage,
203  typename TInputImage::RegionType inputRegion,
204  typename TOutputImage::RegionType outputRegion,
205  typename TInputImage::SizeType Radius,
206  ProgressReporter & progress)
207 {
208  // typedefs
209  typedef TInputImage InputImageType;
210  typedef typename TInputImage::RegionType RegionType;
211  typedef typename TInputImage::SizeType SizeType;
212  typedef typename TInputImage::IndexType IndexType;
213  typedef typename TInputImage::PixelType PixelType;
214  typedef typename TInputImage::OffsetType OffsetType;
215  typedef TOutputImage OutputImageType;
216  typedef typename TOutputImage::PixelType OutputPixelType;
217  typedef typename TInputImage::PixelType InputPixelType;
218  // use the face generator for speed
219  typedef NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< InputImageType > FaceCalculatorType;
220  typedef typename FaceCalculatorType::FaceListType FaceListType;
221  typedef typename FaceCalculatorType::FaceListType::iterator FaceListTypeIt;
222  FaceCalculatorType faceCalculator;
223 
224  FaceListType faceList;
225  FaceListTypeIt fit;
226  ZeroFluxNeumannBoundaryCondition< TInputImage > nbc;
227 
228  // this process is actually slightly asymmetric because we need to
229  // subtract rectangles that are next to our kernel, not overlapping it
230  SizeType kernelSize;
231  SizeType internalRadius;
232  SizeType RegionLimit;
233 
234  IndexType RegionStart = inputRegion.GetIndex();
235  for ( unsigned int i = 0; i < TInputImage::ImageDimension; i++ )
236  {
237  kernelSize[i] = Radius[i] * 2 + 1;
238  internalRadius[i] = Radius[i] + 1;
239  RegionLimit[i] = inputRegion.GetSize()[i] + RegionStart[i] - 1;
240  }
241 
242  typedef typename NumericTraits< OutputPixelType >::RealType AccPixType;
243  // get a set of offsets to corners for a unit hypercube in this image
244  std::vector< OffsetType > UnitCorners = CornerOffsets< TInputImage >(accImage);
245  std::vector< OffsetType > RealCorners;
246  std::vector< AccPixType > Weights;
247  // now compute the weights
248  for ( unsigned int k = 0; k < UnitCorners.size(); k++ )
249  {
250  int prod = 1;
251  OffsetType ThisCorner;
252  for ( unsigned int i = 0; i < TInputImage::ImageDimension; i++ )
253  {
254  prod *= UnitCorners[k][i];
255  if ( UnitCorners[k][i] > 0 )
256  {
257  ThisCorner[i] = Radius[i];
258  }
259  else
260  {
261  ThisCorner[i] = -( Radius[i] + 1 );
262  }
263  }
264  Weights.push_back( (AccPixType)prod );
265  RealCorners.push_back(ThisCorner);
266  }
267 
268  faceList = faceCalculator(accImage, outputRegion, internalRadius);
269  // start with the body region
270  for ( fit = faceList.begin(); fit != faceList.end(); ++fit )
271  {
272  if ( fit == faceList.begin() )
273  {
274  // this is the body region. This is meant to be an optimized
275  // version that doesn't use neighborhood regions
276  // compute the various offsets
277  AccPixType pixelscount = 1;
278  for ( unsigned int i = 0; i < TInputImage::ImageDimension; i++ )
279  {
280  pixelscount *= (AccPixType)( 2 * Radius[i] + 1 );
281  }
282 
283  typedef ImageRegionIterator< OutputImageType > OutputIteratorType;
284  typedef ImageRegionConstIterator< InputImageType > InputIteratorType;
285 
286  typedef std::vector< InputIteratorType > CornerItVecType;
287  CornerItVecType CornerItVec;
288  // set up the iterators for each corner
289  for ( unsigned int k = 0; k < RealCorners.size(); k++ )
290  {
291  typename InputImageType::RegionType tReg = ( *fit );
292  tReg.SetIndex(tReg.GetIndex() + RealCorners[k]);
293  InputIteratorType tempIt(accImage, tReg);
294  tempIt.GoToBegin();
295  CornerItVec.push_back(tempIt);
296  }
297  // set up the output iterator
298  OutputIteratorType oIt(outputImage, *fit);
299  // now do the work
300  for ( oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt )
301  {
302  AccPixType Sum = 0;
303  // check each corner
304  for ( unsigned int k = 0; k < CornerItVec.size(); k++ )
305  {
306  Sum += Weights[k] * CornerItVec[k].Get();
307  // increment each corner iterator
308  ++( CornerItVec[k] );
309  }
310  oIt.Set( static_cast< OutputPixelType >( Sum / pixelscount ) );
311  progress.CompletedPixel();
312  }
313  }
314  else
315  {
316  // now we need to deal with the border regions
317  typedef ImageRegionIteratorWithIndex< OutputImageType > OutputIteratorType;
318  OutputIteratorType oIt(outputImage, *fit);
319  // now do the work
320  for ( oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt )
321  {
322  // figure out the number of pixels in the box by creating an
323  // equivalent region and cropping - this could probably be
324  // included in the loop below.
325  RegionType currentKernelRegion;
326  currentKernelRegion.SetSize(kernelSize);
327  // compute the region's index
328  IndexType kernelRegionIdx = oIt.GetIndex();
329  IndexType CentIndex = kernelRegionIdx;
330  for ( unsigned int i = 0; i < TInputImage::ImageDimension; i++ )
331  {
332  kernelRegionIdx[i] -= Radius[i];
333  }
334  currentKernelRegion.SetIndex(kernelRegionIdx);
335  currentKernelRegion.Crop(inputRegion);
336  OffsetValueType edgepixelscount = currentKernelRegion.GetNumberOfPixels();
337  AccPixType Sum = 0;
338  // rules are : for each corner,
339  // for each dimension
340  // if dimension offset is positive -> this is
341  // a leading edge. Crop if outside the input
342  // region
343  // if dimension offset is negative -> this is
344  // a trailing edge. Ignore if it is outside
345  // image region
346  for ( unsigned int k = 0; k < RealCorners.size(); k++ )
347  {
348  IndexType ThisCorner = CentIndex + RealCorners[k];
349  bool IncludeCorner = true;
350  for ( unsigned int j = 0; j < TInputImage::ImageDimension; j++ )
351  {
352  if ( UnitCorners[k][j] > 0 )
353  {
354  // leading edge - crop it
355  if ( ThisCorner[j] > static_cast< OffsetValueType >( RegionLimit[j] ) )
356  {
357  ThisCorner[j] = static_cast< OffsetValueType >( RegionLimit[j] );
358  }
359  }
360  else
361  {
362  // trailing edge - check bounds
363  if ( ThisCorner[j] < RegionStart[j] )
364  {
365  IncludeCorner = false;
366  break;
367  }
368  }
369  }
370  if ( IncludeCorner )
371  {
372  Sum += accImage->GetPixel(ThisCorner) * Weights[k];
373  }
374  }
375 
376  oIt.Set( static_cast< OutputPixelType >( Sum / (AccPixType)edgepixelscount ) );
377  progress.CompletedPixel();
378  }
379  }
380  }
381 }
382 
383 template< class TInputImage, class TOutputImage >
384 void
385 BoxSigmaCalculatorFunction(const TInputImage *accImage,
386  TOutputImage *outputImage,
387  typename TInputImage::RegionType inputRegion,
388  typename TOutputImage::RegionType outputRegion,
389  typename TInputImage::SizeType Radius,
390  ProgressReporter & progress)
391 {
392  // typedefs
393  typedef TInputImage InputImageType;
394  typedef typename TInputImage::RegionType RegionType;
395  typedef typename TInputImage::SizeType SizeType;
396  typedef typename TInputImage::IndexType IndexType;
397  typedef typename TInputImage::PixelType PixelType;
398  typedef typename TInputImage::OffsetType OffsetType;
399  typedef TOutputImage OutputImageType;
400  typedef typename TOutputImage::PixelType OutputPixelType;
401  typedef typename TInputImage::PixelType InputPixelType;
402  typedef typename InputPixelType::ValueType ValueType;
403  // use the face generator for speed
404  typedef typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator< InputImageType > FaceCalculatorType;
405  typedef typename FaceCalculatorType::FaceListType FaceListType;
406  typedef typename FaceCalculatorType::FaceListType::iterator FaceListTypeIt;
407  FaceCalculatorType faceCalculator;
408 
409  FaceListType faceList;
410  FaceListTypeIt fit;
411  ZeroFluxNeumannBoundaryCondition< TInputImage > nbc;
412 
413  // this process is actually slightly asymmetric because we need to
414  // subtract rectangles that are next to our kernel, not overlapping it
415  SizeType kernelSize;
416  SizeType internalRadius;
417  SizeType RegionLimit;
418  IndexType RegionStart = inputRegion.GetIndex();
419  for ( unsigned int i = 0; i < TInputImage::ImageDimension; i++ )
420  {
421  kernelSize[i] = Radius[i] * 2 + 1;
422  internalRadius[i] = Radius[i] + 1;
423  RegionLimit[i] = inputRegion.GetSize()[i] + RegionStart[i] - 1;
424  }
425 
426  typedef typename NumericTraits< OutputPixelType >::RealType AccPixType;
427  // get a set of offsets to corners for a unit hypercube in this image
428  std::vector< OffsetType > UnitCorners = CornerOffsets< TInputImage >(accImage);
429  std::vector< OffsetType > RealCorners;
430  std::vector< AccPixType > Weights;
431  // now compute the weights
432  for ( unsigned int k = 0; k < UnitCorners.size(); k++ )
433  {
434  int prod = 1;
435  OffsetType ThisCorner;
436  for ( unsigned int i = 0; i < TInputImage::ImageDimension; i++ )
437  {
438  prod *= UnitCorners[k][i];
439  if ( UnitCorners[k][i] > 0 )
440  {
441  ThisCorner[i] = Radius[i];
442  }
443  else
444  {
445  ThisCorner[i] = -( Radius[i] + 1 );
446  }
447  }
448  Weights.push_back( (AccPixType)prod );
449  RealCorners.push_back(ThisCorner);
450  }
451 
452  faceList = faceCalculator(accImage, outputRegion, internalRadius);
453  // start with the body region
454  for ( fit = faceList.begin(); fit != faceList.end(); ++fit )
455  {
456  if ( fit == faceList.begin() )
457  {
458  // this is the body region. This is meant to be an optimized
459  // version that doesn't use neighborhood regions
460  // compute the various offsets
461  AccPixType pixelscount = 1;
462  for ( unsigned int i = 0; i < TInputImage::ImageDimension; i++ )
463  {
464  pixelscount *= (AccPixType)( 2 * Radius[i] + 1 );
465  }
466 
467  typedef ImageRegionIterator< OutputImageType > OutputIteratorType;
468  typedef ImageRegionConstIterator< InputImageType > InputIteratorType;
469 
470  typedef std::vector< InputIteratorType > CornerItVecType;
471  CornerItVecType CornerItVec;
472  // set up the iterators for each corner
473  for ( unsigned int k = 0; k < RealCorners.size(); k++ )
474  {
475  typename InputImageType::RegionType tReg = ( *fit );
476  tReg.SetIndex(tReg.GetIndex() + RealCorners[k]);
477  InputIteratorType tempIt(accImage, tReg);
478  tempIt.GoToBegin();
479  CornerItVec.push_back(tempIt);
480  }
481  // set up the output iterator
482  OutputIteratorType oIt(outputImage, *fit);
483  // now do the work
484  for ( oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt )
485  {
486  AccPixType Sum = 0;
487  AccPixType SquareSum = 0;
488  // check each corner
489  for ( unsigned int k = 0; k < CornerItVec.size(); k++ )
490  {
491  const InputPixelType & i = CornerItVec[k].Get();
492  Sum += Weights[k] * i[0];
493  SquareSum += Weights[k] * i[1];
494  // increment each corner iterator
495  ++( CornerItVec[k] );
496  }
497 
498  oIt.Set( static_cast< OutputPixelType >( vcl_sqrt( ( SquareSum - Sum * Sum / pixelscount ) / ( pixelscount - 1 ) ) ) );
499  progress.CompletedPixel();
500  }
501  }
502  else
503  {
504  // now we need to deal with the border regions
505  typedef ImageRegionIteratorWithIndex< OutputImageType > OutputIteratorType;
506  OutputIteratorType oIt(outputImage, *fit);
507  // now do the work
508  for ( oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt )
509  {
510  // figure out the number of pixels in the box by creating an
511  // equivalent region and cropping - this could probably be
512  // included in the loop below.
513  RegionType currentKernelRegion;
514  currentKernelRegion.SetSize(kernelSize);
515  // compute the region's index
516  IndexType kernelRegionIdx = oIt.GetIndex();
517  IndexType CentIndex = kernelRegionIdx;
518  for ( unsigned int i = 0; i < TInputImage::ImageDimension; i++ )
519  {
520  kernelRegionIdx[i] -= Radius[i];
521  }
522  currentKernelRegion.SetIndex(kernelRegionIdx);
523  currentKernelRegion.Crop(inputRegion);
524  SizeValueType edgepixelscount = currentKernelRegion.GetNumberOfPixels();
525  AccPixType Sum = 0;
526  AccPixType SquareSum = 0;
527  // rules are : for each corner,
528  // for each dimension
529  // if dimension offset is positive -> this is
530  // a leading edge. Crop if outside the input
531  // region
532  // if dimension offset is negative -> this is
533  // a trailing edge. Ignore if it is outside
534  // image region
535  for ( unsigned int k = 0; k < RealCorners.size(); k++ )
536  {
537  IndexType ThisCorner = CentIndex + RealCorners[k];
538  bool IncludeCorner = true;
539  for ( unsigned int j = 0; j < TInputImage::ImageDimension; j++ )
540  {
541  if ( UnitCorners[k][j] > 0 )
542  {
543  // leading edge - crop it
544  if ( ThisCorner[j] > static_cast< OffsetValueType >( RegionLimit[j] ) )
545  {
546  ThisCorner[j] = static_cast< OffsetValueType >( RegionLimit[j] );
547  }
548  }
549  else
550  {
551  // trailing edge - check bounds
552  if ( ThisCorner[j] < RegionStart[j] )
553  {
554  IncludeCorner = false;
555  break;
556  }
557  }
558  }
559  if ( IncludeCorner )
560  {
561  const InputPixelType & i = accImage->GetPixel(ThisCorner);
562  Sum += Weights[k] * i[0];
563  SquareSum += Weights[k] * i[1];
564  }
565  }
566 
567  oIt.Set( static_cast< OutputPixelType >( vcl_sqrt( ( SquareSum - Sum * Sum
568  / edgepixelscount ) / ( edgepixelscount - 1 ) ) ) );
569  progress.CompletedPixel();
570  }
571  }
572  }
573 }
574 
575 template< class TInputImage, class TOutputImage >
576 void
577 BoxSquareAccumulateFunction(const TInputImage *inputImage,
578  TOutputImage *outputImage,
579  typename TInputImage::RegionType inputRegion,
580  typename TOutputImage::RegionType outputRegion,
581  ProgressReporter & progress)
582 {
583  // typedefs
584  typedef TInputImage InputImageType;
585  typedef typename TInputImage::RegionType RegionType;
586  typedef typename TInputImage::SizeType SizeType;
587  typedef typename TInputImage::IndexType IndexType;
588  typedef typename TInputImage::PixelType PixelType;
589  typedef typename TInputImage::OffsetType OffsetType;
590  typedef TOutputImage OutputImageType;
591  typedef typename TOutputImage::PixelType OutputPixelType;
592  typedef typename OutputPixelType::ValueType ValueType;
593  typedef typename TInputImage::PixelType InputPixelType;
594 
595  typedef ImageRegionConstIterator< TInputImage > InputIterator;
596  typedef ImageRegionIterator< TOutputImage > OutputIterator;
597 
598  typedef ShapedNeighborhoodIterator< TOutputImage > NOutputIterator;
599  InputIterator inIt(inputImage, inputRegion);
600  typename TInputImage::SizeType kernelRadius;
601  kernelRadius.Fill(1);
602 
603  NOutputIterator noutIt(kernelRadius, outputImage, outputRegion);
604  // this iterator is fully connected
605  setConnectivityEarlyBox(&noutIt, true);
606 
607  ConstantBoundaryCondition< OutputImageType > oBC;
608  oBC.SetConstant(NumericTraits< OutputPixelType >::Zero);
609  noutIt.OverrideBoundaryCondition(&oBC);
610  // This uses several iterators. An alternative and probably better
611  // approach would be to copy the input to the output and convolve
612  // with the following weights (in 2D)
613  // -(dim - 1) 1
614  // 1 1
615  // The result of each convolution needs to get written back to the
616  // image being convolved so that the accumulation propagates
617  // This should be implementable with neighborhood operators.
618 
619  std::vector< int > Weights;
620  typename NOutputIterator::ConstIterator sIt;
621  for ( typename NOutputIterator::IndexListType::const_iterator idxIt = noutIt.GetActiveIndexList().begin();
622  idxIt != noutIt.GetActiveIndexList().end();
623  idxIt++ )
624  {
625  OffsetType offset = noutIt.GetOffset(*idxIt);
626  int w = -1;
627  for ( unsigned int k = 0; k < InputImageType::ImageDimension; k++ )
628  {
629  if ( offset[k] != 0 )
630  {
631  w *= offset[k];
632  }
633  }
634  Weights.push_back(w);
635  }
636 
637  for ( inIt.GoToBegin(), noutIt.GoToBegin(); !noutIt.IsAtEnd(); ++inIt, ++noutIt )
638  {
639  ValueType Sum = 0;
640  ValueType SquareSum = 0;
641  int k;
642  for ( k = 0, sIt = noutIt.Begin(); !sIt.IsAtEnd(); ++sIt, ++k )
643  {
644  const OutputPixelType & v = sIt.Get();
645  Sum += v[0] * Weights[k];
646  SquareSum += v[1] * Weights[k];
647  }
648  OutputPixelType o;
649  const InputPixelType & i = inIt.Get();
650  o[0] = Sum + i;
651  o[1] = SquareSum + i * i;
652  noutIt.SetCenterPixel(o);
653  progress.CompletedPixel();
654  }
655 }
656 } //namespace itk
657 
658 #endif
659 
Tarballs of the nightly generated Doxygen documentation are available for the html, xml, and tag file.
Generated on Sun Dec 9 2012 00:50:12 for ITK by   doxygen 1.8.2