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

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

(attachments)

v.voronoi.diff (3.75 KB)

Hi Benjamin,
I'm sorry to disappoint You - Your patch does not fix bogous output
from v.delaunay - there's still a garbage output from 46647 input
points. Look at this old bug how to reproduce problem on 32bit system
[1]. Your changes in calling myalloc have no effect, as previously
allocated memory was increased by 4000*sizeof(*sites) chunks and now
is assigned in one shot at beginning. It gives only slight memory saving.

I understand Your (and others) unwillingness to fix this code, still
shipping broken code is also not a good thing. It MUST be
fixed (patch/rewrite) or there must be some additional checks like
if(npoints>MAX_INT) G_error("Blah! Decrease input vector size") to avoid producing bogous results without any warnings.

Maris.

1. http://wald.intevation.org/tracker/index.php?func=detail&aid=335&group_id=21&atid=204

2008/2/1, Benjamin Ducke <benjamin.ducke@ufg.uni-kiel.de>:

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

So maybe the thing to do for now is to mention all these issues
in the HTML man page. The "BUGS" section is silent about them
right now.
Does the v.voronoi part suffer from the same problems?

Benjamin

Maris Nartiss wrote:

Hi Benjamin,
I'm sorry to disappoint You - Your patch does not fix bogous output
from v.delaunay - there's still a garbage output from 46647 input
points. Look at this old bug how to reproduce problem on 32bit system
[1]. Your changes in calling myalloc have no effect, as previously
allocated memory was increased by 4000*sizeof(*sites) chunks and now
is assigned in one shot at beginning. It gives only slight memory saving.

I understand Your (and others) unwillingness to fix this code, still
shipping broken code is also not a good thing. It MUST be
fixed (patch/rewrite) or there must be some additional checks like
if(npoints>MAX_INT) G_error("Blah! Decrease input vector size") to avoid producing bogous results without any warnings.

Maris.

1. http://wald.intevation.org/tracker/index.php?func=detail&aid=335&group_id=21&atid=204

2008/2/1, Benjamin Ducke <benjamin.ducke@ufg.uni-kiel.de>:

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

--
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

I understand Your (and others) unwillingness to fix this code, still
shipping broken code is also not a good thing. It MUST be
fixed (patch/rewrite) or there must be some additional checks like
if(npoints>MAX_INT) G_error("Blah! Decrease input vector size") to avoid producing bogous results without any warnings.

I am not really unwilling to fix this. In fact, I am looking forward to
re-writing the delaunay and voronoi code with full 3D support, as this
would be something to really set GRASS GIS apart.
However, I simply don't have the time to start on this before spring
this year.

So the question is: what are we going to do as a quick fix for the 6.3
release?

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

Hi Benjamin,

On Fri, Feb 1, 2008 at 2:16 PM, Benjamin Ducke
<benjamin.ducke@ufg.uni-kiel.de> wrote:

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.

I tried with nearly 40k points and it worked.

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.

I have fixed it using Vect_tin_get_z() for delaunay (someone needs to
do it for voronoi).

I have submitted the fixed code to SVN:
http://trac.osgeo.org/grass/changeset/31079

For me it now generates a 3D Delaunay TIN.

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.

Forthcoming...:
http://code.google.com/soc/2008/osgeo/appinfo.html?csaid=7599B5B1E8D1F20F

Markus

On Wed, Apr 23, 2008 at 1:07 AM, Markus Neteler <neteler@osgeo.org> wrote:

On Fri, Feb 1, 2008 at 2:16 PM, Benjamin Ducke wrote:

...

> 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.

I have fixed it using Vect_tin_get_z() for delaunay (someone needs to
do it for voronoi).

I have submitted the fixed code to SVN:
http://trac.osgeo.org/grass/changeset/31079

For me it now generates a 3D Delaunay TIN.

I wonder if I could colorize the 3D TIN somehow using v.colors.
So far it does not have a DB connected.

Markus

Markus Neteler wrote:

I wonder if I could colorize the 3D TIN somehow using v.colors.
So far it does not have a DB connected.

if not already, make it into a boundary+centroid map with unique cats for
each centroid (v.type, v.centroids / v.category)

create a DB table for it with columns="int cat, GRASSRGB varchar(11)" and
add new columns for double z1,z2,z3,z_mean for elevation at each triangle
corner and the average of the three.

