OTB  6.1.0
Orfeo Toolbox
otbBinaryFunctorNeighborhoodJoinHistogramImageFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2017 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 otbBinaryFunctorNeighborhoodJoinHistogramImageFilter_txx
22 #define otbBinaryFunctorNeighborhoodJoinHistogramImageFilter_txx
23 
25 #include "itkImageRegionIterator.h"
27 #include "itkProgressReporter.h"
28 
29 namespace otb
30 {
31 
35 template <class TInputImage1, class TInputImage2,
36  class TOutputImage, class TFunction>
39 {
40  this->SetNumberOfRequiredInputs(2);
41  m_Radius = 3;
42  m_UsePaddingValue = false;
43  m_UpperBoundIncreaseFactor = 0.001;
44 }
46 
50 template <class TInputImage1, class TInputImage2,
51  class TOutputImage, class TFunction>
52 void
54 ::SetInput1(const TInputImage1 * image1)
55 {
56  // Process object is not const-correct so the const casting is required.
57  this->SetNthInput(0, const_cast<TInputImage1 *>(image1));
58 }
59 
63 template <class TInputImage1, class TInputImage2,
64  class TOutputImage, class TFunction>
65 void
67 ::SetInput2(const TInputImage2 * image2)
68 {
69  // Process object is not const-correct so the const casting is required.
70  this->SetNthInput(1, const_cast<TInputImage2 *>(image2));
71 }
72 
73 template <class TInputImage1, class TInputImage2,
74  class TOutputImage, class TFunction>
75 const TInputImage1 *
78 {
79  if (this->GetNumberOfInputs() < 1)
80  {
81  return ITK_NULLPTR;
82  }
83  return static_cast<const TInputImage1 *>(this->itk::ProcessObject::GetInput(0));
84 }
85 
86 template <class TInputImage1, class TInputImage2,
87  class TOutputImage, class TFunction>
88 const TInputImage2 *
91 {
92  if (this->GetNumberOfInputs() < 2)
93  {
94  return ITK_NULLPTR;
95  }
96  return static_cast<const TInputImage2 *>(this->itk::ProcessObject::GetInput(1));
97 }
98 
99 template <class TInputImage1, class TInputImage2,
100  class TOutputImage, class TFunction>
101 void
104 {
105  // call the superclass' implementation of this method
106  Superclass::GenerateInputRequestedRegion();
107 
108  // get pointers to the input and output
109  Input1ImagePointer inputPtr1 =
110  const_cast<TInputImage1 *>(this->GetInput1());
111  Input2ImagePointer inputPtr2 =
112  const_cast<TInputImage2 *>(this->GetInput2());
113  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
114 
115  if (!inputPtr1 || !inputPtr2 || !outputPtr)
116  {
117  return;
118  }
119  // get a copy of the input requested region (should equal the output
120  // requested region)
121  typename TInputImage1::RegionType inputRequestedRegion1, inputRequestedRegion2;
122  inputRequestedRegion1 = inputPtr1->GetRequestedRegion();
123 
124  // pad the input requested region by the operator radius
125  inputRequestedRegion1.PadByRadius(m_Radius);
126  inputRequestedRegion2 = inputRequestedRegion1;
127 
128  // crop the input requested region at the input's largest possible region
129  if (inputRequestedRegion1.Crop(inputPtr1->GetLargestPossibleRegion()))
130  {
131  inputPtr1->SetRequestedRegion(inputRequestedRegion1);
132  }
133  else
134  {
135  // Couldn't crop the region (requested region is outside the largest
136  // possible region). Throw an exception.
137 
138  // store what we tried to request (prior to trying to crop)
139  inputPtr1->SetRequestedRegion(inputRequestedRegion1);
140 
141  // build an exception
142  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
143  std::ostringstream msg;
144  msg << this->GetNameOfClass()
145  << "::GenerateInputRequestedRegion()";
146  e.SetLocation(msg.str().c_str());
147  e.SetDescription("Requested region is (at least partially) outside the largest possible region of image 1.");
148  e.SetDataObject(inputPtr1);
149  throw e;
150  }
151  if (inputRequestedRegion2.Crop(inputPtr2->GetLargestPossibleRegion()))
152  {
153  inputPtr2->SetRequestedRegion(inputRequestedRegion2);
154  }
155  else
156  {
157  // Couldn't crop the region (requested region is outside the largest
158  // possible region). Throw an exception.
159 
160  // store what we tried to request (prior to trying to crop)
161  inputPtr2->SetRequestedRegion(inputRequestedRegion2);
162 
163  // build an exception
164  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
165  std::ostringstream msg;
166  msg << this->GetNameOfClass()
167  << "::GenerateInputRequestedRegion()";
168  e.SetLocation(msg.str().c_str());
169  e.SetDescription("Requested region is (at least partially) outside the largest possible region of image 1.");
170  e.SetDataObject(inputPtr2);
171  throw e;
172  }
173  return;
174 }
175 
176 template <class TInputImage1, class TInputImage2,
177  class TOutputImage, class TFunction>
178 void
181 {
182  this->ComputeHistogram();
183 }
184 
188 template <class TInputImage1, class TInputImage2,
189  class TOutputImage, class TFunction>
190 void
193 {
194  // Calculate min and max image values in input1 image.
195  Input1ImageConstPointer pInput1Image
196  = dynamic_cast<const TInputImage1*>(ProcessObjectType::GetInput(0));
197 
198  Input1ImagePixelType minInput1, maxInput1;
200  pInput1Image->GetBufferedRegion());
201  fiIt.GoToBegin();
202  minInput1 = fiIt.Value();
203  maxInput1 = fiIt.Value();
204  ++fiIt;
205  while (!fiIt.IsAtEnd())
206  {
207  Input1ImagePixelType value = fiIt.Value();
208 
209  if (value < minInput1)
210  {
211  minInput1 = value;
212  }
213  else if (value > maxInput1)
214  {
215  maxInput1 = value;
216  }
217 
218  ++fiIt;
219  }
220 
221  // Calculate min and max image values in input2 image.
222  Input2ImageConstPointer pInput2Image
223  = dynamic_cast<const TInputImage2*>(ProcessObjectType::GetInput(1));
224  Input2ImagePixelType minInput2, maxInput2;
226  pInput2Image->GetBufferedRegion());
227  miIt.GoToBegin();
228  minInput2 = miIt.Value();
229  maxInput2 = miIt.Value();
230  ++miIt;
231  while (!miIt.IsAtEnd())
232  {
233  Input2ImagePixelType value = miIt.Value();
234 
235  if (value < minInput2)
236  {
237  minInput2 = value;
238  }
239  else if (value > maxInput2)
240  {
241  maxInput2 = value;
242  }
243  ++miIt;
244  }
245 
246 
247  // Set the size of the upper and lower bounds of the histogram:
248  MeasurementVectorType lowerBound, upperBound;
249  lowerBound.SetSize(2);
250  upperBound.SetSize(2);
251 
252  // Initialize the upper and lower bounds of the histogram.
253  lowerBound[0] = minInput1;
254  lowerBound[1] = minInput2;
255  upperBound[0] =
256  maxInput1 + (maxInput1 - minInput1) * m_UpperBoundIncreaseFactor;
257  upperBound[1] =
258  maxInput2 + (maxInput2 - minInput2) * m_UpperBoundIncreaseFactor;
259 
262 
263  typename Input1ImageType::RegionType input1Region;
264  typename Input1ImageType::RegionType input2Region;
265 
266  input1Region = pInput1Image->GetRequestedRegion();
267  Input1IteratorType ti1(pInput1Image, input1Region);
268 
269  input2Region = pInput2Image->GetRequestedRegion();
270  Input2IteratorType ti2(pInput2Image, input2Region);
271 
272  typename HistogramType::Pointer histogram = HistogramType::New();
273  HistogramSizeType histogramSize(2);
274  histogramSize.Fill(256);
275  histogram->SetMeasurementVectorSize(2);
276  histogram->Initialize(histogramSize, lowerBound, upperBound);
277 
278  ti1.GoToBegin();
279  ti2.GoToBegin();
280  while (!ti1.IsAtEnd() && !ti2.IsAtEnd())
281  {
282 
283  typename HistogramType::IndexType sample(2);
284  sample[0] = ti1.Get();
285  sample[1] = ti2.Get();
287  sample[1] != itk::NumericTraits<Input2ImagePixelType>::Zero) histogram->IncreaseFrequencyOfIndex(sample, 1);
288 
289  ++ti1;
290  ++ti2;
291  }
292 
293  m_Histogram = histogram;
294 }
295 
299 template <class TInputImage1, class TInputImage2, class TOutputImage, class TFunction>
300 void
302 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
303  itk::ThreadIdType threadId)
304 {
307 
308  // We use dynamic_cast since inputs are stored as DataObjects. The
309  // ImageToJoinHistogramImageFilter::GetInput(int) always returns a pointer to a
310  // TInputImage1 so it cannot be used for the second input.
311  Input1ImageConstPointer inputPtr1
312  = dynamic_cast<const TInputImage1*>(ProcessObjectType::GetInput(0));
313  Input2ImageConstPointer inputPtr2
314  = dynamic_cast<const TInputImage2*>(ProcessObjectType::GetInput(1));
315  OutputImagePointer outputPtr = this->GetOutput(0);
316 
317  RadiusType1 r1;
318  r1.Fill(m_Radius);
319  NeighborhoodIteratorType1 neighInputIt1;
320 
321  RadiusType2 r2;
322  r2.Fill(m_Radius);
323  NeighborhoodIteratorType2 neighInputIt2;
324 
326 
327  // Find the data-set boundary "faces"
330  faceList1 = bC1(inputPtr1, outputRegionForThread, r1);
331 
334  faceList2 = bC2(inputPtr2, outputRegionForThread, r2);
335 
338 
339  // support progress methods/callbacks
340  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
341 
342  // Process each of the boundary faces. These are N-d regions which border
343  // the edge of the buffer.
344  for (fit1 = faceList1.begin(), fit2 = faceList2.begin();
345  fit1 != faceList1.end() && fit2 != faceList2.end();
346  ++fit1, ++fit2)
347  {
348  neighInputIt1 = itk::ConstNeighborhoodIterator<TInputImage1>(r1, inputPtr1, *fit1);
349  neighInputIt2 = itk::ConstNeighborhoodIterator<TInputImage2>(r1, inputPtr2, *fit2);
350  outputIt = itk::ImageRegionIterator<TOutputImage>(outputPtr, outputRegionForThread);
351 
352  outputIt = itk::ImageRegionIterator<TOutputImage>(outputPtr, *fit1);
353  neighInputIt1.OverrideBoundaryCondition(&nbc1);
354  neighInputIt1.GoToBegin();
355  neighInputIt2.OverrideBoundaryCondition(&nbc2);
356  neighInputIt2.GoToBegin();
357 
358  while (!outputIt.IsAtEnd())
359  {
360 
361  outputIt.Set(m_Functor(neighInputIt1, neighInputIt2, m_Histogram));
362 
363  ++neighInputIt1;
364  ++neighInputIt2;
365  ++outputIt;
366  progress.CompletedPixel();
367  }
368  }
369 }
370 
371 } // end namespace otb
372 
373 #endif
void SetDataObject(DataObject *dobj)
virtual void SetDescription(const std::string &s)
virtual void SetLocation(const std::string &s)
void Set(const PixelType &value) const
virtual void OverrideBoundaryCondition(const ImageBoundaryConditionPointerType i)
void Fill(TValue const &v)
OutputImageType::Pointer OutputImagePointer
unsigned int ThreadIdType
DataObject * GetInput(const DataObjectIdentifierType &key)
OutputImageType::RegionType OutputImageRegionType
bool IsAtEnd(void) const
void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId) ITK_OVERRIDE