Commit bc7aefa9 authored by Luc Hermitte's avatar Luc Hermitte
Browse files

Merge branch 'master' into 3-lot-of-variations-in-some-places-where-topography-varies-a-lot

parents a2bb4f23 897660f5
#
# Copyright (C) 2005-2020 Centre National d'Etudes Spatiales (CNES)
# Copyright (C) 2005-2022 Centre National d'Etudes Spatiales (CNES)
#
# This file is part of S1Tiling remote module for Orfeo Toolbox
#
......@@ -32,3 +32,9 @@ otb_create_application(
NAME SARComputeLocalIncidenceAngle
SOURCES otbSARComputeLocalIncidenceAngle.cxx
LINK_LIBRARIES ${otb-module} ${${otb-module}_LIBRARIES})
otb_create_application(
NAME SARCartesianMeanEstimation2
SOURCES otbSARCartesianMeanEstimation2.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES})
......@@ -23,7 +23,7 @@
#include "otbRepack.h"
#include "otbPositionHelpers.h"
#include "otbWrapperApplication.h"
#include "otbApplicationWithNoData.h"
#include "otbWrapperApplicationFactory.h"
#include "otbFunctorImageFilter.h"
#include "otbDEMHandler.h"
......@@ -322,7 +322,7 @@ namespace Wrapper
* \ingroup App???
* \ingroup SARCalibrationExtended
*/
class ExtractNormalVector : public Application
class ExtractNormalVector : public ApplicationWithNoData
{
public:
/** Standard class typedefs. */
......@@ -367,10 +367,7 @@ private:
AddParameter(ParameterType_Bool, "mean", "Compute normals from means of cartesian positions on each side");
SetParameterDescription("mean", "Compute normals from means of cartesian positions on each side. By default, use only the positions of the 4 immediate side pixels (right, left, above, below)");
AddParameter(ParameterType_Float, "nodata", "NoData value");
SetParameterDescription("nodata", "DEM empty cells are filled with this value (optional -32768 by default)");
SetDefaultParameterFloat("nodata", -32768);
MandatoryOff("nodata");
DoInit_NoData();
AddRAMParameter();
......@@ -413,6 +410,7 @@ private:
otbLogMacro(Info, <<"Computing normals from "
<< (does_compute_mean_on_sides
? "all pixels on each side" : "a single pixel from each side"));
AddNodataInMetadata(nodata);
auto const spacing = inputImage->GetSignedSpacing();
FloatType const x_spacing = spacing[0];
......
/*
* Copyright (C) 2005-2022 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 "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbSARDEMPolygonsAnalysisImageFilter2.h"
#include "otbSARCartesianMeanFunctor2.h"
#include "otbWrapperOutputImageParameter.h"
#include <string>
#include <sstream>
namespace otb
{
namespace Wrapper
{
class SARCartesianMeanEstimation2 : public Application
{
public:
using Self = SARCartesianMeanEstimation2;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
itkTypeMacro(SARCartesianMeanEstimation2, otb::Wrapper::Application);
// Pixels
using ImageInPixelType = typename FloatVectorImageType::PixelType;
using ImageOutPixelType = typename FloatImageType::PixelType;
// Function
using PolygonsFunctorType = otb::Function::SARPolygonsFunctor<ImageInPixelType, ImageOutPixelType>;
using CartesianMeanFunctorType = otb::Function::SARCartesianMeanFunctor<ImageInPixelType, ImageInPixelType>;
// Filters
using FilterType = otb::SARDEMPolygonsAnalysisImageFilter < FloatVectorImageType, FloatVectorImageType, FloatImageType, ComplexFloatImageType, CartesianMeanFunctorType >;
private:
void DoInit() override
{
SetName("SARCartesianMeanEstimation2");
SetDescription("SAR Cartesian mean estimation thanks to the associated DEM.");
SetDocLongDescription("This application estimates a simulated cartesian mean image thanks to a DEM file.");
//Optional descriptors
SetDocLimitations("Only Sentinel 1 (IW and StripMap mode) and Cosmo products are supported for now.");
SetDocAuthors("OTB-Team");
SetDocSeeAlso(" ");
AddDocTag(Tags::SAR);
AddDocTag("DiapOTB");
//Parameter declarations
AddParameter(ParameterType_InputImage, "indemproj", "Input vector of DEM projected into SAR geometry");
SetParameterDescription("indemproj", "Input vector image for cartesian mean estimation.");
AddParameter(ParameterType_InputImage, "indem", "Input DEM");
SetParameterDescription("indem", "DEM to extract DEM geometry.");
AddParameter(ParameterType_InputImage, "insar", "Input SAR image");
SetParameterDescription("insar", "SAR Image to extract SAR geometry.");
AddParameter(ParameterType_Int, "indirectiondemc", "Range direction for DEM scan");
SetParameterDescription("indirectiondemc", "Range direction for DEM scan.");
SetDefaultParameterInt("indirectiondemc", 1);
MandatoryOff("indirectiondemc");
AddParameter(ParameterType_Int, "indirectiondeml", "Azimut direction for DEM scan");
SetParameterDescription("indirectiondeml", "Azimut direction for DEM scan.");
SetDefaultParameterInt("indirectiondeml", 1);
MandatoryOff("indirectiondeml");
AddParameter(ParameterType_Int, "mlran", "Averaging on distance (output geometry)");
SetParameterDescription("mlran","Averaging on distance (output geometry)");
SetDefaultParameterInt("mlran", 3);
SetMinimumParameterIntValue("mlran", 1);
MandatoryOff("mlran");
AddParameter(ParameterType_Int, "mlazi", "Averaging on azimut (output geometry)");
SetParameterDescription("mlazi","Averaging on azimut (output geometry)");
SetDefaultParameterInt("mlazi", 3);
SetMinimumParameterIntValue("mlazi", 1);
MandatoryOff("mlazi");
AddParameter(ParameterType_OutputImage, "out", "Output cartesian (mean) Image for DEM Projection");
SetParameterDescription("out","Output cartesian (mean) Image for DEM Projection."
" The first 3 bands contain X, Y and Z cartesian coordinates of the point in original SAR geometry."
" The last band is a mask band: 0: XYZ coordinates cannot be computed; 1: coordinates are nominally computed; 2: coordinates are estimated, but point is estimated in shadow"
);
AddRAMParameter();
SetDocExampleParameterValue("indemproj","CLZY_S21E055.tiff");
SetDocExampleParameterValue("indem","S21E055.hgt");
SetDocExampleParameterValue("insar","s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002_SLC.tiff");
SetDocExampleParameterValue("out","s1a-s4-simu-cartMean.tiff");
}
void DoUpdateParameters() override
{
// Nothing to do here : all parameters are independent
}
void DoExecute() override
{
// Get numeric parameters : ML factors and directions for DEM scan
int mlRan = GetParameterInt("mlran");
int mlAzi = GetParameterInt("mlazi");
int DEMScanDirectionC = GetParameterInt("indirectiondemc");
int DEMScanDirectionL = GetParameterInt("indirectiondeml");
otbAppLogINFO(<<"ML Range : "<<mlRan);
otbAppLogINFO(<<"ML Azimut : "<<mlAzi);
otbAppLogINFO(<<"Direction (DEM scan) in range : "<<DEMScanDirectionC);
otbAppLogINFO(<<"Direction (DEM scan) in azimut : "<<DEMScanDirectionL);
// Start the first pipelines (Read SAR and DEM image metedata)
ComplexFloatImageType::Pointer SARPtr = GetParameterComplexFloatImage("insar");
FloatImageType::Pointer DEMPtr = GetParameterFloatImage("indem");
// CartesianMeanEstimation Filter (use SARDEMPolygonsAnalysisImageFilter with SARCartesianMeanEstimationFunctor
// to estimate cartesian coordonates mean image)
FilterType::Pointer filterCartesianMeanEstimation = FilterType::New();
m_Ref.push_back(filterCartesianMeanEstimation.GetPointer());
filterCartesianMeanEstimation->SetSARImageKeyWorList(SARPtr->GetImageKeywordlist());
filterCartesianMeanEstimation->SetSARImagePtr(SARPtr);
filterCartesianMeanEstimation->SetDEMImagePtr(DEMPtr);
filterCartesianMeanEstimation->SetDEMInformation(0, DEMScanDirectionC, DEMScanDirectionL);
filterCartesianMeanEstimation->initializeMarginAndRSTransform();
// Function Type : SARCartesianMeanEstimationFunctor
int nbThreads = filterCartesianMeanEstimation->GetNumberOfThreads();
CartesianMeanFunctorType::Pointer functor = CartesianMeanFunctorType::New(nbThreads, mlRan, mlAzi,
SARPtr->GetOrigin()[0]-0.5);
filterCartesianMeanEstimation->SetFunctorPtr(functor);
// Define the main pipeline (controlled with extended FileName)
std::string origin_FileName = GetParameterString("out");
// Check if FileName is extended (with the ? caracter)
// If not extended then override the FileName
if (origin_FileName.find("?") == std::string::npos && !origin_FileName.empty())
{
std::string extendedFileName = origin_FileName;
// Get the ram value (in MB)
int ram = GetParameterInt("ram");
// Define with the ram value, the number of lines for the streaming
int nbColSAR = SARPtr->GetLargestPossibleRegion().GetSize()[0];
int nbLinesSAR = SARPtr->GetLargestPossibleRegion().GetSize()[1];
// To determine the number of lines :
// nbColSAR * nbLinesStreaming * sizeof(OutputPixel) = RamValue/2 (/2 to be sure)
// OutputPixel are float => sizeof(OutputPixel) = 4 (Bytes)
long int ram_In_KBytes = (ram/2) * 1024;
long int nbLinesStreaming = (ram_In_KBytes / (nbColSAR)) * (1024/4);
// Check the value of nbLinesStreaming
int nbLinesStreamingMax = 10000;
if (nbLinesStreamingMax > nbLinesSAR)
{
nbLinesStreamingMax = nbLinesSAR;
}
if (nbLinesStreaming <= 0 || nbLinesStreaming > nbLinesStreamingMax)
{
nbLinesStreaming = nbLinesStreamingMax;
}
// Construct the extendedPart
std::ostringstream os;
os << "?&streaming:type=stripped&streaming:sizevalue=" << nbLinesStreaming;
// Add the extendedPart
std::string extendedPart = os.str();
extendedFileName.append(extendedPart);
// Set the new FileName with extended options
SetParameterString("out", extendedFileName);
}
// Execute the main Pipeline
filterCartesianMeanEstimation->SetInput(GetParameterImage("indemproj"));
SetParameterOutputImage("out", filterCartesianMeanEstimation->GetOutput());
}
// Vector for filters
std::vector<itk::ProcessObject::Pointer> m_Ref;
};
}
}
OTB_APPLICATION_EXPORT(otb::Wrapper::SARCartesianMeanEstimation2)
......@@ -19,7 +19,7 @@
*/
#include "otbSARComputeLocalIncidenceAngle.h"
#include "otbWrapperApplication.h"
#include "otbApplicationWithNoData.h"
#include "otbWrapperApplicationFactory.h"
// #define DOUBLES_EVERYWHERE
......@@ -47,7 +47,7 @@ namespace Wrapper
* \ingroup App???
* \ingroup SARCalibrationExtended
*/
class SARComputeLocalIncidenceAngle : public Application
class SARComputeLocalIncidenceAngle : public ApplicationWithNoData
{
public:
/** Standard class typedefs. */
......@@ -109,6 +109,7 @@ private:
#endif
SetParameterDescription(k_out_sin, "Output image with abs(sin LIA)");
DoInit_NoData();
AddParameter(ParameterType_Float, "nodata", "NoData value");
SetParameterDescription("nodata", "DEM empty cells are filled with this value (optional -32768 by default)");
SetDefaultParameterFloat("nodata", -32768);
......@@ -139,6 +140,8 @@ private:
// filter->ShallComputeLIASign(true);
auto const nodata = GetParameterFloat("nodata");
AddNodataInMetadata(nodata, k_out_lia);
AddNodataInMetadata(nodata, k_out_sin);
filter->SetNoData(nodata);
SetParameterOutputImage(k_out_sin, filter->GetSinOutputImage());
......
/*
* Copyright (C) 2005-2022 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.
*/
#ifndef equal_range_interval_h
#define equal_range_interval_h
#include <algorithm> // std::lower_bound, std::upper_bound
#include <cassert>
#include <functional> // std::less
#include <utility> // std::pair
namespace otb
{
/**
* Adaptation of `std::equal_range` for finding interval bounds.
* \tparam RandomAccessIterator Iterator type that should model `std::random_access_iterator`.
* \tparam T Type of the interval bounds
* \tparam Comp Comparison function, defaults to `std::less<>`
* \param[in] first Start of the range where bounds are searched
* \param[in] last Past-end of the range where bounds are searched
* \param[in] lo Low bound of the interval
* \param[in] hi High bound of the interval
* \param[in] cmp Comparison predicate
*
* \return a pair of iterators `[it_low, it_high[` such as:
* - `∀ it ∈ [first, it_low[ => cmp(*it, lo)`
* - `∀ it ∈ [it_low, last[ => !cmp(*it, lo)`
* - `∀ it ∈ [first, it_high[ => !cmp(hi, *it)`
* - `∀ it ∈ [it_high, last[ => cmp(hi, *it)`
*
* Somehow, `it_low` is computed with `std::lower_bound`, and `it_hi` is
* computed with `std::upper_bound`.
*
* \pre `cmp(lo, hi)`
* \pre `[first, last[` shall be partitioned according to `cmp` predicate
* around `lo` and around `high`. If the range is sorted according to `cmp`,
* then this also implies it's partitionned.
* \throw None (it depends on `Cmp`...)
*/
template <typename RandomAccessIterator, typename T,
typename Comp = std::less<typename RandomAccessIterator::value_type>>
inline
std::pair<RandomAccessIterator, RandomAccessIterator>
equal_range_interval(
RandomAccessIterator first, RandomAccessIterator last,
T lo, T hi,
Comp cmp = Comp{})
{
// assert(std::is_sorted(first, last, cmp));
auto len = std::distance(first, last);
while (len > 0) {
auto half = len / 2;
auto middle = first + half;
if (cmp(*middle, lo)) {
first = middle + 1;
len -= half + 1;
} else if (cmp(hi, *middle)) {
len = half;
} else {
auto left = std::lower_bound(first, middle, lo, cmp);
auto right = std::upper_bound(middle+1, first+len, hi, cmp);
return {left, right};
}
}
return {first, first};
}
} // namespace otb
#endif // equal_range_interval_h
/*
* Copyright (C) 2005-2022 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.
*/
#ifndef otbApplicationWithNoData_h
#define otbApplicationWithNoData_h
#include "otbWrapperApplication.h"
namespace otb
{
namespace Wrapper
{
/** Helper parent class for application that need to inject `&nodata=` in
* metadata.
*
* This class provides two services:
* - `DoInit_NoData()` meant to be called from `DoInit()` specialisation.
* - `AddNodataInMetadata()` meant to be called from `DoExecute()`
* specialisation.
*
* \author Luc Hermitte
* \copyright CS Group
* \ingroup App???
* \ingroup SARCalibrationExtended
*/
class ApplicationWithNoData : public Application
{
protected:
/**
* Registers a `-nodata` parameter.
*/
void DoInit_NoData()
{
AddParameter(ParameterType_Float, "nodata", "NoData value");
SetParameterDescription("nodata", "DEM empty cells are filled with this value (optional -32768 by default)");
SetDefaultParameterFloat("nodata", -32768);
MandatoryOff("nodata");
}
/**
* Appends `&nodata={value}` at the end of the extended filename to force
* nodata to be added in metadata.
*/
template <typename FloatType>
void AddNodataInMetadata(FloatType nodata, std::string const& out_param = "out")
{
std::string origin_FileName = GetParameterString(out_param);
std::ostringstream oss;
oss << origin_FileName;
// Check if FileName is extended (with the ? caracter)
// If not extended then override the FileName
auto const extension_start = origin_FileName.find('?');
if (extension_start == std::string::npos && !origin_FileName.empty())
oss << '?';
else
{
auto const nodata_start = origin_FileName.find("&nodata=");
if (nodata_start > extension_start)
{
// Let's trust the end-user. Even if the value doesn't match
otbLogMacro(Warning, <<"Trusting the nodata value required in extended filename. User specified " << nodata << " value won't be propagated in image metadata");
return ;
}
}
oss << "&nodata="<<nodata;
// Set the new FileName with extended options
SetParameterString(out_param, oss.str());
}
};
} //end namespace Wrapper
} //end namespace otb
#endif // otbApplicationWithNoData_h
......@@ -134,9 +134,9 @@ auto SimpleNormalMethodForRegularGrid(RealType x_spacing, RealType y_spacing)
* The lamdba used will return the normalized value, or a triplet of `nodata`
* if any of the two vertical or horizontal vector is null.
*/
template <typename PointType>
template <typename PointType, typename FloatType>
inline
auto SimpleNormalMethodForCartesianPoints(auto nodata)
auto SimpleNormalMethodForCartesianPoints(FloatType nodata)
{
auto lambda = [=](
Above<PointType> above, Left <PointType> left,
......
......@@ -23,6 +23,7 @@
#include "otbRepack.h"
#include "otbSarSensorModelAdapter.h"
#include <itkContinuousIndex.h>
namespace otb
{
......@@ -36,6 +37,19 @@ auto TransformIndexToPhysicalPoint(ImageType & im, typename ImageType::IndexType
return p;
}
template <typename IndexRep, typename ImageType>
inline
auto TransformPhysicalPointToContinuousIndex(ImageType & im, typename ImageType::PointType const& p)
{
// There is no ImageType::ContinuousIndex in otb images
// => work around it
// static_assert(std::is_same<typename ImageType::IndexValueType, double>::value, "");
static_assert(ImageType::ImageDimension == 2, "");
itk::ContinuousIndex<IndexRep, ImageType::ImageDimension> index;
im.TransformPhysicalPointToContinuousIndex(p, index);
return index;
}
template <typename Point3DType>
inline
auto WorldToCartesian(Point3DType const& lat_lon_hgt)
......
/*
* Copyright (C) 2005-2022 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.
*/
#ifndef otbSARCartesianMeanFunctor2_h
#define otbSARCartesianMeanFunctor2_h
#include "otbSARPolygonsFunctor2.h"
#include "itkObject.h"
#include "itkObjectFactory.h"
#include <cassert>
namespace otb
{
namespace Function
{
/** \class SARCartesianMeanFunctor
* \brief Catesian mean estimation.
*
* This functor is an implementation of SARPolygonsFunctor. It estimates the Cartesian coordonates mean image.
*
* The characteristics are the following :
* _ Output Geometry : ML (defined by ML factors)
* _ Seven required components : L, C, Y, Z, XCart, YCart and ZCart
* _ Four estimated components : Mean of XCart, YCart, ZCart and an isData mask
* _ Optionnal image : No
*
* \ingroup DiapOTBModule
*/
template <class TInputPixel, class TOutputPixel>
class SARCartesianMeanFunctor : public SARPolygonsFunctor<TInputPixel, TOutputPixel>
{
public:
/** Standard class typedefs */
using Self = SARCartesianMeanFunctor;
using Superclass = SARPolygonsFunctor<TInputPixel, TOutputPixel>;
using Pointer = itk::SmartPointer<Self>;
using ConstPointer = itk::SmartPointer<const Self>;
/** Runtime information */
itkTypeMacro(SARCartesianMeanFunctor, itk::Object);
/**
* Method SetRequiredComponentsToInd (override)
*/
void SetRequiredComponentsToInd(std::map<std::string, int> const& mapForInd) final
{
std::map<std::string, int>::const_iterator it = mapForInd.begin();
while(it!=mapForInd.end())
{
if (it->first == "L")
{
m_indL = it->second;
}
else if (it->first == "C")
{
m_indC = it->second;
}
else if (it->first == "Y")
{
m_indY = it->second;
}
else if (it->first == "Z")
{
m_indZ = it->second;
}
else if (it->first == "XCart")
{
m_indXCart = it->second;
}
else if (it->first == "YCart")
{
m_indYCart = it->second;