Orfeo Toolbox  3.16
otbMeanShiftImageFilter.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 __otbMeanShiftImageFilter_txx
20 #define __otbMeanShiftImageFilter_txx
21 
22 #include "otbMeanShiftImageFilter.h"
23 
24 #include "itkImageRegionConstIteratorWithIndex.h"
25 #include "itkImageRegionIterator.h"
26 #include "itkImageRegionIteratorWithIndex.h"
27 #include "otbMacro.h"
28 
29 #include "msImageProcessor.h"
30 
31 namespace otb
32 {
33 template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
34 MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
35 ::MeanShiftImageFilter()
36 {
37  m_SpatialRadius = 3;
38  m_RangeRadius = 10;
39  m_MinimumRegionSize = 10;
40  m_Scale = 100000.;
41 
42  this->SetNumberOfOutputs(4);
43  this->SetNthOutput(1, OutputImageType::New());
44  this->SetNthOutput(2, LabeledOutputType::New());
45  this->SetNthOutput(3, LabeledOutputType::New());
46 }
47 
48 template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
49 const typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::OutputImageType *
50 MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
51 ::GetClusteredOutput() const
52 {
53  if (this->GetNumberOfOutputs() < 2)
54  {
55  return 0;
56  }
57  return static_cast<const OutputImageType *>(this->itk::ProcessObject::GetOutput(1));
58 }
59 
60 template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
61 typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::OutputImageType *
62 MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
63 ::GetClusteredOutput()
64 {
65  if (this->GetNumberOfOutputs() < 2)
66  {
67  return 0;
68  }
69  return static_cast<OutputImageType *>(this->itk::ProcessObject::GetOutput(1));
70 }
71 
72 template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
73 const typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::LabeledOutputType *
74 MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
75 ::GetLabeledClusteredOutput() const
76 {
77  if (this->GetNumberOfOutputs() < 3)
78  {
79  return 0;
80  }
81  return static_cast<const LabeledOutputType *>(this->itk::ProcessObject::GetOutput(2));
82 }
83 
84 template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
85 typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::LabeledOutputType *
86 MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
87 ::GetLabeledClusteredOutput()
88 {
89  if (this->GetNumberOfOutputs() < 3)
90  {
91  return 0;
92  }
93  return static_cast<LabeledOutputType *>(this->itk::ProcessObject::GetOutput(2));
94 }
95 
96 template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
97 const typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::LabeledOutputType *
98 MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
99 ::GetClusterBoundariesOutput() const
100 {
101  if (this->GetNumberOfOutputs() < 4)
102  {
103  return 0;
104  }
105  return static_cast<const LabeledOutputType *>(this->itk::ProcessObject::GetOutput(3));
106 }
107 
108 template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
109 typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::LabeledOutputType *
110 MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
111 ::GetClusterBoundariesOutput()
112 {
113  if (this->GetNumberOfOutputs() < 4)
114  {
115  return 0;
116  }
117  return static_cast<LabeledOutputType *>(this->itk::ProcessObject::GetOutput(3));
118 }
119 
120 template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
121 void
122 MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
123 ::AllocateOutputs()
124 {
125  typename OutputImageType::Pointer outputPtr = this->GetOutput();
126  typename OutputImageType::Pointer clusteredOutputPtr = this->GetClusteredOutput();
127  typename LabeledOutputType::Pointer labeledClusteredOutputPtr = this->GetLabeledClusteredOutput();
128  typename LabeledOutputType::Pointer clusterBoundariesOutputPtr = this->GetClusterBoundariesOutput();
129 
130  outputPtr->SetBufferedRegion(outputPtr->GetRequestedRegion());
131  outputPtr->Allocate();
132 
133  clusteredOutputPtr->SetBufferedRegion(clusteredOutputPtr->GetRequestedRegion());
134  clusteredOutputPtr->Allocate();
135 
136  labeledClusteredOutputPtr->SetBufferedRegion(labeledClusteredOutputPtr->GetRequestedRegion());
137  labeledClusteredOutputPtr->Allocate();
138 
139  clusterBoundariesOutputPtr->SetBufferedRegion(clusterBoundariesOutputPtr->GetRequestedRegion());
140  clusterBoundariesOutputPtr->Allocate();
141 }
142 
143 
144 template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
145 void
146 MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::EnlargeOutputRequestedRegion( itk::DataObject *output )
147 {
148 
149  // This filter requires all of the output images in the buffer.
150  for ( unsigned int j = 0; j < this->GetNumberOfOutputs(); j++ )
151  {
152  if ( this->itk::ProcessObject::GetOutput(j) )
153  {
154  this->itk::ProcessObject::GetOutput(j)->SetRequestedRegionToLargestPossibleRegion();
155  }
156  }
157 }
158 
159 template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
160 void
161 MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
162 ::GenerateData()
163 {
164  // Allocate the output image
165  this->AllocateOutputs();
166 
167  // Input and output pointers
168  typename InputImageType::ConstPointer inputPtr = this->GetInput();
169  typename OutputImageType::Pointer outputPtr = this->GetOutput();
170 
171 
172  double invScale = 1 / m_Scale;
173 
174  RegionType inputRequestedRegion =inputPtr->GetRequestedRegion();
175  RegionType outputRequestedRegion = inputPtr->GetRequestedRegion();
176 
177  inputRequestedRegion.PadByRadius(m_SpatialRadius);
178  if (!inputRequestedRegion.Crop(inputPtr->GetRequestedRegion()))
179  {
180  itkExceptionMacro(<< "error in requested region cropping");
181  }
182 
183  // Iterators
184  itk::ImageRegionConstIteratorWithIndex<InputImageType> inputIt(inputPtr, inputRequestedRegion);
185  itk::ImageRegionIterator<OutputImageType> outputIt(outputPtr, outputRequestedRegion);
186 
187  //create image processing object
188  msImageProcessor edisonProcessor;
189 
190  float * data = new float[inputRequestedRegion.GetNumberOfPixels() * inputPtr->GetNumberOfComponentsPerPixel()];
191 
192  unsigned int index = 0;
193 
194  for (inputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt)
195  {
196  TBufferConverter::PixelToFloatArray(data, index, inputIt.Get(), m_Scale);
197  index += inputPtr->GetNumberOfComponentsPerPixel();
198  }
199 
200  edisonProcessor.DefineLInput(data,
201  inputRequestedRegion.GetSize()[1],
202  inputRequestedRegion.GetSize()[0],
203  inputPtr->GetNumberOfComponentsPerPixel());
204 
205  //define default kernel paramerters...
206  kernelType k[2] = {Uniform, Uniform};
207  int P[2] = {2, static_cast<int>(inputPtr->GetNumberOfComponentsPerPixel())};
208  float tempH[2] = {1.0, 1.0};
209 
210  edisonProcessor.DefineKernel(k, tempH, P, 2);
211 
212  edisonProcessor.Filter(m_SpatialRadius, m_RangeRadius * m_Scale, MED_SPEEDUP);
213 
214  if (edisonProcessor.ErrorStatus)
215  {
216  itkExceptionMacro(<< "Error while running edison!");
217  }
218 
219  typename OutputImageType::Pointer tmpOutput = OutputImageType::New();
220  tmpOutput->SetRegions(inputRequestedRegion);
221  tmpOutput->SetNumberOfComponentsPerPixel(outputPtr->GetNumberOfComponentsPerPixel());
222  tmpOutput->Allocate();
223  itk::ImageRegionIterator<OutputImageType> tmpIt(tmpOutput, inputRequestedRegion);
224  itk::ImageRegionIterator<OutputImageType> tmp2It(tmpOutput, outputRequestedRegion);
225 
226  edisonProcessor.GetRawData(data);
227 
228  if (edisonProcessor.ErrorStatus)
229  {
230  itkExceptionMacro(<< "Error while running edison!");
231  }
232 
233  index = 0;
234  for (tmpIt.GoToBegin(); !tmpIt.IsAtEnd(); ++tmpIt)
235  {
236  OutputPixelType pixel;
237 
238  TBufferConverter::FloatArrayToPixel(data, index, pixel, outputPtr->GetNumberOfComponentsPerPixel(), invScale);
239  tmpIt.Set(pixel);
240  index += outputPtr->GetNumberOfComponentsPerPixel();
241  }
242 
243  tmp2It.GoToBegin();
244  outputIt.GoToBegin();
245 
246  while (!tmp2It.IsAtEnd() && !outputIt.IsAtEnd())
247  {
248  outputIt.Set(tmp2It.Get());
249  ++tmp2It;
250  ++outputIt;
251  }
252 
253 
254  // typename OutputImageType::Pointer outputPtr = this->GetOutput();
255  typename OutputImageType::Pointer clusteredOutputPtr = this->GetClusteredOutput();
256  typename LabeledOutputType::Pointer labeledClusteredOutputPtr = this->GetLabeledClusteredOutput();
257  typename LabeledOutputType::Pointer clusterBoudariesOutputPtr = this->GetClusterBoundariesOutput();
258  itk::ImageRegionIterator<OutputImageType> clusteredOutputIt(clusteredOutputPtr, outputRequestedRegion);
259 
260  index = 0;
261 
262  for (outputIt.GoToBegin(); !outputIt.IsAtEnd(); ++outputIt)
263  {
264  TBufferConverter::PixelToFloatArray(data, index, outputIt.Get(), m_Scale);
265  index += outputPtr->GetNumberOfComponentsPerPixel();
266  }
267 
268  edisonProcessor.DefineLInput(data,
269  outputRequestedRegion.GetSize()[1],
270  outputRequestedRegion.GetSize()[0],
271  outputPtr->GetNumberOfComponentsPerPixel());
272 
273  // define default kernel paramerters...
274  //kernelType k[2] = {Uniform, Uniform};
275  //int P[2] = {2, outputPtr->GetNumberOfComponentsPerPixel()};
276  //float tempH[2] = {1.0, 1.0};
277 
278  edisonProcessor.DefineKernel(k, tempH, P, 2);
279 
280  edisonProcessor.FuseRegions(m_RangeRadius * m_Scale, m_MinimumRegionSize);
281 
282  if (edisonProcessor.ErrorStatus)
283  {
284  itkExceptionMacro(<< "Error while running edison!");
285  }
286 
287  edisonProcessor.GetRawData(data);
288 
289  if (edisonProcessor.ErrorStatus)
290  {
291  itkExceptionMacro(<< "Error while running edison!");
292  }
293 
294  index = 0;
295  for (clusteredOutputIt.GoToBegin(); !clusteredOutputIt.IsAtEnd(); ++clusteredOutputIt)
296  {
297  OutputPixelType pixel;
298  TBufferConverter::FloatArrayToPixel(data, index, pixel,
299  clusteredOutputPtr->GetNumberOfComponentsPerPixel(), invScale);
300  clusteredOutputIt.Set(pixel);
301  index += clusteredOutputPtr->GetNumberOfComponentsPerPixel();
302  }
303 
304  int * labels = NULL;
305  float * modes = NULL;
306  int * modesPointsCount = NULL;
307 
308  edisonProcessor.GetRegions(&labels, &modes, &modesPointsCount);
309 
310  if (edisonProcessor.ErrorStatus)
311  {
312  itkExceptionMacro(<< "Error while running edison!");
313  }
314 
315  itk::ImageRegionIteratorWithIndex<LabeledOutputType> lcIt(labeledClusteredOutputPtr,
316  labeledClusteredOutputPtr->GetRequestedRegion());
317 
318  index = 0;
319 
320  labeledClusteredOutputPtr->FillBuffer(0);
321 
322  for (lcIt.GoToBegin(); !lcIt.IsAtEnd(); ++lcIt)
323  {
324  lcIt.Set(static_cast<LabelType>(labels[index]));
325  ++index;
326  }
327 
328  clusterBoudariesOutputPtr->FillBuffer(0);
329 
330  //define the boundaries
331  RegionList * regionList = edisonProcessor.GetBoundaries();
332  int * regionIndeces;
333  unsigned int numRegions = regionList->GetNumRegions();
334 
335  typename LabeledOutputType::IndexType boundIndex;
336 
337  for (LabelType label = 0; label < static_cast<LabelType>(numRegions); ++label)
338  {
339  OutputPixelType pixel;
340  TBufferConverter::FloatArrayToPixel(modes,
341  static_cast<unsigned int>(label *
342  clusteredOutputPtr->GetNumberOfComponentsPerPixel()),
343  pixel, clusteredOutputPtr->GetNumberOfComponentsPerPixel(), invScale);
344  // Filling the modes map
345  m_Modes[label] = pixel;
346  regionIndeces = regionList->GetRegionIndeces(static_cast<int>(label));
347  for (int i = 0; i < regionList->GetRegionCount(static_cast<int>(label)); ++i)
348  {
349  boundIndex[0] =
350  (regionIndeces[i] %
351  clusterBoudariesOutputPtr->GetRequestedRegion().GetSize()[0]) +
352  clusterBoudariesOutputPtr->GetRequestedRegion().GetIndex()[0];
353  boundIndex[1] =
354  (regionIndeces[i] /
355  clusterBoudariesOutputPtr->GetRequestedRegion().GetSize()[0]) +
356  clusterBoudariesOutputPtr->GetRequestedRegion().GetIndex()[1];
357  if (clusterBoudariesOutputPtr->GetBufferedRegion().IsInside(boundIndex))
358  {
359  clusterBoudariesOutputPtr->SetPixel(boundIndex, 1);
360  }
361  }
362  }
363 
364  // Free memory
365  delete[] labels;
366  delete[] modes;
367  delete[] modesPointsCount;
368  delete[] data;
369 }
370 
371 
372 template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
373 void
374 MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
375 ::PrintSelf(std::ostream& os, itk::Indent indent) const
376 {
377  Superclass::PrintSelf(os, indent);
378  os << indent << "Spatial radius: " << m_SpatialRadius << std::endl;
379  os << indent << "Range radius: " << m_RangeRadius << std::endl;
380  os << indent << "Minimum region size: " << m_MinimumRegionSize << std::endl;
381  os << indent << "Scale: " << m_Scale << std::endl;
382 }
383 
384 } // end namespace otb
385 
386 #endif
Generated at Sun May 19 2013 00:38:09 for Orfeo Toolbox with doxygen 1.8.3.1