otbMeanShiftSmoothingImageFilter.h 24.5 KB
Newer Older
1
/*
Julien Michel's avatar
Julien Michel committed
2
 * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
20

21 22
#ifndef otbMeanShiftSmoothingImageFilter_h
#define otbMeanShiftSmoothingImageFilter_h
23

24
#include "otbImage.h"
25
#include "otbVectorImage.h"
26
#include "itkImageToImageFilter.h"
27
#include "itkImageRegionConstIterator.h"
28
#include "itkImageRegionConstIteratorWithIndex.h"
29
#include <algorithm>
30

Jonathan Guinet's avatar
Jonathan Guinet committed
31

32 33
namespace otb
{
34 35 36
namespace Meanshift
{

37
template<typename T> inline T simple_pow(T const& v, unsigned int p)
38 39
{
  T res = 1;
40
  for (unsigned int i = 0; i != p; ++i)
41 42 43 44 45
    {
    res *= v;
    }
  return res;
}
46

47 48 49 50 51 52
/** \class SpatialRangeJointDomainTransform
 *
 *
 * Functor returning the joint spatial-range representation of a pixel, i.e. the
 * concatenation of range components and coordinates as a vector. Components are
 * scaled by their respective spatial and range bandwidth.
53 54
 *
 * \ingroup OTBSmoothing
55
 */
56 57 58 59 60 61
template<class TInputImage, class TOutputJointImage>
class SpatialRangeJointDomainTransform
{
public:
  typedef double RealType;

62 63 64 65
  SpatialRangeJointDomainTransform():
      m_ImageDimension(0),
      m_NumberOfComponentsPerPixel(0),
      m_OutputSize(0)
66 67
  {
  }
68

69 70
  typename TOutputJointImage::PixelType operator()(const typename TInputImage::PixelType & inputPixel,
                                                   const typename TInputImage::IndexType & index) const
71
  {
72
    typename TOutputJointImage::PixelType jointPixel(m_ImageDimension + m_NumberOfComponentsPerPixel);
73

74
    for (unsigned int comp = 0; comp < m_ImageDimension; comp++)
75
      {
Julien Michel's avatar
Julien Michel committed
76
      jointPixel[comp] = index[comp] + m_GlobalShift[comp];
77
      }
78
    for (unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++)
79
      {
80
      jointPixel[m_ImageDimension + comp] = inputPixel[comp];
81 82 83 84
      }
    return jointPixel;
  }

Julien Michel's avatar
Julien Michel committed
85
  void Initialize(unsigned int _ImageDimension, unsigned int numberOfComponentsPerPixel_, typename TInputImage::IndexType globalShift_)
86
  {
87
    m_ImageDimension = _ImageDimension;
88
    m_NumberOfComponentsPerPixel = numberOfComponentsPerPixel_;
89
    m_OutputSize = m_ImageDimension + m_NumberOfComponentsPerPixel;
90
    m_GlobalShift = globalShift_;
91 92
  }

93
  unsigned int GetOutputSize() const
94 95 96 97 98
  {
    return m_OutputSize;
  }

private:
99
  unsigned int m_ImageDimension;
100 101
  unsigned int m_NumberOfComponentsPerPixel;
  unsigned int m_OutputSize;
102
  typename TInputImage::IndexType m_GlobalShift;
103
};
104

105 106 107 108 109
class KernelUniform
{
public:
  typedef double RealType;

110 111
  // KernelUniform() {}
  // ~KernelUniform() {}
112

113 114
  RealType operator()(RealType x) const
  {
Julien Michel's avatar
Julien Michel committed
115
    return (x <= 1) ? 1.0 : 0.0;
116 117
  }

118 119
  RealType GetRadius(RealType bandwidth) const
  {
120
    return bandwidth;
121
  }
122
};
123

124 125 126 127 128
class KernelGaussian
{
public:
  typedef double RealType;

129 130
  // KernelGaussian() {}
  // ~KernelGaussian() {}
131

132 133
  RealType operator()(RealType x) const
  {
134
    return std::exp(-0.5 * x);
135 136
  }

137 138 139
  RealType GetRadius(RealType bandwidth) const
  {
    return 3.0 * bandwidth;
140
  }
141 142
};

143 144 145 146
/** \class FastImageRegionConstIterator
 *
 * Iterator for reading pixels over an image region, specialized for faster
 * access to pixels in vector images through the method GetPixelPointer
147 148
 *
 * \ingroup OTBSmoothing
149 150
 */
