Newer
Older
* Copyright (C) 2005-2021 Centre National d'Etudes Spatiales (CNES)
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
*
* 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 otbNormalCompute_h
#define otbNormalCompute_h
#include <vnl/vnl_cross.h>
#include <tuple>
#include <cmath>
namespace otb
{
/**
* Helper class to define strong/opaque types.
* The idea is to prevent easy mismatches whe passing multiples parameters of
* the sam type.
* \tparam T Encapsulated type
* \tparam Tag Tag type used to distinguish different strong types.
*/
template <typename T, typename Tag>
struct StrongType
{
StrongType() = default;
explicit StrongType(T v) : value{v} {}
explicit operator T const&() const noexcept { return value; }
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
private:
T value;
};
namespace Tags
{
struct Above{};
struct Below{};
struct Left {};
struct Right{};
} // namespace Tags
template <typename T> using Above = StrongType<T, Tags::Above>;
template <typename T> using Below = StrongType<T, Tags::Below>;
template <typename T> using Left = StrongType<T, Tags::Left>;
template <typename T> using Right = StrongType<T, Tags::Right>;
/**
* Returns the lambda that'll be used to produce the normal vector at the
* center of the neighbourhood.
*
* \tparam RealType meant to be `float` or `double`
* \param[in] x_spacing spacing in X
* \param[in] y_spacing spacing in Y
* \pre spacings shall be expressed in the same unit as the heights.
*
* This flavour supposes a constant spacing between height _pixels_.
*
* If we note:
* - \f$o\f$ the center
* - \f$a\f$ the pixel above
* - \f$b\f$ the pixel on the right
* - \f$c\f$ the pixel below
* - \f$d\f$ the pixel on the left
* - \f$X\f$ the x-axis spacing
* - \f$Y\f$ the y-axis spacing
* - \f$H_{point}\f$ the height at the given _point_
*
* We can compute the normal vector of
* - \f$\overrightarrow{ob} \wedge \overrightarrow{oa} = \begin{pmatrix}-Y.(H_b - H_o) \\ -X.(H_a - H_o) \\ X.Y\end{pmatrix}\f$
* - \f$\overrightarrow{oc} \wedge \overrightarrow{ob} = \begin{pmatrix}-Y.(H_b - H_o) \\ +X.(H_c - H_o) \\ X.Y\end{pmatrix}\f$
* - \f$\overrightarrow{od} \wedge \overrightarrow{oc} = \begin{pmatrix}+Y.(H_d - H_o) \\ +X.(H_c - H_o) \\ X.Y\end{pmatrix}\f$
* - \f$\overrightarrow{oa} \wedge \overrightarrow{od} = \begin{pmatrix}+Y.(H_d - H_o) \\ -X.(H_a - H_o) \\ X.Y\end{pmatrix}\f$
*
* Their sum is \f$\begin{pmatrix}2.Y(H_d - H_b) \\ 2.X.(H_c-H_a) \\ 4.X.Y\end{pmatrix}\f$. The
* lamdba used will return the normalized value.
*/
template <typename RealType>
inline
auto SimpleNormalMethodForRegularGrid(RealType x_spacing, RealType y_spacing)
{
auto const z = std::abs(RealType{2} * x_spacing * y_spacing);
auto const z2 = z * z;
auto lambda = [=](
Above<RealType> above, Left <RealType> left,
Right<RealType> right, Below<RealType> below)
{
auto const x = y_spacing * (RealType(left) - RealType(right));
auto const y = x_spacing * (RealType(below) - RealType(above));
auto const norm = 1.0f / std::sqrt(x*x + y*y + z2);
return std::make_tuple(x*norm, y*norm, z*norm);
};
return lambda;
}
/**
* Returns the lambda that'll be used to produce the normal vector at the
* center of the neighbourhood.
*
* \tparam PointType Models {X, Y, Z} points
*
* This flavour expects point coordinates.
*
* If we note:
* - \f$o\f$ the center
* - \f$a\f$ the pixel above
* - \f$b\f$ the pixel on the right
* - \f$c\f$ the pixel below
* - \f$d\f$ the pixel on the left
*
* We can compute the sum of all the normal vector to \f$<overrightarrow{ob}, \overrightarrow{oa}>\f$...
* - \f$\overrightarrow{ob} \wedge \overrightarrow{oa}\f$
* - \f$\overrightarrow{oc} \wedge \overrightarrow{ob}\f$
* - \f$\overrightarrow{od} \wedge \overrightarrow{oc}\f$
* - \f$\overrightarrow{oa} \wedge \overrightarrow{od}\f$
*
* Which, given cross product properties, is equal to
* \f$\overrightarrow{ca} \wedge \overrightarrow{bd}\f$
* 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, typename FloatType>
auto SimpleNormalMethodForCartesianPoints(FloatType nodata)
{
auto lambda = [=](
Above<PointType> above, Left <PointType> left,
Right<PointType> right, Below<PointType> below)
{
auto const CA = PointType(above) - PointType(below);
auto const BD = PointType(left) - PointType(right);
auto const n = vnl_cross_3d(CA, BD);
auto const norm = n.two_norm();
if (norm == 0)
{
return PointType{nodata, nodata, nodata};
}
// assert(norm > 0);
auto const inv_norm = static_cast<decltype(CA[0])>(1.0) / norm;
};
return lambda;
}
} // otb namespace
#endif // otbNormalCompute_h