From df49b8d79eb0db3d59299d448964684721e82a22 Mon Sep 17 00:00:00 2001
From: Cyrille Valladeau <cyrille.valladeau@c-s.fr>
Date: Fri, 19 Nov 2010 14:27:23 +0100
Subject: [PATCH] ENH : add test on vnl_minimize
---
.../otbSarParametricMapFunction.txx | 47 +++---
...SarRadiometricCalibrationToImageFilter.txx | 14 +-
Testing/Code/Radiometry/CMakeLists.txt | 5 +-
.../Code/Radiometry/otbRadiometryTests9.cxx | 1 +
.../Code/Radiometry/otbTestVNLminimize.cxx | 140 ++++++++++++++++++
5 files changed, 184 insertions(+), 23 deletions(-)
create mode 100644 Testing/Code/Radiometry/otbTestVNLminimize.cxx
diff --git a/Code/Radiometry/otbSarParametricMapFunction.txx b/Code/Radiometry/otbSarParametricMapFunction.txx
index 98d391ea0b..4e0e44df91 100644
--- a/Code/Radiometry/otbSarParametricMapFunction.txx
+++ b/Code/Radiometry/otbSarParametricMapFunction.txx
@@ -122,16 +122,19 @@ SarParametricMapFunction<TInputImage, TCoordRep>
}
else if (pointSet->GetNumberOfPoints() == 1)
{
+ std::cout<<"EvaluateParam lol1"<<std::endl;
pointSet->GetPointData(0, &pointValue);
m_Coeff(0, 0) = pointValue;
}
else
{
+ std::cout<<"EvaluateParam lol2"<<std::endl;
// Get input region for normalization of coordinates
const itk::MetaDataDictionary& dict = this->GetInputImage()->GetMetaDataDictionary();
ImageKeywordlist imageKeywordlist;
if (dict.HasKey(MetaDataKey::OSSIMKeywordlistKey))
{
+ std::cout<<"EvaluateParam lol2.1"<<std::endl;
itk::ExposeMetaData<ImageKeywordlist>(dict, MetaDataKey::OSSIMKeywordlistKey, imageKeywordlist);
ossimKeywordlist kwl;
imageKeywordlist.convertToOSSIMKeywordlist(kwl);
@@ -142,53 +145,61 @@ SarParametricMapFunction<TInputImage, TCoordRep>
}
else
{
+ std::cout<<"EvaluateParam lol2.2"<<std::endl;
m_ProductHeight = this->GetInputImage()->GetLargestPossibleRegion().GetSize()[0] ;
m_ProductWidth = this->GetInputImage()->GetLargestPossibleRegion().GetSize()[1];
}
-
+std::cout<<"EvaluateParam loli"<<std::endl;
// Perform the plane least square estimation
unsigned int nbRecords = pointSet->GetNumberOfPoints();
- unsigned int nbCoef = m_Coeff.Rows() * m_Coeff.Cols();
+ const unsigned int coeffRows = m_Coeff.Rows();
+ const unsigned int coeffCols = m_Coeff.Cols();
+ unsigned int nbCoef = coeffRows * coeffCols;
+ const double invProductHeight = 1. / m_ProductHeight;
+ const double invProductWidth = 1. / m_ProductWidth;
vnl_sparse_matrix<double> a(nbRecords, nbCoef);
vnl_vector<double> b(nbRecords), bestParams(nbCoef);
b.fill(0);
bestParams.fill(0);
-
+std::cout<<"EvaluateParam lolo"<<std::endl;
// Fill the linear system
for (unsigned int i = 0; i < nbRecords; ++i)
{
this->GetPointSet()->GetPoint(i, &point);
this->GetPointSet()->GetPointData(i, &pointValue);
- b(i) = pointValue;
+ b(i) = pointValue;
//std::cout << "point = " << point << std::endl;
- //std::cout << "b(" << i << ") = " << pointValue << std::endl;
+ std::cout << "b(" << i << ") = " << pointValue << " / "<<point<<std::endl;
- for (unsigned int xcoeff = 0; xcoeff < m_Coeff.Cols(); ++xcoeff)
+ for (unsigned int xcoeff = 0; xcoeff < coeffCols; ++xcoeff)
{
- double xpart = vcl_pow( static_cast<double>(point[0]) / m_ProductWidth, static_cast<double>(xcoeff));
- for (unsigned int ycoeff = 0; ycoeff < m_Coeff.Rows(); ++ycoeff)
+ double xpart = vcl_pow( static_cast<double>(point[0]) + invProductWidth, static_cast<double>(xcoeff));
+ for (unsigned int ycoeff = 0; ycoeff < coeffRows; ++ycoeff)
{
- double ypart = vcl_pow( static_cast<double>(point[1]) / m_ProductHeight, static_cast<double>(ycoeff));
- a(i, xcoeff * m_Coeff.Rows() + ycoeff) = xpart * ypart;
- //std::cout << "a(" << i << "," << xcoeff * m_Coeff.Rows() + ycoeff << ") = " << xpart * ypart << std::endl;
+ double ypart = vcl_pow( static_cast<double>(point[1]) * invProductHeight, static_cast<double>(ycoeff));
+ a(i, xcoeff * coeffRows + ycoeff) = xpart * ypart;
+ std::cout << "a(" << i << "," << xcoeff * coeffRows + ycoeff << ") = " << xpart * ypart << std::endl;
}
}
}
-
+ std::cout<<"EvaluateParam lol"<<std::endl;
// Create the linear system
+ std::cout<<a.cols()<<","<<a.rows()<<" === "<<b.size()<<" === "<<bestParams.size()<<std::endl;
vnl_sparse_matrix_linear_system<double> linearSystem(a, b);
-
+std::cout << "EvaluateParam lol =" << std::endl;
// And solve it
vnl_lsqr linearSystemSolver(linearSystem);
+ std::cout << "EvaluateParam bestParams: " <<bestParams<< std::endl;
linearSystemSolver.minimize(bestParams);
-
- for (unsigned int xcoeff = 0; xcoeff < m_Coeff.Cols(); ++xcoeff)
+
+ std::cout << "m_Coeff bis: " << coeffCols<< " , " << m_Coeff.Rows() << std::endl;
+ for (unsigned int xcoeff = 0; xcoeff < coeffCols; ++xcoeff)
{
- for (unsigned int ycoeff = 0; ycoeff < m_Coeff.Rows(); ++ycoeff)
+ for (unsigned int ycoeff = 0; ycoeff < coeffRows; ++ycoeff)
{
- m_Coeff(ycoeff, xcoeff) = bestParams(xcoeff * m_Coeff.Rows() + ycoeff);
- //std::cout << "m_Coeff(" << ycoeff << "," << xcoeff << ") = " << m_Coeff(ycoeff, xcoeff) << std::endl;
+ m_Coeff(ycoeff, xcoeff) = bestParams(xcoeff * coeffRows + ycoeff);
+ std::cout << "m_Coeff(" << ycoeff << "," << xcoeff << ") = " << m_Coeff(ycoeff, xcoeff) << std::endl;
}
}
}
diff --git a/Code/Radiometry/otbSarRadiometricCalibrationToImageFilter.txx b/Code/Radiometry/otbSarRadiometricCalibrationToImageFilter.txx
index d99e14c8f4..03be545792 100644
--- a/Code/Radiometry/otbSarRadiometricCalibrationToImageFilter.txx
+++ b/Code/Radiometry/otbSarRadiometricCalibrationToImageFilter.txx
@@ -43,6 +43,7 @@ void
SarRadiometricCalibrationToImageFilter<TInputImage, TOutputImage>
::BeforeThreadedGenerateData()
{
+ std::cout<<"Before Starts"<<std::endl;
// will SetInputImage on the function
Superclass::BeforeThreadedGenerateData();
@@ -52,7 +53,7 @@ SarRadiometricCalibrationToImageFilter<TInputImage, TOutputImage>
FunctionPointer function = this->GetFunction();
function->SetScale(imageMetadataInterface->GetRadiometricCalibrationScale());
-
+ std::cout<<"Before Starts 1"<<std::endl;
ParametricFunctionPointer noise;
ParametricFunctionPointer antennaPatternNewGain;
ParametricFunctionPointer antennaPatternOldGain;
@@ -63,7 +64,7 @@ SarRadiometricCalibrationToImageFilter<TInputImage, TOutputImage>
noise->SetPointSet(imageMetadataInterface->GetRadiometricCalibrationNoise());
noise->SetPolynomalSize(imageMetadataInterface->GetRadiometricCalibrationNoisePolynomialDegree());
noise->EvaluateParametricCoefficient();
-
+ std::cout<<"Before Starts 2"<<std::endl;
antennaPatternNewGain = function->GetAntennaPatternNewGain();
antennaPatternNewGain->SetPointSet(imageMetadataInterface->GetRadiometricCalibrationAntennaPatternNewGain());
antennaPatternNewGain->SetPolynomalSize(imageMetadataInterface->GetRadiometricCalibrationAntennaPatternNewGainPolynomialDegree());
@@ -73,16 +74,21 @@ SarRadiometricCalibrationToImageFilter<TInputImage, TOutputImage>
antennaPatternOldGain->SetPointSet(imageMetadataInterface->GetRadiometricCalibrationAntennaPatternOldGain());
antennaPatternOldGain->SetPolynomalSize(imageMetadataInterface->GetRadiometricCalibrationAntennaPatternOldGainPolynomialDegree());
antennaPatternOldGain->EvaluateParametricCoefficient();
-
+ std::cout<<"Before Starts 3"<<std::endl;
incidenceAngle = function->GetIncidenceAngle();
+ std::cout<<function->GetIncidenceAngle()<<std::endl;
incidenceAngle->SetPointSet(imageMetadataInterface->GetRadiometricCalibrationIncidenceAngle());
incidenceAngle->SetPolynomalSize(imageMetadataInterface->GetRadiometricCalibrationIncidenceAnglePolynomialDegree());
+ std::cout<<"================================ Before Starts 3.5"<<std::endl;
+
incidenceAngle->EvaluateParametricCoefficient();
-
+ std::cout<<"Before Starts 4"<<std::endl;
rangeSpreadLoss = function->GetRangeSpreadLoss();
rangeSpreadLoss->SetPointSet(imageMetadataInterface->GetRadiometricCalibrationRangeSpreadLoss());
rangeSpreadLoss->SetPolynomalSize(imageMetadataInterface->GetRadiometricCalibrationRangeSpreadLossPolynomialDegree());
rangeSpreadLoss->EvaluateParametricCoefficient();
+
+ std::cout<<"Before Done"<<std::endl;
}
diff --git a/Testing/Code/Radiometry/CMakeLists.txt b/Testing/Code/Radiometry/CMakeLists.txt
index a8ae745df0..90b3945719 100644
--- a/Testing/Code/Radiometry/CMakeLists.txt
+++ b/Testing/Code/Radiometry/CMakeLists.txt
@@ -1534,7 +1534,9 @@ ADD_TEST(raTvSarBrightnessToImageFilter ${RADIOMETRY_TESTS9}
)
ENDIF(OTB_DATA_USE_LARGEINPUT)
-
+ADD_TEST(raTvTestVNLMinimize ${RADIOMETRY_TESTS9}
+ otbTestVNLMinimize
+)
@@ -1651,6 +1653,7 @@ otbSarRadiometricCalibrationToImageFilterWithRealPixelTest.cxx
otbSarBrightnessFunctor.cxx
otbSarBrightnessFunction.cxx
otbSarBrightnessToImageFilterTest.cxx
+otbTestVNLMinimize.cxx
)
diff --git a/Testing/Code/Radiometry/otbRadiometryTests9.cxx b/Testing/Code/Radiometry/otbRadiometryTests9.cxx
index f96ed4dd19..2feacadd2a 100644
--- a/Testing/Code/Radiometry/otbRadiometryTests9.cxx
+++ b/Testing/Code/Radiometry/otbRadiometryTests9.cxx
@@ -43,4 +43,5 @@ void RegisterTests()
REGISTER_TEST(otbSarBrightnessFunctor);
REGISTER_TEST(otbSarBrightnessFunction);
REGISTER_TEST(otbSarBrightnessToImageFilterTest);
+ REGISTER_TEST(otbTestVNLMinimize);
}
diff --git a/Testing/Code/Radiometry/otbTestVNLminimize.cxx b/Testing/Code/Radiometry/otbTestVNLminimize.cxx
new file mode 100644
index 0000000000..a05547bee0
--- /dev/null
+++ b/Testing/Code/Radiometry/otbTestVNLminimize.cxx
@@ -0,0 +1,140 @@
+/*=========================================================================
+
+ 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 <vnl/algo/vnl_lsqr.h>
+#include <vnl/vnl_sparse_matrix_linear_system.h>
+#include <vnl/vnl_least_squares_function.h>
+#include <iostream>
+
+int otbTestVNLMinimize(int argc, char * argv[])
+{
+ std::cout << "-------------------------- Start with LINUX coeff: "<< std::endl;
+
+ vnl_sparse_matrix<double> a(5, 6);
+ vnl_vector<double> b(5), bestParams(6);
+ b.fill(0);
+ bestParams.fill(0);
+
+
+a(0,0) = 1;
+a(0,1) = 0.500098;
+a(0,2) = 0.5;
+a(0,3) = 0.250049;
+a(0,4) = 0.25;
+a(0,5) = 0.125024;
+
+a(1,0) = 1;
+a(1,1) = 1;
+a(1,2) = 7.93399e-05;
+a(1,3) = 7.93399e-05;
+a(1,4) = 6.29482e-09;
+a(1,5) = 6.29482e-09;
+
+a(2,0) = 1;
+a(2,1) = 1;
+a(2,2) = 0.987544;
+a(2,3) = 0.987544;
+a(2,4) = 0.975242;
+a(2,5) = 0.975242;
+
+a(3,0) = 1;
+a(3,1) = 9.75229e-05;
+a(3,2) = 7.93399e-05;
+a(3,3) = 7.73746e-09;
+a(3,4) = 6.29482e-09;
+a(3,5) = 6.13889e-13;
+
+a(4,0) = 1;
+a(4,1) = 9.75229e-05;
+a(4,2) = 0.987226;
+a(4,3) = 9.62772e-05;
+a(4,4) = 0.974616;
+a(4,5) = 9.50474e-05;
+
+b(0) = 0.451568;
+b(1) = 0.440229;
+b(2) = 0.462236;
+b(3) = 0.44059;
+b(4) = 0.462275;
+
+
+vnl_sparse_matrix_linear_system<double> linearSystem(a, b);
+std::cout << "Satart solve init:" << std::endl;
+ // And solve it
+ vnl_lsqr linearSystemSolver(linearSystem);
+ std::cout << "Solve init done: " <<bestParams<< std::endl;
+ linearSystemSolver.minimize(bestParams);
+ std::cout << "Solve done: " <<bestParams<< std::endl;
+
+std::cout << "-------------------------- Start with visual coeff: "<< std::endl;
+
+ vnl_sparse_matrix<double> a2(5, 6);
+ vnl_vector<double> b2(5), bestParams2(6);
+ b2.fill(0);
+ bestParams2.fill(0);
+
+a2(0,0) = 1;
+a2(0,1) = 0.500098;
+a2(0,2) = 6302;
+a2(0,3) = 3151.61;
+a2(0,4) = 3.97152e+007;
+a2(0,5) = 1.98615e+007;
+
+a2(1,0) = 1;
+a2(1,1) = 1;
+a2(1,2) = 1.00008;
+a2(1,3) = 1.00008;
+a2(1,4) = 1.00016;
+a2(1,5) = 1.00016;
+
+a2(2,0) = 1;
+a2(2,1) = 1;
+a2(2,2) = 12447;
+a2(2,3) = 12447;
+a2(2,4) = 1.54928e+008;
+a2(2,5) = 1.54928e+008;
+
+a2(3,0) = 1;
+a2(3,1) = 9.75229e-005;
+a2(3,2) = 1.00008;
+a2(3,3) = 9.75307e-005;
+a2(3,4) = 1.00016;
+a2(3,5) = 9.75384e-005;
+a2(4,0) = 1;
+a2(4,1) = 9.75229e-005;
+a2(4,2) = 12443;
+a2(4,3) = 1.21348;
+a2(4,4) = 1.54828e+008;
+a2(4,5) = 15099.3;
+
+b2(0) = 0.451568;
+b2(1) = 0.440229;
+b2(2) = 0.462236;
+b2(3) = 0.44059;
+b2(4) = 0.462275;
+
+vnl_sparse_matrix_linear_system<double> linearSystem2(a2, b2);
+std::cout << "Satart solve init:" << std::endl;
+ // And solve it
+ vnl_lsqr linearSystemSolver2(linearSystem2);
+ std::cout << "Solve init done: " <<bestParams2<< std::endl;
+ linearSystemSolver2.minimize(bestParams2);
+ std::cout << "Solve done: " <<bestParams2<< std::endl;
+
+ return EXIT_SUCCESS;
+}
\ No newline at end of file
--
GitLab