From cd32726d47e54ca3f5028464f42268c8795b79c2 Mon Sep 17 00:00:00 2001 From: Julien Michel <julien.michel@c-s.fr> Date: Tue, 9 Sep 2008 11:09:46 +0000 Subject: [PATCH] ENH: Adding a smoothing filter based on the mean shift algorithm. Not fully compliant yet. --- .../otbMeanShiftVectorImageFilter.h | 121 ++++++++ .../otbMeanShiftVectorImageFilter.txx | 263 ++++++++++++++++++ Testing/Code/BasicFilters/CMakeLists.txt | 20 +- .../BasicFilters/otbBasicFiltersTests9.cxx | 4 +- .../otbMeanShiftVectorImageFilter.cxx | 67 +++++ .../otbMeanShiftVectorImageFilterNew.cxx | 32 +++ 6 files changed, 496 insertions(+), 11 deletions(-) create mode 100644 Code/BasicFilters/otbMeanShiftVectorImageFilter.h create mode 100644 Code/BasicFilters/otbMeanShiftVectorImageFilter.txx create mode 100644 Testing/Code/BasicFilters/otbMeanShiftVectorImageFilter.cxx create mode 100644 Testing/Code/BasicFilters/otbMeanShiftVectorImageFilterNew.cxx diff --git a/Code/BasicFilters/otbMeanShiftVectorImageFilter.h b/Code/BasicFilters/otbMeanShiftVectorImageFilter.h new file mode 100644 index 0000000000..cbf621c67a --- /dev/null +++ b/Code/BasicFilters/otbMeanShiftVectorImageFilter.h @@ -0,0 +1,121 @@ +/*========================================================================= + +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 _otbMeanShiftVectorImageFilter_h +#define _otbMeanShiftVectorImageFilter_h + +#include "itkInPlaceImageFilter.h" + +namespace otb +{ + + /** \class MeanShiftVectorImageFilter + * + * + * + * + * \ingroup Streamed + * \ingroup Threaded + */ + +template <class TInputImage, class TOutputImage> +class ITK_EXPORT MeanShiftVectorImageFilter + : public itk::InPlaceImageFilter<TInputImage,TOutputImage> + { + public: + /** Standard class typedef */ + typedef MeanShiftVectorImageFilter Self; + typedef itk::InPlaceImageFilter<TInputImage,TOutputImage> Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + /** New and Type macros */ + itkNewMacro(Self); + itkTypeMacro(MeanShiftVectorImageFilter,InPlaceImageFilter); + + /** Template parameters typedefs */ + typedef TInputImage InputImageType; + typedef typename InputImageType::Pointer InputImagePointerType; + typedef typename InputImageType::PixelType InputPixelType; + typedef typename InputImageType::InternalPixelType InputInternalPixelType; + typedef typename InputImageType::PointType PointType; + typedef TOutputImage OutputImageType; + typedef typename OutputImageType::Pointer OutputImagePointerType; + typedef typename OutputImageType::PixelType OutputPixelType; + typedef typename OutputImageType::InternalPixelType OutputInternalPixelType; + typedef typename OutputImageType::RegionType RegionType; + typedef typename RegionType::SizeType SizeType; + + /** Setters / Getters */ + itkSetMacro(SpatialRadius,double); + itkGetMacro(SpatialRadius,double); + itkSetMacro(RangeRadius,double); + itkGetMacro(RangeRadius,double); + itkSetMacro(MaxNumberOfIterations,unsigned int); + itkGetMacro(MaxNumberOfIterations,unsigned int); + itkGetMacro(ConvergenceDistanceThreshold,double); + itkSetMacro(ConvergenceDistanceThreshold,double); + itkSetMacro(UseImageSpacing,bool); + itkGetMacro(UseImageSpacing,bool); + itkBooleanMacro(UseImageSpacing); + + protected: + /** This filters use a neighborhood around the pixel, so it needs to redfine the + * input requested region */ + virtual void GenerateInputRequestedRegion(); + + /** Threaded generate data */ + virtual void ThreadedGenerateData(const RegionType & outputRegionForThread,int threadId); + + /** Constructor */ + MeanShiftVectorImageFilter(); + + /** destructor */ + ~MeanShiftVectorImageFilter(){}; + + /**PrintSelf method */ + virtual void PrintSelf(std::ostream& os, itk::Indent indent) const; + + private: + MeanShiftVectorImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + /** Spatial radius for mean shift convergence */ + double m_SpatialRadius; + /** Range radius for mean shift convergence */ + double m_RangeRadius; + /** Max number of iterations for convergence */ + unsigned int m_MaxNumberOfIterations; + + /** Internal radius used by the iterator. The value depends on the spatial radius + * and the maximum number of iterations. This value is not intended to be accessed + * by users, thus no getter or setter is provided */ + SizeType m_InternalRadius; + + /** use image spacing */ + bool m_UseImageSpacing; + + /** Distance threshold for convergence in the spatial domain */ + double m_ConvergenceDistanceThreshold; + }; +}// end namespace otb + +#ifndef OTB_MANUAL_INSTANTIATION +#include "otbMeanShiftVectorImageFilter.txx" +#endif + +#endif diff --git a/Code/BasicFilters/otbMeanShiftVectorImageFilter.txx b/Code/BasicFilters/otbMeanShiftVectorImageFilter.txx new file mode 100644 index 0000000000..3f124649bb --- /dev/null +++ b/Code/BasicFilters/otbMeanShiftVectorImageFilter.txx @@ -0,0 +1,263 @@ +/*========================================================================= + +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 _otbMeanShiftVectorImageFilter_txx +#define _otbMeanShiftVectorImageFilter_txx + +#include "otbMeanShiftVectorImageFilter.h" + +#include "itkConstNeighborhoodIterator.h" +#include "itkImageRegionIterator.h" +#include "itkNeighborhoodAlgorithm.h" +#include "itkZeroFluxNeumannBoundaryCondition.h" +#include "itkOffset.h" +#include "itkProgressReporter.h" + +#include "otbMacro.h" + +namespace otb +{ + template <class TInputImage,class TOutputImage> + MeanShiftVectorImageFilter<TInputImage,TOutputImage> + ::MeanShiftVectorImageFilter() + { + m_MaxNumberOfIterations = 20; + m_ConvergenceDistanceThreshold = 0.5; + m_SpatialRadius = 3; + m_RangeRadius = 10; + m_UseImageSpacing = false; + m_InternalRadius.Fill(0); + } + + + template <class TInputImage,class TOutputImage> + void + MeanShiftVectorImageFilter<TInputImage,TOutputImage> + ::GenerateInputRequestedRegion() + { + // call the superclass' implementation of this method + Superclass::GenerateInputRequestedRegion(); + + // get pointers to the input and output + typename Superclass::InputImagePointer inputPtr = + const_cast< TInputImage * >( this->GetInput() ); + typename Superclass::OutputImagePointer outputPtr = this->GetOutput(); + + if ( !inputPtr || !outputPtr ) + { + return; + } + + // get a copy of the input requested region (should equal the output + // requested region) + typename TInputImage::RegionType inputRequestedRegion; + inputRequestedRegion = inputPtr->GetRequestedRegion(); + + //Compute the radius we will need to execute the filter: + if(m_UseImageSpacing) + { + for(unsigned int dimension = 0; dimension < InputImageType::ImageDimension;++dimension) + { + m_InternalRadius[dimension] = static_cast<unsigned long>(vcl_floor(m_MaxNumberOfIterations * m_SpatialRadius / inputPtr->GetSpacing()[0]+0.5)); + } + } + else + { + m_InternalRadius.Fill(static_cast<unsigned long>(vcl_floor(m_MaxNumberOfIterations * m_SpatialRadius+0.5))); + } + std::cout<<"MeanShiftVectorImageFilter: Internal radius = "<<m_InternalRadius<<std::endl; + + + // pad the input requested region by the operator radius + inputRequestedRegion.PadByRadius( m_InternalRadius ); + + // crop the input requested region at the input's largest possible region + if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) ) + { + inputPtr->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) + inputPtr->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(inputPtr); + throw e; + } + } + + template <class TInputImage,class TOutputImage> + void + MeanShiftVectorImageFilter<TInputImage,TOutputImage> + ::ThreadedGenerateData(const RegionType & outputRegionForThread,int threadId) + { + std::cout<<"Call to threaded generate data, threadId: "<<threadId<<", region: "<<outputRegionForThread<<std::endl; + // TODO: update this to handle the UseImageSpacing option. + + // Set up the progress reporter + itk::ProgressReporter progress(this,threadId,outputRegionForThread.GetNumberOfPixels()); + + // Input and output pointers + typename OutputImageType::Pointer outputPtr = this->GetOutput(); + typename InputImageType::ConstPointer inputPtr = this->GetInput(); + + // Iterators + itk::ConstNeighborhoodIterator<InputImageType> inputIt(m_InternalRadius,inputPtr,outputRegionForThread); + itk::ConstNeighborhoodIterator<InputImageType> internalIt; + itk::ImageRegionIterator<OutputImageType> outputIt(outputPtr,outputRegionForThread); + + // Boundary conditions + itk::ZeroFluxNeumannBoundaryCondition<InputImageType> nbc; + inputIt.OverrideBoundaryCondition(&nbc); + + //TODO: Work with a floating point pixel type + + // local declarations + unsigned int nbIterations; + typedef typename itk::ConstNeighborhoodIterator<InputImageType>::OffsetType OffsetType; + + OffsetType offset; + PointType maxDensityPoint,convergencePoint; + InputPixelType maxDensityValue, convergenceValue,diff,current; + unsigned int nbPixelsIntoAccount; + + unsigned int nbComponentsPerPixel = inputPtr->GetNumberOfComponentsPerPixel(); + bool goesOn = true; + long startx,stopx,starty,stopy; + + // intialize iterators + inputIt.GoToBegin(); + outputIt.GoToBegin(); + + // Walk the images + while(!inputIt.IsAtEnd() && !outputIt.IsAtEnd()) + { + nbIterations = 0; + goesOn = true; + convergencePoint.Fill(0); + convergenceValue = inputIt.GetCenterPixel(); + //std::cout<<"Processing pixel: "<<inputIt.GetIndex()<<", value: "<<convergenceValue<<std::endl; + + + // While the max number of iterations has not been reached and the convergence index is still moving + while(nbIterations < m_MaxNumberOfIterations && goesOn) + { + //std::cout<<"Loop: "<<nbIterations<<", goesOn: "<<goesOn<<std::endl; + maxDensityPoint.Fill(0); + maxDensityValue.SetSize(nbComponentsPerPixel); + maxDensityValue.Fill(0); + nbPixelsIntoAccount = 0; + + startx = vcl_floor(convergencePoint[0]-m_SpatialRadius); + stopx = vcl_ceil(convergencePoint[0]+m_SpatialRadius); + starty = vcl_floor(convergencePoint[1]-m_SpatialRadius); + stopy = vcl_ceil(convergencePoint[1]+m_SpatialRadius); + + // loop on the neighborhood + for(long i = startx;i<=stopx;++i) + { + for(long j =starty;j<=stopy;++j) + { + // Setting the offset + offset[0]=i; + offset[1]=j; + + // verify that the current offset is inside or spatial circle + if(vcl_sqrt(vcl_pow(static_cast<double>(convergencePoint[0]-offset[0]),2)+vcl_pow(static_cast<double>(convergencePoint[1]-offset[1]),2))<m_SpatialRadius) + { + //std::cout<<"Point: "<<offset<<" is inside spatial circle"<<std::endl; + current = inputIt.GetPixel(offset); + diff = convergenceValue - current; + //std::cout<<"Difference with current pixel: "<<diff<<std::endl; + // Check if we are inside the spectral - spatial circle + if(diff.GetSquaredNorm()+vcl_pow(static_cast<double>(offset[0]),2)+vcl_pow(static_cast<double>(offset[1]),2)<nbComponentsPerPixel * vcl_pow(m_RangeRadius,2) + vcl_pow(m_SpatialRadius,2)) + { + //std::cout<<"Point: "<<offset<<" is inside the spatial/spectral disk"<<std::endl; + maxDensityPoint[0]+=offset[0]; + maxDensityPoint[1]+=offset[1]; + maxDensityValue+=current; + ++nbPixelsIntoAccount; + } + } + } + } + // If no inside pixels + if(nbPixelsIntoAccount == 0) + { + maxDensityPoint.Fill(0); + maxDensityValue = inputIt.GetCenterPixel(); + } + else + { + maxDensityPoint[0]/=nbPixelsIntoAccount; + maxDensityPoint[1]/=nbPixelsIntoAccount; + maxDensityValue/=nbPixelsIntoAccount; + } + //std::cout<<"MaxDensityValue: "<<maxDensityValue<<std::endl; + convergenceValue = maxDensityValue; + + //std::cout<<"Nb pixels into account: "<<nbPixelsIntoAccount<<", convergence point: "<<convergencePoint<<", maxDensityPoint: "<<maxDensityPoint<<", convergenceValue: "<<convergenceValue<<std::endl; + + // Check if we are still significantly moving + double updatedDistance = vcl_sqrt(vcl_pow(static_cast<double>(maxDensityPoint[0]-convergencePoint[0]),2)+vcl_pow(static_cast<double>(maxDensityPoint[1]-convergencePoint[1]),2)); + if(updatedDistance < m_ConvergenceDistanceThreshold) + { + goesOn = false; + } + else + { + convergencePoint = maxDensityPoint; + } + ++nbIterations; + } + // Set the output value + // TODO: missing value cast here + outputIt.Set(convergenceValue); + // Update progress + progress.CompletedPixel(); + // Increment iterators + ++inputIt; + ++outputIt; + } + } + + template <class TInputImage,class TOutputImage> + void + MeanShiftVectorImageFilter<TInputImage,TOutputImage> + ::PrintSelf(std::ostream& os, itk::Indent indent) const + { + Superclass::PrintSelf(os,indent); + os<<indent<<"Spatial radius: "<<m_SpatialRadius<<std::endl; + os<<indent<<"Range radius: "<<m_RangeRadius<<std::endl; + os<<indent<<"Maximum number of iterations: "<<m_MaxNumberOfIterations<<std::endl; + os<<indent<<"Internal radius: "<<m_InternalRadius<<std::endl; + os<<indent<<"Use image spacing: "<<m_UseImageSpacing<<std::endl; + os<<indent<<"Convergence distance threshold: "<<m_ConvergenceDistanceThreshold<<std::endl; + + } + +} // end namespace otb + +#endif diff --git a/Testing/Code/BasicFilters/CMakeLists.txt b/Testing/Code/BasicFilters/CMakeLists.txt index f532cd7bf9..1d261725ae 100755 --- a/Testing/Code/BasicFilters/CMakeLists.txt +++ b/Testing/Code/BasicFilters/CMakeLists.txt @@ -883,13 +883,15 @@ ADD_TEST(bfTvContinuousMinimumMaximumImageCalculatorTest ${BASICFILTERS_TESTS9} # 2 2 # ) -#ADD_TEST(bfTuMeanShiftImageFilterNew ${BASICFILTERS_TESTS9} -# otbMeanShiftImageFilterNew ) - -#ADD_TEST(bfTvMeanShiftImageFilter ${BASICFILTERS_TESTS9} -# otbMeanShiftImageFilter ) - - +ADD_TEST(bfTuMeanShiftVectorImageFilterNew ${BASICFILTERS_TESTS9} + otbMeanShiftVectorImageFilterNew ) + +ADD_TEST(bfTvMeanShiftVectorImageFilter ${BASICFILTERS_TESTS9} + otbMeanShiftVectorImageFilter + ${INPUTDATA}/qb_RoadExtract2sub200x200.tif + ${TEMP}/bfTvvMeanShiftVectorImageFilterOutput.tif + 3 15 20 0 0.0001 + ) ADD_TEST(bfTvFunctionToImageFilterNew ${BASICFILTERS_TESTS9} otbFunctionToImageFilterNew @@ -1054,8 +1056,8 @@ otbContinuousMinimumMaximumImageCalculatorNew.cxx otbContinuousMinimumMaximumImageCalculatorTest.cxx #MeanShift.cxx #MeanShiftVectorImage.cxx -#otbMeanShiftImageFilterNew.cxx -#otbMeanShiftImageFilter.cxx +otbMeanShiftVectorImageFilterNew.cxx +otbMeanShiftVectorImageFilter.cxx otbFunctionToImageFilterNew.cxx otbFunctionToImageFilter.cxx otbScalarImageTextureFunctorNew.cxx diff --git a/Testing/Code/BasicFilters/otbBasicFiltersTests9.cxx b/Testing/Code/BasicFilters/otbBasicFiltersTests9.cxx index e752c5f9e9..b68119fcf2 100644 --- a/Testing/Code/BasicFilters/otbBasicFiltersTests9.cxx +++ b/Testing/Code/BasicFilters/otbBasicFiltersTests9.cxx @@ -31,8 +31,8 @@ REGISTER_TEST(otbChangeLabelImageFilterNew); REGISTER_TEST(otbChangeLabelImageFilterTest); REGISTER_TEST(otbContinuousMinimumMaximumImageCalculatorNew); REGISTER_TEST(otbContinuousMinimumMaximumImageCalculatorTest); -//REGISTER_TEST(otbMeanShiftImageFilterNew); -//REGISTER_TEST(otbMeanShiftImageFilter); +REGISTER_TEST(otbMeanShiftVectorImageFilterNew); +REGISTER_TEST(otbMeanShiftVectorImageFilter); //REGISTER_TEST(MeanShiftRef); //REGISTER_TEST(MeanShiftVectorImageRef); REGISTER_TEST(otbFunctionToImageFilterNew); diff --git a/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilter.cxx b/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilter.cxx new file mode 100644 index 0000000000..14ac41a5c6 --- /dev/null +++ b/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilter.cxx @@ -0,0 +1,67 @@ +/*========================================================================= + + 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. + +=========================================================================*/ +#include "itkExceptionObject.h" +#include "otbImage.h" +#include "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "otbMeanShiftVectorImageFilter.h" + +int otbMeanShiftVectorImageFilter(int argc, char * argv[]) +{ + if(argc != 8) + { + std::cerr<<"Usage: "<<argv[0]<<" infname outfname spatialRadius rangeRadius maxNbIterations useImageSpacing convergenceDistanceThreshold"<<std::endl; + return EXIT_FAILURE; + } + + const char * infname = argv[1]; + const char * outfname = argv[2]; + const double spatialRadius = atof(argv[3]); + const double rangeRadius = atof(argv[4]); + const unsigned int maxNbIterations = atoi(argv[5]); + const bool useImageSpacing = atoi(argv[6]); + const double convergenceTol = atof(argv[7]); + + const unsigned int Dimension = 2; + typedef double PixelType; + typedef otb::VectorImage<PixelType,Dimension> ImageType; + typedef otb::ImageFileReader<ImageType> ReaderType; + typedef otb::ImageFileWriter<ImageType> WriterType; + typedef otb::MeanShiftVectorImageFilter<ImageType,ImageType> FilterType; + + // Instantiating object + FilterType::Pointer filter = FilterType::New(); + ReaderType::Pointer reader = ReaderType::New(); + WriterType::Pointer writer = WriterType::New(); + + reader->SetFileName(infname); + writer->SetFileName(outfname); + + filter->SetSpatialRadius(spatialRadius); + filter->SetRangeRadius(rangeRadius); + filter->SetMaxNumberOfIterations(maxNbIterations); + filter->SetUseImageSpacing(useImageSpacing); + filter->SetConvergenceDistanceThreshold(convergenceTol); + + filter->SetInput(reader->GetOutput()); + writer->SetInput(filter->GetOutput()); + + writer->Update(); + return EXIT_SUCCESS; +} diff --git a/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilterNew.cxx b/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilterNew.cxx new file mode 100644 index 0000000000..8cb9aa85d7 --- /dev/null +++ b/Testing/Code/BasicFilters/otbMeanShiftVectorImageFilterNew.cxx @@ -0,0 +1,32 @@ +/*========================================================================= + + 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. + +=========================================================================*/ +#include "itkExceptionObject.h" +#include "otbVectorImage.h" +#include "otbMeanShiftVectorImageFilter.h" + +int otbMeanShiftVectorImageFilterNew(int argc, char * argv[]) +{ + const unsigned int Dimension = 2; + typedef double PixelType; + typedef otb::VectorImage<PixelType,Dimension> ImageType; + typedef otb::MeanShiftVectorImageFilter<ImageType,ImageType> FilterType; + + // Instantiating object + FilterType::Pointer object = FilterType::New(); + return EXIT_SUCCESS; +} -- GitLab