# Extracting isosurface from a vtkImageData
**URL:** https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502
**Category:** Support
**Tags:** python
**Created:** [January 19, 2023, 2:30pm UTC](https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502 "2023-01-19T14:30:40Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@Jishnu](https://discourse.vtk.org/u/Jishnu)
#### Post date: [January 19, 2023, 2:30pm UTC](https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502/1 "2023-01-19T14:30:40Z")
I am trying to create an isosurface plot from a vtkImageData. I have a numpy array which has the scalar fields. I kept the scalar field as 0 and 1 for testing. Given below is my pipeline:
```auto
import vtk
from vtk.util import numpy_support
import numpy as np
# Load scalar data from npy file
scalars=np.load("Scalar.npy")
# Convert numpy array to vtk
scalars=numpy_support.numpy_to_vtk(scalars)
scalars.SetName("Scalar")
# Create a vtkImageData
grid =vtk.vtkImageData()
grid.SetDimensions(251, 251, 101)
grid.SetSpacing(0.1, 0.1, 0.1)
grid.SetOrigin(0,0,0)
# Add scalar field to the ImageData
grid.GetPointData().AddArray(scalars)
# Get Isosurface using contour filter
contourFilter = vtk.vtkContourFilter()
contourFilter.SetInputData(grid)
contourFilter.SetArrayComponent (0)
contourFilter.SetValue(0, 1)
contourFilter.Update()
outData=contourFilter.GetOutput()
```
However, I am unable to get the isosurface from the contour filter. The number of points in outData outData. `GetNumberOfPoints()` is giving 0. However, I have value 1 at 211573 places. `np.count_nonzero(scalars[scalars==1])` gives 211573. I have saved the numpy array used in the above code [here](https://wetransfer.com/downloads/a5560af1a8fea58d488955a15349cfb520230119142404/1ca6bc).
Please help me to solve this issue and understand what I am doing wrong here.
---
### Author:  [@will.schroeder](https://discourse.vtk.org/u/will.schroeder)
#### Post date: [January 19, 2023, 2:38pm UTC](https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502/2 "2023-01-19T14:38:19Z")
What happens when you set the isovalue ==0.5 ?
---
### Author:  [@Jishnu](https://discourse.vtk.org/u/Jishnu)
#### Post date: [January 19, 2023, 2:46pm UTC](https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502/3 "2023-01-19T14:46:58Z")
Hi @will.schroeder
isovalue ==0.5 ( `contourFilter.SetValue(0, 0.5)`) is also giving `outData.GetNumberOfPoints()` as 0, I mean contour filter is returning empty polydata.
---
### Author:  [@will.schroeder](https://discourse.vtk.org/u/will.schroeder)
#### Post date: [January 19, 2023, 4:21pm UTC](https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502/4 "2023-01-19T16:21:39Z")
This makes me suspicious that the input data to contourFilter is somehow not correct. We’d have to look at the data to know for sure.
---
### Author:  [@Jishnu](https://discourse.vtk.org/u/Jishnu)
#### Post date: [January 20, 2023, 3:24am UTC](https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502/5 "2023-01-20T03:24:43Z")
Yes, but I am not sure what to check. The size of numpy array which I am assigning as scalar is of the same size as no of points in vtkImageData. I am a bit confused here. Is there any issue in the way I create vtkImageData?
---
### Author:  [@Jishnu](https://discourse.vtk.org/u/Jishnu)
#### Post date: [January 20, 2023, 4:27am UTC](https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502/6 "2023-01-20T04:27:46Z")
Hi @will.schroeder,
I found the mistake. To add the scalar data to vtkImageData, I have to use `grid.GetPointData().SetScalars(scalar)` method rather than `grid.GetPointData().AddArray(scalars)` method. Now I am getting points in the vtkContourFilter output.
However, I am not getting what I am supposed to get. In the above example, I have created a rectangular grid of size (x,y,z) 25,25,25 and grid spacing (dx,dy,dz) of 0.1,0.1,0.25 and assigned value 1 for points which are inside a sphere of radius 5, with center (x,y,z) at 12.5,12.5,12.5. For the rest of the points which are outside the sphere, I am assigned value of 0. So from the contour filter, I am expecting to get the same sphere as the output. But what I am getting is this:
Can you please help me to figure out the right method to solve my problem?
---
### Author:  [@will.schroeder](https://discourse.vtk.org/u/will.schroeder)
#### Post date: [January 20, 2023, 9:44am UTC](https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502/7 "2023-01-20T09:44:03Z")
I’m glad you are making progress. The data is most certainly bad. There are many ways to debug this, one of the simplest ways is to write the data out, and then use ParaView or SLicer to look at it. You can also write out to a legacy VTK files in text-readable form so you can examine it.
---
### Author:  [@marcomusy](https://discourse.vtk.org/u/marcomusy)
#### Post date: [January 20, 2023, 1:48pm UTC](https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502/8 "2023-01-20T13:48:50Z")
i think you just swapped 2 axes 🙂
```python
from vedo import *
scals = np.load("Scalar.npy").reshape(251,251,101)
vol = Volume(scals, spacing=(0.1,0.1,0.25))
vol.isosurface(0.5).show(axes=1)
```
---
### Author:  [@Jishnu](https://discourse.vtk.org/u/Jishnu)
#### Post date: [January 21, 2023, 6:01am UTC](https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502/9 "2023-01-21T06:01:06Z")
Hi @marcomusy, Thank you very much for pointing this out. You are right. I am getting a sphere as contour output as given below:
However, as you see the surface is very rough, is there any way to get a smooth sphere?
@will.schroeder, @marcomusy, based on your experience, am I using the right method here?
---
### Author:  [@will.schroeder](https://discourse.vtk.org/u/will.schroeder)
#### Post date: [January 21, 2023, 10:27am UTC](https://discourse.vtk.org/t/extracting-isosurface-from-a-vtkimagedata/10502/10 "2023-01-21T10:27:28Z")
Smoothing is typically performed with vtkWindowedSincPolyDataFilter. Also, vtkSmoothPolyDataFilter and vtkConstrainedPolyDataFilter can be used.