bresenham_3d绘制3维直线，可设置直线粗细

  1 import matplotlib.pyplot as plt
  2 from mpl_toolkits.mplot3d import Axes3D
  3 import numpy as np
  4 import nibabel as nib
  5 
  6 def bresenham_3d(p0, p1, thickness):
  7     '''
  8     Bresenham's Line Algorithm in 3D with adjustable thickness\n
  9     CONFIRMED
 10     '''
 11 
 12     def add_thickness(points, thickness):
 13         '''Add thickness to the line'''
 14         if thickness <= 1:
 15             return points
 16         thick_points = set()
 17         r = thickness // 2
 18         for x, y, z in points:
 19             for dx in range(-r, r + 1):
 20                 for dy in range(-r, r + 1):
 21                     for dz in range(-r, r + 1):
 22                         if dx ** 2 + dy ** 2 + dz ** 2 <= r ** 2:
 23                             thick_points.add((x + dx, y + dy, z + dz))
 24         return list(thick_points)
 25 
 26     points = []
 27     x0, y0, z0 = [int((num)) for num in p0]
 28     x1, y1, z1 = [int((num)) for num in p1]
 29     dx, dy, dz = np.abs(x1 - x0), np.abs(y1 - y0), np.abs(z1 - z0)
 30     sx, sy, sz = np.sign(x1 - x0), np.sign(y1 - y0), np.sign(z1 - z0)
 31 
 32     if dx >= dy and dx >= dz:
 33         err1, err2 = 2 * dy - dx, 2 * dz - dx
 34         while x0 != x1:
 35             points.append((x0, y0, z0))
 36             if err1 >= 0:
 37                 y0 += sy
 38                 err1 -= 2 * dx
 39             if err2 >= 0:
 40                 z0 += sz
 41                 err2 -= 2 * dx
 42             err1 += 2 * dy
 43             err2 += 2 * dz
 44             x0 += sx
 45     elif dy >= dx and dy >= dz:
 46         err1, err2 = 2 * dx - dy, 2 * dz - dy
 47         while y0 != y1:
 48             points.append((x0, y0, z0))
 49             if err1 >= 0:
 50                 x0 += sx
 51                 err1 -= 2 * dy
 52             if err2 >= 0:
 53                 z0 += sz
 54                 err2 -= 2 * dy
 55             err1 += 2 * dx
 56             err2 += 2 * dz
 57             y0 += sy
 58     else:
 59         err1, err2 = 2 * dy - dz, 2 * dx - dz
 60         while z0 != z1:
 61             points.append((x0, y0, z0))
 62             if err1 >= 0:
 63                 y0 += sy
 64                 err1 -= 2 * dz
 65             if err2 >= 0:
 66                 x0 += sx
 67                 err2 -= 2 * dz
 68             err1 += 2 * dy
 69             err2 += 2 * dx
 70             z0 += sz
 71     points.append((x1, y1, z1))
 72 
 73     # Apply thickness
 74     return add_thickness(points, thickness)
 75 
 76 
 77 # 创建一个图形对象
 78 fig = plt.figure()
 79 
 80 # 在图形中添加一个子图，参数为行数、列数和子图索引
 81 ax = fig.add_subplot(1, 1, 1, projection='3d')
 82 
 83 # 设置坐标轴的标签
 84 ax.set_xlabel('X')
 85 ax.set_ylabel('Y')
 86 ax.set_zlabel('Z')
 87 
 88 data = np.zeros((50, 50, 50))
 89 
 90 # 定义两个三维空间中的点
 91 p0 = (20, 20, 20)  # 起始点
 92 p1 = (30, 30, 30)  # 结束点
 93 
 94 # 定义线宽
 95 thickness = 25
 96 
 97 # 调用函数
 98 line_with_thickness = bresenham_3d(p0, p1, thickness)
 99 
100 # 打印结果
101 for point in line_with_thickness:
102     data[point[0], point[1], point[2]] = 1
103 
104 # 创建一个新的Nifti1Image对象
105 solid_image = nib.Nifti1Image(data, np.eye(4))
106 
107 # 保存为nii.gz文件
108 nib.save(solid_image, 't25.nii.gz')
109 
110 # 绘制体素图
111 ax.voxels(data)
112 
113 # 显示图形
114 plt.show()