[GRASS5] Voxel progress

[benducke at compuserve.de]

As I hear, GRASS-based voxel modelling and its use in
archaeology were a topic at the last GRASS user meeting.
It certainly is a hot topic in archaeology and I am trying
to write some GRASS modules to this end.

It seems that volume visualization in NVIZ still needs some bugfixing.
I would like to use GRASS voxels to represent stratigraphic layers.
For this, I am working on a module that builds a voxel model from a series
of DEMs (each representing a stratigraphic layer surface). 
It uses a sort of 3D flood filling algorithm.

Each DEM is taken as the upper surface of 0 to n input layers.
The module will fill down voxels from each surface and assign each
voxel the number of that surface until it hits the next, underlying surface.

So far, things work fine and if I use r3.out.vtk, I get nice models
that I can visualise in paraview - no problems.
Also, using r3.cross.rast for sampling slices seems to give correct
results.
Also, the output from r3.info is perfectly sane, so it seems the
algorithm is working.

My 3D region settings are:

projection: 0 (x,y)
zone:       0
north:      60.85
south:      58.85
west:       14.9
east:       16.95
top:        12.00000000
bottom:     9.00000000
nsres:      0.05
nsres3:     0.05
ewres:      0.05
ewres3:     0.05
tbres:      0.05
rows:       40
rows3:      40
cols:       41
cols3:      41
depths:     60

The same extents apply to my layer model.

HOWEVER, I fail to visualise such a layer model in NVIZ.
The first problem is that NVIZ seems to be unable to correctly
position volumes that have a bottom > 0.0.

You can test this even with a very primitive volume
(use the 3d region settings as shown above):

r3.mapcalc test=1.0

If you load this into NVIZ (and set the bottom of the 3d region > 0.0), 
it will not display the extents of the volume as it should.
If you start rotating the view around a bit and increase
height enough, the extents will be temporarily visible and you will see 
that the whole model has been
displaced from its original z-range 9.0 to 12.0 and is displayed
at 0.0 to 3.0!

Obviously, it is not possible to display slices or isosurfaces in such a
model, as the NVIZ extent and the real data extent do not match. 
So this seems to be a bug in NVIZ(?).

But here is what really puzzles me: if I adjust 3d region settings and
set the bottom to 0.0, NVIZ shows the correct extent, but has the same
problems!
I.e. I only get a display of the extent as wireframe for as long as
rotate the view, change the height or whatever.

Everything disappears after that and I cannot display slices or isosurfaces.

is it possible that NVIZ has trouble with small 3d regions or high 3d
resolution?
Am I doing something totally silly?

Finally, some points that I could not figure out from any documentation
and am still a bit unsure about. Maybe someone could reaffirm my
conclusions about these:

1. GRASS Z coordinates increase towards the top of the region
2. r3.in.ascii expects the bottom slice first in the ASCII input file

I would so much like to get this working. It would be tremendously
useful for archaeological applications of GRASS GIS!

Best,

Benjamin