Skip to content
Snippets Groups Projects
Commit de72c475 authored by Cyrille Valladeau's avatar Cyrille Valladeau
Browse files

ENH : supress Hermitian Analysis class

parent 78842bb5
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 145 additions and 1389 deletions
Code/SARPolarimetry/otbHermitianEigenAnalysis.h deleted 100644 → 0
+ 0
− 342
View file @ 78842bb5
/*=========================================================================
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 __otbHermitianEigenAnalysis_h
#define __otbHermitianEigenAnalysis_h
#include "itkMacro.h"
#include "itkVector.h"
#include "itkArray2D.h"
namespace otb
{
/** \class HermitianEigenAnalysis
* \brief Find Eigen values of a real 2D symmetric matrix. It
* serves as a thread safe alternative to the class:
* vnl_symmetric_eigensystem, which uses netlib routines.
*
* The class is templated over the input matrix, (which is expected to provide
* access to its elements with the [][] operator), matrix to store eigen
* values, (must provide write operations on its elements with the [] operator),
* EigenMatrix to store eigen vectors (must provide write access to its elements
* with the [][] operator).
*
* The SetOrderEigenValues() method can be used to order eigen values (and their
* corresponding eigen vectors if computed) in ascending order. This is the
* default ordering scheme. Eigen vectors and values can be obtained without
* ordering by calling SetOrderEigenValues(false)
*
* The SetOrderEigenMagnitudes() method can be used to order eigen values (and
* their corresponding eigen vectors if computed) by magnitude in ascending order.
*
* The user of this class is explicitly supposed to set the dimension of the
* 2D matrix using the SetDimension() method.
*
* The class contains routines taken from netlib sources. (www.netlib.org).
* netlib/tql1.c
* netlib/tql2.c
* netlib/tred1.c
* netlib/tred2.c
*
* Reference:
* num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
* wilkinson.
* handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
*
* \ingroup SARPolarimetry
*/
template < typename TMatrix, typename TVector, typename TEigenMatrix=TMatrix >
class HermitianEigenAnalysis
{
public:
typedef enum {
OrderByValue=1,
OrderByMagnitude,
DoNotOrder
}EigenValueOrderType;
typedef itk::Vector<float, 2> ComplexVectorType;
typedef itk::Vector<ComplexVectorType, 3> HermitianVectorType;
typedef itk::Vector<HermitianVectorType, 3> HermitianMatrixType;
HermitianEigenAnalysis():
m_Dimension(0),
m_Order(0),
m_OrderEigenValues(OrderByValue)
{};
HermitianEigenAnalysis( const unsigned int dimension ):
m_Dimension(dimension),
m_Order(dimension),
m_OrderEigenValues(OrderByValue)
{};
virtual ~HermitianEigenAnalysis() {};
typedef TMatrix MatrixType;
typedef TEigenMatrix EigenMatrixType;
typedef TVector VectorType;
/** Compute Eigen values of A
* A is any type that overloads the [][] operator and contains the
* symmetric matrix. In practice only the upper triangle of the
* matrix will be accessed. (Both itk::Matrix and vnl_matrix
* overload [][] operator.)
*
* 'EigenValues' is any type that overloads the [] operator and will contain
* the eigen values.
*
* No size checking is performed. A is expected to be a square matrix of size
* m_Dimension. 'EigenValues' is expected to be of length m_Dimension.
* The matrix is not checked to see if it is symmetric.
*/
unsigned int ComputeEigenValues(
const TMatrix & A,
TVector & EigenValues) const;
/** Compute Eigen values and vectors of A
* A is any type that overloads the [][] operator and contains the
* symmetric matrix. In practice only the upper triangle of the
* matrix will be accessed. (Both itk::Matrix and vnl_matrix
* overload [][] operator.)
*
* 'EigenValues' is any type that overloads the [] operator and will contain
* the eigen values.
*
* 'EigenVectors' is any type that provides access to its elements with the
* [][] operator. It is expected be of size m_Dimension * m_Dimension.
*
* No size checking is performed. A is expected to be a square matrix of size
* m_Dimension. 'EigenValues' is expected to be of length m_Dimension.
* The matrix is not checked to see if it is symmetric.
*
* Each row of the matrix 'EigenVectors' represents an eigen vector. (unlike MATLAB
* where the columns of the [EigenVectors, EigenValues] = eig(A) contains the
* eigenvectors).
*/
unsigned int ComputeEigenValuesAndVectors(
const TMatrix & A,
TVector & EigenValues,
TEigenMatrix & EigenVectors ) const;
/** Matrix order. Defaults to matrix dimension if not set **/
void SetOrder(const unsigned int n)
{
m_Order = n;
}
/** Get the Matrix order. Will be 0 unless explicitly set, or unless a
* call to SetDimension has been made in which case it will be the
* matrix dimension. */
unsigned int GetOrder() const { return m_Order; }
/** Set/Get methods to order the eigen values in ascending order.
* This is the default. ie lambda_1 < lambda_2 < ....
*/
void SetOrderEigenValues( const bool b )
{
if (b) { m_OrderEigenValues = OrderByValue; }
else { m_OrderEigenValues = DoNotOrder; }
}
bool GetOrderEigenValues() const { return (m_OrderEigenValues == OrderByValue); }
/** Set/Get methods to order the eigen value magnitudes in ascending order.
* In other words, |lambda_1| < |lambda_2| < .....
*/
void SetOrderEigenMagnitudes( const bool b )
{
if (b) { m_OrderEigenValues = OrderByMagnitude; }
else { m_OrderEigenValues = DoNotOrder; }
}
bool GetOrderEigenMagnitudes() const { return (m_OrderEigenValues == OrderByMagnitude); }
/** Set the dimension of the input matrix A. A is a square matrix of
* size m_Dimension. */
void SetDimension( const unsigned int n )
{
m_Dimension = n;
if (m_Order == 0 )
{
m_Order = m_Dimension;
}
}
/** Get Matrix dimension, Will be 0 unless explicitly set by a
* call to SetDimension. */
unsigned int GetDimension() const { return m_Dimension; }
private:
unsigned int m_Dimension;
unsigned int m_Order;
EigenValueOrderType m_OrderEigenValues;
/** Reduces a real symmetric matrix to a symmetric tridiagonal matrix using
* orthogonal similarity transformations.
* 'inputMatrix' contains the real symmetric input matrix. Only the lower
* triangle of the matrix need be supplied. The upper triangle is unaltered.
* 'd' contains the diagonal elements of the tridiagonal matrix.
* 'e' contains the subdiagonal elements of the tridiagonal matrix in its
* last n-1 positions. e(1) is set to zero.
* 'e2' contains the squares of the corresponding elements of e.
* 'e2' may coincide with e if the squares are not needed.
* questions and comments should be directed to burton s. garbow.
* mathematics and computer science div, argonne national laboratory
* this version dated august 1983.
*
* Function Adapted from netlib/tred1.c.
* [Changed: remove static vars, enforce const correctness.
* Use vnl routines as necessary].
* Reference:
* num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
* handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). */
void ReduceToTridiagonalMatrix(double *inputMatrix, VectorType &d,
double *e, double *e2) const;
/** Reduces a real symmetric matrix to a symmetric tridiagonal matrix using
* and accumulating orthogonal similarity transformations.
* 'inputMatrix' contains the real symmetric input matrix. Only the lower
* triangle of the matrix need be supplied. The upper triangle is unaltered.
* 'diagonalElements' will contains the diagonal elements of the tridiagonal
* matrix.
* 'subDiagonalElements' will contain the subdiagonal elements of the tridiagonal
* matrix in its last n-1 positions. subDiagonalElements(1) is set to zero.
* 'transformMatrix' contains the orthogonal transformation matrix produced
* in the reduction.
*
* questions and comments should be directed to burton s. garbow.
* mathematics and computer science div, argonne national laboratory
* this version dated august 1983.
*
* Function Adapted from netlib/tred2.c.
* [Changed: remove static vars, enforce const correctness.
* Use vnl routines as necessary].
* Reference:
* num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
* handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). */
void ReduceToTridiagonalMatrixAndGetTransformation(
double *inputMatrix, VectorType &diagonalElements,
double *subDiagonalElements, double *transformMatrix) const;
/* Finds the eigenvalues of a symmetric tridiagonal matrix by the ql method.
*
* On input:
* 'd' contains the diagonal elements of the input matrix.
* 'e' contains the subdiagonal elements of the input matrix
* in its last n-1 positions. e(1) is arbitrary.
* On Output:
* 'd' contains the eigenvalues.
* 'e' has been destroyed.
*
* Returns:
* zero for normal return,
* j if the j-th eigenvalue has not been
* determined after 30 iterations.
*
*
* Reference
* this subroutine is a translation of the algol procedure tql1,
* num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
* wilkinson.
* handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
*
* Questions and comments should be directed to burton s. garbow,
* mathematics and computer science div, argonne national laboratory
* this version dated august 1983.
*
* Function Adapted from netlib/tql1.c.
* [Changed: remove static vars, enforce const correctness.
* Use vnl routines as necessary] */
unsigned int ComputeEigenValuesUsingQL(
VectorType &d, double *e) const;
/* Finds the eigenvalues and eigenvectors of a symmetric tridiagonal matrix
* by the ql method.
*
* On input:
* 'd' contains the diagonal elements of the input matrix.
* 'e' contains the subdiagonal elements of the input matrix
* in its last n-1 positions. e(1) is arbitrary.
* 'z' contains the transformation matrix produced in the reduction by
* ReduceToTridiagonalMatrixAndGetTransformation(), if performed. If the
* eigenvectors of the tridiagonal matrix are desired, z must contain
* the identity matrix.
* On Output:
* 'd' contains the eigenvalues.
* 'e' has been destroyed.
* 'z' contains orthonormal eigenvectors of the symmetric tridiagonal
* (or full) matrix.
*
* Returns:
* zero for normal return,
* j if the j-th eigenvalue has not been
* determined after 1000 iterations.
*
* Reference
* this subroutine is a translation of the algol procedure tql1,
* num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
* wilkinson.
* handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
*
* Questions and comments should be directed to burton s. garbow,
* mathematics and computer science div, argonne national laboratory
* this version dated august 1983.
*
* Function Adapted from netlib/tql2.c.
* [Changed: remove static vars, enforce const correctness.
* Use vnl routines as necessary]
*/
unsigned int ComputeEigenValuesAndVectorsUsingQL(
VectorType &d, double *e, double *z) const;
};
template< typename TMatrix, typename TVector, typename TEigenMatrix >
std::ostream & operator<<(std::ostream& os,
const HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix > &s)
{
os << "[ClassType: HermitianEigenAnalysis]" << std::endl;
os << " Dimension : " << s.GetDimension() << std::endl;
os << " Order : " << s.GetOrder() << std::endl;
os << " OrderEigenValues: " << s.GetOrderEigenValues() << std::endl;
os << " OrderEigenMagnitudes: " << s.GetOrderEigenMagnitudes() << std::endl;
return os;
}
} // end namespace otb
#ifndef ITK_MANUAL_INSTANTIATION
#include "otbHermitianEigenAnalysis.txx"
#endif
#endif
Code/SARPolarimetry/otbHermitianEigenAnalysis.txx deleted 100644 → 0
+ 0
− 983
View file @ 78842bb5
This diff is collapsed.
Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.h
+ 100
− 57
View file @ de72c475
...@@ -20,9 +20,9 @@ ...@@ -20,9 +20,9 @@
#define __ReciprocalHAlphaImageFilter_h #define __ReciprocalHAlphaImageFilter_h
#include "otbUnaryFunctorImageFilter.h" #include "otbUnaryFunctorImageFilter.h"
#include "otbHermitianEigenAnalysis.h"
#include "itkVector.h" #include "itkVector.h"
#include "otbMath.h" #include "otbMath.h"
#include "vnl/algo/vnl_complex_eigensystem.h"
namespace otb namespace otb
{ {
...@@ -44,60 +44,108 @@ template< class TInput, class TOutput> ...@@ -44,60 +44,108 @@ template< class TInput, class TOutput>
class ReciprocalHAlphaFunctor class ReciprocalHAlphaFunctor
{ {
public: public:
typedef typename std::complex <double> ComplexType; typedef typename std::complex<double> ComplexType;
typedef typename TOutput::ValueType OutputValueType; typedef vnl_matrix<ComplexType> VNLMatrixType;
typedef vnl_vector<ComplexType> VNLVectorType;
/** CoherencyMatrix type **/ typedef vnl_vector<double> VNLDoubleVectorType;
typedef itk::Vector<double, 9> CoherencyMatrixType; typedef typename TOutput::ValueType OutputValueType;
/** Vector type used to store eigenvalues. */
typedef itk::Vector<double, 3> EigenvalueType;
/** Matrix type used to store eigenvectors. */
typedef itk::Vector<float, 2> VectorType;
typedef itk::Vector<VectorType, 3> EigenVectorFirstComposantType;
typedef itk::Vector<VectorType, 3> EigenVectorType; // eigenvector type (real part, imaginary part)
typedef itk::Vector<itk::Vector<float, 6>, 3> EigenMatrixType;
typedef itk::Image<EigenVectorType, 2> EigenVectorImageType;
typedef itk::Image<double, 2> EigenValueImageType;
typedef itk::Vector<double, 3> OutputVectorType;
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;
inline TOutput operator()( const TInput & Coherency ) const inline TOutput operator()( const TInput & Coherency ) const
{ {
TOutput result; TOutput result;
result.SetSize(m_NumberOfComponentsPerPixel); result.SetSize(m_NumberOfComponentsPerPixel);
CoherencyMatrixType T; double T0 = static_cast<double>(Coherency[0].real());
EigenvalueType eigenValues; double T1 = static_cast<double>(Coherency[3].real());
EigenMatrixType eigenVectors; double T2 = static_cast<double>(Coherency[5].real());
HermitianAnalysisType HermitianAnalysis; double T3 = static_cast<double>(Coherency[1].real());
double T4 = static_cast<double>(Coherency[1].imag());
T[0] = static_cast<double>(Coherency[0].real()); double T5 = static_cast<double>(Coherency[2].real());
T[1] = static_cast<double>(Coherency[3].real()); double T6 = static_cast<double>(Coherency[2].imag());
T[2] = static_cast<double>(Coherency[5].real()); double T7 = static_cast<double>(Coherency[4].real());
T[3] = static_cast<double>(Coherency[1].real()); double T8 = static_cast<double>(Coherency[4].imag());
T[4] = static_cast<double>(Coherency[1].imag());
T[5] = static_cast<double>(Coherency[2].real()); VNLMatrixType vnlMat(3, 3, 0.);
T[6] = static_cast<double>(Coherency[2].imag()); vnlMat[0][0] = ComplexType(T0, 0.);
T[7] = static_cast<double>(Coherency[4].real()); vnlMat[0][1] = std::conj(ComplexType(Coherency[1]));
T[8] = static_cast<double>(Coherency[4].imag()); vnlMat[0][2] = std::conj(ComplexType(Coherency[2]));
HermitianAnalysis.ComputeEigenValuesAndVectors(T, eigenValues, eigenVectors); vnlMat[1][0] = ComplexType(Coherency[1]);
vnlMat[1][1] = ComplexType(T1, 0.);
vnlMat[1][2] = std::conj(ComplexType(Coherency[4]));
vnlMat[2][0] = ComplexType(Coherency[2]);
vnlMat[2][1] = ComplexType(Coherency[4]);
vnlMat[2][2] = ComplexType(T2, 0.);
// Only compute the left symetry to respect the previous Hermitian Analisys code
vnl_complex_eigensystem syst(vnlMat, false, true);
const VNLMatrixType eigenVectors( syst.L );
const VNLVectorType eigenValues(syst.W);
// Entropy estimation // Entropy estimation
double totalEigenValues(0.0); double totalEigenValues(0.0);
double p[3]; double p[3];
double entropy; double entropy;
double alpha; double alpha;
double anisotropy; double anisotropy;
// Sort eigne values and adapt eigen vectors first component
VNLDoubleVectorType sortedRealEigenValues(3, eigenValues[0].real());
VNLVectorType sortedEigenVectors(3, eigenVectors[0][0]);
if( (eigenValues[1].real() >= eigenValues[0].real()) && (eigenValues[1].real() >= eigenValues[2].real()) )
{
sortedRealEigenValues[0] = eigenValues[1].real();
sortedEigenVectors[0] = eigenVectors[1][0];
if( sortedRealEigenValues[2] >= eigenValues[0].real() )
{
sortedRealEigenValues[1] = eigenValues[2].real();
sortedEigenVectors[1] = eigenVectors[2][0];
}
else
{
sortedRealEigenValues[2] = eigenValues[2].real();
sortedEigenVectors[2] = eigenVectors[2][0];
//init do that sortedRealEigenValues[0] = eigenValues[1].real();
}
}
else if( (eigenValues[2].real() >= eigenValues[0].real()) && (eigenValues[2].real() >= eigenValues[1].real()) )
{
sortedRealEigenValues[0] = eigenValues[2].real();
sortedEigenVectors[0] = eigenVectors[2][0];
if( sortedRealEigenValues[1] >= eigenValues[0].real() )
{
sortedRealEigenValues[1] = eigenValues[1].real();
sortedEigenVectors[1] = eigenVectors[1][0];
}
else
{
sortedRealEigenValues[2] = eigenValues[1].real();
sortedEigenVectors[2] = eigenVectors[1][0];
//init do that sortedRealEigenValues[0] = eigenValues[1].real();
}
}
else
{
if( eigenValues[1].real() >= eigenValues[2].real() )
{
sortedRealEigenValues[1] = eigenValues[1].real();
sortedRealEigenValues[2] = eigenValues[2].real();
sortedEigenVectors[1] = eigenVectors[1][0];
sortedEigenVectors[2] = eigenVectors[2][0];
}
else
{
sortedRealEigenValues[1] = eigenValues[2].real();
sortedRealEigenValues[2] = eigenValues[1].real();
sortedEigenVectors[1] = eigenVectors[2][0];
sortedEigenVectors[2] = eigenVectors[1][0];
}
}
totalEigenValues = static_cast<double>( eigenValues[0] + eigenValues[1] + eigenValues[2]);
totalEigenValues = sortedRealEigenValues[0] + sortedRealEigenValues[1] + sortedRealEigenValues[2];
if (totalEigenValues <m_Epsilon) if (totalEigenValues <m_Epsilon)
{ {
totalEigenValues = m_Epsilon; totalEigenValues = m_Epsilon;
...@@ -105,12 +153,7 @@ public: ...@@ -105,12 +153,7 @@ public:
for (unsigned int k = 0; k < 3; k++) for (unsigned int k = 0; k < 3; k++)
{ {
if (eigenValues[k] <0.) p[k] = std::max(sortedRealEigenValues[k], 0.) / totalEigenValues;
{
eigenValues[k] = 0.;
}
p[k] = static_cast<double>(eigenValues[k]) / totalEigenValues;
} }
if ( (p[0] < m_Epsilon) || (p[1] < m_Epsilon) || (p[2] < m_Epsilon) ) if ( (p[0] < m_Epsilon) || (p[1] < m_Epsilon) || (p[2] < m_Epsilon) )
...@@ -129,7 +172,7 @@ public: ...@@ -129,7 +172,7 @@ public:
for(unsigned int k = 0; k < 3; k++) for(unsigned int k = 0; k < 3; k++)
{ {
p[k] = eigenValues[k] / totalEigenValues; p[k] = sortedRealEigenValues[k] / totalEigenValues;
if (p[k] < 0.) if (p[k] < 0.)
{ {
...@@ -142,13 +185,13 @@ public: ...@@ -142,13 +185,13 @@ public:
} }
} }
val0=sqrt(static_cast<double>(eigenVectors[0][0]*eigenVectors[0][0]) + static_cast<double>(eigenVectors[0][1]*eigenVectors[0][1])); val0 = std::abs(sortedEigenVectors[0]);
a0=acos(vcl_abs(val0)) * CONST_180_PI; a0=acos(vcl_abs(val0)) * CONST_180_PI;
val1=sqrt(static_cast<double>(eigenVectors[0][2]*eigenVectors[0][2]) + static_cast<double>(eigenVectors[0][3]*eigenVectors[0][3])); val1 = std::abs(sortedEigenVectors[1]);
a1=acos(vcl_abs(val1)) * CONST_180_PI; a1=acos(vcl_abs(val1)) * CONST_180_PI;
val2=sqrt(static_cast<double>(eigenVectors[0][4]*eigenVectors[0][4]) + static_cast<double>(eigenVectors[0][5]*eigenVectors[0][5])); val2= std::abs(sortedEigenVectors[2]);
a2=acos(vcl_abs(val2)) * CONST_180_PI; a2=acos(vcl_abs(val2)) * CONST_180_PI;
alpha=p[0]*a0 + p[1]*a1 + p[2]*a2; alpha=p[0]*a0 + p[1]*a1 + p[2]*a2;
...@@ -159,11 +202,11 @@ public: ...@@ -159,11 +202,11 @@ public:
} }
// Anisotropy estimation // Anisotropy estimation
anisotropy=(eigenValues[1] - eigenValues[2])/(eigenValues[1] + eigenValues[2]+m_Epsilon); anisotropy=(sortedRealEigenValues[1] - sortedRealEigenValues[2])/(sortedRealEigenValues[1] + sortedRealEigenValues[2] + m_Epsilon);
result[0] = entropy; result[0] = static_cast<OutputValueType>(entropy);
result[1] = alpha; result[1] = static_cast<OutputValueType>(alpha);
result[2] = anisotropy; result[2] = static_cast<OutputValueType>(anisotropy);
return result; return result;
} }
......
Testing/Code/SARPolarimetry/CMakeLists.txt
+ 0
− 5
View file @ de72c475
...@@ -373,10 +373,6 @@ ADD_TEST(saTvMuellerToCovarianceImageFilter ${SARPOLARIMETRY_TESTS2} ...@@ -373,10 +373,6 @@ ADD_TEST(saTvMuellerToCovarianceImageFilter ${SARPOLARIMETRY_TESTS2}
${BASELINE}/saTvSinclairImageFilter_SinclairToMueller.tif ${BASELINE}/saTvSinclairImageFilter_SinclairToMueller.tif
${TEMP}/saTvMuellerToMLCImageFilter.tif ${TEMP}/saTvMuellerToMLCImageFilter.tif
) )
# Hermitian eigen analysis class
ADD_TEST(saTvHermitianEigenAnalysisTest ${SARPOLARIMETRY_TESTS2}
otbHermitianEigenAnalysisTest
)
# Polarimetric Data # Polarimetric Data
ADD_TEST(saTuPolarimetricDataNew ${SARPOLARIMETRY_TESTS2} ADD_TEST(saTuPolarimetricDataNew ${SARPOLARIMETRY_TESTS2}
...@@ -430,7 +426,6 @@ otbMuellerToPolarisationDegreeAndPowerImageFilterNew.cxx ...@@ -430,7 +426,6 @@ otbMuellerToPolarisationDegreeAndPowerImageFilterNew.cxx
otbMuellerToPolarisationDegreeAndPowerImageFilter.cxx otbMuellerToPolarisationDegreeAndPowerImageFilter.cxx
otbMuellerToCovarianceImageFilterNew.cxx otbMuellerToCovarianceImageFilterNew.cxx
otbMuellerToCovarianceImageFilter.cxx otbMuellerToCovarianceImageFilter.cxx
otbHermitianEigenAnalysisTest.cxx
otbPolarimetricData.cxx otbPolarimetricData.cxx
) )
......
Testing/Code/SARPolarimetry/otbHermitianEigenAnalysisTest.cxx
+ 41
− 0
View file @ de72c475
...@@ -26,6 +26,8 @@ ...@@ -26,6 +26,8 @@
#include "otbVectorImage.h" #include "otbVectorImage.h"
#include "otbHermitianEigenAnalysis.h" #include "otbHermitianEigenAnalysis.h"
#include "itkVariableLengthVector.h" #include "itkVariableLengthVector.h"
#include "vnl/algo/vnl_complex_eigensystem.h"
#include <complex>
int otbHermitianEigenAnalysisTest(int argc, char * argv[]) int otbHermitianEigenAnalysisTest(int argc, char * argv[])
...@@ -37,6 +39,7 @@ int otbHermitianEigenAnalysisTest(int argc, char * argv[]) ...@@ -37,6 +39,7 @@ int otbHermitianEigenAnalysisTest(int argc, char * argv[])
typedef otb::HermitianEigenAnalysis<MatrixType, EigenvalueType, EigenMatrixType> FilterType; typedef otb::HermitianEigenAnalysis<MatrixType, EigenvalueType, EigenMatrixType> FilterType;
EigenMatrixType resEigVal; EigenMatrixType resEigVal;
itk::Vector<float, 6> temp; itk::Vector<float, 6> temp;
temp[0] = 0.162793; temp[0] = 0.162793;
temp[1] = -0.432753; temp[1] = -0.432753;
...@@ -115,5 +118,43 @@ int otbHermitianEigenAnalysisTest(int argc, char * argv[]) ...@@ -115,5 +118,43 @@ int otbHermitianEigenAnalysisTest(int argc, char * argv[])
return EXIT_FAILURE; return EXIT_FAILURE;
} }
std::cout<<"Values: "<<vect<<std::endl;
std::cout<<"Vectors: "<<eigVal<<std::endl;
typedef std::complex<double> ComplexType;
vnl_matrix<ComplexType> vnlMat(3,3,0);
vnlMat[0][0] = ComplexType(0., 0.);
vnlMat[0][1] = ComplexType(1.5, 2.);
vnlMat[0][2] = ComplexType(2.5, 3.);
vnlMat[1][0] = ComplexType(1.5, -2.);
vnlMat[1][1] = ComplexType(0.5, 0.);
vnlMat[1][2] = ComplexType(3.5, 4.);
vnlMat[2][0] = ComplexType(2.5, -3.);
vnlMat[2][1] = ComplexType(3.5, -4.);
vnlMat[2][2] = ComplexType(1., 0.);
std::cout<<"Matrix:: "<<std::endl;
vnlMat.print(std::cout);
vnl_complex_eigensystem syst(vnlMat, true, true);
vnl_matrix< ComplexType > pm = syst.L;
std::cout<<"Left:: "<<std::endl;
pm.print(std::cout);
std::cout<<"Right:: "<<std::endl;
pm = syst.R;
pm.print(std::cout);
std::cout<<"W:: "<<std::endl;
vnl_vector< ComplexType > lol = syst.W;
for(unsigned i=0; i<lol.size(); i++)
{
std::cout<<" "<<lol[i];
}
std::cout<<std::endl;
return EXIT_SUCCESS; return EXIT_SUCCESS;
} }
Testing/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.cxx
+ 4
− 1
View file @ de72c475
...@@ -89,9 +89,10 @@ int otbReciprocalHAlphaImageFilter(int argc, char * argv[]) ...@@ -89,9 +89,10 @@ int otbReciprocalHAlphaImageFilter(int argc, char * argv[])
PerBandMeanFilterType::Pointer perBandMeanFilter = PerBandMeanFilterType::New(); PerBandMeanFilterType::Pointer perBandMeanFilter = PerBandMeanFilterType::New();
perBandMeanFilter->SetInput(sinclairToCoherencyFilter->GetOutput()); perBandMeanFilter->SetInput(sinclairToCoherencyFilter->GetOutput());
FilterType::Pointer filter = FilterType::New(); FilterType::Pointer filter = FilterType::New();
filter->SetInput(perBandMeanFilter->GetOutput()); filter->SetInput(perBandMeanFilter->GetOutput());
filter->SetNumberOfThreads(1);
ExtractType::Pointer extract = ExtractType::New(); ExtractType::Pointer extract = ExtractType::New();
extract->SetInput(filter->GetOutput()); extract->SetInput(filter->GetOutput());
...@@ -102,6 +103,8 @@ int otbReciprocalHAlphaImageFilter(int argc, char * argv[]) ...@@ -102,6 +103,8 @@ int otbReciprocalHAlphaImageFilter(int argc, char * argv[])
writer->SetFileName(outputFilename); writer->SetFileName(outputFilename);
writer->SetInput(extract->GetOutput()); writer->SetInput(extract->GetOutput());
writer->SetNumberOfThreads(1);
writer->SetNumberOfStreamDivisions(1);
writer->Update(); writer->Update();
return EXIT_SUCCESS; return EXIT_SUCCESS;
......
Testing/Code/SARPolarimetry/otbSARPolarimetryTests2.cxx
+ 0
− 1
View file @ de72c475
...@@ -42,7 +42,6 @@ void RegisterTests() ...@@ -42,7 +42,6 @@ void RegisterTests()
REGISTER_TEST(otbMuellerToPolarisationDegreeAndPowerImageFilter); REGISTER_TEST(otbMuellerToPolarisationDegreeAndPowerImageFilter);
REGISTER_TEST(otbMuellerToCovarianceImageFilterNew); REGISTER_TEST(otbMuellerToCovarianceImageFilterNew);
REGISTER_TEST(otbMuellerToCovarianceImageFilter); REGISTER_TEST(otbMuellerToCovarianceImageFilter);
REGISTER_TEST(otbHermitianEigenAnalysisTest);
REGISTER_TEST(otbPolarimetricDataNew); REGISTER_TEST(otbPolarimetricDataNew);
REGISTER_TEST(otbPolarimetricDataTest); REGISTER_TEST(otbPolarimetricDataTest);
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment