sln_polyline.c (9313B)
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 /* Compute the last parameter of the edge equation */ 104 C = -f2_dot(N, p0); 105 106 imax = range[0]+1; 107 err_max = 0; 108 FOR_EACH(ivertex, range[0]+1, range[1]) { 109 float p[2]; 110 float val; 111 float err; 112 113 FETCH_VERTEX(p, ivertex); 114 115 if(N[0]*p[0] + N[1]*p[1] + C < 0) { /* The vertex is above the edge */ 116 has_vertices_above = 1; 117 } 118 119 /* Compute the linear approximation of p */ 120 val = -(N[0]*p[0] + C)/N[1]; 121 122 /* Compute the relative error of the linear approximation of p */ 123 err = absf(val - p[1]) / p[1]; 124 if(err > err_max) { 125 imax = ivertex; 126 err_max = err; 127 } 128 } 129 #undef FETCH_VERTEX 130 131 *out_ivertex = imax; 132 /* To ensure an upper mesh, we cannot delete a vertex above the candidate 133 * edge used to simplify the mesh. We therefore compel ourselves not to 134 * simplify the polyline when such vertices are detected by returning an 135 * infinite error */ 136 if(mesh_type == SLN_MESH_UPPER && has_vertices_above) { 137 *out_err = (float)INF; 138 } else { 139 *out_err = err_max; 140 } 141 } 142 143 static INLINE void 144 check_polyline_vertices 145 (const struct sln_vertex* vertices, 146 const size_t vertices_range[2]) 147 { 148 #ifdef NDEBUG 149 (void)vertices, (void)vertices_range; 150 #else 151 size_t i; 152 ASSERT(vertices); 153 FOR_EACH(i, vertices_range[0]+1, vertices_range[1]+1) { 154 CHK(vertices[i].wavenumber >= vertices[i-1].wavenumber); 155 } 156 #endif 157 } 158 159 /******************************************************************************* 160 * Local function 161 ******************************************************************************/ 162 /* In place simplification of a polyline. Given a curve composed of line 163 * segments, compute a similar curve with fewer points. In the following we 164 * implement the algorithm described in: 165 * 166 * "Algorithms for the reduction of the number of points required to 167 * represent a digitized line or its caricature" - David H. Douglas and Thomas 168 * K Peucker, Cartographica: the international journal for geographic 169 * information and geovisualization - 1973 */ 170 res_T 171 polyline_decimate 172 (struct sln_device* sln, 173 struct sln_vertex* vertices, 174 size_t vertices_range[2], 175 const float err, /* Max relative error */ 176 const enum sln_mesh_type mesh_type) 177 { 178 struct darray_size_t stack; 179 size_t range[2] = {0, 0}; 180 size_t ivtx = 0; 181 size_t nvertices = 0; 182 res_T res = RES_OK; 183 ASSERT(vertices && vertices_range && err >= 0); 184 ASSERT(vertices_range[0] < vertices_range[1]); 185 check_polyline_vertices(vertices, vertices_range); 186 187 darray_size_t_init(sln->allocator, &stack); 188 189 nvertices = vertices_range[1] - vertices_range[0] + 1; 190 if(nvertices <= 2 || err == 0) goto exit; /* Nothing to simplify */ 191 192 /* Helper macros */ 193 #define PUSH(Stack, Val) { \ 194 res = darray_size_t_push_back(&(Stack), &(Val)); \ 195 if(res != RES_OK) goto error; \ 196 } (void)0 197 #define POP(Stack, Val) { \ 198 const size_t sz = darray_size_t_size_get(&(Stack)); \ 199 ASSERT(sz); \ 200 (Val) = darray_size_t_cdata_get(&(Stack))[sz-1]; \ 201 darray_size_t_resize(&(Stack), sz-1); \ 202 } (void)0 203 204 range[0] = vertices_range[0]; 205 range[1] = vertices_range[1]; 206 ivtx = vertices_range[0] + 1; 207 208 /* Push a dummy entry to allow "stack pop" once range[0] reaches range[1] */ 209 PUSH(stack, range[1]); 210 211 while(range[0] != vertices_range[1]) { 212 /* Epsilon below which vertex attributes are considered equal */ 213 const float vtx_eps = 1.e-6f; 214 215 size_t imax; 216 float err_max = -1; 217 218 if(is_polyline_degenerate(vertices, range, vtx_eps)) { 219 /* All the vertices of the polyline are merged. Delete all vertices from 220 * the polyline, i.e., do not save any vertices. Then, move to the next 221 * vertex range */ 222 range[0] = range[1]; 223 POP(stack, range[1]); 224 225 /* Note that this simplification removes vertices (even if they are very 226 * close), and the resulting mesh may therefore not correspond to the 227 * upper limit of the spectrum required by the caller. To try to mitigate 228 * this effect, update the retained vertex by adding the epsilon used to 229 * detect that one vertex is equal to another. 230 * 231 * However, there is no strict guarantee that the mesh constitutes an 232 * upper limit. This should be kept in mind, but it should be assessed in 233 * light of the conditions under which an upward default is possible. The 234 * processed mesh is likely to be superior to the spectrum well beyond the 235 * numerical approximation related to the aforementioned epsilon, which 236 * must nevertheless remain tiny */ 237 if(mesh_type == SLN_MESH_UPPER) { 238 ASSERT(ivtx > 0); 239 vertices[ivtx-1].ka += vtx_eps; 240 } 241 242 } else { 243 if(range[1] - range[0] > 1) { 244 find_falsest_vertex(vertices, range, mesh_type, &imax, &err_max); 245 } 246 247 if(err_max > err) { 248 /* Try to simplify a smaller polyline interval in [range[0], imax] */ 249 PUSH(stack, range[1]); 250 range[1] = imax; 251 } else { 252 /* Remove all vertices in [range[0]+1, range[1]-1]] */ 253 vertices[ivtx] = vertices[range[1]]; 254 ++ivtx; 255 256 /* Setup the next range */ 257 range[0] = range[1]; 258 POP(stack, range[1]); 259 } 260 } 261 } 262 #undef PUSH 263 #undef POP 264 265 vertices_range[1] = ivtx - 1; 266 267 check_polyline_vertices(vertices, vertices_range); 268 269 exit: 270 darray_size_t_release(&stack); 271 return res; 272 error: 273 goto exit; 274 }