VTK 实现对vtkpolyata 膨胀放缩

吃肉夹馍不要夹馍

已于 2024-06-21 10:35:08 修改

阅读量392

点赞数 6

分类专栏： VTK 文章标签： vr python 笔记图像处理图形渲染

于 2024-06-13 14:42:24 首次发布

本文链接：https://blog.csdn.net/sinat_41053216/article/details/139653838

版权

VTK 专栏收录该内容

5 篇文章 0 订阅

订阅专栏

VTK 实现对vtkpolyata 膨胀放缩

实现VTKpolydata 向外扩大，目前了解到的有两种方法

方法1

使用vtkWarpVector 通过设置放缩因子和normals，使得vtkpolyata向外膨胀的效果

但是该种方式更多的是成比例的膨胀，即膨胀多少倍等，很难做到具体的控制膨胀多少距离

另外，注意放缩后可能没那么平滑了，或许需要进行stripper去除表面多边形
在这里插入图片描述

def update_edge_size(input_polydata, scaleFactor, iteration=50):
    """vtkWarpVector 放缩 vtkpolyata"""
    
    
    # 先将polydata表面进行一定平滑处理
    smooth_filter = vtk.vtkSmoothPolyDataFilter()
    smooth_filter.SetInputData(input_polydata)
    smooth_filter.SetNumberOfIterations(iteration)
    smooth_filter.SetRelaxationFactor(0.2)
    smooth_filter.Update()
    
	# 清理polydata 数据中的冗余信息，包括重复的点等，简化数据集
    clean = vtk.vtkCleanPolyData()
    clean.SetInputData(smooth_filter.GetOutput())
    
    normals = vtk.vtkPolyDataNormals()
    normals.SetInputConnection(clean.GetOutputPort())
    normals.SplittingOff()
    # 使用vtkWarpVector 进行放缩
    warp = vtk.vtkWarpVector()
    warp.SetInputConnection(normals.GetOutputPort())
    warp.SetInputArrayToProcess(
        0, 0, 0,
        vtk.vtkDataObject().FIELD_ASSOCIATION_POINTS,
        vtk.vtkDataSetAttributes().NORMALS
    )
    warp.SetScaleFactor(scaleFactor)
    warp.Update()
    warped_polydata = warp.GetOutput()
    return warped_polydata

方法2

使用 vtkImplicitPolyDataDistance 计算一个有向距离场（signed distance field），距离场的数据表示距离polydata表面的距离，在polydata外部值为正数，内部值为负数，表面值为0

具体官方例子可见这里

官方例子中还是以点的形式展示了有向距离场

但是我的目的是通过SDF 实现对polydata按照指定长度的膨胀，就是说膨胀3mm的时候，希望膨胀后的polydata表面距离原先polydata表面的平均距离也是3mm

考虑了一下大致思路是，通过vtkImplicitPolyDataDistance 得到一个有向距离场的矩阵，然后根据希望放大的长度删选除膨胀后的点数据，得到一个mask，再将mask转换为polydata，得到膨胀了指定距离的polydata

在这里插入图片描述

