diff --git a/include/otbSARDEMProjectionImageFilter.h b/include/otbSARDEMProjectionImageFilter.h index d9f11bc615cc9939eb722813e204cfbc25258f13..cd9f0f78e69a783232e622c0d89d4f71dca3151f 100644 --- a/include/otbSARDEMProjectionImageFilter.h +++ b/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 diff --git a/include/otbSARDEMProjectionImageFilter.txx b/include/otbSARDEMProjectionImageFilter.txx index 1d86e34dc042b8940a161ca546bd140f2d4d6f27..da9ba196b899bef948fae568b2b09453f4ead39c 100644 --- a/include/otbSARDEMProjectionImageFilter.txx +++ b/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;