MRG: Merge branch 'develop' into complexImage_integration

.gitlab/issue_templates/feature_request.md

 Short summary of the requested feature
+
+/label ~feature

CONTRIBUTING.md

@@ -89,10 +89,33 @@ Request against the `develop` branch on GitLab using the merge request
 template. The merge request will then be discussed by the community and the core
 OTB team.
 
-* Merge requests can not be merged until all discussions have been resolved (this is enforced by GitLab).
-* Merge requests **must receive at least 2 positives votes from PSC members** before being merged.
-* The merger is responsible for checking that the branch is up-to-date with develop and that the dashboard is ok.
+* Merge requests can not be merged until all discussions have been resolved (this is enforced by GitLab)
+* Merge requests **must receive at least 2 positives votes from PSC members** before being merged
+* The merger is responsible for checking that the branch is up-to-date with develop
 * Merge requests can be merged by anyone (not just PSC or RM) with push access to develop
+* Merge requests can be merged once the dashboard is proven green for this branch
+
+Branches can be registered for dashboard testing by adding one line in `Config/feature_branches.txt` in [otb-devutils repository](https://gitlab.orfeo-toolbox.org/orfeotoolbox/otb-devutils.git).
+
+For branches in the main repository, the syntax is the following:
+
+```
+branch_name [otb-data_branch_name]
+
+```
+The second branch name is optional. It can be set if you need to modify [otb-data](https://gitlab.orfeo-toolbox.org/orfeotoolbox/otb-data.git) according to your changes.
+
+For branches in forks, the syntax is the following:
+```
+user/branch_name [user/otb-data_branch_name]
+```
+Again, the second branch name is optional.
+
+For users without push access to [otb-devutils repository](https://gitlab.orfeo-toolbox.org/orfeotoolbox/otb-devutils.git), the modification can be asked through a merge requests to this repository.
+
+Once the feature branch is registered for testing, it should appear in the *FeatureBranches* section of the [OTB dashboard](https://dash.orfeo-toolbox.org/index.php?project=OTB) next day (remember tests are run on a nighlty basis).
+
+Do not forget to remove the feature branch for testing once it has been merged.
 
 ## Remote modules
 
@@ -104,3 +127,21 @@ OTB source code. Under some conditions (dependencies, official acceptance
 process, etc.), we are also able to distribute your remote module in the
 official standalone binaries. See [the wiki](https://wiki.orfeo-toolbox.org/index.php/Remote_Modules)
 for more information.
+
+## Gitlab guidelines
+
+In order to organize the issues in our Gitlab instance, we use both labels and
+milestones.
+
+The [milestones](https://gitlab.orfeo-toolbox.org/orfeotoolbox/otb/milestones) should be used to track in which release a feature is merged.
+Gitlab can then provide a summary of all features and bugs added to a given release
+version.
+
+Regarding labels, we use the following set:
+* ~story: significant feature to be implemented with a detailed work plan, it can
+  correspond to a Request for Comments that has turned into a development action
+* ~bug: Bug, crash or unexpected behavior, reported by a user or a developer
+* ~feature: Feature request expressed by an OTB user/developer
+* ~"To Do": action is planned
+* ~Doing: work in progress
+* ~api ~app ~documentation ~monteverdi ~packaging ~qgis: optional context information

Modules/Adapters/OSSIMAdapters/include/otbSarSensorModelAdapter.h

@@ -24,6 +24,7 @@
 #include <memory>
 
 #include "otbDEMHandler.h"
+#include "itkPoint.h"
 
 namespace ossimplugins
 {
@@ -62,6 +63,9 @@ public:
 
   typedef std::auto_ptr<ossimplugins::ossimSarSensorModel> InternalModelPointer;
 
+  using Point2DType = itk::Point<double,2>;
+  using Point3DType = itk::Point<double,3>;
+  
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
 
@@ -80,6 +84,21 @@ public:
   /** Deburst metadata if possible and return lines to keep in image file */
   bool Deburst(std::vector<std::pair<unsigned long, unsigned long> > & lines);
 
+  /** Transform world point (lat,lon,hgt) to input image point
+  (col,row) and YZ frame */
+  bool WorldToLineSampleYZ(const Point3DType & inGeoPoint, Point2DType & cr, Point2DType & yz) const;
+
+  /** Transform world point (lat,lon,hgt) to input image point
+  (col,row) */
+  bool WorldToLineSample(const Point3DType & inGEoPOint, Point2DType & cr) const;
+
+/** Transform world point (lat,lon,hgt) to satellite position (x,y,z) and satellite velocity */
+  bool WorldToSatPositionAndVelocity(const Point3DType & inGeoPoint, Point3DType & satelitePosition,  
+				     Point3DType & sateliteVelocity) const;
+
+  /** Transform world point (lat,lon,hgt) to cartesian point (x,y,z) */
+  static bool WorldToCartesian(const Point3DType & inGeoPoint, Point3DType & outCartesianPoint);
+
   static bool ImageLineToDeburstLine(const std::vector<std::pair<unsigned long,unsigned long> >& lines, unsigned long imageLine, unsigned long & deburstLine);
 
   static void DeburstLineToImageLine(const std::vector<std::pair<unsigned long,unsigned long> >& lines, unsigned long deburstLine, unsigned long & imageLine);

Modules/Adapters/OSSIMAdapters/src/otbImageKeywordlist.cxx

@@ -287,60 +287,15 @@ ReadGeometryFromImage(const std::string& filename, bool checkRpcTag)
     }
 
   /**********************************************************/
-  /* Third try : look for external geom file and RPC tags   */
+  /* Third try : look for RPC tags   */
   /**********************************************************/
-  if (!hasMetaData)
+  if (!hasMetaData && checkRpcTag)
     {
-    // If still no metadata, try the ".geom" file
-    ossimFilename ossimGeomFile = ossimFilename(filename).setExtension(".geom");
-    otb_kwl = ReadGeometryFromGEOMFile(ossimGeomFile);
-
-    // also check any RPC tags
-    ImageKeywordlist rpc_kwl;
-    if (checkRpcTag)
-      {
-      rpc_kwl = ReadGeometryFromRPCTag(filename);
-      }
-
-    if (otb_kwl.HasKey("type"))
-      {
-      // external geom has a "type" keyword
-      std::string geomType(otb_kwl.GetMetadataByKey("type"));
-      ossimProjection *testProj = ossimProjectionFactoryRegistry::instance()
-        ->createProjection(ossimString(geomType));
-      if (dynamic_cast<ossimSensorModel*>(testProj))
-        {
-        // "type" keyword corresponds to a sensor model : don't erase it
-        if (rpc_kwl.GetSize() > 0)
-          {
-          rpc_kwl.ClearMetadataByKey("type");
-          }
-
-        ossimKeywordlist ossim_test_kwl;
-        otb_kwl.convertToOSSIMKeywordlist(ossim_test_kwl);
-        testProj = ossimProjectionFactoryRegistry::instance()
-          ->createProjection(ossim_test_kwl);
-        if (testProj)
-          {
-          // external geom contains a valid sensor geometry
-          hasMetaData = true;
-          }
-        }
-      }
+    // check any RPC tags
+    otb_kwl = ReadGeometryFromRPCTag(filename);
 
-    // copy keywords found in RPC tags if the external geom is not valid
-    if (!hasMetaData && rpc_kwl.GetSize() > 0)
+    if (!otb_kwl.Empty())
       {
-      const ImageKeywordlist::KeywordlistMap &kwlMap = rpc_kwl.GetKeywordlist();
-      for (ImageKeywordlist::KeywordlistMap::const_iterator it = kwlMap.begin();
-          it != kwlMap.end();
-          ++it)
-        {
-        if (it->second != "")
-          {
-          otb_kwl.AddKey(it->first , it->second);
-          }
-        }
       hasMetaData = true;
       }
     }

Modules/Adapters/OSSIMAdapters/src/otbSarSensorModelAdapter.cxx

@@ -109,7 +109,112 @@ void SarSensorModelAdapter::DeburstLineToImageLine(const std::vector<std::pair<u
   ossimplugins::ossimSarSensorModel::deburstLineToImageLine(lines,deburstLine,imageLine);
 }
 
+bool SarSensorModelAdapter::WorldToLineSampleYZ(const Point3DType & inGeoPoint, Point2DType & cr, Point2DType & yz) const
+{
+  if(m_SensorModel.get() == ITK_NULLPTR)
+    {
+    return false;
+    }
+
+  ossimGpt inGpt;
+  inGpt.lon = inGeoPoint[0];
+  inGpt.lat = inGeoPoint[1];
+  inGpt.hgt = inGeoPoint[2];
+
+  ossimDpt outDpt;
+
+  double y(0.),z(0.);
+  m_SensorModel->worldToLineSampleYZ(inGpt,outDpt,y,z);
+
+  if(outDpt.isNan())
+    return false;
+
+  cr[0]=outDpt.x;
+  cr[1]=outDpt.y;
+  yz[0]=y;
+  yz[1]=z;
+
+  return true;
+}
 
+bool SarSensorModelAdapter::WorldToLineSample(const Point3DType & inGeoPoint, Point2DType & cr) const
+{
+  if(m_SensorModel.get() == ITK_NULLPTR)
+    {
+    return false;
+    }
+
+  ossimGpt inGpt;
+  inGpt.lon = inGeoPoint[0];
+  inGpt.lat = inGeoPoint[1];
+  inGpt.hgt = inGeoPoint[2];
+
+  ossimDpt outDpt;
+
+  m_SensorModel->worldToLineSample(inGpt,outDpt);
+
+  if(outDpt.isNan())
+    return false;
+
+  cr[0]=outDpt.x;
+  cr[1]=outDpt.y;
+
+  return true;
+}
+
+bool SarSensorModelAdapter::WorldToCartesian(const Point3DType & inGeoPoint, Point3DType & outCartesianPoint)
+{
+  ossimGpt inGpt;
+  inGpt.lon = inGeoPoint[0];
+  inGpt.lat = inGeoPoint[1];
+  inGpt.hgt = inGeoPoint[2];
+
+  ossimEcefPoint outCartesien(inGpt);
+  
+  if(outCartesien.isNan())
+    return false;
+
+  outCartesianPoint[0] = outCartesien.x();
+  outCartesianPoint[1] = outCartesien.y();
+  outCartesianPoint[2] = outCartesien.z();
+
+  return true;
+}
+  
+bool SarSensorModelAdapter::WorldToSatPositionAndVelocity(const Point3DType & inGeoPoint, 
+							  Point3DType & satelitePosition, 
+							  Point3DType & sateliteVelocity) const
+{
+  if(m_SensorModel.get() == ITK_NULLPTR)
+    {
+      return false;
+    }
+
+  ossimGpt inGpt;
+  inGpt.lon = inGeoPoint[0];
+  inGpt.lat = inGeoPoint[1];
+  inGpt.hgt = inGeoPoint[2];
+  
+  ossimplugins::ossimSarSensorModel::TimeType azimuthTime;
+  double rangeTime;
+  ossimEcefPoint sensorPos;
+  ossimEcefVector sensorVel;
+
+  const bool success = m_SensorModel->worldToAzimuthRangeTime(inGpt, azimuthTime, rangeTime,sensorPos,sensorVel);
+  
+  if(sensorPos.isNan() || !success)
+    return false;
+
+  satelitePosition[0] = sensorPos.x();
+  satelitePosition[1] = sensorPos.y();
+  satelitePosition[2] = sensorPos.z();
+
+  sateliteVelocity[0] = sensorVel.x();
+  sateliteVelocity[1] = sensorVel.y();
+  sateliteVelocity[2] = sensorVel.z();    
+
+  return true;
+}
 
 
 } // namespace otb

Modules/Adapters/OSSIMAdapters/test/CMakeLists.txt

@@ -31,6 +31,7 @@ otbGeometricSarSensorModelAdapter.cxx
 otbPlatformPositionAdapter.cxx
 otbDEMHandlerTest.cxx
 otbRPCSolverAdapterTest.cxx
+otbSarSensorModelAdapterTest.cxx
 )
 
 add_executable(otbOSSIMAdaptersTestDriver ${OTBOSSIMAdaptersTests})
@@ -499,3 +500,8 @@ otb_add_test(NAME uaTvRPCSolverAdapterValidationTest COMMAND otbOSSIMAdaptersTes
   ${INPUTDATA}/DEM/egm96.grd
   )
 
+
+otb_add_test(NAME uaTvSarSensorModelAdapter COMMAND otbOSSIMAdaptersTestDriver
+  otbSarSensorModelAdapterTest
+  ${INPUTDATA}/s1a-iw1-slc-vh-amp_xt.geom
+  )

Modules/Adapters/OSSIMAdapters/test/otbOSSIMAdaptersTestDriver.cxx

@@ -33,4 +33,5 @@ void RegisterTests()
   REGISTER_TEST(otbPlatformPositionComputeBaselineTest);
   REGISTER_TEST(otbDEMHandlerTest);
   REGISTER_TEST(otbRPCSolverAdapterTest);
+  REGISTER_TEST(otbSarSensorModelAdapterTest);
 }

Modules/Adapters/OSSIMAdapters/test/otbSarSensorModelAdapterTest.cxx

+/*
+ * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ *     https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+
+#include <fstream>
+#include <iomanip>
+
+#include "otbSarSensorModelAdapter.h"
+#include "otbImageKeywordlist.h"
+
+
+int otbSarSensorModelAdapterTest(int itkNotUsed(argc), char* argv[])
+{
+  std::string infname = argv[1];
+
+  otb::SarSensorModelAdapter::Pointer sensorModel = otb::SarSensorModelAdapter::New();
+
+  auto kwl = otb::ReadGeometryFromGEOMFile(infname);
+
+  bool success = sensorModel->LoadState(kwl);
+
+  if(!success)
+    {
+    std::cerr<<"Could not LoadState() from keyword list read from"<<infname<<std::endl;
+    return EXIT_FAILURE;
+    }
+
+  std::vector<std::pair<unsigned long, unsigned long> > lines;
+  success = sensorModel->Deburst(lines);
+
+  if(!success)
+    {
+    std::cerr<<"Deburst() call failed."<<std::endl;
+    return EXIT_FAILURE;
+    }
+  
+
+  otb::ImageKeywordlist outKwl;
+  success = sensorModel->SaveState(outKwl);
+  if(!success)
+    {
+    std::cerr<<"SaveState() call failed."<<std::endl;
+    return EXIT_FAILURE;
+    }
+  
+
+  otb::SarSensorModelAdapter::Point2DType out1,out2;
+  otb::SarSensorModelAdapter::Point3DType in, out3, out4, out5;
+
+  // GCP 99 from input geom file
+  //support_data.geom.gcp[99].world_pt.hgt:  2.238244926818182e+02
+  //support_data.geom.gcp[99].world_pt.lat:  4.323458093295080e+01
+  //support_data.geom.gcp[99].world_pt.lon:  1.116316013091967e+00
+    
+  in[0] = 4.323458093295080e+01;
+  in[1] = 1.116316013091967e+00;
+  in[2] = 2.238244926818182e+02;
+
+  sensorModel->WorldToLineSample(in,out1);
+  sensorModel->WorldToLineSampleYZ(in,out1,out2);
+
+  sensorModel->WorldToCartesian(in, out5);
+  sensorModel->WorldToSatPositionAndVelocity(in,out3, out4);
+  
+  
+  return EXIT_SUCCESS;
+}

Modules/IO/ImageIO/include/otbImageFileReader.txx

@@ -398,18 +398,33 @@ ImageFileReader<TOutputImage, ConvertPixelTraits>
     {
 
     std::string lFileNameOssimKeywordlist = GetDerivedDatasetSourceFileName(m_FileName);
+    std::string extension = itksys::SystemTools::GetFilenameLastExtension(lFileNameOssimKeywordlist);
+    std::string attachedGeom = lFileNameOssimKeywordlist.substr(
+      0,
+      lFileNameOssimKeywordlist.size() - extension.size()) + std::string(".geom");
 
     // Update otb Keywordlist
     ImageKeywordlist otb_kwl;
-    if (!m_FilenameHelper->ExtGEOMFileNameIsSet())
+
+    // Case 1: external geom supplied through extended filename
+    if (m_FilenameHelper->ExtGEOMFileNameIsSet())
       {
-      otb_kwl = ReadGeometryFromImage(lFileNameOssimKeywordlist,!m_FilenameHelper->GetSkipRpcTag());
-      otbMsgDevMacro(<< "Loading internal kwl");
+      otb_kwl = ReadGeometryFromGEOMFile(m_FilenameHelper->GetExtGEOMFileName());
+      otbMsgDevMacro(<< "Loading external kwl: "<< m_FilenameHelper->GetExtGEOMFileName());
       }
+    // Case 2: attached geom (if present)
+    else if (itksys::SystemTools::FileExists(attachedGeom))
+      {
+      otb_kwl = ReadGeometryFromGEOMFile(attachedGeom);
+      otbMsgDevMacro(<< "Loading attached kwl");
+      }
+    // Case 3: find an ossimPluginProjection
+    // Case 4: find an ossimProjection
+    // Case 5: find RPC tags in TIF
     else
       {
-      otb_kwl = ReadGeometryFromGEOMFile(m_FilenameHelper->GetExtGEOMFileName());
-      otbMsgDevMacro(<< "Loading external kwl");
+      otb_kwl = ReadGeometryFromImage(lFileNameOssimKeywordlist,!m_FilenameHelper->GetSkipRpcTag());
+      otbMsgDevMacro(<< "Loading internal kwl");
       }
 
     // Don't add an empty ossim keyword list

Modules/Learning/Sampling/include/otbPersistentSamplingFilterBase.txx

@@ -265,6 +265,11 @@ PersistentSamplingFilterBase<TInputImage,TMaskImage>
       }
     this->m_InMemoryOutputs.push_back(tmpContainer);
     }
+
+  if (oSRS)
+    {
+    oSRS->Release();
+    }
 }
 
 template <class TInputImage, class TMaskImage>
@@ -734,7 +739,7 @@ PersistentSamplingFilterBase<TInputImage,TMaskImage>
     }
 
   const unsigned int nbFeatThread = std::ceil(inLayer.GetFeatureCount(true) / (float) numberOfThreads);
-  assert(nbFeatThread > 0);
+  //assert(nbFeatThread > 0);
 
   OGRFeatureDefn &layerDefn = inLayer.GetLayerDefn();
   ogr::Layer::const_iterator featIt = inLayer.begin();
@@ -807,6 +812,11 @@ PersistentSamplingFilterBase<TInputImage,TMaskImage>
       OGRFieldDefn fieldDefn(layerDefn.GetFieldDefn(k));
       outLayer.CreateField(fieldDefn);
       }
+
+    if (oSRS)
+      {
+      oSRS->Release();
+      }
     }
 
   // Add new fields

Modules/Learning/Sampling/test/otbImageSampleExtractorFilterTest.cxx

@@ -129,6 +129,7 @@ int otbImageSampleExtractorFilterUpdate(int argc, char* argv[])
   output->CreateLayer( inLayer.GetName(),
                        oSRS,
                        inLayer.GetLayerDefn().GetGeomType());
+  oSRS->Release();
   otb::ogr::Layer dstLayer = output->GetLayer(0);
   OGRFieldDefn labelField(classFieldName.c_str(),OFTString);
   dstLayer.CreateField(labelField, true);

Modules/ThirdParty/OssimPlugins/include/ossim/ossimSarSensorModel.h

@@ -238,13 +238,36 @@ public:
     * \param[in] worldPoint World point to geocode
     * \param[out] azimuthTime Estimated zero-doppler azimuth time
     * \param[out] rangeTime Estimated range time
+    * \param[out] interpSensorPos interpolated ECEF sensor position
+    * \param[out] interpSensorVel interpolated ECEF sensor velocity
     * \return True if success, false otherwise. In this case,
     * azimuthTime and rangeTime will not be modified.
     */
-   /*virtual*/ bool worldToAzimuthRangeTime(const ossimGpt& worldPt, TimeType & azimuthTime, double & rangeTime) const;
+   /*virtual*/ bool worldToAzimuthRangeTime(const ossimGpt& worldPt, TimeType & azimuthTime, double & rangeTime, ossimEcefPoint & interpSensorPos, ossimEcefVector & interpSensorVel) const;
 
-   // TODO: document me
-   /*virtual*/ void lineSampleToAzimuthRangeTime(const ossimDpt & imPt, TimeType & azimuthTime, double & rangeTime) const;
+   /**
+    * This method implement inverse sar geolocation similar to
+    * worldToLineSample, except that it also returns (y,z)
+    * coordinates, defined as follows:
+    * Let n = |sensorPos|,
+    * ps2 = scalar_product(sensorPos,worldPoint)
+    * d = distance(sensorPos,worldPoint)
+    *
+    * z = n - ps2/n
+    * y = sqrt(d*d - z*z)
+    * 
+    * sign of y is furher adapted to be inverted if sensor is left or
+    * right looking 
+    *
+    * \param[in] worldPoint World point to geocode
+    * \param[out] imPt Corresponding estimated image point
+    * \param[out] y 
+    * \param[out] z 
+    * \return True if success, false otherwise. In this case,
+    */
+   void worldToLineSampleYZ(const ossimGpt& worldPt, ossimDpt& imPt, double & y, double & z) const;
+   
+   void lineSampleToAzimuthRangeTime(const ossimDpt & imPt, TimeType & azimuthTime, double & rangeTime) const;
 
    // TODO: document me
    bool autovalidateInverseModelFromGCPs(const double & xtol = 1, const double & ytol = 1, const double azTimeTol = 500, const double &rangeTimeTo=0.0000000001) const;
@@ -419,7 +442,8 @@ protected:
    ProductType                                 theProductType; // GRD/SLC
    DurationType                                theAzimuthTimeOffset; // Offset computed
    double                                      theRangeTimeOffset; // Offset in seconds, computed
-
+   bool                                        theRightLookingFlag;
+   
    static const double C;
 
    static const unsigned int thePluginVersion; // version of the SarSensorModel plugin

Modules/ThirdParty/OssimPlugins/src/ossim/ossimSarSensorModel.cpp

@@ -12,7 +12,7 @@
  *
  * The above copyright notice and this permission notice shall be included in
  * all copies or substantial portions of the Software.
- *
+ 
  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
@@ -137,7 +137,8 @@ namespace ossimplugins
       theRangeResolution(0.),
       theBistaticCorrectionNeeded(false),
       theAzimuthTimeOffset(seconds(0)),
-      theRangeTimeOffset(0.)
+      theRangeTimeOffset(0.),
+      theRightLookingFlag(true)
       {}
 
    ossimSarSensorModel::GCPRecordType const&
@@ -220,7 +221,58 @@ namespace ossimplugins
       TimeType azimuthTime;
       double rangeTime;
 
-      const bool success = worldToAzimuthRangeTime(worldPt, azimuthTime, rangeTime);
+      ossimEcefPoint sensorPos;
+      ossimEcefVector sensorVel;
+
+      const bool success = worldToAzimuthRangeTime(worldPt, azimuthTime, rangeTime, sensorPos, sensorVel);
+
+      if(!success)
+      {
+         imPt.makeNan();
+         return;
+      }
+      // std::clog << "AzimuthTime: " << azimuthTime << "\n";
+      // std::clog << "RangeTime: " << rangeTime << "\n";
+      // std::clog << "GRD: " << isGRD() << "\n";
+
+      // Convert azimuth time to line
+      azimuthTimeToLine(azimuthTime,imPt.y);
+
+      if(isGRD())
+      {
+         // GRD case
+         double groundRange(0);
+         slantRangeToGroundRange(rangeTime*C/2,azimuthTime,groundRange);
+         // std::clog << "GroundRange: " << groundRange << "\n";
+         // std::clog << "TheRangeResolution: " << theRangeResolution << "\n";
+
+         // Eq 32 p. 31
+         // TODO: possible micro-optimization: precompute 1/theRangeResolution, and
+         // use *
+         imPt.x = groundRange/theRangeResolution;
+      }
+      else
+      {
+         // std::clog << "TheNearRangeTime: " << theNearRangeTime << "\n";
+         // std::clog << "TheRangeSamplingRate: " << theRangeSamplingRate << "\n";
+         // SLC case
+         // Eq 23 and 24 p. 28
+         imPt.x = (rangeTime - theNearRangeTime)*theRangeSamplingRate;
+      }
+   }
+
+void ossimSarSensorModel::worldToLineSampleYZ(const ossimGpt& worldPt, ossimDpt & imPt, double & y, double & z) const
+   {
+      // std::clog << "ossimSarSensorModel::worldToLineSample()\n";
+      assert(theRangeResolution>0&&"theRangeResolution is null.");
+
+      // First compute azimuth and range time
+      TimeType azimuthTime;
+      double rangeTime;
+      ossimEcefPoint sensorPos;
+      ossimEcefVector sensorVel;
+
+      const bool success = worldToAzimuthRangeTime(worldPt, azimuthTime, rangeTime,sensorPos,sensorVel);
 
       if(!success)
       {
@@ -255,9 +307,35 @@ namespace ossimplugins
          // Eq 23 and 24 p. 28
          imPt.x = (rangeTime - theNearRangeTime)*theRangeSamplingRate;
       }
+
+      // Now computes Y and Z
+      ossimEcefPoint inputPt(worldPt);
+      double NormeS = sqrt(sensorPos[0]*sensorPos[0] + sensorPos[1]*sensorPos[1] + sensorPos[2]*sensorPos[2]);  /* distance du radar */                                                                           
+      double PS2 = inputPt[0]*sensorPos[0] + inputPt[1]*sensorPos[1] + inputPt[2]*sensorPos[2];
+
+      // TODO check for small NormesS to avoid division by zero ?
+      // Should never happen ...
+      assert(NormeS>1e-6);
+      z = NormeS - PS2/NormeS;
+
+      double distance = sqrt((sensorPos[0]-inputPt[0])*(sensorPos[0]-inputPt[0]) + 
+                             (sensorPos[1]-inputPt[1])*(sensorPos[1]-inputPt[1]) + 
+                             (sensorPos[2]-inputPt[2])*(sensorPos[2]-inputPt[2]));  
+
+      y = sqrt(distance*distance - z*z);
+
+      // Check view side and change sign of Y accordingly
+      if ( (( sensorVel[0] * (sensorPos[1]* inputPt[2] - sensorPos[2]* inputPt[1]) +
+              sensorVel[1] * (sensorPos[2]* inputPt[0] - sensorPos[0]* inputPt[2]) +
+              sensorVel[2] * (sensorPos[0]* inputPt[1] - sensorPos[1]* inputPt[0])) > 0) ^ theRightLookingFlag )
+        {
+        y = -y;
+        }
    }
 
-   bool ossimSarSensorModel::worldToAzimuthRangeTime(const ossimGpt& worldPt, TimeType & azimuthTime, double & rangeTime) const
+
+
+bool ossimSarSensorModel::worldToAzimuthRangeTime(const ossimGpt& worldPt, TimeType & azimuthTime, double & rangeTime, ossimEcefPoint & interpSensorPos, ossimEcefVector & interpSensorVel) const
    {
       // std::clog << "ossimSarSensorModel::worldToAzimuthRangeTime()\n";
       // First convert lat/lon to ECEF
@@ -265,8 +343,6 @@ namespace ossimplugins
 
       // Compute zero doppler time
       TimeType interpTime;
-      ossimEcefPoint interpSensorPos;
-      ossimEcefVector interpSensorVel;
 
       const bool success = zeroDopplerLookup(inputPt,azimuthTime,interpSensorPos,interpSensorVel);
 
@@ -857,7 +933,9 @@ namespace ossimplugins
          double   estimatedRangeTime;
 
          // Estimate times
-         const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime);
+         ossimEcefPoint sensorPos;
+         ossimEcefVector sensorVel;
+         const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime,sensorPos, sensorVel);
          this->worldToLineSample(gcpIt->worldPt,estimatedImPt);
 
          const bool thisSuccess
@@ -964,8 +1042,10 @@ namespace ossimplugins
          TimeType estimatedAzimuthTime;
          double   estimatedRangeTime;
 
+         ossimEcefPoint sensorPos;
+         ossimEcefVector sensorVel;
          // Estimate times
-         const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime);
+         const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime, sensorPos, sensorVel);
 
          if(s1)
          {
@@ -985,8 +1065,11 @@ namespace ossimplugins
          TimeType estimatedAzimuthTime;
          double   estimatedRangeTime;
 
+         ossimEcefPoint sensorPos;
+         ossimEcefVector sensorVel;
+
          // Estimate times
-         const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime);
+         const bool s1 = this->worldToAzimuthRangeTime(gcpIt->worldPt,estimatedAzimuthTime,estimatedRangeTime, sensorPos, sensorVel);
 
          if(s1)
          {