[GRASS-dev] patch for v.voronoi/v.delaunay

[Benjamin Ducke]

Hi all,

As per request, I have patched v.delaunay (source in v.voronoi/)
very carefully.

The memory handling now seems to work better. At least I was able
to triangulate a map with ca 30000 points. Anything more and my
system starts thrashing badly (1 GB, with a lot of other apps running).
Feel free to test with bigger inputs.

The output vector map will now also preserve the 3D coordinates of
the original input points. However, the centroids for the triangles
are all wrong (placed at z=0.0). I don't see a quick way of fixing
this right now.

In fact, I'd prefer not to touch the current code for v.delaunay
anymore. It is in urgent need of a complete overhaul as it is
inefficient, hard to read and still full of site lib artefacts.
Better to re-write from scratch with more time on my (or someone
else's) hands.

Best,

Benjamin

-- 
Benjamin Ducke, M.A.
Archäoinformatik
(Archaeoinformation Science)
Institut für Ur- und Frühgeschichte
(Inst. of Prehistoric and Historic Archaeology)
Christian-Albrechts-Universität zu Kiel
Johanna-Mestorf-Straße 2-6
D 24098 Kiel
Germany

Tel.: ++49 (0)431 880-3378 / -3379
Fax : ++49 (0)431 880-7300
www.uni-kiel.de/ufg

-------------- next part --------------
diff -Naur old/dt_main.c new/dt_main.c
--- old/dt_main.c	2008-02-01 13:04:42.000000000 +0100
+++ new/dt_main.c	2008-02-01 13:05:08.000000000 +0100
@@ -74,12 +74,26 @@
   if ((mapset = G_find_vector2 (in_opt->answer, "")) == NULL) {
       G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
   }
-
+  
   Vect_set_open_level (2);
   Vect_open_old (&In, in_opt->answer, mapset);
 
-  if (0 > Vect_open_new (&Out, out_opt->answer, 0)) {
-    G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
+  /* check if we have a 3D input points map */
+  mode3d = 0;
+  if ( Vect_is_3d ( &In ) ) {
+    	mode3d = 1;
+  }
+
+
+  if ( mode3d ) {
+ 	if (0 > Vect_open_new (&Out, out_opt->answer, 1)) {
+    		G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
+  	}
+  } else {
+ 	if (0 > Vect_open_new (&Out, out_opt->answer, 0)) {
+    		G_fatal_error(_("Unable to create vector map <%s>"), out_opt->answer);
+  	}
+  
   }
 
   Vect_hist_copy (&In, &Out);
diff -Naur old/dt_write.c new/dt_write.c
--- old/dt_write.c	2008-02-01 13:04:42.000000000 +0100
+++ new/dt_write.c	2008-02-01 13:05:08.000000000 +0100
@@ -73,8 +73,13 @@
 
 	    /* Not found, write it */
 	    Vect_reset_line ( Points );
-	    Vect_append_point ( Points, sa->coord.x, sa->coord.y, 0.0 );
-	    Vect_append_point ( Points, sb->coord.x, sb->coord.y, 0.0 );
+	    if ( mode3d ) {
+	    	Vect_append_point ( Points, sa->coord.x, sa->coord.y, sa->coord.z );
+	    	Vect_append_point ( Points, sb->coord.x, sb->coord.y, sb->coord.z );	    
+	    } else {
+	    	Vect_append_point ( Points, sa->coord.x, sa->coord.y, 0.0 );
+	    	Vect_append_point ( Points, sb->coord.x, sb->coord.y, 0.0 );
+	    }
 	    Vect_write_line ( &Out, Type, Points, Cats );
 	}
     }
diff -Naur old/sw_defs.h new/sw_defs.h
--- old/sw_defs.h	2008-02-01 13:04:42.000000000 +0100
+++ new/sw_defs.h	2008-02-01 13:05:08.000000000 +0100
@@ -13,7 +13,7 @@
 
 
 struct Point	{
-    double x,y;
+    double x,y,z;
 };
 
 /* structure used both for sites and for vertices */
@@ -41,7 +41,7 @@
 };
 
 #ifdef MAIN
-int triangulate, sorted, plot, debug;
+int triangulate, sorted, plot, debug, mode3d;
 struct	Site	*sites;
 int		nsites;
 int		siteidx;
@@ -61,7 +61,7 @@
 int PQcount;
 int PQmin;
 #else
-extern int triangulate, sorted, plot, debug;
+extern int triangulate, sorted, plot, debug, mode3d;
 extern struct	Site	*sites;
 extern int		nsites;
 extern int		siteidx;
diff -Naur old/sw_main.c new/sw_main.c
--- old/sw_main.c	2008-02-01 13:04:42.000000000 +0100
+++ new/sw_main.c	2008-02-01 13:05:08.000000000 +0100
@@ -42,12 +42,23 @@
 	i = j = 1;
 	foundDupe = 0;
 	while(i < nsites)
-		if(sites[i].coord.x == sites[i-1].coord.x && sites[i].coord.y == sites[i-1].coord.y)
-			i++;
-		else
-		{
-			if(i != j) sites[j] = sites[i];
-			i++; j++;;
+		if ( mode3d ) {
+			if(sites[i].coord.x == sites[i-1].coord.x && sites[i].coord.y == sites[i-1].coord.y 
+								  && sites[i].coord.z == sites[i-1].coord.z)
+				i++;
+			else
+			{
+				if(i != j) sites[j] = sites[i];
+				i++; j++;;
+			}		
+		} else {
+			if(sites[i].coord.x == sites[i-1].coord.x && sites[i].coord.y == sites[i-1].coord.y)
+				i++;
+			else
+			{
+				if(i != j) sites[j] = sites[i];
+				i++; j++;;
+			}
 		}
 
 	if(j != nsites)
@@ -67,11 +78,11 @@
 
     Points = Vect_new_line_struct ();
 
-    nsites = 0;
-    sites = (struct Site *) myalloc(4000*sizeof(*sites));
-
     nlines = Vect_get_num_lines ( &In );
 
+    nsites = 0;
+    sites = (struct Site *) myalloc(nlines * sizeof(*sites));
+
     for ( line = 1; line <= nlines; line++ ) {
 	int type;
 
@@ -84,6 +95,9 @@
 	
 	sites[nsites].coord.x = Points->x[0];
 	sites[nsites].coord.y = Points->y[0];
+	if ( mode3d ) {
+		sites[nsites].coord.z = Points->z[0];
+	}
 
 	sites[nsites].sitenbr = nsites;
 	sites[nsites].refcnt = 0;