OTB  5.0.0
Orfeo Toolbox
otbMarkovRandomFieldFilter.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 
19 #ifndef __otbMarkovRandomFieldFilter_txx
20 #define __otbMarkovRandomFieldFilter_txx
22 
23 namespace otb
24 {
25 template<class TInputImage, class TClassifiedImage>
28  m_NumberOfClasses(0),
29  m_MaximumNumberOfIterations(50),
30  m_ErrorCounter(0),
31  m_ImageDeltaEnergy(0.0),
32  m_NeighborhoodRadius(1),
33  m_TotalNumberOfValidPixelsInOutputImage(1),
34  m_TotalNumberOfPixelsInInputImage(1),
35  m_ErrorTolerance(0.0),
36  m_SmoothingFactor(1.0),
37  m_NumberOfIterations(0),
38  m_Lambda(1.0),
39  m_ExternalClassificationSet(false),
40  m_StopCondition(MaximumNumberOfIterations)
41 {
42  m_Generator = RandomGeneratorType::GetInstance();
43  m_Generator->SetSeed();
44 
45  this->SetNumberOfRequiredInputs(1);
46  if ((int) InputImageDimension != (int) ClassifiedImageDimension)
47  {
48  std::ostringstream msg;
49  msg << "Input image dimension: " << InputImageDimension << " != output image dimension: " <<
50  ClassifiedImageDimension;
51  throw itk::ExceptionObject(__FILE__, __LINE__, msg.str().c_str(), ITK_LOCATION);
52  }
53  m_InputImageNeighborhoodRadius.Fill(m_NeighborhoodRadius);
54  // m_MRFNeighborhoodWeight.resize(0);
55  // m_NeighborInfluence.resize(0);
56  // m_DummyVector.resize(0);
57  // this->SetMRFNeighborhoodWeight( m_DummyVector );
58  // this->SetDefaultMRFNeighborhoodWeight();
59 
60  //srand((unsigned)time(0));
61 
62 }
63 
64 template<class TInputImage, class TClassifiedImage>
65 void
67 ::SetTrainingInput(const TrainingImageType * trainingImage)
68 {
69  m_ExternalClassificationSet = true;
70  // Process object is not const-correct so the const_cast is required here
71  this->itk::ProcessObject::SetNthInput(1, const_cast<TrainingImageType *>(trainingImage));
72  this->Modified();
73 }
74 
75 template <class TInputImage, class TClassifiedImage>
79 {
80  if (this->GetNumberOfInputs() < 2)
81  {
82  return 0;
83  }
84  return static_cast<const TrainingImageType *>
86 }
87 
88 template<class TInputImage, class TClassifiedImage>
89 void
91 ::PrintSelf(std::ostream& os, itk::Indent indent) const
92 {
93  Superclass::PrintSelf(os, indent);
94 
95  os << indent << " MRF Image filter object " << std::endl;
96 
97  os << indent << " Number of classes: " << m_NumberOfClasses << std::endl;
98 
99  os << indent << " Maximum number of iterations: " <<
100  m_MaximumNumberOfIterations << std::endl;
101 
102  os << indent << " Error tolerance for convergence: " <<
103  m_ErrorTolerance << std::endl;
104 
105  os << indent << " Size of the MRF neighborhood radius:" <<
106  m_InputImageNeighborhoodRadius << std::endl;
107 
108  os << indent << "StopCondition: "
109  << m_StopCondition << std::endl;
110 
111  os << indent << " Number of iterations: " <<
112  m_NumberOfIterations << std::endl;
113 
114  os << indent << " Lambda: " <<
115  m_Lambda << std::endl;
116 } // end PrintSelf
117 
121 template <class TInputImage, class TClassifiedImage>
122 void
125 {
126  // this filter requires the all of the input images
127  // to be at the size of the output requested region
128  InputImagePointer inputPtr =
129  const_cast<InputImageType *>(this->GetInput());
130  OutputImagePointer outputPtr = this->GetOutput();
131  inputPtr->SetRequestedRegion(outputPtr->GetRequestedRegion());
132 }
134 
138 template <class TInputImage, class TClassifiedImage>
139 void
142 {
143  // this filter requires the all of the output image to be in
144  // the buffer
145  TClassifiedImage *imgData;
146  imgData = dynamic_cast<TClassifiedImage*>(output);
147  imgData->SetRequestedRegionToLargestPossibleRegion();
148 }
150 
154 template <class TInputImage, class TClassifiedImage>
155 void
158 {
159  typename TInputImage::ConstPointer input = this->GetInput();
160  typename TClassifiedImage::Pointer output = this->GetOutput();
161  output->SetLargestPossibleRegion(input->GetLargestPossibleRegion());
162 }
164 
165 template<class TInputImage, class TClassifiedImage>
166 void
169 {
170 
171 // InputImageConstPointer inputImage = this->GetInput();
172 
173  //Allocate memory for the labelled images
174  this->Allocate();
175 
176  //Branch the pipeline
177  this->Initialize();
178 
179  //Run the Markov random field
180  this->ApplyMarkovRandomFieldFilter();
181 
182 } // end GenerateData
183 
187 template<class TInputImage, class TClassifiedImage>
188 void
190 ::SetNeighborhoodRadius(const unsigned long radiusValue)
191 {
192  //Set up the neighbor hood
193  NeighborhoodRadiusType radius;
194  for (unsigned int i = 0; i < InputImageDimension; ++i)
195  {
196  radius[i] = radiusValue;
197  }
198  this->SetNeighborhoodRadius(radius);
200 
201 } // end SetNeighborhoodRadius
202 
206 template<class TInputImage, class TClassifiedImage>
207 void
209 ::SetNeighborhoodRadius(const unsigned long *radiusArray)
210 {
211  NeighborhoodRadiusType radius;
212  for (unsigned int i = 0; i < InputImageDimension; ++i)
213  {
214  radius[i] = radiusArray[i];
215  }
216  //Set up the neighbor hood
217  this->SetNeighborhoodRadius(radius);
219 
220 } // end SetNeighborhoodRadius
221 
225 template<class TInputImage, class TClassifiedImage>
226 void
229 {
230  //Set up the neighbor hood
231  for (unsigned int i = 0; i < InputImageDimension; ++i)
232  {
233  m_InputImageNeighborhoodRadius[i] = radius[i];
234  m_LabelledImageNeighborhoodRadius[i] = radius[i];
235  }
237 
238 } // end SetNeighborhoodRadius
239 //-------------------------------------------------------
240 
244 template<class TInputImage, class TClassifiedImage>
245 void
248 {
249  if (m_NumberOfClasses <= 0)
250  {
251  throw itk::ExceptionObject(__FILE__, __LINE__, "NumberOfClasses <= 0.", ITK_LOCATION);
252  }
253 
254  //Set the output labelled and allocate the memory
255  LabelledImagePointer outputPtr = this->GetOutput();
256 
257  //Allocate the output buffer memory
258  outputPtr->SetBufferedRegion(outputPtr->GetRequestedRegion());
259  outputPtr->Allocate();
260 
261  //Copy input data in the output buffer memory or
262  //initialize to random values if not set
264  outImageIt(outputPtr, outputPtr->GetRequestedRegion());
265 
266  if (m_ExternalClassificationSet)
267  {
268  typename TrainingImageType::ConstPointer trainingImage = this->GetTrainingInput();
270  trainingImageIt(trainingImage, outputPtr->GetRequestedRegion());
271 
272  while (!outImageIt.IsAtEnd())
273  {
274  LabelledImagePixelType labelvalue = static_cast<LabelledImagePixelType> (trainingImageIt.Get());
275 
276  outImageIt.Set(labelvalue);
277  ++trainingImageIt;
278  ++outImageIt;
279  } // end while
280  }
281  else //set to random value
282  {
283 // srand((unsigned)time(0));
284  while (!outImageIt.IsAtEnd())
285  {
286  LabelledImagePixelType randomvalue = static_cast<LabelledImagePixelType>(
287  m_Generator->GetIntegerVariate() % static_cast<int>(m_NumberOfClasses)
288  );
289  outImageIt.Set(randomvalue);
290  ++outImageIt;
291  } // end while
292 
293  }
294 
295 } // Allocate
296 
300 template<class TInputImage, class TClassifiedImage>
301 void
303 ::Initialize() throw (itk::ExceptionObject)
304  {
305 
306  m_ImageDeltaEnergy = 0.0;
307 
308  InputImageSizeType inputImageSize =
309  this->GetInput()->GetBufferedRegion().GetSize();
310 
311  //---------------------------------------------------------------------
312  //Get the number of valid pixels in the output MRF image
313  //---------------------------------------------------------------------
314 
315  m_TotalNumberOfPixelsInInputImage = 1;
316  m_TotalNumberOfValidPixelsInOutputImage = 1;
317 
318  for (unsigned int i = 0; i < InputImageDimension; ++i)
319  {
320  m_TotalNumberOfPixelsInInputImage *= static_cast<int>(inputImageSize[i]);
321 
322  m_TotalNumberOfValidPixelsInOutputImage *=
323  (static_cast<int>(inputImageSize[i])
324  - 2 * m_InputImageNeighborhoodRadius[i]);
325  }
326 
327  srand((unsigned) time(0));
328 
329  if (!m_EnergyRegularization)
330  {
331  itkExceptionMacro(<< "EnergyRegularization is not present");
332  }
333 
334  if (!m_EnergyFidelity)
335  {
336  itkExceptionMacro(<< "EnergyFidelity is not present");
337  }
338 
339  if (!m_Optimizer)
340  {
341  itkExceptionMacro(<< "Optimizer is not present");
342  }
343 
344  if (!m_Sampler)
345  {
346  itkExceptionMacro(<< "Sampler is not present");
347  }
348 
349  m_Sampler->SetLambda(m_Lambda);
350  m_Sampler->SetEnergyRegularization(m_EnergyRegularization);
351  m_Sampler->SetEnergyFidelity(m_EnergyFidelity);
352  m_Sampler->SetNumberOfClasses(m_NumberOfClasses);
353  }
354 
358 template<class TInputImage, class TClassifiedImage>
359 void
362 {
363 
364  //Note: error should be defined according to the number of valid pixel in the output
365  int maxNumPixelError = itk::Math::Round<int, double> (m_ErrorTolerance *
366  m_TotalNumberOfPixelsInInputImage);
367 
368  m_NumberOfIterations = 0;
369  m_ErrorCounter = m_TotalNumberOfValidPixelsInOutputImage;
370 
371  while ((m_NumberOfIterations < m_MaximumNumberOfIterations) &&
372  (m_ErrorCounter >= maxNumPixelError))
373  {
374  otbMsgDevMacro(<< "Iteration No." << m_NumberOfIterations);
375 
376  this->MinimizeOnce();
377 
378  otbMsgDevMacro(<< "m_ErrorCounter/m_TotalNumberOfPixelsInInputImage: "
379  << m_ErrorCounter / ((double) (m_TotalNumberOfPixelsInInputImage)));
380  otbMsgDevMacro(<< "m_ImageDeltaEnergy: " << m_ImageDeltaEnergy);
381 
382  ++m_NumberOfIterations;
383 
384  }
385 
386  otbMsgDevMacro(<< "m_NumberOfIterations: " << m_NumberOfIterations);
387  otbMsgDevMacro(<< "m_MaximumNumberOfIterations: " << m_MaximumNumberOfIterations);
388  otbMsgDevMacro(<< "m_ErrorCounter: " << m_ErrorCounter);
389  otbMsgDevMacro(<< "maxNumPixelError: " << maxNumPixelError);
390 
391  //Determine stop condition
392  if (m_NumberOfIterations >= m_MaximumNumberOfIterations)
393  {
394  m_StopCondition = MaximumNumberOfIterations;
395  }
396  else if (m_ErrorCounter <= maxNumPixelError)
397  {
398  m_StopCondition = ErrorTolerance;
399  }
400 
401 } // ApplyMarkovRandomFieldFilter
402 
406 template<class TInputImage, class TClassifiedImage>
407 void
410 {
412  labelledIterator(m_LabelledImageNeighborhoodRadius, this->GetOutput(),
413  this->GetOutput()->GetLargestPossibleRegion());
415  dataIterator(m_InputImageNeighborhoodRadius, this->GetInput(),
416  this->GetInput()->GetLargestPossibleRegion());
417  m_ErrorCounter = 0;
419 
420  for (labelledIterator.GoToBegin(), dataIterator.GoToBegin();
421  !labelledIterator.IsAtEnd();
422  ++labelledIterator, ++dataIterator)
423  {
424 
426  bool changeValueBool;
427  m_Sampler->Compute(dataIterator, labelledIterator);
428  value = m_Sampler->GetValue();
429  changeValueBool = m_Optimizer->Compute(m_Sampler->GetDeltaEnergy());
430  if (changeValueBool)
431  {
432  labelledIterator.SetCenterPixel(value);
433  ++m_ErrorCounter;
434  m_ImageDeltaEnergy += m_Sampler->GetDeltaEnergy();
435  }
436  }
437 
438 }
439 
440 } // namespace otb
441 
442 #endif
virtual bool IsAtEnd() const
void SetNeighborhoodRadius(const NeighborhoodRadiusType &)
virtual void SetTrainingInput(const TrainingImageType *trainingImage)
TClassifiedImage::PixelType LabelledImagePixelType
OutputImageType::Pointer OutputImagePointer
virtual void SetCenterPixel(const PixelType &p)
TInputImage::SizeType NeighborhoodRadiusType
const TrainingImageType * GetTrainingInput(void)
DataObject * GetInput(const DataObjectIdentifierType &key)
TClassifiedImage::Pointer LabelledImagePointer
virtual void EnlargeOutputRequestedRegion(itk::DataObject *)
virtual void SetNthInput(DataObjectPointerArraySizeType num, DataObject *input)
void PrintSelf(std::ostream &os, itk::Indent indent) const
#define otbMsgDevMacro(x)
Definition: otbMacro.h:95