Something went wrong while setting issue due date.
Incorrect indices computation in SARDEMProjection
Indices for the H band, or the SatPos band are incorrectly computed.
The problem is that these numbers are computed incrementally, in a function (GenerateOutputInformation) that is called several times.
An easy fix is the following
diff --git i/include/otbSARDEMProjectionImageFilter.txx w/include/otbSARDEMProjectionImageFilter.txx
index aadb403..9567a9a 100644
--- i/include/otbSARDEMProjectionImageFilter.txx
+++ w/include/otbSARDEMProjectionImageFilter.txx
@@ -152,22 +153,20 @@ SARDEMProjectionImageFilter< TImageIn, TImageOut >
if (m_withXYZ)
{
m_NbComponents += 3;
- m_indH += 3;
- m_indSatPos += 3;
-
inputKwl.AddKey("support_data.band.XCart", std::to_string(4));
inputKwl.AddKey("support_data.band.YCart", std::to_string(5));
inputKwl.AddKey("support_data.band.ZCart", std::to_string(6));
}
if (m_withH)
{
+ m_indH = m_NbComponents;
m_NbComponents += 1;
- m_indSatPos += 1;
inputKwl.AddKey("support_data.band.H", std::to_string(m_indH));
}
if (m_withSatPos)
{
+ m_indSatPos = m_NbComponents;
m_NbComponents += 3;
inputKwl.AddKey("support_data.band.XSatPos", std::to_string(m_indSatPos));
@@ -258,6 +257,8 @@ SARDEMProjectionImageFilter< TImageIn, TImageOut >
Point3DType satPos;
Point3DType satVel;
+ assert(m_NbComponents ==
+ 4 + (m_withXYZ ? 3 : 0) + (m_withH ? 1 : 0) + (m_withSatPos ? 3 : 0));
ImageOutPixelType pixelOutput;
pixelOutput.Reserve(m_NbComponents);
@@ -310,17 +312,20 @@ SARDEMProjectionImageFilter< TImageIn, TImageOut >
if (m_withXYZ)
{
SarSensorModelAdapter::WorldToCartesian(demGeoPoint, xyzCart);
+ assert(pixelOutput.Size() >= 6);
pixelOutput[4] = xyzCart[0];
pixelOutput[5] = xyzCart[1];
pixelOutput[6] = xyzCart[2];
}
if (m_withH)
{
+ assert(pixelOutput.Size() > m_indH);
pixelOutput[m_indH] = demGeoPoint[2];
}
if (m_withSatPos)
{
m_SarSensorModelAdapter->WorldToSatPositionAndVelocity(demGeoPoint, satPos, satVel);
+ assert(pixelOutput.Size() >= m_indSatPos + 3);
pixelOutput[m_indSatPos] = satPos[0];
pixelOutput[m_indSatPos+1] = satPos[1];
pixelOutput[m_indSatPos+2] = satPos[2];PS: it would be nice to be able to completely select what we want to produce, and for instance, ignore the first 4 bands when we don't need them.
PPS: delete nulltpr; is valid.
Activity
This wasn't that complex to implement, but after further investigations this tool isn't the one I was looking for -- as I really need the DEM/XYZ in the original S1 geometry. My current (unsorted) patchs are the following (C++ code modernization, new OTB
DEMHandler(completely untested), code correction on application boolean option handling...)diff --cc include/otbSARDEMProjectionImageFilter.txx index 49d9871,4bc8a5a..0000000 --- i/include/otbSARDEMProjectionImageFilter.txx +++ w/include/otbSARDEMProjectionImageFilter.txx @@@ -49,12 -49,12 +49,9 @@@ namespace ot */ 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_geoidEmg96(NULL) { - DEMHandlerPointerType DEMHandler = DEMHandler::Instance(); + auto& DEMHandler = DEMHandler::GetInstance(); -- m_geoidEmg96 = 0; const std::string isEmg96 = "egm96"; if (!(ConfigurationManager::GetGeoidFile().empty())) @@@ -66,24 -67,28 +63,11 @@@ std::size_t found = ConfigurationManager::GetGeoidFile().find(isEmg96); if (found!=std::string::npos) { - m_geoidEmg96 = new ossimGeoidEgm96(ossimFilename(ConfigurationManager::GetGeoidFile())); - std::cout << "\n\n########\nGeoidFile found !!!\n#########\n" << std::endl; - m_geoidEmg96 = new ossimGeoidEgm96(ossimFilename(ConfigurationManager::GetGeoidFile())); ++ m_geoidEmg96 = std::make_unique<ossimGeoidEgm96>(ossimFilename(ConfigurationManager::GetGeoidFile())); } } - else { - std::cout << "\n\n########\nGeoidFile is enmpty !!!\n#########\n" << std::endl; - } -} - -/** - * Destructor - */ -template <class TImageIn, class TImageOut> -SARDEMProjectionImageFilter<TImageIn ,TImageOut>::~SARDEMProjectionImageFilter() -{ - if (m_geoidEmg96) - { - delete m_geoidEmg96; - m_geoidEmg96 = 0; - } } - /** - * Destructor - */ - template <class TImageIn, class TImageOut> - SARDEMProjectionImageFilter<TImageIn ,TImageOut>::~SARDEMProjectionImageFilter() - { - if (m_geoidEmg96) - { - delete m_geoidEmg96; - m_geoidEmg96 = 0; - } - } - /** * Set Sar Image keyWordList */ @@@ -143,21 -148,18 +127,23 @@@ SARDEMProjectionImageFilter< TImageIn, // m_withSatPos = true => 3 Three components added : Satellite Position into cartesian // Fill KeyWordList -- inputKwl.AddKey("support_data.band.L", std::to_string(1)); -- inputKwl.AddKey("support_data.band.C", std::to_string(0)); -- inputKwl.AddKey("support_data.band.Z", std::to_string(2)); -- inputKwl.AddKey("support_data.band.Y", std::to_string(3)); -- -- m_NbComponents = 4; ++ m_NbComponents = 0; ++ if (m_withCLZY) ++ { ++ m_NbComponents += 4; ++ inputKwl.AddKey("support_data.band.L", "1"); ++ inputKwl.AddKey("support_data.band.C", "2"); ++ inputKwl.AddKey("support_data.band.Z", "3"); ++ inputKwl.AddKey("support_data.band.Y", "4"); ++ } if (m_withXYZ) { ++ m_indXYZ = m_NbComponents; m_NbComponents += 3; - m_indH += 3; - m_indSatPos += 3; - inputKwl.AddKey("support_data.band.XCart", std::to_string(4)); - inputKwl.AddKey("support_data.band.YCart", std::to_string(5)); - inputKwl.AddKey("support_data.band.ZCart", std::to_string(6)); + - inputKwl.AddKey("support_data.band.XCart", std::to_string(4)); - inputKwl.AddKey("support_data.band.YCart", std::to_string(5)); - inputKwl.AddKey("support_data.band.ZCart", std::to_string(6)); ++ inputKwl.AddKey("support_data.band.XCart", std::to_string(m_indXYZ)); ++ inputKwl.AddKey("support_data.band.YCart", std::to_string(m_indXYZ+1)); ++ inputKwl.AddKey("support_data.band.ZCart", std::to_string(m_indXYZ+2)); } if (m_withH) { @@@ -175,7 -178,7 +162,6 @@@ inputKwl.AddKey("support_data.band.ZSatPos", std::to_string(m_indSatPos+2)); } -- outputPtr->SetNumberOfComponentsPerPixel(m_NbComponents); // Create and Initilaze the SarSensorModelAdapter @@@ -258,12 -261,16 +244,12 @@@ SARDEMProjectionImageFilter< TImageIn, Point3DType satPos; Point3DType satVel; + assert(m_NbComponents == - 4 + (m_withXYZ ? 3 : 0) + (m_withH ? 1 : 0) + (m_withSatPos ? 3 : 0)); ++ (m_withCLZY ? 4 : 0) + (m_withXYZ ? 3 : 0) + (m_withH ? 1 : 0) + (m_withSatPos ? 3 : 0)); ImageOutPixelType pixelOutput; pixelOutput.Reserve(m_NbComponents); -- // Elevation from earth geoid -- double h = 0; -- ossimGpt gptPt; - - auto && demHandler = DEMHandler::Instance(); ++ auto && demHandler = DEMHandler::GetInstance(); // For each Line in MNT while ( !InIt.IsAtEnd()) @@@ -283,36 -290,37 +269,40 @@@ demGeoPoint[1] = demLatLonPoint[1]; // Get elevation from earth geoid thanks to DEMHandler or ossimgeoidEmg96 ++ double h; // Elevation from earth geoid if (m_geoidEmg96) { ++ ossimGpt gptPt; gptPt.lon = demLatLonPoint[0]; gptPt.lat = demLatLonPoint[1]; h = m_geoidEmg96->offsetFromEllipsoid(gptPt); } else { - h = DEMHandler::GetInstance().GetHeightAboveEllipsoid(demGeoPoint[0],demGeoPoint[1]); - h = demHandler->GetHeightAboveEllipsoid(demGeoPoint[0],demGeoPoint[1]); ++ h = demHandler.GetHeightAboveEllipsoid(demGeoPoint[0],demGeoPoint[1]); } // Correct the height with earth geoid demGeoPoint[2] = InIt.Get() + h; ++ // Fill the four channels C, L, Z, Y, if required ++ if (m_withCLZY) ++ { // Call the method of sarSensorModelAdapter m_SarSensorModelAdapter->WorldToLineSampleYZ(demGeoPoint, col_row, y_z); -- - // Fill the four channels for output - // Fill the four channels for output pixelOutput[0] = col_row[0]; pixelOutput[1] = col_row[1]; pixelOutput[2] = y_z[1]; pixelOutput[3] = y_z[0]; -- ++ } // Add channels if required if (m_withXYZ) { SarSensorModelAdapter::WorldToCartesian(demGeoPoint, xyzCart); - pixelOutput[4] = xyzCart[0]; - pixelOutput[5] = xyzCart[1]; - pixelOutput[6] = xyzCart[2]; - assert(pixelOutput.Size() >= 6); - pixelOutput[4] = xyzCart[0]; - pixelOutput[5] = xyzCart[1]; - pixelOutput[6] = xyzCart[2]; ++ assert(pixelOutput.Size() >= m_indXYZ + 3); ++ pixelOutput[m_indXYZ] = xyzCart[0]; ++ pixelOutput[m_indXYZ+1] = xyzCart[1]; ++ pixelOutput[m_indXYZ+2] = xyzCart[2]; } if (m_withH) { @@@ -327,7 -337,7 +319,7 @@@ } - // Set the output (without itertor because no VectorImage into iterator template) - // Set the output (without itertor because no VectorImage into iterator template) ++ // Set the output (without iterator because no VectorImage into iterator template) OutIt.Set(pixelOutput); progress.CompletedPixel(); } diff --git i/app/otbSARDEMProjection.cxx w/app/otbSARDEMProjection.cxx index d7f1f97..7c113be 100644 --- i/app/otbSARDEMProjection.cxx +++ w/app/otbSARDEMProjection.cxx @@ -80,6 +80,9 @@ private: MandatoryOff("infilemetadata"); // Empty parameters to allow more components into the projeted DEM + AddParameter(ParameterType_Bool, "withclzy", "Set C, L, Z, Y components into projection -- default: yes"); + SetParameterDescription("withclzy", "If true, then C, L, Z, Y components into projection -- default: yes."); + AddParameter(ParameterType_Bool, "withh", "Set H components into projection"); SetParameterDescription("withh", "If true, then H components into projection."); @@ -167,7 +170,6 @@ void DoExecute() override { otbAppLogWARNING(<<"Wrong Metadata into input file "); } - } preciseMetadataFile.close(); @@ -204,17 +206,29 @@ void DoExecute() override FloatImageType::Pointer inputDEM = GetParameterFloatImage("indem"); filterDEMProjection->SetInput(inputDEM); + if (IsParameterEnabled("withclzy")) + { + bool const enabled = GetParameterInt("withclzy"); + otbAppLogINFO(<<"with CLZY: " << enabled); + filterDEMProjection->SetwithCLZY(enabled); + } if (IsParameterEnabled("withh")) { - filterDEMProjection->SetwithH(true); + bool const enabled = GetParameterInt("withh"); + otbAppLogINFO(<<"with H: " << enabled); + filterDEMProjection->SetwithH(enabled); } if (IsParameterEnabled("withxyz")) { - filterDEMProjection->SetwithXYZ(true); + bool const enabled = GetParameterInt("withxyz"); + otbAppLogINFO(<<"with XYZ: " << enabled); + filterDEMProjection->SetwithXYZ(enabled); } if (IsParameterEnabled("withsatpos")) { - filterDEMProjection->SetwithSatPos(true); + bool const enabled = GetParameterInt("withsatpos"); + otbAppLogINFO(<<"with SatPos: " << enabled); + filterDEMProjection->SetwithSatPos(enabled); } filterDEMProjection->SetNoData(noData); diff --git i/include/otbSARDEMGridImageFilter.txx w/include/otbSARDEMGridImageFilter.txx index 7c429aa..d84fc08 100644 --- i/include/otbSARDEMGridImageFilter.txx +++ w/include/otbSARDEMGridImageFilter.txx @@ -183,8 +183,8 @@ namespace otb // RSTransform m_rsTransform = RSTransformType2D::New(); - m_rsTransform->SetInputImageMetadata(&(m_SarImagePtr->GetImageMetadata()) ); - m_rsTransform->SetOutputImageMetadata(&(DEMProjOnMasterPtr->GetImageMetadata())); + // m_rsTransform->SetInputImageMetadata(&(m_SarImagePtr->GetImageMetadata()) ); + // m_rsTransform->SetOutputImageMetadata(&(DEMProjOnMasterPtr->GetImageMetadata())); m_rsTransform->InstantiateTransform(); // Add GridSteps into KeyWordList diff --git i/include/otbSARDEMPolygonsAnalysisImageFilter.txx w/include/otbSARDEMPolygonsAnalysisImageFilter.txx index 53b88b5..7effd3e 100644 --- i/include/otbSARDEMPolygonsAnalysisImageFilter.txx +++ w/include/otbSARDEMPolygonsAnalysisImageFilter.txx @@ -126,7 +126,7 @@ namespace otb { // Initialize the GenericRSTransform m_RSTransform = RSTransformType2D::New(); - m_RSTransform->SetInputImageMetadata( &(m_SarImagePtr->GetImageMetadata()) ); + // m_RSTransform->SetInputImageMetadata( &(m_SarImagePtr->GetImageMetadata()) ); m_RSTransform->SetInputProjectionRef( m_SarImagePtr->GetProjectionRef() ); m_RSTransform->SetOutputProjectionRef(m_DemImagePtr->GetProjectionRef()); m_RSTransform->InstantiateTransform(); diff --git i/include/otbSARDEMProjectionImageFilter.h w/include/otbSARDEMProjectionImageFilter.h index a239707..306f073 100644 --- i/include/otbSARDEMProjectionImageFilter.h +++ w/include/otbSARDEMProjectionImageFilter.h @@ -40,9 +40,7 @@ # include "ossim/base/ossimGeoidEgm96.h" #endif - - - +#include <memory> namespace otb { @@ -60,10 +58,10 @@ namespace otb public: // Standard class typedefs - typedef SARDEMProjectionImageFilter Self; - typedef itk::ImageToImageFilter<TImageIn,TImageOut> Superclass; - typedef itk::SmartPointer<Self> Pointer; - typedef itk::SmartPointer<const Self> ConstPointer; + using Self = SARDEMProjectionImageFilter; + using Superclass = itk::ImageToImageFilter<TImageIn,TImageOut>; + using Pointer = itk::SmartPointer<Self>; + using ConstPointer = itk::SmartPointer<const Self>; // Method for creation through object factory itkNewMacro(Self); @@ -71,32 +69,32 @@ public: itkTypeMacro(SARDEMProjectionFilter,ImageToImageFilter); /** Typedef to image type */ - typedef TImageIn ImageInType; + using ImageInType = TImageIn; /** Typedef to describe the inout image pointer type. */ - typedef typename ImageInType::Pointer ImageInPointer; - typedef typename ImageInType::ConstPointer ImageInConstPointer; + using ImageInPointer = typename ImageInType::Pointer; + using ImageInConstPointer = typename ImageInType::ConstPointer; /** Typedef to describe the inout image region type. */ - typedef typename ImageInType::RegionType ImageInRegionType; + using ImageInRegionType = typename ImageInType::RegionType; /** Typedef to describe the type of pixel and point for inout image. */ - typedef typename ImageInType::PixelType ImageInPixelType; - typedef typename ImageInType::PointType ImageInPointType; + using ImageInPixelType = typename ImageInType::PixelType; + using ImageInPointType = typename ImageInType::PointType; /** Typedef to describe the image index, size types and spacing for inout image. */ - typedef typename ImageInType::IndexType ImageInIndexType; - typedef typename ImageInType::IndexValueType ImageInIndexValueType; - typedef typename ImageInType::SizeType ImageInSizeType; - typedef typename ImageInType::SizeValueType ImageInSizeValueType; - typedef typename ImageInType::SpacingType ImageInSpacingType; - typedef typename ImageInType::SpacingValueType ImageInSpacingValueType; + using ImageInIndexType = typename ImageInType::IndexType; + using ImageInIndexValueType = typename ImageInType::IndexValueType; + using ImageInSizeType = typename ImageInType::SizeType; + using ImageInSizeValueType = typename ImageInType::SizeValueType; + using ImageInSpacingType = typename ImageInType::SpacingType; + using ImageInSpacingValueType = typename ImageInType::SpacingValueType; - typedef TImageOut ImageOutType; + using ImageOutType = TImageOut; /** Typedef to describe the output image pointer type. */ - typedef typename ImageOutType::Pointer ImageOutPointer; - typedef typename ImageOutType::ConstPointer ImageOutConstPointer; + using ImageOutPointer = typename ImageOutType::Pointer; + using ImageOutConstPointer = typename ImageOutType::ConstPointer; /** Typedef to describe the output image region type. */ - typedef typename ImageOutType::RegionType ImageOutRegionType; + using ImageOutRegionType = typename ImageOutType::RegionType; /** Typedef to describe the type of pixel and point for output image. */ - typedef typename ImageOutType::PixelType ImageOutPixelType; + using ImageOutPixelType = typename ImageOutType::PixelType; typedef typename ImageOutType::PointType ImageOutPointType; /** Typedef to describe the image index, size types and spacing for output image. */ typedef typename ImageOutType::IndexType ImageOutIndexType; @@ -116,12 +114,14 @@ public: // Setter void SetSARImageKeyWorList(ImageKeywordlist sarImageKWL); itkSetMacro(NoData, int); + itkSetMacro(withCLZY, bool); itkSetMacro(withXYZ, bool); itkSetMacro(withH, bool); itkSetMacro(withSatPos, bool); // Getter itkGetMacro(NoData, int); + itkGetMacro(withCLZY, bool); itkGetMacro(withXYZ, bool); itkGetMacro(withH, bool); itkGetMacro(withSatPos, bool); @@ -131,10 +131,10 @@ protected: SARDEMProjectionImageFilter(); // Destructor - virtual ~SARDEMProjectionImageFilter() ITK_OVERRIDE; + ~SARDEMProjectionImageFilter() = default; // Print - void PrintSelf(std::ostream & os, itk::Indent indent) const ITK_OVERRIDE; + void PrintSelf(std::ostream & os, itk::Indent indent) const override; /** SARDEMProjectionImageFilter produces an image which the same size * its input image. The different between output and input images are @@ -143,12 +143,12 @@ protected: * an implementation for GenerateOutputInformation() in order to * inform the pipeline execution model. */ - virtual void GenerateOutputInformation() ITK_OVERRIDE; + void GenerateOutputInformation() override; /** SARDEMProjectionImageFilter needs a input requested region with the same size of the output requested region. * \sa ProcessObject::GenerateInputRequestedRegion() */ - virtual void GenerateInputRequestedRegion() ITK_OVERRIDE; + void GenerateInputRequestedRegion() override; /** SARDEMProjectionImageFilter can be implemented as a multithreaded filter. @@ -161,12 +161,12 @@ protected: * * \sa ImageToImageFilter::ThreadedGenerateData(), * ImageToImageFilter::GenerateData() */ - virtual void ThreadedGenerateData(const ImageOutRegionType& outputRegionForThread, itk::ThreadIdType threadId) ITK_OVERRIDE; + void ThreadedGenerateData(const ImageOutRegionType& outputRegionForThread, itk::ThreadIdType threadId) override; private: - SARDEMProjectionImageFilter(const Self&); // purposely not implemented - void operator=(const Self &); // purposely not + SARDEMProjectionImageFilter(const Self&) = delete; + void operator=(const Self &) = delete; // Instance of SarSensorModelAdapter SarSensorModelAdapter::Pointer m_SarSensorModelAdapter; @@ -175,18 +175,20 @@ protected: ImageKeywordlist m_SarImageKwl; // Value for NoData into DEM - int m_NoData; + int m_NoData = -32768; - // Boleans for more components into DEM projection - bool m_withXYZ; // Cartesian coordonates - bool m_withH; // High - bool m_withSatPos; // Topographic phase + // Flags that tell which components shall be produced + bool m_withCLZY = true; // Cl, L, Z, Y + bool m_withXYZ = false; // Cartesian coordinates + bool m_withH = false; // Height + bool m_withSatPos = false; // Topographic phase - 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 + int m_NbComponents = 0; + int m_indXYZ = -1; // index into output to set the XYZ components + int m_indH = -1; // index into output to set the H component + int m_indSatPos = -1; // index into output to set the SatPos components - ossimGeoidEgm96 * m_geoidEmg96; + std::unique_ptr<ossimGeoidEgm96> m_geoidEmg96; }; } // End namespace otbThanks for your patch. In your case, I think you need SARCartesianEstimationMean applciation with
projeted demas input (aka output of this current application SARDEMProjection)It is not complex, I know. I am just running out of time with a numerous thing in // (not only DiapOTB). To make a change, I have to validate and run some tests => take time.
mentioned in commit 77b30375
closed with merge request !40 (merged)
mentioned in commit ba4aaadc
closed with commit 77b30375
Please register or sign in to reply