template<typename TImage>
151
class FastImageRegionConstIterator: public itk::ImageRegionConstIterator<TImage>
152 153 154
{
public:
  /** Standard class typedef. */
155 156
  typedef FastImageRegionConstIterator<TImage> Self;
  typedef itk::ImageRegionConstIterator<TImage> Superclass;
157

158 159
  typedef typename Superclass::ImageType ImageType;
  typedef typename Superclass::RegionType RegionType;
160

161
  typedef typename TImage::PixelType PixelType;
162 163
  typedef typename TImage::InternalPixelType InternalPixelType;

164
  itkTypeMacro(FastImageRegionConstIterator, ImageRegionConstIterator);
OTB Bot's avatar
STYLE  
OTB Bot committed
165
;
166

167 168 169 170 171 172
  FastImageRegionConstIterator() :
    Superclass()
  {
  }
  FastImageRegionConstIterator(const ImageType *ptr, const RegionType &region) :
    Superclass(ptr, region)
173 174 175 176
  {
    m_NumberOfComponentsPerPixel = ptr->GetNumberOfComponentsPerPixel();
  }

177
  const InternalPixelType * GetPixelPointer() const
178 179 180 181
  {
    return this->m_Buffer + (this->m_Offset * m_NumberOfComponentsPerPixel);
  }

182
private:
183 184 185
  unsigned int m_NumberOfComponentsPerPixel;
};

Jonathan Guinet's avatar
Jonathan Guinet committed
186
#if 0 //disable bucket mode
187
/** \class BucketImage
188 189 190 191 192 193 194 195
 *
 * This class indexes pixels in a N-dimensional image into a N+1-dimensional
 * array of buckets. Each pixel is stored into one bucket depending on its
 * position on the lattice (the width of a bucket is given at construction) and
 * one spectral component (also given at construction by spectralCoordinate).
 *
 * The (spatially and spectrally) neighboring buckets of pixels can then be
 * obtained by using GetNeighborhoodBucketListIndices().
196 197
 *
 * \ingroup OTBSmoothing
198 199
 */
template<class TImage>
200 201 202
class BucketImage
{
public:
203 204 205 206 207 208
  typedef TImage ImageType;
  typedef typename ImageType::ConstPointer ImageConstPointerType;
  typedef typename ImageType::PixelType PixelType;
  typedef typename ImageType::InternalPixelType InternalPixelType;
  typedef typename ImageType::RegionType RegionType;
  typedef typename ImageType::IndexType IndexType;
209

210
  typedef double RealType;
211 212 213 214

  static const unsigned int ImageDimension = ImageType::ImageDimension;

  /** The bucket image has dimension N+1 (ie. usually 3D for most images) */
215
  typedef std::vector<typename ImageType::SizeType::SizeValueType> BucketImageSizeType;
216 217
  typedef std::vector<typename ImageType::IndexType::IndexValueType> BucketImageIndexType;
  //typedef std::vector<long> BucketImageIndexType;
218 219 220


  /** pixel buckets typedefs and declarations */
221 222 223
  typedef const typename ImageType::InternalPixelType * ImageDataPointerType;
  typedef std::vector<ImageDataPointerType> BucketType;
  typedef std::vector<BucketType> BucketListType;
224

225 226 227
  BucketImage()
  {
  }
228 229

