Поиск ячеек файла .stl с отрицательной средней кривизной с использованием VTK в python

Question

У меня есть файл .stl, и я пытаюсь найти координаты ячеек с отрицательной средней кривизной, используя VTK и python.Я написал эти коды, которые отлично работают, чтобы изменить цвета ячеек на основе их средней кривизны, но я хочу достичь, это координаты точных ячеек и треугольников с определенной средней кривизной, например, трехмерные координаты ячеек с наибольшей отрицательной средней кривизной.,Вот коды:

import vtk


def gaussian_curve(fileNameSTL):
    colors = vtk.vtkNamedColors()

    reader = vtk.vtkSTLReader()
    reader.SetFileName(fileNameSTL)
    reader.Update()

    curveGauss = vtk.vtkCurvatures()
    curveGauss.SetInputConnection(reader.GetOutputPort())
    curveGauss.SetCurvatureTypeToGaussian() # SetCurvatureTypeToMean() works better in the case of kidney.

    ctf = vtk.vtkColorTransferFunction()
    ctf.SetColorSpaceToDiverging()
    p1 = [0.0] + list(colors.GetColor3d("MidnightBlue"))
    p2 = [1.0] + list(colors.GetColor3d("DarkRed"))
    ctf.AddRGBPoint(*p1)
    ctf.AddRGBPoint(*p2)
    cc = list()
    for i in range(256):
        cc.append(ctf.GetColor(float(i) / 255.0))

    lut = vtk.vtkLookupTable()
    lut.SetNumberOfColors(256)
    for i, item in enumerate(cc):
        lut.SetTableValue(i, item[0], item[1], item[2], 1.0)
    lut.SetRange(0, 0) # In the case of kidney, the (0, 0) worked better.
    lut.Build()

    cmapper = vtk.vtkPolyDataMapper()
    cmapper.SetInputConnection(curveGauss.GetOutputPort())
    cmapper.SetLookupTable(lut)
    cmapper.SetUseLookupTableScalarRange(1)

    cActor = vtk.vtkActor()
    cActor.SetMapper(cmapper)

    return cActor


def render_scene(my_actor_list):
    renderer = vtk.vtkRenderer()
    for arg in my_actor_list:
        renderer.AddActor(arg)
    namedColors = vtk.vtkNamedColors()
    renderer.SetBackground(namedColors.GetColor3d("SlateGray"))

    window = vtk.vtkRenderWindow()
    window.SetWindowName("Render Window")
    window.AddRenderer(renderer)

    interactor = vtk.vtkRenderWindowInteractor()
    interactor.SetRenderWindow(window)

    # Visualize
    window.Render()
    interactor.Start()


if __name__ == '__main__':
    fileName = "400_tri.stl"
    my_list = list()
    my_list.append(gaussian_curve(fileName))
    render_scene(my_list)

Этот код производит красные ячейки для положительной средней кривизны и синие для отрицательных.

Мне нужен результат (координатыячеек) в виде массивов или чего-то в этом роде.Буду признателен за любые предложения и помощь по этой проблеме.

Answer 1

Возможное решение с vtkplotter :

from vtkplotter import *

torus1 = Torus().addCurvatureScalars().addScalarBar()
print("list of scalars:", torus1.scalars())

torus2 = torus1.clone().addScalarBar()
torus2.threshold("Gauss_Curvature", vmin=-15, vmax=0)

show(torus1, torus2, N=2) # plot on 2 separate renderers

print("vertex coordinates:", len(torus2.coordinates()))
print("cell centers      :", len(torus2.cellCenters()))

посмотрите полученный скриншот здесь

Дополнительный пример здесь .

Надеюсь, это поможет.

Answer 2

Итак, я нашел ответ из kitware weblog , вот код, который отлично работает, используя vtk.numpy_interface и vtk.util.numpy_support, но все равно он не выдает normals_array, и я не знаю почему ??

