Commit 7d24cd40 authored by Cédric Traizet's avatar Cédric Traizet
Browse files

REFAC: use otb::RPCSolver instead of Ossim in GCPsToRPCSensorModelImageFilter

parent 2d6aa87c
Pipeline #6611 failed with stages
in 93 minutes and 3 seconds
Ossim Keyword list:
adjustment_0.adj_param_0.center: 0
adjustment_0.adj_param_0.description: intrack_offset
adjustment_0.adj_param_0.lock_flag: 0
adjustment_0.adj_param_0.parameter: 0
adjustment_0.adj_param_0.sigma: 50
adjustment_0.adj_param_0.units: pixel
adjustment_0.adj_param_1.center: 0
adjustment_0.adj_param_1.description: crtrack_offset
adjustment_0.adj_param_1.lock_flag: 0
adjustment_0.adj_param_1.parameter: 0
adjustment_0.adj_param_1.sigma: 50
adjustment_0.adj_param_1.units: pixel
adjustment_0.adj_param_2.center: 0
adjustment_0.adj_param_2.description: intrack_scale
adjustment_0.adj_param_2.lock_flag: 0
adjustment_0.adj_param_2.parameter: 0
adjustment_0.adj_param_2.sigma: 50
adjustment_0.adj_param_2.units: unknown
adjustment_0.adj_param_3.center: 0
adjustment_0.adj_param_3.description: crtrack_scale
adjustment_0.adj_param_3.lock_flag: 0
adjustment_0.adj_param_3.parameter: 0
adjustment_0.adj_param_3.sigma: 50
adjustment_0.adj_param_3.units: unknown
adjustment_0.adj_param_4.center: 0
adjustment_0.adj_param_4.description: map_rotation
adjustment_0.adj_param_4.lock_flag: 0
adjustment_0.adj_param_4.parameter: 0
adjustment_0.adj_param_4.sigma: 0.1
adjustment_0.adj_param_4.units: degrees
adjustment_0.description: Initial adjustment
adjustment_0.dirty_flag: 0
adjustment_0.number_of_params: 5
current_adjustment: 0
height_off: 0
height_scale: 4.5035996273705e+15
lat_off: 43.608891
lat_scale: 1
line_den_coeff_0: 1
line_den_coeff_1: 0.116647929191405
line_den_coeff_10: 0
line_den_coeff_11: 4.24205291306071e-06
line_den_coeff_12: -1.22170433593066e-05
line_den_coeff_13: 0
line_den_coeff_14: -1.62292152211888e-05
line_den_coeff_15: -1.11751035274519e-05
line_den_coeff_16: 0
line_den_coeff_17: 0
line_den_coeff_18: 0
line_den_coeff_19: 0
line_den_coeff_2: -0.795796952135346
line_den_coeff_3: 0
line_den_coeff_4: -0.00176250458824308
line_den_coeff_5: 0
line_den_coeff_6: 0
line_den_coeff_7: -0.69205662487881
line_den_coeff_8: -0.725158534735126
line_den_coeff_9: 0
line_num_coeff_0: 0.00176507025129761
line_num_coeff_1: 0.469825903760047
line_num_coeff_10: 0
line_num_coeff_11: 4.1371912315397e-05
line_num_coeff_12: -2.66115593034605e-05
line_num_coeff_13: 0
line_num_coeff_14: -0.00378322572708931
line_num_coeff_15: -0.00397219871280679
line_num_coeff_16: 0
line_num_coeff_17: 0
line_num_coeff_18: 0
line_num_coeff_19: 0
line_num_coeff_2: -182.988749105524
line_num_coeff_3: 1.00711799423683e-14
line_num_coeff_4: 0.000167599353017481
line_num_coeff_5: 4.23534328933673e-23
line_num_coeff_6: 0
line_num_coeff_7: -0.000132255684215101
line_num_coeff_8: -0.00299521983182487
line_num_coeff_9: 0
line_off: 2224.5
line_scale: 225.5
long_off: 1.4252414
long_scale: 1
number_of_adjustments: 1
polynomial_format: B
samp_den_coeff_0: 1
samp_den_coeff_1: -0.0140652497583194
samp_den_coeff_10: 0
samp_den_coeff_11: -2.11567008497461e-07
samp_den_coeff_12: -2.52699758468156e-06
samp_den_coeff_13: 0
samp_den_coeff_14: 1.28729228360665e-06
samp_den_coeff_15: 7.0821267053196e-06
samp_den_coeff_16: 0
samp_den_coeff_17: 0
samp_den_coeff_18: 0
samp_den_coeff_19: 0
samp_den_coeff_2: 0.056111938809245
samp_den_coeff_3: 0
samp_den_coeff_4: -0.00165427540398437
samp_den_coeff_5: 0
samp_den_coeff_6: 0
samp_den_coeff_7: -0.720295418717795
samp_den_coeff_8: -0.238046475331491
samp_den_coeff_9: 0
samp_num_coeff_0: -4.03022891014623e-05
samp_num_coeff_1: 140.122696755472
samp_num_coeff_10: 0
samp_num_coeff_11: 0.00514045286544423
samp_num_coeff_12: 0.00169871244049508
samp_num_coeff_13: 0
samp_num_coeff_14: 3.04771969411287e-06
samp_num_coeff_15: 3.0597489842213e-05
samp_num_coeff_16: 0
samp_num_coeff_17: 0
samp_num_coeff_18: 0
samp_num_coeff_19: 0
samp_num_coeff_2: 0.720548641464637
samp_num_coeff_3: 1.72815049746881e-14
samp_num_coeff_4: -0.000307017880169014
samp_num_coeff_5: -1.67843825237436e-25
samp_num_coeff_6: 0
samp_num_coeff_7: 7.0629169037717e-05
samp_num_coeff_8: 0.000393857191344697
samp_num_coeff_9: 0
samp_off: 2224.5
samp_scale: 225.5
type: ossimRpcProjection
Residual ground error: 1.2432282311
{"LineOffset": "2224.5", "SampleOffset": "2224.5", "LatOffset": "43.6089", "LonOffset": "1.42524", "HeightOffset": "0", "LineScale": "225.5", "SampleScale": "225.5", "LatScale": "1", "LonScale": "1", "HeightScale": "1", "LineNum": [ "0.00174891", "0.469879", "-182.989", "-4.38776e-15", "0.000183514", "3.26103e-23", "1.88011e-36", "-0.000133312", "-0.00295522", "0", "0", "4.13724e-05", "-2.66117e-05", "0", "-0.0037833", "-0.00397229", "0", "0", "0", "0", ], "LineDen": [ "1", "0.119555", "-0.788492", "0", "-0.00176245", "0", "0", "-0.69207", "-0.725175", "0", "0", "4.36024e-06", "-1.21576e-05", "0", "-1.60709e-05", "-1.09766e-05", "0", "0", "0", "0", ], "SampleNum": [ "-4.02927e-05", "140.123", "0.720552", "-6.37125e-15", "-0.000306972", "3.63979e-25", "7.44761e-39", "7.06041e-05", "0.000393491", "0", "0", "0.0051405", "0.00169872", "0", "3.04758e-06", "3.05976e-05", "0", "0", "0", "0", ], "SampleDen": [ "1", "-0.0140619", "0.0561058", "0", "-0.00165426", "0", "0", "-0.720302", "-0.238048", "0", "0", "-2.11371e-07", "-2.52599e-06", "0", "1.28711e-06", "7.08205e-06", "0", "0", "0", "0", ], }
Residual ground error: 0.1243213845
......@@ -140,9 +140,6 @@ public:
/** Set the GCP container */
void SetGCPsContainer(const GCPsContainerType& container);
/** Get Keywordlist */
itkGetConstReferenceMacro(Keywordlist, ImageKeywordlist);
/** Add a GCP to the GCPContainer. This version of the AddGCP method
* accepts a 3D ground point and does not use DEM or MeanElevation
* to handle the elevation */
......@@ -205,8 +202,8 @@ private:
/** Container of GCPs */
GCPsContainerType m_GCPsContainer;
/** Keywordlist */
ImageKeywordlist m_Keywordlist;
/** ImageMetadata */
ImageMetadata m_ImageMetadata;
/** Flag to see if model is up-to-date */
mutable bool m_ModelUpToDate;
......
......@@ -219,17 +219,12 @@ void GCPsToRPCSensorModelImageFilter<TImage>::LoadImageGCPs()
template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::ComputeErrors()
{
typedef GenericRSTransform<double, 3, 3> RSTransformType;
using RSTransformType = GenericRSTransform<double, 2, 2>;
RSTransformType::Pointer rsTransform = RSTransformType::New();
  • It be nice to get rid of all these unneeded dynamic allocations around transformations.

  • but the constructor is protected for ITK Object, New() is needed to create the object

  • I know, I known :(

    I've started an experiment in the past to get rid of this, and simplify as much as possible (for instance composite transforms are no longer have final specializations...) but this requires a lot of changes.

    I'll have to see if I can extract my experiment and share it -- along the way I was trying to simplify SARCalibration

Please register or sign in to reply
rsTransform->SetInputImageMetadata(&(this->GetOutput()->GetImageMetadata()));
rsTransform->InstantiateTransform();
ContinuousIndexType idFix, idOut;
Continuous3DIndexType idOut3D, idTrans3D;
rsTransform->SetInputImageMetadata(&m_ImageMetadata);
Point2DType sensorPoint;
Point3DType groundPoint;
rsTransform->InstantiateTransform();
double sum = 0.;
m_MeanError = 0.;
......@@ -240,30 +235,22 @@ void GCPsToRPCSensorModelImageFilter<TImage>::ComputeErrors()
for (unsigned int i = 0; i < m_GCPsContainer.size(); ++i)
{
// GCP value
sensorPoint = m_GCPsContainer[i].first;
groundPoint = m_GCPsContainer[i].second;
const auto & sensorPoint = m_GCPsContainer[i].first;
const auto & groundPoint = m_GCPsContainer[i].second;
// Compute Transform
Point3DType groundPointTemp, sensorPointTemp;
sensorPointTemp[0] = sensorPoint[0];
sensorPointTemp[1] = sensorPoint[1];
sensorPointTemp[2] = groundPoint[2];
groundPointTemp = rsTransform->TransformPoint(sensorPointTemp);
// Compute Euclidian distance
idFix[0] = sensorPoint[0];
idFix[1] = sensorPoint[1];
auto outPoint = rsTransform->TransformPoint(sensorPoint);
idOut3D[0] = groundPoint[0];
idOut3D[1] = groundPoint[1];
idOut3D[2] = groundPoint[2];
Point2DType groundPoint2D;
idTrans3D[0] = groundPointTemp[0];
idTrans3D[1] = groundPointTemp[1];
idTrans3D[2] = groundPointTemp[2];
groundPoint2D[0] = groundPoint[0];
groundPoint2D[1] = groundPoint[1];
double error = idOut3D.EuclideanDistanceTo(idTrans3D);
double error = outPoint.EuclideanDistanceTo(outPoint);
// Add error to the container
m_ErrorsContainer.push_back(error);
......@@ -282,30 +269,28 @@ void GCPsToRPCSensorModelImageFilter<TImage>::GenerateOutputInformation()
Superclass::GenerateOutputInformation();
// First, retrieve the image pointer
typename TImage::Pointer imagePtr = this->GetOutput();
auto imagePtr = this->GetOutput();
if (!m_ModelUpToDate)
{
double rmsError;
ImageKeywordlist otb_kwl;
otb::RPCSolverAdapter::Solve(m_GCPsContainer, rmsError, otb_kwl);
double rmsError;
Projection::RPCParam rpcParams;
otb::RPCSolver::Solve(m_GCPsContainer, rmsError, rpcParams);
// Retrieve the residual ground error
m_RMSGroundError = rmsError;
m_ImageMetadata = imagePtr->GetImageMetadata();
m_ImageMetadata.Add(MDGeom::RPC, rpcParams);
// Compute errors
this->ComputeErrors();
m_Keywordlist = otb_kwl;
m_ModelUpToDate = true;
}
// Update image keywordlist. We need to do it at each stream
// because it will be overwritten by up-stream
// GenerateOutputInformation() calls.
itk::MetaDataDictionary& dict = imagePtr->GetMetaDataDictionary();
itk::EncapsulateMetaData<ImageKeywordlist>(dict, MetaDataKey::OSSIMKeywordlistKey, m_Keywordlist);
imagePtr->SetImageMetadata(m_ImageMetadata);
}
template <class TImage>
......
......@@ -91,9 +91,8 @@ void PhysicalToRPCSensorModelImageFilter<TImage>::GenerateOutputInformation()
otbGenericMsgDebugMacro(<< "RPC model estimated. RMS ground error: " << m_GCPsToSensorModelFilter->GetRMSGroundError()
<< ", Mean error: " << m_GCPsToSensorModelFilter->GetMeanError());
// Encapsulate the keywordlist
itk::MetaDataDictionary& dict = this->GetOutput()->GetMetaDataDictionary();
itk::EncapsulateMetaData<ImageKeywordlist>(dict, MetaDataKey::OSSIMKeywordlistKey, m_GCPsToSensorModelFilter->GetKeywordlist());
this->GetOutput()->SetImageMetadata(m_GCPsToSensorModelFilter->GetOutput()->GetImageMetadata());
// put the flag to true
m_OutputInformationGenerated = true;
......
......@@ -37,9 +37,14 @@ int otbGCPsToRPCSensorModelImageFilterCheckRpcModel(int argc, char* argv[])
nbPoints--; // last argument in not a gcp pairs point
tol = stoi(s_tol.substr(s_tol.find("=") + 1));
}
// Set the DEM Directory
if (std::string(argv[2]).compare("no_output") != 0)
{
otb::DEMHandler::GetInstance().OpenDEMDirectory(argv[2]);
}
// Check if the number of gcp pairs point is consistent
if (nbPoints % 5 != 0)
{
std::cerr << "Inconsistent GCPs description!" << std::endl;
......@@ -51,7 +56,7 @@ int otbGCPsToRPCSensorModelImageFilterCheckRpcModel(int argc, char* argv[])
typedef otb::GCPsToRPCSensorModelImageFilter<ImageType> GCPsToSensorModelFilterType;
typedef GCPsToSensorModelFilterType::Point2DType Point2DType;
typedef GCPsToSensorModelFilterType::Point3DType Point3DType;
typedef otb::GenericRSTransform<double, 3, 3> GenericRSTransformType;
typedef otb::GenericRSTransform<double, 2, 2> GenericRSTransformType;
typedef otb::GeographicalDistance<ImageType::PointType> GeoDistanceType;
ReaderType::Pointer reader = ReaderType::New();
......@@ -90,15 +95,8 @@ int otbGCPsToRPCSensorModelImageFilterCheckRpcModel(int argc, char* argv[])
// geographical coordinates.
GenericRSTransformType::Pointer grsTrasnform = GenericRSTransformType::New();
//TODO OSSIM:Update GCPsToRPCSensorModelImageFilter so it uses ImageMetdata instead of KeywordList
//grsTrasnform->SetInputKeywordList(rpcEstimator->GetKeywordlist());
otbLogMacro(Debug, << rpcEstimator->GetKeywordlist());
grsTrasnform->SetInputImageMetadata(&rpcEstimator->GetOutput()->GetImageMetadata());
grsTrasnform->SetOutputProjectionRef("EPSG:4326");
// Set the DEM Directory
if (std::string(argv[2]).compare("no_output") != 0)
otb::DEMHandler::GetInstance().OpenDEMDirectory(argv[2]);
grsTrasnform->InstantiateTransform();
// Test
......@@ -107,17 +105,11 @@ int otbGCPsToRPCSensorModelImageFilterCheckRpcModel(int argc, char* argv[])
for (unsigned int gcpId = 0; gcpId < nbGCPs; ++gcpId)
{
Point3DType point;
Point2DType point;
point[0] = std::stof(argv[3 + gcpId * 5]);
point[1] = std::stof(argv[4 + gcpId * 5]);
point[2] = std::stof(argv[7 + gcpId * 5]);
Point3DType transformedPoint;
transformedPoint = grsTrasnform->TransformPoint(point);
Point2DType transformedPoint2D;
transformedPoint2D[0] = transformedPoint[0];
transformedPoint2D[1] = transformedPoint[1];
auto transformedPoint = grsTrasnform->TransformPoint(point);
// reference point
Point2DType geoPoint;
......@@ -125,18 +117,18 @@ int otbGCPsToRPCSensorModelImageFilterCheckRpcModel(int argc, char* argv[])
geoPoint[1] = std::stof(argv[6 + gcpId * 5]);
// Search for nans
if (vnl_math_isnan(transformedPoint2D[0]) || vnl_math_isnan(transformedPoint2D[1]))
if (vnl_math_isnan(transformedPoint[0]) || vnl_math_isnan(transformedPoint[1]))
{
otbLogMacro(Warning, << "Reference : " << geoPoint << " --> Result of the reprojection using the estimated RpcModel " << transformedPoint2D);
otbLogMacro(Warning, << "Reference : " << geoPoint << " --> Result of the reprojection using the estimated RpcModel " << transformedPoint);
otbLogMacro(Warning, << "The result of the projection is nan, there is a problem with the estimated RpcModel");
isErrorDetected = true;
}
// Search for wrong projection results
double residual = geoDistance->Evaluate(geoPoint, transformedPoint2D);
double residual = geoDistance->Evaluate(geoPoint, transformedPoint);
if (residual > tol)
{
otbLogMacro(Warning, << "Reference : " << geoPoint << " --> Result of the reprojection using the estimated RpcModel " << transformedPoint2D << std::endl
otbLogMacro(Warning, << "Reference : " << geoPoint << " --> Result of the reprojection using the estimated RpcModel " << transformedPoint << std::endl
<< " Residual [" << residual << "] is higher than the tolerance [" << tol << "], there is a problem with the estimated RpcModel");
isErrorDetected = true;
}
......
......@@ -80,7 +80,9 @@ int otbGCPsToRPCSensorModelImageFilterWithoutDEM(int argc, char* argv[])
ofs.setf(std::ios::fixed, std::ios::floatfield);
ofs.precision(10);
ofs << rpcEstimator->GetOutput()->GetImageKeywordlist() << std::endl;
auto outputRPC = boost::any_cast<otb::Projection::RPCParam>(rpcEstimator->GetOutput()->GetImageMetadata()[otb::MDGeom::RPC]);
ofs << outputRPC.ToJSON() << std::endl;
ofs << "Residual ground error: " << rpcEstimator->GetRMSGroundError() << std::endl;
ofs.close();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment