Commit ab26cbd3 authored by Julien Osman's avatar Julien Osman
Browse files

Merge branch '2160-SarSensorModel' into 2159

parents 5ec85f33 23517628
Pipeline #7519 failed with stages
in 96 minutes and 19 seconds
......@@ -106,8 +106,8 @@ struct OTBMetadata_EXPORT Orbit
struct OTBMetadata_EXPORT BurstRecord
{
MetaData::TimeType azimuthStartTime;
unsigned long startLine;
MetaData::TimeType azimuthStopTime;
unsigned long startLine;
unsigned long endLine;
unsigned long startSample;
unsigned long endSample;
......@@ -127,6 +127,17 @@ struct OTBMetadata_EXPORT GCPTime
double slantRangeTime;
};
/** \struct CoordinateConversionRecord
*
* \brief This structure contains coefficients to convert between coordinates types, e.g.
* from ground range to slant range
*/
struct CoordinateConversionRecord
{
MetaData::TimeType azimuthTime;
double rg0;
std::vector<double> coeffs;
};
/** \struct SARParam
......@@ -145,6 +156,7 @@ struct OTBMetadata_EXPORT SARParam
MetaData::DurationType azimuthTimeInterval;
double nearRangeTime;
double rangeSamplingRate;
double rangeResolution;
/** Doppler centroid estimates */
std::vector<DopplerCentroid> dopplerCentroids;
......@@ -160,6 +172,12 @@ struct OTBMetadata_EXPORT SARParam
/** map between GCP ids and corresponding azimuth and range times */
std::unordered_map<std::string, GCPTime> gcpTimes;
/** Conversion coefficients from slant range to ground range */
std::vector<CoordinateConversionRecord> slantRangeToGroundRangeRecords;
/** Conversion coefficients from ground range to slant range */
std::vector<CoordinateConversionRecord> groundRangeToSlantRangeRecords;
};
/** \struct SARCalib
......
......@@ -126,7 +126,12 @@ protected:
std::vector<SARNoise> GetNoiseVector(const XMLMetadataSupplier&) const;
/* Fetch the burst records */
std::vector<BurstRecord> GetBurstRecords(const XMLMetadataSupplier&) const;
std::vector<BurstRecord> GetBurstRecords(const XMLMetadataSupplier&, const MetaData::DurationType & azimuthTimeInterval) const;
/* Fetch coordinate conversion records (Sr0/Gr0) */
std::vector<CoordinateConversionRecord> GetCoordinateConversionRecord(const XMLMetadataSupplier &xmlMS,
const std::string & rg0_path,
const std::string & coeffs_path) const;
/* Fetch azimuth and range times corresponding to GCPs */
std::unordered_map<std::string, GCPTime> GetGCPTimes(const XMLMetadataSupplier & xmlMS,
......
......@@ -440,8 +440,26 @@ std::vector<Orbit> Sentinel1ImageMetadataInterface::GetOrbits(const XMLMetadataS
std::vector<SARNoise> Sentinel1ImageMetadataInterface::GetNoiseVector(const XMLMetadataSupplier &xmlMS) const
{
std::vector<SARNoise> noiseVector;
std::string rangeNoisePrefix = "noise.noiseVectorList.";
std::string rangeVectorName = "noiseVector_";
std::string rangeLUTName = ".noiseLut";
// Number of entries in the vector
int listCount = xmlMS.GetAs<int>("noise.noiseVectorList.count");
int listCount = xmlMS.GetAs<int>(0, rangeNoisePrefix + "count");
// Starting with IPF 2.9.0, the noiseLut field has been renamed into noiseRangeVectorList, and the
// noiseAzimuthVectorList has been added.
// Ref: MPC-0392 DI-MPC-TN Issue 1.1 2017,Nov.28 Thermal Denoising of Products Generated by the S-1 IPF
if (!listCount)
{
rangeNoisePrefix = "noise.noiseRangeVectorList.";
rangeVectorName = "noiseRangeVector_";
rangeLUTName = ".noiseRangeLut";
listCount = xmlMS.GetAs<int>(rangeNoisePrefix + "count");
}
// This streams wild hold the iteration number
std::ostringstream oss;
for (int listId = 1 ; listId <= listCount ; ++listId)
......@@ -449,7 +467,7 @@ std::vector<SARNoise> Sentinel1ImageMetadataInterface::GetNoiseVector(const XMLM
oss.str("");
oss << listId;
// Base path to the data, that depends on the iteration number
std::string path_root = "noise.noiseVectorList.noiseVector_" + oss.str();
std::string path_root = rangeNoisePrefix + rangeVectorName + oss.str();
SARNoise noiseVect;
std::istringstream(xmlMS.GetAs<std::string>(path_root + ".azimuthTime")) >> noiseVect.azimuthTime;
MetaData::LUT1D noiseLut;
......@@ -457,26 +475,27 @@ std::vector<SARNoise> Sentinel1ImageMetadataInterface::GetNoiseVector(const XMLM
ax1.Size = xmlMS.GetAs<int>(path_root + ".pixel.count");
ax1.Values = xmlMS.GetAsVector<double>(path_root + ".pixel", ' ', ax1.Size);
noiseLut.Axis[0] = ax1;
noiseLut.Array = xmlMS.GetAsVector<double>(path_root + ".noiseLut",
' ', xmlMS.GetAs<int>(path_root + ".noiseLut.count"));
noiseLut.Array = xmlMS.GetAsVector<double>(path_root + rangeLUTName,
' ', xmlMS.GetAs<int>(path_root + rangeLUTName + ".count"));
noiseVect.noiseLut = noiseLut;
noiseVector.push_back(noiseVect);
}
return noiseVector;
}
std::vector<BurstRecord> Sentinel1ImageMetadataInterface::GetBurstRecords(const XMLMetadataSupplier &xmlMS) const
std::vector<BurstRecord> Sentinel1ImageMetadataInterface::GetBurstRecords(const XMLMetadataSupplier &xmlMS, const MetaData::DurationType& azimuthTimeInterval) const
{
std::vector<BurstRecord> burstRecords;
const std::string prefix = "product.swathTiming.";
// number of burst records
int listCount = xmlMS.GetAs<int>(0, "product.swathTiming.burstList.count");
int numberOfBursts = xmlMS.GetAs<int>(0, prefix + "burstList.count");
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));
if (!listCount)
if (!numberOfBursts)
{
BurstRecord record;
......@@ -497,14 +516,101 @@ std::vector<BurstRecord> Sentinel1ImageMetadataInterface::GetBurstRecords(const
}
else
{
// TODO
otbGenericExceptionMacro(MissingMetadataException,
<<"Burst record parsing not implemented yet");
const auto linesPerBurst = xmlMS.GetAs<int>(prefix + "linesPerBurst");
const auto samplesPerBurst = xmlMS.GetAs<int>(prefix + "samplesPerBurst");
for (int i = 0; i<numberOfBursts; i++)
{
const std::string burstPath = prefix + "burstList.burst_" + std::to_string(i+1) + ".";
BurstRecord record;
MetaData::TimeType azimuthTime;
ss << xmlMS.GetAs<std::string>(burstPath + "azimuthTime");
ss >> azimuthTime;
auto firstValidSamples = xmlMS.GetAsVector<int>(burstPath + "firstValidSample");
int firstValidSample = 0, lastValidSample = samplesPerBurst - 1;
int firstValid = 0, lastValid = firstValidSamples.size();
int currentIndex = 0;
bool firstIndexFound = false;
for (const auto & elem : firstValidSamples)
{
if (elem != -1)
{
if (!firstIndexFound)
{
firstIndexFound = true;
firstValid = currentIndex;
}
lastValid = currentIndex;
if (elem < samplesPerBurst && elem > firstValidSample)
{
firstValidSample = elem;
}
}
currentIndex++;
}
for (const auto & elem : xmlMS.GetAsVector<int>(burstPath + "lastValidSample"))
{
if (elem != -1 && elem < lastValidSample)
{
lastValidSample = elem;
}
}
record.startLine = i*linesPerBurst + firstValid;
record.endLine = i*linesPerBurst + lastValid;
record.startSample = firstValidSample;
record.endSample = lastValidSample;
record.azimuthAnxTime = xmlMS.GetAs<double>(burstPath + "azimuthAnxTime");
record.azimuthStartTime = azimuthTime + firstValid * azimuthTimeInterval;
record.azimuthStopTime = azimuthTime + lastValid * azimuthTimeInterval;
burstRecords.push_back(record);
}
}
return burstRecords;
}
std::vector<CoordinateConversionRecord> Sentinel1ImageMetadataInterface::GetCoordinateConversionRecord(const XMLMetadataSupplier &xmlMS,
const std::string & rg0_path,
const std::string & coeffs_path) const
{
std::vector<CoordinateConversionRecord> records;
const std::string prefixPath = "product.coordinateConversion.coordinateConversionList.";
auto listCount = xmlMS.GetAs<int>(prefixPath + "count");
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 i = 0; i < listCount; i++)
{
CoordinateConversionRecord coordinateConversionRecord;
const std::string recordPath = prefixPath +"coordinateConversion_" + std::to_string(i+1) + ".";
//coordinateConversionRecord.azimuthTime;
ss << xmlMS.GetAs<std::string>(recordPath + "azimuthTime");
ss >> coordinateConversionRecord.azimuthTime;
coordinateConversionRecord.rg0 = xmlMS.GetAs<double>(recordPath + rg0_path);
coordinateConversionRecord.coeffs = xmlMS.GetAsVector<double>(recordPath + coeffs_path, ' ');
records.push_back(coordinateConversionRecord);
}
return records;
}
std::unordered_map<std::string, GCPTime> Sentinel1ImageMetadataInterface::GetGCPTimes(const XMLMetadataSupplier &xmlMS,
const Projection::GCPParam & gcps) const
{
......@@ -595,12 +701,18 @@ void Sentinel1ImageMetadataInterface::ParseGdal(ImageMetadata & imd)
sarParam.dopplerCentroids = this->GetDopplerCentroid(AnnotationMS);
sarParam.orbits = this->GetOrbits(AnnotationMS);
sarParam.burstRecords = this->GetBurstRecords(AnnotationMS);
sarParam.slantRangeToGroundRangeRecords = this->GetCoordinateConversionRecord(AnnotationMS, "sr0", "srgrCoefficients");
sarParam.groundRangeToSlantRangeRecords = this->GetCoordinateConversionRecord(AnnotationMS, "gr0", "grsrCoefficients");
sarParam.azimuthTimeInterval = boost::posix_time::precise_duration(
AnnotationMS.GetAs<double>("product.imageAnnotation.imageInformation.azimuthTimeInterval")
* 1e6);
sarParam.burstRecords = this->GetBurstRecords(AnnotationMS, sarParam.azimuthTimeInterval);
sarParam.nearRangeTime = AnnotationMS.GetAs<double>("product.imageAnnotation.imageInformation.slantRangeTime");
sarParam.rangeSamplingRate = AnnotationMS.GetAs<double>("product.generalAnnotation.productInformation.rangeSamplingRate");
sarParam.rangeResolution = AnnotationMS.GetAs<double>("product.imageAnnotation.imageInformation.rangePixelSpacing");
sarParam.gcpTimes = this->GetGCPTimes(AnnotationMS, imd.GetGCPParam());
......
......@@ -958,6 +958,95 @@ void TerraSarImageMetadataInterface::PrintSelf(std::ostream& os, itk::Indent ind
Superclass::PrintSelf(os, indent);
}
void ReadSARSensorModel(const XMLMetadataSupplier & xmlMS,
const XMLMetadataSupplier & geoXmlMS,
const Projection::GCPParam & gcps,
SARParam & param)
{
// Number of entries in the vector
int listCount = xmlMS.GetAs<int>("level1Product.platform.orbit.orbitHeader.numStateVectors");
const auto readFormattedDate = [](const std::string & dateStr, const std::string & format)
{
MetaData::TimeType outputDate;
std::stringstream ss;
auto facet = new boost::posix_time::time_input_facet(format);
ss.imbue(std::locale(std::locale(), facet));
ss << dateStr;
ss >> outputDate;
return outputDate;
};
const std::string dateFormat = "%Y-%m-%dT%H:%M:%S%F";
// This streams wild hold the iteration number
std::ostringstream oss;
for (int listId = 1 ; listId <= listCount ; ++listId)
{
oss.str("");
oss << listId;
// Base path to the data, that depends on the iteration number
std::string path_root = "level1Product.platform.orbit.stateVec_" + oss.str();
Orbit orbit;
orbit.time = readFormattedDate(xmlMS.GetAs<std::string>(path_root + ".timeUTC"), dateFormat);
orbit.position[0] = xmlMS.GetAs<double>(path_root + ".posX");
orbit.position[1] = xmlMS.GetAs<double>(path_root + ".posY");
orbit.position[2] = xmlMS.GetAs<double>(path_root + ".posZ");
orbit.velocity[0] = xmlMS.GetAs<double>(path_root + ".velX");
orbit.velocity[1] = xmlMS.GetAs<double>(path_root + ".velY");
orbit.velocity[2] = xmlMS.GetAs<double>(path_root + ".velZ");
param.orbits.push_back(orbit);
}
param.nearRangeTime = xmlMS.GetAs<double>("level1Product.productInfo.sceneInfo.rangeTime.firstPixel");
param.rangeSamplingRate = xmlMS.GetAs<double>("level1Product.instrument.settings.RSF");
const auto azimuthTimeStart = readFormattedDate(
xmlMS.GetAs<std::string>("level1Product.productInfo.sceneInfo.start.timeUTC"), dateFormat);
const auto azimuthTimeStop = readFormattedDate(
xmlMS.GetAs<std::string>("level1Product.productInfo.sceneInfo.stop.timeUTC"), dateFormat);
const auto td = azimuthTimeStop - azimuthTimeStart;
assert(td > boost::posix_time::seconds(0));
const auto numberOfRows = xmlMS.GetAs<double>("level1Product.productInfo.imageDataInfo.imageRaster.numberOfRows");
param.azimuthTimeInterval = td / numberOfRows;
//For Terrasar-X only 1 burst is supported for now
BurstRecord burstRecord;
burstRecord.azimuthStartTime = azimuthTimeStart;
burstRecord.azimuthStopTime = azimuthTimeStop;
burstRecord.startLine = 0;
burstRecord.endLine = numberOfRows - 1;
burstRecord.startSample = 0;
burstRecord.endSample = xmlMS.GetAs<double>("level1Product.productInfo.imageDataInfo.imageRaster.numberOfColumns") - 1;
burstRecord.azimuthAnxTime = 0.;
param.burstRecords.push_back(burstRecord);
for (const auto & gcp : gcps.GCPs)
{
const std::string prefix = "geoReference.geolocationGrid.gridPoint_" + gcp.m_Id + ".";
GCPTime time;
auto deltaAz = boost::posix_time::precise_duration(geoXmlMS.GetAs<double>(prefix + "t"));
time.azimuthTime = azimuthTimeStart + deltaAz;
time.slantRangeTime = param.nearRangeTime + geoXmlMS.GetAs<double>(prefix + "tau");
param.gcpTimes[gcp.m_Id] = time;
}
}
void TerraSarImageMetadataInterface::ParseGdal(ImageMetadata &imd)
{
// Metadata read by GDAL
......@@ -967,22 +1056,27 @@ void TerraSarImageMetadataInterface::ParseGdal(ImageMetadata &imd)
// Main XML file
std::string MainFilePath = ExtractXMLFiles::GetResourceFile(itksys::SystemTools::GetFilenamePath(m_MetadataSupplierInterface->GetResourceFile("")), "T[S|D]X1_SAR__.*.xml") ;
if (!MainFilePath.empty())
if (MainFilePath.empty())
{
XMLMetadataSupplier MainXMLFileMetadataSupplier(MainFilePath);
imd.Add(MDNum::LineSpacing, MainXMLFileMetadataSupplier.GetAs<double>("level1Product.productSpecific.complexImageInfo.projectedSpacingAzimuth"));
imd.Add(MDNum::PixelSpacing, MainXMLFileMetadataSupplier.GetAs<double>("level1Product.productSpecific.complexImageInfo.projectedSpacingRange.slantRange"));
imd.Add(MDStr::Mission, MainXMLFileMetadataSupplier.GetAs<std::string>("level1Product.generalHeader.mission"));
imd.Add(MDStr::ProductType, MainXMLFileMetadataSupplier.GetAs<std::string>("level1Product.productInfo.productVariantInfo.productType"));
imd.Add(MDStr::SensorID, MainXMLFileMetadataSupplier.GetAs<std::string>("level1Product.productInfo.acquisitionInfo.sensor"));
imd.Add(MDNum::RadarFrequency, MainXMLFileMetadataSupplier.GetAs<double>("level1Product.instrument.radarParameters.centerFrequency"));
imd.Add(MDTime::AcquisitionStartTime, MainXMLFileMetadataSupplier.GetFirstAs<MetaData::Time>("level1Product.productInfo.sceneInfo.start.timeUTC"));
imd.Add(MDTime::AcquisitionStopTime, MainXMLFileMetadataSupplier.GetFirstAs<MetaData::Time>("level1Product.productInfo.sceneInfo.stop.timeUTC"));
imd.Add(MDNum::RangeTimeFirstPixel, MainXMLFileMetadataSupplier.GetFirstAs<double>("level1Product.productInfo.sceneInfo.rangeTime.firstPixel"));
imd.Add(MDNum::RangeTimeLastPixel, MainXMLFileMetadataSupplier.GetFirstAs<double>("level1Product.productInfo.sceneInfo.rangeTime.lastPixel"));
imd.Add(MDNum::PRF, MainXMLFileMetadataSupplier.GetAs<double>("level1Product.productSpecific.complexImageInfo.commonPRF"));
imd.Add(MDNum::RSF, MainXMLFileMetadataSupplier.GetAs<double>("level1Product.productSpecific.complexImageInfo.commonRSF"));
auto numberOfCalFactor = MainXMLFileMetadataSupplier.GetNumberOf("level1Product.calibration.calibrationConstant");
otbGenericExceptionMacro(MissingMetadataException,
<< "Cannot find the TerraSar main XML file");
}
XMLMetadataSupplier MainXMLFileMetadataSupplier(MainFilePath);
imd.Add(MDNum::LineSpacing, MainXMLFileMetadataSupplier.GetAs<double>("level1Product.productSpecific.complexImageInfo.projectedSpacingAzimuth"));
imd.Add(MDNum::PixelSpacing, MainXMLFileMetadataSupplier.GetAs<double>("level1Product.productSpecific.complexImageInfo.projectedSpacingRange.slantRange"));
imd.Add(MDStr::Mission, MainXMLFileMetadataSupplier.GetAs<std::string>("level1Product.generalHeader.mission"));
imd.Add(MDStr::ProductType, MainXMLFileMetadataSupplier.GetAs<std::string>("level1Product.productInfo.productVariantInfo.productType"));
// imd.Add(MDStr::Mode, MainXMLFileMS.GetAs<std::string>("level1Product.productInfo.acquisitionInfo.imagingMode"));
imd.Add(MDStr::SensorID, MainXMLFileMetadataSupplier.GetAs<std::string>("level1Product.productInfo.acquisitionInfo.sensor"));
imd.Add(MDNum::RadarFrequency, MainXMLFileMetadataSupplier.GetAs<double>("level1Product.instrument.radarParameters.centerFrequency"));
imd.Add(MDTime::AcquisitionStartTime, MainXMLFileMetadataSupplier.GetFirstAs<MetaData::Time>("level1Product.productInfo.sceneInfo.start.timeUTC"));
imd.Add(MDTime::AcquisitionStopTime, MainXMLFileMetadataSupplier.GetFirstAs<MetaData::Time>("level1Product.productInfo.sceneInfo.stop.timeUTC"));
imd.Add(MDNum::RangeTimeFirstPixel, MainXMLFileMetadataSupplier.GetFirstAs<double>("level1Product.productInfo.sceneInfo.rangeTime.firstPixel"));
imd.Add(MDNum::RangeTimeLastPixel, MainXMLFileMetadataSupplier.GetFirstAs<double>("level1Product.productInfo.sceneInfo.rangeTime.lastPixel"));
imd.Add(MDNum::PRF, MainXMLFileMetadataSupplier.GetAs<double>("level1Product.productSpecific.complexImageInfo.commonPRF"));
imd.Add(MDNum::RSF, MainXMLFileMetadataSupplier.GetAs<double>("level1Product.productSpecific.complexImageInfo.commonRSF"));
auto numberOfCalFactor = MainXMLFileMetadataSupplier.GetNumberOf("level1Product.calibration.calibrationConstant");
if(numberOfCalFactor == 1)
{
imd.Add(MDNum::CalFactor, MainXMLFileMetadataSupplier.GetAs<double>("level1Product.calibration.calibrationConstant.calFactor"));
......@@ -999,19 +1093,26 @@ void TerraSarImageMetadataInterface::ParseGdal(ImageMetadata &imd)
imd.Add(MDNum::CalScale, MainXMLFileMetadataSupplier.GetAs<double>(path.str()));
}
}
SARCalib sarCalib;
LoadRadiometricCalibrationData(sarCalib, MainXMLFileMetadataSupplier, imd);
imd.Add(MDGeom::SARCalib, sarCalib);
}
SARCalib sarCalib;
LoadRadiometricCalibrationData(sarCalib, MainXMLFileMetadataSupplier, imd);
imd.Add(MDGeom::SARCalib, sarCalib);
assert(m_MetadataSupplierInterface->GetNbBands() == imd.Bands.size());
SARParam sarParam;
// Open the georef file containing GCPs
XMLMetadataSupplier GCPXMLFileMS(m_MetadataSupplierInterface->GetResourceFile()
+ "/ANNOTATION/GEOREF.xml");
ReadSARSensorModel(MainXMLFileMetadataSupplier, GCPXMLFileMS, imd.GetGCPParam(), sarParam);
for (int bandId = 0 ; bandId < m_MetadataSupplierInterface->GetNbBands() ; ++bandId)
{
Fetch(MDStr::Polarization, imd, "POLARIMETRIC_INTERP", bandId);
imd.Bands[bandId].Add(MDGeom::SAR, sarParam);
}
}
void TerraSarImageMetadataInterface::ParseGeom(ImageMetadata & imd)
......@@ -1065,7 +1166,7 @@ void TerraSarImageMetadataInterface::Parse(ImageMetadata & imd)
// Failed to fetch the metadata
else
otbGenericExceptionMacro(MissingMetadataException,
<< "Not a Sentinel1 product");
<< "Not a TerraSar product");
}
......
......@@ -88,6 +88,16 @@ private:
void AzimuthTimeToLine(const TimeType & azimuthTime,
double & line) const;
void SlantRangeToGroundRange(const double & slantRange,
const TimeType & azimuthTime,
double & groundRange) const;
void ApplyCoordinateConversion(const double & in,
const TimeType& azimuthTime,
const std::vector<CoordinateConversionRecord> & records,
double & out) const;
const ImageMetadata & m_Imd;
SARParam m_SarParam;
......
......@@ -93,7 +93,13 @@ void SarSensorModel::WorldToLineSample(const Point3DType& inGeoPoint, Point2DTyp
if (m_IsGrd)
{
//TODO GRD case
// GRD case
double groundRange(0);
SlantRangeToGroundRange(rangeTime*C/2,azimuthTime,groundRange);
// Eq 32 p. 31
outLineSample[0] = groundRange/m_SarParam.rangeResolution;
//imPt.x = groundRange/theRangeResolution;
}
else
{
......@@ -383,4 +389,92 @@ 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
{
ApplyCoordinateConversion(slantRange, azimuthTime, m_SarParam.slantRangeToGroundRangeRecords , groundRange);
}
void SarSensorModel::ApplyCoordinateConversion(const double & in,
const TimeType& azimuthTime,
const std::vector<CoordinateConversionRecord> & records,
double & out) const
{
assert(!records.empty()&&"The records vector is empty.");
// First, we need to find the correct pair of records for interpolation
auto it = records.cbegin();
auto previousRecord = it;
++it;
auto nextRecord = it;
// Look for the correct record
while(it!=records.end())
{
if(azimuthTime >= previousRecord->azimuthTime
&& azimuthTime < nextRecord->azimuthTime)
{
break;
}
previousRecord = nextRecord;
++it;
nextRecord = it;
}
CoordinateConversionRecord srgrRecord;
if(it == records.cend())
{
if(azimuthTime < records.front().azimuthTime)
{
srgrRecord = records.front();
}
else if(azimuthTime >= records.back().azimuthTime)
{
srgrRecord = records.back();
}
}
else
{
assert(nextRecord != records.end());
assert(!previousRecord->coeffs.empty()&&"previousRecord coefficients vector is empty.");
assert(!nextRecord->coeffs.empty()&&"nextRecord coefficients vector is empty.");
// If azimuth time is between 2 records, interpolate
const double interp
= DurationType(azimuthTime - previousRecord->azimuthTime)
/ (nextRecord->azimuthTime - previousRecord->azimuthTime)
;
srgrRecord.rg0 = (1-interp) * previousRecord->rg0 + interp*nextRecord->rg0;
srgrRecord.coeffs.clear();
std::vector<double>::const_iterator pIt = previousRecord->coeffs.cbegin();
std::vector<double>::const_iterator nIt = nextRecord->coeffs.cbegin();
for(;pIt != previousRecord->coeffs.end() && nIt != nextRecord->coeffs.cend();++pIt,++nIt)
{
srgrRecord.coeffs.push_back(interp*(*nIt)+(1-interp)*(*pIt));
}
assert(!srgrRecord.coeffs.empty()&&"Slant range to ground range interpolated coefficients vector is empty.");
}
// Now that we have the interpolated coeffs, compute ground range
// from slant range
const double sr_minus_sr0 = in-srgrRecord.rg0;
assert(!srgrRecord.coeffs.empty()&&"Slant range to ground range coefficients vector is empty.");
out = 0.;
for(auto cIt = srgrRecord.coeffs.crbegin();cIt!=srgrRecord.coeffs.crend();++cIt)
{
out = *cIt + sr_minus_sr0*out;
}
}
} //namespace otb
\ No newline at end of file
......@@ -241,4 +241,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)