1## Copyright (C) 2014 Daniel Kraft <d@domob.eu> 2## GNU Octave level-set package. 3## 4## This program is free software: you can redistribute it and/or modify 5## it under the terms of the GNU General Public License as published by 6## the Free Software Foundation, either version 3 of the License, or 7## (at your option) any later version. 8## 9## This program is distributed in the hope that it will be useful, 10## but WITHOUT ANY WARRANTY; without even the implied warranty of 11## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12## GNU General Public License for more details. 13## 14## You should have received a copy of the GNU General Public License 15## along with this program. If not, see <http://www.gnu.org/licenses/>. 16 17## -*- texinfo -*- 18## @deftypefn {Function File} {@var{d} =} ls_solve_stationary (@var{phi}, @var{f}, @var{h} = 1) 19## @deftypefnx {Function File} {@var{d} =} ls_solve_stationary (@var{phi}, @var{f}, @var{h} = 1, @var{nb}) 20## 21## Solve a generalised Eikonal equation with speeds of arbitrary 22## signs. The equation solved is 23## @tex 24## \begin{equation*} 25## f \left| \nabla d \right| = 1 26## \end{equation*} 27## @end tex 28## @ifnottex 29## 30## @example 31## f | grad d | = 1 32## @end example 33## 34## @end ifnottex 35## with d = 0 on the boundary. The domain is described by the level-set 36## function @var{phi} on a rectangular grid. @var{f} should contain the values 37## of the speed field on the grid points. @var{h} can be given as the 38## grid spacing. In the second form, where the optional @var{nb} is given, 39## it is used to initialise the narrow band with a manual calculation. 40## By default, the result of @code{ls_init_narrowband} is used. Values which 41## are not fixed by the narrow band should be set to @code{NA}. 42## 43## Note that in comparison to @code{fastmarching}, the speed need not be 44## positive. It is the reciprocal of @var{f} in @code{fastmarching}. 45## 46## This is a preparation step, and afterwards, the evolved geometry according 47## to the level-set equation 48## @tex 49## \begin{equation*} 50## \phi_t + f \left| \nabla \phi \right| = 0 51## \end{equation*} 52## @end tex 53## @ifnottex 54## 55## @example 56## d/dt phi + f | grad phi | = 0 57## @end example 58## 59## @end ifnottex 60## can be extracted from @var{d} at arbitrary positive times using 61## the supplemental function @code{ls_extract_solution}. 62## 63## At points where @var{f} is exactly zero, the output will be set 64## to @code{NA}. This case is handled separately in @code{ls_extract_solution}. 65## In the narrow band, the returned distances may actually be negative 66## even for positive @var{f} and vice-versa. This helps to avoid 67## unnecessary errors introduced into the level-set function due to 68## a finite grid-size when the time step is chosen small. 69## 70## @seealso{ls_extract_solution, ls_signed_distance, ls_nb_from_geom} 71## @end deftypefn 72 73function d = ls_solve_stationary (phi, f, h = 1, nb) 74 if (nargin () < 2 || nargin () > 4) 75 print_usage (); 76 endif 77 78 sz = size (phi); 79 if (~all (sz == size (f))) 80 error ("PHI and F must be of the same size"); 81 endif 82 83 if (nargin () < 4) 84 nb = ls_init_narrowband (phi, h); 85 elseif (~all (sz == size (nb))) 86 error ("PHI and NB must be of the same size"); 87 endif 88 89 % Find regions of signs of f. Below for the main calculations, 90 % we will always ignore f = 0 regions and will mostly work with 91 % the absolute value of f, fixing the sign of d at the end. 92 fp = (f > 0); 93 fz = (f == 0); 94 fn = (f < 0); 95 fabs = abs (f); 96 97 % Find regions of the domain. 98 inner = (phi < 0); 99 outer = (phi > 0); 100 101 % Initialise narrow-band distances. Scale them (roughly) according 102 % to our f. 103 d0 = nb; 104 d0(~fz) ./= fabs(~fz); 105 106 % Handle points where the front moves away (depending on the signs of 107 % phi and f). In theory, these distances should be zero according to 108 % the definition in the theoretical paper. Set them to Inf or -Inf 109 % in the numerics for points deeply inside, and leave the narrow-band 110 % values for points near the boundary. Since fastmarching does not 111 % allow -Inf, choose the signs such that the points will be passed with 112 % positive Inf to fastmarching. The sign will be reverted later 113 % on to get the value matching the distance semantics (-Inf inside with 114 % front moving away and Inf outside). 115 % For phi = 0, nothing special needs to be done as d0 is set to zero there 116 % already by internal_init_narrowband. 117 d0(isna (d0) & inner & fp) = Inf; 118 d0(isna (d0) & outer & fn) = -Inf; 119 assert (all (isna (d0(:)) | ~isnan(d0(:)))); 120 121 % Solve for both parts. Invert the negative one's distances. Also for 122 % the "islands" with -Inf there, this is fine, as islands mean that the 123 % respective points are never reached by the front, thus never taken 124 % out of the domain, and thus D < -t should always be true for them 125 % and all arbitrarily large times t. 126 dp = solvePart (d0, fp, h, fabs); 127 % Be careful with -NA ~= NA! 128 invD0 = -d0; 129 invD0(isna (d0)) = NA; 130 dn = -solvePart (invD0, fn, h, fabs); 131 132 % Build everything together. 133 d = NA (size (phi)); 134 d(fp) = dp(fp); 135 d(fn) = dn(fn); 136endfunction 137 138% D = solvePart (D0, REGION, H, FABS) 139% 140% Utility function that performs one of the two solves required for the 141% regions with positive and negative f. D0 should be the initial distances 142% (all positive), REGION a boolean flag of the region we process (f > 0 143% or f < 0), H the grid spacing and FABS the absolute value of f. 144% 145% We solve in the given REGION using fast marching for the distances 146% of points that are eventually (but not initially) reached by the front, 147% and also handle the case of never-reached "islands". 148 149function d = solvePart (d0, region, h, fabs) 150 d = d0; 151 d(~region) = Inf; 152 fval = NA (size (fabs)); 153 fval(region) = h ./ fabs(region); 154 d = fastmarching (d, fval); 155 156 % Set distances from which the front moved away from Inf to -Inf. 157 d(d == Inf & region) = -Inf; 158 159 % If a region is not reached by the fast marching (i. e., NA), this means 160 % that it is separated from the boundary by a sign change of f. This also 161 % means that it will actually never be reached, and thus we can set its 162 % distance to Inf. 163 d(isna (d)) = Inf; 164 165 assert (all (~isnan (d(:)))); 166endfunction 167 168%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 169% Tests. 170 171% Test for error conditions. 172%!error <Invalid call to> 173%! ls_solve_stationary (1) 174%!error <Invalid call to> 175%! ls_solve_stationary (1, 2, 3, 4, 5) 176%!error <PHI and F must be of the same size> 177%! ls_solve_stationary ([1, 2], [3; 4]); 178%!error <PHI and NB must be of the same size> 179%! ls_solve_stationary ([1, 2], [1, 2], :, [3; 4]); 180 181% Check 0D case. 182%!test 183%! assert (ls_solve_stationary ([], []), []); 184 185% Check invariants expected for the output in terms 186% of finite/infinite values and signs. 187%!function checkOutput (phi, f, d) 188%! assert (all (isna (d(f == 0)))); 189%! where = (f ~= 0); 190%! assert (all (isfinite (d(where)) | isinf (d(where)))); 191%! assert (sign (d(where)), sign (phi(where))); 192%!endfunction 193 194% For the test phi's, just check that solving works 195% without any errors. 196%!test 197%! phis = ls_get_tests (); 198%! for i = 1 : length (phis) 199%! f = ones (size (phis{i})); 200%! for j = -1 : 1 201%! phi = ls_normalise (phis{i}); 202%! d = ls_solve_stationary (phi, j * f); 203%! checkOutput (phi, j * f, d); 204%! endfor 205%! endfor 206 207% Test for island handling and the general output conditions 208% (that phi = 0 leads to D = NA and otherwise we only ever get 209% some ordinary values with signs depending on the sign of f 210% or infinities with appropriate signs). 211% 212% The constructed example contains in each half-plane (left / right of the 213% y-axis) an initial circle as domain, and a "ring" around the domain boundary 214% where f is positive or negative. This ensures we have a sign change of f 215% both in the interior and exterior of the initial domain. Depending on the 216% sign of f (and thus the half-plane), one of them will be an "island". 217% 218%!test 219%! n = 101; 220%! x = linspace (-10, 10, n); 221%! h = x(2) - x(1); 222%! [XX, YY] = meshgrid (x, x); 223%! 224%! r1 = (XX - 5).^2 + YY.^2; 225%! r2 = (XX + 5).^2 + YY.^2; 226%! 227%! phi = ls_union (r1 - 2^2, r2 - 2^2); 228%! f1 = max (r1 - 3^2, 1^2 - r1); 229%! f2 = max (r2 - 3^2, 1^2 - r2); 230%! f = min (f1, f2) .* sign (XX); 231%! assert (~all (f(:) ~= 0)); 232%! 233%! d = ls_solve_stationary (phi, f, h); 234%! 235%! checkOutput (phi, f, d); 236%! where = (f ~= 0); 237%! assert (any (isinf (d(where)))); 238%! 239%! %{ 240%! contour (phi, [0, 0]); 241%! colorbar (); 242%! figure (); 243%! imagesc (f); 244%! colorbar (); 245%! figure (); 246%! imagesc (d); 247%! colorbar (); 248%! %} 249 250% Test manual NB initialisation. 251%!test 252%! phi = [3, 1, -1, 1, 3]; 253%! f = ones (size (phi)); 254%! dists1 = ls_solve_stationary (phi, f); 255%! dists2 = ls_solve_stationary (phi, f, :, [NA, NA, -0.1, NA, NA]); 256%! 257%! assert (dists1, [1.5, 0.5, -0.5, 0.5, 1.5], sqrt (eps)); 258%! assert (dists2, [1.9, 0.9, -0.1, 0.9, 1.9]); 259 260%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 261% Demos. 262 263% Visualise refraction and diffraction of a wave front. If we only 264% use positive speeds, then the different times can be plotted nicely 265% just with 'contour'. 266% 267%!demo 268%! n = 100; 269%! x = linspace (-10, 10, n); 270%! h = x(2) - x(1); 271%! [XX, YY] = meshgrid (x, x); 272%! 273%! phi0 = YY - XX + 15; 274%! f = sign (XX) + 2; 275%! 276%! d = ls_solve_stationary (phi0, f, h); 277%! figure (); 278%! contour (XX, YY, d); 279%! line ([0, 0], [-10, 10]); 280%! title ("Refraction"); 281% 282%!demo 283%! n = 101; 284%! x = linspace (-10, 10, n); 285%! h = x(2) - x(1); 286%! [XX, YY] = meshgrid (x, x); 287%! 288%! phi0 = XX + 9; 289%! f = ones (size (phi0)); 290%! f(XX == 0 & abs (YY) > 2) = 0; 291%! 292%! d = ls_solve_stationary (phi0, f, h); 293%! figure (); 294%! contour (XX, YY, d); 295%! line ([0, 0], [-10, -2]); 296%! line ([0, 0], [2, 10]); 297%! title ("Diffraction"); 298 299% XXX: Can we use this for ray-tracing as a nice example??? 300