I’m generating things that are topologically tubes (they aren’t cylinders). I was trying to accomplish this as below - but I get degenerate cells (cell labeled as vtkQuad, but has no area - i.e., it is a line, I think cell 79 is a line):
import numpy as np
from vtkmodules.vtkIOXML import vtkXMLPolyDataWriter
from vtkmodules.vtkCommonCore import vtkPoints
from vtkmodules.vtkCommonDataModel import vtkStructuredGrid
from vtkmodules.vtkFiltersCore import vtkCleanPolyData
from vtkmodules.vtkFiltersGeometry import vtkStructuredGridGeometryFilter
def surface_2(u, v):
x = np.cos(u)
y = np.sin(u)
r = (x**4 + y**4)**(1/4)
x0 = x/2
y0 = y/2
x1 = r * x
y1 = r * y
x = v * x0 + (1 - v) * x1
y = v * y0 + (1 - v) * y1
z = v
return x, y, z
def generate_tube(axial_mesh, radial_mesh):
location_samples = vtkPoints()
for v in axial_mesh:
for u in np.append(radial_mesh, radial_mesh[0]): # periodic in u, append beginning to end
location_samples.InsertNextPoint(*surface_2(u, v))
# location_samples.InsertNextPoint(*surface(u, v))
sg = vtkStructuredGrid()
sg.SetDimensions(radial_mesh.size + 1, axial_mesh.size, 1) # +1 because of periodicity
sg.SetPoints(location_samples)
sgg = vtkStructuredGridGeometryFilter()
sgg.SetInputData(sg)
sgg.SetExtent(0, radial_mesh.size + 1, 0, axial_mesh.size, 0, 0)
cpd = vtkCleanPolyData()
cpd.SetInputConnection(sgg.GetOutputPort())
cpd.Update()
return cpd.GetOutput()
if __name__ == "__main__":
axial_mesh = np.linspace(0, 1, 10)
radial_mesh = np.linspace(0, 1, 40) * 2 * np.pi
t = generate_tube(axial_mesh, radial_mesh)
# ll = t.GetLines()
# print(ll)
# import pdb; pdb.set_trace()
w = vtkXMLPolyDataWriter()
w.SetFileName("tube." + w.GetDefaultFileExtension())
w.SetInputData(t)
w.Write()
That was my original reaction when my coworker pointed out that if you split the screen, open the spreadsheet view, select cell data, and then start cycling around until you find this:
I don’t think think you can fix it be setting tolerances. You could use vtkCellSizeFilter and pick out cells with a very small area e.g. Cell 39 has an area of around 3e-17 then you would have to make sure cells adjacent (38 & 0) have common points, could be very fiddly.
Does it really matter? It’s still a cell, not a line.
Here is another approach which may work for you. Cells 0 and 38 adjoin with no degenerate cell.
It will not be efficient for large datasets.
# See: [Can I avoid or filter out degenerate cells when generating polydata?](https://discourse.vtk.org/t/can-i-avoid-or-filter-out-degenerate-cells-when-generating-polydata/14411)
from decimal import Decimal, ROUND_HALF_UP
import numpy as np
from vtkmodules.vtkCommonCore import vtkPoints
from vtkmodules.vtkCommonDataModel import vtkStructuredGrid
from vtkmodules.vtkFiltersCore import vtkCleanPolyData
from vtkmodules.vtkFiltersGeometry import vtkStructuredGridGeometryFilter
from vtkmodules.vtkFiltersVerdict import vtkCellSizeFilter
from vtkmodules.vtkIOXML import vtkXMLPolyDataWriter
# In order to get rid of cells that look like lines.
# This is Decimal object with an explicit exponent attribute/property.
# It will be interpreted by quantize.
six_places = Decimal('1e-6')
def surface_2(u, v):
x = np.cos(u)
y = np.sin(u)
r = (x ** 4 + y ** 4) ** (1 / 4)
x0 = x / 2
y0 = y / 2
x1 = r * x
y1 = r * y
x = v * x0 + (1 - v) * x1
y = v * y0 + (1 - v) * y1
z = v
# Convert the values to a decimal with a given precision 1.0e-6 in this case.
xx = Decimal(x).quantize(six_places, rounding=ROUND_HALF_UP)
yy = Decimal(y).quantize(six_places, rounding=ROUND_HALF_UP)
zz = Decimal(z).quantize(six_places, rounding=ROUND_HALF_UP)
# These will be automatically converted to float in InsertNextPoint().
return xx, yy, zz
def generate_tube(axial_mesh, radial_mesh):
location_samples = vtkPoints()
for v in axial_mesh:
for u in np.append(radial_mesh, radial_mesh[0]): # periodic in u, append beginning to end
location_samples.InsertNextPoint(tuple(surface_2(u, v)))
sg = vtkStructuredGrid()
sg.SetDimensions(radial_mesh.size + 1, axial_mesh.size, 1) # +1 because of periodicity
sg.SetPoints(location_samples)
sgg = vtkStructuredGridGeometryFilter()
sgg.SetInputData(sg)
sgg.SetExtent(0, radial_mesh.size + 1, 0, axial_mesh.size, 0, 0)
cpd = vtkCleanPolyData()
cpd.SetInputConnection(sgg.GetOutputPort())
# This step is important.
cpd.ConvertPolysToLinesOff()
cpd.Update()
# See what we have, for cells, the area etc.
# This data can be looked at in a ParaView spreadsheet.
cpd1 = vtkCellSizeFilter()
cpd1.SetInputConnection(cpd.GetOutputPort())
cpd1.ComputeAreaOn()
cpd1.ComputeVertexCountOn()
cpd1.Update()
return cpd1.GetOutput()
if __name__ == '__main__':
axial_mesh = np.linspace(0, 1, 10)
radial_mesh = np.linspace(0, 1, 40) * 2 * np.pi
t = generate_tube(axial_mesh, radial_mesh)
w = vtkXMLPolyDataWriter()
w.SetFileName('tube.' + w.GetDefaultFileExtension())
w.SetInputData(t)
w.Write()
