[vtkusers] Confused about vtkImageData + vtkMarchingCubes [without HTML formatting this time]

[Christian Jauvin]

[Sorry to post this message a second time; from this archive:

http://www.vtk.org/pipermail/vtkusers/2012-January/121070.html

it appears that my message was somehow truncated.. maybe due to the
Gmail HTML formatting I used for the fixed-width font?]

First my final goal: given a point cloud, downsample it, and visualize
the resultings voxels (i.e. 1 if containing any point, 0 otherwise) in
an efficient way, using Python-VTK. I already succeeded in doing it
with my own downsampling code, and a voxel visualization mechanism
crafted from vtkUnstructuredGrid, vtkVoxel and vtkPoints. However I
have the strong feeling that there must be a more efficient way of
doing it, and from what I read I gather that a vtkImageData, rendered
with vtkMarchingCubes might work well.

I tried really hard to build an understanding of those objects by
triangulating the many pieces of code and information I gathered
everywhere, but I just can't make it this time. The toy model I'm
trying to first build is a simple 3x3x3 grid where the 3 diagonal
voxels would be "on":

import vtk
n = 3
grid = vtk.vtkImageData()
grid.SetDimensions(n, n, n)
grid.SetOrigin(0, 0, 0)
grid.SetSpacing(1, 1, 1)

>From there, I tried two things: first, to set the desired voxels explicitly:

grid.SetScalarComponentFromFloat(0,0,0,0,1)
grid.SetScalarComponentFromFloat(1,1,1,0,1)
grid.SetScalarComponentFromFloat(2,2,2,0,1)

The other way I tried is to set the scalar field:

scalars = vtk.vtkUnsignedShortArray()
for i in range(n):
    for j in range(n):
        for k in range(n):
            scalars.InsertNextTuple1(1 if i == j == k else 0)
grid.GetPointData().SetScalars(scalars)

But then when I try to render the resulting grid with marching cubes,
it always yields some strange shape, very far from my goal:

surface = vtk.vtkMarchingCubes()
surface.SetInput(grid)
surface.SetValue(0, 0.5)
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(surface.GetOutputPort())
actor = vtk.vtkActor()
actor.SetMapper(mapper)
ren = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
ren.AddActor(actor)
iren.Initialize()
renWin.Render()
iren.Start()

I'd greatly appreciate any help, thanks in advance.

Christian