Skip to content
Snippets Groups Projects
Commit 7f165751 authored by Patrick Imbo's avatar Patrick Imbo
Browse files

ENH: HAlphaImageFilter class inherite from an UnaryFunctorImageFilter

parent 17a782b4
No related branches found
No related tags found
No related merge requests found
/*========================================================================= /*=========================================================================
Program: Insight Segmentation & Registration Toolkit Program: ORFEO Toolbox
Module: $RCSfile: otbHAlphaImageFilter.h,v $
Language: C++ Language: C++
Date: $Date: 2003/09/10 14:28:51 $ Date: $Date$
Version: $Revision: 1.4 $ Version: $Revision$
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 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. PURPOSE. See the above copyright notices for more information.
=========================================================================*/ =========================================================================*/
#ifndef __otbHAlphaImageFilter_h
#define __otbHAlphaImageFilter_h
#include <complex> #ifndef __HAlphaImageFilter_h
#include "math.h" #define __HAlphaImageFilter_h
#include "itkUnaryFunctorImageFilter.h"
#include "otbHermitianEigenAnalysis.h"
#include "itkVector.h"
#include "itkImageToImageFilter.h" namespace otb
#include "itkImage.h" {
#include "itkNumericTraits.h"
#include "itkImageRegionIterator.h"
#include "itkImageRegionConstIterator.h"
#include "itkImageRegionConstIteratorWithIndex.h"
#include "vnl/vnl_matrix_fixed.h"
#include "vnl/algo/vnl_generalized_eigensystem.h"
#include "itkNthElementImageAdaptor.h" namespace Functor {
#include "otbHermitianEigenAnalysis.h"
/** \class otbHAlphaFunctor
* \brief Evaluate the H-Alpha parameters from the coherency matrix image
*
* * Output value are:
* channel #0 : \f$ 0.5 * (S_{hh}+S_{vv}.(S_{hh}+S_{vv})^{*} \f$
* channel #1 : \f$ 0.5 * (S_{hh}+S_{vv}.(S_{hh}-S_{vv})^{*} \f$
* channel #2 : \f$ (S_{hh}+S_{vv}.(S_{hv})^{*} \f$
* channel #3 : \f$ 0.5 * (S_{hh}-S_{vv}.(S_{hh}-S_{vv})^{*} \f$
* channel #4 : \f$ (S_{hh}-S_{vv}.(S_{hv})^{*} \f$
* channel #5 : \f$ 2.0*S_{hv}.S_{hv}^{*} \f$
*
*/
template< class TInput, class TOutput>
class HAlphaFunctor
{
public:
typedef double RealType;
typedef typename std::complex <RealType> ComplexType;
typedef typename TOutput::ValueType OutputValueType;
/** D�claration des constantes */ /** CoherencyMatrix type **/
typedef itk::Vector<RealType,9> CoherencyMatrixType;
#ifndef PI /** Vector type used to store eigenvalues. */
#define PI 3.14159265358979323846 typedef itk::Vector<RealType, 3> EigenvalueType;
#endif
/** Matrix type used to store eigenvectors. */
typedef itk::Vector<float, 2> VectorType;
typedef itk::Vector<VectorType,3> EigenVectorFirstComposantType;
typedef itk::Vector<VectorType,3> EigenVectorType; // type d'un vecteur propre (partie r�elle, partie imaginaire)
typedef itk::Vector<itk::Vector<float, 6>,3> EigenMatrixType;
typedef itk::Image<EigenVectorType,2> EigenVectorImageType;
typedef itk::Image<double,2> EigenValueImageType;
namespace otb typedef itk::Vector<double, 3> OutputVectorType;
{
/** \class otbHAlphaImageFilter
* \brief Applies an averaging filter to an image
*
* Computes an image where a given pixel is the mean value of the
* the pixels in a neighborhood about the corresponding input pixel.
*
* A mean filter is one of the family of linear filters.
*
* \sa Image
* \sa Neighborhood
* \sa NeighborhoodOperator
* \sa NeighborhoodIterator
*
* \ingroup IntensityImageFilters
*/
typedef itk::Vector<float, 2> ComplexVectorType;
typedef itk::Vector<ComplexVectorType, 3> HermitianVectorType;
typedef itk::Vector<HermitianVectorType,3> HermitianMatrixType;
typedef otb::HermitianEigenAnalysis<CoherencyMatrixType,EigenvalueType, EigenMatrixType> HermitianAnalysisType;
template <class TPixel> inline TOutput operator()( const TInput & Coherency ) const
class HAlphaImageFilter : {
public itk::ImageToImageFilter< otb::Image<itk::Vector<TPixel,9>,2>, otb::Image<itk::Vector<TPixel,3>,2> > TOutput result;
result.SetSize(m_NumberOfComponentsPerPixel);
ComplexType T11 = static_cast<ComplexType>(Coherency[0]);
ComplexType T12 = static_cast<ComplexType>(Coherency[1]);
ComplexType T13 = static_cast<ComplexType>(Coherency[2]);
ComplexType T22 = static_cast<ComplexType>(Coherency[3]);
ComplexType T23 = static_cast<ComplexType>(Coherency[4]);
ComplexType T33 = static_cast<ComplexType>(Coherency[5]);
CoherencyMatrixType T;
EigenvalueType eigenValues;
EigenMatrixType eigenVectors;
HermitianAnalysisType HermitianAnalysis;
{
public:
/** Convenient typedefs for simplifying declarations. */
typedef otb::Image<itk::Vector<TPixel,9>,2> InputImageType;
typedef otb::Image<itk::Vector<TPixel,3>,2> OutputImageType;
/** Extract dimension from input and output image. */ T[0] = static_cast<RealType>(T11.real());
itkStaticConstMacro(InputImageDimension, unsigned int, T[1] = static_cast<RealType>(T22.real());
InputImageType::ImageDimension); T[2] = static_cast<RealType>(T33.real());
itkStaticConstMacro(OutputImageDimension, unsigned int, T[3] = static_cast<RealType>(T12.real());
OutputImageType::ImageDimension); T[4] = static_cast<RealType>(T12.imag());
T[5] = static_cast<RealType>(T13.real());
T[6] = static_cast<RealType>(T13.imag());
T[7] = static_cast<RealType>(T23.real());
T[8] = static_cast<RealType>(T23.imag());
HermitianAnalysis.ComputeEigenValuesAndVectors(T,eigenValues,eigenVectors);
// Entropy estimation
RealType totalEigenValues = 0.0;
RealType p[3];
RealType entropy;
totalEigenValues = static_cast<RealType>( eigenValues[0] + eigenValues[1] + eigenValues[2]);
if (totalEigenValues <0.00001)
totalEigenValues = 0.0001;
/** Standard class typedefs. */ for (unsigned int k = 0; k < 3; k++)
typedef HAlphaImageFilter Self; {
typedef itk::ImageToImageFilter< InputImageType, OutputImageType> Superclass; if (eigenValues[k] <0.)
typedef itk::SmartPointer<Self> Pointer; eigenValues[k] = 0.;
typedef itk::SmartPointer<const Self> ConstPointer;
p[k] = static_cast<RealType>(eigenValues[k]) / totalEigenValues;
}
/** Method for creation through the object factory. */ if ( (p[0]<0.0001) || (p[1]<0.0001) || (p[2]<0.0001) )
itkNewMacro(Self); {
entropy =0.0;
}
else
{
entropy = static_cast<RealType>( p[0]*log(p[0]) + p[1]*log(p[1]) + p[2]*log(p[2]) );
entropy /= -log(3.);
}
/** Run-time type information (and related methods). */ // alpha estimation
itkTypeMacro(HAlphaImageFilter, ImageToImageFilter); double val0, val1, val2;
double a0, a1, a2;
/** Image typedef support. */
typedef typename InputImageType::PixelType InputPixelType;
typedef typename OutputImageType::PixelType OutputPixelType;
typedef typename itk::NumericTraits<InputPixelType>::RealType InputRealType;
typedef typename InputImageType::RegionType InputImageRegionType;
typedef typename OutputImageType::RegionType OutputImageRegionType;
typedef typename InputImageType::SizeType InputSizeType; for(unsigned int k = 0; k < 3; k++)
typedef typename InputImageType::IndexType IndexType ; {
p[k] = eigenValues[k] / totalEigenValues;
if (p[k] < 0.) p[k] = 0.;
if (p[k] > 1.) p[k] = 1.;
}
val0=sqrt(eigenVectors[0][0]*eigenVectors[0][0]+eigenvectors[0][1]*eigenvectors[0][1]);
a0=acos(abs(val0)) **180./PI;
/** CoherencyMatrix type **/ val1=sqrt(eigenVectors[0][2]*eigenVectors[0][2]+eigenVectors[0][3]*m_EigenVectors[0][3]);
typedef itk::Vector<float,9> CoherencyMatrixType; a1=acos(abs(val1))*180./PI;
/** Vector type used to store eigenvalues. */ val2=sqrt(eigenVectors[0][4]*eigenVectors[0][4]+eigenvectors[0][5]*eigenvectors[0][5]);
typedef itk::Vector<float, 3> EigenvalueType; a2=acos(abs(val2))*180./PI;
/** Matrix type used to store eigenvectors. */ m_Alpha=p[0]*a0 + p[1]*a1 + p[2]*a2;
// typedef itk::Vector<itk::Vector<double,2>,3> VectorType;// 3 vecteurs propres partie r�elle + parties imaginaire = 6 composantes
typedef itk::Vector<float, 2> VectorType;
typedef itk::Vector<VectorType,3> EigenVectorFirstComposantType;
typedef itk::Vector<VectorType,3> EigenVectorType; // type d'un vecteur propre (partie r�elle, partie imaginaire)
typedef itk::Vector<itk::Vector<float, 6>,3> EigenMatrixType;
typedef itk::Image<EigenVectorType,2> EigenVectorImageType; if (m_Alpha>90) m_Alpha=0.;
typedef itk::Image<double,2> EigenValueImageType;
typedef itk::Vector<double, 3> OutputVectorType; // Anisotropy estimation
anisotropie=(eigenValues[1] - eigenValues[2])/(eigenValues[1] + eigenValues[2]+0.000001);
result[0] = entropy;
result[1] = alpha;
result[2] = anisotropy;
typedef itk::Vector<float, 2> ComplexVectorType; return result;
typedef itk::Vector<ComplexVectorType, 3> HermitianVectorType; }
typedef itk::Vector<HermitianVectorType,3> HermitianMatrixType;
unsigned int GetOutputSize()
{
return m_NumberOfComponentsPerPixel;
}
/** Constructor */
HAlphaFunctor() : m_NumberOfComponentsPerPixel(3) {}
/** Set the radius of the neighborhood used to compute the mean. */ /** Destructor */
itkSetMacro(Radius, InputSizeType); ~HAlphaFunctor() {}
/** Get the radius of the neighborhood used to compute the mean */ private:
itkGetConstReferenceMacro(Radius, InputSizeType); unsigned int m_NumberOfComponentsPerPixel;
};
}
/** \class otbHAlphaImageFilter
/** otbHAlphaImageFilter needs a larger input requested region than * \brief Compute the H-Alpha image (3 channels)
* the output requested region. As such, otbHAlphaImageFilter needs * from the coherency image (6 complex channels)
* to provide an implementation for GenerateInputRequestedRegion() */
* in order to inform the pipeline execution model. template <class TInputImage, class TOutputImage, class TFunction = Functor::HAlphaFunctor<
* ITK_TYPENAME TInputImage::PixelType, ITK_TYPENAME TOutputImage::PixelType> >
* \sa ImageToImageFilter::GenerateInputRequestedRegion() */ class ITK_EXPORT HAlphaImageFilter :
virtual void GenerateInputRequestedRegion() throw(itk::InvalidRequestedRegionError); public itk::UnaryFunctorImageFilter<TInputImage,TOutputImage, TFunction>
{
public:
/** Standard class typedefs. */
typedef HAlphaImageFilter Self;
typedef itk::UnaryFunctorImageFilter<TInputImage,TOutputImage, TFunction> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Runtime information support. */
itkTypeMacro(HAlphaImageFilter,itk::UnaryFunctorImageFilter);
protected: protected:
/**class typedefs for eigenvectors / eigenvalues */ HAlphaImageFilter() {}
/* typedef itk::EigenAnalysis3DImageFilter<InputImageType,EigenValueImageType, EigenVectorImageType> EigenAnalysisType; */
typedef otb::HermitianEigenAnalysis<CoherencyMatrixType,EigenvalueType, EigenMatrixType> HermitianAnalysisType;
HAlphaImageFilter();
virtual ~HAlphaImageFilter() {} virtual ~HAlphaImageFilter() {}
void PrintSelf(std::ostream& os, itk::Indent indent) const;
/** otbHAlphaImageFilter 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: private:
HAlphaImageFilter(const Self&); //purposely not implemented HAlphaImageFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented
InputSizeType m_Radius;
CoherencyMatrixType m_CoherencyMatrix;
EigenvalueType m_Eigenvalues;
EigenVectorFirstComposantType m_Eigenvectors; //r�cup�ration de la premi�re composante de chaque vecteur propre;
float m_Entropie;
float m_Alpha;
float m_Anisotropie;
}; };
} // end namespace otb } // end namespace otb
#ifndef ITK_MANUAL_INSTANTIATION
#include "otbHAlphaImageFilter.txx"
#endif
#endif #endif
...@@ -22,7 +22,6 @@ ...@@ -22,7 +22,6 @@
#include "itkExceptionObject.h" #include "itkExceptionObject.h"
#include <iostream> #include <iostream>
#include "otbImage.h"
#include "otbVectorImage.h" #include "otbVectorImage.h"
#include "otbHAlphaImageFilter.h" #include "otbHAlphaImageFilter.h"
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment