diff --git a/Modules/Adapters/OSSIMAdapters/include/otbSarSensorModelAdapter.h b/Modules/Adapters/OSSIMAdapters/include/otbSarSensorModelAdapter.h
index 692c1574fc7e12b390832546a43154411d6f5e7b..8a713883b6720d0dd9b8e0ce89c03ae924cd75c9 100644
--- a/Modules/Adapters/OSSIMAdapters/include/otbSarSensorModelAdapter.h
+++ b/Modules/Adapters/OSSIMAdapters/include/otbSarSensorModelAdapter.h
@@ -24,6 +24,7 @@
#include <memory>
#include "otbDEMHandler.h"
+#include "itkPoint.h"
namespace ossimplugins
{
@@ -62,6 +63,9 @@ public:
typedef std::auto_ptr<ossimplugins::ossimSarSensorModel> InternalModelPointer;
+ using Point2DType = itk::Point<double,2>;
+ using Point3DType = itk::Point<double,3>;
+
/** Method for creation through the object factory. */
itkNewMacro(Self);
@@ -80,6 +84,21 @@ public:
/** Deburst metadata if possible and return lines to keep in image file */
bool Deburst(std::vector<std::pair<unsigned long, unsigned long> > & lines);
+ /** Transform world point (lat,lon,hgt) to input image point
+ (col,row) and YZ frame */
+ bool WorldToLineSampleYZ(const Point3DType & inGeoPoint, Point2DType & cr, Point2DType & yz) const;
+
+ /** Transform world point (lat,lon,hgt) to input image point
+ (col,row) */
+ bool WorldToLineSample(const Point3DType & inGEoPOint, Point2DType & cr) const;
+
+/** Transform world point (lat,lon,hgt) to satellite position (x,y,z) and satellite velocity */
+ bool WorldToSatPositionAndVelocity(const Point3DType & inGeoPoint, Point3DType & satelitePosition,
+ Point3DType & sateliteVelocity) const;
+
+ /** Transform world point (lat,lon,hgt) to cartesian point (x,y,z) */
+ static bool WorldToCartesian(const Point3DType & inGeoPoint, Point3DType & outCartesianPoint);
+
static bool ImageLineToDeburstLine(const std::vector<std::pair<unsigned long,unsigned long> >& lines, unsigned long imageLine, unsigned long & deburstLine);
static void DeburstLineToImageLine(const std::vector<std::pair<unsigned long,unsigned long> >& lines, unsigned long deburstLine, unsigned long & imageLine);
diff --git a/Modules/Adapters/OSSIMAdapters/src/otbSarSensorModelAdapter.cxx b/Modules/Adapters/OSSIMAdapters/src/otbSarSensorModelAdapter.cxx
index 62d5203af5c4d7dc04f4b9a3d66706cef1fa084a..56bcfbca01f498b1f8cc020b5c8e839ebf809cbb 100644
--- a/Modules/Adapters/OSSIMAdapters/src/otbSarSensorModelAdapter.cxx
+++ b/Modules/Adapters/OSSIMAdapters/src/otbSarSensorModelAdapter.cxx
@@ -109,7 +109,112 @@ void SarSensorModelAdapter::DeburstLineToImageLine(const std::vector<std::pair<u
ossimplugins::ossimSarSensorModel::deburstLineToImageLine(lines,deburstLine,imageLine);
}
+bool SarSensorModelAdapter::WorldToLineSampleYZ(const Point3DType & inGeoPoint, Point2DType & cr, Point2DType & yz) const
+{
+ if(m_SensorModel.get() == ITK_NULLPTR)
+ {
+ return false;
+ }
+
+ ossimGpt inGpt;
+ inGpt.lon = inGeoPoint[0];
+ inGpt.lat = inGeoPoint[1];
+ inGpt.hgt = inGeoPoint[2];
+
+ ossimDpt outDpt;
+
+ double y(0.),z(0.);
+ m_SensorModel->worldToLineSampleYZ(inGpt,outDpt,y,z);
+
+ if(outDpt.isNan())
+ return false;
+
+ cr[0]=outDpt.x;
+ cr[1]=outDpt.y;
+ yz[0]=y;
+ yz[1]=z;
+
+ return true;
+}
+bool SarSensorModelAdapter::WorldToLineSample(const Point3DType & inGeoPoint, Point2DType & cr) const
+{
+ if(m_SensorModel.get() == ITK_NULLPTR)
+ {
+ return false;
+ }
+
+ ossimGpt inGpt;
+ inGpt.lon = inGeoPoint[0];
+ inGpt.lat = inGeoPoint[1];
+ inGpt.hgt = inGeoPoint[2];
+
+ ossimDpt outDpt;
+
+ m_SensorModel->worldToLineSample(inGpt,outDpt);
+
+ if(outDpt.isNan())
+ return false;
+
+ cr[0]=outDpt.x;
+ cr[1]=outDpt.y;
+
+ return true;
+}
+
+bool SarSensorModelAdapter::WorldToCartesian(const Point3DType & inGeoPoint, Point3DType & outCartesianPoint)
+{
+ ossimGpt inGpt;
+ inGpt.lon = inGeoPoint[0];
+ inGpt.lat = inGeoPoint[1];
+ inGpt.hgt = inGeoPoint[2];
+
+ ossimEcefPoint outCartesien(inGpt);
+
+ if(outCartesien.isNan())
+ return false;
+
+ outCartesianPoint[0] = outCartesien.x();
+ outCartesianPoint[1] = outCartesien.y();
+ outCartesianPoint[2] = outCartesien.z();
+
+ return true;
+}
+
+bool SarSensorModelAdapter::WorldToSatPositionAndVelocity(const Point3DType & inGeoPoint,
+ Point3DType & satelitePosition,
+ Point3DType & sateliteVelocity) const
+{
+ if(m_SensorModel.get() == ITK_NULLPTR)
+ {
+ return false;
+ }
+
+ ossimGpt inGpt;
+ inGpt.lon = inGeoPoint[0];
+ inGpt.lat = inGeoPoint[1];
+ inGpt.hgt = inGeoPoint[2];
+
+ ossimplugins::ossimSarSensorModel::TimeType azimuthTime;
+ double rangeTime;
+ ossimEcefPoint sensorPos;
+ ossimEcefVector sensorVel;
+
+ const bool success = m_SensorModel->worldToAzimuthRangeTime(inGpt, azimuthTime, rangeTime,sensorPos,sensorVel);
+
+ if(sensorPos.isNan() || !success)
+ return false;
+
+ satelitePosition[0] = sensorPos.x();
+ satelitePosition[1] = sensorPos.y();
+ satelitePosition[2] = sensorPos.z();
+
+ sateliteVelocity[0] = sensorVel.x();
+ sateliteVelocity[1] = sensorVel.y();
+ sateliteVelocity[2] = sensorVel.z();
+
+ return true;
+}
} // namespace otb
diff --git a/Modules/Adapters/OSSIMAdapters/test/CMakeLists.txt b/Modules/Adapters/OSSIMAdapters/test/CMakeLists.txt
index 54a7117fc93c342e11d0e77b2b5226c70c0104cc..8dadc14ff05f911b3fa0af4cca37e05592cecced 100644
--- a/Modules/Adapters/OSSIMAdapters/test/CMakeLists.txt
+++ b/Modules/Adapters/OSSIMAdapters/test/CMakeLists.txt
@@ -31,6 +31,7 @@ otbGeometricSarSensorModelAdapter.cxx
otbPlatformPositionAdapter.cxx
otbDEMHandlerTest.cxx
otbRPCSolverAdapterTest.cxx
+otbSarSensorModelAdapterTest.cxx
)
add_executable(otbOSSIMAdaptersTestDriver ${OTBOSSIMAdaptersTests})
@@ -499,3 +500,8 @@ otb_add_test(NAME uaTvRPCSolverAdapterValidationTest COMMAND otbOSSIMAdaptersTes
${INPUTDATA}/DEM/egm96.grd
)
+
+otb_add_test(NAME uaTvSarSensorModelAdapter COMMAND otbOSSIMAdaptersTestDriver
+ otbSarSensorModelAdapterTest
+ ${INPUTDATA}/s1a-iw1-slc-vh-amp_xt.geom
+ )
diff --git a/Modules/Adapters/OSSIMAdapters/test/otbOSSIMAdaptersTestDriver.cxx b/Modules/Adapters/OSSIMAdapters/test/otbOSSIMAdaptersTestDriver.cxx
index 911e57e671884e9bb003e399b605b6dfd3dacf29..5f7295d6dffdb208aafde1ab1326ecdd544bd24f 100644
--- a/Modules/Adapters/OSSIMAdapters/test/otbOSSIMAdaptersTestDriver.cxx
+++ b/Modules/Adapters/OSSIMAdapters/test/otbOSSIMAdaptersTestDriver.cxx
@@ -33,4 +33,5 @@ void RegisterTests()
REGISTER_TEST(otbPlatformPositionComputeBaselineTest);
REGISTER_TEST(otbDEMHandlerTest);
REGISTER_TEST(otbRPCSolverAdapterTest);
+ REGISTER_TEST(otbSarSensorModelAdapterTest);
}
diff --git a/Modules/Adapters/OSSIMAdapters/test/otbSarSensorModelAdapterTest.cxx b/Modules/Adapters/OSSIMAdapters/test/otbSarSensorModelAdapterTest.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d48f10066193796158002b55aba127b464814ccc
--- /dev/null
+++ b/Modules/Adapters/OSSIMAdapters/test/otbSarSensorModelAdapterTest.cxx
@@ -0,0 +1,84 @@
+/*
+ * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ * https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+
+#include <fstream>
+#include <iomanip>
+
+#include "otbSarSensorModelAdapter.h"
+#include "otbImageKeywordlist.h"
+
+
+int otbSarSensorModelAdapterTest(int itkNotUsed(argc), char* argv[])
+{
+ std::string infname = argv[1];
+
+ otb::SarSensorModelAdapter::Pointer sensorModel = otb::SarSensorModelAdapter::New();
+
+ auto kwl = otb::ReadGeometryFromGEOMFile(infname);
+
+ bool success = sensorModel->LoadState(kwl);
+
+ if(!success)
+ {
+ std::cerr<<"Could not LoadState() from keyword list read from"<<infname<<std::endl;
+ return EXIT_FAILURE;
+ }
+
+ std::vector<std::pair<unsigned long, unsigned long> > lines;
+ success = sensorModel->Deburst(lines);
+
+ if(!success)
+ {
+ std::cerr<<"Deburst() call failed."<<std::endl;
+ return EXIT_FAILURE;
+ }
+
+
+ otb::ImageKeywordlist outKwl;
+ success = sensorModel->SaveState(outKwl);
+ if(!success)
+ {
+ std::cerr<<"SaveState() call failed."<<std::endl;
+ return EXIT_FAILURE;
+ }
+
+
+ otb::SarSensorModelAdapter::Point2DType out1,out2;
+ otb::SarSensorModelAdapter::Point3DType in, out3, out4, out5;
+
+ // GCP 99 from input geom file
+ //support_data.geom.gcp[99].world_pt.hgt: 2.238244926818182e+02
+ //support_data.geom.gcp[99].world_pt.lat: 4.323458093295080e+01
+ //support_data.geom.gcp[99].world_pt.lon: 1.116316013091967e+00
+
+ in[0] = 4.323458093295080e+01;
+ in[1] = 1.116316013091967e+00;
+ in[2] = 2.238244926818182e+02;
+
+ sensorModel->WorldToLineSample(in,out1);
+ sensorModel->WorldToLineSampleYZ(in,out1,out2);
+
+ sensorModel->WorldToCartesian(in, out5);
+ sensorModel->WorldToSatPositionAndVelocity(in,out3, out4);
+
+
+ return EXIT_SUCCESS;
+}
diff --git a/Modules/ThirdParty/OssimPlugins/include/ossim/ossimSarSensorModel.h b/Modules/ThirdParty/OssimPlugins/include/ossim/ossimSarSensorModel.h
index d76968ecf1929a8f8034e016c8b13fa2d2ff3bfb..542b29e6f8c65dd98cd919bcf0f363cc8072c193 100644
--- a/Modules/ThirdParty/OssimPlugins/include/ossim/ossimSarSensorModel.h
+++ b/Modules/ThirdParty/OssimPlugins/include/ossim/ossimSarSensorModel.h
@@ -238,13 +238,36 @@ public:
* \param[in] worldPoint World point to geocode
* \param[out] azimuthTime Estimated zero-doppler azimuth time
* \param[out] rangeTime Estimated range time
+ * \param[out] interpSensorPos interpolated ECEF sensor position
+ * \param[out] interpSensorVel interpolated ECEF sensor velocity
* \return True if success, false otherwise. In this case,
* azimuthTime and rangeTime will not be modified.
*/
- /*virtual*/ bool worldToAzimuthRangeTime(const ossimGpt& worldPt, TimeType & azimuthTime, double & rangeTime) const;
+ /*virtual*/ bool worldToAzimuthRangeTime(const ossimGpt& worldPt, TimeType & azimuthTime, double & rangeTime, ossimEcefPoint & interpSensorPos, ossimEcefVector & interpSensorVel) const;
- // TODO: document me
- /*virtual*/ void lineSampleToAzimuthRangeTime(const ossimDpt & imPt, TimeType & azimuthTime, double & rangeTime) const;
+ /**
+ * This method implement inverse sar geolocation similar to
+ * worldToLineSample, except that it also returns (y,z)
+ * coordinates, defined as follows:
+ * Let n = |sensorPos|,
+ * ps2 = scalar_product(sensorPos,worldPoint)
+ * d = distance(sensorPos,worldPoint)
+ *
+ * z = n - ps2/n
+ * y = sqrt(d*d - z*z)
+ *
+ * sign of y is furher adapted to be inverted if sensor is left or
+ * right looking
+ *
+ * \param[in] worldPoint World point to geocode
+ * \param[out] imPt Corresponding estimated image point
+ * \param[out] y
+ * \param[out] z
+ * \return True if success, false otherwise. In this case,
+ */
+ void worldToLineSampleYZ(const ossimGpt& worldPt, ossimDpt& imPt, double & y, double & z) const;
+
+ void lineSampleToAzimuthRangeTime(const ossimDpt & imPt, TimeType & azimuthTime, double & rangeTime) const;
// TODO: document me
bool autovalidateInverseModelFromGCPs(const double & xtol = 1, const double & ytol = 1, const double azTimeTol = 500, const double &rangeTimeTo=0.0000000001) const;
@@ -419,7 +442,8 @@ protected:
ProductType theProductType; // GRD/SLC
DurationType theAzimuthTimeOffset; // Offset computed
double theRangeTimeOffset; // Offset in seconds, computed
-
+ bool theRightLookingFlag;
+
static const double C;
static const unsigned int thePluginVersion; // version of the SarSensorModel plugin
diff --git a/Modules/ThirdParty/OssimPlugins/src/ossim/ossimSarSensorModel.cpp b/Modules/ThirdParty/OssimPlugins/src/ossim/ossimSarSensorModel.cpp
index ac695a2a4b790ebce7671fc7e27c62a954e91098..bb8903ad981a1857fbedf1311cf626138c5bd035 100644
--- a/Modules/ThirdParty/OssimPlugins/src/ossim/ossimSarSensorModel.cpp
+++ b/Modules/ThirdParty/OssimPlugins/src/ossim/ossimSarSensorModel.cpp
@@ -12,7 +12,7 @@
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
- *
+
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -137,7 +137,8 @@ namespace ossimplugins
theRangeResolution(0.),
theBistaticCorrectionNeeded(false),
theAzimuthTimeOffset(seconds(0)),
- theRangeTimeOffset(0.)
+ theRangeTimeOffset(0.),
+ theRightLookingFlag(true)
{}
ossimSarSensorModel::GCPRecordType const&
@@ -220,7 +221,58 @@ namespace ossimplugins
TimeType azimuthTime;
double rangeTime;
- const bool success = worldToAzimuthRangeTime(worldPt, azimuthTime, rangeTime);
+ ossimEcefPoint sensorPos;
+ ossimEcefVector sensorVel;
+
+ const bool success = worldToAzimuthRangeTime(worldPt, azimuthTime, rangeTime, sensorPos, sensorVel);
+
+ if(!success)
+ {
+ imPt.makeNan();
+ return;
+ }
+ // std::clog << "AzimuthTime: " << azimuthTime << "\n";
+ // std::clog << "RangeTime: " << rangeTime << "\n";
+ // std::clog << "GRD: " << isGRD() << "\n";
+
+ // Convert azimuth time to line
+ azimuthTimeToLine(azimuthTime,imPt.y);
+
+ if(isGRD())
+ {
+ // GRD case
+ double groundRange(0);
+ slantRangeToGroundRange(rangeTime*C/2,azimuthTime,groundRange);
+ // std::clog << "GroundRange: " << groundRange << "\n";
+ // std::clog << "TheRangeResolution: " << theRangeResolution << "\n";
+
+ // Eq 32 p. 31
+ // TODO: possible micro-optimization: precompute 1/theRangeResolution, and
+ // use *
+ imPt.x = groundRange/theRangeResolution;
+ }
+ else
+ {
+ // std::clog << "TheNearRangeTime: " << theNearRangeTime << "\n";
+ // std::clog << "TheRangeSamplingRate: " << theRangeSamplingRate << "\n";
+ // SLC case
+ // Eq 23 and 24 p. 28
+ imPt.x = (rangeTime - theNearRangeTime)*theRangeSamplingRate;
+ }
+ }
+
+void ossimSarSensorModel::worldToLineSampleYZ(const ossimGpt& worldPt, ossimDpt & imPt, double & y, double & z) const
+ {
+ // std::clog << "ossimSarSensorModel::worldToLineSample()\n";
+ assert(theRangeResolution>0&&"theRangeResolution is null.");
+
+ // First compute azimuth and range time
+ TimeType azimuthTime;
+ double rangeTime;
+ ossimEcefPoint sensorPos;
+ ossimEcefVector sensorVel;
+
+ const bool success = worldToAzimuthRangeTime(worldPt, azimuthTime, rangeTime,sensorPos,sensorVel);
if(!success)
{
@@ -255,9 +307,35 @@ namespace ossimplugins
// Eq 23 and 24 p. 28
imPt.x = (rangeTime - theNearRangeTime)*theRangeSamplingRate;
}
+
+ // Now computes Y and Z
+ ossimEcefPoint inputPt(worldPt);
+ double NormeS = sqrt(sensorPos[0]*sensorPos[0] + sensorPos[1]*sensorPos[1] + sensorPos[2]*sensorPos[2]); /* distance du radar */
+ double PS2 = inputPt[0]*sensorPos[0] + inputPt[1]*sensorPos[1] + inputPt[2]*sensorPos[2];
+
+ // TODO check for small NormesS to avoid division by zero ?
+ // Should never happen ...
+ assert(NormesS>1e-6);
+ z = NormeS - PS2/NormeS;
+
+ double distance = sqrt((sensorPos[0]-inputPt[0])*(sensorPos[0]-inputPt[0]) +
+ (sensorPos[1]-inputPt[1])*(sensorPos[1]-inputPt[1]) +
+ (sensorPos[2]-inputPt[2])*(sensorPos[2]-inputPt[2]));
+
+ y = sqrt(distance*distance - z*z);
+
+ // Check view side and change sign of Y accordingly
+ if ( (( sensorVel[0] * (sensorPos[1]* inputPt[2] - sensorPos[2]* inputPt[1]) +
+ sensorVel[1] * (sensorPos[2]* inputPt[0] - sensorPos[0]* inputPt[2]) +
+ sensorVel[2] * (sensorPos[0]* inputPt[1] - sensorPos[1]* inputPt[0])) > 0) ^ theRightLookingFlag )
+ {
+ y = -y;
+ }
}
- bool ossimSarSensorModel::worldToAzimuthRangeTime(const ossimGpt& worldPt, TimeType & azimuthTime, double & rangeTime) const
+
+
+bool ossimSarSensorModel::worldToAzimuthRangeTime(const ossimGpt& worldPt, TimeType & azimuthTime, double & rangeTime, ossimEcefPoint & interpSensorPos, ossimEcefVector & interpSensorVel) const
{
// std::clog << "ossimSarSensorModel::worldToAzimuthRangeTime()\n";
// First convert lat/lon to ECEF
@@ -265,8 +343,6 @@ namespace ossimplugins
// Compute zero doppler time
TimeType interpTime;
- ossimEcefPoint interpSensorPos;
- ossimEcefVector interpSensorVel;
const bool success = zeroDopplerLookup(inputPt,azimuthTime,interpSensorPos,interpSensorVel);
@@ -857,7 +933,9 @@ namespace ossimplugins
double estimatedRangeTime;
// Estimate times
- const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime);
+ ossimEcefPoint sensorPos;
+ ossimEcefVector sensorVel;
+ const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime,sensorPos, sensorVel);
this->worldToLineSample(gcpIt->worldPt,estimatedImPt);
const bool thisSuccess
@@ -964,8 +1042,10 @@ namespace ossimplugins
TimeType estimatedAzimuthTime;
double estimatedRangeTime;
+ ossimEcefPoint sensorPos;
+ ossimEcefVector sensorVel;
// Estimate times
- const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime);
+ const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime, sensorPos, sensorVel);
if(s1)
{
@@ -985,8 +1065,11 @@ namespace ossimplugins
TimeType estimatedAzimuthTime;
double estimatedRangeTime;
+ ossimEcefPoint sensorPos;
+ ossimEcefVector sensorVel;
+
// Estimate times
- const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime);
+ const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime, sensorPos, sensorVel);
if(s1)
{