From de72c475903cd16d27e838e83a77944464890f08 Mon Sep 17 00:00:00 2001
From: Cyrille Valladeau <cyrille.valladeau@c-s.fr>
Date: Mon, 21 Mar 2011 19:48:32 +0100
Subject: [PATCH] ENH : supress Hermitian Analysis class
---
.../otbHermitianEigenAnalysis.h | 342 ------
.../otbHermitianEigenAnalysis.txx | 983 ------------------
.../otbReciprocalHAlphaImageFilter.h | 157 ++-
Testing/Code/SARPolarimetry/CMakeLists.txt | 5 -
.../otbHermitianEigenAnalysisTest.cxx | 41 +
.../otbReciprocalHAlphaImageFilter.cxx | 5 +-
.../otbSARPolarimetryTests2.cxx | 1 -
7 files changed, 145 insertions(+), 1389 deletions(-)
delete mode 100644 Code/SARPolarimetry/otbHermitianEigenAnalysis.h
delete mode 100644 Code/SARPolarimetry/otbHermitianEigenAnalysis.txx
diff --git a/Code/SARPolarimetry/otbHermitianEigenAnalysis.h b/Code/SARPolarimetry/otbHermitianEigenAnalysis.h
deleted file mode 100644
index 16001824ae..0000000000
--- a/Code/SARPolarimetry/otbHermitianEigenAnalysis.h
+++ /dev/null
@@ -1,342 +0,0 @@
-/*=========================================================================
-
- 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
-
diff --git a/Code/SARPolarimetry/otbHermitianEigenAnalysis.txx b/Code/SARPolarimetry/otbHermitianEigenAnalysis.txx
deleted file mode 100644
index 04d32590c9..0000000000
--- a/Code/SARPolarimetry/otbHermitianEigenAnalysis.txx
+++ /dev/null
@@ -1,983 +0,0 @@
-/*=========================================================================
-
- 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_txx
-#define __otbHermitianEigenAnalysis_txx
-
-#include "otbHermitianEigenAnalysis.h"
-#include "vnl/vnl_math.h"
-
-namespace otb
-{
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-unsigned int
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ComputeEigenValues(const TMatrix & A,
- TVector & D) const
-{
- double *workArea1 = new double[ m_Dimension ];
-
- // Copy the input matrix
- double *inputMatrix = new double [ m_Dimension * m_Dimension ];
-
- unsigned int k = 0;
-
- for( unsigned int row=0; row < m_Dimension; row++ )
- {
- for( unsigned int col=0; col < m_Dimension; col++ )
- {
- inputMatrix[k++] = A(row, col);
- }
- }
-
- ReduceToTridiagonalMatrix( inputMatrix, D, workArea1, workArea1 );
- const unsigned int eigenErrIndex =
- ComputeEigenValuesUsingQL( D, workArea1 );
-
- delete[] workArea1;
- delete[] inputMatrix;
-
- return eigenErrIndex; //index of eigen value that could not be computed
-}
-
-
-/*******************************************************************************
-Routine : diagonalisation
-Authors : Eric POTTIER, Laurent FERRO-FAMIL
-Creation : 01/2002
-Update :
-*-------------------------------------------------------------------------------
-Description : Computes the eigenvectors and eigenvalues of a N*N hermitian
-matrix (with N < 10)
-*-------------------------------------------------------------------------------
-Inputs arguments :
-MatrixDim : Dimension of the Hermitian Matrix (N)
-HermitianMatrix : N*N*2 complex hermitian matrix
-Returned values :
-EigenVect : N*N*2 complex eigenvector matrix
-EigenVal : N elements eigenvalue real vector
-*******************************************************************************/
-
-static void HermitianDiagonalisation( float HM[3][3][2], float EigenVect[3][3][2], float EigenVal[3])
-{
-
- double a[10][10][2], v[10][10][2], d[10], b[10], z[10];
- double w[2], s[2], c[2], titi[2], gc[2], hc[2];
- double sm, tresh, x, toto, e, f, g, h, r, d1, d2;
- int n, p, q, ii, i, j, k;
-
- //n = MatrixDim;
- n=3;
-
-
-
- for (i = 1; i < n + 1; i++)
- {
- for (j = 1; j < n + 1; j++)
- {
- a[i][j][0] = HM[i - 1][j - 1][0];
- a[i][j][1] = HM[i - 1][j - 1][1];
- v[i][j][0] = 0.;
- v[i][j][1] = 0.;
- }
- v[i][i][0] = 1.;
- v[i][i][1] = 0.;
- }
-
- for (p = 1; p < n + 1; p++)
- {
- d[p] = a[p][p][0];
- b[p] = d[p];
- z[p] = 0.;
- }
-
- for (ii = 1; ii < 1000 * n * n; ii++)
- {
- sm = 0.;
- for (p = 1; p < n; p++)
- {
- for (q = p + 1; q < n + 1; q++)
- {
- sm = sm + 2. * sqrt(a[p][q][0] * a[p][q][0] + a[p][q][1] * a[p][q][1]);
- }
- }
- sm = sm / (n * (n - 1));
- if (sm < 1.E-16)
- {
- goto Sortie;
- }
- tresh = 1.E-17;
- if (ii < 4)
- {
- tresh = (long) 0.2 *sm / (n * n);
- }
- x = -1.E-15;
- for (i = 1; i < n; i++)
- {
- for (j = i + 1; j < n + 1; j++)
- {
- toto = sqrt(a[i][j][0] * a[i][j][0] + a[i][j][1] * a[i][j][1]);
- if (x < toto)
- {
- x = toto;
- p = i;
- q = j;
- }
- }
- }
- toto = sqrt(a[p][q][0] * a[p][q][0] + a[p][q][1] * a[p][q][1]);
- if (toto > tresh)
- {
- e = d[p] - d[q];
- w[0] = a[p][q][0];
- w[1] = a[p][q][1];
- g = sqrt(w[0] * w[0] + w[1] * w[1]);
- g = g * g;
- f = sqrt(e * e + 4. * g);
- d1 = e + f;
- d2 = e - f;
- if (fabs(d2) > fabs(d1))
- {
- d1 = d2;
- }
- r = fabs(d1) / sqrt(d1 * d1 + 4. * g);
- s[0] = r;
- s[1] = 0.;
- titi[0] = 2. * r / d1;
- titi[1] = 0.;
- c[0] = titi[0] * w[0] - titi[1] * w[1];
- c[1] = titi[0] * w[1] + titi[1] * w[0];
- r = sqrt(s[0] * s[0] + s[1] * s[1]);
- r = r * r;
- h = (d1 / 2. + 2. * g / d1) * r;
- d[p] = d[p] - h;
- z[p] = z[p] - h;
- d[q] = d[q] + h;
- z[q] = z[q] + h;
- a[p][q][0] = 0.;
- a[p][q][1] = 0.;
-
- for (j = 1; j < p; j++)
- {
- gc[0] = a[j][p][0];
- gc[1] = a[j][p][1];
- hc[0] = a[j][q][0];
- hc[1] = a[j][q][1];
- a[j][p][0] = c[0] * gc[0] - c[1] * gc[1] - s[0] * hc[0] - s[1] * hc[1];
- a[j][p][1] = c[0] * gc[1] + c[1] * gc[0] - s[0] * hc[1] + s[1] * hc[0];
- a[j][q][0] = s[0] * gc[0] - s[1] * gc[1] + c[0] * hc[0] + c[1] * hc[1];
- a[j][q][1] = s[0] * gc[1] + s[1] * gc[0] + c[0] * hc[1] - c[1] * hc[0];
- }
- for (j = p + 1; j < q; j++)
- {
- gc[0] = a[p][j][0];
- gc[1] = a[p][j][1];
- hc[0] = a[j][q][0];
- hc[1] = a[j][q][1];
- a[p][j][0] = c[0] * gc[0] + c[1] * gc[1] - s[0] * hc[0] - s[1] * hc[1];
- a[p][j][1] = c[0] * gc[1] - c[1] * gc[0] + s[0] * hc[1] - s[1] * hc[0];
- a[j][q][0] = s[0] * gc[0] + s[1] * gc[1] + c[0] * hc[0] + c[1] * hc[1];
- a[j][q][1] = -s[0] * gc[1] + s[1] * gc[0] + c[0] * hc[1] - c[1] * hc[0];
- }
- for (j = q + 1; j < n + 1; j++)
- {
- gc[0] = a[p][j][0];
- gc[1] = a[p][j][1];
- hc[0] = a[q][j][0];
- hc[1] = a[q][j][1];
- a[p][j][0] = c[0] * gc[0] + c[1] * gc[1] - s[0] * hc[0] + s[1] * hc[1];
- a[p][j][1] = c[0] * gc[1] - c[1] * gc[0] - s[0] * hc[1] - s[1] * hc[0];
- a[q][j][0] = s[0] * gc[0] + s[1] * gc[1] + c[0] * hc[0] - c[1] * hc[1];
- a[q][j][1] = s[0] * gc[1] - s[1] * gc[0] + c[0] * hc[1] + c[1] * hc[0];
- }
- for (j = 1; j < n + 1; j++)
- {
- gc[0] = v[j][p][0];
- gc[1] = v[j][p][1];
- hc[0] = v[j][q][0];
- hc[1] = v[j][q][1];
- v[j][p][0] = c[0] * gc[0] - c[1] * gc[1] - s[0] * hc[0] - s[1] * hc[1];
- v[j][p][1] = c[0] * gc[1] + c[1] * gc[0] - s[0] * hc[1] + s[1] * hc[0];
- v[j][q][0] = s[0] * gc[0] - s[1] * gc[1] + c[0] * hc[0] + c[1] * hc[1];
- v[j][q][1] = s[0] * gc[1] + s[1] * gc[0] + c[0] * hc[1] - c[1] * hc[0];
- }
- }
- }
-
- Sortie:
-
- for (k = 1; k < n + 1; k++)
- {
- d[k] = 0;
- for (i = 1; i < n + 1; i++)
- {
- for (j = 1; j < n + 1; j++)
- {
- d[k] = d[k] + v[i][k][0] * (HM[i - 1][j - 1][0] * v[j][k][0] - HM[i - 1][j - 1][1] * v[j][k][1]);
- d[k] = d[k] + v[i][k][1] * (HM[i - 1][j - 1][0] * v[j][k][1] + HM[i - 1][j - 1][1] * v[j][k][0]);
- }
- }
- }
-
- for (i = 1; i < n + 1; i++)
- {
- for (j = i + 1; j < n + 1; j++)
- {
- if (d[j] > d[i])
- {
- x = d[i];
- d[i] = d[j];
- d[j] = x;
- for (k = 1; k < n + 1; k++)
- {
- c[0] = v[k][i][0];
- c[1] = v[k][i][1];
- v[k][i][0] = v[k][j][0];
- v[k][i][1] = v[k][j][1];
- v[k][j][0] = c[0];
- v[k][j][1] = c[1];
- }
- }
- }
- }
-
- for (i = 0; i < n; i++)
- {
- EigenVal[i] = d[i + 1];
- for (j = 0; j < n; j++)
- {
- EigenVect[i][j][0] = v[i + 1][j + 1][0];
- EigenVect[i][j][1] = v[i + 1][j + 1][1];
- }
- }
-
-}
-
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-unsigned int
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ComputeEigenValuesAndVectors(
- const TMatrix & T,
- TVector & EigenValues,
- TEigenMatrix & EigenVectors ) const
-{
- // double *workArea1 = new double[ m_Dimension ];
-// double *workArea2 = new double [ m_Dimension * m_Dimension ];
-
- // Copy the input matrix
-
- unsigned int k = 0;
-
- float HM[3][3][2];
-
-
-// Passage de la matrice T T1....T9 une matrice 3*3*2 compatible avec la methode HermitianDiagonalisation
-
-
- HM[0][0][0]=T[0]; HM[0][0][1]=0.;
- HM[0][1][0]=T[3]; HM[0][1][1]=-T[4];
- HM[0][2][0]=T[5]; HM[0][2][1]=-T[6];
- HM[1][0][0]=T[3]; HM[1][0][1]=T[4];
- HM[1][1][0]=T[1]; HM[1][1][1]=0.;
- HM[1][2][0]=T[7]; HM[1][2][1]=-T[8];
- HM[2][0][0]=T[5]; HM[2][0][1]=T[6];
- HM[2][1][0]=T[7]; HM[2][1][1]=T[8];
- HM[2][2][0]=T[2]; HM[2][2][1]=0.;
-
-
-
-
- //Acrobatie pour convectir EigenVectors en un float *** et EigenValues en float*
-
- float eigenVect[3][3][2];
- float eigenVal[3];
-
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++)
- {
- eigenVect[i][j][0]=0.;
- eigenVect[i][j][1]=0.;
- }
-
- for (int i=0; i<3; i++)
- eigenVal[i]=0.;
-
-
-
- HermitianDiagonalisation( HM, eigenVect, eigenVal);
-
-
- //std::cout << HM[0][0][0] << " " << HM[0][1][0]<<"+j"<<HM[0][1][1] << " " << HM[0][2][0]<<"+j"<<HM[0][2][1]<<std::endl;
- //std::cout << HM[1][0][0]<<"+j"<<HM[1][0][1] << " " << HM[1][1][0] << " " << HM[1][2][0]<<"+j"<<HM[1][2][1]<<std::endl;
- //std::cout << HM[2][0][0]<<"+j"<<HM[2][0][1] << " " << HM[2][1][0]<<"+j"<< HM[2][1][1] << " " <<HM[2][2][0]<<std::endl;
- //std::cout << "Valeurs propres : " << eigenVal[0] << "...." << eigenVal[1] << "...." << eigenVal[2] << std::endl<< std::endl;
-
-
- // Recuperation des sorties
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++)
- {
- EigenVectors[i][2*j]=eigenVect[i][j][0];
- EigenVectors[i][2*j+1]=eigenVect[i][j][1];
- }
-
-
- for (int i=0; i<3; i++)
- EigenValues[i]=eigenVal[i];
-
-
-
-
-
- int eigenErrIndex=0;
- return eigenErrIndex; //index of eigen value that could not be computed
-}
-
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-void
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ReduceToTridiagonalMatrix(double * a, VectorType &d,
- double *e, double *e2) const
-{
- double d__1;
-
- /* Local variables */
- double f, g, h;
- int i, j, k, l;
- double scale;
-
- for (i = 0; i < static_cast< int >(m_Order); ++i)
- {
- d[i] = a[m_Order-1 + i * m_Dimension];
- a[m_Order-1 + i * m_Dimension] = a[i + i * m_Dimension];
- }
-
-
- for (i = m_Order-1; i >= 0; --i)
- {
- l = i - 1;
- h = 0.;
- scale = 0.;
-
- /* .......... scale row (algol tol then not needed) .......... */
- for (k = 0; k <= l; ++k)
- {
- scale += vnl_math_abs(d[k]);
- }
- if (scale == 0.)
- {
- for (j = 0; j <= l; ++j)
- {
- d[j] = a[l + j * m_Dimension];
- a[l + j * m_Dimension] = a[i + j * m_Dimension];
- a[i + j * m_Dimension] = 0.;
- }
- e[i] = 0.;
- e2[i] = 0.;
- continue;
- }
- for (k = 0; k <= l; ++k)
- {
- d[k] /= scale;
- h += d[k] * d[k];
- }
-
- e2[i] = scale * scale * h;
- f = d[l];
- d__1 = vcl_sqrt(h);
- g = (-1.0) * vnl_math_sgn0(f) * vnl_math_abs(d__1);
- e[i] = scale * g;
- h -= f * g;
- d[l] = f - g;
- if (l != 0)
- {
-
- /* .......... form a*u .......... */
- for (j = 0; j <= l; ++j)
- {
- e[j] = 0.;
- }
-
- for (j = 0; j <= l; ++j)
- {
- f = d[j];
- g = e[j] + a[j + j * m_Dimension] * f;
-
- for (k = j+1; k <= l; ++k)
- {
- g += a[k + j * m_Dimension] * d[k];
- e[k] += a[k + j * m_Dimension] * f;
- }
- e[j] = g;
- }
-
- /* .......... form p .......... */
- f = 0.;
-
- for (j = 0; j <= l; ++j)
- {
- e[j] /= h;
- f += e[j] * d[j];
- }
-
- h = f / (h + h);
- /* .......... form q .......... */
- for (j = 0; j <= l; ++j)
- {
- e[j] -= h * d[j];
- }
-
- /* .......... form reduced a .......... */
- for (j = 0; j <= l; ++j)
- {
- f = d[j];
- g = e[j];
-
- for (k = j; k <= l; ++k)
- {
- a[k + j * m_Dimension] = a[k + j * m_Dimension] - f * e[k] - g * d[k];
- }
- }
- }
-
- for (j = 0; j <= l; ++j)
- {
- f = d[j];
- d[j] = a[l + j * m_Dimension];
- a[l + j * m_Dimension] = a[i + j * m_Dimension];
- a[i + j * m_Dimension] = f * scale;
- }
- }
-}
-
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-void
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ReduceToTridiagonalMatrixAndGetTransformation(double * a, VectorType &d,
- double *e, double *z) const
-{
- double d__1;
-
- /* Local variables */
- double f, g, h;
- unsigned int i, j, k, l;
- double scale, hh;
-
- for (i = 0; i < m_Order; ++i)
- {
- for (j = i; j < m_Order; ++j)
- {
- z[j + i * m_Dimension] = a[j + i * m_Dimension];
- }
- d[i] = a[m_Order -1 + i * m_Dimension];
- }
-
- for (i = m_Order-1; i > 0; --i)
- {
- l = i - 1;
- h = 0.0;
- scale = 0.0;
-
- /* .......... scale row (algol tol then not needed) .......... */
- for (k = 0; k <= l; ++k)
- {
- scale += vnl_math_abs(d[k]);
- }
- if (scale == 0.0)
- {
- e[i] = d[l];
-
- for (j = 0; j <= l; ++j)
- {
- d[j] = z[l + j * m_Dimension];
- z[i + j * m_Dimension] = 0.0;
- z[j + i * m_Dimension] = 0.0;
- }
- }
- else
- {
- for (k = 0; k <= l; ++k)
- {
- d[k] /= scale;
- h += d[k] * d[k];
- }
-
- f = d[l];
- d__1 = vcl_sqrt(h);
- g = (-1.0) * vnl_math_sgn0(f) * vnl_math_abs(d__1);
- e[i] = scale * g;
- h -= f * g;
- d[l] = f - g;
-
- /* .......... form a*u .......... */
- for (j = 0; j <= l; ++j)
- {
- e[j] = 0.0;
- }
-
- for (j = 0; j <= l; ++j)
- {
- f = d[j];
- z[j + i * m_Dimension] = f;
- g = e[j] + z[j + j * m_Dimension] * f;
-
- for (k = j+1; k <= l; ++k)
- {
- g += z[k + j * m_Dimension] * d[k];
- e[k] += z[k + j * m_Dimension] * f;
- }
-
- e[j] = g;
- }
-
- /* .......... form p .......... */
- f = 0.0;
-
- for (j = 0; j <= l; ++j)
- {
- e[j] /= h;
- f += e[j] * d[j];
- }
-
- hh = f / (h + h);
-
- /* .......... form q .......... */
- for (j = 0; j <= l; ++j)
- {
- e[j] -= hh * d[j];
- }
-
- /* .......... form reduced a .......... */
- for (j = 0; j <= l; ++j)
- {
- f = d[j];
- g = e[j];
-
- for (k = j; k <= l; ++k)
- {
- z[k + j * m_Dimension] = z[k + j * m_Dimension] - f * e[k] - g * d[k];
- }
-
- d[j] = z[l + j * m_Dimension];
- z[i + j * m_Dimension] = 0.0;
- }
- }
-
- d[i] = h;
- }
-
- /* .......... accumulation of transformation matrices .......... */
- for (i = 1; i < m_Order; ++i)
- {
- l = i - 1;
- z[m_Order-1 + l * m_Dimension] = z[l + l * m_Dimension];
- z[l + l * m_Dimension] = 1.0;
- h = d[i];
- if (h != 0.0)
- {
- for (k = 0; k <= l; ++k)
- {
- d[k] = z[k + i * m_Dimension] / h;
- }
-
- for (j = 0; j <= l; ++j)
- {
- g = 0.0;
-
- for (k = 0; k <= l; ++k)
- {
- g += z[k + i * m_Dimension] * z[k + j * m_Dimension];
- }
-
- for (k = 0; k <= l; ++k)
- {
- z[k + j * m_Dimension] -= g * d[k];
- }
- }
- }
-
- for (k = 0; k <= l; ++k)
- {
- z[k + i * m_Dimension] = 0.0;
- }
-
- }
-
- for (i = 0; i < m_Order; ++i)
- {
- d[i] = z[m_Order-1 + i * m_Dimension];
- z[m_Order-1 + i * m_Dimension] = 0.0;
- }
-
- z[m_Order-1 + (m_Order-1) * m_Dimension] = 1.0;
- e[0] = 0.0;
-
-}
-
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-unsigned int
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ComputeEigenValuesUsingQL(VectorType &d, double *e) const
-{
-
- const double c_b10 = 1.0;
-
- /* Local variables */
- double c, f, g, h;
- unsigned int i, j, l, m;
- double p, r, s, c2, c3=0.0;
- double s2=0.0;
- double dl1, el1;
- double tst1, tst2;
-
- unsigned int ierr = 0;
- if (m_Order == 1) {
- return 1;
- }
-
- for (i = 1; i < m_Order; ++i) {
- e[i-1] = e[i];
- }
-
- f = 0.;
- tst1 = 0.;
- e[m_Order-1] = 0.;
-
- for (l = 0; l < m_Order; ++l)
- {
- j = 0;
- h = vnl_math_abs(d[l]) + vnl_math_abs(e[l]);
- if (tst1 < h)
- {
- tst1 = h;
- }
- /* .......... look for small sub-diagonal element .......... */
- for (m = l; m < m_Order; ++m)
- {
- tst2 = tst1 + vnl_math_abs(e[m]);
- if (tst2 == tst1)
- {
- break;
- }
- /* .......... e(n) is always zero, so there is no exit */
- /* through the bottom of the loop .......... */
- }
-
- if (m != l)
- {
-
- do
- {
- if (j == 30)
- {
- /* .......... set error -- no convergence to an */
- /* eigenvalue after 30 iterations .......... */
- ierr = l+1;
- return ierr;
- }
- ++j;
- /* .......... form shift .......... */
- g = d[l];
- p = (d[l+1] - g) / (e[l] * 2.);
- r = vnl_math_hypot(p, c_b10);
- d[l] = e[l] / (p + vnl_math_sgn0(p) * vnl_math_abs(r));
- d[l+1] = e[l] * (p + vnl_math_sgn0(p) * vnl_math_abs(r));
- dl1 = d[l+1];
- h = g - d[l];
-
- for (i = l+2; i < m_Order; ++i)
- {
- d[i] -= h;
- }
-
- f += h;
- /* .......... ql transformation .......... */
- p = d[m];
- c = 1.;
- c2 = c;
- el1 = e[l+1];
- s = 0.;
- for (i = m-1; i >= l; --i)
- {
- c3 = c2;
- c2 = c;
- s2 = s;
- g = c * e[i];
- h = c * p;
- r = vnl_math_hypot(p, e[i]);
- e[i+1] = s * r;
- s = e[i] / r;
- c = p / r;
- p = c * d[i] - s * g;
- d[i+1] = h + s * (c * g + s * d[i]);
- if( i == l )
- {
- break;
- }
- }
-
- p = -s * s2 * c3 * el1 * e[l] / dl1;
- e[l] = s * p;
- d[l] = c * p;
- tst2 = tst1 + vnl_math_abs(e[l]);
- } while (tst2 > tst1);
- }
-
- p = d[l] + f;
-
- if( m_OrderEigenValues == OrderByValue )
- {
- // Order by value
- for (i = l; i > 0; --i)
- {
- if (p >= d[i-1])
- break;
- d[i] = d[i-1];
- }
- d[i] = p;
- }
- else if( m_OrderEigenValues == OrderByMagnitude )
- {
- // Order by magnitude.. make eigen values positive
- for (i = l; i > 0; --i)
- {
- if (vnl_math_abs(p) >= vnl_math_abs(d[i-1]))
- break;
- d[i] = vnl_math_abs(d[i-1]);
- }
- d[i] = vnl_math_abs(p);
- }
- else
- {
- d[l] = p;
- }
- }
-
- return ierr; //ierr'th eigen value that couldn't be computed
-
-}
-
-
-template< class TMatrix, class TVector, class TEigenMatrix >
-unsigned int
-HermitianEigenAnalysis< TMatrix, TVector, TEigenMatrix >::
-ComputeEigenValuesAndVectorsUsingQL(VectorType &d, double *e, double *z) const
-{
-
- const double c_b10 = 1.0;
-
- /* Local variables */
- double c, f, g, h;
- unsigned int i, j, k, l, m;
- double p, r, s, c2, c3=0.0;
- double s2=0.0;
- double dl1, el1;
- double tst1, tst2;
-
- unsigned int ierr = 0;
- if (m_Order == 1)
- {
- return 1;
- }
-
- for (i = 1; i < m_Order; ++i)
- {
- e[i - 1] = e[i];
- }
-
- f = 0.0;
- tst1 = 0.;
- e[m_Order-1] = 0.;
-
- for (l = 0; l < m_Order; ++l)
- {
- j = 0;
- h = vnl_math_abs(d[l]) + vnl_math_abs(e[l]);
- if (tst1 < h)
- {
- tst1 = h;
- }
-
- /* .......... look for small sub-diagonal element .......... */
- for (m = l; m < m_Order; ++m)
- {
- tst2 = tst1 + vnl_math_abs(e[m]);
- if (tst2 == tst1)
- {
- break;
- }
-
- /* .......... e(n) is always zero, so there is no exit */
- /* through the bottom of the loop .......... */
- }
-
- if (m != l)
- {
- do
- {
-
- if (j == 1000)
- {
- /* .......... set error -- no convergence to an */
- /* eigenvalue after 1000 iterations .......... */
- ierr = l+1;
- return ierr;
- }
- ++j;
- /* .......... form shift .......... */
- g = d[l];
- p = (d[l+1] - g) / (e[l] * 2.);
- r = vnl_math_hypot(p, c_b10);
- d[l] = e[l] / (p + vnl_math_sgn0(p) * vnl_math_abs(r));
- d[l+1] = e[l] * (p + vnl_math_sgn0(p) * vnl_math_abs(r));
- dl1 = d[l+1];
- h = g - d[l];
-
- for (i = l+2; i < m_Order; ++i)
- {
- d[i] -= h;
- }
-
- f += h;
- /* .......... ql transformation .......... */
- p = d[m];
- c = 1.0;
- c2 = c;
- el1 = e[l+1];
- s = 0.;
-
- for (i = m-1; i >= l; --i)
- {
- c3 = c2;
- c2 = c;
- s2 = s;
- g = c * e[i];
- h = c * p;
- r = vnl_math_hypot(p, e[i]);
- e[i + 1] = s * r;
- s = e[i] / r;
- c = p / r;
- p = c * d[i] - s * g;
- d[i + 1] = h + s * (c * g + s * d[i]);
-
- /* .......... form vector .......... */
- for (k = 0; k < m_Order; ++k)
- {
- h = z[k + (i + 1) * m_Dimension];
- z[k + (i + 1) * m_Dimension] = s * z[k + i * m_Dimension] + c * h;
- z[k + i * m_Dimension] = c * z[k + i * m_Dimension] - s * h;
- }
- if( i == l )
- {
- break;
- }
- }
-
- p = -s * s2 * c3 * el1 * e[l] / dl1;
- e[l] = s * p;
- d[l] = c * p;
- tst2 = tst1 + vnl_math_abs(e[l]);
- } while (tst2 > tst1);
-
- }
-
- d[l] += f;
- }
-
- /* .......... order eigenvalues and eigenvectors .......... */
- if( m_OrderEigenValues == OrderByValue )
- {
- // Order by value
- for (i = 0; i < m_Order-1; ++i)
- {
- k = i;
- p = d[i];
-
- for (j = i+1; j < m_Order; ++j)
- {
- if (d[j] >= p)
- {
- continue;
- }
- k = j;
- p = d[j];
- }
-
- if (k == i)
- {
- continue;
- }
- d[k] = d[i];
- d[i] = p;
-
- for (j = 0; j < m_Order; ++j)
- {
- p = z[j + i * m_Dimension];
- z[j + i * m_Dimension] = z[j + k * m_Dimension];
- z[j + k * m_Dimension] = p;
- }
- }
- }
- else if( m_OrderEigenValues == OrderByMagnitude )
- {
- // Order by magnitude
- for (i = 0; i < m_Order-1; ++i)
- {
- k = i;
- p = d[i];
-
- for (j = i+1; j < m_Order; ++j)
- {
- if (vnl_math_abs(d[j]) >= vnl_math_abs(p))
- {
- continue;
- }
- k = j;
- p = d[j];
- }
-
- if (k == i)
- {
- continue;
- }
- d[k] = vnl_math_abs(d[i]);
- d[i] = vnl_math_abs(p);
-
- for (j = 0; j < m_Order; ++j)
- {
- p = z[j + i * m_Dimension];
- z[j + i * m_Dimension] = z[j + k * m_Dimension];
- z[j + k * m_Dimension] = p;
- }
- }
- }
-
-
- return ierr;
-}
-
-
-} // end namespace otb
-
-#endif
-
diff --git a/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.h b/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.h
index 18292eca3b..5a18300e3d 100644
--- a/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.h
+++ b/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.h
@@ -20,9 +20,9 @@
#define __ReciprocalHAlphaImageFilter_h
#include "otbUnaryFunctorImageFilter.h"
-#include "otbHermitianEigenAnalysis.h"
#include "itkVector.h"
#include "otbMath.h"
+#include "vnl/algo/vnl_complex_eigensystem.h"
namespace otb
{
@@ -44,60 +44,108 @@ template< class TInput, class TOutput>
class ReciprocalHAlphaFunctor
{
public:
- typedef typename std::complex <double> ComplexType;
- typedef typename TOutput::ValueType OutputValueType;
-
- /** CoherencyMatrix type **/
- typedef itk::Vector<double, 9> CoherencyMatrixType;
-
- /** 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;
-
-
+ typedef typename std::complex<double> ComplexType;
+ typedef vnl_matrix<ComplexType> VNLMatrixType;
+ typedef vnl_vector<ComplexType> VNLVectorType;
+ typedef vnl_vector<double> VNLDoubleVectorType;
+ typedef typename TOutput::ValueType OutputValueType;
+
+
inline TOutput operator()( const TInput & Coherency ) const
{
TOutput result;
result.SetSize(m_NumberOfComponentsPerPixel);
- CoherencyMatrixType T;
- EigenvalueType eigenValues;
- EigenMatrixType eigenVectors;
- HermitianAnalysisType HermitianAnalysis;
-
- T[0] = static_cast<double>(Coherency[0].real());
- T[1] = static_cast<double>(Coherency[3].real());
- T[2] = static_cast<double>(Coherency[5].real());
- T[3] = static_cast<double>(Coherency[1].real());
- T[4] = static_cast<double>(Coherency[1].imag());
- T[5] = static_cast<double>(Coherency[2].real());
- T[6] = static_cast<double>(Coherency[2].imag());
- T[7] = static_cast<double>(Coherency[4].real());
- T[8] = static_cast<double>(Coherency[4].imag());
- HermitianAnalysis.ComputeEigenValuesAndVectors(T, eigenValues, eigenVectors);
-
+ double T0 = static_cast<double>(Coherency[0].real());
+ double T1 = static_cast<double>(Coherency[3].real());
+ double T2 = static_cast<double>(Coherency[5].real());
+ double T3 = static_cast<double>(Coherency[1].real());
+ double T4 = static_cast<double>(Coherency[1].imag());
+ double T5 = static_cast<double>(Coherency[2].real());
+ double T6 = static_cast<double>(Coherency[2].imag());
+ double T7 = static_cast<double>(Coherency[4].real());
+ double T8 = static_cast<double>(Coherency[4].imag());
+
+ VNLMatrixType vnlMat(3, 3, 0.);
+ vnlMat[0][0] = ComplexType(T0, 0.);
+ vnlMat[0][1] = std::conj(ComplexType(Coherency[1]));
+ vnlMat[0][2] = std::conj(ComplexType(Coherency[2]));
+ 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
double totalEigenValues(0.0);
double p[3];
double entropy;
double alpha;
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)
{
totalEigenValues = m_Epsilon;
@@ -105,12 +153,7 @@ public:
for (unsigned int k = 0; k < 3; k++)
{
- if (eigenValues[k] <0.)
- {
- eigenValues[k] = 0.;
- }
-
- p[k] = static_cast<double>(eigenValues[k]) / totalEigenValues;
+ p[k] = std::max(sortedRealEigenValues[k], 0.) / totalEigenValues;
}
if ( (p[0] < m_Epsilon) || (p[1] < m_Epsilon) || (p[2] < m_Epsilon) )
@@ -129,7 +172,7 @@ public:
for(unsigned int k = 0; k < 3; k++)
{
- p[k] = eigenValues[k] / totalEigenValues;
+ p[k] = sortedRealEigenValues[k] / totalEigenValues;
if (p[k] < 0.)
{
@@ -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;
- 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;
- 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;
alpha=p[0]*a0 + p[1]*a1 + p[2]*a2;
@@ -159,11 +202,11 @@ public:
}
// 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[1] = alpha;
- result[2] = anisotropy;
+ result[0] = static_cast<OutputValueType>(entropy);
+ result[1] = static_cast<OutputValueType>(alpha);
+ result[2] = static_cast<OutputValueType>(anisotropy);
return result;
}
diff --git a/Testing/Code/SARPolarimetry/CMakeLists.txt b/Testing/Code/SARPolarimetry/CMakeLists.txt
index 654b91338e..c291daf45d 100644
--- a/Testing/Code/SARPolarimetry/CMakeLists.txt
+++ b/Testing/Code/SARPolarimetry/CMakeLists.txt
@@ -373,10 +373,6 @@ ADD_TEST(saTvMuellerToCovarianceImageFilter ${SARPOLARIMETRY_TESTS2}
${BASELINE}/saTvSinclairImageFilter_SinclairToMueller.tif
${TEMP}/saTvMuellerToMLCImageFilter.tif
)
-# Hermitian eigen analysis class
-ADD_TEST(saTvHermitianEigenAnalysisTest ${SARPOLARIMETRY_TESTS2}
- otbHermitianEigenAnalysisTest
-)
# Polarimetric Data
ADD_TEST(saTuPolarimetricDataNew ${SARPOLARIMETRY_TESTS2}
@@ -430,7 +426,6 @@ otbMuellerToPolarisationDegreeAndPowerImageFilterNew.cxx
otbMuellerToPolarisationDegreeAndPowerImageFilter.cxx
otbMuellerToCovarianceImageFilterNew.cxx
otbMuellerToCovarianceImageFilter.cxx
-otbHermitianEigenAnalysisTest.cxx
otbPolarimetricData.cxx
)
diff --git a/Testing/Code/SARPolarimetry/otbHermitianEigenAnalysisTest.cxx b/Testing/Code/SARPolarimetry/otbHermitianEigenAnalysisTest.cxx
index d21ca4c3dd..62a756f674 100644
--- a/Testing/Code/SARPolarimetry/otbHermitianEigenAnalysisTest.cxx
+++ b/Testing/Code/SARPolarimetry/otbHermitianEigenAnalysisTest.cxx
@@ -26,6 +26,8 @@
#include "otbVectorImage.h"
#include "otbHermitianEigenAnalysis.h"
#include "itkVariableLengthVector.h"
+#include "vnl/algo/vnl_complex_eigensystem.h"
+#include <complex>
int otbHermitianEigenAnalysisTest(int argc, char * argv[])
@@ -37,6 +39,7 @@ int otbHermitianEigenAnalysisTest(int argc, char * argv[])
typedef otb::HermitianEigenAnalysis<MatrixType, EigenvalueType, EigenMatrixType> FilterType;
EigenMatrixType resEigVal;
+
itk::Vector<float, 6> temp;
temp[0] = 0.162793;
temp[1] = -0.432753;
@@ -115,5 +118,43 @@ int otbHermitianEigenAnalysisTest(int argc, char * argv[])
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;
}
diff --git a/Testing/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.cxx b/Testing/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.cxx
index ed970d3191..79c312c5cb 100644
--- a/Testing/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.cxx
+++ b/Testing/Code/SARPolarimetry/otbReciprocalHAlphaImageFilter.cxx
@@ -89,9 +89,10 @@ int otbReciprocalHAlphaImageFilter(int argc, char * argv[])
PerBandMeanFilterType::Pointer perBandMeanFilter = PerBandMeanFilterType::New();
perBandMeanFilter->SetInput(sinclairToCoherencyFilter->GetOutput());
-
+
FilterType::Pointer filter = FilterType::New();
filter->SetInput(perBandMeanFilter->GetOutput());
+ filter->SetNumberOfThreads(1);
ExtractType::Pointer extract = ExtractType::New();
extract->SetInput(filter->GetOutput());
@@ -102,6 +103,8 @@ int otbReciprocalHAlphaImageFilter(int argc, char * argv[])
writer->SetFileName(outputFilename);
writer->SetInput(extract->GetOutput());
+ writer->SetNumberOfThreads(1);
+writer->SetNumberOfStreamDivisions(1);
writer->Update();
return EXIT_SUCCESS;
diff --git a/Testing/Code/SARPolarimetry/otbSARPolarimetryTests2.cxx b/Testing/Code/SARPolarimetry/otbSARPolarimetryTests2.cxx
index bdda5c3240..97a6220613 100644
--- a/Testing/Code/SARPolarimetry/otbSARPolarimetryTests2.cxx
+++ b/Testing/Code/SARPolarimetry/otbSARPolarimetryTests2.cxx
@@ -42,7 +42,6 @@ void RegisterTests()
REGISTER_TEST(otbMuellerToPolarisationDegreeAndPowerImageFilter);
REGISTER_TEST(otbMuellerToCovarianceImageFilterNew);
REGISTER_TEST(otbMuellerToCovarianceImageFilter);
- REGISTER_TEST(otbHermitianEigenAnalysisTest);
REGISTER_TEST(otbPolarimetricDataNew);
REGISTER_TEST(otbPolarimetricDataTest);
}
--
GitLab