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