  /** Constructor for the bucket image. It operates on the specified
230 231 232 233 234 235 236 237 238 239
   * region.
   * spatialRadius specifies the width of a bucket in pixels.
   * rangeRadius is the spectral width for the specified spectral coordinate in
   * one bucket.
   * spectralCoordinate is the index of the pixel used for classification in buckets.
   */
  BucketImage(ImageConstPointerType image, const RegionType & region, RealType spatialRadius, RealType rangeRadius,
              unsigned int spectralCoordinate) :
    m_Image(image), m_Region(region), m_SpatialRadius(spatialRadius), m_RangeRadius(rangeRadius),
        m_SpectralCoordinate(spectralCoordinate)
240 241 242 243 244 245 246 247
  {

    // Find max and min of the used spectral band
    itk::ImageRegionConstIterator<ImageType> inputIt(m_Image, m_Region);
    inputIt.GoToBegin();
    InternalPixelType minValue = inputIt.Get()[spectralCoordinate];
    InternalPixelType maxValue = minValue;
    ++inputIt;
248
    while (!inputIt.IsAtEnd())
249 250
      {
      const PixelType &p = inputIt.Get();
251 252
      minValue = std::min(minValue, p[m_SpectralCoordinate]);
      maxValue = std::max(maxValue, p[m_SpectralCoordinate]);
253 254 255 256 257
      ++inputIt;
      }

    m_MinValue = minValue;
    m_MaxValue = maxValue;
258 259 260

    // Compute bucket image dimensions. Note: empty buckets are at each border
    // to simplify image border issues
261 262
    m_DimensionVector.resize(ImageDimension + 1); // NB: pays for a 0-innit
    for (unsigned int dim = 0; dim < ImageDimension; ++dim)
263
      {
264
      m_DimensionVector[dim] = m_Region.GetSize()[dim] / m_SpatialRadius + 3;
265
      }
266
    m_DimensionVector[ImageDimension] = (unsigned int) ((maxValue - minValue) / m_RangeRadius) + 3;
267 268

    unsigned int numBuckets = m_DimensionVector[0];
269
    for (unsigned int dim = 1; dim <= ImageDimension; ++dim)
270 271 272 273 274 275 276 277 278 279
      numBuckets *= m_DimensionVector[dim];

    m_BucketList.resize(numBuckets);
    // Build buckets
    itk::ImageRegionConstIteratorWithIndex<ImageType> it(m_Image, m_Region);
    it.GoToBegin();
    // this iterator is only used to get the pixel data pointer
    FastImageRegionConstIterator<ImageType> fastIt(m_Image, m_Region);
    fastIt.GoToBegin();

280
    while (!it.IsAtEnd())
281 282 283 284 285
      {
      const IndexType & index = it.GetIndex();
      const PixelType & pixel = it.Get();

      // Find which bucket this pixel belongs to
286
      const BucketImageIndexType bucketIndex = GetBucketIndex(pixel, index);
287 288 289 290 291 292 293

      unsigned int bucketListIndex = BucketIndexToBucketListIndex(bucketIndex);
      assert(bucketListIndex < numBuckets);
      m_BucketList[bucketListIndex].push_back(fastIt.GetPixelPointer());
      ++it;
      ++fastIt;
      }
294 295

    // Prepare neighborhood offset vector
296
    // BucketImageIndexType zeroOffsetIndex(ImageDimension+1);
297
    std::vector<BucketImageIndexType> neighborsIndexList;
298 299
    neighborsIndexList.reserve(simple_pow(3, ImageDimension + 1));
    neighborsIndexList.resize(1, BucketImageIndexType(ImageDimension + 1)); // zeroOffsetIndex
300
    // neighborsIndexList.push_back(zeroOffsetIndex);
301
    for (unsigned dim = 0; dim <= ImageDimension; ++dim)
302 303 304
      {
      // take all neighbors already in the list and add their direct neighbor
      // along the current dim
305
      const unsigned int curSize = neighborsIndexList.size();
306
      for (unsigned int i = 0; i < curSize; ++i)
307 308 309 310 311 312 313 314 315
        {
        BucketImageIndexType index = neighborsIndexList[i];
        index[dim]--;
        neighborsIndexList.push_back(index);
        index[dim] += 2;
        neighborsIndexList.push_back(index);
        }
      }
    // Convert all neighbors n-dimensional indices to bucket list 1D indices
316 317
    const unsigned int neighborhoodOffsetVectorSize = neighborsIndexList.size();
    m_NeighborhoodOffsetVector.reserve(neighborhoodOffsetVectorSize);
318
    for (unsigned int i = 0; i < neighborhoodOffsetVectorSize; ++i)
319
      {
320
      const int listIndex = BucketIndexToBucketListIndex(neighborsIndexList[i]);
321 322
      m_NeighborhoodOffsetVector.push_back(listIndex);
      }
323 324
  }

325 326 327
  ~BucketImage()
  {
  }
328

329
  /** Returns the N+1-dimensional bucket index for a given pixel value at the given index */
330 331
  BucketImageIndexType GetBucketIndex(const PixelType & pixel, const IndexType & index)
  {
332 333
    BucketImageIndexType bucketIndex(ImageDimension + 1);
    for (unsigned int dim = 0; dim < ImageDimension; ++dim)
334
      {
335
      bucketIndex[dim] = (index[dim] - m_Region.GetIndex()[dim]) / m_SpatialRadius + 1;
336
      }
337
    bucketIndex[ImageDimension] = (pixel[m_SpectralCoordinate] - m_MinValue) / m_RangeRadius + 1;
338 339 340
    return bucketIndex;
  }

341
  /** Converts a N+1-dimensional bucket index into the 1D list index usable by
342
   GetBucket() */
343
  int BucketIndexToBucketListIndex(const BucketImageIndexType & bucketIndex) const
344
  {
345
    int bucketListIndex = bucketIndex[0];
346
    for (unsigned int dim = 1; dim <= ImageDimension; ++dim)
347 348 349 350 351 352
      {
      bucketListIndex = bucketListIndex * m_DimensionVector[dim] + bucketIndex[dim];
      }
    return bucketListIndex;
  }

353
  /** Retrieves the list of all buckets in the neighborhood of the given bucket */
354
  std::vector<unsigned int> GetNeighborhoodBucketListIndices(int bucketIndex) const
355
  {
356 357
    const unsigned int neighborhoodOffsetVectorSize = m_NeighborhoodOffsetVector.size();
    std::vector<unsigned int> indices(neighborhoodOffsetVectorSize);
358

359
    for (unsigned int i = 0; i < neighborhoodOffsetVectorSize; ++i)
360
      {
361
      indices[i] = bucketIndex + m_NeighborhoodOffsetVector[i];
362 363
      }
    return indices;
364
  }
365

366
  /* Returns the list of pixels (actually pointer to pixel data) contained in a bucket */
367
  const BucketType & GetBucket(unsigned int index) const
368 369 370 371
  {
    return m_BucketList[index];
  }

372
  unsigned int GetNumberOfNeighborBuckets() const
373
  {
374
    return m_NeighborhoodOffsetVector.size();
375 376
  }

377 378 379 380 381 382 383 384 385 386 387
private:
  /** Input image */
  ImageConstPointerType m_Image;
  /** Processed region */
  RegionType m_Region;

