Mercurial > octave

comparison libinterp/corefcn/tsearch.cc @ 31277:185799b2a566

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
tsearch.cc: Minor performance improvements tsearch.cc: replaced if-else with conditional operator in inline functions, replaced redeclarations of local variables with one set of declarations and initializations, added more comments. For searching 1e3 points through 1e6 triangles, these minor changes reduced 1.55-1.75 seconds to 1.35-1.45 seconds, so roughly 15% faster. Still O(M*N) for M points and N triangles, but acceptably fast for common use.
author Arun Giridhar <arungiridhar@gmail.com>
date Sat, 08 Oct 2022 13:14:02 -0400
parents 796f54d4ddbf
children e88a07dec498
comparison
equal deleted inserted replaced
31276:773a361bc476 31277:185799b2a566
37 37
38 OCTAVE_NAMESPACE_BEGIN 38 OCTAVE_NAMESPACE_BEGIN
39 39
40 inline double max (double a, double b, double c) 40 inline double max (double a, double b, double c)
41 { 41 {
42 if (a < b) 42 return (a > b) ? (a > c ? a : c) : (b > c ? b : c);
43 return (b < c ? c : b);
44 else
45 return (a < c ? c : a);
46 } 43 }
47 44
48 inline double min (double a, double b, double c) 45 inline double min (double a, double b, double c)
49 { 46 {
50 if (a > b) 47 return (a < b) ? (a < c ? a : c) : (b < c ? b : c);
51 return (b > c ? c : b);
52 else
53 return (a > c ? c : a);
54 } 48 }
55 49
56 #define REF(x,k,i) x(static_cast<octave_idx_type> (elem((k), (i))) - 1) 50 #define REF(x,k,i) x(static_cast<octave_idx_type> (elem((k), (i))) - 1)
57 51
58 // for large data set the algorithm is very slow 52 // for large data set the algorithm is very slow
100 } 94 }
101 95
102 const octave_idx_type np = xi.numel (); 96 const octave_idx_type np = xi.numel ();
103 ColumnVector values (np); 97 ColumnVector values (np);
104 98
105 double x0, y0, a11, a12, a21, a22, det; 99 double x0 = 0.0, y0 = 0.0;
106 x0 = y0 = 0.0; 100 double a11 = 0.0, a12 = 0.0, a21 = 0.0, a22 = 0.0, det = 0.0;
107 a11 = a12 = a21 = a22 = 0.0; 101 double xt = 0.0, yt = 0.0;
108 det = 0.0; 102 double dx1 = 0.0, dx2 = 0.0, c1 = 0.0, c2 = 0.0;
109 103
110 octave_idx_type k = nelem; // k is a counter of elements 104 octave_idx_type k = nelem; // k is more than just an index variable.
111 for (octave_idx_type kp = 0; kp < np; kp++) 105
106 for (octave_idx_type kp = 0; kp < np; kp++) // for each point
112 { 107 {
113 const double xt = xi(kp); 108 xt = xi (kp);
114 const double yt = yi(kp); 109 yt = yi (kp);
115 110
116 // check if last triangle contains the next point 111 // Check if point (xt,yt) is in the triangle that was last examined.
117 if (k < nelem) 112 // This is for inputs where points are in contiguous order,
113 // like when the points are sampled from a continuous path.
114 if (k < nelem) // This check will be false for the very first point.
118 { 115 {
119 const double dx1 = xt - x0; 116 // If we are here, then x0, y0, det all exist from before.
120 const double dx2 = yt - y0; 117 dx1 = xt - x0;
121 const double c1 = (a22 * dx1 - a21 * dx2) / det; 118 dx2 = yt - y0;
122 const double c2 = (-a12 * dx1 + a11 * dx2) / det; 119 c1 = (a22 * dx1 - a21 * dx2) / det;
123 if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= (1 + eps)) 120 c2 = (-a12 * dx1 + a11 * dx2) / det;
121 if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= 1 + eps)
124 { 122 {
125 values(kp) = double(k+1); 123 values (kp) = k+1;
126 continue; 124 continue;
127 } 125 }
128 } 126 }
129 127
130 // it doesn't, so go through all elements 128 // The point is not in the same triangle, so go through all triangles.
131 for (k = 0; k < nelem; k++) 129 for (k = 0; k < nelem; k++)
132 { 130 {
133 octave_quit ();
134
135 if (xt >= minx(k) && xt <= maxx(k) && yt >= miny(k) && yt <= maxy(k)) 131 if (xt >= minx(k) && xt <= maxx(k) && yt >= miny(k) && yt <= maxy(k))
136 { 132 {
137 // element inside the minimum rectangle: examine it closely 133 // Point is inside the triangle's bounding rectangle:
134 // See if it's inside the triangle itself.
138 x0 = REF (x, k, 0); 135 x0 = REF (x, k, 0);
139 y0 = REF (y, k, 0); 136 y0 = REF (y, k, 0);
140 a11 = REF (x, k, 1) - x0; 137 a11 = REF (x, k, 1) - x0;
141 a12 = REF (y, k, 1) - y0; 138 a12 = REF (y, k, 1) - y0;
142 a21 = REF (x, k, 2) - x0; 139 a21 = REF (x, k, 2) - x0;
143 a22 = REF (y, k, 2) - y0; 140 a22 = REF (y, k, 2) - y0;
144 det = a11 * a22 - a21 * a12; 141 det = a11 * a22 - a21 * a12;
145 142
146 // solve the system 143 // solve the system
147 const double dx1 = xt - x0; 144 dx1 = xt - x0;
148 const double dx2 = yt - y0; 145 dx2 = yt - y0;
149 const double c1 = (a22 * dx1 - a21 * dx2) / det; 146 c1 = (a22 * dx1 - a21 * dx2) / det;
150 const double c2 = (-a12 * dx1 + a11 * dx2) / det; 147 c2 = (-a12 * dx1 + a11 * dx2) / det;
151 if ((c1 >= -eps) && (c2 >= -eps) && ((c1 + c2) <= (1 + eps))) 148 if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= 1 + eps)
152 { 149 {
153 values(kp) = double(k+1); 150 values (kp) = k+1;
154 break; 151 break;
155 } 152 }
156 } //endif # examine this element closely 153 } //end see if it's inside the triangle itself
157 } //endfor # each element 154 } //end for each triangle
158 155
159 if (k == nelem) 156 if (k == nelem)
160 values(kp) = lo_ieee_nan_value (); 157 values (kp) = lo_ieee_nan_value ();
161 158
162 } //endfor # kp 159 } //end for each point
163 160
164 return ovl (values); 161 return ovl (values);
165 } 162 }
166 163
167 /* 164 /*