sln_polyline.c (10063B)
1 /* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2026 Université de Lorraine 3 * Copyright (C) 2022 Centre National de la Recherche Scientifique 4 * Copyright (C) 2022 Université Paul Sabatier 5 * 6 * This program is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 18 19 #include "sln.h" 20 #include "sln_device_c.h" 21 #include "sln_polyline.h" 22 23 #include <rsys/dynamic_array_size_t.h> 24 #include <rsys/float2.h> 25 #include <rsys/math.h> 26 27 /******************************************************************************* 28 * Helper function 29 ******************************************************************************/ 30 static INLINE int 31 vtx_eq_eps 32 (const struct sln_vertex* v0, 33 const struct sln_vertex* v1, 34 const float eps) 35 { 36 ASSERT(v0 && v1); 37 38 return eq_eps(v0->wavenumber, v1->wavenumber, eps) 39 && eq_eps(v0->ka, v1->ka, eps); 40 } 41 42 /* Defines whether a polyline is degenerate, i.e., whether it is actually a 43 * point */ 44 static int 45 is_polyline_degenerate 46 (const struct sln_vertex* vertices, 47 const size_t range[2], /* Bounds are inclusive */ 48 const float epsilon) 49 { 50 size_t i = 0; 51 ASSERT(range[0] < range[1] && epsilon >= 0); 52 53 /* Find the first peak that is not (approximately) equal to the first one */ 54 FOR_EACH(i, range[0]+1, range[1]+1/*Inclusive bound*/) { 55 if(!vtx_eq_eps(vertices+range[0], vertices + i, epsilon)) break; 56 } 57 58 /* If all vertices are verified without interruption, then all vertices are 59 * (approximately) equal. The polyline is therefore degenerate. */ 60 return i > range[1]; 61 } 62 63 /* Given a set of points in [range[0], range[1]], find the point in 64 * [range[0]+1, range[0]-1] whose maximize the error regarding its linear 65 * approximation by the line (range[0], range[1]). Let K the real value of the point 66 * and Kl its linear approximation, the error is computed as: 67 * 68 * err = |K-Kl|/K */ 69 static void 70 find_falsest_vertex 71 (const struct sln_vertex* vertices, 72 const size_t range[2], 73 const enum sln_mesh_type mesh_type, 74 size_t* out_ivertex, /* Index of the falsest vertex */ 75 float* out_err) /* Error of the falsest vertex */ 76 { 77 float p0[2], p1[2]; /* 1st and last point of the submitted range */ 78 float N[2], C; /* edge equation N.p + C = 0 */ 79 float len; 80 size_t ivertex; 81 size_t imax; /* Index toward the falsest vertex */ 82 float err_max; 83 int has_vertices_above = 0; 84 ASSERT(vertices && range && range[0] < range[1]-1 && out_ivertex && out_err); 85 ASSERT((unsigned)mesh_type < SLN_MESH_TYPES_COUNT__); 86 (void)len; 87 88 #define FETCH_VERTEX(Dst, Id) { \ 89 (Dst)[0] = vertices[(Id)].wavenumber; \ 90 (Dst)[1] = vertices[(Id)].ka; \ 91 } (void)0 92 FETCH_VERTEX(p0, range[0]); 93 FETCH_VERTEX(p1, range[1]); 94 95 /* Compute the normal of the edge [p0, p1] 96 * N[0] = (p1 - p0).y 97 * N[1] =-(p1 - p0).x */ 98 N[0] = p1[1] - p0[1]; 99 N[1] = p0[0] - p1[0]; 100 len = f2_normalize(N, N); 101 ASSERT(len > 0); 102 103 104 if(eq_eps(N[1], 0, 1e-6)) { 105 /* A normal with a value of zero in Y means a vertical segment. In fact, 106 * this is not supposed to happen because a line cannot have 2 different 107 * values for a single wave number. But this can happen due to numerical 108 * errors. In such a situation, all intermediate vertices may be considered 109 * to be located on the same vertical segment. And so, the error of using it 110 * to represent all intermediate vertices can be assumed to be zero. */ 111 *out_ivertex = range[0]; 112 *out_err = 0; 113 return; 114 } 115 116 /* Compute the last parameter of the edge equation */ 117 C = -f2_dot(N, p0); 118 119 imax = range[0]+1; 120 err_max = 0; 121 FOR_EACH(ivertex, range[0]+1, range[1]) { 122 float p[2]; 123 float val; 124 float err; 125 126 FETCH_VERTEX(p, ivertex); 127 128 if(N[0]*p[0] + N[1]*p[1] + C < 0) { /* The vertex is above the edge */ 129 has_vertices_above = 1; 130 } 131 132 /* Compute the linear approximation of p */ 133 val = -(N[0]*p[0] + C)/N[1]; 134 135 /* Compute the relative error of the linear approximation of p */ 136 err = absf(val - p[1]); 137 if(err != 0) { 138 /* Divide by the maximum between the real and the approximated value to 139 * avoid a zero division */ 140 err /= MMAX(p[1], val); 141 } 142 ASSERT(!IS_NaN(err)); 143 144 if(err > err_max) { 145 imax = ivertex; 146 err_max = err; 147 } 148 } 149 #undef FETCH_VERTEX 150 151 *out_ivertex = imax; 152 /* To ensure an upper mesh, we cannot delete a vertex above the candidate 153 * edge used to simplify the mesh. We therefore compel ourselves not to 154 * simplify the polyline when such vertices are detected by returning an 155 * infinite error */ 156 if(mesh_type == SLN_MESH_UPPER && has_vertices_above) { 157 *out_err = (float)INF; 158 } else { 159 *out_err = err_max; 160 } 161 } 162 163 static INLINE void 164 check_polyline_vertices 165 (const struct sln_vertex* vertices, 166 const size_t vertices_range[2]) 167 { 168 #ifdef NDEBUG 169 (void)vertices, (void)vertices_range; 170 #else 171 size_t i; 172 ASSERT(vertices); 173 FOR_EACH(i, vertices_range[0]+1, vertices_range[1]+1) { 174 CHK(vertices[i].wavenumber >= vertices[i-1].wavenumber); 175 } 176 #endif 177 } 178 179 /******************************************************************************* 180 * Local function 181 ******************************************************************************/ 182 /* In place simplification of a polyline. Given a curve composed of line 183 * segments, compute a similar curve with fewer points. In the following we 184 * implement the algorithm described in: 185 * 186 * "Algorithms for the reduction of the number of points required to 187 * represent a digitized line or its caricature" - David H. Douglas and Thomas 188 * K Peucker, Cartographica: the international journal for geographic 189 * information and geovisualization - 1973 */ 190 res_T 191 polyline_decimate 192 (struct sln_device* sln, 193 struct sln_vertex* vertices, 194 size_t vertices_range[2], 195 const float err, /* Max relative error */ 196 const enum sln_mesh_type mesh_type) 197 { 198 struct darray_size_t stack; 199 size_t range[2] = {0, 0}; 200 size_t ivtx = 0; 201 size_t nvertices = 0; 202 res_T res = RES_OK; 203 ASSERT(vertices && vertices_range && err >= 0); 204 ASSERT(vertices_range[0] < vertices_range[1]); 205 check_polyline_vertices(vertices, vertices_range); 206 207 darray_size_t_init(sln->allocator, &stack); 208 209 nvertices = vertices_range[1] - vertices_range[0] + 1; 210 if(nvertices <= 2 || err == 0) goto exit; /* Nothing to simplify */ 211 212 /* Helper macros */ 213 #define PUSH(Stack, Val) { \ 214 res = darray_size_t_push_back(&(Stack), &(Val)); \ 215 if(res != RES_OK) goto error; \ 216 } (void)0 217 #define POP(Stack, Val) { \ 218 const size_t sz = darray_size_t_size_get(&(Stack)); \ 219 ASSERT(sz); \ 220 (Val) = darray_size_t_cdata_get(&(Stack))[sz-1]; \ 221 darray_size_t_resize(&(Stack), sz-1); \ 222 } (void)0 223 224 range[0] = vertices_range[0]; 225 range[1] = vertices_range[1]; 226 ivtx = vertices_range[0] + 1; 227 228 /* Push a dummy entry to allow "stack pop" once range[0] reaches range[1] */ 229 PUSH(stack, range[1]); 230 231 while(range[0] != vertices_range[1]) { 232 /* Epsilon below which vertex attributes are considered equal */ 233 const float vtx_eps = 1.e-6f; 234 235 size_t imax; 236 float err_max = -1; 237 238 if(is_polyline_degenerate(vertices, range, vtx_eps)) { 239 /* All the vertices of the polyline are merged. Delete all vertices from 240 * the polyline, i.e., do not save any vertices. Then, move to the next 241 * vertex range */ 242 range[0] = range[1]; 243 POP(stack, range[1]); 244 245 /* Note that this simplification removes vertices (even if they are very 246 * close), and the resulting mesh may therefore not correspond to the 247 * upper limit of the spectrum required by the caller. To try to mitigate 248 * this effect, update the retained vertex by adding the epsilon used to 249 * detect that one vertex is equal to another. 250 * 251 * However, there is no strict guarantee that the mesh constitutes an 252 * upper limit. This should be kept in mind, but it should be assessed in 253 * light of the conditions under which an upward default is possible. The 254 * processed mesh is likely to be superior to the spectrum well beyond the 255 * numerical approximation related to the aforementioned epsilon, which 256 * must nevertheless remain tiny */ 257 if(mesh_type == SLN_MESH_UPPER) { 258 ASSERT(ivtx > 0); 259 vertices[ivtx-1].ka += vtx_eps; 260 } 261 262 } else { 263 if(range[1] - range[0] > 1) { 264 find_falsest_vertex(vertices, range, mesh_type, &imax, &err_max); 265 } 266 267 if(err_max > err) { 268 /* Try to simplify a smaller polyline interval in [range[0], imax] */ 269 PUSH(stack, range[1]); 270 range[1] = imax; 271 } else { 272 /* Remove all vertices in [range[0]+1, range[1]-1]] */ 273 vertices[ivtx] = vertices[range[1]]; 274 ++ivtx; 275 276 /* Setup the next range */ 277 range[0] = range[1]; 278 POP(stack, range[1]); 279 } 280 } 281 } 282 #undef PUSH 283 #undef POP 284 285 vertices_range[1] = ivtx - 1; 286 287 check_polyline_vertices(vertices, vertices_range); 288 289 exit: 290 darray_size_t_release(&stack); 291 return res; 292 error: 293 goto exit; 294 }