  /** Spatial radius of one bucket of pixels */
  RealType m_SpatialRadius;
  /** Range radius (at a single dimension) of one bucket of pixels */
  RealType m_RangeRadius;
  /** pixels are separated in buckets depending on their spatial position and
388
   also their value at one coordinate */
389 390 391 392 393 394 395 396 397
  unsigned int m_SpectralCoordinate;
  /** Min and Max of selected spectral coordinate */
  InternalPixelType m_MinValue;
  InternalPixelType m_MaxValue;

  /** the buckets are stored in this list */
  BucketListType m_BucketList;
  /** This vector holds the dimensions of the 3D (ND?) bucket image */
  BucketImageSizeType m_DimensionVector;
398
  /** Vector of offsets in the buckets list to get all buckets in the
399 400
   * neighborhood
   */
401
  std::vector<int> m_NeighborhoodOffsetVector;
402
};
403
#endif
404

405 406
} // end namespace Meanshift

407
/** \class MeanShiftSmoothingImageFilter
408 409
 *
 *
410 411 412 413 414
 * Mean shift is an edge-preserving smoothing algorithm often used in image
 * processing and segmentation. It will iteratively smooth a given pixel with
 * its neighbors that are within a spatial distance (set using
 * SetSpatialBandwidth()) and within a spectral range (set using
 * SetRangeBandwidth()). The resulting filtered image can be retrieved by
415
 * GetOutput() or GetRangeOutput(). Parameter SetRangeBandwidthRamp()
416
 * allows linearly adapting the range bandwidth to the intensity of
417
 * each channel if set greater than 0 (default value is 0).
418
 *
419 420 421 422 423
 * There are additional output images, as explained below.
 * Internally, the algorithm will iteratively update a pixel both in position
 * and spectral value, with respect to its neighbors, until convergence to a
 * local mode. The map of the distance traveled by pixels is obtained by
 * GetSpatialOutput(). A map of detected local modes is also available in
424 425
 * GetLabelOutput() if mode search is set to true (default set to false)
 * and can be seen as a first segmentation of the input image,
426 427 428
 * although usually highly oversegmented.
 * Finally, GetIterationOutput() will return the number of algorithm iterations
 * for each pixel.
429
 *
430 431 432 433 434 435 436 437 438 439
 * The class template parameter TKernel allows one to choose how pixels in the
 * spatial and spectral neighborhood of a given pixel participate in the
 * smoothed result. By default, a uniform kernel is used (KernelUniform), giving
 * an equal weight to all neighbor pixels. KernelGaussian can also be used,
 * although the computation time is significantly higher. The TKernel class
 * should define operator(), taking a squared norm as parameter and returning a
 * real value between 0 and 1. It should also define GetRadius(), converting the
 * spatial bandwidth parameter to the spatial radius defining how many pixels
 * are in the processing window local to a pixel.
 *
440
 * MeanShifVector squared norm is compared with Threshold (set using Get/Set accessor) to define pixel convergence (1e-3 by default).
Jonathan Guinet's avatar
Jonathan Guinet committed
441
 * MaxIterationNumber defines maximum iteration number for each pixel convergence (set using Get/Set accessor). Set to 4 by default.
442
 * ModeSearch is a boolean value, to choose between optimized and non optimized algorithm. If set to true (by default), assign mode value to each pixel on a path covered in convergence steps.
443 444 445 446 447 448 449 450 451 452
 *
 * For more information on mean shift techniques, one might consider reading the following article:
 *
 * D. Comaniciu, P. Meer, "Mean Shift: A Robust Approach Toward Feature Space Analysis," IEEE Transactions on
 * Pattern Analysis and Machine Intelligence, vol. 24, no. 5, pp. 603-619, May, 2002
 * D. Comaniciu, P. Meer, "Robust analysis of feature spaces: color image segmentation," cvpr, p. 750, 1997
 * IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR'97), 1997
 * D. Comaniciu, P. Meer, "Mean Shift Analysis and Applications," iccv, p. 1197, Seventh International Conference
 * on Computer Vision (ICCV'99) - Volume 2, 1999
 *
453
 * \sa MeanShiftSegmentationFilter
454
 *
Jonathan Guinet's avatar
Jonathan Guinet committed
455
 * \ingroup ImageSegmentation
Jonathan Guinet's avatar
Jonathan Guinet committed
456
 * \ingroup ImageEnhancement
457 458
 *
 * \ingroup OTBSmoothing
459
 */
