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

ENH: add support for GRD and multi-burst products

parent d4f3ba79
Pipeline #7471 canceled with stages
in 4 minutes and 52 seconds
......@@ -126,8 +126,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;
......@@ -147,6 +147,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
......@@ -171,6 +182,7 @@ struct OTBMetadata_EXPORT SARParam
MetaData::DurationType azimuthTimeInterval;
double nearRangeTime;
double rangeSamplingRate;
double rangeResolution;
/** Doppler centroid estimates */
std::vector<DopplerCentroid> dopplerCentroids;
......@@ -186,6 +198,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
......
......@@ -129,7 +129,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,
......
......@@ -540,18 +540,19 @@ std::vector<SARNoise> Sentinel1ImageMetadataInterface::GetNoiseVector(const XMLM
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;
......@@ -572,14 +573,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
{
......@@ -670,12 +758,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());
......
......@@ -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
{
  • All these loops can be replaced with something more concise (and I guess the other as well). In SRGR case I already did the work in an unpublished branch

    Note that we don't need to copy the coef vector ( srgrRecord = records.back(); x 2), but as I left this optimization in comments I guess the difference wasn't perceptible -- I my benchs most of the time was lost in the use of MAT library (too many allocations)

       void ossimSarSensorModel::applyCoordinateConversion(const double & in, const TimeType& azimuthTime, const std::vector<CoordinateConversionRecordType> & records, double & out) const
       {
          assert(!records.empty()&&"The records vector is empty.");
          // std::clog << "conv coord(" << in << ", az="<<azimuthTime<<")\n";
    
          // First, we need to find the correct pair of records for interpolation
          auto it = records.begin();
    
          CoordinateConversionRecordType srgrRecord;
          auto nextRecord = std::find_if(std::begin(records), std::end(records),
                [azimuthTime](auto const& record) { return azimuthTime < record.azimuthTime; }
                );
          // Note: I did try with upper_bound instead of find_if (as the elements are sorted, but there aren't enough elements to make a perceptible difference
    
    
          assert(nextRecord == records.end() || nextRecord->azimuthTime > azimuthTime);
          if (nextRecord == records.end())
          {
             srgrRecord = records.back(); // ### don't copy the coefs vector
             // auto const& ref = records.back();
             // srgrRecord.azimuthTime = ref.azimuthTime;
             // srgrRecord.rg0         = ref.rg0        ;
          }
          else if (nextRecord == records.begin())
          {
             assert(azimuthTime < nextRecord->azimuthTime);
             srgrRecord = records.front(); // ### don't copy the coefs vector
             // auto const& ref = records.front();
             // srgrRecord.azimuthTime = ref.azimuthTime;
             // srgrRecord.rg0         = ref.rg0        ;
          }
          else
          {
             auto previousRecord = nextRecord - 1;
             assert(previousRecord->azimuthTime <= azimuthTime);
             assert(azimuthTime < nextRecord->azimuthTime);
             assert(nextRecord != records.end());
             assert(!previousRecord->coefs.empty()&&"previousRecord coefficients vector is empty.");
             assert(!nextRecord->coefs.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)
                ;
             // std::clog << "interp: " << interp << " ="
             // << " (" << azimuthTime             << " - " << previousRecord->azimuthTime << " (="<< (azimuthTime             - previousRecord->azimuthTime)<< ") )"
             // << "/(" << nextRecord->azimuthTime << " - " << previousRecord->azimuthTime << " (="<< (nextRecord->azimuthTime - previousRecord->azimuthTime)<< ") )"
             // << "\n";
    
             srgrRecord.rg0 = (1.0-interp) * previousRecord->rg0 + interp*nextRecord->rg0;
    
             srgrRecord.coefs.clear();
             auto pIt = previousRecord->coefs.begin();
             auto nIt = nextRecord->coefs.begin();
    
             srgrRecord.coefs.reserve(std::min(previousRecord->coefs.size(), nextRecord->coefs.size()));
             for(;pIt != previousRecord->coefs.end() && nIt != nextRecord->coefs.end();++pIt,++nIt)
             {
                srgrRecord.coefs.push_back(interp*(*nIt)+(1.0-interp)*(*pIt));
             }
    
             assert(!srgrRecord.coefs.empty()&&"Slant range to ground range interpolated coefficients vector is empty.");
          }
    
          // Now that we have the interpolated coefs, compute ground range
          // from slant range
          const double sr_minus_sr0 =  in-srgrRecord.rg0;
    
          assert(!srgrRecord.coefs.empty()&&"Slant range to ground range coefficients vector is empty.");
    
          double rout = 0.0; // avoid aliasing in loop
    
          for(auto cIt = srgrRecord.coefs.rbegin();cIt!=srgrRecord.coefs.rend();++cIt)
          {
             rout = *cIt + sr_minus_sr0*rout;
          }
          out = rout;
       }
  • I my benchs most of the time was lost in the use of MAT library (too many allocations)

    I plan on using itk fixed size matrices (itk::Matrix) and vectors (itk::Vector) to replace the NEWMAT objets in the inverse SAR transform. This should reduce the number of allocations.

  • If I'm not mistaken, ITK have started using eigen. It'll be nice to rely on this library for all the transformations.

  • Eigen3 has been integrated into ITK release 5.0.0. Quoting the ITK release notes:

    Eigen, a C++ library for linear algebra, was integrated in ITK and gradual transition from VXL to Eigen for Linear algebra is planned

    It looks like a long term goal of ITK is to replace uses of VXL by Eigen (e.g. use vectors from Eigen in itk::VariableLengthVector). See for example https://discourse.itk.org/t/add-internal-third-party-module-eigen3/1429

    Not that in current release (5.2.1), vnl objects are still used as internal type of ITK vector and matrix objects.

Please register or sign in to reply
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
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