OTB  9.0.0
Orfeo Toolbox
otbMeanShiftSmoothingImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2022 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 otbMeanShiftSmoothingImageFilter_hxx
23 #define otbMeanShiftSmoothingImageFilter_hxx
24 
25 #include "otbMeanShiftSmoothingImageFilter.h"
26 #include "itkImageRegionIterator.h"
27 #include "otbUnaryFunctorWithIndexWithOutputSizeImageFilter.h"
28 #include "otbMacro.h"
29 
30 #include "itkProgressReporter.h"
31 
32 
33 namespace otb
34 {
35 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
36 MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::MeanShiftSmoothingImageFilter()
37  : m_RangeBandwidth(16.),
38  m_RangeBandwidthRamp(0),
39  m_SpatialBandwidth(3)
40  // , m_SpatialRadius(???)
41  ,
42  m_Threshold(1e-3),
43  m_MaxIterationNumber(10)
44  // , m_Kernel(...)
45  ,
46  m_NumberOfComponentsPerPixel(0)
47  // , m_JointImage(0)
48  // , m_ModeTable(0)
49  ,
50  m_ModeSearch(false),
51  m_ThreadIdNumberOfBits(0)
52 #if 0
53  , m_BucketOptimization(false)
54 #endif
55 {
56  this->SetNumberOfRequiredOutputs(4);
57  this->SetNthOutput(0, OutputImageType::New());
58  this->SetNthOutput(1, OutputSpatialImageType::New());
59  this->SetNthOutput(2, OutputIterationImageType::New());
60  this->SetNthOutput(3, OutputLabelImageType::New());
61  m_GlobalShift.Fill(0);
62 }
63 
64 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
65 MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::~MeanShiftSmoothingImageFilter()
66 {
67 }
68 
69 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
70 const typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputSpatialImageType*
71 MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetSpatialOutput() const
72 {
73  return static_cast<const OutputSpatialImageType*>(this->itk::ProcessObject::GetOutput(1));
74 }
75 
76 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
77 typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputSpatialImageType*
78 MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetSpatialOutput()
79 {
80  return static_cast<OutputSpatialImageType*>(this->itk::ProcessObject::GetOutput(1));
81 }
82 
83 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
84 const typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputImageType*
85 MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetRangeOutput() const
86 {
87  return static_cast<const OutputImageType*>(this->itk::ProcessObject::GetOutput(0));
88 }
89 
90 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
91 typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputImageType*
92 MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetRangeOutput()
93 {
94  return static_cast<OutputImageType*>(this->itk::ProcessObject::GetOutput(0));
95 }
96 
97 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
98 typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputIterationImageType*
99 MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetIterationOutput()
100 {
101  return static_cast<OutputIterationImageType*>(this->itk::ProcessObject::GetOutput(2));
102 }
103 
104 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
105 const typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputIterationImageType*
106 MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetIterationOutput() const
107 {
108  return static_cast<OutputIterationImageType*>(this->itk::ProcessObject::GetOutput(2));
109 }
110 
111 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
112 typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputLabelImageType*
113 MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetLabelOutput()
114 {
115  return static_cast<OutputLabelImageType*>(this->itk::ProcessObject::GetOutput(3));
116 }
117 
118 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
119 const typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputLabelImageType*
120 MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetLabelOutput() const
121 {
122  return static_cast<OutputLabelImageType*>(this->itk::ProcessObject::GetOutput(3));
123 }
124 
125 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
126 void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::AllocateOutputs()
127 {
128  typename OutputSpatialImageType::Pointer spatialOutputPtr = this->GetSpatialOutput();
129  typename OutputImageType::Pointer rangeOutputPtr = this->GetRangeOutput();
130  typename OutputIterationImageType::Pointer iterationOutputPtr = this->GetIterationOutput();
131  typename OutputLabelImageType::Pointer labelOutputPtr = this->GetLabelOutput();
132 
133  spatialOutputPtr->SetBufferedRegion(spatialOutputPtr->GetRequestedRegion());
134  spatialOutputPtr->Allocate();
135 
136  rangeOutputPtr->SetBufferedRegion(rangeOutputPtr->GetRequestedRegion());
137  rangeOutputPtr->Allocate();
138 
139  iterationOutputPtr->SetBufferedRegion(iterationOutputPtr->GetRequestedRegion());
140  iterationOutputPtr->Allocate();
141 
142  labelOutputPtr->SetBufferedRegion(labelOutputPtr->GetRequestedRegion());
143  labelOutputPtr->Allocate();
144 }
145 
146 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
147 void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GenerateOutputInformation()
148 {
149  Superclass::GenerateOutputInformation();
150 
151  m_NumberOfComponentsPerPixel = this->GetInput()->GetNumberOfComponentsPerPixel();
152 
153  if (this->GetSpatialOutput())
154  {
155  this->GetSpatialOutput()->SetNumberOfComponentsPerPixel(ImageDimension); // image lattice
156  }
157  if (this->GetRangeOutput())
158  {
159  this->GetRangeOutput()->SetNumberOfComponentsPerPixel(m_NumberOfComponentsPerPixel);
160  }
161 }
162 
163 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
164 void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GenerateInputRequestedRegion()
165 {
166  // Call superclass implementation
167  Superclass::GenerateInputRequestedRegion();
168 
169  // Retrieve input pointers
170  InputImagePointerType inPtr = const_cast<TInputImage*>(this->GetInput());
171  OutputImagePointerType outRangePtr = this->GetRangeOutput();
172 
173  // Check pointers before using them
174  if (!inPtr || !outRangePtr)
175  {
176  return;
177  }
178 
179  // Retrieve requested region (TODO: check if we need to handle
180  // region for outHDispPtr)
181  RegionType outputRequestedRegion = outRangePtr->GetRequestedRegion();
182 
183  // Pad by the appropriate radius
184  RegionType inputRequestedRegion = outputRequestedRegion;
185 
186  // Initializes the spatial radius from kernel bandwidth
187  m_SpatialRadius.Fill(m_Kernel.GetRadius(m_SpatialBandwidth));
188 
189  InputSizeType margin;
190 
191  for (unsigned int comp = 0; comp < ImageDimension; ++comp)
192  {
193  margin[comp] = (m_MaxIterationNumber * m_SpatialRadius[comp]) + 1;
194  }
195 
196  inputRequestedRegion.PadByRadius(margin);
197 
198  // Crop the input requested region at the input's largest possible region
199  if (inputRequestedRegion.Crop(inPtr->GetLargestPossibleRegion()))
200  {
201  inPtr->SetRequestedRegion(inputRequestedRegion);
202  return;
203  }
204  else
205  {
206  // Couldn't crop the region (requested region is outside the largest
207  // possible region). Throw an exception.
208 
209  // store what we tried to request (prior to trying to crop)
210  inPtr->SetRequestedRegion(inputRequestedRegion);
211 
212  // build an exception
213  itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
214  e.SetLocation(ITK_LOCATION);
215  e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
216  e.SetDataObject(inPtr);
217  throw e;
218  }
219 }
220 
221 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
222 void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::BeforeThreadedGenerateData()
223 {
224  // typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputIteratorWithIndexType;
225  // typedef itk::ImageRegionIterator<RealVectorImageType> JointImageIteratorType;
226 
227  OutputSpatialImagePointerType outSpatialPtr = this->GetSpatialOutput();
228  OutputImagePointerType outRangePtr = this->GetRangeOutput();
229  typename InputImageType::ConstPointer inputPtr = this->GetInput();
230  typename OutputIterationImageType::Pointer iterationOutput = this->GetIterationOutput();
231  typename OutputSpatialImageType::Pointer spatialOutput = this->GetSpatialOutput();
232 
233  // InputIndexType index;
234 
235 
236  m_SpatialRadius.Fill(m_Kernel.GetRadius(m_SpatialBandwidth));
237 
238  m_NumberOfComponentsPerPixel = this->GetInput()->GetNumberOfComponentsPerPixel();
239 
240  // Allocate output images
241  this->AllocateOutputs();
242 
243  // Initialize output images to zero
244  iterationOutput->FillBuffer(0);
245  OutputSpatialPixelType zero(spatialOutput->GetNumberOfComponentsPerPixel());
246  zero.Fill(0);
247  spatialOutput->FillBuffer(zero);
248 
249  // m_JointImage is the input data expressed in the joint spatial-range
250  // domain, i.e. spatial coordinates are concatenated to the range values.
251  // Moreover, pixel components in this image are normalized by their respective
252  // (spatial or range) bandwidth.
253  typedef Meanshift::SpatialRangeJointDomainTransform<InputImageType, RealVectorImageType> FunctionType;
254  typedef otb::UnaryFunctorWithIndexWithOutputSizeImageFilter<InputImageType, RealVectorImageType, FunctionType> JointImageFunctorType;
255 
256  typename JointImageFunctorType::Pointer jointImageFunctor = JointImageFunctorType::New();
257 
258  jointImageFunctor->SetInput(inputPtr);
259  jointImageFunctor->GetFunctor().Initialize(ImageDimension, m_NumberOfComponentsPerPixel, m_GlobalShift);
260  jointImageFunctor->GetOutput()->SetRequestedRegion(this->GetInput()->GetBufferedRegion());
261  jointImageFunctor->Update();
262  m_JointImage = jointImageFunctor->GetOutput();
263 
264 #if 0
265  if (m_BucketOptimization)
266  {
267  // Create bucket image
268  // Note: because values in the input m_JointImage are normalized, the
269  // rangeRadius argument is just 1
270  m_BucketImage = BucketImageType(static_cast<typename RealVectorImageType::ConstPointer> (m_JointImage),
271  m_JointImage->GetRequestedRegion(), m_Kernel.GetRadius(m_SpatialBandwidth), 1,
272  ImageDimension);
273  }
274 #endif
275  /*
276  // Allocate the joint domain image
277  m_JointImage = RealVectorImageType::New();
278  m_JointImage->SetNumberOfComponentsPerPixel(ImageDimension + m_NumberOfComponentsPerPixel);
279  m_JointImage->SetRegions(inputPtr->GetRequestedRegion());
280  m_JointImage->Allocate();
281 
282  InputIteratorWithIndexType inputIt(inputPtr, inputPtr->GetRequestedRegion());
283  JointImageIteratorType jointIt(m_JointImage, inputPtr->GetRequestedRegion());
284 
285  // Initialize the joint image with scaled values
286  inputIt.GoToBegin();
287  jointIt.GoToBegin();
288 
289  while (!inputIt.IsAtEnd())
290  {
291  typename InputImageType::PixelType const& inputPixel = inputIt.Get();
292  index = inputIt.GetIndex();
293 
294  RealVector & jointPixel = jointIt.Get();
295  for(unsigned int comp = 0; comp < ImageDimension; comp++)
296  {
297  jointPixel[comp] = index[comp] / m_SpatialBandwidth;
298  }
299  for(unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++)
300  {
301  jointPixel[ImageDimension + comp] = inputPixel[comp] / m_RangeBandwidth;
302  }
303  // jointIt.Set(jointPixel);
304 
305  ++inputIt;
306  ++jointIt;
307  }
308  */
309 
310  // TODO don't create mode table iterator when ModeSearch is set to false
311  m_ModeTable = ModeTableImageType::New();
312  m_ModeTable->SetRegions(inputPtr->GetRequestedRegion());
313  m_ModeTable->Allocate();
314  m_ModeTable->FillBuffer(0);
315 
316  if (m_ModeSearch)
317  {
318  // Image to store the status at each pixel:
319  // 0 : no mode has been found yet
320  // 1 : a mode has been assigned to this pixel
321  // 2 : a mode will be assigned to this pixel
322 
323 
324  // Initialize counters for mode (also used for mode labeling)
325  // Most significant bits of label counters are used to identify the thread
326  // Id.
327  unsigned int numThreads;
328 
329  numThreads = this->GetNumberOfThreads();
330  m_ThreadIdNumberOfBits = -1;
331  unsigned int n = numThreads;
332  while (n != 0)
333  {
334  n >>= 1;
335  m_ThreadIdNumberOfBits++;
336  }
337  if (m_ThreadIdNumberOfBits == 0)
338  m_ThreadIdNumberOfBits = 1; // minimum 1 bit
339  m_NumLabels.SetSize(numThreads);
340  for (unsigned int i = 0; i < numThreads; i++)
341  {
342  m_NumLabels[i] = static_cast<LabelType>(i) << (sizeof(LabelType) * 8 - m_ThreadIdNumberOfBits);
343  }
344  }
345 }
346 
347 // Calculates the mean shift vector at the position given by jointPixel
348 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
349 void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::CalculateMeanShiftVector(
350  const typename RealVectorImageType::Pointer jointImage, const RealVector& jointPixel, const OutputRegionType& outputRegion, const RealVector& bandwidth,
351  RealVector& meanShiftVector)
352 {
353  const unsigned int jointDimension = ImageDimension + m_NumberOfComponentsPerPixel;
354 
355  InputIndexType inputIndex;
356  InputIndexType regionIndex;
357  InputSizeType regionSize;
358 
359  assert(meanShiftVector.GetSize() == jointDimension);
360  meanShiftVector.Fill(0);
361 
362  // Calculates current pixel neighborhood region, restricted to the output image region
363  for (unsigned int comp = 0; comp < ImageDimension; ++comp)
364  {
365  inputIndex[comp] = std::floor(jointPixel[comp] + 0.5) - m_GlobalShift[comp];
366 
367  regionIndex[comp] =
368  std::max(static_cast<long int>(outputRegion.GetIndex().GetElement(comp)), static_cast<long int>(inputIndex[comp] - m_SpatialRadius[comp] - 1));
369  const long int indexRight = std::min(static_cast<long int>(outputRegion.GetIndex().GetElement(comp) + outputRegion.GetSize().GetElement(comp) - 1),
370  static_cast<long int>(inputIndex[comp] + m_SpatialRadius[comp] + 1));
371 
372  regionSize[comp] = std::max(0l, indexRight - static_cast<long int>(regionIndex[comp]) + 1);
373  }
374 
375  RegionType neighborhoodRegion;
376  neighborhoodRegion.SetIndex(regionIndex);
377  neighborhoodRegion.SetSize(regionSize);
378 
379  RealType weightSum = 0;
380  RealVector shifts(jointDimension);
381 
382  // An iterator on the neighborhood of the current pixel (in joint
383  // spatial-range domain)
384  otb::Meanshift::FastImageRegionConstIterator<RealVectorImageType> it(jointImage, neighborhoodRegion);
385  // itk::ImageRegionConstIterator<RealVectorImageType> it(jointImage, neighborhoodRegion);
386 
387  it.GoToBegin();
388  while (!it.IsAtEnd())
389  {
390  const RealType* jointNeighbor = it.GetPixelPointer();
391 
392  // Compute the squared norm of the difference
393  // This is the L2 norm, TODO: replace by the templated norm
394  RealType norm2 = 0;
395  for (unsigned int comp = 0; comp < jointDimension; comp++)
396  {
397  shifts[comp] = jointNeighbor[comp] - jointPixel[comp];
398  double d = shifts[comp] / bandwidth[comp];
399  norm2 += d * d;
400  }
401 
402  // Compute pixel weight from kernel
403  const RealType weight = m_Kernel(norm2);
404  /*
405  // The following code is an alternative way to compute norm2 and weight
406  // It separates the norms of spatial and range elements
407  RealType spatialNorm2;
408  RealType rangeNorm2;
409  spatialNorm2 = 0;
410  for (unsigned int comp = 0; comp < ImageDimension; comp++)
411  {
412  RealType d;
413  d = jointNeighbor[comp] - jointPixel[comp];
414  spatialNorm2 += d*d;
415  }
416 
417  if(spatialNorm2 >= 1.0)
418  {
419  weight = 0;
420  }
421  else
422  {
423  rangeNorm2 = 0;
424  for (unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++)
425  {
426  RealType d;
427  d = jointNeighbor[ImageDimension + comp] - jointPixel[ImageDimension + comp];
428  rangeNorm2 += d*d;
429  }
430 
431  weight = (rangeNorm2 <= 1.0)? 1.0 : 0.0;
432  }
433  */
434 
435  // Update sum of weights
436  weightSum += weight;
437 
438  // Update mean shift vector
439  for (unsigned int comp = 0; comp < jointDimension; comp++)
440  {
441  meanShiftVector[comp] += weight * shifts[comp];
442  }
443 
444  ++it;
445  }
446 
447  if (weightSum > 0)
448  {
449  for (unsigned int comp = 0; comp < jointDimension; comp++)
450  {
451  meanShiftVector[comp] = meanShiftVector[comp] / weightSum;
452  }
453  }
454 }
455 
456 #if 0
457 // Calculates the mean shift vector at the position given by jointPixel
458 template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
459 void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::CalculateMeanShiftVectorBucket(
460  const RealVector& jointPixel,
461  RealVector& meanShiftVector)
462 {
463  const unsigned int jointDimension = ImageDimension + m_NumberOfComponentsPerPixel;
464 
465  RealType weightSum = 0;
466 
467  for (unsigned int comp = 0; comp < jointDimension; comp++)
468  {
469  meanShiftVector[comp] = 0;
470  }
471 
472  RealVector jointNeighbor(ImageDimension + m_NumberOfComponentsPerPixel);
473 
474  InputIndexType index;
475  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
476  {
477  index[dim] = jointPixel[dim] * m_SpatialBandwidth + 0.5;
478  }
479 
480  const std::vector<unsigned int>
481  neighborBuckets = m_BucketImage.GetNeighborhoodBucketListIndices(
482  m_BucketImage.BucketIndexToBucketListIndex(
483  m_BucketImage.GetBucketIndex(
484  jointPixel,
485  index)));
486 
487  unsigned int numNeighbors = m_BucketImage.GetNumberOfNeighborBuckets();
488  for (unsigned int neighborIndex = 0; neighborIndex < numNeighbors; ++neighborIndex)
489  {
490  const typename BucketImageType::BucketType & bucket = m_BucketImage.GetBucket(neighborBuckets[neighborIndex]);
491  if (bucket.empty()) continue;
492  typename BucketImageType::BucketType::const_iterator it = bucket.begin();
493  while (it != bucket.end())
494  {
495  jointNeighbor.SetData(const_cast<RealType*> (*it));
496 
497  // Compute the squared norm of the difference
498  // This is the L2 norm, TODO: replace by the templated norm
499  RealType norm2 = 0;
500  for (unsigned int comp = 0; comp < jointDimension; comp++)
501  {
502  const RealType d = jointNeighbor[comp] - jointPixel[comp];
503  norm2 += d * d;
504  }
505 
506  // Compute pixel weight from kernel
507  const RealType weight = m_Kernel(norm2);
508 
509  // Update sum of weights
510  weightSum += weight;
511 
512  // Update mean shift vector
513  for (unsigned int comp = 0; comp < jointDimension; comp++)
514  {
515  meanShiftVector[comp] += weight * jointNeighbor[comp];
516  }
517 
518  ++it;
519  }
520  }
521 
522  if (weightSum > 0)
523  {
524  for (unsigned int comp = 0; comp < jointDimension; comp++)
525  {
526  meanShiftVector[comp] = meanShiftVector[comp] / weightSum - jointPixel[comp];
527  }
528  }
529 }
530 #endif
531 
532 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
533 void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::ThreadedGenerateData(
534  const OutputRegionType& outputRegionForThread, itk::ThreadIdType threadId)
535 {
536  // at the first iteration
537 
538 
539  // Retrieve output images pointers
540  typename OutputSpatialImageType::Pointer spatialOutput = this->GetSpatialOutput();
541  typename OutputImageType::Pointer rangeOutput = this->GetRangeOutput();
542  typename OutputIterationImageType::Pointer iterationOutput = this->GetIterationOutput();
543  typename OutputLabelImageType::Pointer labelOutput = this->GetLabelOutput();
544 
545  // Get input image pointer
546  typename InputImageType::ConstPointer input = this->GetInput();
547 
548  // defines input and output iterators
549  typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType;
550  typedef itk::ImageRegionIterator<OutputSpatialImageType> OutputSpatialIteratorType;
551  typedef itk::ImageRegionIterator<OutputIterationImageType> OutputIterationIteratorType;
552  typedef itk::ImageRegionIterator<OutputLabelImageType> OutputLabelIteratorType;
553 
554  const unsigned int jointDimension = ImageDimension + m_NumberOfComponentsPerPixel;
555 
556  typename OutputImageType::PixelType rangePixel(m_NumberOfComponentsPerPixel);
557  typename OutputSpatialImageType::PixelType spatialPixel(ImageDimension);
558 
559  RealVector jointPixel(jointDimension);
560 
561  RealVector bandwidth(jointDimension);
562  for (unsigned int comp = 0; comp < ImageDimension; comp++)
563  bandwidth[comp] = m_SpatialBandwidth;
564 
565  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
566 
567  RegionType const& requestedRegion = input->GetRequestedRegion();
568 
569  typedef itk::ImageRegionConstIteratorWithIndex<RealVectorImageType> JointImageIteratorType;
570  JointImageIteratorType jointIt(m_JointImage, outputRegionForThread);
571 
572  OutputIteratorType rangeIt(rangeOutput, outputRegionForThread);
573  OutputSpatialIteratorType spatialIt(spatialOutput, outputRegionForThread);
574  OutputIterationIteratorType iterationIt(iterationOutput, outputRegionForThread);
575  OutputLabelIteratorType labelIt(labelOutput, outputRegionForThread);
576 
577  typedef itk::ImageRegionIterator<ModeTableImageType> ModeTableImageIteratorType;
578  ModeTableImageIteratorType modeTableIt(m_ModeTable, outputRegionForThread);
579 
580  jointIt.GoToBegin();
581  rangeIt.GoToBegin();
582  spatialIt.GoToBegin();
583  iterationIt.GoToBegin();
584  modeTableIt.GoToBegin();
585  labelIt.GoToBegin();
586 
587  unsigned int iteration = 0;
588 
589  // Mean shift vector, updating the joint pixel at each iteration
590  RealVector meanShiftVector(jointDimension);
591 
592  // Variables used by mode search optimization
593  // List of indices where the current pixel passes through
594  std::vector<InputIndexType> pointList;
595  if (m_ModeSearch)
596  pointList.resize(m_MaxIterationNumber);
597  // Number of times an already processed candidate pixel is encountered, resulting in no
598  // further computation (Used for statistics only)
599  unsigned int numBreaks = 0;
600  // index of the current pixel updated during the mean shift loop
601  InputIndexType modeCandidate;
602 
603  for (; !jointIt.IsAtEnd(); ++jointIt, ++rangeIt, ++spatialIt, ++iterationIt, ++modeTableIt, ++labelIt, progress.CompletedPixel())
604  {
605 
606  // if pixel has been already processed (by mode search optimization), skip
607  typename ModeTableImageType::InternalPixelType const& currentPixelMode = modeTableIt.Get();
608  if (m_ModeSearch && currentPixelMode == 1)
609  {
610  numBreaks++;
611  continue;
612  }
613 
614  bool hasConverged = false;
615 
616  // get input pixel in the joint spatial-range domain (with components
617  // normalized by bandwidth)
618  const RealVector& jointPixelVal = jointIt.Get(); // Pixel in the joint spatial-range domain
619  for (unsigned int comp = 0; comp < jointDimension; comp++)
620  jointPixel[comp] = jointPixelVal[comp];
621 
622  for (unsigned int comp = ImageDimension; comp < jointDimension; comp++)
623  bandwidth[comp] = m_RangeBandwidthRamp * jointPixel[comp] + m_RangeBandwidth;
624 
625  // index of the currently processed output pixel
626  InputIndexType currentIndex = jointIt.GetIndex();
627 
628  // Number of points currently in the pointList
629  unsigned int pointCount = 0; // Note: used only in mode search optimization
630  iteration = 0;
631  while ((iteration < m_MaxIterationNumber) && (!hasConverged))
632  {
633 
634  if (m_ModeSearch)
635  {
636  // Find index of the pixel closest to the current jointPixel (not normalized by bandwidth)
637  for (unsigned int comp = 0; comp < ImageDimension; comp++)
638  {
639  modeCandidate[comp] = std::floor(jointPixel[comp] - m_GlobalShift[comp] + 0.5);
640  }
641  // Check status of candidate mode
642 
643  // If pixel candidate has status 0 (no mode assigned) or 1 (mode assigned)
644  // but not 2 (pixel in current search path), and pixel has actually moved
645  // from its initial position, and pixel candidate is inside the output
646  // region, then perform optimization tasks
647  if (modeCandidate != currentIndex && m_ModeTable->GetPixel(modeCandidate) != 2 && outputRegionForThread.IsInside(modeCandidate))
648  {
649  // Obtain the data point to see if it close to jointPixel
650  RealType diff = 0;
651  RealVector const& candidatePixel = m_JointImage->GetPixel(modeCandidate);
652  for (unsigned int comp = ImageDimension; comp < jointDimension; comp++)
653  {
654  const RealType d = (candidatePixel[comp] - jointPixel[comp]) / bandwidth[comp];
655  diff += d * d;
656  }
657 
658  if (diff < 0.5) // Spectral value is close enough
659  {
660  // If no mode has been associated to the candidate pixel then
661  // associate it to the upcoming mode
662  if (m_ModeTable->GetPixel(modeCandidate) == 0)
663  {
664  // Add the candidate to the list of pixels that will be assigned the
665  // finally calculated mode value
666  pointList[pointCount++] = modeCandidate;
667  m_ModeTable->SetPixel(modeCandidate, 2);
668  }
669  else // == 1
670  {
671  // The candidate pixel has already been assigned to a mode
672  // Assign the same value
673  rangePixel = rangeOutput->GetPixel(modeCandidate);
674  for (unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++)
675  {
676  jointPixel[ImageDimension + comp] = rangePixel[comp];
677  }
678  // Update the mode table because pixel will be assigned just now
679  modeTableIt.Set(2); // m_ModeTable->SetPixel(currentIndex, 2);
680  // bypass further calculation
681  numBreaks++;
682  break;
683  }
684  }
685  }
686  } // end if (m_ModeSearch)
687 
688 // Calculate meanShiftVector
689 #if 0
690  if (m_BucketOptimization)
691  {
692  this->CalculateMeanShiftVectorBucket(jointPixel, meanShiftVector);
693  }
694  else
695  {
696 #endif
697  this->CalculateMeanShiftVector(m_JointImage, jointPixel, requestedRegion, bandwidth, meanShiftVector);
698 
699 #if 0
700  }
701 #endif
702 
703  // Compute mean shift vector squared norm (not normalized by bandwidth)
704  // and add mean shift vector to current joint pixel
705  double meanShiftVectorSqNorm = 0;
706  for (unsigned int comp = 0; comp < jointDimension; comp++)
707  {
708  const double v = meanShiftVector[comp];
709  meanShiftVectorSqNorm += v * v;
710  jointPixel[comp] += meanShiftVector[comp];
711  }
712 
713  // TODO replace SSD Test with templated metric
714  hasConverged = meanShiftVectorSqNorm < m_Threshold;
715  iteration++;
716  }
717 
718  for (unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++)
719  {
720  rangePixel[comp] = jointPixel[ImageDimension + comp];
721  }
722 
723  for (unsigned int comp = 0; comp < ImageDimension; comp++)
724  {
725  spatialPixel[comp] = jointPixel[comp] - currentIndex[comp] - m_GlobalShift[comp];
726  }
727 
728  rangeIt.Set(rangePixel);
729  spatialIt.Set(spatialPixel);
730 
731  const typename OutputIterationImageType::PixelType iterationPixel = iteration;
732  iterationIt.Set(iterationPixel);
733 
734  if (m_ModeSearch)
735  {
736  // Update the mode table now that the current pixel has been assigned
737  modeTableIt.Set(1); // m_ModeTable->SetPixel(currentIndex, 1);
738 
739  // If the loop exited with hasConverged or too many iterations, then we have a new mode
740  LabelType label;
741  if (hasConverged || iteration == m_MaxIterationNumber)
742  {
743  m_NumLabels[threadId]++;
744  label = m_NumLabels[threadId];
745  }
746  else // the loop exited through a break. Use the already assigned mode label
747  {
748  label = labelOutput->GetPixel(modeCandidate);
749  }
750  labelIt.Set(label);
751 
752  // Also assign all points in the list to the same mode
753  for (unsigned int i = 0; i < pointCount; i++)
754  {
755  rangeOutput->SetPixel(pointList[i], rangePixel);
756  m_ModeTable->SetPixel(pointList[i], 1);
757  labelOutput->SetPixel(pointList[i], label);
758  }
759  }
760  else // if ModeSearch is not set LabelOutput can't be generated
761  {
762  LabelType labelZero = 0;
763  labelIt.Set(labelZero);
764  }
765  }
766  // std::cout << "numBreaks: " << numBreaks << " Break ratio: " << numBreaks / (RealType)outputRegionForThread.GetNumberOfPixels() << std::endl;
767 }
768 
769 /* after threaded convergence test */
770 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
771 void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::AfterThreadedGenerateData()
772 {
773  typename OutputLabelImageType::Pointer labelOutput = this->GetLabelOutput();
774  typedef itk::ImageRegionIterator<OutputLabelImageType> OutputLabelIteratorType;
775  OutputLabelIteratorType labelIt(labelOutput, labelOutput->GetRequestedRegion());
776 
777  // Reassign mode labels
778  // Note: Labels are only computed when mode search optimization is enabled
779  if (m_ModeSearch)
780  {
781  // New labels will be consecutive. The following vector contains the new
782  // start label for each thread.
783  itk::VariableLengthVector<LabelType> newLabelOffset;
784  newLabelOffset.SetSize(this->GetNumberOfThreads());
785  newLabelOffset[0] = 0;
786  for (itk::ThreadIdType i = 1; i < this->GetNumberOfThreads(); i++)
787  {
788  // Retrieve the number of labels in the thread by removing the threadId
789  // from the most significant bits
790  LabelType localNumLabel =
791  m_NumLabels[i - 1] & ((static_cast<LabelType>(1) << (sizeof(LabelType) * 8 - m_ThreadIdNumberOfBits)) - static_cast<LabelType>(1));
792  newLabelOffset[i] = localNumLabel + newLabelOffset[i - 1];
793  }
794 
795  labelIt.GoToBegin();
796 
797  while (!labelIt.IsAtEnd())
798  {
799  LabelType const label = labelIt.Get();
800 
801  // Get threadId from most significant bits
802  const itk::ThreadIdType threadId = label >> (sizeof(LabelType) * 8 - m_ThreadIdNumberOfBits);
803 
804  // Relabeling
805  // First get the label number by removing the threadId bits
806  // Then add the label offset specific to the threadId
807  LabelType newLabel = label & ((static_cast<LabelType>(1) << (sizeof(LabelType) * 8 - m_ThreadIdNumberOfBits)) - static_cast<LabelType>(1));
808  newLabel += newLabelOffset[threadId];
809 
810  labelIt.Set(newLabel);
811  ++labelIt;
812  }
813  }
814 }
815 
816 template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage>
817 void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
818 {
819  Superclass::PrintSelf(os, indent);
820  os << indent << "Spatial bandwidth: " << m_SpatialBandwidth << std::endl;
821  os << indent << "Range bandwidth: " << m_RangeBandwidth << std::endl;
822 }
823 
824 } // end namespace otb
825 
826 #endif
otb::MeanShiftSmoothingImageFilter::OutputRegionType
OutputImageType::RegionType OutputRegionType
Definition: otbMeanShiftSmoothingImageFilter.h:499
otb::MeanShiftSmoothingImageFilter::InputSizeType
InputImageType::SizeType InputSizeType
Definition: otbMeanShiftSmoothingImageFilter.h:490
otb::MeanShiftSmoothingImageFilter::OutputImageType
TOutputImage OutputImageType
Definition: otbMeanShiftSmoothingImageFilter.h:496
otb::MeanShiftSmoothingImageFilter::m_GlobalShift
InputIndexType m_GlobalShift
Definition: otbMeanShiftSmoothingImageFilter.h:698
otb::VectorImage::Pointer
itk::SmartPointer< Self > Pointer
Definition: otbVectorImage.h:53
otb::MeanShiftSmoothingImageFilter::GetIterationOutput
const OutputIterationImageType * GetIterationOutput() const
Definition: otbMeanShiftSmoothingImageFilter.hxx:99
otb::MeanShiftSmoothingImageFilter::GenerateInputRequestedRegion
void GenerateInputRequestedRegion() override
Definition: otbMeanShiftSmoothingImageFilter.hxx:164
otb::MeanShiftSmoothingImageFilter::RealVector
itk::VariableLengthVector< RealType > RealVector
Definition: otbMeanShiftSmoothingImageFilter.h:514
otb::VectorImage::ConstPointer
itk::SmartPointer< const Self > ConstPointer
Definition: otbVectorImage.h:54
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::Image
Creation of an "otb" image which contains metadata.
Definition: otbImage.h:89
otb::MeanShiftSmoothingImageFilter::InputIndexType
InputImageType::IndexType InputIndexType
Definition: otbMeanShiftSmoothingImageFilter.h:489
otb::MeanShiftSmoothingImageFilter::BeforeThreadedGenerateData
void BeforeThreadedGenerateData() override
Definition: otbMeanShiftSmoothingImageFilter.hxx:222
otb::MeanShiftSmoothingImageFilter::ThreadedGenerateData
void ThreadedGenerateData(const OutputRegionType &outputRegionForThread, itk::ThreadIdType threadId) override
Definition: otbMeanShiftSmoothingImageFilter.hxx:533
otbMacro.h
otb::MeanShiftSmoothingImageFilter::GetRangeOutput
const OutputImageType * GetRangeOutput() const
Definition: otbMeanShiftSmoothingImageFilter.hxx:85
otb::UnaryFunctorWithIndexWithOutputSizeImageFilter
Implements neighborhood-wise generic operation on image.
Definition: otbUnaryFunctorWithIndexWithOutputSizeImageFilter.h:42
otb::MeanShiftSmoothingImageFilter::OutputSpatialPixelType
OutputSpatialImageType::PixelType OutputSpatialPixelType
Definition: otbMeanShiftSmoothingImageFilter.h:508
otb::MeanShiftSmoothingImageFilter::OutputIterationImageType
TOutputIterationImage OutputIterationImageType
Definition: otbMeanShiftSmoothingImageFilter.h:501
otb::MeanShiftSmoothingImageFilter::MeanShiftSmoothingImageFilter
MeanShiftSmoothingImageFilter()
Definition: otbMeanShiftSmoothingImageFilter.hxx:36
otb::MeanShiftSmoothingImageFilter::OutputImagePointerType
OutputImageType::Pointer OutputImagePointerType
Definition: otbMeanShiftSmoothingImageFilter.h:497
otb::MeanShiftSmoothingImageFilter::GetLabelOutput
const OutputLabelImageType * GetLabelOutput() const
Definition: otbMeanShiftSmoothingImageFilter.hxx:113
otb::Meanshift::FastImageRegionConstIterator::GetPixelPointer
const InternalPixelType * GetPixelPointer() const
Definition: otbMeanShiftSmoothingImageFilter.h:172
otb::Image::InternalPixelType
Superclass::InternalPixelType InternalPixelType
Definition: otbImage.h:120
otb::Image::New
static Pointer New()
otb::MeanShiftSmoothingImageFilter::~MeanShiftSmoothingImageFilter
~MeanShiftSmoothingImageFilter() override
Definition: otbMeanShiftSmoothingImageFilter.hxx:65
otb::MeanShiftSmoothingImageFilter::RegionType
InputImageType::RegionType RegionType
Definition: otbMeanShiftSmoothingImageFilter.h:493
otb::MeanShiftSmoothingImageFilter::GenerateOutputInformation
void GenerateOutputInformation(void) override
Definition: otbMeanShiftSmoothingImageFilter.hxx:147
otb::VectorImage::PixelType
Superclass::PixelType PixelType
Definition: otbVectorImage.h:63
otb::MeanShiftSmoothingImageFilter::OutputSpatialImagePointerType
OutputSpatialImageType::Pointer OutputSpatialImagePointerType
Definition: otbMeanShiftSmoothingImageFilter.h:507
otb::MeanShiftSmoothingImageFilter::GetSpatialOutput
const OutputSpatialImageType * GetSpatialOutput() const
Definition: otbMeanShiftSmoothingImageFilter.hxx:71
otb::MeanShiftSmoothingImageFilter::AfterThreadedGenerateData
void AfterThreadedGenerateData() override
Definition: otbMeanShiftSmoothingImageFilter.hxx:771
otb::VectorImage::New
static Pointer New()
otb::MeanShiftSmoothingImageFilter
Definition: otbMeanShiftSmoothingImageFilter.h:467
otb::MeanShiftSmoothingImageFilter::LabelType
unsigned long LabelType
Definition: otbMeanShiftSmoothingImageFilter.h:503
otb::MeanShiftSmoothingImageFilter::InputImagePointerType
InputImageType::Pointer InputImagePointerType
Definition: otbMeanShiftSmoothingImageFilter.h:487
otb::MeanShiftSmoothingImageFilter::RealType
double RealType
Definition: otbMeanShiftSmoothingImageFilter.h:475
otbMeanShiftSmoothingImageFilter.h
otb::MeanShiftSmoothingImageFilter::AllocateOutputs
void AllocateOutputs() override
Definition: otbMeanShiftSmoothingImageFilter.hxx:126
otb::Image::Pointer
itk::SmartPointer< Self > Pointer
Definition: otbImage.h:97
otb::Meanshift::SpatialRangeJointDomainTransform
Definition: otbMeanShiftSmoothingImageFilter.h:58
otb::Meanshift::FastImageRegionConstIterator
Definition: otbMeanShiftSmoothingImageFilter.h:148
otb::MeanShiftSmoothingImageFilter::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbMeanShiftSmoothingImageFilter.hxx:817
otb::VectorImage
Creation of an "otb" vector image which contains metadata.
Definition: otbVectorImage.h:45
otb::MeanShiftSmoothingImageFilter::CalculateMeanShiftVector
virtual void CalculateMeanShiftVector(const typename RealVectorImageType::Pointer inputImagePtr, const RealVector &jointPixel, const OutputRegionType &outputRegion, const RealVector &bandwidth, RealVector &meanShiftVector)
Definition: otbMeanShiftSmoothingImageFilter.hxx:349
otb::MeanShiftSmoothingImageFilter::OutputLabelImageType
otb::Image< LabelType, InputImageType::ImageDimension > OutputLabelImageType
Definition: otbMeanShiftSmoothingImageFilter.h:504
otbUnaryFunctorWithIndexWithOutputSizeImageFilter.h
Generated at Wed Feb 28 2024 15:23:03 for OTB with   doxygen 1.8.17