summaryrefslogtreecommitdiff
path: root/meme/slice.c
diff options
context:
space:
mode:
authorWerner Almesberger <werner@almesberger.net>2015-03-14 21:35:46 (GMT)
committerWerner Almesberger <werner@almesberger.net>2015-03-14 21:35:46 (GMT)
commit83e28d001844f160bfe73a35dd03868a9971756a (patch)
tree313080fd5746915113b97ab34109c804a9ebb163 /meme/slice.c
parent6b7491bd5298dd31c09c63bf08adbc1f7013345d (diff)
downloadmisc-83e28d001844f160bfe73a35dd03868a9971756a.zip
misc-83e28d001844f160bfe73a35dd03868a9971756a.tar.gz
misc-83e28d001844f160bfe73a35dd03868a9971756a.tar.bz2
meme/slice.c (slice_plane): use proper linear regression for finding the plane
Luckily, GSL has all the math we need: https://www.gnu.org/software/gsl/manual/html_node/Robust-linear-regression.html#Robust-linear-regression
Diffstat (limited to 'meme/slice.c')
-rw-r--r--meme/slice.c49
1 files changed, 46 insertions, 3 deletions
diff --git a/meme/slice.c b/meme/slice.c
index 2ce8657..9618249 100644
--- a/meme/slice.c
+++ b/meme/slice.c
@@ -16,6 +16,10 @@
#include <stdio.h>
#include <math.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_multifit.h>
+
#include "util.h"
#include "mesh.h"
#include "slice.h"
@@ -41,9 +45,48 @@ static double z0, cx, cy; /* z = cx * x + cy * y + z0 */
static void slice_plane(void)
{
- /* @@@ */
- cx = cy = 0;
- z0 = points->z;
+ const struct point *p;
+ unsigned n = 0;
+ gsl_multifit_robust_workspace *work;
+ gsl_vector *y, *c;
+ gsl_matrix *x, *cov;
+
+ for (p = points; p; p = p->next)
+ n++;
+
+ if (n < 3) {
+ cx = cy = 0;
+ z0 = points->z;
+ return;
+ }
+
+ x = gsl_matrix_alloc(n, 3);
+ y = gsl_vector_alloc(n);
+ n = 0;
+ for (p = points; p; p = p->next) {
+ gsl_vector_set(y, n, p->z);
+ gsl_matrix_set(x, n, 0, p->x);
+ gsl_matrix_set(x, n, 1, p->y);
+ gsl_matrix_set(x, n, 2, 1);
+ n++;
+ }
+
+ c = gsl_vector_alloc(3);
+ cov = gsl_matrix_alloc(3, 3);
+
+ work = gsl_multifit_robust_alloc(gsl_multifit_robust_default, n, 3);
+ gsl_multifit_robust(x, y, c, cov, work);
+ gsl_multifit_robust_free(work);
+
+ cx = gsl_vector_get(c, 0);
+ cy = gsl_vector_get(c, 1);
+ z0 = gsl_vector_get(c, 2);
+
+ gsl_matrix_free(x);
+ gsl_vector_free(y);
+ gsl_vector_free(c);
+
+fprintf(stderr, "cx %g cy %g z0 %g\n", cx, cy, z0);
}