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