460 461 462
template<class TInputImage, class TOutputImage, class TKernel = Meanshift::KernelUniform,
    class TOutputIterationImage = otb::Image<unsigned int, TInputImage::ImageDimension> >
class ITK_EXPORT MeanShiftSmoothingImageFilter: public itk::ImageToImageFilter<TInputImage, TOutputImage>
463 464 465
{
public:
  /** Standard class typedef */
466
  typedef MeanShiftSmoothingImageFilter Self;
467
  typedef itk::ImageToImageFilter<TInputImage, TOutputImage> Superclass;
468 469 470
  typedef itk::SmartPointer<Self> Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;
  typedef double RealType;
471

472
  /** Type macro */
473 474
  itkTypeMacro(MeanShiftSmoothingImageFilter, ImageToImageFilter);
; itkNewMacro(Self);
OTB Bot's avatar
STYLE  
OTB Bot committed
475
;
476 477

  /** Template parameters typedefs */
478

479 480 481 482 483
  typedef TInputImage InputImageType;
  typedef typename InputImageType::Pointer InputImagePointerType;
  typedef typename InputImageType::PixelType InputPixelType;
  typedef typename InputImageType::IndexType InputIndexType;
  typedef typename InputImageType::SizeType InputSizeType;
484
  typedef typename InputImageType::IndexValueType InputIndexValueType;
485 486 487
  typedef typename InputImageType::PointType PointType;
  typedef typename InputImageType::RegionType RegionType;
  typedef typename InputImageType::SizeType SizeType;
488

489 490 491 492
  typedef TOutputImage OutputImageType;
  typedef typename OutputImageType::Pointer OutputImagePointerType;
  typedef typename OutputImageType::PixelType OutputPixelType;
  typedef typename OutputImageType::RegionType OutputRegionType;
493

494
  typedef TOutputIterationImage OutputIterationImageType;
495

496
  typedef unsigned long LabelType;
497 498
  typedef otb::Image<LabelType, InputImageType::ImageDimension> OutputLabelImageType;

499
  typedef otb::VectorImage<RealType, InputImageType::ImageDimension> OutputSpatialImageType;
500 501
  typedef typename OutputSpatialImageType::Pointer OutputSpatialImagePointerType;
  typedef typename OutputSpatialImageType::PixelType OutputSpatialPixelType;
502

503
  typedef TKernel KernelType;
504

505 506
  itkStaticConstMacro(ImageDimension, unsigned int, InputImageType::ImageDimension);

507
  typedef itk::VariableLengthVector<RealType> RealVector;
508
  typedef otb::VectorImage<RealType, InputImageType::ImageDimension> RealVectorImageType;
509
  typedef otb::Image<unsigned short, InputImageType::ImageDimension> ModeTableImageType;
510

511
  /** Sets the spatial bandwidth (or radius in the case of a uniform kernel)
512 513
   * of the neighborhood for each pixel
   */
514 515
  itkSetMacro(SpatialBandwidth, RealType);
  itkGetConstReferenceMacro(SpatialBandwidth, RealType);
516 517

  /** Sets the spectral bandwidth (or radius for a uniform kernel) for pixels
518 519
   * to be included in the same mode
   */
520 521 522 523 524 525 526 527
  itkSetMacro(RangeBandwidth, RealType);
  itkGetConstReferenceMacro(RangeBandwidth, RealType);

   /** Sets the range bandwidth ramp. If > 0, the range bandwidth
    * will be y = RangeBandwidthRamp * x + RangeBandwidth, where x is
    * the band value. */
  itkSetMacro(RangeBandwidthRamp, RealType);
  itkGetConstReferenceMacro(RangeBandwidthRamp,RealType);
528

529
  /** Sets the maximum number of algorithm iterations */
530 531
  itkGetConstReferenceMacro(MaxIterationNumber, unsigned int);
  itkSetMacro(MaxIterationNumber, unsigned int);
532

533
  /** Sets the threshold value for the mean shift vector's squared norm,
534 535
   * under which convergence is assumed
   */
536 537
  itkGetConstReferenceMacro(Threshold, double);
  itkSetMacro(Threshold, double);
538

539
  /** Toggle mode search, which is enabled by default.
540 541 542
   * When off, the output label image is not available
   * Be careful, with this option, the result will slightly depend on thread number.
   */
543 544
  itkSetMacro(ModeSearch, bool);
  itkGetConstReferenceMacro(ModeSearch, bool);
545

Jonathan Guinet's avatar
Jonathan Guinet committed
546
#if 0
547
  /** Toggle bucket optimization, which is disabled by default.
548
   */
549 550
  itkSetMacro(BucketOptimization, bool);
  itkGetConstReferenceMacro(BucketOptimization, bool);
551
#endif
552

553
  /** Global shift allows tackling down numerical instabilities by
Julien Michel's avatar
Julien Michel committed
554
  aligning pixel indices when performing tile processing */
555 556
  itkSetMacro(GlobalShift,InputIndexType);

Jonathan Guinet's avatar
Jonathan Guinet committed
557
  /** Returns the const spatial image output,spatial image output is a displacement map (pixel position after convergence minus pixel index)  */
558
  const OutputSpatialImageType * GetSpatialOutput() const;
Jonathan Guinet's avatar
Jonathan Guinet committed
559
  /** Returns the const spectral image output */
560
  const OutputImageType * GetRangeOutput() const;
Jonathan Guinet's avatar
Jonathan Guinet committed
561
  /** Returns the const number of iterations map. */
562
  const OutputIterationImageType * GetIterationOutput() const;
Jonathan Guinet's avatar
Jonathan Guinet committed
563
  /** Returns the const image of region labels. This output does not have sense without mode search optimization (each label codes for one mode)*/
564 565
  const OutputLabelImageType * GetLabelOutput() const;

Jonathan Guinet's avatar
Jonathan Guinet committed
566
  /** Returns the spatial image output,spatial image output is a displacement map (pixel position after convergence minus pixel index)  */
567
  OutputSpatialImageType * GetSpatialOutput();
568
  /** Returns the spectral image output */
569
  OutputImageType * GetRangeOutput();
570 571
  /** Returns the number of iterations done at each pixel */
  OutputIterationImageType * GetIterationOutput();
Jonathan Guinet's avatar
Jonathan Guinet committed
572
  /** Returns the image of region labels. This output does not have sense without mode search optimization (each label codes for one mode) */
573
  OutputLabelImageType * GetLabelOutput();
574

575 576 577
protected:

  /** GenerateOutputInformation
578 579 580
   *  Define output pixel size
   *
   **/
581
  void GenerateOutputInformation(void) override;
582

583
  void GenerateInputRequestedRegion() override;
584

585
  void BeforeThreadedGenerateData() override;
586 587 588 589 590 591 592 593 594 595 596

  /** MeanShiftFilter can be implemented as a multithreaded filter.
   * Therefore, this implementation provides a ThreadedGenerateData()
   * routine which is called for each processing thread. The output
   * image data is allocated automatically by the superclass prior to
   * calling ThreadedGenerateData().  ThreadedGenerateData can only
   * write to the portion of the output image specified by the
   * parameter "outputRegionForThread"
   *
   * \sa ImageToImageFilter::ThreadedGenerateData(),
   *     ImageToImageFilter::GenerateData() */
597
  void ThreadedGenerateData(const OutputRegionType& outputRegionForThread, itk::ThreadIdType threadId) override;
598

599
  void AfterThreadedGenerateData() override;
600

601
  /** Allocates the outputs (need to be reimplemented since outputs have different type) */
602
  void AllocateOutputs() override;
603 604

  /** Constructor */
605
  MeanShiftSmoothingImageFilter();
606

607
  /** Destructor */
608
  ~MeanShiftSmoothingImageFilter() override;
609

610
  /** PrintSelf method */
611
  void PrintSelf(std::ostream& os, itk::Indent indent) const override;
612

613 614
  virtual void CalculateMeanShiftVector(const typename RealVectorImageType::Pointer inputImagePtr,
                                        const RealVector& jointPixel, const OutputRegionType& outputRegion,
615
                                        const RealVector& bandwidth,
616
                                        RealVector& meanShiftVector);
Jonathan Guinet's avatar
Jonathan Guinet committed
617
#if 0
618
  virtual void CalculateMeanShiftVectorBucket(const RealVector& jointPixel, RealVector& meanShiftVector);
619 620
#endif

621
private:
622 623
  MeanShiftSmoothingImageFilter(const Self &) = delete;
  void operator =(const Self&) = delete;
624

625
  /** Range bandwidth */
626
  RealType m_RangeBandwidth;
627

628 629 630
  /** Coefficient */
  RealType m_RangeBandwidthRamp;

631
  /** Spatial bandwidth */
632
  RealType m_SpatialBandwidth;
633

634
  /** Radius of pixel neighborhood, determined by the kernel from the spatial bandwidth  */
635
  InputSizeType m_SpatialRadius;
636

637
  /** Threshold on the squared norm of the mean shift vector to decide when to stop iterating **/
638
  double m_Threshold;
639

640
  /** Maximum number of iterations **/
641 642
  unsigned int m_MaxIterationNumber;

643
  /** Kernel object, implementing operator() which returns a weight between 0 and 1
644
   * depending on the squared norm given in parameter **/
645
  KernelType m_Kernel;
646

647
  /** Number of components per pixel in the input image */
648
  unsigned int m_NumberOfComponentsPerPixel;
649 650 651 652

  /** Input data in the joint spatial-range domain, scaled by the bandwidths */
  typename RealVectorImageType::Pointer m_JointImage;

653 654 655 656 657 658
  /** Image to store the status at each pixel:
   * 0 : no mode has been found yet
   * 1 : a mode has been assigned to this pixel
   * 2 : pixel is in the path of the currently processed pixel and a mode will
   *     be assigned to it
   */
659
  typename ModeTableImageType::Pointer m_ModeTable;
660

661 662
  /** Boolean to enable mode search  */
  bool m_ModeSearch;
663

Jonathan Guinet's avatar
Jonathan Guinet committed
664
#if 0
665 666
  /** Boolean to enable bucket optimization */
  bool m_BucketOptimization;
667
#endif
668 669 670 671

  /** Mode counters (local to each thread) */
  itk::VariableLengthVector<LabelType> m_NumLabels;
  /** Number of bits used to represent the threadId in the most significant bits
672
   of labels */
673
  unsigned int m_ThreadIdNumberOfBits;
674

Jonathan Guinet's avatar
Jonathan Guinet committed
675
#if 0
676
  typedef Meanshift::BucketImage<RealVectorImageType> BucketImageType;
677
  BucketImageType m_BucketImage;
678
#endif
679

680 681
  InputIndexType m_GlobalShift;

682
};
683

684 685 686
} // end namespace otb

#ifndef OTB_MANUAL_INSTANTIATION
687
#include "otbMeanShiftSmoothingImageFilter.hxx"
688 689 690
#endif

#endif