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