Commit a376e3ff authored by Guillaume Pasero's avatar Guillaume Pasero
Browse files

Remove EdisonMeanShift module

parent 2de54084
project(OTBEdisonMeanShift)
otb_module_impl()
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
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.
=========================================================================*/
#ifndef __otbMeanShiftImageFilter_h
#define __otbMeanShiftImageFilter_h
#include "vcl_deprecated_header.h"
#include "itkImageToImageFilter.h"
#include "itkVariableLengthVector.h"
#include "otbImage.h"
#include "otbObjectList.h"
#include "otbPolygon.h"
#include "otbLabeledOutputAccessor.h"
namespace otb
{
namespace MeanShift
{
/** \class ScalarBufferConverter
* \brief Small utilities class used for buffer conversion with EDISON standard.
* This class handles buffer conversion when pixel type is scalar.
* \deprecated
*
*
* \ingroup OTBEdisonMeanShift
**/
class ScalarBufferConverter
{
public:
/**
* Fill the pixel with the float array values at index, using nbComp values scaled by scale.
*/
template <class TPixel>
static inline void
FloatArrayToPixel(const float * data, unsigned int index, TPixel& pixel,
const unsigned int /*nbComp*/, double scale)
{
pixel = static_cast<TPixel>(scale * data[index]);
}
/**
* Fill the float array data with the pixel values at index, using nbComp values scaled by scale.
*/
template <class TPixel>
static inline void
PixelToFloatArray(float * data, unsigned int index, const TPixel& pixel, double scale)
{
data[index] = static_cast<float>(scale * pixel);
}
};
}
/** \class MeanShiftImageFilter
*
*
* Mean shift is a data clustering algorithm often used in image processing and segmentation.
* For a given pixel, the mean shift will build a set of neighboring pixels within a given spatial
* radius (can be set using SetSpatialRadius()) and a color range (can be set using SetRangeRadius()).
* The spatial and color center of this set is then computed and the algorithm iterates with this new spatial
* and color center.
*
* Mean shift can be used for edge-preserving smoothing, or for clustering. The GetOutput() method will allow
* you to get the smoothed image, whereas the
* GetClusteredOutput() methods returns the clustered output. The GetLabeledClusteredOutput() returns
* a labeled clustered image, and the GetClusterBoundariesOutput()
* an image of the cluster boundaries.
*
* The MinimumRegionSize parameter allows you to prune small clustered regions.
*
* This filter uses the Edison mean shift algorithm implementation. Please note that data whose precision
* is more than float are casted to float before processing.
*
* The Scale parameter allows you to stretch the data dynamic
*
* 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
*
* \deprecated use MeanShiftSegmentationFilter instead
*
* \sa MeanShiftSmootingImageFilter, MeanShiftSegmentationFilter
*
* \ingroup ImageEnhancement
*
* \ingroup OTBEdisonMeanShift
*/
template <class TInputImage, class TOutputImage,
class TLabeledOutput = otb::Image<unsigned short, 2>,
class TBufferConverter = MeanShift::ScalarBufferConverter>
class ITK_EXPORT MeanShiftImageFilter
: public itk::ImageToImageFilter<TInputImage, TOutputImage>
{
public:
/** Standard class typedef */
typedef MeanShiftImageFilter Self;
typedef itk::ImageToImageFilter<TInputImage, TOutputImage> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Type macro */
itkTypeMacro(MeanShiftImageFilter, ImageToImageFilter);
itkNewMacro(Self);
/** Template parameters typedefs */
typedef TInputImage InputImageType;
typedef typename InputImageType::Pointer InputImagePointerType;
typedef typename InputImageType::PixelType InputPixelType;
typedef typename InputImageType::PointType PointType;
typedef TOutputImage OutputImageType;
typedef typename OutputImageType::Pointer OutputImagePointerType;
typedef typename OutputImageType::PixelType OutputPixelType;
typedef typename OutputImageType::RegionType RegionType;
typedef TLabeledOutput LabeledOutputType;
typedef typename LabeledOutputType::Pointer LabeledOutputPointerType;
typedef typename LabeledOutputType::PixelType LabelType;
/** Typedefs for vectorized output */
typedef otb::Polygon<InputPixelType> PolygonType;
typedef typename PolygonType::Pointer PolygonPointerType;
typedef otb::ObjectList<PolygonType> PolygonListType;
typedef typename PolygonListType::Pointer PolygonListPointerType;
/** Typedef for mean-shift modes map */
typedef std::map<LabelType, InputPixelType> ModeMapType;
/** Setters / Getters */
itkSetMacro(SpatialRadius, unsigned int);
itkGetMacro(SpatialRadius, unsigned int);
itkSetMacro(RangeRadius, double);
itkGetMacro(RangeRadius, double);
itkSetMacro(MinimumRegionSize, unsigned int);
itkGetMacro(MinimumRegionSize, unsigned int);
itkSetMacro(Scale, double);
itkGetMacro(Scale, double);
/** Return the const clustered image output */
const OutputImageType * GetClusteredOutput() const;
/** Return the clustered image output */
OutputImageType * GetClusteredOutput();
/** Return the const labeled clustered image output */
const LabeledOutputType * GetLabeledClusteredOutput() const;
/** Return the labeled clustered image output */
LabeledOutputType * GetLabeledClusteredOutput();
/** Return the const cluster boundaries image output */
const LabeledOutputType * GetClusterBoundariesOutput() const;
/** Return the cluster boundaries image output */
LabeledOutputType * GetClusterBoundariesOutput();
/** Return the mean-shift mode by label */
const ModeMapType& GetModes()
{
return m_Modes;
}
protected:
virtual void EnlargeOutputRequestedRegion( itk::DataObject *output );
virtual void GenerateData();
/** Allocate the outputs (need to be reimplemented since outputs have differents type) */
virtual void AllocateOutputs();
/** If modified, we have to reset the list of modes */
virtual void Modified() const
{
Superclass::Modified();
m_Modes.clear();
}
/** Constructor */
MeanShiftImageFilter();
/** destructor */
virtual ~MeanShiftImageFilter() {}
/**PrintSelf method */
virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
private:
MeanShiftImageFilter(const Self &); //purposely not implemented
void operator =(const Self&); //purposely not implemented
/** Spatial radius for mean shift convergence */
unsigned int m_SpatialRadius;
/** Range radius for mean shift convergence */
double m_RangeRadius;
/** Minimum region size in pixels for clustering */
unsigned int m_MinimumRegionSize;
/** Data scale (used to stretch data range) */
double m_Scale;
/** A map of the different modes by segmented regions */
mutable ModeMapType m_Modes;
};
/**
* \class LabeledOutputAccessor
* \brief Specialized class to get the index of the labeled output image in mean shift filter.
*
* \ingroup OTBEdisonMeanShift
*/
template <class TInputImage, class TOutputImage, class TLabeledImage, class TBufferConverter>
class LabeledOutputAccessor<MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledImage, TBufferConverter> >
{
public:
typedef typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledImage, TBufferConverter>::LabeledOutputType LabelImageType;
itkStaticConstMacro(LabeledOutputIndex, unsigned int, 2);
};
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbMeanShiftImageFilter.txx"
#endif
#endif
/*=========================================================================
Program: ORFEO Toolbox
Language: C++
Date: $Date$
Version: $Revision$
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.
=========================================================================*/
#ifndef __otbMeanShiftImageFilter_txx
#define __otbMeanShiftImageFilter_txx
#include "otbMeanShiftImageFilter.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "otbMacro.h"
#include "msImageProcessor.h"
namespace otb
{
template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
::MeanShiftImageFilter()
{
m_SpatialRadius = 3;
m_RangeRadius = 10;
m_MinimumRegionSize = 10;
m_Scale = 100000.;
this->SetNumberOfRequiredOutputs(4);
this->SetNthOutput(1, OutputImageType::New());
this->SetNthOutput(2, LabeledOutputType::New());
this->SetNthOutput(3, LabeledOutputType::New());
}
template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
const typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::OutputImageType *
MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
::GetClusteredOutput() const
{
if (this->GetNumberOfOutputs() < 2)
{
return 0;
}
return static_cast<const OutputImageType *>(this->itk::ProcessObject::GetOutput(1));
}
template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::OutputImageType *
MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
::GetClusteredOutput()
{
if (this->GetNumberOfOutputs() < 2)
{
return 0;
}
return static_cast<OutputImageType *>(this->itk::ProcessObject::GetOutput(1));
}
template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
const typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::LabeledOutputType *
MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
::GetLabeledClusteredOutput() const
{
if (this->GetNumberOfOutputs() < 3)
{
return 0;
}
return static_cast<const LabeledOutputType *>(this->itk::ProcessObject::GetOutput(2));
}
template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::LabeledOutputType *
MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
::GetLabeledClusteredOutput()
{
if (this->GetNumberOfOutputs() < 3)
{
return 0;
}
return static_cast<LabeledOutputType *>(this->itk::ProcessObject::GetOutput(2));
}
template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
const typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::LabeledOutputType *
MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
::GetClusterBoundariesOutput() const
{
if (this->GetNumberOfOutputs() < 4)
{
return 0;
}
return static_cast<const LabeledOutputType *>(this->itk::ProcessObject::GetOutput(3));
}
template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
typename MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::LabeledOutputType *
MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
::GetClusterBoundariesOutput()
{
if (this->GetNumberOfOutputs() < 4)
{
return 0;
}
return static_cast<LabeledOutputType *>(this->itk::ProcessObject::GetOutput(3));
}
template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
void
MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
::AllocateOutputs()
{
typename OutputImageType::Pointer outputPtr = this->GetOutput();
typename OutputImageType::Pointer clusteredOutputPtr = this->GetClusteredOutput();
typename LabeledOutputType::Pointer labeledClusteredOutputPtr = this->GetLabeledClusteredOutput();
typename LabeledOutputType::Pointer clusterBoundariesOutputPtr = this->GetClusterBoundariesOutput();
outputPtr->SetBufferedRegion(outputPtr->GetRequestedRegion());
outputPtr->Allocate();
clusteredOutputPtr->SetBufferedRegion(clusteredOutputPtr->GetRequestedRegion());
clusteredOutputPtr->Allocate();
labeledClusteredOutputPtr->SetBufferedRegion(labeledClusteredOutputPtr->GetRequestedRegion());
labeledClusteredOutputPtr->Allocate();
clusterBoundariesOutputPtr->SetBufferedRegion(clusterBoundariesOutputPtr->GetRequestedRegion());
clusterBoundariesOutputPtr->Allocate();
}
template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
void
MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>::EnlargeOutputRequestedRegion( itk::DataObject *itkNotUsed(output) )
{
// This filter requires all of the output images in the buffer.
for ( unsigned int j = 0; j < this->GetNumberOfOutputs(); j++ )
{
if ( this->itk::ProcessObject::GetOutput(j) )
{
this->itk::ProcessObject::GetOutput(j)->SetRequestedRegionToLargestPossibleRegion();
}
}
}
template <class TInputImage, class TOutputImage, class TLabeledOutput, class TBufferConverter>
void
MeanShiftImageFilter<TInputImage, TOutputImage, TLabeledOutput, TBufferConverter>
::GenerateData()
{
// Allocate the output image
this->AllocateOutputs();
// Input and output pointers
typename InputImageType::ConstPointer inputPtr = this->GetInput();
typename OutputImageType::Pointer outputPtr = this->GetOutput();
double invScale = 1 / m_Scale;
RegionType inputRequestedRegion =inputPtr->GetRequestedRegion();
RegionType outputRequestedRegion = inputPtr->GetRequestedRegion();
inputRequestedRegion.PadByRadius(m_SpatialRadius);
if (!inputRequestedRegion.Crop(inputPtr->GetRequestedRegion()))
{
itkExceptionMacro(<< "error in requested region cropping");
}
// Iterators
itk::ImageRegionConstIteratorWithIndex<InputImageType> inputIt(inputPtr, inputRequestedRegion);
itk::ImageRegionIterator<OutputImageType> outputIt(outputPtr, outputRequestedRegion);
//create image processing object
msImageProcessor edisonProcessor;
float * data = new float[inputRequestedRegion.GetNumberOfPixels() * inputPtr->GetNumberOfComponentsPerPixel()];
unsigned int index = 0;
for (inputIt.GoToBegin(); !inputIt.IsAtEnd(); ++inputIt)
{
TBufferConverter::PixelToFloatArray(data, index, inputIt.Get(), m_Scale);
index += inputPtr->GetNumberOfComponentsPerPixel();
}
edisonProcessor.DefineLInput(data,
inputRequestedRegion.GetSize()[1],
inputRequestedRegion.GetSize()[0],
inputPtr->GetNumberOfComponentsPerPixel());
//define default kernel paramerters...
kernelType k[2] = {Uniform, Uniform};
int P[2] = {2, static_cast<int>(inputPtr->GetNumberOfComponentsPerPixel())};
float tempH[2] = {1.0, 1.0};
edisonProcessor.DefineKernel(k, tempH, P, 2);
edisonProcessor.Filter(m_SpatialRadius, m_RangeRadius * m_Scale, MED_SPEEDUP);
if (edisonProcessor.ErrorStatus)
{
itkExceptionMacro(<< "Error while running edison!");
}
typename OutputImageType::Pointer tmpOutput = OutputImageType::New();
tmpOutput->SetRegions(inputRequestedRegion);
tmpOutput->SetNumberOfComponentsPerPixel(outputPtr->GetNumberOfComponentsPerPixel());
tmpOutput->Allocate();
itk::ImageRegionIterator<OutputImageType> tmpIt(tmpOutput, inputRequestedRegion);
itk::ImageRegionIterator<OutputImageType> tmp2It(tmpOutput, outputRequestedRegion);
edisonProcessor.GetRawData(data);
if (edisonProcessor.ErrorStatus)
{
itkExceptionMacro(<< "Error while running edison!");
}
index = 0;
for (tmpIt.GoToBegin(); !tmpIt.IsAtEnd(); ++tmpIt)
{
OutputPixelType pixel;
TBufferConverter::FloatArrayToPixel(data, index, pixel, outputPtr->GetNumberOfComponentsPerPixel(), invScale);
tmpIt.Set(pixel);
index += outputPtr->GetNumberOfComponentsPerPixel();
}
tmp2It.GoToBegin();
outputIt.GoToBegin();
while (!tmp2It.IsAtEnd() && !outputIt.IsAtEnd())
{
outputIt.Set(tmp2It.Get());
++tmp2It;
++outputIt;
}
// typename OutputImageType::Pointer outputPtr = this->GetOutput();
typename OutputImageType::Pointer clusteredOutputPtr = this->GetClusteredOutput();
typename LabeledOutputType::Pointer labeledClusteredOutputPtr = this->GetLabeledClusteredOutput();
typename LabeledOutputType::Pointer clusterBoudariesOutputPtr = this->GetClusterBoundariesOutput();
itk::ImageRegionIterator<OutputImageType> clusteredOutputIt(clusteredOutputPtr, outputRequestedRegion);
index = 0;
for (outputIt.GoToBegin(); !outputIt.IsAtEnd(); ++outputIt)
{
TBufferConverter::PixelToFloatArray(data, index, outputIt.Get(), m_Scale);
index += outputPtr->GetNumberOfComponentsPerPixel();
}
edisonProcessor.DefineLInput(data,
outputRequestedRegion.GetSize()[1],
outputRequestedRegion.GetSize()[0],
outputPtr->GetNumberOfComponentsPerPixel());
// define default kernel paramerters...
//kernelType k[2] = {Uniform, Uniform};
//int P[2] = {2, outputPtr->GetNumberOfComponentsPerPixel()};
//float tempH[2] = {1.0, 1.0};
edisonProcessor.DefineKernel(k, tempH, P, 2);
edisonProcessor.FuseRegions(m_RangeRadius * m_Scale, m_MinimumRegionSize);
if (edisonProcessor.ErrorStatus)
{
itkExceptionMacro(<< "Error while running edison!");
}
edisonProcessor.GetRawData(data);
if (edisonProcessor.ErrorStatus)
{
itkExceptionMacro(<< "Error while running edison!");
}
index = 0;
for (clusteredOutputIt.GoToBegin(); !clusteredOutputIt.IsAtEnd(); ++clusteredOutputIt)
{
OutputPixelType pixel;
TBufferConverter::FloatArrayToPixel(data, index, pixel,
clusteredOutputPtr->GetNumberOfComponentsPerPixel(), invScale);
clusteredOutputIt.Set(pixel);
index += clusteredOutputPtr->GetNumberOfComponentsPerPixel();
}
int * labels = NULL;
float * modes = NULL;
int * modesPointsCount = NULL;
edisonProcessor.GetRegions(&labels, &modes, &modesPointsCount);
if (edisonProcessor.ErrorStatus)
{
itkExceptionMacro(<< "Error while running edison!");
}
itk::ImageRegionIteratorWithIndex<LabeledOutputType> lcIt(labeledClusteredOutputPtr,
labeledClusteredOutputPtr->GetRequestedRegion());
index = 0;
labeledClusteredOutputPtr->FillBuffer(0);
for (lcIt.GoToBegin(); !lcIt.IsAtEnd(); ++lcIt)
{
lcIt.Set(static_cast<LabelType>(labels[index]));
++index;
}
clusterBoudariesOutputPtr->FillBuffer(0);
//define the boundaries
RegionList * regionList = edisonProcessor.GetBoundaries();
int * regionIndeces;
unsigned int numRegions = regionList->GetNumRegions();
typename LabeledOutputType::IndexType boundIndex;
for (LabelType label = 0; label < static_cast<LabelType>(numRegions); ++label)
{
OutputPixelType pixel;
TBufferConverter::FloatArrayToPixel(modes,
static_cast<unsigned int>(label *
clusteredOutputPtr->GetNumberOfComponentsPerPixel()),
pixel, clusteredOutputPtr->GetNumberOfComponentsPerPixel(), invScale);
// Filling the modes map
m_Modes[label] = pixel;
regionIndeces = regionList->GetRegionIndeces(static_cast<int>(label));
for (int i = 0; i < regionList->GetRegionCount(static_cast<int>(label)); ++i)
{
boundIndex[0] =
(regionIndeces[i] %
clusterBoudariesOutputPtr->GetRequestedRegion().GetSize()[0]) +
clusterBoudariesOutputPtr->GetRequestedRegion().GetIndex()[0];
boundIndex[1] =
(regionIndeces[i] /
clusterBoudariesOutputPtr->GetRequestedRegion().GetSize()[0]) +
clusterBoudariesOutputPtr->GetRequestedRegion().GetIndex()[1];
if (clusterBoudariesOutputPtr->GetBufferedRegion().IsInside(boundIndex))
{
clusterBoudariesOutputPtr->SetPixel(boundIndex, 1);
}
}
}
// Free memory