segfault with MarchingCubes

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

segfault with MarchingCubes

Clément D.
Hello all,

I am trying to use the Marching Cubes algorithm to show a 3D surface of a 3D image given by a Python (numpy array).
I thought I had it all mapped out, things were going great with arrays of small size (generated and rendered 3D isosurfaces), but now that I'm working with arrays of slightly bigger size, the programm crashes with Python error 11 SIGSEGV which means seg fault. What worries me a bit is that the array I'm using is not that big, 150*150*150, and I have a fairly recent computer who shouldn't have any problem handling it.

Code MWE is below (in Python)

=========================================== Code snippet =======================

import numpy as np
import vtk

# Creating a numpy array representing a 3D image
npa=np.zeros((150,150,150)).astype(np.uint8)
npa[100:102,100:102,100:102]=255

# Importing it as  object vtk (based on vtk example http://www.vtk.org/Wiki/VTK/Examples/Python/vtkWithNumpy)
dataImporter = vtk.vtkImageImport()
dataString = npa.tostring()
[hh, ww, dd] = npa.shape
dataImporter.CopyImportVoidPointer(dataString, len(dataString))
dataImporter.SetDataScalarTypeToUnsignedChar()
dataImporter.SetNumberOfScalarComponents(3)
dataImporter.SetDataExtent(0, dd - 1, 0, ww - 1, 0, hh - 1)
dataImporter.SetWholeExtent(0, dd - 1, 0, ww - 1, 0, hh - 1)
dataImporter.Update()
vol = dataImporter.GetOutputPort()
# (print(dataImporter.GetOutput().GetClassName()) yields "vtk.ImageData")

# Use MarchingCubes to generate iso-surface
surface = vtk.vtkMarchingCubes();
surface.SetInputConnection(vol);
surface.ComputeNormalsOn();
surface.ComputeScalarsOn();
isoValue = 128
surface.SetValue(0, isoValue);

#This line causes segfault if npa is "too big"
surface.Update();

# .... rest of pipeline
# decimator = vtk.vtkDecimatePro(); ....
# smoother = vtk.vtkSmoothPolyDataFilter(); ....
# connectivityFilter = vtk.vtkPolyDataConnectivityFilter(); ....
# mesh = vtk.vtkPolyData(); ....
# mapper = vtk.vtkPolyDataMapper(); ....
# actor = vtk.vtkActor(); ....


=========================================== End of code snippet =======================



So I have a couple of questions here :
1/ Is this a problem of input for MarchingCubes being too big ? Or is this segfault indicating something else ?
- arguments for the former : Reducing the size of npa in this code makes it work
- arguments for the latter : This code I'm showing is a MWE for my code which is more complicated, and in my code, this code portion fails with a large array, but succeeds if I reduce in size *another* array which has nothing to do with MarchingCubes() (again, that other array is not that big, so I don't think I'm reaching the limits of Python memory or anything like that...)

2/ Am I doing the whole import numpy thing correctly ? I've based myself on the vtk example, but I am unsure of whether it is the recommended way for my case or if it was case specific


Thanks everyone,
Clément

_______________________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html

Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ

Search the list archives at: http://markmail.org/search/?q=vtkusers

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/vtkusers
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: segfault with MarchingCubes

David Gobbi
Hi Clément,

The "number of scalar components" should be 1 for your data:

  dataImporter.SetNumberOfScalarComponents(1)

There are a couple places where you can improve efficiency:

  # directly create an array of the desired type
  npa=np.zeros((150,150,150), np.uint8)

  # directly pass the array to VTK (without tostring())
  dataImporter.CopyImportVoidPointer(npa, npa.size)

Cheers,
 - David


On Mon, Jul 10, 2017 at 2:29 AM, Clément D. <[hidden email]> wrote:
Hello all,

I am trying to use the Marching Cubes algorithm to show a 3D surface of a 3D image given by a Python (numpy array).
I thought I had it all mapped out, things were going great with arrays of small size (generated and rendered 3D isosurfaces), but now that I'm working with arrays of slightly bigger size, the programm crashes with Python error 11 SIGSEGV which means seg fault. What worries me a bit is that the array I'm using is not that big, 150*150*150, and I have a fairly recent computer who shouldn't have any problem handling it.

Code MWE is below (in Python)

=========================================== Code snippet =======================

import numpy as np
import vtk

# Creating a numpy array representing a 3D image
npa=np.zeros((150,150,150)).astype(np.uint8)
npa[100:102,100:102,100:102]=255

# Importing it as  object vtk (based on vtk example http://www.vtk.org/Wiki/VTK/Examples/Python/vtkWithNumpy)
dataImporter = vtk.vtkImageImport()
dataString = npa.tostring()
[hh, ww, dd] = npa.shape
dataImporter.CopyImportVoidPointer(dataString, len(dataString))
dataImporter.SetDataScalarTypeToUnsignedChar()
dataImporter.SetNumberOfScalarComponents(3)
dataImporter.SetDataExtent(0, dd - 1, 0, ww - 1, 0, hh - 1)
dataImporter.SetWholeExtent(0, dd - 1, 0, ww - 1, 0, hh - 1)
dataImporter.Update()
vol = dataImporter.GetOutputPort()
# (print(dataImporter.GetOutput().GetClassName()) yields "vtk.ImageData")

