Orfeo Toolbox  4.2
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 
964  // Faces iterations
965  typename NeighborhoodIteratorType::RadiusType radiusMax;
966  for (unsigned int idx = 0; idx < OutputImageDimension; ++idx)
967  {
968  radiusMax[idx] = lowPassOperator.GetRadius(idx);
969  if (radiusMax[idx] < highPassOperator.GetRadius(idx)) radiusMax[idx] = highPassOperator.GetRadius(idx);
970  }
971 
972  // The multiresolution case requires a SubsampleImageFilter step
973  if (m_SubsampleImageFactor > 1)
974  {
975  for (unsigned int i = 0; i < this->GetNumberOfInputs(); i += 2)
976  {
977  InputImagePointerType imgLowPass = const_cast<InputImageType*>(this->GetInput(i));
978  InputImagePointerType imgHighPass = const_cast<InputImageType*>(this->GetInput(i + 1));
979 
980  OutputImagePointerType outputImage = this->GetOutput();
981  if (dir != InputImageDimension - 1)
982  {
983  outputImage = m_InternalImages[0][i / 2];
984  }
985 
987  typename FilterType::InputImageIndexType delta;
988  delta.Fill(1);
989  delta[dir] = this->GetSubsampleImageFactor();
990 
991  InputImagePointerType cropedLowPass = InputImageType::New();
992  cropedLowPass->SetRegions(inputRegionForThread);
993  cropedLowPass->Allocate();
994  cropedLowPass->FillBuffer(0.);
995  itk::ImageRegionIterator<InputImageType> cropedLowPassIt(cropedLowPass, inputRegionForThread);
996  itk::ImageRegionIterator<InputImageType> imgLowPassIt(imgLowPass, inputRegionForThread);
997  for (cropedLowPassIt.GoToBegin(), imgLowPassIt.GoToBegin();
998  !cropedLowPassIt.IsAtEnd() && !imgLowPassIt.IsAtEnd();
999  ++cropedLowPassIt, ++imgLowPassIt)
1000  {
1001  cropedLowPassIt.Set(imgLowPassIt.Get());
1002  }
1003 
1004  typename FilterType::Pointer overSampledLowPass = FilterType::New();
1005  overSampledLowPass->SetInput(cropedLowPass);
1006  overSampledLowPass->SetSubsampleFactor(delta);
1007  overSampledLowPass->Update();
1008 
1009  InputImagePointerType cropedHighPass = InputImageType::New();
1010  cropedHighPass->SetRegions(inputRegionForThread);
1011  cropedHighPass->Allocate();
1012  cropedHighPass->FillBuffer(0.);
1013  itk::ImageRegionIterator<InputImageType> cropedHighPassIt(cropedHighPass, inputRegionForThread);
1014  itk::ImageRegionIterator<InputImageType> imgHighPassIt(imgHighPass, inputRegionForThread);
1015  for (cropedHighPassIt.GoToBegin(), imgHighPassIt.GoToBegin();
1016  !cropedHighPassIt.IsAtEnd() && !imgHighPassIt.IsAtEnd();
1017  ++cropedHighPassIt, ++imgHighPassIt)
1018  {
1019  cropedHighPassIt.Set(imgHighPassIt.Get());
1020  }
1021 
1022  typename FilterType::Pointer overSampledHighPass = FilterType::New();
1023  overSampledHighPass->SetInput(cropedHighPass);
1024  overSampledHighPass->SetSubsampleFactor(delta);
1025  overSampledHighPass->Update();
1026 
1027  InnerProductType innerProduct;
1028 
1030  (outputImage, overSampledLowPass->GetOutput()->GetRequestedRegion());
1031 
1032  NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(),
1033  overSampledLowPass->GetOutput(),
1034  overSampledLowPass->GetOutput()->GetRequestedRegion());
1036  lowIter.OverrideBoundaryCondition(&boundaryCondition);
1037 
1038  NeighborhoodIteratorType highIter(highPassOperator.GetRadius(),
1039  overSampledHighPass->GetOutput(),
1040  overSampledHighPass->GetOutput()->GetRequestedRegion());
1041  highIter.OverrideBoundaryCondition(&boundaryCondition);
1042 
1043  lowIter.GoToBegin();
1044  highIter.GoToBegin();
1045  out.GoToBegin();
1046 
1047  while (!out.IsAtEnd())
1048  {
1049  out.Set(innerProduct(lowIter, lowPassOperator)
1050  + innerProduct(highIter, highPassOperator));
1051 
1052  ++lowIter;
1053  ++highIter;
1054  ++out;
1055 
1056  reporter.CompletedPixel();
1057  }
1058  } // end for each overSampledLowPass/overSampledhighPass pair of entry
1059  }
1060  else // multiscale case
1061  {
1062  for (unsigned int i = 0; i < this->GetNumberOfInputs(); i += 2)
1063  {
1064  InputImagePointerType imgLowPass = const_cast<InputImageType*>(this->GetInput(i));
1065  InputImagePointerType imgHighPass = const_cast<InputImageType*>(this->GetInput(i + 1));
1066 
1067  OutputImagePointerType outputImage = this->GetOutput();
1068  if (dir != InputImageDimension - 1)
1069  {
1070  outputImage = m_InternalImages[0][i / 2];
1071  }
1072 
1073  InnerProductType innerProduct;
1074 
1075  itk::ImageRegionIterator<OutputImageType> out(outputImage, imgLowPass->GetRequestedRegion());
1076 
1077  NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(), imgLowPass, imgLowPass->GetRequestedRegion());
1079  lowIter.OverrideBoundaryCondition(&boundaryCondition);
1080 
1081  NeighborhoodIteratorType highIter(highPassOperator.GetRadius(), imgHighPass, imgLowPass->GetRequestedRegion());
1082  highIter.OverrideBoundaryCondition(&boundaryCondition);
1083 
1084  lowIter.GoToBegin();
1085  highIter.GoToBegin();
1086  out.GoToBegin();
1087 
1088  while (!out.IsAtEnd())
1089  {
1090  out.Set((innerProduct(lowIter, lowPassOperator)
1091  + innerProduct(highIter, highPassOperator)) / 2.);
1092 
1093  ++lowIter;
1094  ++highIter;
1095  ++out;
1096 
1097  reporter.CompletedPixel();
1098  }
1099  } // end for each imgLowPass/imghighPass pair of entry
1100  } // end multiscale case
1101 
1102  if (dir != InputImageDimension - 1)
1103  {
1104  // Note that outputImageRegion correspond to the actual region of (local) input !
1105  OutputImageRegionType outputImageRegion;
1106  this->CallCopyInputRegionToOutputRegion(dir, outputImageRegion, inputRegionForThread);
1107 
1108  ThreadedGenerateDataAtDimensionN(dir + 1, reporter, outputImageRegion, threadId);
1109  }
1110 
1111 }
1112 
1113 template <class TInputImage, class TOutputImage, class TWaveletOperator>
1114 void
1117  (unsigned int direction,
1118  itk::ProgressReporter& reporter,
1119  const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId)
1120 {
1121  OutputImageRegionType outputImageRegion;
1122  this->CallCopyInputRegionToOutputRegion(direction, outputImageRegion, outputRegionForThread);
1123 
1124  // Low pass part calculation
1125  LowPassOperatorType lowPassOperator;
1126  lowPassOperator.SetDirection(direction);
1127  lowPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
1128  lowPassOperator.CreateDirectional();
1129 
1130  // High pass part calculation
1131  HighPassOperatorType highPassOperator;
1132  highPassOperator.SetDirection(direction);
1133  highPassOperator.SetUpSampleFactor(this->GetUpSampleFilterFactor());
1134  highPassOperator.CreateDirectional();
1135 
1136  // typedef for the iterations over the input image
1137  typedef itk::ConstNeighborhoodIterator<OutputImageType> NeighborhoodIteratorType;
1138  typedef itk::NeighborhoodInnerProduct<OutputImageType> InnerProductType;
1139  // Faces iterations
1140  typename NeighborhoodIteratorType::RadiusType radiusMax;
1141  for (unsigned int i = 0; i < InputImageDimension; ++i)
1142  {
1143  radiusMax[i] = lowPassOperator.GetRadius(i);
1144  if (radiusMax[i] < highPassOperator.GetRadius(i)) radiusMax[i] = highPassOperator.GetRadius(i);
1145  }
1146 
1147  // The multiresolution case requires a SubsampleImageFilter step
1148  if (m_SubsampleImageFactor > 1)
1149  {
1150  for (unsigned int i = 0; i < m_InternalImages[direction - 1].size(); i += 2)
1151  {
1152  InputImagePointerType imgLowPass = m_InternalImages[direction - 1][i];
1153  InputImagePointerType imgHighPass = m_InternalImages[direction - 1][i + 1];
1154 
1155  OutputImagePointerType outputImage = this->GetOutput();
1156  if (direction < InputImageDimension - 1)
1157  {
1158  outputImage = m_InternalImages[direction][i / 2];
1159  }
1160 
1162  typename FilterType::InputImageIndexType delta;
1163  delta.Fill(1);
1164  delta[direction] = this->GetSubsampleImageFactor();
1165 
1166  InputImagePointerType cropedLowPass = InputImageType::New();
1167  cropedLowPass->SetRegions(outputRegionForThread);
1168  cropedLowPass->Allocate();
1169  cropedLowPass->FillBuffer(0.);
1170  itk::ImageRegionIterator<InputImageType> cropedLowPassIt(cropedLowPass, outputRegionForThread);
1171  itk::ImageRegionIterator<InputImageType> imgLowPassIt(imgLowPass, outputRegionForThread);
1172  for (cropedLowPassIt.GoToBegin(), imgLowPassIt.GoToBegin();
1173  !cropedLowPassIt.IsAtEnd() && !imgLowPassIt.IsAtEnd();
1174  ++cropedLowPassIt, ++imgLowPassIt)
1175  {
1176  cropedLowPassIt.Set(imgLowPassIt.Get());
1177  }
1178 
1179  typename FilterType::Pointer overSampledLowPass = FilterType::New();
1180  overSampledLowPass->SetInput(cropedLowPass);
1181  overSampledLowPass->SetSubsampleFactor(delta);
1182  overSampledLowPass->Update();
1183 
1184  InputImagePointerType cropedHighPass = InputImageType::New();
1185  cropedHighPass->SetRegions(outputRegionForThread);
1186  cropedHighPass->Allocate();
1187  cropedHighPass->FillBuffer(0.);
1188  itk::ImageRegionIterator<InputImageType> cropedHighPassIt(cropedHighPass, outputRegionForThread);
1189  itk::ImageRegionIterator<InputImageType> imgHighPassIt(imgHighPass, outputRegionForThread);
1190  for (cropedHighPassIt.GoToBegin(), imgHighPassIt.GoToBegin();
1191  !cropedHighPassIt.IsAtEnd() && !imgHighPassIt.IsAtEnd();
1192  ++cropedHighPassIt, ++imgHighPassIt)
1193  {
1194  cropedHighPassIt.Set(imgHighPassIt.Get());
1195  }
1196 
1197  typename FilterType::Pointer overSampledHighPass = FilterType::New();
1198  overSampledHighPass->SetInput(cropedHighPass);
1199  overSampledHighPass->SetSubsampleFactor(delta);
1200  overSampledHighPass->Update();
1201 
1202  InnerProductType innerProduct;
1203 
1205  overSampledLowPass->GetOutput()->GetRequestedRegion());
1206 
1207  // TODO: This might be the cause of the multithreading bug : we use a neighborhood iterator on cropped data
1208  // Are we sure that we have cropped enough data to access the neighborhood ?
1209  NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(),
1210  overSampledLowPass->GetOutput(),
1211  overSampledLowPass->GetOutput()->GetRequestedRegion());
1213  lowIter.OverrideBoundaryCondition(&boundaryCondition);
1214 
1215  NeighborhoodIteratorType highIter(highPassOperator.GetRadius(),
1216  overSampledHighPass->GetOutput(),
1217  overSampledHighPass->GetOutput()->GetRequestedRegion());
1218  highIter.OverrideBoundaryCondition(&boundaryCondition);
1219 
1220  lowIter.GoToBegin();
1221  highIter.GoToBegin();
1222  out.GoToBegin();
1223 
1224  while (!out.IsAtEnd())
1225  {
1226  out.Set(innerProduct(lowIter, lowPassOperator)
1227  + innerProduct(highIter, highPassOperator));
1228 
1229  ++lowIter;
1230  ++highIter;
1231  ++out;
1232 
1233  reporter.CompletedPixel();
1234  }
1235  } // end for each overSampledLowPass/overSampledhighPass pair of entry
1236  }
1237  else // multiscale case
1238  {
1239  for (unsigned int i = 0; i < m_InternalImages[direction - 1].size(); i += 2)
1240  {
1241  InputImagePointerType imgLowPass = m_InternalImages[direction - 1][i];
1242  InputImagePointerType imgHighPass = m_InternalImages[direction - 1][i + 1];
1243 
1244  OutputImagePointerType outputImage = this->GetOutput();
1245  if (direction < InputImageDimension - 1)
1246  {
1247  outputImage = m_InternalImages[direction][i / 2];
1248  }
1249 
1250  InnerProductType innerProduct;
1251 
1252  itk::ImageRegionIterator<OutputImageType> out(outputImage, imgLowPass->GetRequestedRegion());
1253 
1254  NeighborhoodIteratorType lowIter(lowPassOperator.GetRadius(), imgLowPass, imgLowPass->GetRequestedRegion());
1256  lowIter.OverrideBoundaryCondition(&boundaryCondition);
1257 
1258  NeighborhoodIteratorType highIter(highPassOperator.GetRadius(), imgHighPass, imgLowPass->GetRequestedRegion());
1259  highIter.OverrideBoundaryCondition(&boundaryCondition);
1260 
1261  lowIter.GoToBegin();
1262  highIter.GoToBegin();
1263  out.GoToBegin();
1264 
1265  while (!out.IsAtEnd())
1266  {
1267  out.Set((innerProduct(lowIter, lowPassOperator)
1268  + innerProduct(highIter, highPassOperator)) / 2.);
1269 
1270  ++lowIter;
1271  ++highIter;
1272  ++out;
1273 
1274  reporter.CompletedPixel();
1275  }
1276  } // end for each imgLowPass/imghighPass pair of entry
1277  }
1278 
1279  if (direction < InputImageDimension - 1)
1280  {
1281  ThreadedGenerateDataAtDimensionN(direction + 1, reporter, outputImageRegion, threadId);
1282  }
1283 }
1284 
1285 } // end of namespace
1286 
1287 #endif

Generated at Sat Aug 30 2014 16:33:18 for Orfeo Toolbox with doxygen 1.8.3.1