Mercurial > octave

comparison scripts/geometry/tsearchn.m @ 31543:72ef3d097059

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
tsearchn.m: Speed up performance (bug #63376) tsearchn.m: Eliminate a subfunction and inline its contents, eliminate ones() and replace with broadcasting, add input validation, pass to a faster compiled function tsearch() where possible, expand documentation. Cumulative speedup is some 20%.
author Arun Giridhar <arungiridhar@gmail.com>
date Fri, 25 Nov 2022 11:31:47 -0500
parents 796f54d4ddbf
children c664627d601e
comparison
equal deleted inserted replaced
31540:ce5b4a00b022 31543:72ef3d097059
24 ######################################################################## 24 ########################################################################
25 25
26 ## -*- texinfo -*- 26 ## -*- texinfo -*-
27 ## @deftypefn {} {@var{idx} =} tsearchn (@var{x}, @var{t}, @var{xi}) 27 ## @deftypefn {} {@var{idx} =} tsearchn (@var{x}, @var{t}, @var{xi})
28 ## @deftypefnx {} {[@var{idx}, @var{p}] =} tsearchn (@var{x}, @var{t}, @var{xi}) 28 ## @deftypefnx {} {[@var{idx}, @var{p}] =} tsearchn (@var{x}, @var{t}, @var{xi})
29 ## Search for the enclosing Delaunay convex hull. 29 ## Find the simplexes enclosing the given points.
30 ## 30 ##
31 ## For @code{@var{t} = delaunayn (@var{x})}, finds the index in @var{t} 31 ## @code{tsearchn} is typically used with @code{delaunayn}:
32 ## containing the points @var{xi}. For points outside the convex hull, 32 ## @code{@var{t} = delaunayn (@var{x})} returns a set of simplexes @code{t},
33 ## @var{idx} is NaN. 33 ## then @code{tsearchn} returns the row index of @var{t} containing each point
34 ## of @var{xi}. For points outside the convex hull, @var{idx} is NaN.
34 ## 35 ##
35 ## If requested @code{tsearchn} also returns the Barycentric coordinates 36 ## If requested, @code{tsearchn} also returns the barycentric coordinates
36 ## @var{p} of the enclosing triangles. 37 ## @var{p} of the enclosing simplexes.
37 ## @seealso{delaunay, delaunayn} 38 ##
39 ## @seealso{delaunay, delaunayn, tsearch}
38 ## @end deftypefn 40 ## @end deftypefn
39 41
40 function [idx, p] = tsearchn (x, t, xi) 42 function [idx, p] = tsearchn (x, t, xi)
41 43
42 if (nargin != 3) 44 if (nargin != 3)
43 print_usage (); 45 print_usage ();
44 endif 46 endif
45 47
48 if (columns (x) != columns (xi))
49 error ("columns (x) should equal columns (xi)")
50 end
51
52 if (max (t(:)) > rows (x))
53 error ("triangles should only access points in x")
54 end
55
56 if (nargout <= 1 && columns (x) == 2) # pass to the faster tsearch.cc
57 idx = tsearch (x(:,1), x(:,2), t, xi(:,1), xi(:,2));
58 return
59 endif
60
46 nt = rows (t); 61 nt = rows (t);
47 [m, n] = size (x); 62 [m, n] = size (x);
48 mi = rows (xi); 63 mi = rows (xi);
49 idx = NaN (mi, 1); 64 idx = NaN (mi, 1);
50 p = NaN (mi, n + 1); 65 p = NaN (mi, n + 1);
66 ni = [1:mi].';
51 67
52 ni = [1:mi].'; 68 for i = 1 : nt # each simplex in turn
53 for i = 1 : nt
54 ## Only calculate the Barycentric coordinates for points that have not
55 ## already been found in a triangle.
56 b = cart2bary (x (t (i, :), :), xi(ni,:));
57 69
58 ## Our points xi are in the current triangle if 70 T = x(t(i, :), :); # T is the current simplex
59 ## (all (b >= 0) && all (b <= 1)). However as we impose that 71 P = xi(ni, :); # P is the set of points left to calculate
60 ## sum (b,2) == 1 we only need to test all(b>=0). Note need to add 72
61 ## a small margin for rounding errors 73 ## Convert to barycentric coords: these are used to express a point P
62 intri = all (b >= -1e-12, 2); 74 ## as P = Beta * T
63 idx(ni(intri)) = i; 75 ## where T is a simplex.
64 p(ni(intri),:) = b(intri, :); 76 ##
65 ni(intri) = []; 77 ## If 0 <= Beta <= 1, then the linear combination is also convex,
78 ## and the point P is inside the simplex T, otherwise it is outside.
79 ## Since the equation system is underdetermined, we apply the constraint
80 ## sum (Beta) == 1 to make it unique up to scaling.
81 ##
82 ## Note that the code below is vectorized over P, one point per row.
83
84 b = (P - T(end,:)) / (T(1:end-1,:) - T(end,:));
85 b(:, end+1) = 1 - sum (b, 2);
86
87 ## The points xi are inside the current simplex if
88 ## (all (b >= 0) && all (b <= 1)). As sum (b,2) == 1, we only need to
89 ## test all(b>=0).
90 inside = all (b >= -1e-12, 2); # -1e-12 instead of 0 for rounding errors
91 idx (ni (inside)) = i;
92 p(ni(inside), :) = b(inside, :);
93 ni = ni (~inside);
66 endfor 94 endfor
67 95
68 endfunction 96 endfunction
69
70 function Beta = cart2bary (T, P)
71
72 ## Conversion of Cartesian to Barycentric coordinates.
73 ## Given a reference simplex in N dimensions represented by an
74 ## N+1-by-N matrix, an arbitrary point P in Cartesian coordinates,
75 ## represented by an N-by-1 column vector can be written as
76 ##
77 ## P = Beta * T
78 ##
79 ## Where Beta is an N+1 vector of the barycentric coordinates. A criteria
80 ## on Beta is that
81 ##
82 ## sum (Beta) == 1
83 ##
84 ## and therefore we can write the above as
85 ##
86 ## P - T(end, :) = Beta(1:end-1) * (T(1:end-1,:) - ones (N,1) * T(end,:))
87 ##
88 ## and then we can solve for Beta as
89 ##
90 ## Beta(1:end-1) = (P - T(end,:)) / (T(1:end-1,:) - ones (N,1) * T(end,:))
91 ## Beta(end) = sum (Beta)
92 ##
93 ## Note code below is generalized for multiple values of P, one per row.
94 [M, N] = size (P);
95 Beta = (P - ones (M,1) * T(end,:)) / (T(1:end-1,:) - ones (N,1) * T(end,:));
96 Beta (:,end+1) = 1 - sum (Beta, 2);
97
98 endfunction
99
100 97
101 %!shared x, tri 98 %!shared x, tri
102 %! x = [-1,-1;-1,1;1,-1]; 99 %! x = [-1,-1;-1,1;1,-1];
103 %! tri = [1, 2, 3]; 100 %! tri = [1, 2, 3];
104 %!test 101 %!test