[GRASS-SVN] r64538 - in grass-addons/grass7/vector: . v.ellipse
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Feb 10 14:06:23 PST 2015
Author: martinl
Date: 2015-02-10 14:06:23 -0800 (Tue, 10 Feb 2015)
New Revision: 64538
Added:
grass-addons/grass7/vector/v.ellipse/
grass-addons/grass7/vector/v.ellipse/Makefile
grass-addons/grass7/vector/v.ellipse/ellipse.c
grass-addons/grass7/vector/v.ellipse/main.c
grass-addons/grass7/vector/v.ellipse/proto.h
grass-addons/grass7/vector/v.ellipse/read.c
grass-addons/grass7/vector/v.ellipse/v.ellipse.html
grass-addons/grass7/vector/v.ellipse/v_ellipse.png
grass-addons/grass7/vector/v.ellipse/write.c
Log:
v.ellipse: Computes the best-fitting ellipse for given vector data
(author: Tereza Fiedlerova, OSGeoREL, CTU in Prague)
Added: grass-addons/grass7/vector/v.ellipse/Makefile
===================================================================
--- grass-addons/grass7/vector/v.ellipse/Makefile (rev 0)
+++ grass-addons/grass7/vector/v.ellipse/Makefile 2015-02-10 22:06:23 UTC (rev 64538)
@@ -0,0 +1,14 @@
+# fix this relative to include/
+# or use absolute path to the GRASS source code
+MODULE_TOPDIR = ../..
+
+PGM = v.ellipse
+
+LIBES = $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB) $(GMATHLIB)
+DEPENDENCIES = $(VECTORDEP) $(DBMIDEP) $(GISDEP) $(GMATHDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
Property changes on: grass-addons/grass7/vector/v.ellipse/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.ellipse/ellipse.c
===================================================================
--- grass-addons/grass7/vector/v.ellipse/ellipse.c (rev 0)
+++ grass-addons/grass7/vector/v.ellipse/ellipse.c 2015-02-10 22:06:23 UTC (rev 64538)
@@ -0,0 +1,220 @@
+#include <math.h>
+
+#include <grass/vector.h>
+#include <grass/gmath.h>
+#include <grass/glocale.h>
+
+#include "proto.h"
+
+void create_ellipse(struct Parameters *par, struct line_pnts *Ellipse,
+ struct line_cats *Cats, double stepsize)
+{
+ Vect_reset_line(Ellipse);
+
+ /* angle */
+ double t = 0;
+
+ /* number of points */
+ if (stepsize < 0.01) /* too small step */
+ stepsize = 0.01;
+ int n = 360.0 / stepsize;
+
+ /* compute n points on ellipse */
+ double *x;
+ double *y;
+
+ x = (double *)G_malloc((n + 1) * sizeof(double));
+ y = (double *)G_malloc((n + 1) * sizeof(double));
+
+ int i;
+
+ for (i = 0; i < n; i++) {
+ x[i] = par->a * cos(t) * cos(par->angle) - par->b * sin(t) * sin(par->angle) + par->xc; /* parametric equations of an ellipse (non-tilted) */
+ y[i] =
+ par->a * cos(t) * sin(par->angle) +
+ par->b * sin(t) * cos(par->angle) + par->yc;
+ Vect_cat_set(Cats, 1, i);
+ t += 2 * M_PI / n;
+ }
+
+ x[n] = x[0];
+ y[n] = y[0];
+ Vect_cat_set(Cats, 1, n);
+
+ Vect_copy_xyz_to_pnts(Ellipse, x, y, 0, n + 1);
+
+ G_free(x);
+ G_free(y);
+}
+
+int fitting_ellipse(struct line_pnts *Points, struct Parameters *pars)
+{
+ /* variables */
+ int i, j, k, n;
+ double meanx, meany;
+ double *x, *y, *nn;
+ double **AA, **NN;
+ int singular;
+
+ /* LSM for conic representation of an ellipse */
+ /* a*x^2 + b*x*y + c*y^2 + d*x + e*y = f */
+
+ /* allocate memory */
+ n = Points->n_points;
+ x = G_alloc_vector(n); /* x-coordinates */
+ y = G_alloc_vector(n); /* y-coordinates */
+ AA = G_alloc_matrix(n, 5); /* design matrix [x.^2, y.^2, x.*y, y.^2, x, y] */
+ NN = G_alloc_matrix(5, 5); /* N = A'*A */
+ nn = G_alloc_vector(5); /* n = A'*l, vector of unknown
+ * parameters of conic representation of
+ * ellipse [a, b, c, d, e] */
+
+ /* remove bias of the ellipse (to be more accurate for inversion) */
+ /* compute mean */
+ meanx = 0;
+ meany = 0;
+ /* TODO: do this more effective */
+ for (i = 0; i < n; i++) {
+ meanx += Points->x[i] / n;
+ meany += Points->y[i] / n;
+ }
+ /* change coordinates */
+ for (i = 0; i < n; i++) {
+ x[i] = Points->x[i] - meanx;
+ y[i] = Points->y[i] - meany;
+ }
+
+ /* compute design matrix */
+ for (i = 0; i < n; i++) {
+ AA[i][0] = x[i] * x[i];
+ AA[i][1] = x[i] * y[i];
+ AA[i][2] = y[i] * y[i];
+ AA[i][3] = x[i];
+ AA[i][4] = y[i];
+ }
+
+ /* compute N=A'*A */
+ /* TODO: use fn - mat_struct *G_matrix_product (mat_struct *mt1, mat_struct *mt2) */
+ for (i = 0; i < 5; i++) {
+ for (j = 0; j < 5; j++) {
+ NN[i][j] = 0;
+
+ if (i > j) {
+ NN[i][j] = NN[j][i]; // matrix is symmetric
+ }
+ else {
+ for (k = 0; k < n; k++) {
+ NN[i][j] += AA[k][i] * AA[k][j];
+ }
+ }
+ }
+ }
+
+ /* compute n = A'*l (as column sum) */
+ for (i = 0; i < 5; i++) {
+ nn[i] = 0;
+ for (k = 0; k < n; k++) {
+ nn[i] += AA[k][i];
+ }
+ }
+
+ /* solve x = N^(-1)*n */
+ singular = G_math_solv(NN, nn, 5); /* TODO: ensure use of the most effective method from G_math_sol... */
+
+ /* free memory */
+ G_free_vector(x);
+ G_free_vector(y);
+ G_free_matrix(AA);
+ G_free_matrix(NN);
+
+ if (singular == -1) {
+ G_free_vector(nn);
+ G_warning(_("Singular matrix failure"));
+ return 0;
+ }
+
+ /* extract parameters of ellipse from conic representation */
+
+ double a, b, c, d, e, f;
+ double test, tol, cos_phi, sin_phi, x0, y0, sqrta, sqrtb; /*, sqrta, sqrtb; */
+
+ tol = 1e-3;
+
+ /* parameters of conic representation */
+ a = nn[0];
+ b = nn[1];
+ c = nn[2];
+ d = nn[3];
+ e = nn[4];
+ /* f = -1; */
+
+ G_free_vector(nn);
+
+ /* check if it is an ellipse */
+ test = a * c;
+ if (test == 0) {
+ G_warning(_("Parabola found. Exiting."));
+ return 0;
+ }
+ else if (test < 0) {
+ G_warning(_("Hyperbola found. Exiting."));
+ return 0;
+ }
+
+ /* orientation (and its removing for further computing) */
+ if ((abs(b / a) > tol) && (abs(b / c) > tol)) {
+ pars->angle = 0.5 * atan2(b, (c - a));
+ cos_phi = cos(pars->angle);
+ sin_phi = sin(pars->angle);
+ a = a * cos_phi * cos_phi - b * cos_phi * sin_phi +
+ c * sin_phi * sin_phi;
+ b = 0;
+ c = a * sin_phi * sin_phi - b * cos_phi * sin_phi +
+ c * cos_phi * cos_phi;
+ d = d * cos_phi - e * sin_phi;
+ e = d * sin_phi + e * cos_phi;
+ }
+ else {
+ pars->angle = 0;
+ cos_phi = 1;
+ sin_phi = 0;
+ }
+
+ /* make sure coefficients are positive */
+ if (a < 0) {
+ a *= -1;
+ b *= -1;
+ c *= -1;
+ d *= -1;
+ e *= -1;
+ }
+
+ /* center of non-tilted ellipse */
+ x0 = meanx - d / 2 / a;
+ y0 = meany - e / 2 / c;
+ /* center of tilted ellipse */
+ pars->xc = cos_phi * x0 + sin_phi * y0;
+ pars->yc = cos_phi * y0 - sin_phi * x0;
+
+ /* semi-axis */
+ f = 1 + d * d / 4 / a + e * e / 4 / c;
+ sqrta = sqrt(f / a);
+ sqrtb = sqrt(f / c);
+ if (sqrta > sqrtb) {
+ pars->a = sqrta;
+ pars->b = sqrtb;
+ }
+ else {
+ pars->b = sqrta;
+ pars->a = sqrtb;
+ }
+
+ G_verbose_message("%s: %f*x^2+%f*x*y+%f*y^2+%f*x+%f*y-1=0",
+ _("conic representation"), a, b, c, d, e);
+ G_verbose_message
+ ("parameters of ellipse:\na = %f \nb = %f \nangle = %f \n"
+ "xc = %f \nyc = %f", pars->a, pars->b, pars->angle, pars->xc,
+ pars->yc);
+
+ return 1;
+}
Property changes on: grass-addons/grass7/vector/v.ellipse/ellipse.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.ellipse/main.c
===================================================================
--- grass-addons/grass7/vector/v.ellipse/main.c (rev 0)
+++ grass-addons/grass7/vector/v.ellipse/main.c 2015-02-10 22:06:23 UTC (rev 64538)
@@ -0,0 +1,109 @@
+/***************************************************************
+ *
+ * MODULE: v.ellipse
+ *
+ * AUTHOR(S): Tereza Fiedlerova
+ *
+ * PURPOSE: Computes the best-fitting ellipse for
+ * given vector data
+ *
+ * COPYRIGHT: (C) 2015 by Tereza Fiedlerova, and the GRASS Development Team
+ *
+ * This program is free software under the GNU General
+ * Public License (>=v2). Read the file COPYING that
+ * comes with GRASS for details.
+ **************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include <grass/glocale.h>
+
+#include "proto.h"
+
+int main(int argc, char *argv[])
+{
+ /* variables */
+ struct Map_info In, Out;
+ static struct line_pnts *Points, *Ellipse;
+ struct line_cats *Cats, *Cats2;
+ struct GModule *module;
+ struct Option *old, *new, *stepopt;
+ struct Parameters *pars;
+ double stepsize;
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]);
+
+ /* initialize module */
+ module = G_define_module();
+ G_add_keyword(_("vector"));
+ G_add_keyword(_("geometry"));
+ G_add_keyword(_("best-fitting ellipse"));
+ module->description =
+ _("Computes the best-fitting ellipse for given vector data.");
+
+ /* define options */
+ old = G_define_standard_option(G_OPT_V_INPUT);
+ new = G_define_standard_option(G_OPT_V_OUTPUT);
+
+ stepopt = G_define_option();
+ stepopt->type = TYPE_DOUBLE;
+ stepopt->key = "step";
+ stepopt->description = "Step size in degrees";
+ stepopt->required = NO;
+ stepopt->answer = "4";
+
+ /* call parser */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ /* read step size */
+ stepsize = atof(stepopt->answer); /* conversion to double stepopt->answer */
+
+ G_debug(1, "step size: %f degrees", stepsize);
+
+ /* test if input and output are different */
+ Vect_check_input_output_name(old->answer, new->answer, G_FATAL_EXIT);
+
+ /* open input vector map */
+ Vect_set_open_level(2);
+ if (1 > Vect_open_old(&In, old->answer, ""))
+ G_fatal_error(_("Unable to open vector map <%s> on topological level"),
+ old->answer);
+
+ /* open output vector map */
+ if (0 > Vect_open_new(&Out, new->answer, WITHOUT_Z)) {
+ Vect_close(&In);
+ G_fatal_error(_("Unable to create vector map <%s>"), new->answer);
+ }
+ Vect_set_error_handler_io(&In, &Out);
+
+ /* create and initializes structs where to store points, categories */
+ Cats = Vect_new_cats_struct();
+ Cats2 = Vect_new_cats_struct();
+ Points = Vect_new_line_struct();
+ Ellipse = Vect_new_line_struct();
+
+ /* allocate memory */
+ pars = (struct Parameters *)G_malloc(sizeof(struct Parameters));
+
+ int ok;
+
+ load_points(&In, Points, Cats);
+ ok = fitting_ellipse(Points, pars);
+ G_debug(1, "exit %d", ok);
+ if (ok) {
+ create_ellipse(pars, Ellipse, Cats2, stepsize);
+ write_ellipse(&Out, Ellipse, Cats2);
+ }
+
+ Vect_destroy_line_struct(Points);
+
+ G_free(pars);
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/grass7/vector/v.ellipse/main.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.ellipse/proto.h
===================================================================
--- grass-addons/grass7/vector/v.ellipse/proto.h (rev 0)
+++ grass-addons/grass7/vector/v.ellipse/proto.h 2015-02-10 22:06:23 UTC (rev 64538)
@@ -0,0 +1,24 @@
+#ifndef __V_ELLIPSE_PROTO__
+#define __V_ELLIPSE_PROTO__
+
+/* struct with ellipse parameters */
+struct Parameters
+{
+ double a; /* main semi-axis */
+ double b; /* second semi-axis */
+ double xc; /* x of center of ellipse */
+ double yc; /* y of center of ellipse */
+ double angle; /* rotation */
+};
+
+/* read.c */
+void load_points(struct Map_info *, struct line_pnts *, struct line_cats *);
+
+/* write.c */
+void write_ellipse(struct Map_info *, struct line_pnts *, struct line_cats *);
+
+/* ellipse.c */
+void create_ellipse(struct Parameters *, struct line_pnts *, struct line_cats *, double);
+int fitting_ellipse(struct line_pnts *, struct Parameters *);
+
+#endif
Property changes on: grass-addons/grass7/vector/v.ellipse/proto.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.ellipse/read.c
===================================================================
--- grass-addons/grass7/vector/v.ellipse/read.c (rev 0)
+++ grass-addons/grass7/vector/v.ellipse/read.c 2015-02-10 22:06:23 UTC (rev 64538)
@@ -0,0 +1,34 @@
+#include <grass/vector.h>
+
+#include "proto.h"
+
+void load_points(struct Map_info *In, struct line_pnts *Points,
+ struct line_cats *Cats)
+{
+ /* variables */
+ struct line_pnts *OPoints;
+ int type, i, cat;
+
+ /* create and initializes structs where to store points, categories */
+ OPoints = Vect_new_line_struct();
+
+ Vect_reset_cats(Cats);
+
+ /* load points */
+ while ((type = Vect_read_next_line(In, OPoints, Cats)) > 0) {
+ if (type == GV_LINE || type == GV_POINT || type == GV_CENTROID) {
+ if (Vect_cat_get(Cats, 1, &cat) == 0) {
+ Vect_cat_set(Cats, 1, i);
+ i++;
+ }
+ }
+
+ Vect_append_points(Points, OPoints, GV_FORWARD);
+ Vect_reset_line(OPoints);
+ }
+
+ Vect_destroy_line_struct(OPoints);
+
+ /* close map */
+ Vect_close(In);
+}
Property changes on: grass-addons/grass7/vector/v.ellipse/read.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.ellipse/v.ellipse.html
===================================================================
--- grass-addons/grass7/vector/v.ellipse/v.ellipse.html (rev 0)
+++ grass-addons/grass7/vector/v.ellipse/v.ellipse.html 2015-02-10 22:06:23 UTC (rev 64538)
@@ -0,0 +1,45 @@
+<h2>DESCRIPTION</h2>
+
+<em>v.ellipse</em> computes the best-fitting ellipse for <b>input</b>
+vector map and creates new <b>output</b> vector map with
+ellipse. Input vector data might be 2D points, lines, or areas.
+
+The parameters of ellipse are printed on output if <b>--verbose</b>
+flag is given.
+
+<center>
+<img src="v_ellipse.png" border="1"><br>
+<i>Fig: Fitting ellipse created with v.ellipse</i>
+</center>
+
+<h2>EXAMPLE</h2>
+
+Example of <em>v.ellipse</em> created around set of points (using data
+<i>points_of_interest</i>, North Carolina sample data set). Ellipse is
+is approximated by linestring with point distance 1 degree
+(<b>step</b>).
+
+<div class="code"><pre>
+v.ellipse input=points_of_interest output=ellipse step=1
+</pre></div>
+
+<h2>REFERENCES</h2>
+
+<ul>
+ <li><a href="https://www.cs.cornell.edu/cv/OtherPdf/Ellipse.pdf">Charles
+ F. Van Loan: Using the Ellipse to Fit and Enclose Data
+ Points.</a></li>
+</ul>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="v.hull.html">v.hull</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+Tereza Fiedlerova, OSGeoREL, Czech Technical University in Prague
+
+<p>
+<i>Last changed: $Date$</i>
Property changes on: grass-addons/grass7/vector/v.ellipse/v.ellipse.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/vector/v.ellipse/v_ellipse.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/vector/v.ellipse/v_ellipse.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: grass-addons/grass7/vector/v.ellipse/write.c
===================================================================
--- grass-addons/grass7/vector/v.ellipse/write.c (rev 0)
+++ grass-addons/grass7/vector/v.ellipse/write.c 2015-02-10 22:06:23 UTC (rev 64538)
@@ -0,0 +1,36 @@
+#include <grass/vector.h>
+
+#include "proto.h"
+
+void write_ellipse(struct Map_info *Out, struct line_pnts *Ellipse,
+ struct line_cats *Cats)
+{
+ /* variables */
+ double xc, yc;
+ int i, n;
+
+ /* write boundary */
+ Vect_write_line(Out, GV_BOUNDARY, Ellipse, Cats);
+
+ /* write centroid */
+ xc = 0;
+ yc = 0;
+ n = Ellipse->n_points;
+
+ for (i = 0; i < n; i++) {
+ xc += Ellipse->x[i];
+ yc += Ellipse->y[i];
+ }
+
+ xc /= n;
+ yc /= n;
+ Vect_reset_line(Ellipse);
+ Vect_append_point(Ellipse, xc, yc, 0);
+ Vect_cat_set(Cats, 1, 1);
+ Vect_write_line(Out, GV_CENTROID, Ellipse, Cats);
+ Vect_destroy_line_struct(Ellipse);
+
+ /* build topology and close map */
+ Vect_build(Out);
+ Vect_close(Out);
+}
Property changes on: grass-addons/grass7/vector/v.ellipse/write.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list