Orfeo Toolbox  4.0
otbWaveletFilterBank.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12  Copyright (c) Institut Telecom / Telecom Bretagne. All rights reserved.
13  See ITCopyright.txt for details.
14 
15  This software is distributed WITHOUT ANY WARRANTY; without even
16  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
17  PURPOSE. See the above copyright notices for more information.
18 
19 =========================================================================*/
20 
21 #ifndef __otbWaveletFilterBank_txx
22 #define __otbWaveletFilterBank_txx
23 #include "otbWaveletFilterBank.h"
24 
25 #include "otbMacro.h"
28 
31 #include "itkProgressReporter.h"
32 
33 // FIXME
34 #define __myDebug__ 0
35 #include "otbImageFileWriter.h"
36 
37 namespace otb {
38 
43 template <class TInputImage, class TOutputImage, class TWaveletOperator>
46 {
47  this->SetNumberOfRequiredInputs(1);
48  this->SetNumberOfRequiredInputs(1);
49 
50  unsigned int numOfOutputs = 1 << InputImageDimension;
51 
52  this->SetNumberOfRequiredOutputs(numOfOutputs);
53  for (unsigned int i = 0; i < numOfOutputs; ++i)
54  {
55  this->SetNthOutput(i, OutputImageType::New());
56  }
57 
58  m_UpSampleFilterFactor = 0;
59  m_SubsampleImageFactor = 1;
60 
61  //this->SetNumberOfThreads(1);
62 }
63 
64 template <class TInputImage, class TOutputImage, class TWaveletOperator>
65 void
68 {
69  Superclass::GenerateOutputInformation();
70 
71  if (GetSubsampleImageFactor() == 1) return;
72 
73 #if __myDebug__
74  otbGenericMsgDebugMacro(<< " down sampling output regions by a factor of " << GetSubsampleImageFactor());
75  otbGenericMsgDebugMacro(<< "initial region " << this->GetInput()->GetLargestPossibleRegion().GetSize()[0]
76  << "," << this->GetInput()->GetLargestPossibleRegion().GetSize()[1]);
77 #endif
78 
79  OutputImageRegionType newRegion;
80  this->CallCopyInputRegionToOutputRegion(newRegion, this->GetInput()->GetLargestPossibleRegion());
81 
82  for (unsigned int i = 0; i < this->GetNumberOfOutputs(); ++i)
83  {
84  this->GetOutput(i)->SetRegions(newRegion);
85  }
86 
87 #if __myDebug__
88  otbGenericMsgDebugMacro(<< "new region output " << newRegion.GetSize()[0] << "," << newRegion.GetSize()[1]);
89 #endif
90 }
91 
92 template <class TInputImage, class TOutputImage, class TWaveletOperator>
93 void
96 throw (itk::InvalidRequestedRegionError)
97  {
98  Superclass::GenerateInputRequestedRegion();
99 
100  // Filter length calculation
101  LowPassOperatorType lowPassOperator;
102  lowPassOperator.SetDirection(0);
103  lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
104  lowPassOperator.CreateDirectional();
105 
106  unsigned int radius = lowPassOperator.GetRadius()[0];
107 
108  HighPassOperatorType highPassOperator;
109  highPassOperator.SetDirection(0);
110  highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
111  highPassOperator.CreateDirectional();
112 
113  if (radius < highPassOperator.GetRadius()[0]) radius = highPassOperator.GetRadius()[0];
114 
115  // Get the requested region and pad it
116  InputImagePointerType input = const_cast<InputImageType*>(this->GetInput());
117  InputImageRegionType inputRegion = input->GetRequestedRegion();
118  inputRegion.PadByRadius(radius);
119 
120  if (inputRegion.Crop(input->GetLargestPossibleRegion()))
121  {
122  input->SetRequestedRegion(inputRegion);
123  }
124  else
125  {
126  input->SetRequestedRegion(inputRegion);
127  itk::InvalidRequestedRegionError err(__FILE__, __LINE__);
128  err.SetLocation(ITK_LOCATION);
129  err.SetDescription("Requested region is (at least partially) outside the largest possible region.");
130  err.SetDataObject(input);
131  throw err;
132  }
133  }
134 
135 template <class TInputImage, class TOutputImage, class TWaveletOperator>
136 void
139 {
140  if (m_SubsampleImageFactor > 1)
141  {
142  // Check the dimension
143  for (unsigned int i = 0; i < InputImageDimension; ++i)
144  {
145  if ((m_SubsampleImageFactor
146  * (this->GetInput()->GetRequestedRegion().GetSize()[i] / m_SubsampleImageFactor))
147  != this->GetInput()->GetRequestedRegion().GetSize()[i])
148  {
149  itk::InvalidRequestedRegionError err(__FILE__, __LINE__);
150  err.SetLocation(ITK_LOCATION);
151  err.SetDescription("Requested region dimension cannot be used in multiresolution analysis (crop it).");
152  err.SetDataObject(const_cast<InputImageType*>(this->GetInput()));
153  throw err;
154  }
155  }
156 
157  if (InputImageDimension > 1)
158  {
159  // Internal images will be used only if m_SubsampledInputImages != 1
160  m_InternalImages.resize(InputImageDimension - 1);
161  for (unsigned int i = 0; i < m_InternalImages.size(); ++i)
162  {
163  // the size is linked to the SubsampleImageFactor that is assume to be 2!!!
164  m_InternalImages[InputImageDimension - 2 - i].resize(1 << (i + 1));
165  }
166 
167  OutputImageRegionType intermediateRegion;
168  this->Superclass::CallCopyInputRegionToOutputRegion(intermediateRegion,
169  this->GetInput()->GetLargestPossibleRegion());
170 
171  AllocateInternalData(intermediateRegion);
172  }
173  }
174 }
175 
176 template <class TInputImage, class TOutputImage, class TWaveletOperator>
177 void
180  (const OutputImageRegionType& outputRegion)
181 {
182  OutputImageRegionType smallerRegion;
183  OutputImageRegionType largerRegion = outputRegion;
184 
185  for (unsigned int direction = 0; direction < InputImageDimension - 1; direction++)
186  {
187  this->CallCopyInputRegionToOutputRegion(InputImageDimension - 1 - direction,
188  smallerRegion, largerRegion);
189 
190  for (unsigned int i = 0; i < m_InternalImages[direction].size(); ++i)
191  {
192  m_InternalImages[InputImageDimension - 2 - direction][i] = OutputImageType::New();
193  m_InternalImages[InputImageDimension - 2 - direction][i]->SetRegions(smallerRegion);
194  m_InternalImages[InputImageDimension - 2 - direction][i]->Allocate();
195  m_InternalImages[InputImageDimension - 2 - direction][i]->FillBuffer(0);
196  }
197 
198  largerRegion = smallerRegion;
199  }
200 }
201 
202 template <class TInputImage, class TOutputImage, class TWaveletOperator>
203 void
206 {
207  if (m_SubsampleImageFactor > 1 && InputImageDimension > 1)
208  {
209  m_InternalImages.clear();
210  }
211 }
212 
213 template <class TInputImage, class TOutputImage, class TWaveletOperator>
214 void
217  (InputImageRegionType& destRegion, const OutputImageRegionType& srcRegion)
218 {
219  Superclass::CallCopyOutputRegionToInputRegion(destRegion, srcRegion);
220 
221  if (GetSubsampleImageFactor() > 1)
222  {
223  OutputIndexType srcIndex = srcRegion.GetIndex();
224  OutputSizeType srcSize = srcRegion.GetSize();
225 
226  InputIndexType destIndex;
227  InputSizeType destSize;
228 
229  for (unsigned int i = 0; i < InputImageDimension; ++i)
230  {
231  destIndex[i] = srcIndex[i] * GetSubsampleImageFactor();
232  destSize[i] = srcSize[i] * GetSubsampleImageFactor();
233  }
234 
235  destRegion.SetIndex(destIndex);
236  destRegion.SetSize(destSize);
237 
238 #if 0
239  // Contrairement a Wavelet::INVERSE, ici ca ne sera a rien apparemment...
240 
241  // Region Padding
242  LowPassOperatorType lowPassOperator;
243  lowPassOperator.SetDirection(0);
244  lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
245  lowPassOperator.CreateDirectional();
246 
247  unsigned long radius[InputImageDimension];
248  radius[0] = lowPassOperator.GetRadius()[0];
249 
250  HighPassOperatorType highPassOperator;
251  highPassOperator.SetDirection(0);
252  highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
253  highPassOperator.CreateDirectional();
254 
255  if (radius[0] < highPassOperator.GetRadius()[0]) radius[0] = highPassOperator.GetRadius()[0];
256 
257  for (unsigned int i = 1; i < InputImageDimension; ++i)
258  radius[i] = 0;
259 
260  InputImageRegionType paddedRegion = destRegion;
261  paddedRegion.PadByRadius(radius);
262 
263  if (paddedRegion.Crop(this->GetInput()->GetLargestPossibleRegion()))
264  {
265  destRegion = paddedRegion;
266  }
267 #endif
268 
269  }
270 }
271 
272 template <class TInputImage, class TOutputImage, class TWaveletOperator>
273 void
276  (unsigned int direction,
277  InputImageRegionType& destRegion, const OutputImageRegionType& srcRegion)
278 {
279  Superclass::CallCopyOutputRegionToInputRegion(destRegion, srcRegion);
280 
281  if (GetSubsampleImageFactor() > 1)
282  {
283  OutputIndexType srcIndex = srcRegion.GetIndex();
284  OutputSizeType srcSize = srcRegion.GetSize();
285 
286  InputIndexType destIndex;
287  InputSizeType destSize;
288 
289  for (unsigned int i = 0; i < InputImageDimension; ++i)
290  {
291  if (i == direction)
292  {
293  destIndex[i] = srcIndex[i] * GetSubsampleImageFactor();
294  destSize[i] = srcSize[i] * GetSubsampleImageFactor();
295  }
296  else
297  {
298  destIndex[i] = srcIndex[i];
299  destSize[i] = srcSize[i];
300  }
301  }
302 
303  destRegion.SetIndex(destIndex);
304  destRegion.SetSize(destSize);
305  }
306 }
307 
308 template <class TInputImage, class TOutputImage, class TWaveletOperator>
309 void
312  (OutputImageRegionType& destRegion, const InputImageRegionType& srcRegion)
313 {
314  Superclass::CallCopyInputRegionToOutputRegion(destRegion, srcRegion);
315 
316  if (GetSubsampleImageFactor() > 1)
317  {
318  typename InputImageRegionType::IndexType srcIndex = srcRegion.GetIndex();
319  typename InputImageRegionType::SizeType srcSize = srcRegion.GetSize();
320 
321  typename OutputImageRegionType::IndexType destIndex;
322  typename OutputImageRegionType::SizeType destSize;
323 
324  for (unsigned int i = 0; i < InputImageDimension; ++i)
325  {
326  // TODO: This seems not right in odd index cases
327  destIndex[i] = srcIndex[i] / GetSubsampleImageFactor();
328  destSize[i] = srcSize[i] / GetSubsampleImageFactor();
329  }
330 
331  destRegion.SetIndex(destIndex);
332  destRegion.SetSize(destSize);
333  }
334 }
335 
336 template <class TInputImage, class TOutputImage, class TWaveletOperator>
337 void
340  (unsigned int direction,
341  OutputImageRegionType& destRegion, const InputImageRegionType& srcRegion)
342 {
343  Superclass::CallCopyInputRegionToOutputRegion(destRegion, srcRegion);
344 
345  if (GetSubsampleImageFactor() > 1)
346  {
347  typename InputImageRegionType::IndexType srcIndex = srcRegion.GetIndex();
348  typename InputImageRegionType::SizeType srcSize = srcRegion.GetSize();
349 
350  typename OutputImageRegionType::IndexType destIndex;
351  typename OutputImageRegionType::SizeType destSize;
352 
353  for (unsigned int i = 0; i < InputImageDimension; ++i)
354  {
355  if (i == direction)
356  {
357  // TODO: This seems not right in odd index cases
358  destIndex[i] = srcIndex[i] / GetSubsampleImageFactor();
359  destSize[i] = srcSize[i] / GetSubsampleImageFactor();
360  }
361  else
362  {
363  destIndex[i] = srcIndex[i];
364  destSize[i] = srcSize[i];
365  }
366  }
367 
368  destRegion.SetIndex(destIndex);
369  destRegion.SetSize(destSize);
370  }
371 }
372 
373 template <class TInputImage, class TOutputImage, class TWaveletOperator>
374 void
377  (const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
378 {
379  unsigned int dir = InputImageDimension - 1;
380 
381  if ((1 << dir) >= static_cast<int>(this->GetNumberOfOutputs()))
382  {
383  std::ostringstream msg;
384  msg << "Output number 1<<" << dir << " = " << (1 << dir) << " not allocated\n";
385  msg << "Number of excpected outputs " << this->GetNumberOfOutputs() << "\n";
386  throw itk::ExceptionObject(__FILE__, __LINE__, msg.str().c_str(), ITK_LOCATION);
387  }
388 
389  itk::ProgressReporter reporter(this, threadId,
390  outputRegionForThread.GetNumberOfPixels() * this->GetNumberOfOutputs() * 2);
391 
392  const InputImageType * input = this->GetInput();
393  InputImageRegionType inputRegionForThread;
394  this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread);
395 
396  OutputImagePointerType outputLowPass = this->GetOutput(0);
397  OutputImagePointerType outputHighPass = this->GetOutput(1 << dir);
398 
399  // On multidimensional case, if m_SubsampleImageFactor != 1, we need internal images of different size
400  if (dir != 0 && m_SubsampleImageFactor > 1)
401  {
402  outputLowPass = m_InternalImages[dir - 1][0];
403  outputHighPass = m_InternalImages[dir - 1][1];
404  }
405 
406  // typedef for the iterations over the input image
407  typedef itk::ConstNeighborhoodIterator<InputImageType> NeighborhoodIteratorType;
408  typedef itk::NeighborhoodInnerProduct<InputImageType> InnerProductType;
409  typedef itk::ImageRegionIterator<OutputImageType> IteratorType;
410 
411  // Prepare the subsampling image factor, if any.
412  typedef SubsampledImageRegionConstIterator<InputImageType> SubsampleIteratorType;
413  typename SubsampleIteratorType::IndexType subsampling;
414  subsampling.Fill(1);
415  subsampling[dir] = GetSubsampleImageFactor();
416 
417  // Inner product
418  InnerProductType innerProduct;
419 
420  // High pass part calculation
421  HighPassOperatorType highPassOperator;
422  highPassOperator.SetDirection(dir);
423  highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
424  highPassOperator.CreateDirectional();
425 
426  SubsampleIteratorType subItHighPass(input, inputRegionForThread);
427  subItHighPass.SetSubsampleFactor(subsampling);
428  subItHighPass.GoToBegin();
429 
430  NeighborhoodIteratorType itHighPass(highPassOperator.GetRadius(), input, inputRegionForThread);
432  itHighPass.OverrideBoundaryCondition(&boundaryCondition);
433 
434  IteratorType outHighPass(outputHighPass, subItHighPass.GenerateOutputInformation());
435  outHighPass.GoToBegin();
436 
437  while (!subItHighPass.IsAtEnd() && !outHighPass.IsAtEnd())
438  {
439  itHighPass.SetLocation(subItHighPass.GetIndex());
440  outHighPass.Set(innerProduct(itHighPass, highPassOperator));
441 
442  ++subItHighPass;
443  ++outHighPass;
444 
445  reporter.CompletedPixel();
446  }
447 
448  // Low pass part calculation
449  LowPassOperatorType lowPassOperator;
450  lowPassOperator.SetDirection(dir);
451  lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
452  lowPassOperator.CreateDirectional();
453 
454  SubsampleIteratorType subItLowPass(input, inputRegionForThread);
455  subItLowPass.SetSubsampleFactor(subsampling);
456  subItLowPass.GoToBegin();
457 
458  NeighborhoodIteratorType itLowPass(lowPassOperator.GetRadius(), input, inputRegionForThread);
459  itLowPass.OverrideBoundaryCondition(&boundaryCondition);
460 
461  IteratorType outLowPass(outputLowPass, subItLowPass.GenerateOutputInformation());
462  outLowPass.GoToBegin();
463 
464  while (!subItLowPass.IsAtEnd() && !outLowPass.IsAtEnd())
465  {
466  itLowPass.SetLocation(subItLowPass.GetIndex());
467  outLowPass.Set(innerProduct(itLowPass, lowPassOperator));
468 
469  ++subItLowPass;
470  ++outLowPass;
471 
472  reporter.CompletedPixel();
473  }
474 
475  if (dir > 0)
476  {
477  // Note that outputImageRegion correspond to the actual region of (local) input !
478  OutputImageRegionType outputImageRegion;
479  this->CallCopyInputRegionToOutputRegion(dir, outputImageRegion, inputRegionForThread);
480 
481  ThreadedGenerateDataAtDimensionN(0, dir - 1, reporter, outputImageRegion, threadId);
482  ThreadedGenerateDataAtDimensionN(1 << dir, dir - 1, reporter, outputImageRegion, threadId);
483  }
484 }
485 
486 template <class TInputImage, class TOutputImage, class TWaveletOperator>
487 void
490  (unsigned int idx, unsigned int direction, itk::ProgressReporter& reporter,
491  const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
492 {
493  // Note that outputRegionForThread correspond to the actual region of input !
494  OutputImagePointerType input = this->GetOutput(idx);
495  OutputImagePointerType outputHighPass = this->GetOutput(idx + (1 << direction));
496  OutputImagePointerType outputLowPass = OutputImageType::New();
497 
498  OutputImageRegionType outputImageRegion;
499  this->CallCopyInputRegionToOutputRegion(direction, outputImageRegion, outputRegionForThread);
500 
501  if (m_SubsampleImageFactor > 1)
502  {
503  input = m_InternalImages[direction][idx / (1 << (direction + 1))];
504 
505  if (direction != 0)
506  {
507  outputLowPass = m_InternalImages[direction - 1][idx / (1 << direction)];
508  outputHighPass = m_InternalImages[direction - 1][idx / (1 << direction) + 1];
509  }
510  }
511 
512  if (direction == 0)
513  {
514  // The output image has to be allocated
515  // May be not valid on multithreaded process ???
516  outputLowPass->SetRegions(outputImageRegion);
517  outputLowPass->Allocate(); // FIXME Check this line...
518  outputLowPass->FillBuffer(0);
519  }
520 
522  typedef itk::ConstNeighborhoodIterator<OutputImageType> NeighborhoodIteratorType;
523  typedef itk::NeighborhoodInnerProduct<OutputImageType> InnerProductType;
524  typedef itk::ImageRegionIterator<OutputImageType> IteratorType;
525 
526  // Prepare the subsampling image factor, if any.
527  typedef SubsampledImageRegionConstIterator<InputImageType> SubsampleIteratorType;
528  typename SubsampleIteratorType::IndexType subsampling;
529  subsampling.Fill(1);
530  subsampling[direction] = GetSubsampleImageFactor();
531 
532  // Inner products
533  InnerProductType innerProduct;
534 
535  // High pass part calculation
536  HighPassOperatorType highPassOperator;
537  highPassOperator.SetDirection(direction);
538  highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
539  highPassOperator.CreateDirectional();
540 
541  SubsampleIteratorType subItHighPass(input, outputRegionForThread);
542  subItHighPass.SetSubsampleFactor(subsampling);
543  subItHighPass.GoToBegin();
544 
545  NeighborhoodIteratorType itHighPass(highPassOperator.GetRadius(), input, outputRegionForThread);
547  itHighPass.OverrideBoundaryCondition(&boundaryCondition);
548 
549  IteratorType outHighPass(outputHighPass, subItHighPass.GenerateOutputInformation());
550  outHighPass.GoToBegin();
551 
552  while (!subItHighPass.IsAtEnd() && !outHighPass.IsAtEnd())
553  {
554  itHighPass.SetLocation(subItHighPass.GetIndex());
555  outHighPass.Set(innerProduct(itHighPass, highPassOperator));
556 
557  ++subItHighPass;
558  ++outHighPass;
559 
560  reporter.CompletedPixel();
561  }
562 
563  // Low pass part calculation
564  LowPassOperatorType lowPassOperator;
565  lowPassOperator.SetDirection(direction);
566  lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
567  lowPassOperator.CreateDirectional();
568 
569  SubsampleIteratorType subItLowPass(input, outputRegionForThread);
570  subItLowPass.SetSubsampleFactor(subsampling);
571  subItLowPass.GoToBegin();
572 
573  NeighborhoodIteratorType itLowPass(lowPassOperator.GetRadius(), input, outputRegionForThread);
574  itLowPass.OverrideBoundaryCondition(&boundaryCondition);
575 
576  IteratorType outLowPass(outputLowPass, subItLowPass.GenerateOutputInformation());
577  outLowPass.GoToBegin();
578 
579  while (!subItLowPass.IsAtEnd() && !outLowPass.IsAtEnd())
580  {
581  itLowPass.SetLocation(subItLowPass.GetIndex());
582  outLowPass.Set(innerProduct(itLowPass, lowPassOperator));
583 
584  ++subItLowPass;
585  ++outLowPass;
586 
587  reporter.CompletedPixel();
588  }
589 
590  // Swap input and lowPassOutput
591  itk::ImageRegionConstIterator<OutputImageType> inIt(outputLowPass, outputImageRegion);
592  IteratorType outIt((direction != 0 && m_SubsampleImageFactor > 1) ?
593  static_cast<OutputImageType*>(m_InternalImages[direction - 2][idx / (1 << (direction - 1))])
594  : this->GetOutput(idx),
595  outputImageRegion);
596 
597  for (inIt.GoToBegin(), outIt.GoToBegin(); !inIt.IsAtEnd(); ++inIt, ++outIt)
598  {
599  outIt.Set(inIt.Get());
600  }
601 
602  if (direction != 0)
603  {
604  ThreadedGenerateDataAtDimensionN(idx, direction - 1, reporter, outputImageRegion, threadId);
605  ThreadedGenerateDataAtDimensionN(idx + (1 << direction), direction - 1, reporter, outputImageRegion, threadId);
606  }
607 }
608 
613 template <class TInputImage, class TOutputImage, class TWaveletOperator>
616 {
617  this->SetNumberOfRequiredInputs(1 << InputImageDimension);
618 
619  m_UpSampleFilterFactor = 0;
620  m_SubsampleImageFactor = 1;
621 
622  // TODO: For now, we force the number threads to 1 because there is a bug with multithreading in INVERSE transform
623  // Resulting in discontinuities in the reconstructed images
624  this->SetNumberOfThreads(1);
625 }
626 
627 template <class TInputImage, class TOutputImage, class TWaveletOperator>
628 void
631 {
632  Superclass::GenerateOutputInformation();
633 
634  for (unsigned int i = 1; i < this->GetNumberOfInputs(); ++i)
635  {
636  for (unsigned int dim = 0; dim < InputImageDimension; dim++)
637  {
638  if (this->GetInput(0)->GetLargestPossibleRegion().GetSize()[dim]
639  != this->GetInput(i)->GetLargestPossibleRegion().GetSize()[dim])
640  {
641  throw itk::ExceptionObject(__FILE__, __LINE__,
642  "Input images are not of the same dimension", ITK_LOCATION);
643  }
644  }
645  }
646 
647 #if __myDebug__
648  otbGenericMsgDebugMacro(<< " up sampling output regions by a factor of " << GetSubsampleImageFactor());
649 
650  otbGenericMsgDebugMacro(<< "initial region "
651  << this->GetInput(0)->GetLargestPossibleRegion().GetSize()[0]
652  << "," << this->GetInput(0)->GetLargestPossibleRegion().GetSize()[1]);
653 #endif
654 
655  OutputImageRegionType newRegion;
656  this->CallCopyInputRegionToOutputRegion(newRegion, this->GetInput(0)->GetLargestPossibleRegion());
657  this->GetOutput()->SetRegions(newRegion);
658 
659 #if __myDebug__
660  otbGenericMsgDebugMacro(<< "new region output " << newRegion.GetSize()[0] << "," << newRegion.GetSize()[1]);
661 #endif
662 }
663 
664 template <class TInputImage, class TOutputImage, class TWaveletOperator>
665 void
668 throw (itk::InvalidRequestedRegionError)
669  {
670  Superclass::GenerateInputRequestedRegion();
671 
672  // Filter length calculation
673  LowPassOperatorType lowPassOperator;
674  lowPassOperator.SetDirection(0);
675  lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
676  lowPassOperator.CreateDirectional();
677 
678  unsigned int radius = lowPassOperator.GetRadius()[0];
679 
680  HighPassOperatorType highPassOperator;
681  highPassOperator.SetDirection(0);
682  highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
683  highPassOperator.CreateDirectional();
684 
685  if (radius < highPassOperator.GetRadius()[0]) radius = highPassOperator.GetRadius()[0];
686 
687  // Get the requested regionand pad it
688  for (unsigned int idx = 0; idx < this->GetNumberOfInputs(); ++idx)
689  {
690  InputImagePointerType input = const_cast<InputImageType*>(this->GetInput(idx));
691  InputImageRegionType inputRegion = input->GetRequestedRegion();
692  inputRegion.PadByRadius(radius);
693 
694  if (inputRegion.Crop(input->GetLargestPossibleRegion()))
695  {
696  input->SetRequestedRegion(inputRegion);
697  }
698  else
699  {
700  input->SetRequestedRegion(inputRegion);
701  itk::InvalidRequestedRegionError err(__FILE__, __LINE__);
702  err.SetLocation(ITK_LOCATION);
703  err.SetDescription("Requested region is (at least partially) outside the largest possible region.");
704  err.SetDataObject(input);
705  throw err;
706  }
707  }
708  }
709 
710 template <class TInputImage, class TOutputImage, class TWaveletOperator>
711 void
714 {
715  if (InputImageDimension > 1)
716  {
717  // Internal images will be used only if m_SubsampleImageFactor != 1
718  m_InternalImages.resize(InputImageDimension - 1);
719  for (unsigned int i = 0; i < m_InternalImages.size(); ++i)
720  {
721  // the size is linked to the SubsampleImageFactor that is assume to be 2!!!
722  m_InternalImages[i].resize(1 << (i + 1));
723  }
724 
725  OutputImageRegionType intermediateRegion;
726  Superclass::CallCopyInputRegionToOutputRegion(intermediateRegion,
727  this->GetInput(0)->GetLargestPossibleRegion());
728 
729  AllocateInternalData(intermediateRegion);
730  }
731 }
732 
733 template <class TInputImage, class TOutputImage, class TWaveletOperator>
734 void
737  (const OutputImageRegionType& outputRegion)
738 {
739  OutputImageRegionType largerRegion;
740  OutputImageRegionType smallerRegion = outputRegion;
741 
742  for (unsigned int direction = 0; direction < InputImageDimension - 1; direction++)
743  {
744  this->CallCopyInputRegionToOutputRegion(direction,
745  largerRegion, smallerRegion);
746 
747  for (unsigned int i = 0; i < m_InternalImages[direction].size(); ++i)
748  {
749  m_InternalImages[direction][i] = OutputImageType::New();
750  m_InternalImages[direction][i]->SetRegions(largerRegion);
751  m_InternalImages[direction][i]->Allocate();
752  m_InternalImages[direction][i]->FillBuffer(0);
753  }
754 
755  smallerRegion = largerRegion;
756  }
757 }
758 
759 template <class TInputImage, class TOutputImage, class TWaveletOperator>
760 void
763 {
764  if (m_SubsampleImageFactor > 1 && InputImageDimension > 1)
765  {
766  m_InternalImages.clear();
767  }
768 }
769 
770 template <class TInputImage, class TOutputImage, class TWaveletOperator>
771 void
774  (InputImageRegionType& destRegion, const OutputImageRegionType& srcRegion)
775 {
776  Superclass::CallCopyOutputRegionToInputRegion(destRegion, srcRegion);
777 
778  if (GetSubsampleImageFactor() > 1)
779  {
780  OutputIndexType srcIndex = srcRegion.GetIndex();
781  OutputSizeType srcSize = srcRegion.GetSize();
782 
783  InputIndexType destIndex;
784  InputSizeType destSize;
785 
786  for (unsigned int i = 0; i < InputImageDimension; ++i)
787  {
788  // TODO: This seems not right in odd index cases
789  destIndex[i] = srcIndex[i] / GetSubsampleImageFactor();
790  destSize[i] = srcSize[i] / GetSubsampleImageFactor();
791  }
792 
793  destRegion.SetIndex(destIndex);
794  destRegion.SetSize(destSize);
795 
796 #if 1
797  // Region Padding
798  LowPassOperatorType lowPassOperator;
799  lowPassOperator.SetDirection(0);
800  lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
801  lowPassOperator.CreateDirectional();
802 
803  typename InputImageRegionType::SizeType radius;
804  radius[0] = lowPassOperator.GetRadius()[0];
805 
806  HighPassOperatorType highPassOperator;
807  highPassOperator.SetDirection(0);
808  highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
809  highPassOperator.CreateDirectional();
810 
811  if (radius[0] < highPassOperator.GetRadius()[0]) radius[0] = highPassOperator.GetRadius()[0];
812 
813  for (unsigned int i = 1; i < InputImageDimension; ++i)
814  radius[i] = 0;
815 
816 // for ( unsigned int i = 0; i < InputImageDimension; ++i )
817 // {
818 // radius[i] = lowPassOperator.GetRadius()[i];
819 // if ( radius[i] < highPassOperator.GetRadius()[i] )
820 // radius[i] = highPassOperator.GetRadius()[i];
821 // }
822 
823  InputImageRegionType paddedRegion = destRegion;
824  paddedRegion.PadByRadius(radius);
825 
826  if (paddedRegion.Crop(this->GetInput(0)->GetLargestPossibleRegion()))
827  {
828  destRegion = paddedRegion;
829  }
830 #endif
831  }
832 }
833 
834 template <class TInputImage, class TOutputImage, class TWaveletOperator>
835 void
838  (OutputImageRegionType& destRegion, const InputImageRegionType& srcRegion)
839 {
840  Superclass::CallCopyInputRegionToOutputRegion(destRegion, srcRegion);
841 
842  if (GetSubsampleImageFactor() > 1)
843  {
844  OutputIndexType srcIndex = srcRegion.GetIndex();
845  OutputSizeType srcSize = srcRegion.GetSize();
846 
847  InputIndexType destIndex;
848  InputSizeType destSize;
849 
850  for (unsigned int i = 0; i < InputImageDimension; ++i)
851  {
852  destIndex[i] = srcIndex[i] * GetSubsampleImageFactor();
853  destSize[i] = srcSize[i] * GetSubsampleImageFactor();
854  }
855 
856  destRegion.SetIndex(destIndex);
857  destRegion.SetSize(destSize);
858  }
859 }
860 
861 template <class TInputImage, class TOutputImage, class TWaveletOperator>
862 void
865  (unsigned int direction,
866  InputImageRegionType& destRegion, const OutputImageRegionType& srcRegion)
867 {
868  Superclass::CallCopyOutputRegionToInputRegion(destRegion, srcRegion);
869 
870  if (GetSubsampleImageFactor() > 1)
871  {
872  OutputIndexType srcIndex = srcRegion.GetIndex();
873  OutputSizeType srcSize = srcRegion.GetSize();
874 
875  InputIndexType destIndex;
876  InputSizeType destSize;
877 
878  for (unsigned int i = 0; i < InputImageDimension; ++i)
879  {
880  if (i == direction)
881  {
882  // TODO: This seems not right in odd index cases
883  destIndex[i] = srcIndex[i] / GetSubsampleImageFactor();
884  destSize[i] = srcSize[i] / GetSubsampleImageFactor();
885  }
886  else
887  {
888  destIndex[i] = srcIndex[i];
889  destSize[i] = srcSize[i];
890  }
891  }
892 
893  destRegion.SetIndex(destIndex);
894  destRegion.SetSize(destSize);
895  }
896 }
897 
898 template <class TInputImage, class TOutputImage, class TWaveletOperator>
899 void
902  (unsigned int direction,
903  OutputImageRegionType& destRegion, const InputImageRegionType& srcRegion)
904 {
905  Superclass::CallCopyInputRegionToOutputRegion(destRegion, srcRegion);
906 
907  if (GetSubsampleImageFactor() > 1)
908  {
909  typename InputImageRegionType::IndexType srcIndex = srcRegion.GetIndex();
910  typename InputImageRegionType::SizeType srcSize = srcRegion.GetSize();
911 
912  typename OutputImageRegionType::IndexType destIndex;
913  typename OutputImageRegionType::SizeType destSize;
914 
915  for (unsigned int i = 0; i < InputImageDimension; ++i)
916  {
917  if (i == direction)
918  {
919  destIndex[i] = srcIndex[i] * GetSubsampleImageFactor();
920  destSize[i] = srcSize[i] * GetSubsampleImageFactor();
921  }
922  else
923  {
924  destIndex[i] = srcIndex[i];
925  destSize[i] = srcSize[i];
926  }
927  }
928 
929  destRegion.SetIndex(destIndex);
930  destRegion.SetSize(destSize);
931  }
932 }
933 
934 template <class TInputImage, class TOutputImage, class TWaveletOperator>
935 void
938  (const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
939 {
940  itk::ProgressReporter reporter(this, threadId,
941  outputRegionForThread.GetNumberOfPixels() * this->GetNumberOfInputs());
942 
943  InputImageRegionType inputRegionForThread;
944  this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread);
945 
946  unsigned int dir = 0;
947 
948  // Low pass part calculation
949  LowPassOperatorType lowPassOperator;
950  lowPassOperator.SetDirection(dir);
951  lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
952  lowPassOperator.CreateDirectional();
953 
954  // High pass part calculation
955  HighPassOperatorType highPassOperator;
956  highPassOperator.SetDirection(dir);
957  highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
958  highPassOperator.CreateDirectional();
959 
960  // typedef for the iterations over the input image
961  typedef itk::ConstNeighborhoodIterator<OutputImageType> NeighborhoodIteratorType;
962  typedef itk::NeighborhoodInnerProduct<OutputImageType> InnerProductType;
963  typedef itk::ImageRegionIterator<OutputImageType> IteratorType;
964 
965  // Faces iterations
966  typename NeighborhoodIteratorType::RadiusType radiusMax;
967  for (unsigned int idx = 0; idx < OutputImageDimension; ++idx)
968  {
969  radiusMax[idx] = lowPassOperator.GetRadius(idx);
970  if (radiusMax[idx] < highPassOperator.GetRadius(idx)) radiusMax[idx] = highPassOperator.GetRadius(idx);
971  }
972 
973  // The multiresolution case requires a SubsampleImageFilter step
974  if (m_SubsampleImageFactor > 1)
975  {
976  for (unsigned int i = 0; i < this->GetNumberOfInputs(); i += 2)
977  {
978  InputImagePointerType imgLowPass = const_cast<InputImageType*>(this->GetInput(i));
979  InputImagePointerType imgHighPass = const_cast<InputImageType*>(this->GetInput(i + 1));
980 
981  OutputImagePointerType outputImage = this->GetOutput();
982  if (dir != InputImageDimension - 1)
983  {
984  outputImage = m_InternalImages[0][i / 2];
985  }
986 
988  typename FilterType::InputImageIndexType delta;
989  delta.Fill(1);
990  delta[dir] = this->GetSubsampleImageFactor();
991 
992  InputImagePointerType cropedLowPass = InputImageType::New();
993  cropedLowPass->SetRegions(inputRegionForThread);
994  cropedLowPass->Allocate();
995  cropedLowPass->FillBuffer(0.);
996  itk::ImageRegionIterator<InputImageType> cropedLowPassIt(cropedLowPass, inputRegionForThread);
997  itk::ImageRegionIterator<InputImageType> imgLowPassIt(imgLowPass, inputRegionForThread);
998  for (cropedLowPassIt.GoToBegin(), imgLowPassIt.GoToBegin();
999  !cropedLowPassIt.IsAtEnd() && !imgLowPassIt.IsAtEnd();
1000  ++cropedLowPassIt, ++imgLowPassIt)
1001  {
1002  cropedLowPassIt.Set(imgLowPassIt.Get());
1003  }
1004 
1005  typename FilterType::Pointer overSampledLowPass = FilterType::New();
1006  overSampledLowPass->SetInput(cropedLowPass);
1007  overSampledLowPass->SetSubsampleFactor(delta);
1008  overSampledLowPass->Update();
1009 
1010  InputImagePointerType cropedHighPass = InputImageType::New();
1011  cropedHighPass->SetRegions(inputRegionForThread);
1012  cropedHighPass->Allocate();
1013  cropedHighPass->FillBuffer(0.);
1014  itk::ImageRegionIterator<InputImageType> cropedHighPassIt(cropedHighPass, inputRegionForThread);
1015  itk::ImageRegionIterator<InputImageType> imgHighPassIt(imgHighPass, inputRegionForThread);
1016  for (cropedHighPassIt.GoToBegin(), imgHighPassIt.GoToBegin();
1017  !cropedHighPassIt.IsAtEnd() && !imgHighPassIt.IsAtEnd();
1018  ++cropedHighPassIt, ++imgHighPassIt)
1019  {
1020  cropedHighPassIt.Set(imgHighPassIt.Get());
1021  }
1022 
1023  typename FilterType::Pointer overSampledHighPass = FilterType::New();
1024  overSampledHighPass->SetInput(cropedHighPass);
1025  overSampledHighPass->SetSubsampleFactor(delta);
1026  overSampledHighPass->Update();
1027 
1028  InnerProductType innerProduct;
1029 
1031  (outputImage, overSampledLowPass->GetOutput()->GetRequestedRegion());
1032 
1033  NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(),
1034  overSampledLowPass->GetOutput(),
1035  overSampledLowPass->GetOutput()->GetRequestedRegion());
1037  lowIter.OverrideBoundaryCondition(&boundaryCondition);
1038 
1039  NeighborhoodIteratorType highIter(highPassOperator.GetRadius(),
1040  overSampledHighPass->GetOutput(),
1041  overSampledHighPass->GetOutput()->GetRequestedRegion());
1042  highIter.OverrideBoundaryCondition(&boundaryCondition);
1043 
1044  lowIter.GoToBegin();
1045  highIter.GoToBegin();
1046  out.GoToBegin();
1047 
1048  while (!out.IsAtEnd())
1049  {
1050  out.Set(innerProduct(lowIter, lowPassOperator)
1051  + innerProduct(highIter, highPassOperator));
1052 
1053  ++lowIter;
1054  ++highIter;
1055  ++out;
1056 
1057  reporter.CompletedPixel();
1058  }
1059  } // end for each overSampledLowPass/overSampledhighPass pair of entry
1060  }
1061  else // multiscale case
1062  {
1063  for (unsigned int i = 0; i < this->GetNumberOfInputs(); i += 2)
1064  {
1065  InputImagePointerType imgLowPass = const_cast<InputImageType*>(this->GetInput(i));
1066  InputImagePointerType imgHighPass = const_cast<InputImageType*>(this->GetInput(i + 1));
1067 
1068  OutputImagePointerType outputImage = this->GetOutput();
1069  if (dir != InputImageDimension - 1)
1070  {
1071  outputImage = m_InternalImages[0][i / 2];
1072  }
1073 
1074  InnerProductType innerProduct;
1075 
1076  itk::ImageRegionIterator<OutputImageType> out(outputImage, imgLowPass->GetRequestedRegion());
1077 
1078  NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(), imgLowPass, imgLowPass->GetRequestedRegion());
1080  lowIter.OverrideBoundaryCondition(&boundaryCondition);
1081 
1082  NeighborhoodIteratorType highIter(highPassOperator.GetRadius(), imgHighPass, imgLowPass->GetRequestedRegion());
1083  highIter.OverrideBoundaryCondition(&boundaryCondition);
1084 
1085  lowIter.GoToBegin();
1086  highIter.GoToBegin();
1087  out.GoToBegin();
1088 
1089  while (!out.IsAtEnd())
1090  {
1091  out.Set((innerProduct(lowIter, lowPassOperator)
1092  + innerProduct(highIter, highPassOperator)) / 2.);
1093 
1094  ++lowIter;
1095  ++highIter;
1096  ++out;
1097 
1098  reporter.CompletedPixel();
1099  }
1100  } // end for each imgLowPass/imghighPass pair of entry
1101  } // end multiscale case
1102 
1103  if (dir != InputImageDimension - 1)
1104  {
1105  // Note that outputImageRegion correspond to the actual region of (local) input !
1106  OutputImageRegionType outputImageRegion;
1107  this->CallCopyInputRegionToOutputRegion(dir, outputImageRegion, inputRegionForThread);
1108 
1109  ThreadedGenerateDataAtDimensionN(dir + 1, reporter, outputImageRegion, threadId);
1110  }
1111 
1112 }
1113 
1114 template <class TInputImage, class TOutputImage, class TWaveletOperator>
1115 void
1118  (unsigned int direction,
1119  itk::ProgressReporter& reporter,
1120  const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
1121 {
1122  OutputImageRegionType outputImageRegion;
1123  this->CallCopyInputRegionToOutputRegion(direction, outputImageRegion, outputRegionForThread);
1124 
1125  // Low pass part calculation
1126  LowPassOperatorType lowPassOperator;
1127  lowPassOperator.SetDirection(direction);
1128  lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
1129  lowPassOperator.CreateDirectional();
1130 
1131  // High pass part calculation
1132  HighPassOperatorType highPassOperator;
1133  highPassOperator.SetDirection(direction);
1134  highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
1135  highPassOperator.CreateDirectional();
1136 
1137  // typedef for the iterations over the input image
1138  typedef itk::ConstNeighborhoodIterator<OutputImageType> NeighborhoodIteratorType;
1139  typedef itk::NeighborhoodInnerProduct<OutputImageType> InnerProductType;
1140  typedef itk::ImageRegionIterator<OutputImageType> IteratorType;
1141 
1142  // Faces iterations
1143  typename NeighborhoodIteratorType::RadiusType radiusMax;
1144  for (unsigned int i = 0; i < InputImageDimension; ++i)
1145  {
1146  radiusMax[i] = lowPassOperator.GetRadius(i);
1147  if (radiusMax[i] < highPassOperator.GetRadius(i)) radiusMax[i] = highPassOperator.GetRadius(i);
1148  }
1149 
1150  // The multiresolution case requires a SubsampleImageFilter step
1151  if (m_SubsampleImageFactor > 1)
1152  {
1153  for (unsigned int i = 0; i < m_InternalImages[direction - 1].size(); i += 2)
1154  {
1155  InputImagePointerType imgLowPass = m_InternalImages[direction - 1][i];
1156  InputImagePointerType imgHighPass = m_InternalImages[direction - 1][i + 1];
1157 
1158  OutputImagePointerType outputImage = this->GetOutput();
1159  if (direction < InputImageDimension - 1)
1160  {
1161  outputImage = m_InternalImages[direction][i / 2];
1162  }
1163 
1165  typename FilterType::InputImageIndexType delta;
1166  delta.Fill(1);
1167  delta[direction] = this->GetSubsampleImageFactor();
1168 
1169  InputImagePointerType cropedLowPass = InputImageType::New();
1170  cropedLowPass->SetRegions(outputRegionForThread);
1171  cropedLowPass->Allocate();
1172  cropedLowPass->FillBuffer(0.);
1173  itk::ImageRegionIterator<InputImageType> cropedLowPassIt(cropedLowPass, outputRegionForThread);
1174  itk::ImageRegionIterator<InputImageType> imgLowPassIt(imgLowPass, outputRegionForThread);
1175  for (cropedLowPassIt.GoToBegin(), imgLowPassIt.GoToBegin();
1176  !cropedLowPassIt.IsAtEnd() && !imgLowPassIt.IsAtEnd();
1177  ++cropedLowPassIt, ++imgLowPassIt)
1178  {
1179  cropedLowPassIt.Set(imgLowPassIt.Get());
1180  }
1181 
1182  typename FilterType::Pointer overSampledLowPass = FilterType::New();
1183  overSampledLowPass->SetInput(cropedLowPass);
1184  overSampledLowPass->SetSubsampleFactor(delta);
1185  overSampledLowPass->Update();
1186 
1187  InputImagePointerType cropedHighPass = InputImageType::New();
1188  cropedHighPass->SetRegions(outputRegionForThread);
1189  cropedHighPass->Allocate();
1190  cropedHighPass->FillBuffer(0.);
1191  itk::ImageRegionIterator<InputImageType> cropedHighPassIt(cropedHighPass, outputRegionForThread);
1192  itk::ImageRegionIterator<InputImageType> imgHighPassIt(imgHighPass, outputRegionForThread);
1193  for (cropedHighPassIt.GoToBegin(), imgHighPassIt.GoToBegin();
1194  !cropedHighPassIt.IsAtEnd() && !imgHighPassIt.IsAtEnd();
1195  ++cropedHighPassIt, ++imgHighPassIt)
1196  {
1197  cropedHighPassIt.Set(imgHighPassIt.Get());
1198  }
1199 
1200  typename FilterType::Pointer overSampledHighPass = FilterType::New();
1201  overSampledHighPass->SetInput(cropedHighPass);
1202  overSampledHighPass->SetSubsampleFactor(delta);
1203  overSampledHighPass->Update();
1204 
1205  InnerProductType innerProduct;
1206 
1208  overSampledLowPass->GetOutput()->GetRequestedRegion());
1209 
1210  // TODO: This might be the cause of the multithreading bug : we use a neighborhood iterator on cropped data
1211  // Are we sure that we have cropped enough data to access the neighborhood ?
1212  NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(),
1213  overSampledLowPass->GetOutput(),
1214  overSampledLowPass->GetOutput()->GetRequestedRegion());
1216  lowIter.OverrideBoundaryCondition(&boundaryCondition);
1217 
1218  NeighborhoodIteratorType highIter(highPassOperator.GetRadius(),
1219  overSampledHighPass->GetOutput(),
1220  overSampledHighPass->GetOutput()->GetRequestedRegion());
1221  highIter.OverrideBoundaryCondition(&boundaryCondition);
1222 
1223  lowIter.GoToBegin();
1224  highIter.GoToBegin();
1225  out.GoToBegin();
1226 
1227  while (!out.IsAtEnd())
1228  {
1229  out.Set(innerProduct(lowIter, lowPassOperator)
1230  + innerProduct(highIter, highPassOperator));
1231 
1232  ++lowIter;
1233  ++highIter;
1234  ++out;
1235 
1236  reporter.CompletedPixel();
1237  }
1238  } // end for each overSampledLowPass/overSampledhighPass pair of entry
1239  }
1240  else // multiscale case
1241  {
1242  for (unsigned int i = 0; i < m_InternalImages[direction - 1].size(); i += 2)
1243  {
1244  InputImagePointerType imgLowPass = m_InternalImages[direction - 1][i];
1245  InputImagePointerType imgHighPass = m_InternalImages[direction - 1][i + 1];
1246 
1247  OutputImagePointerType outputImage = this->GetOutput();
1248  if (direction < InputImageDimension - 1)
1249  {
1250  outputImage = m_InternalImages[direction][i / 2];
1251  }
1252 
1253  InnerProductType innerProduct;
1254 
1255  itk::ImageRegionIterator<OutputImageType> out(outputImage, imgLowPass->GetRequestedRegion());
1256 
1257  NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(), imgLowPass, imgLowPass->GetRequestedRegion());
1259  lowIter.OverrideBoundaryCondition(&boundaryCondition);
1260 
1261  NeighborhoodIteratorType highIter(highPassOperator.GetRadius(), imgHighPass, imgLowPass->GetRequestedRegion());
1262  highIter.OverrideBoundaryCondition(&boundaryCondition);
1263 
1264  lowIter.GoToBegin();
1265  highIter.GoToBegin();
1266  out.GoToBegin();
1267 
1268  while (!out.IsAtEnd())
1269  {
1270  out.Set((innerProduct(lowIter, lowPassOperator)
1271  + innerProduct(highIter, highPassOperator)) / 2.);
1272 
1273  ++lowIter;
1274  ++highIter;
1275  ++out;
1276 
1277  reporter.CompletedPixel();
1278  }
1279  } // end for each imgLowPass/imghighPass pair of entry
1280  }
1281 
1282  if (direction < InputImageDimension - 1)
1283  {
1284  ThreadedGenerateDataAtDimensionN(direction + 1, reporter, outputImageRegion, threadId);
1285  }
1286 }
1287 
1288 } // end of namespace
1289 
1290 #endif

Generated at Sat Mar 8 2014 16:25:57 for Orfeo Toolbox with doxygen 1.8.3.1