1 /*
2 GNU Octave level-set package.
3 Copyright (C) 2015 Daniel Kraft <d@domob.eu>
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 /* Implement an approximation to |grad phi| using an appropriate
20 upwind scheme. The scheme implemented is basically (6.17)
21 from Sethian's book "Level Set Methods and Fast Marching Methods",
22 Cambridge University Press, second edition, 1999.
23
24 The discretisation could probably be discretised also in Octave code
25 with proper vector operations, but this would to harder-to-read code.
26 The direct, point-wise C++ formulation is closer to the formulas
27 that are derived in the book. */
28
29 #include <octave/oct.h>
30 #include <octave/lo-ieee.h>
31
32 #include <algorithm>
33 #include <cmath>
34 #include <cstddef>
35 #include <stdexcept>
36
37 /** Type used to index n-dimensional Octave arrays. */
38 typedef Array<octave_idx_type> GridIndex;
39 /** Type used to denote number of dimensions. */
40 typedef size_t dimensionT;
41
42 /**
43 * Calculate the gradient norm in a given grid point.
44 * @param idx The current point's index.
45 * @param phi The function whose gradient is requested.
46 * @param F The sign gives the upwind direction.
47 * @param size The total domain size.
48 * @param D The number of dimensions.
49 * @return The approximate gradient norm in the point, not yet divided by h.
50 */
51 static double
gradientNormInPoint(const GridIndex & idx,const NDArray & phi,double F,const dim_vector & size,dimensionT D)52 gradientNormInPoint (const GridIndex& idx, const NDArray& phi, double F,
53 const dim_vector& size, dimensionT D)
54 {
55 double res = 0.0;
56 for (dimensionT d = 0; d < D; ++d)
57 {
58 GridIndex neighbour(idx);
59 double forward = 0.0;
60 double backward = 0.0;
61
62 if (idx(d) > 0)
63 {
64 neighbour(d) = idx(d) - 1;
65 backward = phi(idx) - phi(neighbour);
66 }
67 if (idx(d) + 1 < size(d))
68 {
69 neighbour(d) = idx(d) + 1;
70 forward = phi(neighbour) - phi(idx);
71 }
72
73 if (F < 0.0)
74 std::swap (forward, backward);
75
76 if (backward > 0.0)
77 res += std::pow (backward, 2);
78 if (forward < 0.0)
79 res += std::pow (forward, 2);
80 }
81
82 return std::sqrt (res);
83 }
84
85 /**
86 * Recursively (since the number of dimensions and thus required loops
87 * is not known) iterate the whole array, calling gradientNormInPoint
88 * at each grid point.
89 * @param idx The current index.
90 * @param depth The current depth in the recursion.
91 * @param phi The function whose gradient is requested.
92 * @param F The sign defines the upwind direction.
93 * @param size Size of the array.
94 * @param D Number of dimensions.
95 * @param grad Return the gradient norm here.
96 */
97 static void
iterate(GridIndex & idx,dimensionT depth,const NDArray & phi,const NDArray & F,const dim_vector & size,dimensionT D,NDArray & grad)98 iterate (GridIndex& idx, dimensionT depth, const NDArray& phi, const NDArray& F,
99 const dim_vector& size, dimensionT D, NDArray& grad)
100 {
101 if (depth == D)
102 {
103 grad(idx) = gradientNormInPoint (idx, phi, F(idx), size, D);
104 return;
105 }
106
107 for (idx(depth) = 0; idx(depth) < size(depth); ++idx(depth))
108 iterate (idx, depth + 1, phi, F, size, D, grad);
109 }
110
111 /* Octave routine interfacing to fast marching code. */
112 DEFUN_DLD (__levelset_upwindGrad, args, nargout,
113 "GNORM = upwindGrad (PHI, F, H)\n\n"
114 "Internal routine to calculate the upwind-discretised gradient norm\n"
115 "of PHI. F is the speed field, and only the sign of F is used to decide\n"
116 "about the upwind direction. H is the grid size as usual.\n\n"
117 "Returned is an array GNORM of the same size as PHI, with |grad phi| in\n"
118 "the entries.\n")
119 {
120 try
121 {
122 if (args.length () != 3 || nargout != 1)
123 throw std::invalid_argument ("invalid argument counts");
124
125 const NDArray phi = args(0).array_value ();
126 const NDArray F = args(1).array_value ();
127 const double h = args(2).double_value ();
128
129 const dim_vector size = phi.dims ();
130 const dimensionT D = size.length ();
131
132 if (size != F.dims ())
133 throw std::invalid_argument ("argument dimensions mismatch");
134
135 NDArray grad(size);
136 grad.fill (octave_NA);
137
138 GridIndex idx(dim_vector (D, 1));
139 iterate (idx, 0, phi, F, size, D, grad);
140 grad /= h;
141
142 octave_value_list res;
143 res(0) = grad;
144
145 return res;
146 }
147 catch (const std::exception& exc)
148 {
149 error (exc.what ());
150 return octave_value_list ();
151 }
152 }
153