From 910711f1da1a25169ae0e5abeaabe2d3cc59de07 Mon Sep 17 00:00:00 2001
From: Julien Michel <julien.michel@c-s.fr>
Date: Wed, 18 Apr 2007 13:37:12 +0000
Subject: [PATCH] Ajout du filtre Spectral Angle.
---
.../otbSpectralAngleDistanceImageFilter.h | 110 +++++++++++++++++
.../otbSpectralAngleDistanceImageFilter.txx | 112 ++++++++++++++++++
Testing/Code/BasicFilters/CMakeLists.txt | 18 +++
.../BasicFilters/otbBasicFiltersTests.cxx | 2 +
.../otbSpectralAngleDistanceImageFilter.cxx | 71 +++++++++++
...otbSpectralAngleDistanceImageFilterNew.cxx | 50 ++++++++
6 files changed, 363 insertions(+)
create mode 100644 Code/BasicFilters/otbSpectralAngleDistanceImageFilter.h
create mode 100644 Code/BasicFilters/otbSpectralAngleDistanceImageFilter.txx
create mode 100644 Testing/Code/BasicFilters/otbSpectralAngleDistanceImageFilter.cxx
create mode 100644 Testing/Code/BasicFilters/otbSpectralAngleDistanceImageFilterNew.cxx
diff --git a/Code/BasicFilters/otbSpectralAngleDistanceImageFilter.h b/Code/BasicFilters/otbSpectralAngleDistanceImageFilter.h
new file mode 100644
index 0000000000..dd6a86899a
--- /dev/null
+++ b/Code/BasicFilters/otbSpectralAngleDistanceImageFilter.h
@@ -0,0 +1,110 @@
+/*=========================================================================
+
+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 _otbSpectralAngleDistanceImageFilter_h
+#define _otbSpectralAngleDistanceImageFilter_h
+
+#include "itkImageToImageFilter.h"
+
+namespace otb
+{
+/** \class SpectralAngleDistanceImageFilter
+ * \brief This filter implements the computation of the spectral angle
+ * distance with respect to a reference pixel.
+ *
+ * Spectral angle distance is an efficient way to convert spectral information to
+ * scalar data with respect to a reference pixel. The spectral angle distance
+ * is defined as:
+ *
+ * \f[SA = cos^{-1}\Big(\frac{{\sum_{b=1}^{n_{b}}}r(b)\cdot p(b)}{\sqrt{{\sum_{b=1}^{n_{b}}}r(b)^{2}{\sum_{b=1}^{n_{b}}p(b)^{2}}}}\Big) \f]
+ * with \f$b\f$ being the spectral band, \f$r\f$ the reference pixel and \f$p\f$
+ * the current pixel.
+ *
+ * Since the spectral angle deals with multi-bands image, the InputImage pixels are suposed to
+ * support the [] operator, and the input image to support the GetNumberOfComponentsPerPixel() method.
+ *
+ * \sa VectorImage
+ *
+ * This filter is implemented as a multithreaded filter.
+ *
+ * \ingroup IntensityImageFilters
+ * \ingroup Threading
+ * \ingroup Streamed
+ */
+template <class TInputImage,class TOutputImage>
+class ITK_EXPORT SpectralAngleDistanceImageFilter
+ : public itk::ImageToImageFilter<TInputImage,TOutputImage>
+{
+ public:
+ /** Standard typedefs */
+ typedef SpectralAngleDistanceImageFilter Self;
+ typedef itk::ImageToImageFilter<TInputImage,TOutputImage> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Type macro */
+ itkNewMacro(Self);
+
+ /** Creation through object factory macro */
+ itkTypeMacro(SpectralAngleDistanceImageFilter,itk::ImageToImageFilter);
+
+ /** Template parameters typedefs */
+ typedef TInputImage InputImageType;
+ typedef typename InputImageType::ConstPointer InputImageConstPointerType;
+ typedef typename InputImageType::RegionType InputImageRegionType;
+ typedef typename InputImageType::PixelType InputPixelType;
+ typedef TOutputImage OutputImageType;
+ typedef typename OutputImageType::Pointer OutputImagePointerType;
+ typedef typename OutputImageType::RegionType OutputImageRegionType;
+ typedef typename OutputImageType::PixelType OutputPixelType;
+ /** Get/Set the reference pixel */
+ itkGetConstReferenceMacro(ReferencePixel,InputPixelType);
+ itkSetMacro(ReferencePixel,InputPixelType);
+
+protected:
+ /** Constructor */
+ SpectralAngleDistanceImageFilter();
+ /** Destructor */
+ virtual ~SpectralAngleDistanceImageFilter() {};
+ /**PrintSelf method */
+ virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
+ /** SpectralAngleDistanceImageFilter 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 OutputImageRegionType& outputRegionForThread,
+ int threadId );
+
+private:
+ SpectralAngleDistanceImageFilter(const Self&); //purposely not implemented
+ void operator=(const Self&); //purposely not implemented
+
+ /** The reference pixel */
+ InputPixelType m_ReferencePixel;
+};
+}// End namespace otb
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSpectralAngleDistanceImageFilter.txx"
+#endif
+
+#endif
diff --git a/Code/BasicFilters/otbSpectralAngleDistanceImageFilter.txx b/Code/BasicFilters/otbSpectralAngleDistanceImageFilter.txx
new file mode 100644
index 0000000000..1acc68aeac
--- /dev/null
+++ b/Code/BasicFilters/otbSpectralAngleDistanceImageFilter.txx
@@ -0,0 +1,112 @@
+/*=========================================================================
+
+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 _otbSpectralAngleDistanceImageFilter_txx
+#define _otbSpectralAngleDistanceImageFilter_txx
+
+#include "otbSpectralAngleDistanceImageFilter.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkProgressReporter.h"
+#include "otbMacro.h"
+
+namespace otb
+{
+/**
+ * Constructor
+ */
+template <class TInputImage, class TOutputImage>
+SpectralAngleDistanceImageFilter<TInputImage,TOutputImage>
+::SpectralAngleDistanceImageFilter()
+{
+ m_ReferencePixel=0;
+}
+template <class TInputImage, class TOutputImage>
+void
+SpectralAngleDistanceImageFilter<TInputImage,TOutputImage>
+::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,int threadId)
+{
+
+ if (m_ReferencePixel == 0)
+ {
+ itkExceptionMacro(<<"Reference pixel is not set!");
+ }
+
+ InputImageConstPointerType inputPtr = this->GetInput();
+ OutputImagePointerType outputPtr = this->GetOutput();
+
+
+ // Define the portion of the input to walk for this thread, using
+ // the CallCopyOutputRegionToInputRegion method allows for the input
+ // and output images to be different dimensions
+ InputImageRegionType inputRegionForThread;
+ this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread);
+
+ // Define the iterators
+ itk::ImageRegionConstIterator<InputImageType> inputIt(inputPtr, inputRegionForThread);
+ itk::ImageRegionIterator<OutputImageType> outputIt(outputPtr, outputRegionForThread);
+ itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
+
+ inputIt.GoToBegin();
+ outputIt.GoToBegin();
+
+ while(!inputIt.IsAtEnd() && !outputIt.IsAtEnd())
+ {
+ double dist=0.0;
+ double scalarProd=0.0;
+ double normProd=0.0;
+ double normProd1=0.0;
+ double normProd2=0.0;
+ InputPixelType pixel = inputIt.Get();
+ for (int i=0; i<pixel.Size(); i++)
+ {
+ scalarProd += pixel[i]*m_ReferencePixel[i];
+ normProd1 += pixel[i]*pixel[i];
+ normProd2 += m_ReferencePixel[i]*m_ReferencePixel[i];
+ }
+ normProd = normProd1 * normProd2;
+
+ if ( normProd == 0.0)
+ {
+ dist = 0.0;
+ }
+ else
+ {
+ dist = vcl_acos(scalarProd/vcl_sqrt(normProd));
+ }
+ // Spectral angle normalisation
+ dist = dist/(M_PI/2);
+ //ponderation par un carre
+ dist = sqrt(dist)*255;
+ outputIt.Set(static_cast<OutputPixelType>(dist));
+ ++inputIt;
+ ++outputIt;
+ progress.CompletedPixel(); // potential exception thrown here
+ }
+}
+/**
+ * PrintSelf Method
+ */
+template <class TInputImage, class TOutputImage>
+void
+SpectralAngleDistanceImageFilter<TInputImage,TOutputImage>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ Superclass::PrintSelf(os, indent);
+}
+} // End namespace otb
+#endif
diff --git a/Testing/Code/BasicFilters/CMakeLists.txt b/Testing/Code/BasicFilters/CMakeLists.txt
index 2515b6d46c..bb2a347c01 100755
--- a/Testing/Code/BasicFilters/CMakeLists.txt
+++ b/Testing/Code/BasicFilters/CMakeLists.txt
@@ -198,6 +198,22 @@ ADD_TEST(bfTvStreamingShrinkImageFilterQBMUL ${BASICFILTERS_TESTS}
${TEMP}/bfTvImportGeoInformationOutput.img
)
+# ------- otb::SpectralAngleDistanceImageFilter ----------------------------
+
+ ADD_TEST(bfTuSpectralAngleDistanceImageFilterNew ${BASICFILTERS_TESTS}
+ otbSpectralAngleDistanceImageFilterNew)
+
+
+ ADD_TEST(bfTvSpectralAngleDistanceImageFilter ${BASICFILTERS_TESTS}
+ --compare-image ${TOL}
+ ${BASELINE}/bfTvSpectralAngleDistanceOutput.hdr
+ ${TEMP}/bfTvSpectralAngleDistanceOutput.hdr
+ otbSpectralAngleDistanceImageFilter
+ ${INPUTDATA}/poupees.jpg
+ ${TEMP}/bfTvSpectralAngleDistanceOutput.hdr
+ 210 432
+)
+
# A enrichir
SET(BasicFilters_SRCS
otbLeeFilter.cxx
@@ -223,6 +239,8 @@ otbSpatialObjectToImageDrawingFilterNew.cxx
otbSpatialObjectToImageDrawingFilter.cxx
otbImportGeoInformationImageFilterNew.cxx
otbImportGeoInformationImageFilter.cxx
+otbSpectralAngleDistanceImageFilterNew.cxx
+otbSpectralAngleDistanceImageFilter.cxx
)
diff --git a/Testing/Code/BasicFilters/otbBasicFiltersTests.cxx b/Testing/Code/BasicFilters/otbBasicFiltersTests.cxx
index ad491acd74..3e6495c129 100755
--- a/Testing/Code/BasicFilters/otbBasicFiltersTests.cxx
+++ b/Testing/Code/BasicFilters/otbBasicFiltersTests.cxx
@@ -50,4 +50,6 @@ REGISTER_TEST(otbSpatialObjectToImageDrawingFilterNew);
REGISTER_TEST(otbSpatialObjectToImageDrawingFilter);
REGISTER_TEST(otbImportGeoInformationImageFilterNew);
REGISTER_TEST(otbImportGeoInformationImageFilter);
+REGISTER_TEST(otbSpectralAngleDistanceImageFilterNew);
+REGISTER_TEST(otbSpectralAngleDistanceImageFilter);
}
diff --git a/Testing/Code/BasicFilters/otbSpectralAngleDistanceImageFilter.cxx b/Testing/Code/BasicFilters/otbSpectralAngleDistanceImageFilter.cxx
new file mode 100644
index 0000000000..fdc5d90e1a
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbSpectralAngleDistanceImageFilter.cxx
@@ -0,0 +1,71 @@
+/*=========================================================================
+
+ 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 "otbSpectralAngleDistanceImageFilter.h"
+
+int otbSpectralAngleDistanceImageFilter(int argc, char * argv[])
+{
+ try
+ {
+ const unsigned int Dimension = 2;
+ typedef double PixelType;
+ typedef otb::VectorImage<PixelType,Dimension> InputImageType;
+ typedef otb::Image<PixelType,Dimension> OutputImageType;
+ typedef otb::ImageFileReader<InputImageType> ReaderType;
+ typedef otb::ImageFileWriter<OutputImageType> WriterType;
+ typedef otb::SpectralAngleDistanceImageFilter<InputImageType,OutputImageType> DistanceFilterType;
+
+ // Instantiating object
+ DistanceFilterType::Pointer filter = DistanceFilterType::New();
+ ReaderType::Pointer reader = ReaderType::New();
+ WriterType::Pointer writer = WriterType::New();
+
+ reader->SetFileName(argv[1]);
+ writer->SetFileName(argv[2]);
+
+ InputImageType::IndexType refPixelIndex;
+ refPixelIndex[0]=atoi(argv[3]);
+ refPixelIndex[1]=atoi(argv[4]);
+
+ reader->Update();
+ filter->SetInput(reader->GetOutput());
+ filter->SetReferencePixel(reader->GetOutput()->GetPixel(refPixelIndex));
+
+ writer->SetInput(filter->GetOutput());
+
+ writer->Update();
+ }
+
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject thrown !" << std::endl;
+ std::cout << err << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ catch( ... )
+ {
+ std::cout << "Unknown exception thrown !" << std::endl;
+ return EXIT_FAILURE;
+ }
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/BasicFilters/otbSpectralAngleDistanceImageFilterNew.cxx b/Testing/Code/BasicFilters/otbSpectralAngleDistanceImageFilterNew.cxx
new file mode 100644
index 0000000000..02a9da54c6
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbSpectralAngleDistanceImageFilterNew.cxx
@@ -0,0 +1,50 @@
+/*=========================================================================
+
+ 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 "otbSpectralAngleDistanceImageFilter.h"
+
+int otbSpectralAngleDistanceImageFilterNew(int argc, char * argv[])
+{
+ try
+ {
+ const unsigned int Dimension = 2;
+ typedef double PixelType;
+ typedef otb::VectorImage<PixelType,Dimension> InputImageType;
+ typedef otb::Image<PixelType,Dimension> OutputImageType;
+ typedef otb::SpectralAngleDistanceImageFilter<InputImageType,OutputImageType> DistanceFilterType;
+
+ // Instantiating object
+ DistanceFilterType::Pointer object = DistanceFilterType::New();
+ }
+
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject thrown !" << std::endl;
+ std::cout << err << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ catch( ... )
+ {
+ std::cout << "Unknown exception thrown !" << std::endl;
+ return EXIT_FAILURE;
+ }
+ return EXIT_SUCCESS;
+}
--
GitLab