Newer
Older
* Copyright (C) 2005-2021 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 "otbNormalCompute.h"
#include "otbRepack.h"
#include "otbPositionHelpers.h"
#include "otbApplicationWithNoData.h"
#include "otbWrapperApplicationFactory.h"
#include "otbFunctorImageFilter.h"
#include "otbDEMHandler.h"
#include "otbSarSensorModelAdapter.h"
// Elevation handler
#include "otbWrapperElevationParametersHandler.h"
#include "otbStringHelpers.h"
#include "itkFixedArray.h"
#include <array>
#include <cmath>
namespace // anonymous
{
using FloatType = float;
/**
* Returns the lambda that'll be used to produce the normal vector at the
* center of the neighbourhood.
*
* \note As a simplification, this method ignores Geiod correction on
* neighbourhood pixels. They are all supposed identical.
* \todo Handle cases with nodata => assume the data from the center; 0 if the
* center is nodata
*/
template <typename OutputPixelType, typename InputImageType>
inline
auto AppSimpleNormalMethod(
FloatType x_spacing,
FloatType y_spacing,
FloatType nodata,
using InputPixelType = typename InputImageType::PixelType;
using RealType = typename itk::PixelTraits<InputPixelType>::ValueType;
using InputV3 = vnl_vector_fixed<RealType, 3>;
using NeighborhoodIteratorType = itk::ConstNeighborhoodIterator<InputImageType>;
auto const kwl = image.GetImageKeywordlist();
auto const xband = value_or_unless(kwl, otb::bands::k_id_X, "extracting number of X band", 0);
auto const yband = value_or_unless(kwl, otb::bands::k_id_Y, "extracting number of Y band", 1);
auto const zband = value_or_unless(kwl, otb::bands::k_id_Z, "extracting number of Z band", 2);
otbLogMacro(Info, <<"Using XYZ cartesian coordinates in bands " << xband << ", " << yband << ", " << zband);
auto const XYZ = [xband, yband, zband](auto const& pixel) {
return otb::repack<InputV3>(pixel[xband], pixel[yband], pixel[zband]);
// x_spacing > 0 and y_spacing > 0 ==>
// 6 7 8
// 3 4 5
// 0 1 2
// x_spacing > 0 and y_spacing < 0 ==>
// 0 1 2
// 3 4 5
// 6 7 8
// x_spacing < 0 and y_spacing > 0 ==>
// 8 7 6
// 5 4 3
// 2 1 0
auto const above_idx = y_spacing > 0 ? 7 : 1;
auto const below_idx = y_spacing > 0 ? 1 : 7;
auto const left_idx = x_spacing > 0 ? 3 : 5;
auto const right_idx = x_spacing > 0 ? 5 : 3;
auto const center_idx = 4;
// TODO: Can we be sure this is enough to guarantee dot(res, center) > 0?
auto const above = [above_idx, XYZ](NeighborhoodIteratorType const& n) { return otb::Above<InputV3>(XYZ(n.GetPixel(above_idx)));};
auto const below = [below_idx, XYZ](NeighborhoodIteratorType const& n) { return otb::Below<InputV3>(XYZ(n.GetPixel(below_idx)));};
auto const left = [left_idx, XYZ](NeighborhoodIteratorType const& n) { return otb::Left <InputV3>(XYZ(n.GetPixel(left_idx)));};
auto const right = [right_idx, XYZ](NeighborhoodIteratorType const& n) { return otb::Right<InputV3>(XYZ(n.GetPixel(right_idx)));};
auto const center = [center_idx, XYZ](NeighborhoodIteratorType const& n) { return XYZ(n.GetPixel(center_idx));};
auto l0 = otb::SimpleNormalMethodForCartesianPoints<vnl_vector_fixed<FloatType, 3>>(nodata);
auto const nodata_res = otb::repack<OutputPixelType>(nodata, nodata, nodata);
auto lambda = [=](NeighborhoodIteratorType const& in)
{
// Compute n-> at current point
assert(in.Size() == 9);
auto const is_invalid = [nodata](auto v)
{
InputV3 const& w = v; // Remove the StrongType<> layer
return
// SARCartesianMeanEstimation will return {0, 0, 0}
w == InputV3{0., 0., 0.}
// ossim may produce nan
|| std::isnan(w[0]) || std::isnan(w[1]) || std::isnan(w[2])
#if 0
// other tools may produce nodata;
// but we only support with SARDEMProjection => #if 0
|| w[0] == nodata || w[1] == nodata || w[2] == nodata
auto const a = above(in);
auto const l = left(in);
auto const r = right(in);
auto const b = below(in);
auto const c = center(in);
// TODO: Try to produce something if any left/right/above/below pixels is
// valid => do a mean on the pixels of a side
if (is_invalid(a) || is_invalid(l) || is_invalid(r) || is_invalid(b) ||
is_invalid(c))
{
return nodata_res;
}
InputV3 n = l0(a, l, r, b);
if (InputV3{0., 0., 0.} == n) // only way to be invalid => {0, 0, 0}
// if (is_invalid(n))
{
return nodata_res;
}
assert(! std::isnan(dot_product(n, c)));
if (dot_product(n, c) < 0)
{
// TODO: we could expect this is the case for every single vector
// computed for this image => factorize
n = -n;
}
assert(dot_product(n, c) >= 0);
std::cout << "-----@" << in.GetIndex() << " from " <<
InputV3(above(in)) << " ; " << InputV3(left(in)) << " ; " << InputV3(right(in)) << " ; " << InputV3(below(in))
<< " --> " << n
<<std::endl;
auto const res = otb::repack<OutputPixelType>(n);
return res;
};
return lambda;
}
} // anonymous namespace.
namespace otb
{
namespace Wrapper
{
/**
* This application takes cares of building the normal vector from a DEM image.
*
* \author Luc Hermitte
* \copyright CS Group
* \ingroup App???
* \ingroup SARCalibrationExtended
*/
class ExtractNormalVector : public ApplicationWithNoData
{
public:
/** Standard class typedefs. */
using Self = ExtractNormalVector;
// using InputPixelType = itk::FixedArray<FloatType, 3>;
/** Pixel Type of input image: VLV.
* Some XYZ producers may provide mode than 3 bands... */
using InputPixelType = itk::VariableLengthVector<FloatType>;
// using InputImageType = otb::Image<InputPixelType, 2>;
using InputImageType = otb::VectorImage<FloatType, 2>;
using OutputPixelType = itk::VariableLengthVector<FloatType>;
// using InputImageType = otb::Image<OutputPixelType, 2>;
// using WriterType = otb::ImageFileWriter<InputImageType>;
/** Standard macro */
itkNewMacro(Self);
itkTypeMacro(ExtractNormalVector, unused);
private:
void DoInit() override
{
SetName("ExtractNormalVector");
SetDescription("This application produces normal vectors from a DEM image");
// Documentation
SetDocLongDescription("This application produces normal vectors from a DEM image" );
SetDocLimitations("None");
SetDocAuthors("Luc Hermitte (CS Group)");
AddDocTag(Tags::Calibration);
AddParameter(ParameterType_InputImage, "xyz", "filename for the file containing XYZ cartesian coordinates aligned of their file geometry");
SetParameterDescription("xyz", "filename for the file containing XYZ cartesian coordinates aligned of their file geometry");
AddParameter(ParameterType_OutputImage, "out", "Image of normal vectors");
SetParameterDescription("out", "Image of normal vectors");
AddRAMParameter();
SetDocExampleParameterValue("xyz", "xyz.tif");
SetDocExampleParameterValue("out", "normal.tif");
SetOfficialDocLink();
}
void DoUpdateParameters() override
{
}
void DoExecute() override
{
auto* inputImage = GetParameterImage<InputImageType>("xyz");
assert(inputImage);
auto const nodata = GetParameterFloat("nodata");
// std::cout << "nodata => " << nodata << std::endl;
AddNodataInMetadata(nodata);
auto const spacing = inputImage->GetSignedSpacing();
FloatType const x_spacing = spacing[0];
FloatType const y_spacing = spacing[1];
auto const lambda = AppSimpleNormalMethod<OutputPixelType>(x_spacing, y_spacing, nodata, *inputImage);
auto filter = otb::NewFunctorFilter(lambda, 3, {{1,1}});
filter->SetInput(inputImage);
SetParameterOutputImage("out", filter->GetOutput());
// SetParameterOutputImagePixelType("out",ImagePixelType_float);
// Call register pipeline to allow streaming and garbage collection
RegisterPipeline();
}
};
} //end namespace Wrapper
} //end namespace otb
OTB_APPLICATION_EXPORT(otb::Wrapper::ExtractNormalVector)