sln_polyline.c (7005B)
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 /* Given a set of points in [range[0], range[1]], find the point in 31 * [range[0]+1, range[0]-1] whose maximize the error regarding its linear 32 * approximation by the line (range[0], range[1]). Let K the real value of the point 33 * and Kl its linear approximation, the error is computed as: 34 * 35 * err = |K-Kl|/K */ 36 static void 37 find_falsest_vertex 38 (const struct sln_vertex* vertices, 39 const size_t range[2], 40 const enum sln_mesh_type mesh_type, 41 size_t* out_ivertex, /* Index of the falsest vertex */ 42 float* out_err) /* Error of the falsest vertex */ 43 { 44 float p0[2], p1[2]; /* 1st and last point of the submitted range */ 45 float N[2], C; /* edge equation N.p + C = 0 */ 46 float len; 47 size_t ivertex; 48 size_t imax; /* Index toward the falsest vertex */ 49 float err_max; 50 int has_vertices_above = 0; 51 ASSERT(vertices && range && range[0] < range[1]-1 && out_ivertex && out_err); 52 ASSERT((unsigned)mesh_type < SLN_MESH_TYPES_COUNT__); 53 (void)len; 54 55 #define FETCH_VERTEX(Dst, Id) { \ 56 (Dst)[0] = vertices[(Id)].wavenumber; \ 57 (Dst)[1] = vertices[(Id)].ka; \ 58 } (void)0 59 FETCH_VERTEX(p0, range[0]); 60 FETCH_VERTEX(p1, range[1]); 61 62 /* Compute the normal of the edge [p0, p1] 63 * N[0] = (p1 - p0).y 64 * N[1] =-(p1 - p0).x */ 65 N[0] = p1[1] - p0[1]; 66 N[1] = p0[0] - p1[0]; 67 len = f2_normalize(N, N); 68 ASSERT(len > 0); 69 70 /* Compute the last parameter of the edge equation */ 71 C = -f2_dot(N, p0); 72 73 imax = range[0]+1; 74 err_max = 0; 75 FOR_EACH(ivertex, range[0]+1, range[1]) { 76 float p[2]; 77 float val; 78 float err; 79 80 FETCH_VERTEX(p, ivertex); 81 82 if(N[0]*p[0] + N[1]*p[1] + C < 0) { /* The vertex is above the edge */ 83 has_vertices_above = 1; 84 } 85 86 /* Compute the linear approximation of p */ 87 val = -(N[0]*p[0] + C)/N[1]; 88 89 /* Compute the relative error of the linear approximation of p */ 90 err = absf(val - p[1]) / p[1]; 91 if(err > err_max) { 92 imax = ivertex; 93 err_max = err; 94 } 95 } 96 #undef FETCH_VERTEX 97 98 *out_ivertex = imax; 99 /* To ensure an upper mesh, we cannot delete a vertex above the candidate 100 * edge used to simplify the mesh. We therefore compel ourselves not to 101 * simplify the polyline when such vertices are detected by returning an 102 * infinite error */ 103 if(mesh_type == SLN_MESH_UPPER && has_vertices_above) { 104 *out_err = (float)INF; 105 } else { 106 *out_err = err_max; 107 } 108 } 109 110 static INLINE void 111 check_polyline_vertices 112 (const struct sln_vertex* vertices, 113 const size_t vertices_range[2]) 114 { 115 #ifdef NDEBUG 116 (void)vertices, (void)vertices_range; 117 #else 118 size_t i; 119 ASSERT(vertices); 120 FOR_EACH(i, vertices_range[0]+1, vertices_range[1]+1) { 121 CHK(vertices[i].wavenumber >= vertices[i-1].wavenumber); 122 } 123 #endif 124 } 125 126 /******************************************************************************* 127 * Local function 128 ******************************************************************************/ 129 /* In place simplification of a polyline. Given a curve composed of line 130 * segments, compute a similar curve with fewer points. In the following we 131 * implement the algorithm described in: 132 * 133 * "Algorithms for the reduction of the number of points required to 134 * represent a digitized line or its caricature" - David H. Douglas and Thomas 135 * K Peucker, Cartographica: the international journal for geographic 136 * information and geovisualization - 1973 */ 137 res_T 138 polyline_decimate 139 (struct sln_device* sln, 140 struct sln_vertex* vertices, 141 size_t vertices_range[2], 142 const float err, /* Max relative error */ 143 const enum sln_mesh_type mesh_type) 144 { 145 struct darray_size_t stack; 146 size_t range[2] = {0, 0}; 147 size_t ivtx = 0; 148 size_t nvertices = 0; 149 res_T res = RES_OK; 150 ASSERT(vertices && vertices_range && err >= 0); 151 ASSERT(vertices_range[0] < vertices_range[1]); 152 check_polyline_vertices(vertices, vertices_range); 153 154 darray_size_t_init(sln->allocator, &stack); 155 156 nvertices = vertices_range[1] - vertices_range[0] + 1; 157 if(nvertices <= 2 || err == 0) goto exit; /* Nothing to simplify */ 158 159 /* Helper macros */ 160 #define PUSH(Stack, Val) { \ 161 res = darray_size_t_push_back(&(Stack), &(Val)); \ 162 if(res != RES_OK) goto error; \ 163 } (void)0 164 #define POP(Stack, Val) { \ 165 const size_t sz = darray_size_t_size_get(&(Stack)); \ 166 ASSERT(sz); \ 167 (Val) = darray_size_t_cdata_get(&(Stack))[sz-1]; \ 168 darray_size_t_resize(&(Stack), sz-1); \ 169 } (void)0 170 171 range[0] = vertices_range[0]; 172 range[1] = vertices_range[1]; 173 ivtx = vertices_range[0] + 1; 174 175 /* Push a dummy entry to allow "stack pop" once range[0] reaches range[1] */ 176 PUSH(stack, range[1]); 177 178 while(range[0] != vertices_range[1]) { 179 size_t imax; 180 float err_max = -1; 181 182 if(range[1] - range[0] > 1) { 183 find_falsest_vertex(vertices, range, mesh_type, &imax, &err_max); 184 } 185 186 if(err_max > err) { 187 /* Try to simplify a smaller polyline interval in [range[0], imax] */ 188 PUSH(stack, range[1]); 189 range[1] = imax; 190 } else { 191 /* Remove all vertices in [range[0]+1, range[1]-1]] */ 192 vertices[ivtx++] = vertices[range[1]]; 193 194 /* Setup the next range */ 195 range[0] = range[1]; 196 POP(stack, range[1]); 197 } 198 } 199 #undef PUSH 200 #undef POP 201 202 vertices_range[1] = ivtx - 1; 203 204 check_polyline_vertices(vertices, vertices_range); 205 206 exit: 207 darray_size_t_release(&stack); 208 return res; 209 error: 210 goto exit; 211 }