somehow for each triangle extract the 3 elev values (v.extract +
v.out.ascii in a loop would work but probably there is a better way,
maybe v.to.db or like v.what.vect). do either SQL magic or awk magic or
..?. to get an average value for the triangle

upload the mean to the centroid with v.db.update or similar, then run
v.colors..

Hamish

      ____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile. Try it now. http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ

On Wed, Apr 23, 2008 at 8:14 AM, Hamish <hamish_b@yahoo.com> wrote:

Markus Neteler wrote:
> I wonder if I could colorize the 3D TIN somehow using v.colors.
> So far it does not have a DB connected.

...

create a DB table for it with columns="int cat, GRASSRGB varchar(11)" and

v.db.addtable tin col="east double precision, north double precision, \
      height double precision, GRASSRGB varchar(11)"

[... not needed ... ]

somehow for each triangle extract the 3 elev values (v.extract +
v.out.ascii in a loop would work but probably there is a better way,
maybe v.to.db or like v.what.vect). do either SQL magic or awk magic or
..?. to get an average value for the triangle

That's now automatically done with the recent 3D extentions to v.delaunay
(the Vect_tin_get_z() was sleeping for some years, now used).

upload the mean to the centroid with v.db.update or similar, then run
v.colors..

# transfer geometry for colorizing
v.to.db tin opt=coor col="east,north,height"
v.db.select tin

v.colors tin column=height rgb_column=GRASSRGB color=rainbow

Should work... but it fails due to a precision problem in the output
of v.db.select:

v.colors tin column=height rgb_column=GRASSRGB color=rainbow
...
Looking up colors ...
+ v.db.select map=tin layer=1 column=height
+ sort -n
+ grep '^[-0-9]'
+ uniq
+ r.what.color -i in=tmp_colr_22337
+ sed -e 's/: /|/'
+ grep -v '|\*$'
+ read LINE
++ echo '1.353553|255:255:0'
++ cut -f1 '-d|'
+ VALUE=1.353553
++ echo '1.353553|255:255:0'
++ cut -f2 '-d|'
+ COLR=255:255:0
+ echo 'UPDATE tin SET GRASSRGB = '\''255:255:0'\'' WHERE height = 1.353553;'
...

[ for comparison:
v.db.select output:
...
D3/3: fetch row = 4
D3/3: col 0, litetype 2, sqltype 6: val = '4.0'
D3/3: Row fetched
4
D3/3: fetch row = 5
D3/3: col 0, litetype 2, sqltype 6: val = '2.35355339059327'
D3/3: Row fetched
2.353553
...
]

sqlite3 /home/neteler/grassdata/eth_utm32/PERMANENT/sqlite.db
SQLite version 3.4.2
Enter ".help" for instructions
sqlite> select * from tin;
1|3.53553390593274|6.46446609406726|1.35355339059327|
2|6.46446609406726|3.53553390593274|1.35355339059327|
3|13.5355339059327|6.46446609406726|2.64644660940673|
4|6.46446609406726|13.5355339059327|2.64644660940673|
5|16.4644660940673|3.53553390593274|2.06066017177982|
6|16.4644660940673|13.5355339059327|4.0|255:0:1
7|3.53553390593274|16.4644660940673|2.35355339059327|
8|13.5355339059327|16.4644660940673|4.0|255:0:1
sqlite>

Should we fix the v.db.select output (new flag for full precision output)?

Markus

On Fri, May 2, 2008 at 11:15 PM, Markus Neteler <neteler@osgeo.org> wrote:
...

v.db.select output:
...
D3/3: fetch row = 4
D3/3: col 0, litetype 2, sqltype 6: val = '4.0'
D3/3: Row fetched
4
D3/3: fetch row = 5
D3/3: col 0, litetype 2, sqltype 6: val = '2.35355339059327'
D3/3: Row fetched
2.353553

It seems that the precision gets somewhere lost in
db_convert_column_value_to_string()

i.e. lib/db/dbmi_base/columnfmt.c

but I don't find where is is rounded.

Markus

Markus Neteler wrote:

It seems that the precision gets somewhere lost in
db_convert_column_value_to_string()

i.e. lib/db/dbmi_base/columnfmt.c

but I don't find where is is rounded.

