diff --git a/Code/BasicFilters/otbMeanShiftSmoothingImageFilter.h b/Code/BasicFilters/otbMeanShiftSmoothingImageFilter.h index a333247e5d043fe4505fac9768c90bc3254e3e53..8caa39d518c269f6a34abc1ae4efea26a8e0603c 100644 --- a/Code/BasicFilters/otbMeanShiftSmoothingImageFilter.h +++ b/Code/BasicFilters/otbMeanShiftSmoothingImageFilter.h @@ -1,20 +1,20 @@ /*========================================================================= - Program: ORFEO Toolbox - Language: C++ - Date: $Date$ - Version: $Revision$ + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ - Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. - See OTBCopyright.txt for details. + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. - This software is distributed WITHOUT ANY WARRANTY; without even - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the above copyright notices for more information. + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. -=========================================================================*/ + =========================================================================*/ #ifndef __otbMeanShiftSmoothingImageFilter_h #define __otbMeanShiftSmoothingImageFilter_h @@ -25,23 +25,22 @@ #include "itkImageRegionConstIteratorWithIndex.h" #include <vcl_algorithm.h> - +#define USEBUCKET 0 // don't use bucket optimization namespace otb { namespace Meanshift { -template <typename T> inline T simple_pow(T const& v, unsigned int p) +template<typename T> inline T simple_pow(T const& v, unsigned int p) { T res = 1; - for (unsigned int i=0; i!=p; ++i) + for (unsigned int i = 0; i != p; ++i) { res *= v; } return res; } - /** \class SpatialRangeJointDomainTransform * * @@ -55,33 +54,35 @@ class SpatialRangeJointDomainTransform public: typedef double RealType; - SpatialRangeJointDomainTransform() {} + SpatialRangeJointDomainTransform() + { + } // ~SpatialRangeJointDomainTransform() {} - typename TOutputJointImage::PixelType operator() - (const typename TInputImage::PixelType & inputPixel, const typename TInputImage::IndexType & index) const + typename TOutputJointImage::PixelType operator()(const typename TInputImage::PixelType & inputPixel, + const typename TInputImage::IndexType & index) const { typename TOutputJointImage::PixelType jointPixel(m_ImageDimension + m_NumberOfComponentsPerPixel); - for(unsigned int comp = 0; comp < m_ImageDimension; comp++) + for (unsigned int comp = 0; comp < m_ImageDimension; comp++) { jointPixel[comp] = index[comp] / m_SpatialBandwidth; } - for(unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++) + for (unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++) { jointPixel[m_ImageDimension + comp] = inputPixel[comp] / m_RangeBandwidth; } return jointPixel; } - void Initialize(unsigned int _ImageDimension, unsigned int numberOfComponentsPerPixel_, - RealType spatialBandwidth_, RealType rangeBandwidth_) + void Initialize(unsigned int _ImageDimension, unsigned int numberOfComponentsPerPixel_, RealType spatialBandwidth_, + RealType rangeBandwidth_) { - m_ImageDimension = _ImageDimension; + m_ImageDimension = _ImageDimension; m_NumberOfComponentsPerPixel = numberOfComponentsPerPixel_; - m_SpatialBandwidth = spatialBandwidth_; - m_RangeBandwidth = rangeBandwidth_; - m_OutputSize = m_ImageDimension + m_NumberOfComponentsPerPixel; + m_SpatialBandwidth = spatialBandwidth_; + m_RangeBandwidth = rangeBandwidth_; + m_OutputSize = m_ImageDimension + m_NumberOfComponentsPerPixel; } unsigned int GetOutputSize() const @@ -93,11 +94,10 @@ private: unsigned int m_ImageDimension; unsigned int m_NumberOfComponentsPerPixel; unsigned int m_OutputSize; - RealType m_SpatialBandwidth; - RealType m_RangeBandwidth; + RealType m_SpatialBandwidth; + RealType m_RangeBandwidth; }; - class KernelUniform { public: @@ -106,11 +106,13 @@ public: // KernelUniform() {} // ~KernelUniform() {} - RealType operator() (RealType x) const{ + RealType operator()(RealType x) const + { return (x <= 1) ? 1.0 : 0.0; } - RealType GetRadius(RealType bandwidth) const { + RealType GetRadius(RealType bandwidth) const + { return bandwidth; } }; @@ -123,12 +125,14 @@ public: // KernelGaussian() {} // ~KernelGaussian() {} - RealType operator() (RealType x) const { - return vcl_exp(-0.5*x); + RealType operator()(RealType x) const + { + return vcl_exp(-0.5 * x); } - RealType GetRadius(RealType bandwidth) const { - return 3.0*bandwidth; + RealType GetRadius(RealType bandwidth) const + { + return 3.0 * bandwidth; } }; @@ -138,25 +142,28 @@ public: * access to pixels in vector images through the method GetPixelPointer */ template<typename TImage> -class FastImageRegionConstIterator - : public itk::ImageRegionConstIterator<TImage> +class FastImageRegionConstIterator: public itk::ImageRegionConstIterator<TImage> { public: /** Standard class typedef. */ - typedef FastImageRegionConstIterator<TImage> Self; - typedef itk::ImageRegionConstIterator<TImage> Superclass; + typedef FastImageRegionConstIterator<TImage> Self; + typedef itk::ImageRegionConstIterator<TImage> Superclass; - typedef typename Superclass::ImageType ImageType; - typedef typename Superclass::RegionType RegionType; + typedef typename Superclass::ImageType ImageType; + typedef typename Superclass::RegionType RegionType; - typedef typename TImage::PixelType PixelType; + typedef typename TImage::PixelType PixelType; typedef typename TImage::InternalPixelType InternalPixelType; - itkTypeMacro(FastImageRegionConstIterator, ImageRegionConstIterator); + itkTypeMacro(FastImageRegionConstIterator, ImageRegionConstIterator) + ; - FastImageRegionConstIterator() : Superclass() {} - FastImageRegionConstIterator(const ImageType *ptr, const RegionType ®ion) - : Superclass(ptr, region) + FastImageRegionConstIterator() : + Superclass() + { + } + FastImageRegionConstIterator(const ImageType *ptr, const RegionType ®ion) : + Superclass(ptr, region) { m_NumberOfComponentsPerPixel = ptr->GetNumberOfComponentsPerPixel(); } @@ -170,58 +177,58 @@ private: unsigned int m_NumberOfComponentsPerPixel; }; - +#ifdef USEBUCKET //disable bucket mode /** \class BucketImage - * - * 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(). - */ -template <class TImage> + * + * 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(). + */ +template<class TImage> class BucketImage { public: - 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; + 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; - typedef double RealType; + typedef double RealType; static const unsigned int ImageDimension = ImageType::ImageDimension; /** The bucket image has dimension N+1 (ie. usually 3D for most images) */ - typedef std::vector<typename ImageType::SizeType::SizeValueType> BucketImageSizeType; + typedef std::vector<typename ImageType::SizeType::SizeValueType> BucketImageSizeType; typedef std::vector<typename ImageType::IndexType::IndexValueType> BucketImageIndexType; //typedef std::vector<long> BucketImageIndexType; /** pixel buckets typedefs and declarations */ - typedef const typename ImageType::InternalPixelType * ImageDataPointerType; - typedef std::vector<ImageDataPointerType> BucketType; - typedef std::vector<BucketType> BucketListType; + typedef const typename ImageType::InternalPixelType * ImageDataPointerType; + typedef std::vector<ImageDataPointerType> BucketType; + typedef std::vector<BucketType> BucketListType; - BucketImage() {} + BucketImage() + { + } /** Constructor for the bucket image. It operates on the specified - * 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 ) + * 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) { // Find max and min of the used spectral band @@ -230,7 +237,7 @@ public: InternalPixelType minValue = inputIt.Get()[spectralCoordinate]; InternalPixelType maxValue = minValue; ++inputIt; - while( !inputIt.IsAtEnd() ) + while (!inputIt.IsAtEnd()) { const PixelType &p = inputIt.Get(); minValue = vcl_min(minValue, p[m_SpectralCoordinate]); @@ -243,15 +250,15 @@ public: // Compute bucket image dimensions. Note: empty buckets are at each border // to simplify image border issues - m_DimensionVector.resize(ImageDimension+1); // NB: pays for a 0-innit - for(unsigned int dim = 0; dim < ImageDimension; ++dim) + m_DimensionVector.resize(ImageDimension + 1); // NB: pays for a 0-innit + for (unsigned int dim = 0; dim < ImageDimension; ++dim) { m_DimensionVector[dim] = m_Region.GetSize()[dim] / m_SpatialRadius + 3; } - m_DimensionVector[ImageDimension] = (unsigned int)((maxValue - minValue) / m_RangeRadius) + 3; + m_DimensionVector[ImageDimension] = (unsigned int) ((maxValue - minValue) / m_RangeRadius) + 3; unsigned int numBuckets = m_DimensionVector[0]; - for(unsigned int dim = 1; dim <= ImageDimension; ++dim) + for (unsigned int dim = 1; dim <= ImageDimension; ++dim) numBuckets *= m_DimensionVector[dim]; m_BucketList.resize(numBuckets); @@ -262,7 +269,7 @@ public: FastImageRegionConstIterator<ImageType> fastIt(m_Image, m_Region); fastIt.GoToBegin(); - while( !it.IsAtEnd() ) + while (!it.IsAtEnd()) { const IndexType & index = it.GetIndex(); const PixelType & pixel = it.Get(); @@ -280,15 +287,15 @@ public: // Prepare neighborhood offset vector // BucketImageIndexType zeroOffsetIndex(ImageDimension+1); std::vector<BucketImageIndexType> neighborsIndexList; - neighborsIndexList.reserve(simple_pow(3, ImageDimension+1)); - neighborsIndexList.resize(1, BucketImageIndexType(ImageDimension+1)); // zeroOffsetIndex + neighborsIndexList.reserve(simple_pow(3, ImageDimension + 1)); + neighborsIndexList.resize(1, BucketImageIndexType(ImageDimension + 1)); // zeroOffsetIndex // neighborsIndexList.push_back(zeroOffsetIndex); - for(unsigned dim = 0; dim <= ImageDimension; ++dim) + for (unsigned dim = 0; dim <= ImageDimension; ++dim) { // take all neighbors already in the list and add their direct neighbor // along the current dim const unsigned int curSize = neighborsIndexList.size(); - for(unsigned int i = 0; i < curSize; ++i) + for (unsigned int i = 0; i < curSize; ++i) { BucketImageIndexType index = neighborsIndexList[i]; index[dim]--; @@ -300,33 +307,35 @@ public: // Convert all neighbors n-dimensional indices to bucket list 1D indices const unsigned int neighborhoodOffsetVectorSize = neighborsIndexList.size(); m_NeighborhoodOffsetVector.reserve(neighborhoodOffsetVectorSize); - for(unsigned int i = 0; i < neighborhoodOffsetVectorSize; ++i) + for (unsigned int i = 0; i < neighborhoodOffsetVectorSize; ++i) { const int listIndex = BucketIndexToBucketListIndex(neighborsIndexList[i]); m_NeighborhoodOffsetVector.push_back(listIndex); } } - ~BucketImage() {} + ~BucketImage() + { + } /** Returns the N+1-dimensional bucket index for a given pixel value at the given index */ BucketImageIndexType GetBucketIndex(const PixelType & pixel, const IndexType & index) { - BucketImageIndexType bucketIndex(ImageDimension+1); - for(unsigned int dim = 0; dim < ImageDimension; ++dim) + BucketImageIndexType bucketIndex(ImageDimension + 1); + for (unsigned int dim = 0; dim < ImageDimension; ++dim) { bucketIndex[dim] = (index[dim] - m_Region.GetIndex()[dim]) / m_SpatialRadius + 1; - } + } bucketIndex[ImageDimension] = (pixel[m_SpectralCoordinate] - m_MinValue) / m_RangeRadius + 1; return bucketIndex; } /** Converts a N+1-dimensional bucket index into the 1D list index useable by - GetBucket() */ + GetBucket() */ int BucketIndexToBucketListIndex(const BucketImageIndexType & bucketIndex) const { int bucketListIndex = bucketIndex[0]; - for(unsigned int dim = 1; dim <= ImageDimension; ++dim) + for (unsigned int dim = 1; dim <= ImageDimension; ++dim) { bucketListIndex = bucketListIndex * m_DimensionVector[dim] + bucketIndex[dim]; } @@ -335,16 +344,16 @@ public: /** Retrieves the list of all buckets in the neighborhood of the given bucket */ std::vector<unsigned int> GetNeighborhoodBucketListIndices(int bucketIndex) const - { + { const unsigned int neighborhoodOffsetVectorSize = m_NeighborhoodOffsetVector.size(); std::vector<unsigned int> indices(neighborhoodOffsetVectorSize); - for(unsigned int i = 0; i < neighborhoodOffsetVectorSize; ++i) + for (unsigned int i = 0; i < neighborhoodOffsetVectorSize; ++i) { indices[i] = bucketIndex + m_NeighborhoodOffsetVector[i]; } return indices; - } + } /* Returns the list of pixels (actually pointer to pixel data) contained in a bucket */ const BucketType & GetBucket(unsigned int index) const @@ -368,7 +377,7 @@ private: /** 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 - also their value at one coordinate */ + also their value at one coordinate */ unsigned int m_SpectralCoordinate; /** Min and Max of selected spectral coordinate */ InternalPixelType m_MinValue; @@ -379,10 +388,11 @@ private: /** This vector holds the dimensions of the 3D (ND?) bucket image */ BucketImageSizeType m_DimensionVector; /** Vector of offsets in the buckets list to get all buckets in the - * neighborhood - */ + * neighborhood + */ std::vector<int> m_NeighborhoodOffsetVector; }; +#endif } // end namespace Meanshift @@ -434,89 +444,98 @@ private: * \ingroup ImageSegmentation * \ingroup ImageEnhancement */ -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> +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> { public: /** Standard class typedef */ - typedef MeanShiftSmoothingImageFilter Self; + typedef MeanShiftSmoothingImageFilter Self; typedef itk::ImageToImageFilter<TInputImage, TOutputImage> Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; - typedef double RealType; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + typedef double RealType; /** Type macro */ - itkTypeMacro(MeanShiftSmoothingImageFilter, ImageToImageFilter); - itkNewMacro(Self); + itkTypeMacro(MeanShiftSmoothingImageFilter, ImageToImageFilter) + ;itkNewMacro(Self) + ; /** Template parameters typedefs */ - typedef TInputImage InputImageType; - typedef typename InputImageType::Pointer InputImagePointerType; - typedef typename InputImageType::PixelType InputPixelType; - typedef typename InputImageType::IndexType InputIndexType; - typedef typename InputImageType::SizeType InputSizeType; + typedef TInputImage InputImageType; + typedef typename InputImageType::Pointer InputImagePointerType; + typedef typename InputImageType::PixelType InputPixelType; + typedef typename InputImageType::IndexType InputIndexType; + typedef typename InputImageType::SizeType InputSizeType; typedef typename InputImageType::IndexValueType InputIndexValueType; - typedef typename InputImageType::PointType PointType; - typedef typename InputImageType::RegionType RegionType; - typedef typename InputImageType::SizeType SizeType; + typedef typename InputImageType::PointType PointType; + typedef typename InputImageType::RegionType RegionType; + typedef typename InputImageType::SizeType SizeType; - typedef TOutputImage OutputImageType; - typedef typename OutputImageType::Pointer OutputImagePointerType; - typedef typename OutputImageType::PixelType OutputPixelType; - typedef typename OutputImageType::RegionType OutputRegionType; + typedef TOutputImage OutputImageType; + typedef typename OutputImageType::Pointer OutputImagePointerType; + typedef typename OutputImageType::PixelType OutputPixelType; + typedef typename OutputImageType::RegionType OutputRegionType; - typedef TOutputIterationImage OutputIterationImageType; + typedef TOutputIterationImage OutputIterationImageType; - typedef unsigned long LabelType; + typedef unsigned long LabelType; typedef otb::Image<LabelType, InputImageType::ImageDimension> OutputLabelImageType; typedef otb::VectorImage<RealType, InputImageType::ImageDimension> OutputSpatialImageType; - typedef typename OutputSpatialImageType::Pointer OutputSpatialImagePointerType; - typedef typename OutputSpatialImageType::PixelType OutputSpatialPixelType; + typedef typename OutputSpatialImageType::Pointer OutputSpatialImagePointerType; + typedef typename OutputSpatialImageType::PixelType OutputSpatialPixelType; typedef TKernel KernelType; itkStaticConstMacro(ImageDimension, unsigned int, InputImageType::ImageDimension); - typedef itk::VariableLengthVector<RealType> RealVector; + typedef itk::VariableLengthVector<RealType> RealVector; typedef otb::VectorImage<RealType, InputImageType::ImageDimension> RealVectorImageType; typedef otb::Image<unsigned short, InputImageType::ImageDimension> ModeTableImageType; /** Sets the spatial bandwidth (or radius in the case of a uniform kernel) - * of the neighborhood for each pixel - */ - itkSetMacro(SpatialBandwidth, RealType); - itkGetConstReferenceMacro(SpatialBandwidth, RealType); + * of the neighborhood for each pixel + */ + itkSetMacro(SpatialBandwidth, RealType) + ;itkGetConstReferenceMacro(SpatialBandwidth, RealType) + ; /** Sets the spectral bandwidth (or radius for a uniform kernel) for pixels - * to be included in the same mode - */ - itkSetMacro(RangeBandwidth, RealType); - itkGetConstReferenceMacro(RangeBandwidth, RealType); + * to be included in the same mode + */ + itkSetMacro(RangeBandwidth, RealType) + ;itkGetConstReferenceMacro(RangeBandwidth, RealType) + ; /** Sets the maximum number of algorithm iterations */ - itkGetConstReferenceMacro(MaxIterationNumber, unsigned int); - itkSetMacro(MaxIterationNumber, unsigned int); + itkGetConstReferenceMacro(MaxIterationNumber, unsigned int) + ;itkSetMacro(MaxIterationNumber, unsigned int) + ; /** Sets the threshold value for the mean shift vector's squared norm, - * under which convergence is assumed - */ - itkGetConstReferenceMacro(Threshold, double); - itkSetMacro(Threshold, double); + * under which convergence is assumed + */ + itkGetConstReferenceMacro(Threshold, double) + ;itkSetMacro(Threshold, double) + ; /** Toggle mode search, which is enabled by default. - * When off, the output label image is not available - * Be careful, with this option, the result will slightly depend on thread number. - */ - itkSetMacro(ModeSearch, bool); - itkGetConstReferenceMacro(ModeSearch, bool); + * When off, the output label image is not available + * Be careful, with this option, the result will slightly depend on thread number. + */ + itkSetMacro(ModeSearch, bool) + ;itkGetConstReferenceMacro(ModeSearch, bool) + ; +#ifdef USEBUCKET /** Toggle bucket optimization, which is disabled by default. - */ - itkSetMacro(BucketOptimization, bool); - itkGetConstReferenceMacro(BucketOptimization, bool); + */ + itkSetMacro(BucketOptimization, bool) + ;itkGetConstReferenceMacro(BucketOptimization, bool) + ; +#endif /** Returns the const spatial image output,spatial image output is a displacement map (pixel position after convergence minus pixel index) */ const OutputSpatialImageType * GetSpatialOutput() const; @@ -529,7 +548,7 @@ public: /** Returns the spatial image output,spatial image output is a displacement map (pixel position after convergence minus pixel index) */ OutputSpatialImageType * GetSpatialOutput(); - /** Returns the spectral image output */ + /** Returns the spectral image output */ OutputImageType * GetRangeOutput(); /** Returns the number of iterations done at each pixel */ OutputIterationImageType * GetIterationOutput(); @@ -538,29 +557,27 @@ public: protected: - /** GenerateOutputInformation - * Define output pixel size - * - **/ - virtual void GenerateOutputInformation(void); - - virtual void GenerateInputRequestedRegion(); - - virtual void BeforeThreadedGenerateData(); - - /** 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() */ - void ThreadedGenerateData(const OutputRegionType& outputRegionForThread, - int threadId ); + * Define output pixel size + * + **/ + virtual void GenerateOutputInformation(void); + + virtual void GenerateInputRequestedRegion(); + + virtual void BeforeThreadedGenerateData(); + + /** 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() */ + void ThreadedGenerateData(const OutputRegionType& outputRegionForThread, int threadId); virtual void AfterThreadedGenerateData(); @@ -579,19 +596,22 @@ protected: virtual void CalculateMeanShiftVector(const typename RealVectorImageType::Pointer inputImagePtr, const RealVector& jointPixel, const OutputRegionType& outputRegion, RealVector& meanShiftVector); +#ifdef USEBUCKET virtual void CalculateMeanShiftVectorBucket(const RealVector& jointPixel, RealVector& meanShiftVector); +#endif + private: MeanShiftSmoothingImageFilter(const Self &); //purposely not implemented - void operator =(const Self&); //purposely not implemented + void operator =(const Self&); //purposely not implemented /** Range bandwidth */ - RealType m_RangeBandwidth; + RealType m_RangeBandwidth; /** Spatial bandwidth */ - RealType m_SpatialBandwidth; + RealType m_SpatialBandwidth; /** Radius of pixel neighborhood, determined by the kernel from the spatial bandwidth */ - InputSizeType m_SpatialRadius; + InputSizeType m_SpatialRadius; /** Threshold on the squared norm of the mean shift vector to decide when to stop iterating **/ double m_Threshold; @@ -600,7 +620,7 @@ private: unsigned int m_MaxIterationNumber; /** Kernel object, implementing operator() which returns a weight between 0 and 1 -* depending on the squared norm given in parameter **/ + * depending on the squared norm given in parameter **/ KernelType m_Kernel; /** Number of components per pixel in the input image */ @@ -619,17 +639,22 @@ private: /** Boolean to enable mode search */ bool m_ModeSearch; + +#ifdef USEBUCKET /** Boolean to enable bucket optimization */ bool m_BucketOptimization; +#endif /** Mode counters (local to each thread) */ itk::VariableLengthVector<LabelType> m_NumLabels; /** Number of bits used to represent the threadId in the most significant bits - of labels */ + of labels */ unsigned int m_ThreadIdNumberOfBits; +#ifdef USEBUCKET typedef Meanshift::BucketImage<RealVectorImageType> BucketImageType; BucketImageType m_BucketImage; +#endif }; diff --git a/Code/BasicFilters/otbMeanShiftSmoothingImageFilter.txx b/Code/BasicFilters/otbMeanShiftSmoothingImageFilter.txx index 5dbbbb48342d0ee2557b20974c6aeb46c12cdd7c..32406f67901f25b8582f2a87f762b79f437f5ba4 100644 --- a/Code/BasicFilters/otbMeanShiftSmoothingImageFilter.txx +++ b/Code/BasicFilters/otbMeanShiftSmoothingImageFilter.txx @@ -30,20 +30,19 @@ namespace otb { -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::MeanShiftSmoothingImageFilter() -: m_RangeBandwidth( 16. ) - , m_SpatialBandwidth( 3 ) +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::MeanShiftSmoothingImageFilter() : + m_RangeBandwidth(16.), m_SpatialBandwidth(3) // , m_SpatialRadius(???) - , m_Threshold( 1e-3 ) - , m_MaxIterationNumber( 10 ) - // , m_Kernel(...) - // , m_NumberOfComponentsPerPixel(...) - // , m_JointImage(0) - // , m_ModeTable(0) - , m_ModeSearch( true ) - , m_BucketOptimization( false ) + , m_Threshold(1e-3), m_MaxIterationNumber(10) + // , m_Kernel(...) + // , m_NumberOfComponentsPerPixel(...) + // , m_JointImage(0) + // , m_ModeTable(0) + , m_ModeSearch(true) +#ifdef USEBUCKET + , m_BucketOptimization(false) +#endif { this->SetNumberOfOutputs(4); this->SetNthOutput(0, OutputImageType::New()); @@ -52,89 +51,74 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati this->SetNthOutput(3, OutputLabelImageType::New()); } - -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::~ MeanShiftSmoothingImageFilter() +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::~MeanShiftSmoothingImageFilter() { } -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> const typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputSpatialImageType * -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::GetSpatialOutput() const +MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetSpatialOutput() const { - return static_cast<const OutputSpatialImageType *>(this->itk::ProcessObject::GetOutput(1)); + return static_cast<const OutputSpatialImageType *> (this->itk::ProcessObject::GetOutput(1)); } -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputSpatialImageType * -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::GetSpatialOutput() +MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetSpatialOutput() { - return static_cast<OutputSpatialImageType *>(this->itk::ProcessObject::GetOutput(1)); + return static_cast<OutputSpatialImageType *> (this->itk::ProcessObject::GetOutput(1)); } - -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> const typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputImageType * -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::GetRangeOutput() const +MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetRangeOutput() const { - return static_cast<const OutputImageType *>(this->itk::ProcessObject::GetOutput(0)); + return static_cast<const OutputImageType *> (this->itk::ProcessObject::GetOutput(0)); } -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputImageType * -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::GetRangeOutput() +MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetRangeOutput() { - return static_cast<OutputImageType *>(this->itk::ProcessObject::GetOutput(0)); + return static_cast<OutputImageType *> (this->itk::ProcessObject::GetOutput(0)); } - -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputIterationImageType * -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::GetIterationOutput() +MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetIterationOutput() { - return static_cast<OutputIterationImageType *>(this->itk::ProcessObject::GetOutput(2)); + return static_cast<OutputIterationImageType *> (this->itk::ProcessObject::GetOutput(2)); } -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> const typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputIterationImageType * -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::GetIterationOutput() const +MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetIterationOutput() const { - return static_cast<OutputIterationImageType *>(this->itk::ProcessObject::GetOutput(2)); + return static_cast<OutputIterationImageType *> (this->itk::ProcessObject::GetOutput(2)); } -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputLabelImageType * -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::GetLabelOutput() +MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetLabelOutput() { - return static_cast<OutputLabelImageType *>(this->itk::ProcessObject::GetOutput(3)); + return static_cast<OutputLabelImageType *> (this->itk::ProcessObject::GetOutput(3)); } -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> const typename MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::OutputLabelImageType * -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::GetLabelOutput() const +MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GetLabelOutput() const { - return static_cast<OutputLabelImageType *>(this->itk::ProcessObject::GetOutput(3)); + return static_cast<OutputLabelImageType *> (this->itk::ProcessObject::GetOutput(3)); } - -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> -void -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::AllocateOutputs() +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::AllocateOutputs() { - typename OutputSpatialImageType::Pointer spatialOutputPtr = this->GetSpatialOutput(); - typename OutputImageType::Pointer rangeOutputPtr = this->GetRangeOutput(); + typename OutputSpatialImageType::Pointer spatialOutputPtr = this->GetSpatialOutput(); + typename OutputImageType::Pointer rangeOutputPtr = this->GetRangeOutput(); typename OutputIterationImageType::Pointer iterationOutputPtr = this->GetIterationOutput(); - typename OutputLabelImageType::Pointer labelOutputPtr = this->GetLabelOutput(); + typename OutputLabelImageType::Pointer labelOutputPtr = this->GetLabelOutput(); spatialOutputPtr->SetBufferedRegion(spatialOutputPtr->GetRequestedRegion()); spatialOutputPtr->Allocate(); @@ -147,13 +131,10 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati labelOutputPtr->SetBufferedRegion(labelOutputPtr->GetRequestedRegion()); labelOutputPtr->Allocate(); - } - +} -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> -void -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::GenerateOutputInformation() +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GenerateOutputInformation() { Superclass::GenerateOutputInformation(); @@ -169,32 +150,28 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati } } - -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> -void -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::GenerateInputRequestedRegion() +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::GenerateInputRequestedRegion() { // Call superclass implementation Superclass::GenerateInputRequestedRegion(); // Retrieve input pointers - InputImagePointerType inPtr = const_cast<TInputImage *>(this->GetInput()); + InputImagePointerType inPtr = const_cast<TInputImage *> (this->GetInput()); OutputImagePointerType outRangePtr = this->GetRangeOutput(); // Check pointers before using them - if ( !inPtr || !outRangePtr ) + if (!inPtr || !outRangePtr) { return; } - // Retrieve requested region (TODO: check if we need to handle // region for outHDispPtr) RegionType outputRequestedRegion = outRangePtr->GetRequestedRegion(); // Pad by the appropriate radius - RegionType inputRequestedRegion = outputRequestedRegion; + RegionType inputRequestedRegion = outputRequestedRegion; // Initializes the spatial radius from kernel bandwidth m_SpatialRadius.Fill(m_Kernel.GetRadius(m_SpatialBandwidth)); @@ -202,43 +179,40 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati inputRequestedRegion.PadByRadius(m_SpatialRadius); // Crop the input requested region at the input's largest possible region - if ( inputRequestedRegion.Crop(inPtr->GetLargestPossibleRegion()) ) - { - inPtr->SetRequestedRegion( inputRequestedRegion ); - return; - } - else - { - // Couldn't crop the region (requested region is outside the largest - // possible region). Throw an exception. - - // store what we tried to request (prior to trying to crop) - inPtr->SetRequestedRegion( inputRequestedRegion ); - - // build an exception - itk::InvalidRequestedRegionError e(__FILE__, __LINE__); - e.SetLocation(ITK_LOCATION); - e.SetDescription("Requested region is (at least partially) outside the largest possible region."); - e.SetDataObject(inPtr); - throw e; - } + if (inputRequestedRegion.Crop(inPtr->GetLargestPossibleRegion())) + { + inPtr->SetRequestedRegion(inputRequestedRegion); + return; + } + else + { + // Couldn't crop the region (requested region is outside the largest + // possible region). Throw an exception. + + // store what we tried to request (prior to trying to crop) + inPtr->SetRequestedRegion(inputRequestedRegion); + + // build an exception + itk::InvalidRequestedRegionError e(__FILE__, __LINE__); + e.SetLocation(ITK_LOCATION); + e.SetDescription("Requested region is (at least partially) outside the largest possible region."); + e.SetDataObject(inPtr); + throw e; + } } - -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> -void -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::BeforeThreadedGenerateData() +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::BeforeThreadedGenerateData() { // typedef itk::ImageRegionConstIteratorWithIndex<InputImageType> InputIteratorWithIndexType; // typedef itk::ImageRegionIterator<RealVectorImageType> JointImageIteratorType; - OutputSpatialImagePointerType outSpatialPtr = this->GetSpatialOutput(); - OutputImagePointerType outRangePtr = this->GetRangeOutput(); - typename InputImageType::ConstPointer inputPtr = this->GetInput(); + OutputSpatialImagePointerType outSpatialPtr = this->GetSpatialOutput(); + OutputImagePointerType outRangePtr = this->GetRangeOutput(); + typename InputImageType::ConstPointer inputPtr = this->GetInput(); typename OutputIterationImageType::Pointer iterationOutput = this->GetIterationOutput(); - typename OutputSpatialImageType::Pointer spatialOutput = this->GetSpatialOutput(); + typename OutputSpatialImageType::Pointer spatialOutput = this->GetSpatialOutput(); //InputIndexType index; @@ -247,7 +221,6 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati m_NumberOfComponentsPerPixel = this->GetInput()->GetNumberOfComponentsPerPixel(); - // Allocate output images this->AllocateOutputs(); @@ -257,64 +230,67 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati zero.Fill(0); spatialOutput->FillBuffer(zero); - // m_JointImage is the input data expressed in the joint spatial-range // domain, i.e. spatial coordinates are concatenated to the range values. // Moreover, pixel components in this image are normalized by their respective // (spatial or range) bandwith. typedef Meanshift::SpatialRangeJointDomainTransform<InputImageType, RealVectorImageType> FunctionType; - typedef otb::UnaryFunctorWithIndexWithOutputSizeImageFilter<InputImageType, RealVectorImageType, FunctionType> JointImageFunctorType; + typedef otb::UnaryFunctorWithIndexWithOutputSizeImageFilter<InputImageType, RealVectorImageType, FunctionType> + JointImageFunctorType; typename JointImageFunctorType::Pointer jointImageFunctor = JointImageFunctorType::New(); jointImageFunctor->SetInput(inputPtr); - jointImageFunctor->GetFunctor().Initialize(ImageDimension, m_NumberOfComponentsPerPixel, m_SpatialBandwidth, m_RangeBandwidth); + jointImageFunctor->GetFunctor().Initialize(ImageDimension, m_NumberOfComponentsPerPixel, m_SpatialBandwidth, + m_RangeBandwidth); jointImageFunctor->Update(); m_JointImage = jointImageFunctor->GetOutput(); - if(m_BucketOptimization) +#ifdef USEBUCKET + if (m_BucketOptimization) { - // Create bucket image - // Note: because values in the input m_JointImage are normalized, the - // rangeRadius argument is just 1 - m_BucketImage = BucketImageType(static_cast<typename RealVectorImageType::ConstPointer>(m_JointImage), - m_JointImage->GetRequestedRegion(), - m_Kernel.GetRadius(m_SpatialBandwidth), 1, ImageDimension); + // Create bucket image + // Note: because values in the input m_JointImage are normalized, the + // rangeRadius argument is just 1 + m_BucketImage = BucketImageType(static_cast<typename RealVectorImageType::ConstPointer> (m_JointImage), + m_JointImage->GetRequestedRegion(), m_Kernel.GetRadius(m_SpatialBandwidth), 1, + ImageDimension); } -/* - // Allocate the joint domain image - m_JointImage = RealVectorImageType::New(); - m_JointImage->SetNumberOfComponentsPerPixel(ImageDimension + m_NumberOfComponentsPerPixel); - m_JointImage->SetRegions(inputPtr->GetRequestedRegion()); - m_JointImage->Allocate(); - - InputIteratorWithIndexType inputIt(inputPtr, inputPtr->GetRequestedRegion()); - JointImageIteratorType jointIt(m_JointImage, inputPtr->GetRequestedRegion()); - - // Initialize the joint image with scaled values - inputIt.GoToBegin(); - jointIt.GoToBegin(); - - while (!inputIt.IsAtEnd()) - { - typename InputImageType::PixelType const& inputPixel = inputIt.Get(); - index = inputIt.GetIndex(); - - RealVector & jointPixel = jointIt.Get(); - for(unsigned int comp = 0; comp < ImageDimension; comp++) - { - jointPixel[comp] = index[comp] / m_SpatialBandwidth; - } - for(unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++) - { - jointPixel[ImageDimension + comp] = inputPixel[comp] / m_RangeBandwidth; - } - // jointIt.Set(jointPixel); - - ++inputIt; - ++jointIt; - } -*/ +#endif + /* + // Allocate the joint domain image + m_JointImage = RealVectorImageType::New(); + m_JointImage->SetNumberOfComponentsPerPixel(ImageDimension + m_NumberOfComponentsPerPixel); + m_JointImage->SetRegions(inputPtr->GetRequestedRegion()); + m_JointImage->Allocate(); + + InputIteratorWithIndexType inputIt(inputPtr, inputPtr->GetRequestedRegion()); + JointImageIteratorType jointIt(m_JointImage, inputPtr->GetRequestedRegion()); + + // Initialize the joint image with scaled values + inputIt.GoToBegin(); + jointIt.GoToBegin(); + + while (!inputIt.IsAtEnd()) + { + typename InputImageType::PixelType const& inputPixel = inputIt.Get(); + index = inputIt.GetIndex(); + + RealVector & jointPixel = jointIt.Get(); + for(unsigned int comp = 0; comp < ImageDimension; comp++) + { + jointPixel[comp] = index[comp] / m_SpatialBandwidth; + } + for(unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++) + { + jointPixel[ImageDimension + comp] = inputPixel[comp] / m_RangeBandwidth; + } + // jointIt.Set(jointPixel); + + ++inputIt; + ++jointIt; + } + */ //TODO don't create mode table iterator when ModeSearch is set to false m_ModeTable = ModeTableImageType::New(); @@ -343,48 +319,47 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati n >>= 1; m_ThreadIdNumberOfBits++; } - if(m_ThreadIdNumberOfBits == 0) m_ThreadIdNumberOfBits = 1; // minimum 1 bit + if (m_ThreadIdNumberOfBits == 0) m_ThreadIdNumberOfBits = 1; // minimum 1 bit m_NumLabels.SetSize(numThreads); - for(unsigned int i = 0; i < numThreads; i++) + for (unsigned int i = 0; i < numThreads; i++) { - m_NumLabels[i] = static_cast<LabelType>(i) << (sizeof(LabelType)*8 - m_ThreadIdNumberOfBits); + m_NumLabels[i] = static_cast<LabelType> (i) << (sizeof(LabelType) * 8 - m_ThreadIdNumberOfBits); } } } - // Calculates the mean shift vector at the position given by jointPixel -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> -void -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::CalculateMeanShiftVector( - const typename RealVectorImageType::Pointer jointImage, const RealVector& jointPixel, - const OutputRegionType& outputRegion, RealVector& meanShiftVector) - { +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::CalculateMeanShiftVector( + const typename RealVectorImageType::Pointer jointImage, + const RealVector& jointPixel, + const OutputRegionType& outputRegion, + RealVector& meanShiftVector) +{ const unsigned int jointDimension = ImageDimension + m_NumberOfComponentsPerPixel; InputIndexType inputIndex; InputIndexType regionIndex; - InputSizeType regionSize; + InputSizeType regionSize; assert(meanShiftVector.GetSize() == jointDimension); meanShiftVector.Fill(0); // Calculates current pixel neighborhood region, restricted to the output image region - for(unsigned int comp = 0; comp < ImageDimension; ++comp) + for (unsigned int comp = 0; comp < ImageDimension; ++comp) { inputIndex[comp] = jointPixel[comp] * m_SpatialBandwidth; - regionIndex[comp] = vcl_max( - static_cast<long int>(outputRegion.GetIndex().GetElement(comp)), - static_cast<long int>(inputIndex[comp] - m_SpatialRadius[comp])); + regionIndex[comp] = vcl_max(static_cast<long int> (outputRegion.GetIndex().GetElement(comp)), + static_cast<long int> (inputIndex[comp] - m_SpatialRadius[comp])); const long int indexRight = vcl_min( - static_cast<long int>(outputRegion.GetIndex().GetElement(comp) + outputRegion.GetSize().GetElement(comp) - 1), - static_cast<long int>(inputIndex[comp] + m_SpatialRadius[comp])); + static_cast<long int> (outputRegion.GetIndex().GetElement(comp) + + outputRegion.GetSize().GetElement(comp) - 1), + static_cast<long int> (inputIndex[comp] + m_SpatialRadius[comp])); - regionSize[comp] = vcl_max(0l, indexRight - static_cast<long int>(regionIndex[comp]) + 1); + regionSize[comp] = vcl_max(0l, indexRight - static_cast<long int> (regionIndex[comp]) + 1); } RegionType neighborhoodRegion; @@ -400,85 +375,85 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati //itk::ImageRegionConstIterator<RealVectorImageType> it(jointImage, neighborhoodRegion); it.GoToBegin(); - while(!it.IsAtEnd()) + while (!it.IsAtEnd()) { - jointNeighbor.SetData(const_cast<RealType*>(it.GetPixelPointer())); + jointNeighbor.SetData(const_cast<RealType*> (it.GetPixelPointer())); // Compute the squared norm of the difference // This is the L2 norm, TODO: replace by the templated norm RealType norm2 = 0; - for(unsigned int comp = 0; comp < jointDimension; comp++) + for (unsigned int comp = 0; comp < jointDimension; comp++) { const RealType d = jointNeighbor[comp] - jointPixel[comp]; - norm2 += d*d; + norm2 += d * d; } // Compute pixel weight from kernel const RealType weight = m_Kernel(norm2); -/* - // The following code is an alternative way to compute norm2 and weight - // It separates the norms of spatial and range elements - RealType spatialNorm2; - RealType rangeNorm2; - spatialNorm2 = 0; - for (unsigned int comp = 0; comp < ImageDimension; comp++) - { - RealType d; - d = jointNeighbor[comp] - jointPixel[comp]; - spatialNorm2 += d*d; - } + /* + // The following code is an alternative way to compute norm2 and weight + // It separates the norms of spatial and range elements + RealType spatialNorm2; + RealType rangeNorm2; + spatialNorm2 = 0; + for (unsigned int comp = 0; comp < ImageDimension; comp++) + { + RealType d; + d = jointNeighbor[comp] - jointPixel[comp]; + spatialNorm2 += d*d; + } - if(spatialNorm2 >= 1.0) - { - weight = 0; - } - else - { - rangeNorm2 = 0; - for (unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++) - { - RealType d; - d = jointNeighbor[ImageDimension + comp] - jointPixel[ImageDimension + comp]; - rangeNorm2 += d*d; - } + if(spatialNorm2 >= 1.0) + { + weight = 0; + } + else + { + rangeNorm2 = 0; + for (unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++) + { + RealType d; + d = jointNeighbor[ImageDimension + comp] - jointPixel[ImageDimension + comp]; + rangeNorm2 += d*d; + } - weight = (rangeNorm2 <= 1.0)? 1.0 : 0.0; - } -*/ + weight = (rangeNorm2 <= 1.0)? 1.0 : 0.0; + } + */ // Update sum of weights weightSum += weight; // Update mean shift vector - for(unsigned int comp = 0; comp < jointDimension; comp++) + for (unsigned int comp = 0; comp < jointDimension; comp++) { meanShiftVector[comp] += weight * jointNeighbor[comp]; } - ++it; } - if(weightSum > 0) + if (weightSum > 0) { - for(unsigned int comp = 0; comp < jointDimension; comp++) + for (unsigned int comp = 0; comp < jointDimension; comp++) { meanShiftVector[comp] = meanShiftVector[comp] / weightSum - jointPixel[comp]; } } - } +} +#ifdef USEBUCKET // Calculates the mean shift vector at the position given by jointPixel -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> -void -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::CalculateMeanShiftVectorBucket(const RealVector& jointPixel, RealVector& meanShiftVector) - { +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::CalculateMeanShiftVectorBucket( + const RealVector& jointPixel, + RealVector& meanShiftVector) +{ const unsigned int jointDimension = ImageDimension + m_NumberOfComponentsPerPixel; RealType weightSum = 0; - for(unsigned int comp = 0; comp < jointDimension; comp++) + for (unsigned int comp = 0; comp < jointDimension; comp++) { meanShiftVector[comp] = 0; } @@ -486,34 +461,35 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati RealVector jointNeighbor(ImageDimension + m_NumberOfComponentsPerPixel); InputIndexType index; - for(unsigned int dim = 0; dim < ImageDimension; ++dim) + for (unsigned int dim = 0; dim < ImageDimension; ++dim) { index[dim] = jointPixel[dim] * m_SpatialBandwidth + 0.5; } - const std::vector<unsigned int> neighborBuckets = m_BucketImage.GetNeighborhoodBucketListIndices( - m_BucketImage.BucketIndexToBucketListIndex( - m_BucketImage.GetBucketIndex(jointPixel, index) - ) - ); + const std::vector<unsigned int> + neighborBuckets = m_BucketImage.GetNeighborhoodBucketListIndices( + m_BucketImage.BucketIndexToBucketListIndex( + m_BucketImage.GetBucketIndex( + jointPixel, + index))); unsigned int numNeighbors = m_BucketImage.GetNumberOfNeighborBuckets(); - for(unsigned int neighborIndex = 0; neighborIndex < numNeighbors; ++neighborIndex) + for (unsigned int neighborIndex = 0; neighborIndex < numNeighbors; ++neighborIndex) { const typename BucketImageType::BucketType & bucket = m_BucketImage.GetBucket(neighborBuckets[neighborIndex]); - if(bucket.empty()) continue; + if (bucket.empty()) continue; typename BucketImageType::BucketType::const_iterator it = bucket.begin(); - while(it != bucket.end()) + while (it != bucket.end()) { - jointNeighbor.SetData(const_cast<RealType*>(*it)); + jointNeighbor.SetData(const_cast<RealType*> (*it)); // Compute the squared norm of the difference // This is the L2 norm, TODO: replace by the templated norm RealType norm2 = 0; - for(unsigned int comp = 0; comp < jointDimension; comp++) + for (unsigned int comp = 0; comp < jointDimension; comp++) { const RealType d = jointNeighbor[comp] - jointPixel[comp]; - norm2 += d*d; + norm2 += d * d; } // Compute pixel weight from kernel @@ -523,7 +499,7 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati weightSum += weight; // Update mean shift vector - for(unsigned int comp = 0; comp < jointDimension; comp++) + for (unsigned int comp = 0; comp < jointDimension; comp++) { meanShiftVector[comp] += weight * jointNeighbor[comp]; } @@ -532,57 +508,56 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati } } - if(weightSum > 0) + if (weightSum > 0) { - for(unsigned int comp = 0; comp < jointDimension; comp++) + for (unsigned int comp = 0; comp < jointDimension; comp++) { meanShiftVector[comp] = meanShiftVector[comp] / weightSum - jointPixel[comp]; } } - } - +} +#endif -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> -void -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::ThreadedGenerateData(const OutputRegionType& outputRegionForThread, int threadId) +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::ThreadedGenerateData( + const OutputRegionType& outputRegionForThread, + int threadId) { // at the first iteration // Retrieve output images pointers - typename OutputSpatialImageType::Pointer spatialOutput = this->GetSpatialOutput(); - typename OutputImageType::Pointer rangeOutput = this->GetRangeOutput(); + typename OutputSpatialImageType::Pointer spatialOutput = this->GetSpatialOutput(); + typename OutputImageType::Pointer rangeOutput = this->GetRangeOutput(); typename OutputIterationImageType::Pointer iterationOutput = this->GetIterationOutput(); - typename OutputLabelImageType::Pointer labelOutput = this->GetLabelOutput(); + typename OutputLabelImageType::Pointer labelOutput = this->GetLabelOutput(); // Get input image pointer typename InputImageType::ConstPointer input = this->GetInput(); // defines input and output iterators - typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType; - typedef itk::ImageRegionIterator<OutputSpatialImageType> OutputSpatialIteratorType; - typedef itk::ImageRegionIterator<OutputIterationImageType> OutputIterationIteratorType; - typedef itk::ImageRegionIterator<OutputLabelImageType> OutputLabelIteratorType; - + typedef itk::ImageRegionIterator<OutputImageType> OutputIteratorType; + typedef itk::ImageRegionIterator<OutputSpatialImageType> OutputSpatialIteratorType; + typedef itk::ImageRegionIterator<OutputIterationImageType> OutputIterationIteratorType; + typedef itk::ImageRegionIterator<OutputLabelImageType> OutputLabelIteratorType; const unsigned int jointDimension = ImageDimension + m_NumberOfComponentsPerPixel; - typename OutputImageType::PixelType rangePixel(m_NumberOfComponentsPerPixel); - typename OutputSpatialImageType::PixelType spatialPixel(ImageDimension); + typename OutputImageType::PixelType rangePixel(m_NumberOfComponentsPerPixel); + typename OutputSpatialImageType::PixelType spatialPixel(ImageDimension); RealVector jointPixel; - RealVector bandwidth(jointDimension); - for (unsigned int comp = 0; comp < ImageDimension; comp++) bandwidth[comp] = m_SpatialBandwidth; - for (unsigned int comp = ImageDimension; comp < jointDimension; comp++) bandwidth[comp] = m_RangeBandwidth; + for (unsigned int comp = 0; comp < ImageDimension; comp++) + bandwidth[comp] = m_SpatialBandwidth; + for (unsigned int comp = ImageDimension; comp < jointDimension; comp++) + bandwidth[comp] = m_RangeBandwidth; itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); RegionType const& requestedRegion = input->GetRequestedRegion(); - typedef itk::ImageRegionConstIteratorWithIndex<RealVectorImageType> JointImageIteratorType; JointImageIteratorType jointIt(m_JointImage, outputRegionForThread); @@ -594,7 +569,6 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati typedef itk::ImageRegionIterator<ModeTableImageType> ModeTableImageIteratorType; ModeTableImageIteratorType modeTableIt(m_ModeTable, outputRegionForThread); - jointIt.GoToBegin(); rangeIt.GoToBegin(); spatialIt.GoToBegin(); @@ -610,21 +584,19 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati // Variables used by mode search optimization // List of indices where the current pixel passes through std::vector<InputIndexType> pointList; - if(m_ModeSearch) pointList.resize(m_MaxIterationNumber); + if (m_ModeSearch) pointList.resize(m_MaxIterationNumber); // Number of times an already processed candidate pixel is encountered, resulting in no // further computation (Used for statistics only) unsigned int numBreaks = 0; // index of the current pixel updated during the mean shift loop InputIndexType modeCandidate; - for (; !jointIt.IsAtEnd(); - ++jointIt, ++rangeIt, ++spatialIt, ++iterationIt, - ++modeTableIt, ++labelIt, progress.CompletedPixel()) + for (; !jointIt.IsAtEnd(); ++jointIt, ++rangeIt, ++spatialIt, ++iterationIt, ++modeTableIt, ++labelIt, progress.CompletedPixel()) { // if pixel has been already processed (by mode search optimization), skip typename ModeTableImageType::InternalPixelType const& currentPixelMode = modeTableIt.Get(); - if(m_ModeSearch && currentPixelMode == 1) + if (m_ModeSearch && currentPixelMode == 1) { numBreaks++; continue; @@ -658,7 +630,8 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati // but not 2 (pixel in current search path), and pixel has actually moved // from its initial position, and pixel candidate is inside the output // region, then perform optimization tasks - if (modeCandidate != currentIndex && m_ModeTable->GetPixel(modeCandidate) != 2 && outputRegionForThread.IsInside(modeCandidate)) + if (modeCandidate != currentIndex && m_ModeTable->GetPixel(modeCandidate) != 2 + && outputRegionForThread.IsInside(modeCandidate)) { // Obtain the data point to see if it close to jointPixel RealType diff = 0; @@ -666,20 +639,21 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati for (unsigned int comp = ImageDimension; comp < jointDimension; comp++) { const RealType d = candidatePixel[comp] - jointPixel[comp]; - diff += d*d; + diff += d * d; } if (diff < 0.5) // Spectral value is close enough { // If no mode has been associated to the candidate pixel then // associate it to the upcoming mode - if( m_ModeTable->GetPixel(modeCandidate) == 0) + if (m_ModeTable->GetPixel(modeCandidate) == 0) { // Add the candidate to the list of pixels that will be assigned the // finally calculated mode value pointList[pointCount++] = modeCandidate; m_ModeTable->SetPixel(modeCandidate, 2); - } else // == 1 + } + else // == 1 { // The candidate pixel has already been assigned to a mode // Assign the same value @@ -700,14 +674,19 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati } // end if (m_ModeSearch) //Calculate meanShiftVector - if(m_BucketOptimization) +#ifdef USEBUCKET + if (m_BucketOptimization) { this->CalculateMeanShiftVectorBucket(jointPixel, meanShiftVector); } else { +#endif this->CalculateMeanShiftVector(m_JointImage, jointPixel, requestedRegion, meanShiftVector); + +#ifdef USEBUCKET } +#endif // Compute mean shift vector squared norm (not normalized by bandwidth) // and add mean shift vector to current joint pixel @@ -715,7 +694,7 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati for (unsigned int comp = 0; comp < jointDimension; comp++) { const double v = meanShiftVector[comp] * bandwidth[comp]; - meanShiftVectorSqNorm += v*v; + meanShiftVectorSqNorm += v * v; jointPixel[comp] += meanShiftVector[comp]; } @@ -724,12 +703,12 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati iteration++; } - for(unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++) + for (unsigned int comp = 0; comp < m_NumberOfComponentsPerPixel; comp++) { rangePixel[comp] = jointPixel[ImageDimension + comp] * m_RangeBandwidth; } - for(unsigned int comp = 0; comp < ImageDimension; comp++) + for (unsigned int comp = 0; comp < ImageDimension; comp++) { spatialPixel[comp] = jointPixel[comp] * m_SpatialBandwidth - currentIndex[comp]; } @@ -751,7 +730,8 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati { m_NumLabels[threadId]++; label = m_NumLabels[threadId]; - } else // the loop exited through a break. Use the already assigned mode label + } + else // the loop exited through a break. Use the already assigned mode label { label = labelOutput->GetPixel(modeCandidate); } @@ -767,20 +747,17 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati } else // if ModeSearch is not set LabelOutput can't be generated { - LabelType labelZero=0; - labelIt.Set(labelZero); + LabelType labelZero = 0; + labelIt.Set(labelZero); } } // std::cout << "numBreaks: " << numBreaks << " Break ratio: " << numBreaks / (RealType)outputRegionForThread.GetNumberOfPixels() << std::endl; } - /* after threaded convergence test */ -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> -void -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::AfterThreadedGenerateData() +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::AfterThreadedGenerateData() { typename OutputLabelImageType::Pointer labelOutput = this->GetLabelOutput(); typedef itk::ImageRegionIterator<OutputLabelImageType> OutputLabelIteratorType; @@ -799,23 +776,25 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati { // Retrieve the number of labels in the thread by removing the threadId // from the most significant bits - LabelType localNumLabel = m_NumLabels[i-1] & ( (static_cast<LabelType>(1) << (sizeof(LabelType)*8 - m_ThreadIdNumberOfBits)) - static_cast<LabelType>(1) ); - newLabelOffset[i] = localNumLabel + newLabelOffset[i-1]; + LabelType localNumLabel = m_NumLabels[i - 1] & ((static_cast<LabelType> (1) << (sizeof(LabelType) * 8 + - m_ThreadIdNumberOfBits)) - static_cast<LabelType> (1)); + newLabelOffset[i] = localNumLabel + newLabelOffset[i - 1]; } labelIt.GoToBegin(); - while ( !labelIt.IsAtEnd() ) + while (!labelIt.IsAtEnd()) { LabelType const label = labelIt.Get(); // Get threadId from most significant bits - const unsigned int threadId = label >> (sizeof(LabelType)*8 - m_ThreadIdNumberOfBits); + const unsigned int threadId = label >> (sizeof(LabelType) * 8 - m_ThreadIdNumberOfBits); // Relabeling // First get the label number by removing the threadId bits // Then add the label offset specific to the threadId - LabelType newLabel = label & ( (static_cast<LabelType>(1) << (sizeof(LabelType)*8 - m_ThreadIdNumberOfBits)) - static_cast<LabelType>(1) ); + LabelType newLabel = label & ((static_cast<LabelType> (1) << (sizeof(LabelType) * 8 - m_ThreadIdNumberOfBits)) + - static_cast<LabelType> (1)); newLabel += newLabelOffset[threadId]; labelIt.Set(newLabel); @@ -825,16 +804,15 @@ MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterati } - -template <class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> -void -MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage> -::PrintSelf(std::ostream& os, itk::Indent indent) const +template<class TInputImage, class TOutputImage, class TKernel, class TOutputIterationImage> +void MeanShiftSmoothingImageFilter<TInputImage, TOutputImage, TKernel, TOutputIterationImage>::PrintSelf( + std::ostream& os, + itk::Indent indent) const { Superclass::PrintSelf(os, indent); - os << indent << "Spatial bandwidth: " << m_SpatialBandwidth << std::endl; - os << indent << "Range bandwidth: " << m_RangeBandwidth << std::endl; - } + os << indent << "Spatial bandwidth: " << m_SpatialBandwidth << std::endl; + os << indent << "Range bandwidth: " << m_RangeBandwidth << std::endl; +} } // end namespace otb