ossimSarSensorModel.cpp 75.6 KB
Newer Older
1
/*
Julien Michel's avatar
Julien Michel committed
2
 * Copyright (C) 2005-2019 by Centre National d'Etudes Spatiales (CNES)
3 4 5 6 7 8 9 10 11 12 13 14
 *
 * This file is licensed under MIT license:
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
15
 
16 17 18 19 20 21 22 23 24
 * 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
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

25

26 27 28 29 30
#include "ossim/ossimSarSensorModel.h"
#include "ossim/ossimKeyWordListUtilities.h"
#include "ossim/ossimTraceHelpers.h"
#include "ossim/ossimRangeUtilities.h"
#include "ossim/ossimSarSensorModelPathsAndKeys.h"
31
#include <ossim/base/ossimRegExp.h>
32
#include <ossim/base/ossimLsrSpace.h>
33
#include <ossim/base/ossimKeywordNames.h>
34 35 36 37 38 39 40
#include <boost/static_assert.hpp>
#include <iostream>
#include <vector>
#include <algorithm>

#if defined(USE_BOOST_TIME)
using boost::posix_time::microseconds;
41
using boost::posix_time::seconds;
42 43
#else
using ossimplugins::time::microseconds;
44
using ossimplugins::time::seconds;
45
#endif
46 47

namespace {// Anonymous namespace
48
   const bool         k_verbose = false; // global verbose constant; TODO: use an option
49 50 51 52 53 54 55 56

   // Sometimes, we don't need to compare the actual distance, its square value is
   // more than enough.
   inline double squareDistance(ossimDpt const& lhs, ossimDpt const& rhs) {
      const double dx = lhs.x - rhs.x;
      const double dy = lhs.y - rhs.y;
      return dx*dx + dy*dy;
   }
57

58 59 60 61 62
   inline double squareDistance(ossimGpt const& lhs, ossimGpt const& rhs) {
      const ossimEcefPoint l(lhs);
      const ossimEcefPoint r(rhs);
      return (l-r).norm2();
   }
63

64
   template <typename Container>
65
      inline void unzip(Container const& in, Container& out_even, Container& out_odd)
66
      {
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
         typedef typename Container::const_iterator const_iterator;
         typedef typename Container::size_type      size_type;

         const size_type size                         = in.size();
         const bool      has_a_odd_number_of_elements = size % 2 == 1;

         out_even.reserve(size/2+1);
         out_odd.reserve(size/2);

         const_iterator end = in.end();
         if (has_a_odd_number_of_elements)
         {
            std::advance(end, -1);
         }
         for (const_iterator it=in.begin(); it != end ; )
         {
            out_even.push_back(*it++);
            out_odd.push_back(*it++);
         }
         if (has_a_odd_number_of_elements)
         {
            assert(end != in.end());
            out_even.push_back(*end);
         }
         assert(out_even.size() >= out_odd.size());
         assert(out_odd.capacity() == size/2); // The correct number of element have been reserved
         assert(out_even.capacity() == size/2+1); // The correct number of element have been reserved
         assert(out_odd.size() + out_even.size() == size);
95 96 97 98
      }

   ossimTrace traceExec  ("ossimSarSensorModel:exec");
   ossimTrace traceDebug ("ossimSarSensorModel:debug");
99

100 101
   typedef char const* const* strings_iterator;
   static char const* const PRODUCTTYPE_STRINGS[] = { "SLC", "GRD", "MGD", "GEC", "EEC" };
102 103 104 105
}// Anonymous namespace

namespace ossimplugins
{
106 107
   RTTI_DEF1(ossimSarSensorModel, "ossimSarSensorModel", ossimSensorModel);

108
   const double ossimSarSensorModel::C = 299792458;
109

110 111
   const unsigned int ossimSarSensorModel::thePluginVersion = 3;
  const unsigned int ossimSarSensorModel::thePluginVersionMin = 2;
112

113 114 115 116 117 118 119 120 121 122
   ossimSarSensorModel::ProductType::ProductType(string_view const& s)
   {
      using ossimplugins::begin;
      using ossimplugins::end;
      strings_iterator const ProductType_it = std::find(begin(::PRODUCTTYPE_STRINGS),end(::PRODUCTTYPE_STRINGS), s);
      if (ProductType_it == end(::PRODUCTTYPE_STRINGS))  {
         throw std::runtime_error("Invalid Sar Sensor Product type: `"+s+"'");
      }
      m_value = Type(std::distance(begin(::PRODUCTTYPE_STRINGS), ProductType_it));
      assert(m_value < MAX__);
123 124
   }

125 126 127 128 129 130 131
   string_view ossimSarSensorModel::ProductType::ToString() const
   {
      BOOST_STATIC_ASSERT((MAX__ == array_size(::PRODUCTTYPE_STRINGS)));
      assert(m_value != UNDEFINED__); // Yes, I know UNDEFINED__ > MAX__
      assert(m_value < MAX__);
      return PRODUCTTYPE_STRINGS[m_value];
   }
132

133 134
   ossimSarSensorModel::ossimSarSensorModel()
      : theRadarFrequency(0.),
135
      theAzimuthTimeInterval(seconds(0)),
136 137 138 139
      theNearRangeTime(0.),
      theRangeSamplingRate(0.),
      theRangeResolution(0.),
      theBistaticCorrectionNeeded(false),
140
      theAzimuthTimeOffset(seconds(0)),
141
      theRangeTimeOffset(0.),
142 143
      theRightLookingFlag(true),
      redaptMedataAfterDeburst(false)
144 145 146 147 148 149
      {}

   ossimSarSensorModel::GCPRecordType const&
      ossimSarSensorModel::findClosestGCP(ossimDpt const& imPt) const
      {
         assert(!theGCPRecords.empty()&&"theGCPRecords is empty.");
150

151 152
         // Find the closest GCP
         double distance2 = squareDistance(imPt, theGCPRecords.front().imPt);
153

154 155 156 157 158 159
         std::vector<GCPRecordType>::const_iterator refGcp = theGCPRecords.begin();
         std::vector<GCPRecordType>::const_iterator gcpIt  = theGCPRecords.begin();
         std::vector<GCPRecordType>::const_iterator gcpEnd = theGCPRecords.end();
         for(++gcpIt ; gcpIt!=gcpEnd ; ++gcpIt)
         {
            const double currentDistance2 = squareDistance(imPt, gcpIt->imPt);
160

161 162 163 164 165 166 167 168 169
            if(currentDistance2 < distance2)
            {
               distance2 = currentDistance2;
               refGcp = gcpIt;
            }
         }

         assert(refGcp != theGCPRecords.end() && "A GCP record shall have been found!");
         return *refGcp;
170 171
      }

172 173
   void ossimSarSensorModel::lineSampleHeightToWorld(const ossimDpt& imPt, const double & heightAboveEllipsoid, ossimGpt& worldPt) const
   {
174
      // std::clog << "ossimSarSensorModel::lineSampleHeightToWorld()\n";
175
      assert(!theGCPRecords.empty()&&"theGCPRecords is empty.");
176

177
      GCPRecordType const& refGcp = findClosestGCP(imPt);
178

179 180 181 182 183 184 185 186 187 188 189
      // Set the height reference
      ossim_float64 hgtSet;
      if ( ossim::isnan(heightAboveEllipsoid) )
      {
         hgtSet = refGcp.worldPt.height();
      }
      else
      {
         hgtSet = heightAboveEllipsoid;
      }
      const ossimHgtRef hgtRef(AT_HGT, hgtSet);
190

191
      ossimEcefPoint ellPt;
192

193 194
      // Simple iterative inversion of inverse model starting at closest gcp
      projToSurface(refGcp,imPt,hgtRef,ellPt);
195

196 197
      worldPt = ossimGpt(ellPt);
   }
198

199 200
   void ossimSarSensorModel::lineSampleToWorld(ossimDpt const& imPt, ossimGpt& worldPt) const
   {
201
      // std::clog << "ossimSarSensorModel::lineSampleToWorld()\n";
202
      assert(!theGCPRecords.empty()&&"theGCPRecords is empty.");
203

204 205
      GCPRecordType const& refGcp = findClosestGCP(imPt);
      ossimGpt      const& refPt = refGcp.worldPt;
206

207
      const ossimHgtRef hgtRef(AT_DEM);
208

209
      ossimEcefPoint ellPt;
210

211 212
      // Simple iterative inversion of inverse model starting at closest gcp
      projToSurface(refGcp,imPt,hgtRef,ellPt);
213

214 215
      worldPt = ossimGpt(ellPt);
   }
216

217 218
   void ossimSarSensorModel::worldToLineSample(const ossimGpt& worldPt, ossimDpt & imPt) const
   {
219
      // std::clog << "ossimSarSensorModel::worldToLineSample()\n";
220
      assert(theRangeResolution>0&&"theRangeResolution is null.");
221

222 223 224
      // First compute azimuth and range time
      TimeType azimuthTime;
      double rangeTime;
225

226 227 228 229
      ossimEcefPoint sensorPos;
      ossimEcefVector sensorVel;

      const bool success = worldToAzimuthRangeTime(worldPt, azimuthTime, rangeTime, sensorPos, sensorVel);
230

231 232 233 234 235
      if(!success)
      {
         imPt.makeNan();
         return;
      }
236 237
      // std::clog << "AzimuthTime: " << azimuthTime << "\n";
      // std::clog << "RangeTime: " << rangeTime << "\n";
238
      // std::clog << "GRD: " << isGRD() << "\n";
239

240 241
      // Convert azimuth time to line
      azimuthTimeToLine(azimuthTime,imPt.y);
242

243 244 245 246 247
      if(isGRD())
      {
         // GRD case
         double groundRange(0);
         slantRangeToGroundRange(rangeTime*C/2,azimuthTime,groundRange);
248 249
         // std::clog << "GroundRange: " << groundRange << "\n";
         // std::clog << "TheRangeResolution: " << theRangeResolution << "\n";
250 251 252 253 254 255 256 257

         // Eq 32 p. 31
         // TODO: possible micro-optimization: precompute 1/theRangeResolution, and
         // use *
         imPt.x = groundRange/theRangeResolution;
      }
      else
      {
258 259
         // std::clog << "TheNearRangeTime: " << theNearRangeTime << "\n";
         // std::clog << "TheRangeSamplingRate: " << theRangeSamplingRate << "\n";
260 261 262 263 264
         // SLC case
         // Eq 23 and 24 p. 28
         imPt.x = (rangeTime - theNearRangeTime)*theRangeSamplingRate;
      }
   }
265

266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319
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)
      {
         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;
      }

      // 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 ...
320
      assert(NormeS>1e-6);
321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340
      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, ossimEcefPoint & interpSensorPos, ossimEcefVector & interpSensorVel) const
341
   {
342
      // std::clog << "ossimSarSensorModel::worldToAzimuthRangeTime()\n";
343 344
      // First convert lat/lon to ECEF
      ossimEcefPoint inputPt(worldPt);
345

346 347
      // Compute zero doppler time
      TimeType interpTime;
348

349
      const bool success = zeroDopplerLookup(inputPt,azimuthTime,interpSensorPos,interpSensorVel);
350

351 352 353 354 355
      if(!success)
      {
         // TODO: check whether we could throw instead
         return false;
      }
356

357 358 359 360 361
      if(theBistaticCorrectionNeeded)
      {
         // Compute bistatic correction if needed
         DurationType bistaticCorrection;
         computeBistaticCorrection(inputPt,interpSensorPos,bistaticCorrection);
362

363 364
         // Update interpolated azimuth time
         azimuthTime += bistaticCorrection;
365

366 367 368
         // Update sensor position and velocity
         interpolateSensorPosVel(interpTime,interpSensorPos,interpSensorVel);
      }
369

370 371 372
      // rangeTime is the round-tripping time to target
      const double rangeDistance = (interpSensorPos-inputPt).magnitude();
      rangeTime = theRangeTimeOffset + 2*rangeDistance/C;
373

374 375
      return true;
   }
376

377 378
   void ossimSarSensorModel::lineSampleToAzimuthRangeTime(const ossimDpt & imPt, TimeType & azimuthTime, double & rangeTime) const
   {
379
      // std::clog << "ossimSarSensorModel::lineSampleToAzimuthRangeTime()\n";
380 381
      // First compute azimuth time here
      lineToAzimuthTime(imPt.y,azimuthTime);
382

383 384 385 386 387 388 389 390 391 392 393 394 395
      // Then compute range time
      if(isGRD())
      {
         // Handle grd case here
         double slantRange;
         groundRangeToSlantRange(imPt.x*theRangeResolution,azimuthTime, slantRange);
         rangeTime = theRangeTimeOffset + 2*slantRange/C;
      }
      else
      {
         rangeTime = theRangeTimeOffset + theNearRangeTime + imPt.x*(1/theRangeSamplingRate);
      }
   }
396

397 398
   void ossimSarSensorModel::computeRangeDoppler(const ossimEcefPoint & inputPt, const ossimEcefPoint & sensorPos, const ossimEcefVector sensorVel, double & range, double & doppler) const
   {
399
      // std::clog << "ossimSarSensorModel::computeRangeDoppler()\n";
400
      assert(theRadarFrequency>0&&"theRadarFrequency is null");
401

402 403
      // eq. 19, p. 25
      const ossimEcefVector s2gVec = inputPt - sensorPos;
404

405
      range = s2gVec.magnitude();
406

407
      const double coef = -2*C/(theRadarFrequency*range);
408

409 410
      doppler = coef * sensorVel.dot(s2gVec);
   }
411

412 413 414
   void ossimSarSensorModel::interpolateSensorPosVel(const TimeType & azimuthTime, ossimEcefPoint& sensorPos, ossimEcefVector& sensorVel, unsigned int deg) const
   {
      assert(!theOrbitRecords.empty()&&"The orbit records vector is empty");
415

416
      // Lagrangian interpolation of sensor position and velocity
417

418
      unsigned int nBegin(0), nEnd(0);
419

420 421 422
      sensorPos[0] = 0;
      sensorPos[1] = 0;
      sensorPos[2] = 0;
423

424 425 426
      sensorVel[0] = 0;
      sensorVel[1] = 0;
      sensorVel[2] = 0;
427

428 429
      // First, we search for the correct set of record to use during
      // interpolation
430

431 432
      // If there are less records than degrees, use them all
      if(theOrbitRecords.size()<deg)
433
      {
434
         nEnd = theOrbitRecords.size()-1;
435
      }
436
      else
437
      {
438 439
         // Search for the deg number of records around the azimuth time
         unsigned int t_min_idx = 0;
440
         DurationType t_min = abs(azimuthTime - theOrbitRecords.front().azimuthTime);
441

442
         unsigned int count = 0;
443

444 445
         for(std::vector<OrbitRecordType>::const_iterator it = theOrbitRecords.begin();it!=theOrbitRecords.end();++it,++count)
         {
446
            const DurationType current_time = abs(azimuthTime-it->azimuthTime);
447

448 449 450 451 452 453
            if(t_min > current_time)
            {
               t_min_idx = count;
               t_min = current_time;
            }
         }
454
         // TODO: see if these expressions can be simplified
455 456 457 458
         nBegin = std::max((int)t_min_idx-(int)deg/2+1,(int)0);
         nEnd = std::min(nBegin+deg-1,(unsigned int)theOrbitRecords.size());
         nBegin = nEnd<theOrbitRecords.size()-1 ? nBegin : nEnd-deg+1;
      }
459

460 461 462 463 464 465 466 467
      // Compute lagrangian interpolation using records from nBegin to nEnd
      for(unsigned int i = nBegin; i < nEnd; ++i)
      {
         double w = 1.;

         unsigned int j = nBegin;
         for( ; j != i ; ++j)
         {
468 469 470 471
            const DurationType td1 = azimuthTime                    - theOrbitRecords[j].azimuthTime;
            const DurationType td2 = theOrbitRecords[i].azimuthTime - theOrbitRecords[j].azimuthTime;
            const double f = td1 / td2;
            w *= f;
472 473 474 475
         }
         ++j;
         for( ; j < nEnd; ++j)
         {
476 477 478 479
            const DurationType td1 = azimuthTime                    - theOrbitRecords[j].azimuthTime;
            const DurationType td2 = theOrbitRecords[i].azimuthTime - theOrbitRecords[j].azimuthTime;
            const double f = td1 / td2;
            w *= f;
480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504
         }

         sensorPos[0]+=w*theOrbitRecords[i].position[0];
         sensorPos[1]+=w*theOrbitRecords[i].position[1];
         sensorPos[2]+=w*theOrbitRecords[i].position[2];

         sensorVel[0]+=w*theOrbitRecords[i].velocity[0];
         sensorVel[1]+=w*theOrbitRecords[i].velocity[1];
         sensorVel[2]+=w*theOrbitRecords[i].velocity[2];
      }
   }

   void ossimSarSensorModel::slantRangeToGroundRange(const double & slantRange, const TimeType & azimuthTime, double & groundRange) const
   {
      applyCoordinateConversion(slantRange,azimuthTime,theSlantRangeToGroundRangeRecords,groundRange);
   }

   void ossimSarSensorModel::groundRangeToSlantRange(const double & groundRange, const TimeType & azimuthTime, double & slantRange) const
   {
      applyCoordinateConversion(groundRange,azimuthTime,theGroundRangeToSlantRangeRecords,slantRange);
   }

   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.");
505
      // std::clog << "conv coord(" << in << ", az="<<azimuthTime<<")\n";
506

507 508
      // First, we need to find the correct pair of records for interpolation
      std::vector<CoordinateConversionRecordType>::const_iterator it = records.begin();
509

510
      CoordinateConversionRecordType srgrRecord;
511

512 513
      std::vector<CoordinateConversionRecordType>::const_iterator  previousRecord = it;
      ++it;
514

515
      std::vector<CoordinateConversionRecordType>::const_iterator nextRecord = it;
516

517
      // Look for the correct record
518
      // std::clog << "Looking for " << azimuthTime << " within records:\n";
519 520
      while(it!=records.end())
      {
521
         // std::clog << "- record: " << it->azimuthTime << "...";
522 523 524 525 526
         // nextRecord = it;

         if(azimuthTime >= previousRecord->azimuthTime
               && azimuthTime < nextRecord->azimuthTime)
         {
527
            // std::clog << " found!\n";
528
            break;
529 530 531
         }
         else
         {
532 533 534
            previousRecord = nextRecord;
            ++it;
            nextRecord = it;
535
            // std::clog << " NOT found => next!\n";
536 537 538 539 540 541 542
         }
      }
      assert(nextRecord == it);
      if(it == records.end())
      {
         if(azimuthTime < records.front().azimuthTime)
         {
543
            srgrRecord = records.front();
544
            // std::clog << "Not found, but before first => srgrRecord: " << srgrRecord.azimuthTime << "\n";
545 546 547
         }
         else if(azimuthTime >= records.back().azimuthTime)
         {
548
            srgrRecord = records.back();
549
            // std::clog << "Not found, but after last => srgrRecord: " << srgrRecord.azimuthTime << "\n";
550 551 552 553 554 555 556 557 558 559
         }
      }
      else
      {
         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
560
            = DurationType(azimuthTime             - previousRecord->azimuthTime)
561 562
            / (nextRecord->azimuthTime - previousRecord->azimuthTime)
            ;
563
         // std::clog << "interp: " << interp << " ="
564 565 566
         // << " (" << azimuthTime             << " - " << previousRecord->azimuthTime << " (="<< (azimuthTime             - previousRecord->azimuthTime)<< ") )"
         // << "/(" << nextRecord->azimuthTime << " - " << previousRecord->azimuthTime << " (="<< (nextRecord->azimuthTime - previousRecord->azimuthTime)<< ") )"
         // << "\n";
567

568
         srgrRecord.rg0 = (1-interp) * previousRecord->rg0 + interp*nextRecord->rg0;
569

570 571 572
         srgrRecord.coefs.clear();
         std::vector<double>::const_iterator pIt = previousRecord->coefs.begin();
         std::vector<double>::const_iterator nIt = nextRecord->coefs.begin();
573

574 575
         for(;pIt != previousRecord->coefs.end() && nIt != nextRecord->coefs.end();++pIt,++nIt)
         {
576
            srgrRecord.coefs.push_back(interp*(*nIt)+(1-interp)*(*pIt));
577
         }
578

579 580
         assert(!srgrRecord.coefs.empty()&&"Slant range to ground range interpolated coefficients vector is empty.");
      }
581

582 583 584
      // Now that we have the interpolated coefs, compute ground range
      // from slant range
      const double sr_minus_sr0 =  in-srgrRecord.rg0;
585

586
      assert(!srgrRecord.coefs.empty()&&"Slant range to ground range coefficients vector is empty.");
587

588
      out = 0;
589

590 591 592 593 594
      for(std::vector<double>::const_reverse_iterator cIt = srgrRecord.coefs.rbegin();cIt!=srgrRecord.coefs.rend();++cIt)
      {
         out = *cIt + sr_minus_sr0*out;
      }
   }
595 596


597 598
   bool ossimSarSensorModel::zeroDopplerLookup(const ossimEcefPoint & inputPt, TimeType & interpAzimuthTime, ossimEcefPoint & interpSensorPos, ossimEcefVector & interpSensorVel) const
   {
599
      assert((theOrbitRecords.size()>=2) && "Orbit records vector contains less than 2 elements");
600

601
      std::vector<OrbitRecordType>::const_iterator it = theOrbitRecords.begin();
602

603
      double doppler2(0.);
604

605 606 607
      // Compute range and doppler of first record
      // NOTE: here we only use the scalar product with vel and discard
      // the constant coef as it has no impact on doppler sign
608

609
      double doppler1 = (inputPt-it->position).dot(it->velocity);
610 611


612
      bool dopplerSign1 = doppler1 < 0;
613

614
      ++it; // -> it != begin
615

616 617 618 619 620 621
      // Look for the consecutive records where doppler freq changes sign
      // Note: implementing a bisection algorithm here might be faster
      for ( ; it!=theOrbitRecords.end() ; ++it)
      {
         // compute range and doppler of current record
         doppler2 = (inputPt-it->position).dot(it->velocity);
622

623
         const bool dopplerSign2 = doppler2 <0;
624

625 626 627 628 629 630 631 632 633
         // If a change of sign is detected
         if(dopplerSign1 != dopplerSign2)
         {
            break;
         }
         else
         {
            doppler1 = doppler2;
         }
634
      }
635

636
      // In this case, we need to extrapolate
637
      if(it == theOrbitRecords.end())
638
      {
639 640
         std::vector<OrbitRecordType>::const_iterator record1 = theOrbitRecords.begin();
         std::vector<OrbitRecordType>::const_iterator record2 = record1 + theOrbitRecords.size()-1;
641 642 643 644
         doppler1 = (inputPt-record1->position).dot(record1->velocity);
         doppler2 = (inputPt-record2->position).dot(record2->velocity);
         const DurationType delta_td = record2->azimuthTime - record1->azimuthTime;
         interpAzimuthTime = record1->azimuthTime - doppler1 / (doppler2 - doppler1) * delta_td;
645
      }
646
      else
647
      {
648 649 650 651
         assert(it != theOrbitRecords.begin());
         assert(it != theOrbitRecords.end());
         std::vector<OrbitRecordType>::const_iterator record2 = it;
         std::vector<OrbitRecordType>::const_iterator record1 = --it;
652 653 654
         // now interpolate time and sensor position
         const double abs_doppler1 = std::abs(doppler1);
         const double interpDenom = abs_doppler1+std::abs(doppler2);
655

656
         assert(interpDenom>0&&"Both doppler frequency are null in interpolation weight computation");
657

658
         const double interp = abs_doppler1/interpDenom;
659
         // std::clog << "interp: " << interp << "\n";
660

661
         const DurationType delta_td = record2->azimuthTime - record1->azimuthTime;
662
         // std::clog << "delta_td: " << delta_td << " = " << record2->azimuthTime <<" - " <<record1->azimuthTime<< "\n";
663

664
         // Compute interpolated time offset wrt record1
665
         // (No need for that many computations (day-frac -> ms -> day frac))
666
         const DurationType td     = delta_td * interp;
667
         // std::clog << "td: " << td  << "(" << td.total_microseconds() << "us)\n";
668
         // Compute interpolated azimuth time
669
         interpAzimuthTime = record1->azimuthTime + td + theAzimuthTimeOffset;
670
      }
671

672
      // std::clog << "interpAzimuthTime: " << interpAzimuthTime << "\n";
673

674 675
      // Interpolate sensor position and velocity
      interpolateSensorPosVel(interpAzimuthTime,interpSensorPos, interpSensorVel);
676

677 678
      return true;
   }
679

680 681 682
   void ossimSarSensorModel::computeBistaticCorrection(const ossimEcefPoint & inputPt, const ossimEcefPoint & sensorPos, DurationType & bistaticCorrection) const
   {
      // Bistatic correction (eq 25, p 28)
683
      double halftrange = 1000000. * (sensorPos-inputPt).magnitude()/C;
684 685
      bistaticCorrection= microseconds(static_cast<unsigned long>(floor(halftrange+0.5)));
   }
686 687


688 689 690
   void ossimSarSensorModel::azimuthTimeToLine(const TimeType & azimuthTime, double & line) const
   {
      assert(!theBurstRecords.empty()&&"Burst records are empty (at least one burst should be available)");
691

692
      std::vector<BurstRecordType>::const_iterator currentBurst = theBurstRecords.begin();
693

694 695 696 697 698 699 700 701 702 703 704 705 706
      // Look for the correct burst. In most cases the number of burst
      // records will be 1 (except for TOPSAR Sentinel1 products)
      std::vector<BurstRecordType>::const_iterator it = theBurstRecords.begin();
      std::vector<BurstRecordType>::const_iterator itend = theBurstRecords.end();
      for(; it!= itend ; ++it)
      {
         if(azimuthTime >= it->azimuthStartTime
               && azimuthTime < it->azimuthStopTime)
         {
            currentBurst = it;
            break;
         }
      }
707

708 709 710
      // If no burst is found, we will use the first (resp. last burst to
      // extrapolate line
      if(it == itend)
711
      {
712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727
         if(! theBurstRecords.empty())
         {
            if(azimuthTime < theBurstRecords.front().azimuthStartTime)
            {
               currentBurst = theBurstRecords.begin();
            }
            else if (azimuthTime > theBurstRecords.back().azimuthStopTime)
            {
               currentBurst = theBurstRecords.end()-1;
            }
         }
         else
         {
            // Fall back to the only record
            currentBurst = theBurstRecords.begin();
         }
728 729
      }

730
      const DurationType timeSinceStart = azimuthTime - currentBurst->azimuthStartTime;
731

732
      // Eq 22 p 27
733
      line = (timeSinceStart/theAzimuthTimeInterval) + currentBurst->startLine;
734
      // std::clog << "line = " << line << " <- " << timeSinceStart << "/" << theAzimuthTimeInterval << "+" << currentBurst->startLine << "\n";
735
   }
736

737 738 739
   void ossimSarSensorModel::lineToAzimuthTime(const double & line, TimeType & azimuthTime) const
   {
      assert(!theBurstRecords.empty()&&"Burst records are empty (at least one burst should be available)");
740

741
      std::vector<BurstRecordType>::const_iterator currentBurst = theBurstRecords.begin();
742

743 744 745 746 747 748 749 750
      if(theBurstRecords.size() != 1)
      {
         // Look for the correct burst. In most cases the number of burst
         // records will be 1 (except for TOPSAR Sentinel1 products)
         std::vector<BurstRecordType>::const_iterator it = theBurstRecords.begin();
         std::vector<BurstRecordType>::const_iterator itend = theBurstRecords.end();
         for( ; it!= itend; ++it)
         {
751 752
            if(line >= it->startLine && line < it->endLine)
            {
753 754
               currentBurst = it;
               break;
755
            }
756
         }
757

758 759
         if(it == itend)
         {
760 761
            if(line < theBurstRecords.front().startLine)
            {
762
               currentBurst = theBurstRecords.begin();
763 764 765
            }
            else if (line >= theBurstRecords.back().endLine)
            {
766
               currentBurst = theBurstRecords.end()-1;
767
            }
768
         }
769

770
      }
771

772 773
      const DurationType timeSinceStart = (line - currentBurst->startLine)*theAzimuthTimeInterval;
      // std::clog << "timeSinceStart: " << timeSinceStart.total_microseconds() << "us\n";
774

775
      // Eq 22 p 27
776 777
      azimuthTime = currentBurst->azimuthStartTime + timeSinceStart + theAzimuthTimeOffset;
      // std::clog << "offset: "         << theAzimuthTimeOffset << "\n";
778
      // std::clog << "->azimuthTime: "  << azimuthTime << "\n";
779
   }
780 781 782



783 784
   bool ossimSarSensorModel::projToSurface(const GCPRecordType & initGcp, const ossimDpt & target, const ossimHgtRef & hgtRef, ossimEcefPoint & ellPt) const
   {
785
     
786 787
      // Initialize current estimation
      ossimEcefPoint currentEstimation(initGcp.worldPt);
788

789
      // Compute corresponding image position
790
      // std::clog << "initGCP: " << initGcp.imPt << "\n";
791
      ossimDpt currentImPoint(initGcp.imPt);
792

793 794
      ossim_float64 currentImSquareResidual = squareDistance(target,currentImPoint);
      double currentHeightResidual = initGcp.worldPt.height() - hgtRef.getRefHeight(initGcp.worldPt);
795

796
      bool init = true;
797

798
      unsigned int iter = 0;
799

800 801 802 803 804 805 806
      // TODO: Every time the function is called, an allocation (+a free) is done.
      // This is not efficient. => Find a static matrix of 3x3 elements
      // Moreover, NEWMAT implies a lot of objet creations, hence allocations
      NEWMAT::SymmetricMatrix BtB(3);
      NEWMAT::ColumnVector BtF(3);
      NEWMAT::ColumnVector F(3);
      NEWMAT::ColumnVector dR(3);
807

808 809 810 811 812 813
      // Stop condition: img residual < 1e-2 pixels, height residual² <
      // 0.01² m, nb iter < 50. init ensure that loop runs at least once.
      while((init || (currentImSquareResidual > (0.01*0.01) || std::abs(currentHeightResidual) > 0.01))  && iter < 50)
      {
         if(init)
            init =false;
814

815
         // std::clog<<"Iter: "<<iter<<", Res: im="<<currentImSquareResidual<<", hgt="<<currentHeightResidual<<'\n';
816

817 818 819 820
         // compute residuals
         F(1) = target.x - currentImPoint.x;
         F(2) = target.y - currentImPoint.y;
         F(3) = currentHeightResidual;
821

822
         // std::clog<<"F("<<iter<<")="<<F<<'\n';
823

824 825
         // Delta use for partial derivatives estimation (in meters)
         const double d = 10.;
826

827 828 829
         // Compute partial derivatives
         ossimEcefVector p_fx, p_fy, p_fh,dx(d,0,0),dy(0,d,0),dz(0,0,d);
         ossimDpt tmpImPt;
830

831
         ossim_float64 rdx(0.0),rdy(0.0),rdz(0.0), fdx(0.0),fdy(0.0),fdz(0.0);
832

833 834 835
         ossimGpt currentEstimationWorld(currentEstimation);
         ossimGpt tmpGpt = ossimGpt(currentEstimation+dx);
         worldToLineSample(tmpGpt,tmpImPt);
836 837 838 839
         // std::clog << "currentEstimationWorld: " << currentEstimationWorld << "\n";
         // std::clog << "currentEstimation: " << currentEstimation << "\n";
         // std::clog << "tmpGpt: " << tmpGpt << "\n";
         // std::clog << "tmpImPt: " << tmpImPt << "\n";
840 841 842
         p_fx[0] = (currentImPoint.x-tmpImPt.x)/d;
         p_fy[0] = (currentImPoint.y-tmpImPt.y)/d;
         p_fh[0] = (currentEstimationWorld.height()-tmpGpt.height())/d;
843

844 845 846 847 848
         tmpGpt = ossimGpt(currentEstimation+dy);
         worldToLineSample(tmpGpt,tmpImPt);
         p_fx[1] = (currentImPoint.x-tmpImPt.x)/d;
         p_fy[1] = (currentImPoint.y-tmpImPt.y)/d;
         p_fh[1] = (currentEstimationWorld.height()-tmpGpt.height())/d;
849

850 851 852 853 854
         tmpGpt = ossimGpt(currentEstimation+dz);
         worldToLineSample(tmpGpt,tmpImPt);
         p_fx[2] = (currentImPoint.x-tmpImPt.x)/d;
         p_fy[2] = (currentImPoint.y-tmpImPt.y)/d;
         p_fh[2] = (currentEstimationWorld.height()-tmpGpt.height())/d;
855

856 857 858 859
         // Form B-matrix
         NEWMAT::Matrix B = ossimMatrix3x3::create(p_fx[0], p_fx[1], p_fx[2],
               p_fy[0], p_fy[1], p_fy[2],
               p_fh[0], p_fh[1], p_fh[2]);
860

861
         // std::clog<<"B: "<<B<<'\n';
862

863 864 865
         // Invert system
         try {
            dR = B.i() * F;
866
         } catch (NEWMAT::SingularException const& e) {
867 868 869 870 871 872
         ellPt = currentEstimation;
         
         ossimNotify(ossimNotifyLevel_WARN) <<"ossim::SarSensorModel::projToSurface(): singular matrix can not be inverted. Returning best estimation so far ("<<ellPt<<") for output point ("<<target<<")\n";
         std::clog << "initGCP: " << initGcp.imPt <<", "<<initGcp.worldPt<< "\n";
         
         return true;
873
         }
874

875 876 877 878 879
         // Update estimate
         for (ossim_int32 k=0; k<3; k++)
         {
            currentEstimation[k] -= dR(k+1);
         }
880

881
         // std::clog<<"dR: "<<dR<<'\n';
882

883
         currentEstimationWorld=ossimGpt(currentEstimation);
884

885 886 887
         // Update residuals
         const ossim_float64 atHgt = hgtRef.getRefHeight(currentEstimationWorld);
         currentHeightResidual = atHgt - currentEstimationWorld.height();
888

889
         worldToLineSample(currentEstimationWorld,currentImPoint);
890

891
         // std::clog<<currentImPoint<<'\n';
892

893
         currentImSquareResidual = squareDistance(currentImPoint,target);
894

895 896
         ++iter;
      }
897

898
      // std::clog<<"Iter: "<<iter<<", Res: im="<<currentImSquareResidual<<", hgt="<<currentHeightResidual<<'\n';
899

900 901 902
      ellPt = currentEstimation;
      return true;
   }
903 904 905



906 907 908 909 910 911 912 913 914 915 916
   //*************************************************************************************************
   // Infamous DUP
   //*************************************************************************************************
   ossimObject* ossimSarSensorModel::dup() const
   {
      return new ossimSarSensorModel(*this);
   }
   bool ossimSarSensorModel::useForward() const
   {
      return false;
   }
917

918 919
   bool ossimSarSensorModel::autovalidateInverseModelFromGCPs(const double & xtol, const double & ytol, const double azTimeTol, const double & rangeTimeTol) const
   {
920
      std::clog << "ossimSarSensorModel::autovalidateInverseModelFromGCPs()\n";
921 922 923 924
      if(theGCPRecords.empty())
      {
         return false;
      }
925

926
      bool success = true;
927

928
      unsigned int gcpId = 1;
929

930
      std::clog << theGCPRecords.size() << " GCPS\n";
931 932 933 934 935 936 937
      for(std::vector<GCPRecordType>::const_iterator gcpIt = theGCPRecords.begin(); gcpIt!=theGCPRecords.end();++gcpIt,++gcpId)
      {
         ossimDpt estimatedImPt;
         TimeType estimatedAzimuthTime;
         double   estimatedRangeTime;

         // Estimate times
938 939 940
         ossimEcefPoint sensorPos;
         ossimEcefVector sensorVel;
         const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime,sensorPos, sensorVel);
941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956
         this->worldToLineSample(gcpIt->worldPt,estimatedImPt);

         const bool thisSuccess
            =  s1
            && (std::abs(estimatedImPt.x - gcpIt->imPt.x) <= xtol)
            && (std::abs(estimatedImPt.y - gcpIt->imPt.y) <= ytol)
            && (std::abs((estimatedAzimuthTime-gcpIt->azimuthTime).total_microseconds()) <= azTimeTol)
            && (std::abs(estimatedRangeTime - gcpIt->slantRangeTime) <= rangeTimeTol)
            ;

         const bool verbose = k_verbose;

         success = success && thisSuccess;

         if(verbose)
         {
957
            std::clog<<"GCP #"<<gcpId<< (thisSuccess ? "succeeded" : "failed") << '\n';
958 959 960 961
            std::clog<<"Azimuth time: ref="<<gcpIt->azimuthTime<<", predicted: "<<estimatedAzimuthTime<<", res="<<to_simple_string(estimatedAzimuthTime-gcpIt->azimuthTime)<<'\n';
            std::clog<<"Slant range time: ref="<<gcpIt->slantRangeTime<<", predicted: "<<estimatedRangeTime<<", res="<<std::abs(estimatedRangeTime - gcpIt->slantRangeTime)<<'\n';
            std::clog<<"Image point: ref="<<gcpIt->imPt<<", predicted="<<estimatedImPt<<", res="<<estimatedImPt-gcpIt->imPt<<'\n';
            std::clog<<'\n';
962 963
         }
      }
964

965
      if(success)
966
      {
967
         std::cout<<"All GCPs within "<<ytol <<" azimuth pixel, "<<xtol<<" range pixel, "<<azTimeTol<<" microseconds of azimuth time, "<<rangeTimeTol<<" of range time\n";
968 969
      }

970 971 972 973 974
      return success;
   }

   bool ossimSarSensorModel::autovalidateForwardModelFromGCPs(double resTol)
   {
975
      // std::clog << "ossimSarSensorModel::autovalidateForwardModelFromGCPs()\n";
976
      resTol *= resTol; // as internally we won't be using sqrt on norms
977

978 979 980
      // First, split half of the gcps to serve as tests, and remove them
      // temporarily from theGCPRecord.
      std::vector<GCPRecordType> gcpRecordSave, testGcps;
981

982
      gcpRecordSave.swap(theGCPRecords); // steal the data; `gcpRecordSave = move(theGCPRecords);` in C++11
983

984 985 986 987 988 989 990
      unsigned int count = 0;

      unzip(gcpRecordSave, theGCPRecords, testGcps);
      assert(theGCPRecords.size() >= testGcps.size());

      bool success = true;
      const bool verbose = k_verbose;
991

992
      unsigned int gcpId = 1;
993

994
      // std::clog << testGcps.size() << " GCPS\n";
995
      for(std::vector<GCPRecordType>::const_iterator gcpIt = testGcps.begin(); gcpIt!=testGcps.end();++gcpIt,++gcpId)
996
      {
997 998 999 1000 1001 1002 1003 1004 1005
         ossimGpt estimatedWorldPt;
         ossimGpt const& refPt = gcpIt->worldPt;

         double estimatedRangeTime;
         TimeType estimatedAzimuthTime;

         lineSampleToAzimuthRangeTime(gcpIt->imPt,estimatedAzimuthTime,estimatedRangeTime);

         lineSampleHeightToWorld(gcpIt->imPt,refPt.height(),estimatedWorldPt);
1006

1007 1008 1009 1010 1011 1012 1013 1014
         const double res = squareDistance(refPt, estimatedWorldPt);

         if(res>resTol || estimatedWorldPt.hasNans())
         {
            success = false;

            if(verbose)
            {
1015 1016 1017 1018 1019 1020
               std::clog<<"GCP #"<<gcpId<<'\n';
               std::clog<<"Azimuth time: ref="<<gcpIt->azimuthTime<<", predicted: "<<estimatedAzimuthTime<<", res="<<to_simple_string(estimatedAzimuthTime-gcpIt->azimuthTime)<<'\n';
               std::clog<<"Slant range time: ref="<<gcpIt->slantRangeTime<<", predicted: "<<estimatedRangeTime<<", res="<<std::abs(estimatedRangeTime - gcpIt->slantRangeTime)<<'\n';
               std::clog<<"Im point: "<<gcpIt->imPt<<'\n';
               std::clog<<"World point: ref="<<refPt<<", predicted="<<estimatedWorldPt<<", res="<<sqrt(res)<<" m\n";
               std::clog<<'\n';
1021 1022 1023
            }
         }
      }
1024

1025
      theGCPRecords.swap(gcpRecordSave);
1026

1027 1028
      return success;
   }
1029

1030 1031
   void ossimSarSensorModel::optimizeTimeOffsetsFromGcps()
   {
1032
      // std::clog << "ossimSarSensorModel::optimizeTimeOffsetsFromGcps()\n";
1033
      DurationType cumulAzimuthTime(seconds(0));
1034
      double cumulRangeTime(0);
1035
      unsigned int count=0;
1036 1037 1038
      // reset offsets before optimisation
      theAzimuthTimeOffset = seconds(0);
      theRangeTimeOffset = 0.0;
1039

1040 1041
      // First, fix the azimuth time
      for(std::vector<GCPRecordType>::const_iterator gcpIt = theGCPRecords.begin(); gcpIt!=theGCPRecords.end();++gcpIt)
1042
      {
1043 1044 1045 1046
         ossimDpt estimatedImPt;
         TimeType estimatedAzimuthTime;
         double   estimatedRangeTime;

1047 1048
         ossimEcefPoint sensorPos;
         ossimEcefVector sensorVel;
1049
         // Estimate times
1050
         const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime, sensorPos, sensorVel);
1051 1052 1053

         if(s1)
         {
1054
            cumulAzimuthTime -= (estimatedAzimuthTime-gcpIt->azimuthTime);
1055 1056
            ++count;
         }
1057 1058
      }

1059
      theAzimuthTimeOffset = count > 0 ? cumulAzimuthTime / count : DurationType(0);
1060

1061 1062 1063 1064 1065 1066 1067 1068 1069
      // Then, fix the range time
      count=0;

      for(std::vector<GCPRecordType>::const_iterator gcpIt = theGCPRecords.begin(); gcpIt!=theGCPRecords.end();++gcpIt)
      {
         ossimDpt estimatedImPt;
         TimeType estimatedAzimuthTime;
         double   estimatedRangeTime;

1070 1071 1072
         ossimEcefPoint sensorPos;
         ossimEcefVector sensorVel;

1073
         // Estimate times
1074
         const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime, sensorPos, sensorVel);
1075 1076 1077 1078 1079 1080 1081 1082

         if(s1)
         {
            cumulRangeTime+=-estimatedRangeTime+gcpIt->slantRangeTime;
            ++count;
         }
      }

1083
      theRangeTimeOffset = count > 0 ? cumulRangeTime/count : 0;
1084 1085
   }

1086 1087 1088 1089 1090
   void get(
         ossimKeywordlist                             const& kwl,
         std::vector<ossimSarSensorModel::OrbitRecordType> & orbitRecords)
   {
      char orbit_prefix_[256];
1091
      std::size_t nbOrbits(0);
1092 1093 1094
      try {
         get(kwl, "orbitList.nb_orbits", nbOrbits);
      } catch (kw_runtime_error const& e) {
1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108
         nbOrbits = 0;
         ossimRegExp regExp;
         regExp.compile("orbitList\\.orbit\\[.*\\]\\.time");
         ossimKeywordlist::KeywordMap::const_iterator i =
           kwl.getMap().begin();
         for(; i != kwl.getMap().end(); ++i)
         {
            if(regExp.find( (*i).first.c_str()))
            {
               ++nbOrbits;
            }
         }
         // Method getNumberOfKeysThatMatch not available in ossim 1.8.16
         //nbOrbits = kwl.getNumberOfKeysThatMatch("orbitList\\.orbit\\[.*\\]\\.time");
1109 1110 1111 1112 1113 1114
         ossimNotify(ossimNotifyLevel_WARN)
            << "WARNING: " << e.what()
            << "\n\tNumber of orbits manually counted to " << nbOrbits
            << ".\n\tPlease update your geom file.\n";
      }

1115
      for (std::size_t i=0; i!=nbOrbits ; ++i) {
1116
         const int pos = s_printf(orbit_prefix_, "orbitList.orbit[%d].", int(i));
1117 1118 1119 1120 1121 1122 1123 1124
         assert(pos > 0 && pos < 256);
         const std::string orbit_prefix(orbit_prefix_, pos);

         ossimSarSensorModel::OrbitRecordType orbitRecord;
         get(kwl, orbit_prefix + keyTime, orbitRecord.azimuthTime);
         get(kwl, orbit_prefix + keyPosX, orbitRecord.position[0]);
         get(kwl, orbit_prefix + keyPosY, orbitRecord.position[1]);
         get(kwl, orbit_prefix + keyPosZ, orbitRecord.position[2]);
1125
         get(kwl, orbit_prefix + keyVelX, orbitRecord.velocity[0]);
1126 1127 1128 1129
         get(kwl, orbit_prefix + keyVelY, orbitRecord.velocity[1]);
         get(kwl, orbit_prefix + keyVelZ, orbitRecord.velocity[2]);
         orbitRecords.push_back(orbitRecord);
      }
1130 1131
   }

1132 1133 1134 1135 1136 1137
   void add(
       ossimKeywordlist & kwl,
       const std::vector<ossimSarSensorModel::OrbitRecordType> & orbitRecords)
   {
     char orbit_prefix_[256];

1138
     add(kwl,ORBIT_NUMBER_KEY,(ossim_uint32)orbitRecords.size());
1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155

     for(std::size_t i = 0; i!=orbitRecords.size();++i)
       {
       const int pos = s_printf(orbit_prefix_, "%s[%d].", ORBIT_PREFIX.c_str(), int(i));
       assert(pos > 0 && pos < 256);
       const std::string orbit_prefix(orbit_prefix_, pos);

       add(kwl, orbit_prefix + keyTime, orbitRecords[i].azimuthTime);
       add(kwl, orbit_prefix + keyPosX, orbitRecords[i].position[0]);
       add(kwl, orbit_prefix + keyPosY, orbitRecords[i].position[1]);
       add(kwl, orbit_prefix + keyPosZ, orbitRecords[i].position[2]);
       add(kwl, orbit_prefix + keyVelX, orbitRecords[i].velocity[0]);
       add(kwl, orbit_prefix + keyVelY, orbitRecords[i].velocity[1]);
       add(kwl, orbit_prefix + keyVelZ, orbitRecords[i].velocity[2]);
       }
   }

1156 1157 1158
   void get(
         ossimKeywordlist                             const& kwl,
         std::vector<ossimSarSensorModel::BurstRecordType> & burstRecords)
1159
   {
1160
      char burstPrefix_[1024];
1161
      std::size_t nbBursts(0);
1162 1163
      get(kwl, BURST_NUMBER_KEY, nbBursts);
      for (std::size_t burstId=0; burstId!=nbBursts ; ++burstId) {
1164
         const int pos = s_printf(burstPrefix_, "%s[%d].", BURST_PREFIX.c_str(), burstId);
1165
         assert(pos > 0 && pos < sizeof(burstPrefix_));
1166
         const std::string burstPrefix(burstPrefix_, pos);
1167 1168 1169 1170

         ossimSarSensorModel::BurstRecordType burstRecord;
         get(kwl, burstPrefix + keyStartLine,        burstRecord.startLine);
         get(kwl, burstPrefix + keyEndLine,          burstRecord.endLine);
1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186
	 
	 try {
            unsigned int version;
            get(kwl, HEADER_PREFIX, "version", version);
	    // startSample and endSample since version 3
            if (version >= 3) 
	      {
                get(kwl, burstPrefix + keyStartSample,        burstRecord.startSample);
		get(kwl, burstPrefix + keyEndSample,          burstRecord.endSample);
	      }
         } 
	 catch (...) {
            throw std::runtime_error("Geom file generated with previous version of ossim plugins");
         }
        
	 get(kwl, burstPrefix + keyAzimuthStartTime, burstRecord.azimuthStartTime);
1187 1188 1189
         get(kwl, burstPrefix + keyAzimuthStopTime,  burstRecord.azimuthStopTime);
         burstRecords.push_back(burstRecord);
      }
1190 1191
   }

1192 1193 1194 1195 1196
   void add(
     ossimKeywordlist                                        & kwl,
     const std::vector<ossimSarSensorModel::BurstRecordType> & burstRecords)
   {
     char burstPrefix_[1024];
1197
     add(kwl, BURST_NUMBER_KEY, (ossim_uint32)burstRecords.size());
1198 1199 1200 1201 1202
     for (std::size_t burstId=0; burstId!=burstRecords.size() ; ++burstId) {
     const int pos = s_printf(burstPrefix_, "%s[%d].", BURST_PREFIX.c_str(), burstId);
     assert(pos > 0 && pos < sizeof(burstPrefix_));
     const std::string burstPrefix(burstPrefix_, pos);

1203 1204
     add(kwl, burstPrefix + keyStartLine, (ossim_uint32) burstRecords[burstId].startLine);
     add(kwl, burstPrefix + keyEndLine, (ossim_uint32) burstRecords[burstId].endLine);
1205 1206
     add(kwl, burstPrefix + keyStartSample, (ossim_uint32) burstRecords[burstId].startSample);
     add(kwl, burstPrefix + keyEndSample, (ossim_uint32) burstRecords[burstId].endSample);
1207 1208 1209 1210 1211
     add(kwl, burstPrefix + keyAzimuthStartTime, burstRecords[burstId].azimuthStartTime);
     add(kwl, burstPrefix + keyAzimuthStopTime,  burstRecords[burstId].azimuthStopTime);
     }
   }

1212 1213 1214 1215 1216
   void get(
         ossimKeywordlist                             const& kwl,
         std::vector<ossimSarSensorModel::GCPRecordType> & gcpRecords)
   {
      char prefix_[1024];
1217
      std::size_t nbGCPs(0);
1218 1219
      get(kwl, GCP_NUMBER_KEY, nbGCPs);
      for (std::size_t gcpId=0; gcpId