import vtk
from vtk.numpy_interface import dataset_adapter as dsa
from vtk.util.numpy_support import vtk_to_numpy


def curvature_to_numpy(fileNameSTL, curve_type='Mean'):
    colors = vtk.vtkNamedColors()

    reader = vtk.vtkSTLReader()
    reader.SetFileName(fileNameSTL)
    reader.Update()
    # Defining the curvature type.
    curve = vtk.vtkCurvatures()
    curve.SetInputConnection(reader.GetOutputPort())
    if curve_type == "Mean":
        curve.SetCurvatureTypeToMean()
    else:
        curve.SetCurvatureTypeToGaussian()
    curve.Update()

    # Applying color lookup table.
    ctf = vtk.vtkColorTransferFunction()
    ctf.SetColorSpaceToDiverging()
    p1 = [0.0] + list(colors.GetColor3d("MidnightBlue"))
    p2 = [1.0] + list(colors.GetColor3d("DarkOrange"))
    ctf.AddRGBPoint(*p1)
    ctf.AddRGBPoint(*p2)
    cc = list()
    for i in range(256):
        cc.append(ctf.GetColor(float(i) / 255.0))

    lut = vtk.vtkLookupTable()
    lut.SetNumberOfColors(256)
    for i, item in enumerate(cc):
        lut.SetTableValue(i, item[0], item[1], item[2], 1.0)
    lut.SetRange(0, 0)  # In the case of kidney, the (0, 0) worked better.
    lut.Build()

    # Creating Mappers and Actors.
    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInputConnection(curve.GetOutputPort())
    mapper.SetLookupTable(lut)
    mapper.SetUseLookupTableScalarRange(1)

    actor = vtk.vtkActor()
    actor.SetMapper(mapper)

    # Scalar values to numpy array. (Curvature).
    dataObject = dsa.WrapDataObject(curve.GetOutput())
    normals_array = dataObject.PointData['Normals'] # Output array.
    curvature_array = dataObject.PointData['Mean_Curvature'] # output array.

    # Node values to numpy array.
    nodes = curve.GetOutput().GetPoints().GetData()
    nodes_array = vtk_to_numpy(nodes)

    # Creating a report file (.vtk file).
    writer = vtk.vtkPolyDataWriter()
    writer.SetFileName('vtk_file_generic.vtk')
    writer.SetInputConnection(curve.GetOutputPort())
    writer.Write()

   #  EDIT:

   # Creating the point normal array using vtkPolyDataNormals().
    normals = vtk.vtkPolyDataNormals()
    normals.SetInputConnection(reader.GetOutputPort())  # Here "curve" could be replaced by "reader".
    normals.ComputePointNormalsOn()
    normals.SplittingOff()
    normals.Update()
    dataNormals = dsa.WrapDataObject(normals.GetOutput())
    normals_array = dataNormals.PointData["Normals"]



    return actor, normals_array, curvature_array, nodes_array


def render_scene(my_actor_list):
    renderer = vtk.vtkRenderer()
    for arg in my_actor_list:
        renderer.AddActor(arg)
    namedColors = vtk.vtkNamedColors()
    renderer.SetBackground(namedColors.GetColor3d("SlateGray"))

    window = vtk.vtkRenderWindow()
    window.SetWindowName("Render Window")
    window.AddRenderer(renderer)

    interactor = vtk.vtkRenderWindowInteractor()
    interactor.SetRenderWindow(window)

    # Visualize
    window.Render()
    interactor.Start()


if __name__ == '__main__':
    filename = "400_tri.stl"
    my_list = list()
    my_actor, my_normals, my_curve, my_nodes = curvature_to_numpy(filename, curve_type="Mean")
    my_list.append(my_actor)
    render_scene(my_list) # Visualization.
    print(my_nodes) # Data points.
    print(my_normals) # Normal vectors.
    print(my_curve) # Mean curvatures.