From d1214e56251b76ab2f5f312d795d87cc767c5e44 Mon Sep 17 00:00:00 2001 From: Gregoire Mercier <gregoire.mercier@telecom-bretagne.eu> Date: Wed, 2 Mar 2011 15:45:34 +0100 Subject: [PATCH] TEST: test NAPCA --- Code/Hyperspectral/otbNAPCAImageFilter.txx | 18 +-- Testing/CMakeLists.txt | 12 ++ Testing/otbHyperTests1.cxx | 2 + Testing/otbNAPCAImageFilter.cxx | 155 +++++++++++++++++++++ 4 files changed, 178 insertions(+), 9 deletions(-) create mode 100644 Testing/otbNAPCAImageFilter.cxx diff --git a/Code/Hyperspectral/otbNAPCAImageFilter.txx b/Code/Hyperspectral/otbNAPCAImageFilter.txx index 53f65f6bad..eb91bc27fd 100644 --- a/Code/Hyperspectral/otbNAPCAImageFilter.txx +++ b/Code/Hyperspectral/otbNAPCAImageFilter.txx @@ -35,8 +35,8 @@ void NAPCAImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTransformation > ::GetTransformationMatrixFromCovarianceMatrix () { - InternalMatrixType An = m_NoiseCovarianceMatrix.GetVnlMatrix(); - InternalMatrixType Ax = m_CovarianceMatrix.GetVnlMatrix(); + InternalMatrixType An = this->GetNoiseCovarianceMatrix().GetVnlMatrix(); + InternalMatrixType Ax = this->GetCovarianceMatrix().GetVnlMatrix(); vnl_svd< MatrixElementType > An_solver ( An ); InternalMatrixType U_cholesky = An_solver.U(); @@ -44,24 +44,24 @@ NAPCAImageFilter< TInputImage, TOutputImage, TNoiseImageFilter, TDirectionOfTran InternalMatrixType C = U.transpose() * Ax * U; InternalMatrixType Id ( C.rows(), C.cols(), vnl_matrix_identity ); - vnl_generalized_eigensystem solver ( C, I ); + vnl_generalized_eigensystem solver ( C, Id ); InternalMatrixType transf = solver.V; transf *= U.transpose(); transf.fliplr(); - if ( m_NumberOfPrincipalComponentsRequired + if ( this->GetNumberOfPrincipalComponentsRequired() != this->GetInput()->GetNumberOfComponentsPerPixel() ) - m_TransformationMatrix = transf.get_n_rows( 0, m_NumberOfPrincipalComponentsRequired ); + this->m_TransformationMatrix = transf.get_n_rows( 0, this->GetNumberOfPrincipalComponentsRequired() ); else - m_TransformationMatrix = transf; + this->m_TransformationMatrix = transf; vnl_vector< double > valP = solver.D.diagonal(); valP.flip(); - m_EigenValues.SetSize( m_NumberOfPrincipalComponentsRequired ); - for ( unsigned int i = 0; i < m_NumberOfPrincipalComponentsRequired; i++ ) - m_EigenValues[i] = static_cast< RealType >( valP[i] ); + this->m_EigenValues.SetSize( this->GetNumberOfPrincipalComponentsRequired() ); + for ( unsigned int i = 0; i < this->GetNumberOfPrincipalComponentsRequired(); i++ ) + this->m_EigenValues[i] = static_cast< RealType >( valP[i] ); } diff --git a/Testing/CMakeLists.txt b/Testing/CMakeLists.txt index 3d35698b87..a63ba85d9c 100644 --- a/Testing/CMakeLists.txt +++ b/Testing/CMakeLists.txt @@ -26,6 +26,7 @@ SET( otbHyperTests1_SRC otbNormalizeVectorImageFilter.cxx otbPCAImageFilter.cxx otbMNFImageFilter.cxx + otbNAPCAImageFilter.cxx ) ADD_EXECUTABLE(otbHyperTests1 ${otbHyperTests1_SRC} otbHyperTests1.cxx) @@ -421,6 +422,17 @@ ADD_TEST(hyTvMNFImageFilter -inv ${TEMP}/cupriteMNFinv.hdr -out ${TEMP}/cupriteMNF.hdr) +ADD_TEST(hyTuNAPCAImageFilterNew + ${TESTEXE_DIR}/otbHyperTests1 + otbNAPCAImageFilterNewTest) + +ADD_TEST(hyTvNAPCAImageFilter + ${TESTEXE_DIR}/otbHyperTests1 + otbNAPCAImageFilterTest + -in ${DATA}/CupriteSubHsi/cupriteSubHsi.tif + -inv ${TEMP}/cupriteNAPCAinv.hdr + -out ${TEMP}/cupriteNAPCA.hdr) + ADD_TEST(hyTuEigenvalueLikelihoodMaximizationNew ${TESTEXE_DIR}/otbHyperTests1 otbEigenvalueLikelihoodMaximizationNewTest) diff --git a/Testing/otbHyperTests1.cxx b/Testing/otbHyperTests1.cxx index bdebe01a67..aa1d71282f 100644 --- a/Testing/otbHyperTests1.cxx +++ b/Testing/otbHyperTests1.cxx @@ -63,4 +63,6 @@ void RegisterTests() REGISTER_TEST(otbPCAImageFilterTest); REGISTER_TEST(otbMNFImageFilterNewTest); REGISTER_TEST(otbMNFImageFilterTest); + REGISTER_TEST(otbNAPCAImageFilterNewTest); + REGISTER_TEST(otbNAPCAImageFilterTest); } diff --git a/Testing/otbNAPCAImageFilter.cxx b/Testing/otbNAPCAImageFilter.cxx new file mode 100644 index 0000000000..4c69d68f23 --- /dev/null +++ b/Testing/otbNAPCAImageFilter.cxx @@ -0,0 +1,155 @@ +/*========================================================================= + + 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 "otbVectorImage.h" +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "otbCommandProgressUpdate.h" +#include "otbCommandLineArgumentParser.h" + +#include "otbNAPCAImageFilter.h" + +#include "otbLocalActivityVectorImageFilter.h" + +int otbNAPCAImageFilterNewTest ( int argc, char* argv[] ) +{ + const unsigned int Dimension = 2; + typedef double PixelType; + typedef otb::VectorImage< PixelType, Dimension > ImageType; + + typedef otb::LocalActivityVectorImageFilter< ImageType, ImageType > NoiseFilterType; + + typedef otb::NAPCAImageFilter< ImageType, ImageType, + NoiseFilterType, otb::Transform::FORWARD > FilterType; + FilterType::Pointer filter = FilterType::New(); + + typedef otb::NAPCAImageFilter< ImageType, ImageType, + NoiseFilterType, otb::Transform::INVERSE > InvFilterType; + InvFilterType::Pointer invFilter = InvFilterType::New(); + + return EXIT_SUCCESS; +} + +int otbNAPCAImageFilterTest ( int argc, char* argv[] ) +{ + typedef otb::CommandLineArgumentParser ParserType; + ParserType::Pointer parser = ParserType::New(); + + parser->AddInputImage(); + parser->AddOption( "--NumComponents", "Number of components to keep for output", "-n", 1, false ); + parser->AddOption( "--Inverse", "Performs also the inverse transformation (give the output name)", "-inv", 1, false ); + parser->AddOption( "--Radius", "Set the radius of the sliding window (def.1)", "-r", 2, false ); + parser->AddOption( "--Normalize", "center and reduce data before NAPCA", "-norm", 0, false ); + parser->AddOutputImage(); + + typedef otb::CommandLineArgumentParseResult ParserResultType; + ParserResultType::Pointer parseResult = ParserResultType::New(); + + try + { + parser->ParseCommandLine( argc, argv, parseResult ); + } + catch( itk::ExceptionObject & err ) + { + std::cerr << argv[0] << " applies NAPCA transformations\n"; + std::string descriptionException = err.GetDescription(); + if ( descriptionException.find("ParseCommandLine(): Help Parser") + != std::string::npos ) + return EXIT_SUCCESS; + if(descriptionException.find("ParseCommandLine(): Version Parser") + != std::string::npos ) + return EXIT_SUCCESS; + return EXIT_FAILURE; + } + + const char * inputImageName = parseResult->GetInputImage().c_str(); + const char * outputImageName = parseResult->GetOutputImage().c_str(); + const unsigned int nbComponents = parseResult->IsOptionPresent("--NumComponents") ? + parseResult->GetParameterUInt("--NumComponents") : 0; + unsigned int radiusX = 1; + unsigned int radiusY = 1; + if ( parseResult->IsOptionPresent("--Radius") ) + { + radiusX = parseResult->GetParameterUInt("--Radius",0); + radiusY = parseResult->GetParameterUInt("--Radius",1); + } + const bool normalization = parseResult->IsOptionPresent("--Normalize"); + + // Main type definition + const unsigned int Dimension = 2; + typedef double PixelType; + typedef otb::VectorImage< PixelType, Dimension > ImageType; + + // Reading input images + typedef otb::ImageFileReader<ImageType> ReaderType; + ReaderType::Pointer reader = ReaderType::New(); + reader->SetFileName(inputImageName); + + // Noise filtering + typedef otb::LocalActivityVectorImageFilter< ImageType, ImageType > NoiseFilterType; + NoiseFilterType::RadiusType radius = {{ radiusX, radiusY }}; + + // Image filtering + typedef otb::NAPCAImageFilter< ImageType, ImageType, + NoiseFilterType, otb::Transform::FORWARD > FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput( reader->GetOutput() ); + filter->SetNumberOfPrincipalComponentsRequired( nbComponents ); + filter->SetUseNormalization( normalization ); + filter->GetNoiseImageFilter()->SetRadius( radius ); + + typedef otb::CommandProgressUpdate< FilterType > CommandType; + CommandType::Pointer observer = CommandType::New(); + filter->AddObserver( itk::ProgressEvent(), observer ); + + // Writing + typedef otb::ImageFileWriter< ImageType > ImageWriterType; + ImageWriterType::Pointer writer = ImageWriterType::New(); + writer->SetFileName( outputImageName ); + writer->SetInput( filter->GetOutput() ); + writer->Update(); + + + if ( parseResult->IsOptionPresent("--Inverse") ) + { + typedef otb::NAPCAImageFilter< ImageType, ImageType, + NoiseFilterType, otb::Transform::INVERSE > InvFilterType; + InvFilterType::Pointer invFilter = InvFilterType::New(); + invFilter->SetInput( filter->GetOutput() ); + invFilter->SetMeanValues( filter->GetMeanValues() ); + if ( normalization ) + invFilter->SetStdDevValues( filter->GetStdDevValues() ); + invFilter->SetTransformationMatrix( filter->GetTransformationMatrix() ); + + typedef otb::CommandProgressUpdate< InvFilterType > CommandType2; + CommandType2::Pointer invObserver = CommandType2::New(); + invFilter->AddObserver( itk::ProgressEvent(), invObserver ); + + ImageWriterType::Pointer invWriter = ImageWriterType::New(); + invWriter->SetFileName( parseResult->GetParameterString("--Inverse") ); + invWriter->SetInput( invFilter->GetOutput() ); + invWriter->Update(); + } + + return EXIT_SUCCESS; +} + + + + + -- GitLab