Commit e35280b2 authored by Luc Hermitte's avatar Luc Hermitte
Browse files

ENH: merge code from github for New SAR Sensor Model (OTB-999)

parent 9c6917e9
......@@ -2,6 +2,7 @@ set(DOCUMENTATION "This module provides Ossim Plugins")
otb_module(OTBOssimPlugins
DEPENDS
OTBBoost
OTBGeoTIFF
OTBGDAL
OTBOssim
......
......@@ -50,7 +50,9 @@ set(ossimplugins_SOURCES
)
add_library(otbossimplugins ${ossimplugins_LIBTYPE} ${ossimplugins_SOURCES})
message(STATUS "lib for ossim: ${OTBBoost_LIBRARIES}")
target_link_libraries(otbossimplugins
${OTBBoost_LIBRARIES}
${OTBGDAL_LIBRARIES}
${OTBOssim_LIBRARIES}
${OTBGeoTIFF_LIBRARIES}
......
//----------------------------------------------------------------------------
//
// "Copyright Centre National d'Etudes Spatiales"
//
// License: LGPL-2
//
// See LICENSE.txt file in the top level directory for more details.
//
//----------------------------------------------------------------------------
// $Id$
#ifndef ossimSarSensorModel_HEADER
#define ossimSarSensorModel_HEADER
#include <boost/config.hpp>
#if defined(__GNUC__) || defined(__clang__)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-parameter"
// #pragma GCC diagnostic ignored "-Woverloaded-virtual"
// #pragma GCC diagnostic ignored "-Wshadow"
#include <ossim/projection/ossimSensorModel.h>
#include <ossim/elevation/ossimHgtRef.h>
#pragma GCC diagnostic pop
#else
#include <ossim/projection/ossimSensorModel.h>
#include <ossim/elevation/ossimHgtRef.h>
#endif
#include <boost/date_time/posix_time/posix_time.hpp>
namespace ossimplugins
{
class ossimSarSensorModel : public ossimSensorModel
{
public:
typedef boost::posix_time::ptime TimeType;
typedef boost::posix_time::time_duration DurationType;
struct OrbitRecordType
{
TimeType azimuthTime;
ossimEcefPoint position;
ossimEcefVector velocity;
};
struct GCPRecordType
{
TimeType azimuthTime;
double slantRangeTime;
ossimDpt imPt;
ossimGpt worldPt;
};
struct BurstRecordType
{
TimeType azimuthStartTime;
unsigned long startLine;
TimeType azimuthStopTime;
unsigned long endLine;
};
struct CoordinateConversionRecordType
{
TimeType azimuthTime;
double rg0;
std::vector<double> coefs;
};
/** Constructor */
ossimSarSensorModel();
#if ! (defined(BOOST_NO_DEFAULTED_FUNCTIONS) || defined(BOOST_NO_CXX1_DEFAULTED_FUNCTIONS))
/** Copy constructor */
ossimSarSensorModel(ossimSarSensorModel const& m) = default;
/** Move constructor */
ossimSarSensorModel(ossimSarSensorModel && m) = default;
/** Destructor */
virtual ~ossimSarSensorModel() = default;
#endif
virtual void lineSampleHeightToWorld(const ossimDpt& imPt, const double & heightEllipsoid, ossimGpt& worldPt) const;
virtual void lineSampleToWorld(const ossimDpt& imPt, ossimGpt& worldPt) const;
/** This method implement inverse sar geolocation using method found
* in ESA document "Guide to ASAR geocoding" (ref
* RSL-ASAR-GC-AD). Equation numbers can be found in source code
* comments.
*
* \param[in] worldPt World point to geocode
* \param[out] imPt Corresponding estimated image point
*/
virtual void worldToLineSample(const ossimGpt& worldPt, ossimDpt & imPt) const;
/**
* Sub-routine of lineSampleToWorld that computes azimuthTime and
* slant range time from worldPoint
*
* \param[in] worldPoint World point to geocode
* \param[out] azimuthTime Estimated zero-doppler azimuth time
* \param[out] rangeTime Estimated range time
* \return True if sucess, false otherwise. In this case,
* azimuthTime and rangeTime will not be modified.
*/
/*virtual*/ bool worldToAzimuthRangeTime(const ossimGpt& worldPt, TimeType & azimuthTime, double & rangeTime) const;
// TODO: document me
/*virtual*/ 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;
// TODO: document me
bool autovalidateForwardModelFromGCPs(double resTol = 10);
//Pure virtual in base class
bool useForward() const;
void optimizeTimeOffsetsFromGcps();
/**
* Returns pointer to a new instance, copy of this.
*/
virtual ossimObject* dup() const;
//TODO: Add virtual method readAnnotationFile?
protected:
/**
* Compute range and doppler frequency from an input point, sensor
* position and velocity.
*
* \param[in] inputPt The target point
* \param[in] sensorPos The sensor position
* \param[in] sensorvel The sensor velocity
* \param[out] range Estimated range
* \param[out] doppler Estimated doppler frequency
*/
virtual void computeRangeDoppler(const ossimEcefPoint & inputPt, const ossimEcefPoint & sensorPos, const ossimEcefVector sensorVel, double & range, double & doppler) const;
/**
* Interpolate sensor position and velocity at given azimuth time
* using lagragian interpolation of orbital records.
*
* \param[in] azimuthTime The time at which to interpolate
* \param[out] sensorPos Interpolated sensor position
* \param[out] sensorvel Interpolated sensor velocity
* \param[in] deg Degree of lagragian interpolation
*/
/*virtual*/ void interpolateSensorPosVel(const TimeType & azimuthTime, ossimEcefPoint& sensorPos, ossimEcefVector& sensorVel, unsigned int deg = 8) const;
/**
* Convert slant range to ground range by interpolating slant range
* to ground range coefficients.
*
* \param[in] slantRangeTime The slantRange to convert (meters)
* \param[in] azimuthTime The corresponding azimuth time
* \param[out] groundRange The estimated ground range (meters)
*/
/*virtual*/ void slantRangeToGroundRange(const double & slantRange, const TimeType & azimuthTime, double & groundRange) const;
// TODO: Document me
/*virtual*/ void groundRangeToSlantRange(const double & groundRange, const TimeType & azimuthTime, double & slantRange) const;
// TODO: Document me
/*virtual*/ void applyCoordinateConversion(const double & in, const TimeType& azimuthTime, const std::vector<CoordinateConversionRecordType> & records, double & out) const;
/**
* Estimate the zero-doppler azimuth time and corresponding sensor
* position and velocity from the inputPt.
*
* \param[in] inputPt The point to estimated zero-doppler time on
* \param[out] interpAzimuthTime Interpolated azimuth time
* \param[out] interpSensorPos Interpolated sensor position
* \param[out] interpSensorVel Interpolated sensor velocity
* \return True if success, false otherwise. In this case, output
* parameters are left untouched.
*/
/*virtual*/ bool zeroDopplerLookup(const ossimEcefPoint & inputPt, TimeType & interpAzimuthTime, ossimEcefPoint & interpSensorPos, ossimEcefVector & interpSensorVel) const;
/**
* Compute the bistatic correction to apply.
*
* \param[in] inputPt The point to compute bistatic correction on
* \param[in] sensorPos The corresponding sensor position
* \param[out] bistaticCorrection The estimated bistatic correction
*/
/*virtual*/ void computeBistaticCorrection(const ossimEcefPoint & inputPt, const ossimEcefPoint & sensorPos, DurationType & bistaticCorrection) const;
/**
* Convert azimuth time to fractional line.
*
* \param[in] azimuthTime The azimuth time to convert
* \param[out] The estimated fractional line
*/
/*virtual*/ // TODO: check why virtual
void azimuthTimeToLine(const TimeType & azimuthTime, double & line) const;
// TODO: document me
/*virtual*/ // TODO: check why virtual
void lineToAzimuthTime(const double & line, TimeType & azimuthTime) const;
// TODO: document me
/*virtual*/ // TODO: check why virtual
bool projToSurface(const GCPRecordType & initGcp, const ossimDpt & target, const ossimHgtRef & hgtRef, ossimEcefPoint & ellpt) const;
/**
* Finds closest GCP.
*
* \param[in] imPt «imPt-explanations»
* \return the closest GCP record to \c imPt.
* \throw None
* \pre `theGCPRecords` shall not be empty.
*/
GCPRecordType const& findClosestGCP(ossimDpt const& imPt) const;
std::vector<OrbitRecordType> theOrbitRecords;
std::vector<GCPRecordType> theGCPRecords;
std::vector<BurstRecordType> theBurstRecords;
std::vector<CoordinateConversionRecordType> theSlantRangeToGroundRangeRecords;
std::vector<CoordinateConversionRecordType> theGroundRangeToSlantRangeRecords;
double theRadarFrequency; // in Hz
double theAzimuthTimeInterval; // in microseconds
double theNearRangeTime; // in seconds
double theRangeSamplingRate; // in Hz
double theRangeResolution; // in meters
bool theBistaticCorrectionNeeded; // Do we need to compute
// bistatic correction ?
bool isGRD; // True if the product is GRD. False if it is SLC
double theAzimuthTimeOffset; // Offset in microseconds
double theRangeTimeOffset; // Offset in seconds;
static const double C;
private:
/** Disabled assignment operator. */
ossimSarSensorModel& operator=(ossimSarSensorModel const& rhs);
};
}
#endif
//----------------------------------------------------------------------------
//
// "Copyright Centre National d'Etudes Spatiales"
//
// License: LGPL-2
//
// See LICENSE.txt file in the top level directory for more details.
//
//----------------------------------------------------------------------------
// $Id$
#include <ossimSentinel1SarSensorModel.h>
#include <ossim/base/ossimXmlDocument.h>
#include <ossimXmlTools.h>
namespace {// Anonymous namespace
const ossimString attAzimuthTime = "azimuthTime";
const ossimString attFirstValidSample = "firstValidSample";
const ossimString attGr0 = "gr0";
const ossimString attGrsrCoefficients = "grsrCoefficients";
const ossimString attHeight = "height";
const ossimString attLatitude = "latitude";
const ossimString attLine = "line";
const ossimString attLongitude = "longitude";
const ossimString attPixel = "pixel";
const ossimString attPosition = "position";
const ossimString attSlantRangeTime = "slantRangeTime";
const ossimString attSr0 = "sr0";
const ossimString attSrgrCoefficients = "srgrCoefficients";
const ossimString attTime = "time";
const ossimString attVelocity = "velocity";
const ossimString attX = "x";
const ossimString attY = "y";
const ossimString attZ = "z";
}// Anonymous namespace
void ossimplugins::ossimSentinel1SarSensorModel::readCoordinates(
ossimXmlDocument const& xmlDoc, ossimString const& xpath,
ossimString const& rg0_xpath, ossimString const& coeffs_xpath,
std::vector<CoordinateConversionRecordType> & outputRecords)
{
std::vector<ossimRefPtr<ossimXmlNode> > xnodes;
xmlDoc.findNodes(xpath, xnodes);
for(std::vector<ossimRefPtr<ossimXmlNode> >::iterator itNode = xnodes.begin(); itNode!=xnodes.end();++itNode)
{
CoordinateConversionRecordType coordRecord;
coordRecord.azimuthTime = getTimeFromFirstNode(**itNode, attAzimuthTime);
coordRecord.rg0 = getDoubleFromFirstNode(**itNode, rg0_xpath);;
ossimString const& s = getTextFromFirstNode(**itNode, coeffs_xpath);
std::vector<ossimString> ssplit = s.split(" ");
if (ssplit.empty())
{
throw std::runtime_error("The "+rg0_xpath+" record has an empty coef vector");
}
for (std::vector<ossimString>::const_iterator cIt = ssplit.begin(), e = ssplit.end()
; cIt != e
; ++cIt
)
{
coordRecord.coefs.push_back(cIt->toDouble());
}
assert(!coordRecord.coefs.empty()&&"The rg0 record has empty coefs vector.");
outputRecords.push_back(coordRecord);
}
}
namespace ossimplugins
{
void ossimSentinel1SarSensorModel::readAnnotationFile(const std::string & annotationXml)
{
ossimRefPtr<ossimXmlDocument> xmlDoc = new ossimXmlDocument(annotationXml);
const ossimXmlNode & xmlRoot = *xmlDoc->getRoot();
//Parse specific metadata for Sentinel1
//TODO add as members to the Sentinel1SarSensorModel
const std::string & product_type = getTextFromFirstNode(xmlRoot, "adsHeader/productType");
const std::string & mode = getTextFromFirstNode(xmlRoot, "adsHeader/mode");
const std::string & swath = getTextFromFirstNode(xmlRoot, "adsHeader/swath");
const std::string & polarisation = getTextFromFirstNode(xmlRoot, "adsHeader/polarisation");
isGRD = (product_type == "GRD");
// First, lookup position/velocity records
std::vector<ossimRefPtr<ossimXmlNode> > xnodes;
xmlDoc->findNodes("/product/generalAnnotation/orbitList/orbit",xnodes);
//TODO uncomment and adapt following code from s1_inverse to fill
//SarSensorModel structure
for(std::vector<ossimRefPtr<ossimXmlNode> >::iterator itNode = xnodes.begin(); itNode!=xnodes.end();++itNode)
{
OrbitRecordType orbitRecord;
// Retrieve acquisition time
orbitRecord.azimuthTime = getTimeFromFirstNode(**itNode, attTime);
// Retrieve ECEF position
ossimXmlNode const& positionNode = getExpectedFirstNode(**itNode, attPosition);
orbitRecord.position[0] = getDoubleFromFirstNode(positionNode, attX);
orbitRecord.position[1] = getDoubleFromFirstNode(positionNode, attY);
orbitRecord.position[2] = getDoubleFromFirstNode(positionNode, attZ);
// Retrieve ECEF velocity
ossimXmlNode const& velocityNode = getExpectedFirstNode(**itNode, attVelocity);
orbitRecord.velocity[0] = getDoubleFromFirstNode(velocityNode, attX);
orbitRecord.velocity[1] = getDoubleFromFirstNode(velocityNode, attY);
orbitRecord.velocity[2] = getDoubleFromFirstNode(velocityNode, attZ);
//Add one orbits record
theOrbitRecords.push_back(orbitRecord);
}
//Parse the near range time (in seconds)
theNearRangeTime = getDoubleFromFirstNode(xmlRoot, "imageAnnotation/imageInformation/slantRangeTime");;
//Parse the range sampling rate
theRangeSamplingRate = getDoubleFromFirstNode(xmlRoot, "generalAnnotation/productInformation/rangeSamplingRate");;
//Parse the range resolution
theRangeResolution = getDoubleFromFirstNode(xmlRoot, "imageAnnotation/imageInformation/rangePixelSpacing");;
//Parse the radar frequency
theRadarFrequency = getDoubleFromFirstNode(xmlRoot, "generalAnnotation/productInformation/radarFrequency");;
//Parse azimuth time interval
theAzimuthTimeInterval = getDoubleFromFirstNode(xmlRoot, "imageAnnotation/imageInformation/azimuthTimeInterval")*1000000;
// Now read burst records as well
xnodes.clear();
xmlDoc->findNodes("/product/swathTiming/burstList/burst",xnodes);
if(xnodes.empty())
{
BurstRecordType burstRecord;
burstRecord.startLine = 0;
burstRecord.azimuthStartTime = getTimeFromFirstNode(xmlRoot,"imageAnnotation/imageInformation/productFirstLineUtcTime");
std::cout<< burstRecord.azimuthStartTime<<'\n';
burstRecord.azimuthStopTime = getTimeFromFirstNode(xmlRoot,"imageAnnotation/imageInformation/productLastLineUtcTime");
burstRecord.endLine = getTextFromFirstNode(xmlRoot, "imageAnnotation/imageInformation/numberOfLines").toUInt16()-1;
theBurstRecords.push_back(burstRecord);
}
else
{
unsigned int linesPerBurst = xmlRoot.findFirstNode("swathTiming/linesPerBurst")->getText().toUInt16();
unsigned int burstId(0);
for(std::vector<ossimRefPtr<ossimXmlNode> >::iterator itNode = xnodes.begin(); itNode!=xnodes.end();++itNode,++burstId)
{
BurstRecordType burstRecord;
const ossimSarSensorModel::TimeType azTime = getTimeFromFirstNode(**itNode, attAzimuthTime);
ossimString const& s = getTextFromFirstNode(**itNode, attFirstValidSample);
long first_valid(0), last_valid(0);
bool begin_found(false), end_found(false);
std::vector<ossimString> ssp = s.split(" ");
for (std::vector<ossimString>::const_iterator sIt = ssp.begin(), e = ssp.end()
; sIt != e && !end_found
; ++sIt
)
{
if(!begin_found)
{
if(*sIt!="-1")
{
begin_found = true;
}
else
{
++first_valid;
}
++last_valid;
}
else
{
if(!end_found && *sIt=="-1")
{
end_found = true;
}
else
{
++last_valid;
}
}
}
burstRecord.startLine = burstId*linesPerBurst + first_valid;
burstRecord.endLine = burstId*linesPerBurst + last_valid;
burstRecord.azimuthStartTime = azTime + boost::posix_time::microseconds(first_valid*theAzimuthTimeInterval);
burstRecord.azimuthStopTime = azTime + boost::posix_time::microseconds(last_valid*theAzimuthTimeInterval);
theBurstRecords.push_back(burstRecord);
}
}
if(isGRD)
{
readCoordinates(*xmlDoc,
"/product/coordinateConversion/coordinateConversionList/coordinateConversion",
attSr0, attSrgrCoefficients,
theSlantRangeToGroundRangeRecords);
readCoordinates(*xmlDoc,
"/product/coordinateConversion/coordinateConversionList/coordinateConversion",
attGr0, attGrsrCoefficients,
theGroundRangeToSlantRangeRecords);
}
xnodes.clear();
xmlDoc->findNodes("/product/geolocationGrid/geolocationGridPointList/geolocationGridPoint",xnodes);
for(std::vector<ossimRefPtr<ossimXmlNode> >::iterator itNode = xnodes.begin(); itNode!=xnodes.end();++itNode)
{
GCPRecordType gcpRecord;
// Retrieve acquisition time
gcpRecord.azimuthTime = getTimeFromFirstNode(**itNode, attAzimuthTime);
gcpRecord.slantRangeTime = getDoubleFromFirstNode(**itNode, attSlantRangeTime);
gcpRecord.imPt.x = getDoubleFromFirstNode(**itNode, attPixel);
// In TOPSAR products, GCPs are weird (they fall in black lines
// between burst. This code allows to move them to a valid area of
// the image.
if(theBurstRecords.size()>2)
{
ossimSarSensorModel::TimeType acqStart;
bool burstFound(false);
unsigned long acqStartLine(0);
for(std::vector<BurstRecordType>::reverse_iterator bIt = theBurstRecords.rbegin();bIt!=theBurstRecords.rend() && !burstFound;++bIt)
{
if(gcpRecord.azimuthTime >= bIt->azimuthStartTime && gcpRecord.azimuthTime < bIt->azimuthStopTime)
{
burstFound = true;
acqStart = bIt->azimuthStartTime;
acqStartLine = bIt->startLine;
}
}
if(!burstFound)
{
if(gcpRecord.azimuthTime < theBurstRecords.front().azimuthStartTime)
{
acqStart = theBurstRecords.front().azimuthStartTime;
acqStartLine = theBurstRecords.front().startLine;
}
else if (gcpRecord.azimuthTime >= theBurstRecords.front().azimuthStopTime)
{
acqStart = theBurstRecords.back().azimuthStartTime;
acqStartLine = theBurstRecords.back().startLine;
}
}
boost::posix_time::time_duration timeSinceStart = gcpRecord.azimuthTime - acqStart;
const double timeSinceStartInMicroSeconds = timeSinceStart.total_microseconds();
gcpRecord.imPt.y= timeSinceStartInMicroSeconds/theAzimuthTimeInterval + acqStartLine;
}
else
{
gcpRecord.imPt.y = getDoubleFromFirstNode(**itNode, attLine);;
}
ossimGpt geoPoint;
gcpRecord.worldPt.lat = getDoubleFromFirstNode(**itNode, attLatitude);
gcpRecord.worldPt.lon = getDoubleFromFirstNode(**itNode, attLongitude);
gcpRecord.worldPt.hgt = getDoubleFromFirstNode(**itNode, attHeight);
theGCPRecords.push_back(gcpRecord);
}
this->optimizeTimeOffsetsFromGcps();
}
} // namespace ossimplugins
//----------------------------------------------------------------------------
//
// "Copyright Centre National d'Etudes Spatiales"
//
// License: LGPL-2
//
// See LICENSE.txt file in the top level directory for more details.
//
//----------------------------------------------------------------------------
// $Id$
#ifndef ossimSentinel1SarSensorModel_HEADER
#define ossimSentinel1SarSensorModel_HEADER