def sdf_polydata(polydata:vtk.vtkPolyData, volume_sitk_img: sitk.Image,target_size=20):
    """
		通过SDF 膨胀 polydata
		params:
			polydata: 输入待膨胀的polydata
			volume_sitk_img： sitk.Image 相应的体数据，主要是为了得到矩阵大小以及进行RAS坐标和IJK世界坐标的转换，如果有指定的矩阵大小和RAS坐标和IJK坐标转换则不需要
			target_size： 希望膨胀的大小
        """
    if not target_size or target_size <= 1:
        return None

    volume_arr = sitk.GetArrayFromImage(volume_sitk_img)
    sdf_arr = np.zeros_like(volume_arr, np.float16)
    origin = volume_sitk_img.GetOrigin()
    spacing = volume_sitk_img.GetSpacing()

    volume_shape = volume_arr.shape

    # 膨胀后 mask的包围盒
    polydata_bounds = [0] * 6
    polydata.GetBounds(polydata_bounds)
    ras_x_min, ras_x_max, ras_y_min, ras_y_max, ras_z_min, ras_z_max = polydata_bounds

    ras_x_min, ras_y_min, ras_z_min = (
        ras_x_min - target_size - 1,
        ras_y_min - target_size - 1,
        ras_z_min - target_size - 1
    )
    ras_x_max, ras_y_max, ras_z_max = (
        ras_x_max + target_size + 1,
        ras_y_max + target_size + 1,
        ras_z_max + target_size + 1
    )

    impl_distance = vtk.vtkImplicitPolyDataDistance()
    impl_distance.SetInput(polydata)

    min_ijk_point = volume_sitk_img.TransformPhysicalPointToIndex([ras_x_min, ras_y_min, ras_z_min])
    max_ijk_point = volume_sitk_img.TransformPhysicalPointToIndex([ras_x_max, ras_y_max, ras_z_max])

    # 这部分比较耗时
    # 遍历包围盒，计算包围盒中每个值保存的是当前点距离polydata表面的距离
    for x in range(min_ijk_point[0], max_ijk_point[0]):
        for y in range(min_ijk_point[1], max_ijk_point[1]):
            for z in range(min_ijk_point[2], max_ijk_point[2]):
                tmp_ijk_point = [x, y, z]
                tmp_ras_point = volume_sitk_img.TransformIndexToPhysicalPoint(tmp_ijk_point)

                if not 0 <= tmp_ijk_point[0] <= volume_shape[2] or not 0 <= tmp_ijk_point[1] <= volume_shape[
                    1] or not 0 <= tmp_ijk_point[2] <= volume_shape[0]:
                    print(f"point {tmp_ijk_point} is over rand {volume_shape}")
                    continue
				# 计算当前点到polydata表面的距离， 由于计算需要ras坐标，所以把索引下标再转回了ras
                distance = impl_distance.EvaluateFunction(tmp_ras_point)
                sdf_arr[tmp_ijk_point[2], tmp_ijk_point[1], tmp_ijk_point[0]] = round(distance, 2)


    # 得到膨胀后mask
    new_mask = np.zeros_like(sdf_arr, np.uint8)
    logic = np.logical_and(sdf_arr <= target_size, sdf_arr != 0)
    new_mask[logic] = 1
	
    # 进行开运算，去掉膨胀后mask中间的一些空洞
    for i in range(0, new_mask.shape[0]):
        new_mask[i] = open_compute(new_mask[i], 3)

    # 得到膨胀后的polydatamask之后需要将其转为polydata
    # 分两步走。 1，mask转为imagedata; 2, imagedata转为polydata

    # 将膨胀后mask转为imagedata
    tmp_imagedata = vtk.vtkImageData()
    tmp_imagedata.SetDimensions(*volume_sitk_img.GetSize())
    tmp_imagedata.SetOrigin(origin)
    tmp_imagedata.SetSpacing(spacing)

    vtk_array = numpy_support.numpy_to_vtk(new_mask.flatten(), array_type=vtk.VTK_INT)
    tmp_imagedata.GetPointData().SetScalars(vtk_array)
	
    # imagedata转为polydata
    marching_cubes = vtk.vtkDiscreteMarchingCubes()
    marching_cubes.SetInputData(tmp_imagedata)
    marching_cubes.SetValue(0, 1)
    # 设置法向量计算开启
    marching_cubes.ComputeNormalsOn()
    marching_cubes.ComputeScalarsOn()
    marching_cubes.Update()
    res_polydata = marching_cubes.GetOutput()
    return res_polydata


def open_compute(img_arr, kernel_size=3, iterations=1):
	# 开运算
    kernel = np.ones((kernel_size, kernel_size), np.int8)
    dilated_arr = cv2.dilate(img_arr, kernel, iterations=iterations)
    erode_arr = cv2.erode(dilated_arr, kernel, iterations=iterations)
    return erode_arr

补充：

使用vtkImplicitPolyDataDistance 计算时候可能会受到 polydata的normals的影响，当法线指向内部时候，外部的距离变为负数，内部变为正数。所以需要法线指向外部（更方便计算处理）

修正polydata 法线：

def updatepolydata_normals(polydata):
	normals = vtk.vtkPolyDataNormals()
    normal_filter.SetInputData(polydata)
    normal_filter.SetComputePointNormals(0)
    normal_filter.SetComputeCellNormals(1)
    normal_filter.SetAutoOrientNormals(1)
    normal_filter.Update()
    polydata = normal_filter.GetOutput()
    return polydata