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