From a2456748f2f99ad312a7c8b60d1f46914f49177c Mon Sep 17 00:00:00 2001
From: Julien Michel <julien.michel@orfeo-toolbox.org>
Date: Mon, 3 Dec 2012 15:50:32 +0100
Subject: [PATCH] TEST: Adding tests for the new RPCSolverAdapter class
---
Testing/Code/UtilitiesAdapters/CMakeLists.txt | 37 ++++
.../otbRPCSolverAdapterTest.cxx | 180 ++++++++++++++++++
.../otbUtilitiesAdaptersTests1.cxx | 1 +
3 files changed, 218 insertions(+)
create mode 100644 Testing/Code/UtilitiesAdapters/otbRPCSolverAdapterTest.cxx
diff --git a/Testing/Code/UtilitiesAdapters/CMakeLists.txt b/Testing/Code/UtilitiesAdapters/CMakeLists.txt
index cd155d1d3b..e7e06978d7 100644
--- a/Testing/Code/UtilitiesAdapters/CMakeLists.txt
+++ b/Testing/Code/UtilitiesAdapters/CMakeLists.txt
@@ -32,12 +32,49 @@ ADD_TEST(uaTvGeometricSarSensorModelAdapterTest
otbGeometricSarSensorModelAdapterTest
${LARGEINPUT}/TERRASARX/2008-03-10_GrandCanyon_SSC/TSX1_SAR__SSC______SM_S_SRA_20080310T133220_20080310T133228/TSX1_SAR__SSC______SM_S_SRA_20080310T133220_20080310T133228.xml
${TEMP}/uaTvGeometricSarSensorModelAdapterTest.txt)
+
+ADD_TEST(uaTvRPCSolverAdapterValidationTest
+ ${UtilitiesAdapters_TESTS1}
+ otbRPCSolverAdapterTest
+ ${LARGEINPUT}/QUICKBIRD/TOULOUSE/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF
+ 10 0.25 0.35
+ ${INPUTDATA}/DEM/srtm_directory/
+ ${INPUTDATA}/DEM/egm96.grd
+)
+
+ADD_TEST(uaTvRPCSolverAdapterOutGeomTest
+ ${UtilitiesAdapters_TESTS1}
+ --compare-ascii ${EPSILON_9}
+ ${BASELINE_FILES}/uaTvRPCSolverAdapterOutGeomTest.geom
+ ${TEMP}/uaTvRPCSolverAdapterOutGeomTest.geom
+ otbRPCSolverAdapterTest
+ ${LARGEINPUT}/QUICKBIRD/TOULOUSE/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF
+ 10 0.25 0.35
+ ${INPUTDATA}/DEM/srtm_directory/
+ ${INPUTDATA}/DEM/egm96.grd
+ ${TEMP}/uaTvRPCSolverAdapterOutGeomTest.geom
+)
+
+ADD_TEST(uaTvRPCSolverAdapterNotEnoughPointsTest
+ ${UtilitiesAdapters_TESTS1}
+ otbRPCSolverAdapterTest
+ ${LARGEINPUT}/QUICKBIRD/TOULOUSE/000000128955_01_P001_PAN/02APR01105228-P1BS-000000128955_01_P001.TIF
+ 4 0.25 0.35
+ ${INPUTDATA}/DEM/srtm_directory/
+ ${INPUTDATA}/DEM/egm96.grd
+ )
+SET_TESTS_PROPERTIES(uaTvRPCSolverAdapterNotEnoughPointsTest PROPERTIES WILL_FAIL TRUE)
+
+
+
+
ENDIF(OTB_DATA_USE_LARGEINPUT)
SET(UtilitiesAdapters_SRCS1
otbUtilitiesAdaptersTests1.cxx
otbPlatformPositionAdapter.cxx
otbGeometricSarSensorModelAdapter.cxx
+otbRPCSolverAdapterTest.cxx
)
diff --git a/Testing/Code/UtilitiesAdapters/otbRPCSolverAdapterTest.cxx b/Testing/Code/UtilitiesAdapters/otbRPCSolverAdapterTest.cxx
new file mode 100644
index 0000000000..61e948f3e9
--- /dev/null
+++ b/Testing/Code/UtilitiesAdapters/otbRPCSolverAdapterTest.cxx
@@ -0,0 +1,180 @@
+/*=========================================================================
+
+ 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 "otbMacro.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include "otbGenericRSTransform.h"
+#include "otbGeographicalDistance.h"
+#include "itkEuclideanDistance.h"
+#include "otbDEMHandler.h"
+#include "otbSensorModelAdapter.h"
+#include "otbRPCSolverAdapter.h"
+
+typedef otb::Image<double> ImageType;
+typedef otb::ImageFileReader<ImageType> ReaderType;
+typedef otb::GenericRSTransform<> RSTranformType;
+typedef otb::RPCSolverAdapter::Point2DType Point2DType;
+typedef otb::RPCSolverAdapter::Point3DType Point3DType;
+typedef otb::GenericRSTransform<double,3,3> RSTranform3dType;
+typedef itk::Statistics::EuclideanDistance<Point2DType> EuclideanDistanceType;
+typedef otb::GeographicalDistance<Point2DType> GeoDistanceType;
+
+int otbRPCSolverAdapterTest(int argc, char* argv[])
+{
+ // This test takes a sensor model (possibly a rpc one), use it to
+ // generate gcps and estimate a rpc model. It then checks the
+ // precision of both forward and inverse transform
+ std::string infname = argv[1];
+ const unsigned int gridSize = atoi(argv[2]);
+ const double geoTol = atof(argv[3]);
+ const double imgTol = atof(argv[4]);
+ const std::string demdir = argv[5];
+ const std::string geoid = argv[6];
+
+ std::cout<<"GeoTol: "<<geoTol<<" meters"<<std::endl;
+ std::cout<<"ImgTol: "<<imgTol<<" pixels"<<std::endl;
+
+ otb::DEMHandler::Instance()->OpenDEMDirectory(demdir);
+ otb::DEMHandler::Instance()->OpenGeoidFile(geoid);
+
+ ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName(infname);
+ reader->UpdateOutputInformation();
+
+ RSTranformType::Pointer fwd2dTransform = RSTranformType::New();
+ fwd2dTransform->SetInputKeywordList(reader->GetOutput()->GetImageKeywordlist());
+ fwd2dTransform->InstanciateTransform();
+
+ ImageType::SizeType size = reader->GetOutput()->GetLargestPossibleRegion().GetSize();
+
+ unsigned int stepx = size[0]/gridSize;
+ unsigned int stepy = size[1]/gridSize;
+
+ otb::RPCSolverAdapter::GCPsContainerType gcps;
+
+ // Generate gcps
+ for(unsigned int i = 0; i < gridSize;++i)
+ {
+ for(unsigned int j = 0; j < gridSize;++j)
+ {
+ ImageType::IndexType currentIndex;
+ currentIndex[0] = i*stepx;
+ currentIndex[1] = j*stepy;
+
+ Point2DType currentPoint, currentWgs84Point;
+ reader->GetOutput()->TransformIndexToPhysicalPoint(currentIndex,currentPoint);
+
+ currentWgs84Point = fwd2dTransform->TransformPoint(currentPoint);
+
+ double height = otb::DEMHandler::Instance()->GetHeightAboveEllipsoid(currentWgs84Point);
+
+ Point3DType current3DWgs84Point;
+ current3DWgs84Point[0] = currentWgs84Point[0];
+ current3DWgs84Point[1] = currentWgs84Point[1];
+ current3DWgs84Point[2] = height;
+
+ gcps.push_back(std::make_pair(currentPoint,current3DWgs84Point));
+
+ std::cout<<"Adding gcp: "<<currentPoint<<" <-> " << current3DWgs84Point<<std::endl;
+ }
+ }
+
+ // Solve rpc
+ otb::ImageKeywordlist rpcKwl;
+ double rmse;
+
+ Point2DType ul = reader->GetOutput()->GetOrigin();
+ Point2DType lr;
+ lr[0] = ul[0] + size[0]*reader->GetOutput()->GetSpacing()[0];
+ lr[1] = ul[1] + size[1]*reader->GetOutput()->GetSpacing()[1];
+
+ // Call solver: either write geom and exit, or evaluate model precision
+ if(argc==8)
+ {
+ const std::string outgeom = argv[7];
+ bool success = otb::RPCSolverAdapter::Solve(gcps,ul,lr,rmse,outgeom);
+ if(success)
+ return EXIT_SUCCESS;
+ else
+ return EXIT_FAILURE;
+ }
+
+ otb::RPCSolverAdapter::Solve(gcps,ul,lr,rmse,rpcKwl);
+
+ std::cout<<"Optimization done, RMSE="<<rmse<<std::endl;
+
+ // Build forward and inverse rpc transform
+ RSTranform3dType::Pointer rpcFwdTransform = RSTranform3dType::New();
+ rpcFwdTransform->SetInputKeywordList(rpcKwl);
+ rpcFwdTransform->InstanciateTransform();
+ RSTranformType::Pointer rpcInvTransform = RSTranformType::New();
+ rpcInvTransform->SetOutputKeywordList(rpcKwl);
+ rpcInvTransform->InstanciateTransform();
+
+ EuclideanDistanceType::Pointer euclideanDistance = EuclideanDistanceType::New();
+ GeoDistanceType::Pointer geoDistance = GeoDistanceType::New();
+
+ bool fail = false;
+
+ for(otb::RPCSolverAdapter::GCPsContainerType::iterator it = gcps.begin();
+ it != gcps.end();++it)
+ {
+ Point2DType imgPoint,groundPoint, groundPoint2dRef;
+ Point3DType imgPoint3D,groundPoint3D;
+
+ groundPoint2dRef[0] = it->second[0];
+ groundPoint2dRef[1] = it->second[1];
+
+ // Check forward transform
+ imgPoint3D[0] = it->first[0];
+ imgPoint3D[1] = it->first[1];
+ imgPoint3D[2] = it->second[2];
+
+ groundPoint3D = rpcFwdTransform->TransformPoint(imgPoint3D);
+
+ groundPoint[0] = groundPoint3D[0];
+ groundPoint[1] = groundPoint3D[1];
+
+ double groundRes = geoDistance->Evaluate(groundPoint,groundPoint2dRef);
+
+ if(groundRes>geoTol)
+ {
+ fail = true;
+ std::cerr<<"Imprecise result with forward estimated model: fwd("<<it->first<<") = "<<groundPoint<<", but reference is "<<groundPoint2dRef<<" error: "<<groundRes<<" meters)"<<std::endl;
+ }
+
+ // Check inverse transform
+ imgPoint = rpcInvTransform->TransformPoint(groundPoint2dRef);
+
+ double imgRes = euclideanDistance->Evaluate(imgPoint,it->first);
+
+ if(imgRes>imgTol)
+ {
+ fail = true;
+ std::cerr<<"Imprecise result with inverse estimated model: inv("<<groundPoint2dRef<<") = "<<imgPoint<<", but reference is "<<it->first<<" error: "<<imgRes<<" pixels)"<<std::endl;
+ }
+ }
+
+
+ if(fail)
+ {
+ return EXIT_FAILURE;
+ }
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/UtilitiesAdapters/otbUtilitiesAdaptersTests1.cxx b/Testing/Code/UtilitiesAdapters/otbUtilitiesAdaptersTests1.cxx
index 27df316b36..8e1eb88c51 100644
--- a/Testing/Code/UtilitiesAdapters/otbUtilitiesAdaptersTests1.cxx
+++ b/Testing/Code/UtilitiesAdapters/otbUtilitiesAdaptersTests1.cxx
@@ -27,4 +27,5 @@ void RegisterTests()
REGISTER_TEST(otbPlatformPositionComputeBaselineTest);
REGISTER_TEST(otbGeometricSarSensorModelAdapterNewTest);
REGISTER_TEST(otbGeometricSarSensorModelAdapterTest);
+ REGISTER_TEST(otbRPCSolverAdapterTest);
}
--
GitLab