OTB  6.7.0
Orfeo Toolbox
otbOverlapSaveConvolutionImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbOverlapSaveConvolutionImageFilter_hxx
22 #define otbOverlapSaveConvolutionImageFilter_hxx
23 
24 #include "itkConfigure.h"
25 
27 
31 #include "itkOffset.h"
32 #include "itkProgressReporter.h"
33 
34 //debug
35 #include "itkImageRegionIterator.h"
36 #include "otbMath.h"
37 
38 #ifdef ITK_USE_FFTWD
39 #include "itkFFTWCommon.h"
40 #endif
41 
42 namespace otb
43 {
44 
45 template <class TInputImage, class TOutputImage, class TBoundaryCondition>
48 {
49  m_Radius.Fill(1);
50  m_Filter.SetSize(3 * 3);
51  m_Filter.Fill(1);
52  m_NormalizeFilter = false;
53 }
54 
55 template <class TInputImage, class TOutputImage, class TBoundaryCondition>
56 void
58 ::GenerateInputRequestedRegion() throw (itk::InvalidRequestedRegionError)
59  {
60 #if defined ITK_USE_FFTWD
61  // call the superclass' implementation of this method
62  Superclass::GenerateInputRequestedRegion();
63 
64  // get pointers to the input and output
65  typename Superclass::InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
66  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
67 
68  if (!inputPtr || !outputPtr)
69  {
70  return;
71  }
72 
73  // get a copy of the input requested region (should equal the output
74  // requested region)
75  typename TInputImage::RegionType inputRequestedRegion;
76  inputRequestedRegion = inputPtr->GetRequestedRegion();
77 
78  // Pad by filter radius
79  inputRequestedRegion.PadByRadius(m_Radius);
80 
81  // crop the input requested region at the input's largest possible region
82  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
83  {
84  inputPtr->SetRequestedRegion(inputRequestedRegion);
85  return;
86  }
87  else
88  {
89  // Couldn't crop the region (requested region is outside the largest
90  // possible region). Throw an exception.
91 
92  // store what we tried to request (prior to trying to crop)
93  inputPtr->SetRequestedRegion(inputRequestedRegion);
94 
95  // build an exception
96  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
97  e.SetLocation(ITK_LOCATION);
98  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
99  e.SetDataObject(inputPtr);
100  throw e;
101  }
102 #else
103  itkGenericExceptionMacro(
104  <<
105  "The OverlapSaveConvolutionImageFilter can not operate without the FFTW library (double implementation). Please build ITK with USE_FFTD set to ON, and rebuild OTB.");
106 #endif
107  }
108 
109 template<class TInputImage, class TOutputImage, class TBoundaryCondition>
110 void
112 /* TODO commented out since multi-threading is not supported for the moment
113  * ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, itk::ThreadIdType threadId) */
114 ::GenerateData()
115 {
116 #if defined ITK_USE_FFTWD
117  // Input/Output pointers
118  typename OutputImageType::Pointer output = this->GetOutput();
119  typename InputImageType::ConstPointer input = this->GetInput();
120 
121  /* TODO: This is a patch to switch from GenerateData() to ThreadedGenerateData(). Remove these two lines
122  once multi-threading problem is solved */
123  this->AllocateOutputs();
124  OutputImageRegionType outputRegionForThread = output->GetRequestedRegion();
125 
126  // Size of the filter
127  typename InputImageType::SizeType sizeOfFilter;
128  sizeOfFilter[0] = 2 * m_Radius[0] + 1;
129  sizeOfFilter[1] = 2 * m_Radius[1] + 1;
130 
131  // Filter normalization
132  InputRealType norm;
133 
134  // Compute the input region for the given thread
135  OutputImageRegionType inputRegionForThread = outputRegionForThread;
136  inputRegionForThread.PadByRadius(m_Radius);
137 
138  // Compute the piece region for the input thread. Piece region is different
139  // from input region on boundaries
140  typename InputImageType::RegionType pieceRegion = inputRegionForThread;
141  typename InputImageType::SizeType pieceSize = pieceRegion.GetSize();
142  typename InputImageType::IndexType pieceIndex = pieceRegion.GetIndex();
143 
144  // Compute the size of the FFT and the size of the piece
145  unsigned int pieceNbOfPixel = pieceRegion.GetNumberOfPixels();
146  unsigned int sizeFFT = (pieceSize[0] / 2 + 1) * pieceSize[1];
147 
148  // Achieve the computation of the inputRegionForThread
149  inputRegionForThread.Crop(input->GetLargestPossibleRegion());
150  typename InputImageType::IndexType inputIndex = inputRegionForThread.GetIndex();
151  typename InputImageType::SizeType inputSize = inputRegionForThread.GetSize();
152 
153  //Iterator of input image
155  inputIt(m_Radius, input, inputRegionForThread);
156  inputIt.GoToBegin();
157 
158  //Iterator of output image
160  outputIt = itk::ImageRegionIteratorWithIndex<OutputImageType>(output, outputRegionForThread);
161 
162  //variables for loops
163  unsigned int i, j, k, l;
164 
165  // ITK proxy to the fftw library
166  typedef typename itk::fftw::Proxy<double> FFTWProxyType;
167 
168  //memory allocation
169  InputPixelType* resampledFilterPiece;
170  resampledFilterPiece = static_cast<FFTWProxyType::PixelType*>(fftw_malloc(pieceNbOfPixel * sizeof(InputPixelType)));
171 
172  FFTWProxyType::ComplexType* filterPieceFFT;
173  filterPieceFFT = static_cast<FFTWProxyType::ComplexType*>(fftw_malloc(sizeFFT * sizeof(FFTWProxyType::ComplexType)));
174 
175  InputPixelType *inputPiece;
176  inputPiece = static_cast<FFTWProxyType::PixelType*>(fftw_malloc(pieceNbOfPixel * sizeof(InputPixelType)));
177 
178  FFTWProxyType::ComplexType* inputPieceFFT;
179  inputPieceFFT = static_cast<FFTWProxyType::ComplexType*>(fftw_malloc(sizeFFT * sizeof(FFTWProxyType::ComplexType)));
180 
181  // Image piece FFT
182  FFTWProxyType::PlanType inputPlan = FFTWProxyType::Plan_dft_r2c_2d(pieceSize[1],
183  pieceSize[0],
184  inputPiece,
185  inputPieceFFT,
186  FFTW_MEASURE);
187 
188  // left zero padding
189  unsigned int leftskip = static_cast<unsigned int>(std::max(
190  (typename InputImageType::IndexValueType) 0, inputIndex[0] - pieceIndex[0]));
191  unsigned int topskip = pieceSize[0] * static_cast<unsigned int>(std::max(
192  (typename InputImageType::IndexValueType) 0, inputIndex[1] - pieceIndex[1]));
193 
194  // zero filling
195  memset(inputPiece,0,pieceNbOfPixel * sizeof(InputPixelType));
196 
197  // Filling the buffer with image values
198  for (l = 0; l < inputSize[1]; ++l)
199  {
200  for (k = 0; k < inputSize[0]; ++k)
201  {
202  inputPiece[topskip + pieceSize[0] * l + k + leftskip] = inputIt.GetCenterPixel();
203  ++inputIt;
204  }
205  }
206 
207  FFTWProxyType::Execute(inputPlan);
208 
209  // Resampled filter FFT
210  FFTWProxyType::PlanType filterPlan = FFTWProxyType::Plan_dft_r2c_2d(pieceSize[1],
211  pieceSize[0],
212  resampledFilterPiece,
213  filterPieceFFT,
214  FFTW_MEASURE);
215 
216  // zero filling
217  memset(resampledFilterPiece,0,pieceNbOfPixel * sizeof(InputPixelType));
218 
219  k = 0;
220  // Filling the buffer with filter values
221  for (j = 0; j < sizeOfFilter[1]; ++j)
222  {
223  for (i = 0; i < sizeOfFilter[0]; ++i)
224  {
225  resampledFilterPiece[i + j * pieceSize[0]] = m_Filter.GetElement(k); //Copy values
226  ++k;
227  }
228  }
229 
230  FFTWProxyType::Execute(filterPlan);
231 
232  // memory allocation for inverse FFT
233  FFTWProxyType::ComplexType* multipliedFFTarray;
234  multipliedFFTarray = static_cast<FFTWProxyType::ComplexType*>(fftw_malloc(sizeFFT * sizeof(FFTWProxyType::ComplexType)));
235 
236  FFTWProxyType::PixelType* inverseFFTpiece;
237  inverseFFTpiece = static_cast<FFTWProxyType::PixelType*>(fftw_malloc(pieceNbOfPixel * sizeof(FFTWProxyType::PixelType)));
238 
239  // Inverse FFT of the product of FFT (actually do filtering here)
240  FFTWProxyType::PlanType outputPlan = FFTWProxyType::Plan_dft_c2r_2d(pieceSize[1],
241  pieceSize[0],
242  multipliedFFTarray,
243  inverseFFTpiece,
244  FFTW_MEASURE);
245 
246  // Filling the buffer with complex product values
247  for (k = 0; k < sizeFFT; ++k)
248  {
249  //complex mutiplication
250  multipliedFFTarray[k][0] = inputPieceFFT[k][0] * filterPieceFFT[k][0] - inputPieceFFT[k][1] * filterPieceFFT[k][1];
251  multipliedFFTarray[k][1] = inputPieceFFT[k][0] * filterPieceFFT[k][1] + inputPieceFFT[k][1] * filterPieceFFT[k][0];
252  }
253 
254  FFTWProxyType::Execute(outputPlan);
255 
256  // Computing the filter normalization
257 
258  if (m_NormalizeFilter)
259  {
261  for (i = 0; i < sizeOfFilter[0] * sizeOfFilter[1]; ++i)
262  {
263  norm += static_cast<InputRealType>(m_Filter(i));
264  }
265  if (norm == 0.0)
266  {
267  norm = 1.0;
268  }
269  else
270  {
271  norm = 1 / norm;
272  }
273  }
274  else
275  {
276  norm = 1.0;
277  }
278 
279  // Fill the output image
280  outputIt.GoToBegin();
281  while (!outputIt.IsAtEnd())
282  {
283  typename InputImageType::IndexType index = outputIt.GetIndex();
284  unsigned int linearIndex =
285  (index[1] + sizeOfFilter[1] - 1 -
286  outputRegionForThread.GetIndex()[1]) * pieceSize[0] - 1 + index[0] + sizeOfFilter[0] -
287  outputRegionForThread.GetIndex()[0];
288  outputIt.Set(static_cast<OutputPixelType>((inverseFFTpiece[linearIndex] /
289  pieceNbOfPixel) * static_cast<double>(norm)));
290  ++outputIt;
291  }
292 
293  // destroy the FFT plans
294  FFTWProxyType::DestroyPlan(inputPlan);
295  FFTWProxyType::DestroyPlan(filterPlan);
296  FFTWProxyType::DestroyPlan(outputPlan);
297 
298  //frees memory
299  fftw_free(resampledFilterPiece);
300  fftw_free(inputPiece);
301  fftw_free(filterPieceFFT);
302  fftw_free(inputPieceFFT);
303  fftw_free(multipliedFFTarray);
304  fftw_free(inverseFFTpiece);
305 #else
306  itkGenericExceptionMacro(
307  <<
308  "The OverlapSaveConvolutionImageFilter can not operate without the FFTW library (double implementation). Please build ITK with USE_FFTD set to ON, and rebuild OTB.");
309 #endif
310 }
311 
313 template <class TInputImage, class TOutput, class TBoundaryCondition>
314 void
316 ::PrintSelf(std::ostream& os, itk::Indent indent) const
317 {
318  Superclass::PrintSelf(os, indent);
319  os << indent << "Radius: " << m_Radius << std::endl;
320  os << indent << "Normalize filter: " << m_NormalizeFilter << std::endl;
321 }
322 } // end namespace otb
324 
325 #endif
itk::Size< Monteverdi_DIMENSION > SizeType
Definition: mvdTypes.h:137
void SetDataObject(DataObject *dobj)
virtual void SetDescription(const std::string &s)
virtual void SetLocation(const std::string &s)
void Set(const PixelType &value) const
const IndexType & GetIndex() const
itk::Index< Monteverdi_DIMENSION > IndexType
Definition: mvdTypes.h:133
PixelType GetCenterPixel() const
itk::NumericTraits< InputPixelType >::RealType InputRealType
void PrintSelf(std::ostream &os, itk::Indent indent) const override
OutputImageType::RegionType OutputImageRegionType
const SizeValueType * GetSize() const