Orfeo ToolBox  4.2
Orfeo ToolBox is not a black box
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 *>
85  (this->itk::ProcessObject::GetInput(1));
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 }
133 
137 template <class TInputImage, class TClassifiedImage>
138 void
141 {
142  // this filter requires the all of the output image to be in
143  // the buffer
144  TClassifiedImage *imgData;
145  imgData = dynamic_cast<TClassifiedImage*>(output);
146  imgData->SetRequestedRegionToLargestPossibleRegion();
147 }
148 
152 template <class TInputImage, class TClassifiedImage>
153 void
156 {
157  typename TInputImage::ConstPointer input = this->GetInput();
158  typename TClassifiedImage::Pointer output = this->GetOutput();
159  output->SetLargestPossibleRegion(input->GetLargestPossibleRegion());
160 }
161 
162 template<class TInputImage, class TClassifiedImage>
163 void
166 {
167 
168 // InputImageConstPointer inputImage = this->GetInput();
169 
170  //Allocate memory for the labelled images
171  this->Allocate();
172 
173  //Branch the pipeline
174  this->Initialize();
175 
176  //Run the Markov random field
177  this->ApplyMarkovRandomFieldFilter();
178 
179 } // end GenerateData
180 
184 template<class TInputImage, class TClassifiedImage>
185 void
187 ::SetNeighborhoodRadius(const unsigned long radiusValue)
188 {
189  //Set up the neighbor hood
190  NeighborhoodRadiusType radius;
191  for (unsigned int i = 0; i < InputImageDimension; ++i)
192  {
193  radius[i] = radiusValue;
194  }
195  this->SetNeighborhoodRadius(radius);
196 
197 } // end SetNeighborhoodRadius
198 
202 template<class TInputImage, class TClassifiedImage>
203 void
205 ::SetNeighborhoodRadius(const unsigned long *radiusArray)
206 {
207  NeighborhoodRadiusType radius;
208  for (unsigned int i = 0; i < InputImageDimension; ++i)
209  {
210  radius[i] = radiusArray[i];
211  }
212  //Set up the neighbor hood
213  this->SetNeighborhoodRadius(radius);
214 
215 } // end SetNeighborhoodRadius
216 
220 template<class TInputImage, class TClassifiedImage>
221 void
224 {
225  //Set up the neighbor hood
226  for (unsigned int i = 0; i < InputImageDimension; ++i)
227  {
228  m_InputImageNeighborhoodRadius[i] = radius[i];
229  m_LabelledImageNeighborhoodRadius[i] = radius[i];
230  }
231 
232 } // end SetNeighborhoodRadius
233 //-------------------------------------------------------
234 
238 template<class TInputImage, class TClassifiedImage>
239 void
242 {
243  if (m_NumberOfClasses <= 0)
244  {
245  throw itk::ExceptionObject(__FILE__, __LINE__, "NumberOfClasses <= 0.", ITK_LOCATION);
246  }
247 
248  //Set the output labelled and allocate the memory
249  LabelledImagePointer outputPtr = this->GetOutput();
250 
251  //Allocate the output buffer memory
252  outputPtr->SetBufferedRegion(outputPtr->GetRequestedRegion());
253  outputPtr->Allocate();
254 
255  //Copy input data in the output buffer memory or
256  //initialize to random values if not set
258  outImageIt(outputPtr, outputPtr->GetRequestedRegion());
259 
260  if (m_ExternalClassificationSet)
261  {
262  typename TrainingImageType::ConstPointer trainingImage = this->GetTrainingInput();
264  trainingImageIt(trainingImage, outputPtr->GetRequestedRegion());
265 
266  while (!outImageIt.IsAtEnd())
267  {
268  LabelledImagePixelType labelvalue = static_cast<LabelledImagePixelType> (trainingImageIt.Get());
269 
270  outImageIt.Set(labelvalue);
271  ++trainingImageIt;
272  ++outImageIt;
273  } // end while
274  }
275  else //set to random value
276  {
277 // srand((unsigned)time(0));
278  while (!outImageIt.IsAtEnd())
279  {
280  LabelledImagePixelType randomvalue = static_cast<LabelledImagePixelType>(
281  m_Generator->GetIntegerVariate() % static_cast<int>(m_NumberOfClasses)
282  );
283  outImageIt.Set(randomvalue);
284  ++outImageIt;
285  } // end while
286 
287  }
288 
289 } // Allocate
290 
294 template<class TInputImage, class TClassifiedImage>
295 void
297 ::Initialize() throw (itk::ExceptionObject)
298  {
299 
300  m_ImageDeltaEnergy = 0.0;
301 
302  InputImageSizeType inputImageSize =
303  this->GetInput()->GetBufferedRegion().GetSize();
304 
305  //---------------------------------------------------------------------
306  //Get the number of valid pixels in the output MRF image
307  //---------------------------------------------------------------------
308 
309  m_TotalNumberOfPixelsInInputImage = 1;
310  m_TotalNumberOfValidPixelsInOutputImage = 1;
311 
312  for (unsigned int i = 0; i < InputImageDimension; ++i)
313  {
314  m_TotalNumberOfPixelsInInputImage *= static_cast<int>(inputImageSize[i]);
315 
316  m_TotalNumberOfValidPixelsInOutputImage *=
317  (static_cast<int>(inputImageSize[i])
318  - 2 * m_InputImageNeighborhoodRadius[i]);
319  }
320 
321  srand((unsigned) time(0));
322 
323  if (!m_EnergyRegularization)
324  {
325  itkExceptionMacro(<< "EnergyRegularization is not present");
326  }
327 
328  if (!m_EnergyFidelity)
329  {
330  itkExceptionMacro(<< "EnergyFidelity is not present");
331  }
332 
333  if (!m_Optimizer)
334  {
335  itkExceptionMacro(<< "Optimizer is not present");
336  }
337 
338  if (!m_Sampler)
339  {
340  itkExceptionMacro(<< "Sampler is not present");
341  }
342 
343  m_Sampler->SetLambda(m_Lambda);
344  m_Sampler->SetEnergyRegularization(m_EnergyRegularization);
345  m_Sampler->SetEnergyFidelity(m_EnergyFidelity);
346  m_Sampler->SetNumberOfClasses(m_NumberOfClasses);
347  }
348 
352 template<class TInputImage, class TClassifiedImage>
353 void
356 {
357 
358  //Note: error should be defined according to the number of valid pixel in the output
359  int maxNumPixelError = itk::Math::Round<int, double> (m_ErrorTolerance *
360  m_TotalNumberOfPixelsInInputImage);
361 
362  m_NumberOfIterations = 0;
363  m_ErrorCounter = m_TotalNumberOfValidPixelsInOutputImage;
364 
365  while ((m_NumberOfIterations < m_MaximumNumberOfIterations) &&
366  (m_ErrorCounter >= maxNumPixelError))
367  {
368  otbMsgDevMacro(<< "Iteration No." << m_NumberOfIterations);
369 
370  this->MinimizeOnce();
371 
372  otbMsgDevMacro(<< "m_ErrorCounter/m_TotalNumberOfPixelsInInputImage: "
373  << m_ErrorCounter / ((double) (m_TotalNumberOfPixelsInInputImage)));
374  otbMsgDevMacro(<< "m_ImageDeltaEnergy: " << m_ImageDeltaEnergy);
375 
376  ++m_NumberOfIterations;
377 
378  }
379 
380  otbMsgDevMacro(<< "m_NumberOfIterations: " << m_NumberOfIterations);
381  otbMsgDevMacro(<< "m_MaximumNumberOfIterations: " << m_MaximumNumberOfIterations);
382  otbMsgDevMacro(<< "m_ErrorCounter: " << m_ErrorCounter);
383  otbMsgDevMacro(<< "maxNumPixelError: " << maxNumPixelError);
384 
385  //Determine stop condition
386  if (m_NumberOfIterations >= m_MaximumNumberOfIterations)
387  {
388  m_StopCondition = MaximumNumberOfIterations;
389  }
390  else if (m_ErrorCounter <= maxNumPixelError)
391  {
392  m_StopCondition = ErrorTolerance;
393  }
394 
395 } // ApplyMarkovRandomFieldFilter
396 
400 template<class TInputImage, class TClassifiedImage>
401 void
404 {
406  labelledIterator(m_LabelledImageNeighborhoodRadius, this->GetOutput(),
407  this->GetOutput()->GetLargestPossibleRegion());
409  dataIterator(m_InputImageNeighborhoodRadius, this->GetInput(),
410  this->GetInput()->GetLargestPossibleRegion());
411  m_ErrorCounter = 0;
412 
413  for (labelledIterator.GoToBegin(), dataIterator.GoToBegin();
414  !labelledIterator.IsAtEnd();
415  ++labelledIterator, ++dataIterator)
416  {
417 
419  bool changeValueBool;
420  m_Sampler->Compute(dataIterator, labelledIterator);
421  value = m_Sampler->GetValue();
422  changeValueBool = m_Optimizer->Compute(m_Sampler->GetDeltaEnergy());
423  if (changeValueBool)
424  {
425  labelledIterator.SetCenterPixel(value);
426  ++m_ErrorCounter;
427  m_ImageDeltaEnergy += m_Sampler->GetDeltaEnergy();
428  }
429  }
430 
431 }
432 
433 } // namespace otb
434 
435 #endif
virtual bool IsAtEnd() const
void SetNeighborhoodRadius(const NeighborhoodRadiusType &)
virtual void SetTrainingInput(const TrainingImageType *trainingImage)
Superclass::OutputImagePointer OutputImagePointer
TClassifiedImage::PixelType LabelledImagePixelType
TInputImage::SizeType NeighborhoodRadiusType
const TrainingImageType * GetTrainingInput(void)
TClassifiedImage::Pointer LabelledImagePointer
virtual void EnlargeOutputRequestedRegion(itk::DataObject *)
void PrintSelf(std::ostream &os, itk::Indent indent) const
#define otbMsgDevMacro(x)
Definition: otbMacro.h:94