Python在3D数据中根据点的位置以及法线方向提取切片

Python在3D数据中根据点以及法线方向提取切片

在处理3D图像数据（CT图像，MRI图像）中，有时需要对其进行提取切片的操作，一般都是按照X,Y,Z三个轴来提取切片，然而在一些情况下需要根据某点的法线方向提取过这点的切片，在这里根据的 matlab代码 改写成python的版本。

原理

在matlab的代码上主要是分为几步来获得斜切平面的
（1）初始化切平面，在这里切平面可以当成一把刀，初始化的切平面是平行于XY平面的。
（2）通过绕一点以及一个方向旋转改切平面，这样的话就知道往哪切了，只要将初始切平面的法线方向变换至切平面法线的方向，同时对初始平面上的点进行同样的变换，便可得到最后的切平面。
（3）因为旋转之后的点会产生小数，所以还需要取整操作。

代码

import numpy as np
import math
import scipy.linalg as linalg
import nibabel as nib
import matplotlib.pyplot as plt

def sphere(shape, radius, position):
    semisizes = (radius,) * 3
    grid = [slice(-x0, dim - x0) for x0, dim in zip(position, shape)]
    position = np.ogrid[grid]
    arr = np.zeros(shape, dtype=float)
    for x_i, semisize in zip(position, semisizes):
     arr += (np.abs(x_i/semisize) ** 2)
    return arr <= 1.0


def loc_convert(loc, axis, radian):
    '''
	实现点按某个轴旋转一定度数后得到新点坐标
    :param loc:原始坐标
    :param axis: 绕旋转的轴（点）
    :param radian: 角度
    :return: 新坐标
    '''
    radian = np.deg2rad(radian)
    rot_matrix = linalg.expm(np.cross(np.eye(3), axis / linalg.norm(axis) * radian))
    new_loc = np.dot(rot_matrix, loc)
    return new_loc

def extract_slice(img, c, v, radius):
    '''
    :param V:3d 图像
    :param center-c: 中心（x,y,z）
    :param normal-v: 法向量（v1,v2,v3）
    :param radius: 半径，即边长的一半
    :return:
    slicer：得到的2d切片
    loc: 得到切片对应的原3d坐标
    '''
    # 设置初始面
    epsilon = 1e-12
    x = np.arange(-radius, radius, 1)
    y = np.arange(-radius, radius, 1)
    X, Y = np.meshgrid(x, y)
    Z = np.zeros_like(X)
    loc = np.array([X.flatten(), Y.flatten(), Z.flatten()])

    # 设置初始平面，垂直于Z轴，将向量变为单位向量
    hspInitialVector = np.array([0, 0, 1])
    h_norm = np.linalg.norm(hspInitialVector)
    h_v = hspInitialVector / h_norm
    h_v[h_v == 0] = epsilon
    v = v / np.linalg.norm(v)
    v[v == 0] = epsilon

    # 计算初始法线与最后法线的角度
    hspVecXvec = np.cross(h_v, v) / np.linalg.norm(np.cross(h_v, v))
    acosineVal = np.arccos(np.dot(h_v, v))
    hspVecXvec[np.isnan(hspVecXvec)] = epsilon
    acosineVal = epsilon if np.isnan(acosineVal) else acosineVal

    # 得到旋转后的坐标
    loc = loc_convert(loc, hspVecXvec, 180 * acosineVal / math.pi)
    sub_loc = loc + np.reshape(c, (3, 1))
    loc = np.round(sub_loc)
    loc = np.reshape(loc, (3, X.shape[0], X.shape[1]))

    # 生成初始切片，以及对应的索引值
    sliceInd = np.zeros_like(X, dtype=np.float)
    sliceInd[sliceInd == 0] = np.nan
    slicer = np.copy(sliceInd)

    # 将3D图像对应的像素值以及对应的坐标赋值给对应的切片
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            if loc[0, i, j] >= 0 and loc[0, i, j] < img.shape[0] and loc[1, i, j] >= 0 and loc[1, i, j] < img.shape[1] and loc[2, i, j] >= 0 and loc[2, i, j] < img.shape[2]:
                slicer[i, j] = img[
                    loc[0, i, j].astype(np.int), loc[1, i, j].astype(np.int), loc[2, i, j].astype(np.int)]
    slicer[np.isnan(slicer)]=0
    return slicer, sub_loc,loc

if __name__ == '__main__':

    arr = sphere((256, 256, 256), 100, (127, 127, 127))
    print(arr.shape)
    c = [150,150,150]
    n = [1, 2, 3]
    slicer, sub_loc, loc = extract_slice(arr, c, n, 100)
    print(slicer)
    plt.imshow(slicer)
    plt.show()