lib/db/dbmi_base/valuefmt.c db_convert_value_to_string() does:
        case DB_C_TYPE_DOUBLE:
            sprintf (buf, "%lf",db_get_value_double(value));
            G_trim_decimal(buf);

the %lf is rounding to %.6f?

maybe
- "%lf"
+ "%.14f"

then let G_trim_decimal() remove any extra .00000s.

but maybe such noisy output is not always wanted. :-/
(e.g. pretty output for d.what.vect tcl form)
?

FWIW v.out.ascii defaults to dp=8.

for the v.colors script, SQL "select * WHERE dbl_col = <FP number>" is
always going to be problematic.

the C way:
if( abs(test - value) < GRASS_EPSILON )

Hamish

____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile. Try it now. http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ

      ____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile. Try it now. http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ

On Sat, May 3, 2008 at 5:33 AM, Hamish <hamish_b@yahoo.com> wrote:

Markus Neteler wrote:
> It seems that the precision gets somewhere lost in
> db_convert_column_value_to_string()
>
> i.e. lib/db/dbmi_base/columnfmt.c
>
> but I don't find where is is rounded.

lib/db/dbmi_base/valuefmt.c db_convert_value_to_string() does:
        case DB_C_TYPE_DOUBLE:
            sprintf (buf, "%lf",db_get_value_double(value));
            G_trim_decimal(buf);

the %lf is rounding to %.6f?

maybe
- "%lf"
+ "%.14f"

then let G_trim_decimal() remove any extra .00000s.

but maybe such noisy output is not always wanted. :-/

Let's do the fix because a library should keep the precision.
We can reduce noise in the modules then.

(e.g. pretty output for d.what.vect tcl form)
?

Yes, but trim that at module level.

FWIW v.out.ascii defaults to dp=8.

So we should consequently use the same length everywhere
(adjustable) and full precision at library level.

for the v.colors script, SQL "select * WHERE dbl_col = <FP number>" is
always going to be problematic.

the C way:
  if( abs(test - value) < GRASS_EPSILON )

On Sat, May 3, 2008 at 7:41 AM, Hamish <hamish_b@yahoo.com> wrote:

In the case of v.colors, the SQL set loop could be rewritten to avoid
this problem by matching to the cat number, not trying to logically match
the FP value.

This would be good.

Markus

> Markus wrote:
> > It seems that the precision gets somewhere lost in
> > db_convert_column_value_to_string()
> >
> > i.e. lib/db/dbmi_base/columnfmt.c
> >
> > but I don't find where is is rounded.
>

Hamish:

> lib/db/dbmi_base/valuefmt.c db_convert_value_to_string() does:
> case DB_C_TYPE_DOUBLE:
> sprintf (buf, "%lf",db_get_value_double(value));
> G_trim_decimal(buf);
>
> the %lf is rounding to %.6f?
>
> maybe
> - "%lf"
> + "%.14f"
>
> then let G_trim_decimal() remove any extra .00000s.
>
> but maybe such noisy output is not always wanted. :-/

Let's do the fix because a library should keep the precision.
We can reduce noise in the modules then.

understand that you will then change the way the DB attribute library is
currently designed, it's not as simple as a small patch or new lib fn
(beyond adjusting %lf to show more digits). For me, I don't feel
qualified to mess around with Radim's DB lib code much.

output from computer stored IEEE FP representation to human-readable will
always need some form of %f simplification, and associated loss of
exactness. Life is a compromise, as they say.

[v.]db.select is the important one, as that is specifically data out.

> (e.g. pretty output for d.what.vect tcl form)
> ?

Yes, but trim that at module level.

I think it might be annoying for every module to have their own random
%.[6-16]f choice, I think it is better to write it into one place, or
have two lib fns, one standard and the other full precision. But having
two would be a bit annoying and doesn't address the real problem.

> FWIW v.out.ascii defaults to dp=8.

So we should consequently use the same length everywhere
(adjustable) and full precision at library level.

I expect SQL "a * b" works at the DB level, and so doesn't suffer the
loss of precision?

> In the case of v.colors, the SQL set loop could be rewritten to
> avoid this problem by matching to the cat number, not trying to
> logically match the FP value.

This would be good.

Yes. I feel some responsibility to fix that, but unfortunately I have
very little time to work on that currently.

Hamish

____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile. Try it now. http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ

      ____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile. Try it now. http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