Orfeo ToolBox  4.2
Orfeo ToolBox is not a black box
otbBinaryFunctorNeighborhoodJoinHistogramImageFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef __otbBinaryFunctorNeighborhoodJoinHistogramImageFilter_txx
19 #define __otbBinaryFunctorNeighborhoodJoinHistogramImageFilter_txx
20 
22 #include "itkImageRegionIterator.h"
25 #include "itkProgressReporter.h"
26 
27 namespace otb
28 {
29 
33 template <class TInputImage1, class TInputImage2,
34  class TOutputImage, class TFunction>
37 {
38  this->SetNumberOfRequiredInputs(2);
39  m_Radius = 3;
40  m_UsePaddingValue = false;
41  m_UpperBoundIncreaseFactor = 0.001;
42 }
43 
47 template <class TInputImage1, class TInputImage2,
48  class TOutputImage, class TFunction>
49 void
51 ::SetInput1(const TInputImage1 * image1)
52 {
53  // Process object is not const-correct so the const casting is required.
54  this->SetNthInput(0, const_cast<TInputImage1 *>(image1));
55 }
56 
60 template <class TInputImage1, class TInputImage2,
61  class TOutputImage, class TFunction>
62 void
64 ::SetInput2(const TInputImage2 * image2)
65 {
66  // Process object is not const-correct so the const casting is required.
67  this->SetNthInput(1, const_cast<TInputImage2 *>(image2));
68 }
69 
70 template <class TInputImage1, class TInputImage2,
71  class TOutputImage, class TFunction>
72 const TInputImage1 *
75 {
76  if (this->GetNumberOfInputs() < 1)
77  {
78  return 0;
79  }
80  return static_cast<const TInputImage1 *>(this->itk::ProcessObject::GetInput(0));
81 }
82 
83 template <class TInputImage1, class TInputImage2,
84  class TOutputImage, class TFunction>
85 const TInputImage2 *
88 {
89  if (this->GetNumberOfInputs() < 2)
90  {
91  return 0;
92  }
93  return static_cast<const TInputImage2 *>(this->itk::ProcessObject::GetInput(1));
94 }
95 
96 template <class TInputImage1, class TInputImage2,
97  class TOutputImage, class TFunction>
98 void
101 {
102  // call the superclass' implementation of this method
103  Superclass::GenerateInputRequestedRegion();
104 
105  // get pointers to the input and output
106  Input1ImagePointer inputPtr1 =
107  const_cast<TInputImage1 *>(this->GetInput1());
108  Input2ImagePointer inputPtr2 =
109  const_cast<TInputImage2 *>(this->GetInput2());
110  typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
111 
112  if (!inputPtr1 || !inputPtr2 || !outputPtr)
113  {
114  return;
115  }
116  // get a copy of the input requested region (should equal the output
117  // requested region)
118  typename TInputImage1::RegionType inputRequestedRegion1, inputRequestedRegion2;
119  inputRequestedRegion1 = inputPtr1->GetRequestedRegion();
120 
121  // pad the input requested region by the operator radius
122  inputRequestedRegion1.PadByRadius(m_Radius);
123  inputRequestedRegion2 = inputRequestedRegion1;
124 
125  // crop the input requested region at the input's largest possible region
126  if (inputRequestedRegion1.Crop(inputPtr1->GetLargestPossibleRegion()))
127  {
128  inputPtr1->SetRequestedRegion(inputRequestedRegion1);
129  }
130  else
131  {
132  // Couldn't crop the region (requested region is outside the largest
133  // possible region). Throw an exception.
134 
135  // store what we tried to request (prior to trying to crop)
136  inputPtr1->SetRequestedRegion(inputRequestedRegion1);
137 
138  // build an exception
139  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
140  std::ostringstream msg;
141  msg << this->GetNameOfClass()
142  << "::GenerateInputRequestedRegion()";
143  e.SetLocation(msg.str().c_str());
144  e.SetDescription("Requested region is (at least partially) outside the largest possible region of image 1.");
145  e.SetDataObject(inputPtr1);
146  throw e;
147  }
148  if (inputRequestedRegion2.Crop(inputPtr2->GetLargestPossibleRegion()))
149  {
150  inputPtr2->SetRequestedRegion(inputRequestedRegion2);
151  }
152  else
153  {
154  // Couldn't crop the region (requested region is outside the largest
155  // possible region). Throw an exception.
156 
157  // store what we tried to request (prior to trying to crop)
158  inputPtr2->SetRequestedRegion(inputRequestedRegion2);
159 
160  // build an exception
161  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
162  std::ostringstream msg;
163  msg << this->GetNameOfClass()
164  << "::GenerateInputRequestedRegion()";
165  e.SetLocation(msg.str().c_str());
166  e.SetDescription("Requested region is (at least partially) outside the largest possible region of image 1.");
167  e.SetDataObject(inputPtr2);
168  throw e;
169  }
170  return;
171 }
172 
173 template <class TInputImage1, class TInputImage2,
174  class TOutputImage, class TFunction>
175 void
178 {
179  this->ComputeHistogram();
180 }
181 
185 template <class TInputImage1, class TInputImage2,
186  class TOutputImage, class TFunction>
187 void
190 {
191  // Calculate min and max image values in input1 image.
192  Input1ImageConstPointer pInput1Image
193  = dynamic_cast<const TInputImage1*>(ProcessObjectType::GetInput(0));
194 
195  Input1ImagePixelType minInput1, maxInput1;
197  pInput1Image->GetBufferedRegion());
198  fiIt.GoToBegin();
199  minInput1 = fiIt.Value();
200  maxInput1 = fiIt.Value();
201  ++fiIt;
202  while (!fiIt.IsAtEnd())
203  {
204  Input1ImagePixelType value = fiIt.Value();
205 
206  if (value < minInput1)
207  {
208  minInput1 = value;
209  }
210  else if (value > maxInput1)
211  {
212  maxInput1 = value;
213  }
214 
215  ++fiIt;
216  }
217 
218  // Calculate min and max image values in input2 image.
219  Input2ImageConstPointer pInput2Image
220  = dynamic_cast<const TInputImage2*>(ProcessObjectType::GetInput(1));
221  Input2ImagePixelType minInput2, maxInput2;
223  pInput2Image->GetBufferedRegion());
224  miIt.GoToBegin();
225  minInput2 = miIt.Value();
226  maxInput2 = miIt.Value();
227  ++miIt;
228  while (!miIt.IsAtEnd())
229  {
230  Input2ImagePixelType value = miIt.Value();
231 
232  if (value < minInput2)
233  {
234  minInput2 = value;
235  }
236  else if (value > maxInput2)
237  {
238  maxInput2 = value;
239  }
240  ++miIt;
241  }
242 
243 
244  // Set the size of the upper and lower bounds of the histogram:
245  MeasurementVectorType lowerBound, upperBound;
246  lowerBound.SetSize(2);
247  upperBound.SetSize(2);
248 
249  // Initialize the upper and lower bounds of the histogram.
250  lowerBound[0] = minInput1;
251  lowerBound[1] = minInput2;
252  upperBound[0] =
253  maxInput1 + (maxInput1 - minInput1) * m_UpperBoundIncreaseFactor;
254  upperBound[1] =
255  maxInput2 + (maxInput2 - minInput2) * m_UpperBoundIncreaseFactor;
256 
259 
260  typename Input1ImageType::RegionType input1Region;
261  typename Input1ImageType::RegionType input2Region;
262 
263  input1Region = pInput1Image->GetRequestedRegion();
264  Input1IteratorType ti1(pInput1Image, input1Region);
265 
266  input2Region = pInput2Image->GetRequestedRegion();
267  Input2IteratorType ti2(pInput2Image, input2Region);
268 
269  typename HistogramType::Pointer histogram = HistogramType::New();
270  HistogramSizeType histogramSize(2);
271  histogramSize.Fill(256);
272  histogram->SetMeasurementVectorSize(2);
273  histogram->Initialize(histogramSize, lowerBound, upperBound);
274 
275  ti1.GoToBegin();
276  ti2.GoToBegin();
277  while (!ti1.IsAtEnd() && !ti2.IsAtEnd())
278  {
279 
280  typename HistogramType::IndexType sample(2);
281  sample[0] = ti1.Get();
282  sample[1] = ti2.Get();
284  sample[1] != itk::NumericTraits<Input2ImagePixelType>::Zero) histogram->IncreaseFrequencyOfIndex(sample, 1);
285 
286  ++ti1;
287  ++ti2;
288  }
289 
290  m_Histogram = histogram;
291 }
292 
296 template <class TInputImage1, class TInputImage2, class TOutputImage, class TFunction>
297 void
299 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
300  itk::ThreadIdType threadId)
301 {
304 
305  // We use dynamic_cast since inputs are stored as DataObjects. The
306  // ImageToJoinHistogramImageFilter::GetInput(int) always returns a pointer to a
307  // TInputImage1 so it cannot be used for the second input.
308  Input1ImageConstPointer inputPtr1
309  = dynamic_cast<const TInputImage1*>(ProcessObjectType::GetInput(0));
310  Input2ImageConstPointer inputPtr2
311  = dynamic_cast<const TInputImage2*>(ProcessObjectType::GetInput(1));
312  OutputImagePointer outputPtr = this->GetOutput(0);
313 
314  RadiusType1 r1;
315  r1.Fill(m_Radius);
316  NeighborhoodIteratorType1 neighInputIt1;
317 
318  RadiusType2 r2;
319  r2.Fill(m_Radius);
320  NeighborhoodIteratorType2 neighInputIt2;
321 
323 
324  // Find the data-set boundary "faces"
327  faceList1 = bC1(inputPtr1, outputRegionForThread, r1);
328 
331  faceList2 = bC2(inputPtr2, outputRegionForThread, r2);
332 
335 
336  // support progress methods/callbacks
337  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
338 
339  // Process each of the boundary faces. These are N-d regions which border
340  // the edge of the buffer.
341  for (fit1 = faceList1.begin(), fit2 = faceList2.begin();
342  fit1 != faceList1.end() && fit2 != faceList2.end();
343  ++fit1, ++fit2)
344  {
345  neighInputIt1 = itk::ConstNeighborhoodIterator<TInputImage1>(r1, inputPtr1, *fit1);
346  neighInputIt2 = itk::ConstNeighborhoodIterator<TInputImage2>(r1, inputPtr2, *fit2);
347  outputIt = itk::ImageRegionIterator<TOutputImage>(outputPtr, outputRegionForThread);
348 
349  outputIt = itk::ImageRegionIterator<TOutputImage>(outputPtr, *fit1);
350  neighInputIt1.OverrideBoundaryCondition(&nbc1);
351  neighInputIt1.GoToBegin();
352  neighInputIt2.OverrideBoundaryCondition(&nbc2);
353  neighInputIt2.GoToBegin();
354 
355  while (!outputIt.IsAtEnd())
356  {
357 
358  outputIt.Set(m_Functor(neighInputIt1, neighInputIt2, m_Histogram));
359 
360  ++neighInputIt1;
361  ++neighInputIt2;
362  ++outputIt;
363  progress.CompletedPixel();
364  }
365  }
366 }
367 
368 } // end namespace otb
369 
370 #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)
virtual void ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, itk::ThreadIdType threadId)
void Fill(TValue const &v)
bool IsAtEnd(void) const
unsigned int ThreadIdType