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

Merge remote-tracking branch 'origin/2160-SarSensorModel' into 2159

parents 8b32a001 08a2be01
......@@ -115,6 +115,7 @@ public:
std::vector<AzimuthFmRate> GetAzimuthFmRateGeom() const;
std::vector<DopplerCentroid> GetDopplerCentroidGeom() const;
std::vector<Orbit> GetOrbitsGeom() const;
std::vector<BurstRecord> GetBurstRecordsGeom() const;
std::vector<SARNoise> GetNoiseVectorGeom() const;
/**
......
......@@ -236,8 +236,13 @@ std::vector<Orbit> SarImageMetadataInterface::GetOrbitsGeom() const
std::vector<Orbit> orbitVector;
// Number of entries in the vector
int listCount = m_MetadataSupplierInterface->GetAs<int>("orbitList.nb_orbits");
// This streams wild hold the iteration number
// This streams will hold the iteration number
std::ostringstream oss;
std::stringstream ss;
auto facet = new boost::posix_time::time_input_facet("%Y-%m-%dT%H:%M:%S%F");
ss.imbue(std::locale(std::locale(), facet));
for (int listId = 0 ; listId <= listCount - 1 ; ++listId)
{
oss.str("");
......@@ -245,18 +250,56 @@ std::vector<Orbit> SarImageMetadataInterface::GetOrbitsGeom() const
// Base path to the data, that depends on the iteration number
std::string path_root = "orbitList.orbit[" + oss.str() + "]";
Orbit orbit;
std::istringstream(m_MetadataSupplierInterface->GetAs<std::string>(path_root + ".time")) >> orbit.time;
orbit.position = m_MetadataSupplierInterface->GetAs<double>(path_root + ".x_pos");
orbit.position = m_MetadataSupplierInterface->GetAs<double>(path_root + ".y_pos");
orbit.position = m_MetadataSupplierInterface->GetAs<double>(path_root + ".z_pos");
orbit.velocity = m_MetadataSupplierInterface->GetAs<double>(path_root + ".x_vel");
orbit.velocity = m_MetadataSupplierInterface->GetAs<double>(path_root + ".y_vel");
orbit.velocity = m_MetadataSupplierInterface->GetAs<double>(path_root + ".z_vel");
ss << m_MetadataSupplierInterface->GetAs<std::string>(path_root + ".time");
ss >> orbit.time;
orbit.position[0] = m_MetadataSupplierInterface->GetAs<double>(path_root + ".x_pos");
orbit.position[1] = m_MetadataSupplierInterface->GetAs<double>(path_root + ".y_pos");
orbit.position[2] = m_MetadataSupplierInterface->GetAs<double>(path_root + ".z_pos");
orbit.velocity[0] = m_MetadataSupplierInterface->GetAs<double>(path_root + ".x_vel");
orbit.velocity[1] = m_MetadataSupplierInterface->GetAs<double>(path_root + ".y_vel");
orbit.velocity[2] = m_MetadataSupplierInterface->GetAs<double>(path_root + ".z_vel");
orbitVector.push_back(std::move(orbit));
}
return orbitVector;
}
std::vector<BurstRecord> SarImageMetadataInterface::GetBurstRecordsGeom() const
{
const std::string prefix = "support_data.";
std::vector<BurstRecord> burstRecords;
int listCount = m_MetadataSupplierInterface->GetAs<int>(prefix + "geom.bursts.number");
std::stringstream ss;
auto facet = new boost::posix_time::time_input_facet("%Y-%m-%dT%H:%M:%S%F");
ss.imbue(std::locale(std::locale(), facet));
for (int listId = 0 ; listId <= listCount - 1 ; ++listId)
{
const std::string burstName = prefix + "geom.bursts.burst[" + std::to_string(listId) + "].";
BurstRecord record;
ss << m_MetadataSupplierInterface->GetAs<std::string>(burstName + "azimuth_start_time");
ss >> record.azimuthStartTime;
ss << m_MetadataSupplierInterface->GetAs<std::string>(burstName + "azimuth_stop_time");
ss >> record.azimuthStopTime;
record.azimuthAnxTime = m_MetadataSupplierInterface->GetAs<double>(burstName + "azimuth_anx_time");
record.startLine = m_MetadataSupplierInterface->GetAs<int>(burstName + "start_line");
record.endLine = m_MetadataSupplierInterface->GetAs<int>(burstName + "end_line");
record.startSample = m_MetadataSupplierInterface->GetAs<int>(burstName + "start_sample");
record.endSample = m_MetadataSupplierInterface->GetAs<int>(burstName + "end_sample");
burstRecords.push_back(std::move(record));
}
return burstRecords;
}
std::vector<SARNoise> SarImageMetadataInterface::GetNoiseVectorGeom() const
{
std::vector<SARNoise> noiseVector;
......@@ -294,6 +337,21 @@ bool SarImageMetadataInterface::GetSAR(SARParam & sarParam) const
sarParam.dopplerCentroids = this->GetDopplerCentroidGeom();
sarParam.orbits = this->GetOrbitsGeom();
sarParam.noiseVector = this->GetNoiseVectorGeom();
sarParam.burstRecords = this->GetBurstRecordsGeom();
const std::string supportDataPrefix = "support_data.";
sarParam.rangeSamplingRate = m_MetadataSupplierInterface->GetAs<double>(
supportDataPrefix + "range_sampling_rate");
sarParam.nearRangeTime = m_MetadataSupplierInterface->GetAs<double>(
supportDataPrefix + "slant_range_to_first_pixel");
sarParam.rangeResolution = m_MetadataSupplierInterface->GetAs<double>(
supportDataPrefix + "range_spacing");
sarParam.azimuthTimeInterval = boost::posix_time::precise_duration(m_MetadataSupplierInterface->GetAs<double>(
supportDataPrefix + "line_time_interval") * 1e6);
return true;
}
......
......@@ -35,6 +35,9 @@ public:
SarSensorModel(const ImageMetadata & imd);
virtual ~SarSensorModel() = default;
SarSensorModel(const SarSensorModel&) = delete; // non construction-copyable
SarSensorModel& operator=(const SarSensorModel&) = delete; // non copyable
using Point2DType = itk::Point<double, 2>;
using Point3DType = itk::Point<double, 3>;
......@@ -45,7 +48,8 @@ public:
/** Transform world point (lat,lon,hgt) to input image point
(col,row) */
void WorldToLineSample(const Point3DType& inGeoPoint, Point2DType& outLineSample) const;
void WorldToLineSample(const Point3DType& inGeoPoint,
Point2DType& outLineSample) const;
bool WorldToAzimuthRangeTime(const Point3DType& inGeoPoint,
TimeType & azimuthTime,
......@@ -53,6 +57,10 @@ public:
Point3DType& sensorPos,
Vector3DType& sensorVel) const;
void LineSampleHeightToWorld(const Point2DType& imPt,
double heightAboveEllipsoid,
Point3DType& worldPt) const;
protected:
private:
......@@ -88,16 +96,23 @@ private:
void AzimuthTimeToLine(const TimeType & azimuthTime,
double & line) const;
void SlantRangeToGroundRange(const double & slantRange,
void SlantRangeToGroundRange(double slantRange,
const TimeType & azimuthTime,
double & groundRange) const;
void ApplyCoordinateConversion(const double & in,
void ApplyCoordinateConversion(double in,
const TimeType& azimuthTime,
const std::vector<CoordinateConversionRecord> & records,
double & out) const;
const GCP & findClosestGCP(const Point2DType& imPt, const Projection::GCPParam & gcpParam) const;
Point3DType projToSurface(const GCP & gcp,
const Point2DType & imPt,
double heightAboveEllipsoid) const;
const ImageMetadata & m_Imd;
SARParam m_SarParam;
......@@ -105,7 +120,7 @@ private:
double m_RangeTimeOffset; // Offset in seconds
// Speed of light
const double C = 299792458;
static constexpr double C = 299792458;
// True if the input product is a ground product
bool m_IsGrd;
......
......@@ -134,6 +134,19 @@ bool SarSensorModel::WorldToAzimuthRangeTime(const Point3DType& inGeoPoint,
}
void SarSensorModel::LineSampleHeightToWorld(const Point2DType& imPt,
double heightAboveEllipsoid,
Point3DType& worldPt) const
{
assert(m_Imd.Has(MDGeom::GCP));
const auto& gcp = findClosestGCP(imPt, m_Imd.GetGCPParam());
Point3DType ecefPoint;
// Simple iterative inversion of inverse model starting at closest gcp
worldPt = EcefToWorld(projToSurface(gcp, imPt, heightAboveEllipsoid));
}
bool SarSensorModel::ZeroDopplerLookup(const Point3DType& inEcefPoint,
TimeType & azimuthTime,
......@@ -302,10 +315,8 @@ void SarSensorModel::OptimizeTimeOffsetsFromGcps()
}
DurationType cumulAzimuthTime(seconds(0));
double cumulRangeTime(0);
unsigned int count=0;
// reset offsets before optimisation
m_AzimuthTimeOffset = seconds(0);
auto gcpParam = m_Imd.GetGCPParam();
......@@ -345,7 +356,6 @@ void SarSensorModel::AzimuthTimeToLine(const TimeType & azimuthTime, double & li
auto currentBurst = m_SarParam.burstRecords.cbegin();
// Look for the correct burst. In most cases the number of burst
// records will be 1 (except for TOPSAR Sentinel1 products)
auto it = m_SarParam.burstRecords.cbegin();
......@@ -389,12 +399,12 @@ void SarSensorModel::AzimuthTimeToLine(const TimeType & azimuthTime, double & li
line = (timeSinceStart/m_SarParam.azimuthTimeInterval) + currentBurst->startLine;
}
void SarSensorModel::SlantRangeToGroundRange(const double & slantRange, const TimeType & azimuthTime, double & groundRange) const
void SarSensorModel::SlantRangeToGroundRange(double slantRange, const TimeType & azimuthTime, double & groundRange) const
{
ApplyCoordinateConversion(slantRange, azimuthTime, m_SarParam.slantRangeToGroundRangeRecords , groundRange);
}
void SarSensorModel::ApplyCoordinateConversion(const double & in,
void SarSensorModel::ApplyCoordinateConversion(double in,
const TimeType& azimuthTime,
const std::vector<CoordinateConversionRecord> & records,
double & out) const
......@@ -477,4 +487,152 @@ void SarSensorModel::ApplyCoordinateConversion(const double & in,
}
}
const GCP & SarSensorModel::findClosestGCP(const Point2DType& imPt, const Projection::GCPParam & gcpParam) const
{
assert(gcpParam.GCPs.size() > 0);
const GCP* closest;
// Squared distance between a Point and a gcp
auto squaredDistance = [](const Point2DType & imPt, const GCP & gcp)
{
double dx = imPt[0] - gcp.m_GCPCol;
double dy = imPt[1] - gcp.m_GCPRow;
return dx*dx + dy*dy;
};
double minDistance = std::numeric_limits<double>::max();
for (const auto & elem : gcpParam.GCPs)
{
auto currentDistance = squaredDistance(imPt, elem);
if (currentDistance < minDistance)
{
minDistance = currentDistance;
closest = &elem;
}
}
return *closest;
}
SarSensorModel::Point3DType SarSensorModel::projToSurface(const GCP & gcp, const Point2DType& imPt, double heightAboveEllipsoid) const
{
// Stop condition: img residual < 1e-2 pixels, height residual² <
// 0.01² m, nb iter < 50. the loop runs at least once.
const int maxIter = 50; //50
const double imgResidual = 1e-2;
const double heightResidual = 1e-2;
using MatrixType = itk::Matrix<double, 3, 3>;
using VectorType = itk::Vector<double, 3>;
Point3DType groundGcp;
groundGcp[0] = gcp.m_GCPX;
groundGcp[1] = gcp.m_GCPY;
groundGcp[2] = gcp.m_GCPZ;
// Initialize current estimation
auto currentEstimation = WorldToEcef(groundGcp);
// Corresponding image position
Point2DType currentImPoint;
currentImPoint[0] = gcp.m_GCPCol;
currentImPoint[1] = gcp.m_GCPRow;
auto currentImSquareResidual = imPt.SquaredEuclideanDistanceTo(currentImPoint);
double currentHeightResidual = 0.; //TODO
MatrixType B, BtB;
VectorType BtF;
VectorType F;
VectorType dR;
// Delta use for partial derivatives estimation (in meters)
const double d = 10.;
for (int i = 0; i < maxIter; i++)
{
// compute residuals
F[0] = imPt[0] - currentImPoint[0];
F[1] = imPt[1] - currentImPoint[1];
F[2] = currentHeightResidual;
// Compute partial derivatives
VectorType p_fx, p_fy, p_fh, dx(0.) ,dy(0.), dz(0.);
dx[0] = d;
dy[1] = d;
dz[2] = d;
Point2DType tmpImPt;
auto currentEstimationWorld = EcefToWorld(currentEstimation);
auto tmpGpt = EcefToWorld(currentEstimation + dx);
WorldToLineSample(tmpGpt, tmpImPt);
p_fx[0] = (currentImPoint[0] - tmpImPt[0])/d;
p_fy[0] = (currentImPoint[1] - tmpImPt[1])/d;
p_fh[0] = (currentEstimationWorld[2] - tmpGpt[2])/d;
tmpGpt = EcefToWorld(currentEstimation + dy);
WorldToLineSample(tmpGpt,tmpImPt);
p_fx[1] = (currentImPoint[0] - tmpImPt[0])/d;
p_fy[1] = (currentImPoint[1] - tmpImPt[1])/d;
p_fh[1] = (currentEstimationWorld[2] - tmpGpt[2])/d;
tmpGpt = EcefToWorld(currentEstimation + dz);
WorldToLineSample(tmpGpt,tmpImPt);
p_fx[2] = (currentImPoint[0] - tmpImPt[0])/d;
p_fy[2] = (currentImPoint[1] - tmpImPt[1])/d;
p_fh[2] = (currentEstimationWorld[2] - tmpGpt[2])/d;
B(0,0) = p_fx[0];
B(0,1) = p_fx[1];
B(0,2) = p_fx[2];
B(1,0) = p_fy[0];
B(1,1) = p_fy[1];
B(1,2) = p_fy[2];
B(2,0) = p_fh[0];
B(2,1) = p_fh[1];
B(2,2) = p_fh[2];
// Invert system
try
{
dR = MatrixType(B.GetInverse()) * F;
}
catch (itk::ExceptionObject const& e)
{
otbLogMacro(Warning, << "SarSensorModel::projToSurface(): singular matrix can not be inverted. Returning best estimation so far("
<< currentEstimation
<< ") for output point ("
<< imPt << ")" << std::endl);
return currentEstimation;
}
currentEstimation -= dR;
currentEstimationWorld = EcefToWorld(currentEstimation);
currentImSquareResidual = imPt.SquaredEuclideanDistanceTo(currentImPoint);
/* TODO manage elevation
const ossim_float64 atHgt = hgtRef.getRefHeight(currentEstimationWorld);
currentHeightResidual = atHgt - currentEstimationWorld.height();
*/
WorldToLineSample(currentEstimationWorld, currentImPoint);
if (currentImSquareResidual <= imgResidual && currentHeightResidual <= heightResidual)
{
break;
}
}
return currentEstimation;
}
} //namespace otb
\ No newline at end of file
......@@ -226,13 +226,19 @@ if (Boost_UNIT_TEST_FRAMEWORK_FOUND)
${OTBTransform-Test_LIBRARIES}
${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} )
otb_add_test(NAME prTvSarSensorModelSentinel1
otb_add_test(NAME prTvSarSensorModelSentinel1_SLC
COMMAND otbSarSensorModelTest --
#LARGEINPUT{SENTINEL1/S1A_S6_SLC__1SSV_20150619T195043/measurement/s1a-s6-slc-vv-20150619t195043-20150619t195101-006447-00887d-001.tiff}
LARGEINPUT{SENTINEL1/S1A_S6_SLC__1SSV_20150619T195043/manifest.safe}
-23.9431732671333
15.7637929806278
0)
otb_add_test(NAME prTvSarSensorModelSentinel1_geom
COMMAND otbSarSensorModelTest --
LARGEINPUT{SENTINEL1/S1A_S6_SLC__1SSV_20150619T195043/measurement/s1a-s6-slc-vv-20150619t195043-20150619t195101-006447-00887d-001.tiff}?&geom=${INPUTDATA}/sentinel1/sentinel1-2.geom
-23.9431732671333
15.7637929806278
0)
otb_add_test(NAME prTvSarSensorModelCosmoSkyMed
COMMAND otbSarSensorModelTest --
......@@ -241,11 +247,11 @@ if (Boost_UNIT_TEST_FRAMEWORK_FOUND)
43.5726757080402
0)
otb_add_test(NAME prTvSarSensorModelTerraSarX
COMMAND otbSarSensorModelTest --
LARGEINPUT{TERRASARX/2008-03-10_GrandCanyon_SSC/TSX1_SAR__SSC______SM_S_SRA_20080310T133220_20080310T133228}
-111.664738184032
36.2692959995967
0)
# otb_add_test(NAME prTvSarSensorModelTerraSarX
# COMMAND otbSarSensorModelTest --
# LARGEINPUT{TERRASARX/2008-03-10_GrandCanyon_SSC/TSX1_SAR__SSC______SM_S_SRA_20080310T133220_20080310T133228}
# -111.664738184032
# 36.2692959995967
# 0)
endif()
\ No newline at end of file
......@@ -37,6 +37,8 @@
//TODO: replace by baselines
#include "otbSarSensorModelAdapter.h"
#include "ossim/ossimSarSensorModel.h"
using namespace boost::unit_test;
......@@ -52,7 +54,7 @@ BOOST_AUTO_TEST_CASE(SARSensorModel_parameters)
BOOST_AUTO_TEST_CASE(SARSensorModel_WorldToLineSample)
{
using ImageType = otb::Image<unsigned int, 2>;
using ImageType = otb::VectorImage<unsigned int, 2>;
using ReaderType = otb::ImageFileReader<ImageType>;
auto reader = ReaderType::New();
......@@ -70,21 +72,9 @@ BOOST_AUTO_TEST_CASE(SARSensorModel_WorldToLineSample)
otb::SarSensorModel::Point2DType outLineSampleOssim;
std::cout << "test with ossim" << std::endl;
auto ossimModel = otb::SarSensorModelAdapter::New();
ossimModel->LoadState(reader->GetOutput()->GetImageKeywordlist());
ossimModel->WorldToLineSample(inWorldPoint, outLineSampleOssim);
std::cout << inWorldPoint << " " << outLineSampleOssim << std::endl;
std::cout << "test with otb" << std::endl;
std::cout << "Direct test with otb" << std::endl;
otb::SarSensorModel::Point2DType outLineSampleOtb;
model.WorldToLineSample(inWorldPoint, outLineSampleOtb);
std::cout << inWorldPoint << " " << outLineSampleOtb << std::endl;
//TODO compare results with a baseline
}
......@@ -93,7 +83,7 @@ BOOST_AUTO_TEST_CASE(SARSensorModel_auto_validate_inverse_transform )
double lineTol = 1.;
double sampleTol = 1.;
using ImageType = otb::Image<unsigned int, 2>;
using ImageType = otb::VectorImage<unsigned int, 2>;
using ReaderType = otb::ImageFileReader<ImageType>;
auto reader = ReaderType::New();
......@@ -128,5 +118,69 @@ BOOST_AUTO_TEST_CASE(SARSensorModel_auto_validate_inverse_transform )
BOOST_TEST(std::abs(lineSample[0] - lineSampleBaseline[0]) < lineTol);
BOOST_TEST(std::abs(lineSample[1] - lineSampleBaseline[1]) < sampleTol);
}
}
BOOST_AUTO_TEST_CASE(SARSensorModel_auto_validate_forward_transform)
{
double sqResTol = 25;
using ImageType = otb::VectorImage<unsigned int, 2>;
using ReaderType = otb::ImageFileReader<ImageType>;
auto reader = ReaderType::New();
reader->SetFileName(framework::master_test_suite().argv[1]);
reader->GenerateOutputInformation();
if (reader->GetOutput()->GetGCPCount() == 0)
{
otbLogMacro(Info, << "Input product has no gcp, skipping gcp forward transform validation.");
return;
}
auto & imd = reader->GetOutput()->m_Imd;
auto & GCPParam = imd.GetGCPParam();
std::vector<otb::GCP> testGCPs;
otb::Projection::GCPParam productGCPs;
bool odd = false;
for (auto gcp: GCPParam.GCPs)
{
gcp.m_GCPZ = 0.;
if (odd)
{
productGCPs.GCPs.push_back(gcp);
}
else
{
testGCPs.push_back(gcp);
}
odd = !odd;
}
imd.Add(otb::MDGeom::GCP, productGCPs);
otb::SarSensorModel model(imd);
for (const auto & gcp : testGCPs)
{
itk::Point<double, 2> sensorPoint;
sensorPoint[0] = gcp.m_GCPCol;
sensorPoint[1] = gcp.m_GCPRow;
itk::Point<double, 3> geoPointBaseline;
geoPointBaseline[0] = gcp.m_GCPX;
geoPointBaseline[1] = gcp.m_GCPY;
geoPointBaseline[2] = gcp.m_GCPZ;
itk::Point<double, 3> geoPoint;
model.LineSampleHeightToWorld(sensorPoint, gcp.m_GCPZ, geoPoint);
BOOST_TEST(geoPoint.SquaredEuclideanDistanceTo(geoPointBaseline) < sqResTol);
}
}
\ No newline at end of file
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