python中的vtk包的部分代码

Please refresh the page if equations are not rendered correctly.
---------------------------------------------------------------

Pyvista应该是原生的vtk包的一个很好替代，以下代码目前基本不再使用，作为记录以防以后需要。

然而，单纯的通过PyVista来调用vtk，却可能使某些功能的运行速度下降百倍以上，这是由于PyVista为了兼顾易用性，将不少函数整合在了一个函数中，在确定调用哪一个时需要不断检查输入的参数进行判断，并相应的去检查数据格式，消耗不少资源。另外，wrapper中还可能包含直接调用时可以避免的循环结构。

def vtu_loader(path, volume_file):
    """
    Load vtu (Unstructured mesh) file.

    Parameters
    ----------
    path: str
        path to vtu file
    volume_file: str
        vtu file name

    Returns
    -------
    volume: vtkUnstructuredGrid
        vtkUnstructuredGrid object
    """
    volume = vtk.vtkXMLUnstructuredGridReader()
    volume.SetFileName(os.path.join(path, volume_file))
    volume.Update()
    volume = volume.GetOutput()
    volume.BuildLinks()
    return volume

import numpy as np
import vtk
import pyvista as pv
from collections import OrderedDict


def get_boundary_faces(volume):
    """ Extract boundary faces from a vtkUnstructuredGrid object.
    :param volume: vtkUnstructuredGrid object
    :return: boundary_faces: a list of boundary faces
    :return surf: a vtkPolyData object containing all boundary faces
    """
    surface_filter = vtk.vtkDataSetSurfaceFilter()
    surface_filter.SetInputData(volume)
    surface_filter.SetPassThroughPointIds(True)
    surface_filter.SetPassThroughCellIds(True)
    surface_filter.Update()

    # convert the vtkPolyData to a pyvista object
    surf = pv.wrap(surface_filter.GetOutputDataObject(0))
    vtkOriginalPointIds = surf['vtkOriginalPointIds']

    for i in range(surf.n_cells):
        cell_points_id = surf.GetCell(i).GetPointIds()
        cell_points_id = pv.vtk_id_list_to_array(cell_points_id)
        try:
            boundary_face = np.vstack([boundary_face, vtkOriginalPointIds[cell_points_id]])
        except NameError:
            boundary_face = vtkOriginalPointIds[cell_points_id]
    return boundary_face, surf

def get_midpoint(face, scale=1.0):
    """
    Get the midpoint of a face.

    Parameters
    ----------
    face : vtk.vtkQuad
        A vtkQuad object.
    scale : float, optional
        The scale factor to convert the unit of points. The default is 1.0.

    Returns
    -------
    midpoint : numpy.ndarray
        The midpoint of the face.
    pts : numpy.ndarray
        The coordinates of the face vertices.
    """
    np = face.GetNumberOfPoints()
    vtkpts = face.GetPoints()
    pts = numpy.array([vtkpts.GetPoint(i) for i in range(np)]) / scale
    midpoint = mean(pts, axis=0)
    return midpoint, pts

def get_angle_position(coordinates):
    """
    Get the angle position of a point in local coordinate system
    :param coordinate: numpy.ndarray, shape(n,2)
    :return angle_position: float in degree
    """
    x = coordinates[:, 0]
    y = coordinates[:, 1]
    angle_position = arctan2(y, x) * 180 / pi
    angle_position[angle_position < 0] += 360
    return angle_position

def find_nearest(array, value):
    """
    # find the closest value in an array
    :param array: array to search (n,)
    :param value: value, scalar
    :return: index and value of the closest value
    """
    idx = (np.abs(array - value)).argmin()
    return idx, array[idx]