# Use MarchingCubes to generate iso-surface
surface = vtk.vtkMarchingCubes();
surface.SetInputConnection(vol);
surface.ComputeNormalsOn();
surface.ComputeScalarsOn();
isoValue = 128
surface.SetValue(0, isoValue);

#This line causes segfault if npa is "too big"
surface.Update();

# .... rest of pipeline
# decimator = vtk.vtkDecimatePro(); ....
# smoother = vtk.vtkSmoothPolyDataFilter(); ....
# connectivityFilter = vtk.vtkPolyDataConnectivityFilter(); ....
# mesh = vtk.vtkPolyData(); ....
# mapper = vtk.vtkPolyDataMapper(); ....
# actor = vtk.vtkActor(); ....


=========================================== End of code snippet =======================



So I have a couple of questions here :
1/ Is this a problem of input for MarchingCubes being too big ? Or is this segfault indicating something else ?
- arguments for the former : Reducing the size of npa in this code makes it work
- arguments for the latter : This code I'm showing is a MWE for my code which is more complicated, and in my code, this code portion fails with a large array, but succeeds if I reduce in size *another* array which has nothing to do with MarchingCubes() (again, that other array is not that big, so I don't think I'm reaching the limits of Python memory or anything like that...)

2/ Am I doing the whole import numpy thing correctly ? I've based myself on the vtk example, but I am unsure of whether it is the recommended way for my case or if it was case specific


Thanks everyone,
Clément

_______________________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html

Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ

Search the list archives at: http://markmail.org/search/?q=vtkusers

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/vtkusers
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: segfault with MarchingCubes

Clément D.
Great David, thanks again.
The SetNumberOfScalarComponents(3) was indeed a residual of a not-careful-enough copy-paste from another part of the code.
Everything's working now, and thanks for the optimisation tips !


2017-07-10 14:48 GMT+02:00 David Gobbi <[hidden email]>:
Hi Clément,

The "number of scalar components" should be 1 for your data:

  dataImporter.SetNumberOfScalarComponents(1)

There are a couple places where you can improve efficiency:

  # directly create an array of the desired type
  npa=np.zeros((150,150,150), np.uint8)

  # directly pass the array to VTK (without tostring())
  dataImporter.CopyImportVoidPointer(npa, npa.size)

Cheers,
 - David


On Mon, Jul 10, 2017 at 2:29 AM, Clément D. <[hidden email]> wrote:
Hello all,

I am trying to use the Marching Cubes algorithm to show a 3D surface of a 3D image given by a Python (numpy array).
I thought I had it all mapped out, things were going great with arrays of small size (generated and rendered 3D isosurfaces), but now that I'm working with arrays of slightly bigger size, the programm crashes with Python error 11 SIGSEGV which means seg fault. What worries me a bit is that the array I'm using is not that big, 150*150*150, and I have a fairly recent computer who shouldn't have any problem handling it.

Code MWE is below (in Python)

=========================================== Code snippet =======================

import numpy as np
import vtk

# Creating a numpy array representing a 3D image
npa=np.zeros((150,150,150)).astype(np.uint8)
npa[100:102,100:102,100:102]=255

# Importing it as  object vtk (based on vtk example http://www.vtk.org/Wiki/VTK/Examples/Python/vtkWithNumpy)
dataImporter = vtk.vtkImageImport()
dataString = npa.tostring()
[hh, ww, dd] = npa.shape
dataImporter.CopyImportVoidPointer(dataString, len(dataString))
dataImporter.SetDataScalarTypeToUnsignedChar()
dataImporter.SetNumberOfScalarComponents(3)
dataImporter.SetDataExtent(0, dd - 1, 0, ww - 1, 0, hh - 1)
dataImporter.SetWholeExtent(0, dd - 1, 0, ww - 1, 0, hh - 1)
dataImporter.Update()
vol = dataImporter.GetOutputPort()
# (print(dataImporter.GetOutput().GetClassName()) yields "vtk.ImageData")

# Use MarchingCubes to generate iso-surface
surface = vtk.vtkMarchingCubes();
surface.SetInputConnection(vol);
surface.ComputeNormalsOn();
surface.ComputeScalarsOn();
isoValue = 128
surface.SetValue(0, isoValue);

#This line causes segfault if npa is "too big"
surface.Update();

# .... rest of pipeline
# decimator = vtk.vtkDecimatePro(); ....
# smoother = vtk.vtkSmoothPolyDataFilter(); ....
# connectivityFilter = vtk.vtkPolyDataConnectivityFilter(); ....
# mesh = vtk.vtkPolyData(); ....
# mapper = vtk.vtkPolyDataMapper(); ....
# actor = vtk.vtkActor(); ....


=========================================== End of code snippet =======================



So I have a couple of questions here :
1/ Is this a problem of input for MarchingCubes being too big ? Or is this segfault indicating something else ?
- arguments for the former : Reducing the size of npa in this code makes it work
- arguments for the latter : This code I'm showing is a MWE for my code which is more complicated, and in my code, this code portion fails with a large array, but succeeds if I reduce in size *another* array which has nothing to do with MarchingCubes() (again, that other array is not that big, so I don't think I'm reaching the limits of Python memory or anything like that...)

2/ Am I doing the whole import numpy thing correctly ? I've based myself on the vtk example, but I am unsure of whether it is the recommended way for my case or if it was case specific


Thanks everyone,
Clément


_______________________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html

Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ

Search the list archives at: http://markmail.org/search/?q=vtkusers

Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/vtkusers
Loading...