Merge branch 'default_geoid_emg96' into 'master'

ENH : Default Geoid handler for emg96 file

See merge request remote_modules/diapotb!11

include/otbSARDEMProjectionImageFilter.h

@@ -28,6 +28,8 @@
 #include "otbSarSensorModelAdapter.h"
 #include "otbImageKeywordlist.h"
 
+#include "ossim/base/ossimGeoidEgm96.h"
+
 namespace otb
 {
 /** \class SARDEMProjectionImageFilter 
@@ -116,7 +118,7 @@ protected:
   SARDEMProjectionImageFilter();
 
   // Destructor
-  virtual ~SARDEMProjectionImageFilter() ITK_OVERRIDE {};
+  virtual ~SARDEMProjectionImageFilter() ITK_OVERRIDE;
 
   // Print
   void PrintSelf(std::ostream & os, itk::Indent indent) const ITK_OVERRIDE;
@@ -170,6 +172,8 @@ protected:
   int m_NbComponents;
   int m_indH; // indice into output to set the H component
   int m_indSatPos; // indice  into output to set the SatPos components
+
+  ossimGeoidEgm96 * m_geoidEmg96;
 };
 
 } // End namespace otb

include/otbSARDEMProjectionImageFilter.txx

@@ -29,6 +29,9 @@
 #include "itkProgressReporter.h"
 #include "itkNumericTraitsPointPixel.h"
 
+#include "ossim/base/ossimFilename.h"
+#include "ossim/base/ossimGpt.h"
+
 #include <cmath>
 
 namespace otb
@@ -39,16 +42,39 @@ namespace otb
   template <class TImageIn, class TImageOut> 
   SARDEMProjectionImageFilter<TImageIn ,TImageOut>::SARDEMProjectionImageFilter()
     : m_SarSensorModelAdapter(ITK_NULLPTR), m_NoData(-32768), m_withXYZ(false), m_withH(false), 
-      m_withSatPos(false), m_NbComponents(4), m_indH(4), m_indSatPos(4)
+      m_withSatPos(false), m_NbComponents(4), m_indH(4), m_indSatPos(4), m_geoidEmg96(NULL)
   {
     std::cout << "ConfigurationManager::GetGeoidFile() = " << ConfigurationManager::GetGeoidFile() << std::endl;
    
     DEMHandlerPointerType DEMHandler = DEMHandler::Instance();
+    
+    m_geoidEmg96 = 0;
+    const std::string isEmg96 =  "egm96";
 
     if (!(ConfigurationManager::GetGeoidFile().empty()))
       {
 	// DEMHandler instance to specify the geoid file
 	DEMHandler->OpenGeoidFile(ConfigurationManager::GetGeoidFile()); 
+
+	// If Geoid by default (emg96) instanciate directly a ossimGeoidEgm96 (increase performance)
+	std::size_t found = ConfigurationManager::GetGeoidFile().find(isEmg96);
+	if (found!=std::string::npos)
+	  {
+	    m_geoidEmg96 = new ossimGeoidEgm96(ossimFilename(ConfigurationManager::GetGeoidFile()));
+	  }
+      }
+  }
+
+  /** 
+   * Destructor
+   */
+  template <class TImageIn, class TImageOut> 
+  SARDEMProjectionImageFilter<TImageIn ,TImageOut>::~SARDEMProjectionImageFilter()
+  {
+    if (m_geoidEmg96)
+      {
+	delete m_geoidEmg96;
+	m_geoidEmg96 = 0;
       }
   }
 
@@ -233,6 +259,7 @@ SARDEMProjectionImageFilter< TImageIn, TImageOut >
   
   // Elevation from earth geoid
   double h = 0;
+  ossimGpt gptPt;
 
   // For each Line in MNT
   while ( !InIt.IsAtEnd())
@@ -251,8 +278,17 @@ SARDEMProjectionImageFilter< TImageIn, TImageOut >
 	      demGeoPoint[0] = demLatLonPoint[0];
 	      demGeoPoint[1] = demLatLonPoint[1];
 
-	      // Get elevation from earth geoid thanks to DEMHandler
-	      h = DEMHandler::Instance()->GetHeightAboveEllipsoid(demGeoPoint[0],demGeoPoint[1]);
+	      // Get elevation from earth geoid thanks to DEMHandler or ossimgeoidEmg96
+	      if (m_geoidEmg96)
+		{
+		  gptPt.lon =  demLatLonPoint[0];
+		  gptPt.lat =  demLatLonPoint[1];
+		  h = m_geoidEmg96->offsetFromEllipsoid(gptPt);
+		}
+	      else
+		{
+		  h = DEMHandler::Instance()->GetHeightAboveEllipsoid(demGeoPoint[0],demGeoPoint[1]); 
+		}
 
 	      // Correct the height with earth geoid
 	      demGeoPoint[2] = InIt.Get() + h;