Open3d计算点云切片

computePlane()给出三个点计算平面方程参数
设通过P1，P2，P3三点的平面方程为A(x - x1) + B(y - y1) + C(z - z1) = 0 。化简为一般式：Ax + By + Cz + D = 0。
将P1(x1,y1,z1)点数值代入方程Ax + By + Cz + D = 0。
即可得到：Ax1 + By 1+ Cz1 + D = 0。
化简得D = -(A * x1 + B * y1 + C * z1)。
则可以根据P1(x1,y1,z1)，P2(x2,y2,z2)，P3(x3,y3,z3)三点坐标分别求得A、B、C的值，如下：
A = (y3 - y1)(z3 - z1) - (z2 -z1)(y3 - y1);
B = (x3 - x1)(z2 - z1) - (x2 - x1)(z3 - z1);
C = (x2 - x1)(y3 - y1) - (x3 - x1)(y2 - y1);
又D = -(A * x1 + B * y1 + C * z1)，所以可以求得D的值。

import open3d as o3d
import  numpy as np
#给定三个空间点，计算平面表达式Ax+By+Cz+D=0
def computePlane(a,b,c):
        x1 = a[0]
        x2 = b[0]
        x3 = c[0]
        y1 = a[1]
        y2 = b[1]
        y3 = c[1]
        z1 = a[2]
        z2 = b[2]
        z3 = c[2]
        A = y1 * (z2 - z3) + y2 * (z3 - z1) + y3 * (z1 - z2)
        B = z1 * (x2 - x3) + z2 * (x3 - x1) + z3 * (x1 - x2)
        C = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
        D = -(x1*(y2*z3 - y3 * z2) + x2 * (y3*z1 - y1 * z3) + x3 * (y1*z2 - y2 * z1))
        return [A,B,C,D]
if __name__ == '__main__':
   ply=o3d.io.read_point_cloud("data/bunny.ply")
   #点个数
   num_points = len(ply.points)
   #计算平面参数
   A, B, C, D = computePlane([0,0,0], [0,0,1], [0,1,0])
   
   #两个平面
   d1=D+0.01
   d2=D-0.01
   idx=[]
   for i,point in enumerate(ply.points):
       #处在两个平面之间 (Ax+By+Cz+D1)  *   （Ax+By+Cz+D2）  <    0
        if (A*point[0]+B*point[1]+C*point[2]+d1)*(A*point[0]+B*point[1]+C*point[2]+d2)<0:#符号相反
                idx.append(i)
    #提取切片的点
   slice= ply.select_by_index(idx)
   slice.paint_uniform_color([1,0,0])
   #提取剩余的点
   rest_points= ply.select_by_index(idx,invert=True)
   rest_points.paint_uniform_color([0,0,1])
   o3d.visualization.draw_geometries([slice,rest_points])
   o3d.io.write_point_cloud("4.ply",slice)
   o3d.io.write_point_cloud("5.ply",rest_points)