commit 8478852f51c71b9fe3cfa1d74e9cafbcedc07b78
parent 80ccb818abcc9d71d43deb824176265d9f44dcce
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Fri, 6 May 2022 16:40:35 +0200
Decimate the line mesh
Implement the Douglas & Peucker polyline decimation algorithm.
Diffstat:
7 files changed, 254 insertions(+), 8 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -52,13 +52,15 @@ set(SLN_FILES_SRC
sln_device.c
sln_line.c
sln_log.c
+ sln_polyline.c
sln_tree.c
sln_tree_build.c
sln_voigt.c)
set(SLN_FILES_INC
sln_device_c.h
sln_line.h
- sln_log.h)
+ sln_log.h
+ sln_polyline.h)
set(SLN_FILES_INC_API
sln.h)
diff --git a/src/sln.h b/src/sln.h
@@ -108,11 +108,15 @@ struct sln_tree_create_args {
/* Hint on the number of vertices around the the line center */
size_t nvertices_hint;
+ /* Relative error used to simplify the spectrum mesh. The larger it is, the
+ * coarser the mesh */
+ float mesh_decimation_err;
+
enum sln_line_profile line_profile;
};
#define SLN_TREE_CREATE_ARGS_DEFAULT__ \
{NULL, NULL, {SLN_MOLECULE_NULL__}, 0, {0, DBL_MAX}, 0, 0, 1, \
- 16, SLN_LINE_PROFILE_VOIGT}
+ 16, 0.01f, SLN_LINE_PROFILE_VOIGT}
static const struct sln_tree_create_args SLN_TREE_CREATE_ARGS_DEFAULT =
SLN_TREE_CREATE_ARGS_DEFAULT__;
diff --git a/src/sln_polyline.c b/src/sln_polyline.c
@@ -0,0 +1,176 @@
+/* Copyright (C) 2022 CNRS - LMD
+ * Copyright (C) 2022 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2022 Université Paul Sabatier - IRIT
+ * Copyright (C) 2022 Université Paul Sabatier - Laplace
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "sln.h"
+#include "sln_device_c.h"
+#include "sln_polyline.h"
+
+#include <rsys/dynamic_array_size_t.h>
+#include <rsys/float2.h>
+#include <rsys/math.h>
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+/* Given a set of points in [range[0], range[1]], find the point in
+ * [range[0]+1, range[0]-1] whose maximize the error regarding its linear
+ * approximation by the line (range[0], range[1]). Let K the real value of the point
+ * and Kl its linear approximation, the error is computed as:
+ *
+ * err = |K-K|/min(K, Kl) */
+static void
+find_falsest_vertex
+ (const struct sln_vertex* vertices,
+ const size_t range[2],
+ size_t* out_ivertex, /* Index of the falsest vertex */
+ float* out_err) /* Error of the falsest vertex */
+{
+ float p0[2], p1[2]; /* 1st and last point of the submitted range */
+ float N[2], C; /* edge equation N.p + C = 0 */
+ float len;
+
+ size_t ivertex;
+ size_t imax; /* Index toward the falsest vertex */
+ float err_max;
+ ASSERT(vertices && range && range[0] < range[1]-1 && out_ivertex && out_err);
+ (void)len;
+
+ #define FETCH_VERTEX(Dst, Id) { \
+ (Dst)[0] = vertices[(Id)].wavenumber; \
+ (Dst)[1] = vertices[(Id)].ka; \
+ } (void)0
+ FETCH_VERTEX(p0, range[0]);
+ FETCH_VERTEX(p1, range[1]);
+
+ /* Compute the normal of the edge [p0, p1]
+ * N[0] = (p1 - p0).y
+ * N[1] =-(p1 - p0).x */
+ N[0] = p1[1] - p0[1];
+ N[1] = p0[0] - p1[0];
+ len = f2_normalize(N, N);
+ ASSERT(len > 0);
+
+ /* Compute the last parameter of the edge equation */
+ C = -f2_dot(N, p0);
+
+ imax = 0, err_max = 0;
+ FOR_EACH(ivertex, range[0]+1, range[1]) {
+ float p[2];
+ float val;
+ float err;
+
+ FETCH_VERTEX(p, ivertex);
+
+ /* Compute the linear approximation of p */
+ val = -(N[0]*p[0] + C)/N[1];
+
+ /* Compute the relative error of the linear approximation of p */
+ err = absf(val - p[1]) / MMIN(val, p[1]);
+ if(err > err_max) {
+ imax = ivertex;
+ err_max = err;
+ }
+ }
+ #undef FETCH_VERTEX
+
+ *out_ivertex = imax;
+ *out_err = err_max;
+}
+
+/*******************************************************************************
+ * Local function
+ ******************************************************************************/
+/* In place simplification of a polyline. Given a curve composed of line
+ * segments, compute a similar curve with fewer points. In the following we
+ * implement the algorithm described in:
+ *
+ * "Algorithms for the reduction of the number of points required to
+ * represent a digitized line or its caricature" - David H. Douglas and Thomas
+ * K Peucker, Cartographica: the international journal for geographic
+ * information and geovisualization - 1973 */
+res_T
+polyline_decimate
+ (struct sln_device* sln,
+ struct sln_vertex* vertices,
+ size_t vertices_range[2],
+ const float err) /* Max relative error */
+{
+ struct darray_size_t stack;
+ size_t range[2] = {0, 0};
+ size_t ivtx = 0;
+ size_t nvertices = 0;
+ res_T res = RES_OK;
+ ASSERT(vertices && vertices_range && err >= 0);
+ ASSERT(vertices_range[0] < vertices_range[1]);
+
+ darray_size_t_init(sln->allocator, &stack);
+
+ nvertices = vertices_range[1] - vertices_range[0] + 1;
+ if(nvertices <= 2 || err == 0) goto exit; /* Nothing to simplify */
+
+ /* Helper macros */
+ #define PUSH(Stack, Val) { \
+ res = darray_size_t_push_back(&(Stack), &(Val)); \
+ if(res != RES_OK) goto error; \
+ } (void)0
+ #define POP(Stack, Val) { \
+ const size_t sz = darray_size_t_size_get(&(Stack)); \
+ ASSERT(sz); \
+ (Val) = darray_size_t_cdata_get(&(Stack))[sz-1]; \
+ darray_size_t_resize(&(Stack), sz-1); \
+ } (void)0
+
+ range[0] = vertices_range[0];
+ range[1] = vertices_range[1];
+ ivtx = vertices_range[0] + 1;
+
+ /* Push a dummy entry to allow "stack pop" once range[0] reaches range[1] */
+ PUSH(stack, range[1]);
+
+ while(range[0] != vertices_range[1]) {
+ size_t imax;
+ float err_max = -1;
+
+ if(range[1] - range[0] > 1) {
+ find_falsest_vertex(vertices, range, &imax, &err_max);
+ }
+
+ if(err_max > err) {
+ /* Try to simplify a smaller polyline interval in [range[0], imax] */
+ PUSH(stack, range[1]);
+ range[1] = imax;
+ } else {
+ /* Remove all vertices in [range[0]+1, range[1]-1]] */
+ vertices[ivtx++] = vertices[range[1]];
+
+ /* Setup the next range */
+ range[0] = range[1];
+ POP(stack, range[1]);
+ }
+ }
+ #undef PUSH
+ #undef POP
+
+ vertices_range[1] = ivtx - 1;
+
+exit:
+ darray_size_t_release(&stack);
+ return res;
+error:
+ goto exit;
+}
diff --git a/src/sln_polyline.h b/src/sln_polyline.h
@@ -0,0 +1,37 @@
+/* Copyright (C) 2022 CNRS - LMD
+ * Copyright (C) 2022 |Meso|Star> (contact@meso-star.com)
+ * Copyright (C) 2022 Université Paul Sabatier - IRIT
+ * Copyright (C) 2022 Université Paul Sabatier - Laplace
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef SLN_POLYLINE_H
+#define SLN_POLYLINE_H
+
+#include <rsys/rsys.h>
+
+/* Forward declaration */
+struct sln_device;
+struct sln_vertex;
+
+/* In place simplification of a polyline. Given a curve composed of line
+ * segments, compute a similar curve with fewer points. */
+extern LOCAL_SYM res_T
+polyline_decimate
+ (struct sln_device* sln,
+ struct sln_vertex* vertices,
+ size_t vertices_range[2],
+ const float err); /* Relative error */
+
+#endif /* SLN_POLYLINE_H */
diff --git a/src/sln_tree.c b/src/sln_tree.c
@@ -201,6 +201,12 @@ check_sln_tree_create_args
return RES_BAD_ARG;
}
+ if(args->mesh_decimation_err < 0) {
+ log_err(sln, "%s: invalid decimation error %g.\n",
+ func_name, args->mesh_decimation_err);
+ return RES_BAD_ARG;
+ }
+
if((unsigned)args->line_profile >= SLN_LINE_PROFILES_COUNT__) {
log_err(sln, "%s: invalid line profile %d.\n",
func_name, args->line_profile);
diff --git a/src/sln_tree_build.c b/src/sln_tree_build.c
@@ -17,9 +17,14 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include "sln.h"
+#include "sln_log.h"
+#include "sln_polyline.h"
#include "sln_tree_c.h"
+
#include <star/shtr.h>
+#include <rsys/cstr.h>
+
#define STACK_SIZE 64
/*******************************************************************************
@@ -39,20 +44,32 @@ setup_leaf
const struct sln_tree_create_args* args,
struct sln_node* leaf)
{
+ size_t vertices_range[2];
res_T res = RES_OK;
/* Currently assume that we have only one line per leaf */
ASSERT(leaf->range[0] == leaf->range[1]);
- leaf->offset = 0;
-
- leaf->ivertices = darray_vertex_size_get(&tree->vertices);
-
+ /* Line meshing */
+ vertices_range[0] = darray_vertex_size_get(&tree->vertices);
res = line_mesh(tree, leaf->range[0], args->nvertices_hint, &tree->vertices);
if(res != RES_OK) goto error;
+ vertices_range[1] = darray_vertex_size_get(&tree->vertices) - 1;
+
+ /* Decimating the line mesh */
+ res = polyline_decimate(tree->sln, darray_vertex_data_get(&tree->vertices),
+ vertices_range, args->mesh_decimation_err);
+ if(res != RES_OK) {
+ log_err(tree->sln, "Error while decimating the line mesh -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+ darray_vertex_resize(&tree->vertices, vertices_range[1] + 1);
- leaf->nvertices = ui64_to_ui32
- (darray_vertex_size_get(&tree->vertices) - leaf->ivertices);
+ /* Setup leaf data */
+ leaf->offset = 0;
+ leaf->ivertices = vertices_range[0];
+ leaf->nvertices = ui64_to_ui32(vertices_range[1] - vertices_range[0] + 1);
exit:
return res;
diff --git a/src/test_sln_tree.c b/src/test_sln_tree.c
@@ -430,6 +430,10 @@ test_tree
CHK(sln_tree_create(sln, &tree_args, &tree) == RES_BAD_ARG);
tree_args.max_nlines_per_leaf = 1;
+ tree_args.mesh_decimation_err = -1;
+ CHK(sln_tree_create(sln, &tree_args, &tree) == RES_BAD_ARG);
+
+ tree_args.mesh_decimation_err = 1e-1f;
CHK(sln_tree_create(sln, &tree_args, &tree) == RES_OK);
CHK(sln_tree_ref_put(tree) == RES_OK);
