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