Exporting python polar plot into VTK

ParaView Support
1

I have python numpy arrays of radius, theta, and values that will generate the following image
(radius: 1d, theta: 1d, values: 2d [dim_theta, dim_radius])

fig, ax = plt.subplots(subplot_kw=dict(projection='polar'))
meshRadius, meshTheta = np.meshgrid(radius, theta)
ax.contourf(meshTheta, meshRadius, values)

image

I tried to export this into VTK with the following

    import vtk
    from vtk.util import numpy_support

    npts = np.shape(values)
    Points = vtk.vtkPoints()
    Vertices = vtk.vtkCellArray()
    values_vtk = []
    for i_radius in range(npts[1]):
        _r = radius[i_radius]
        for i_theta in range(npts[0]):
            _theta = theta[i_theta]
            _id = Points.InsertNextPoint(_r*np.cos(_theta), _r*np.sin(_theta), 0.)
            Vertices.InsertNextCell(1)
            Vertices.InsertCellPoint(_id)
            values_vtk.append(values[i_theta,i_radius])

    values_vtk = np.array(values_vtk)

    polyData = vtk.vtkPolyData()
    polyData.SetPoints(Points)
    polyData.SetVerts(Vertices)

    pointData = polyData.GetPointData()
    values_vtk2 = numpy_support.numpy_to_vtk(values_vtk)
    values_vtk2.SetName("InvDFT")
    pointData.AddArray(values_vtk2)

    writer = vtk.vtkPolyDataWriter()
    writer.SetFileName("vtk.vtk")
    writer.SetInputData(polyData)
    writer.Write()

image
image1520×874 115 KB

I think I am almost there, but somehow I should specify the connectivities sort of things to make this ‘surface’ plot.

Could anyone help a bit here?

2

This is more of a VTK question, but we can follow up here.

That’s correct. You have inserted 0-dimensional (no connectivity) cells, one for each point. Instead, you will need to insert a 2-dimensional cell, such as quads or triangles, to produce a mesh like you produce with your call to np.meshgrid.

There isn’t a convenient call like that in VTK. Instead, I recommend figuring our the connectivity among neighboring points and setting up cells with those point IDs manually. You’ll probably want to insert the points first in a loop and the insert the cells in a later loop. You don’t need to know the point positions to define the cells, just the point IDs and how they are arranged in space to determine the polygon you want to insert (an exercise left to the poster).

1 Like
3

Thanks to the comments from @cory.quammen , I get it to work.

image
image1523×873 157 KB

The point hindered me the most was that my wrong usage of ‘SetVerts’ instead of ‘SetPolys’, which is commented in the following modified version of the Python script.

    import vtk
    from vtk.util import numpy_support

    npts = np.shape(values)
    Points = vtk.vtkPoints()
    Vertices = vtk.vtkCellArray()
    values_vtk = []
    # Modification 1) Now I am storing ids array to make connectivity construction easier
    ids = np.zeros((npts[1],npts[0]))
    for i_radius in range(npts[1]):
        _r = radius[i_radius]
        for i_theta in range(npts[0]):
            _theta = theta[i_theta]
            _id = Points.InsertNextPoint(_r*np.cos(_theta), _r*np.sin(_theta), 0.)
            ids[i_theta,i_radius] = _id
            values_vtk.append(values[i_theta,i_radius])

    values_vtk = np.array(values_vtk)

    # Modification 2) Separating the connectivity loop
    for i_radius in range(npts[1]-1):  # Note: be aware of '-1' following npts[0]
        for i_theta in range(npts[0]-1):  # Note: be aware of '-1' following npts[0]
            # Modification 3) Now I am using '4' instead of '1' for InsertNextCell
            Vertices.InsertNextCell(4)
            Vertices.InsertCellPoint(int(ids[i_theta,i_radius]))
            Vertices.InsertCellPoint(int(ids[i_theta,i_radius+1]))
            Vertices.InsertCellPoint(int(ids[i_theta+1,i_radius+1]))
            Vertices.InsertCellPoint(int(ids[i_theta+1,i_radius]))

    polyData = vtk.vtkPolyData()
    polyData.SetPoints(Points)
    # Modiciation 4) SetVerts have been replaced by SetPolys
    #polyData.SetVerts(Vertices)
    polyData.SetPolys(Vertices)

    pointData = polyData.GetPointData()
    values_vtk2 = numpy_support.numpy_to_vtk(values_vtk)
    values_vtk2.SetName("InvDFT")
    pointData.AddArray(values_vtk2)

    writer = vtk.vtkPolyDataWriter()
    writer.SetFileName("vtk.vtk")
    writer.SetInputData(polyData)
    writer.Write()
4

Awesome! Nice job!