init

crop-demo-0.py

@@ -0,0 +1,24 @@
+# ----------------------------------------------------------------------------
+# -                        Open3D: www.open3d.org                            -
+# ----------------------------------------------------------------------------
+# Copyright (c) 2018-2023 www.open3d.org
+# SPDX-License-Identifier: MIT
+# ----------------------------------------------------------------------------
+
+import open3d as o3d
+
+if __name__ == "__main__":
+    print("Load a ply point cloud, crop it, and render it")
+    sample_ply_data = o3d.data.DemoCropPointCloud()
+    pcd = o3d.io.read_point_cloud(sample_ply_data.point_cloud_path)
+    vol = o3d.visualization.read_selection_polygon_volume(
+        sample_ply_data.cropped_json_path)
+    chair = vol.crop_point_cloud(pcd)
+    # Flip the pointclouds, otherwise they will be upside down.
+    pcd.transform([[1, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]])
+    chair.transform([[1, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]])
+
+    print("Displaying original pointcloud ...")
+    o3d.visualization.draw([pcd])
+    print("Displaying cropped pointcloud")
+    o3d.visualization.draw([chair])

crop-demo-1.py

@@ -0,0 +1,56 @@
+import json
+import numpy as np
+import open3d as o3d
+
+CUBOID_EXTENT_METERS = 200
+
+METERS_BELOW_START = 5
+METERS_ABOVE_START = 30
+
+def main():
+  ## Point Cloud
+  points = np.array([
+    ## These points lie inside the cuboid
+    [-2770.94365061042, 722.0595600050154, -20.004812609192445],
+    [-2755.94365061042, 710.0595600050154, -20.004812609192445],
+    [-2755.94365061042, 710.0595600050154, -15.004812609192445],
+
+    ## These points lie outside the cuboid
+    [-2755.94365061042 + CUBOID_EXTENT_METERS, 710.0595600050154, -15.004812609192445],
+    [-2755.94365061042, 710.0595600050154 + CUBOID_EXTENT_METERS, -15.004812609192445],
+  ]).reshape([-1, 3])
+
+  point_cloud = o3d.geometry.PointCloud()
+  point_cloud.points = o3d.utility.Vector3dVector(points)
+
+  ## Start point here corresponds to an ego vehicle position start in a point cloud
+  start_position = {'x': -2755.94365061042, 'y': 722.0595600050154, 'z': -20.004812609192445}
+  cuboid_points = getCuboidPoints(start_position)
+
+  points = o3d.utility.Vector3dVector(cuboid_points)
+  oriented_bounding_box = o3d.geometry.OrientedBoundingBox.create_from_points(points)
+  point_cloud_crop = point_cloud.crop(oriented_bounding_box)
+
+  # View original point cloud with the cuboid, all 5 points present
+  o3d.visualization.draw_geometries([point_cloud, oriented_bounding_box])
+
+  # View cropped point cloud with the cuboid, only 3 points present
+  o3d.visualization.draw_geometries([point_cloud_crop, oriented_bounding_box])
+
+def getCuboidPoints(start_position):
+  return np.array([
+    # Vertices Polygon1
+    [start_position['x'] + (CUBOID_EXTENT_METERS / 2), start_position['y'] + (CUBOID_EXTENT_METERS / 2), start_position['z'] + METERS_ABOVE_START], # face-topright
+    [start_position['x'] - (CUBOID_EXTENT_METERS / 2), start_position['y'] + (CUBOID_EXTENT_METERS / 2), start_position['z'] + METERS_ABOVE_START], # face-topleft
+    [start_position['x'] - (CUBOID_EXTENT_METERS / 2), start_position['y'] - (CUBOID_EXTENT_METERS / 2), start_position['z'] + METERS_ABOVE_START], # rear-topleft
+    [start_position['x'] + (CUBOID_EXTENT_METERS / 2), start_position['y'] - (CUBOID_EXTENT_METERS / 2), start_position['z'] + METERS_ABOVE_START], # rear-topright
+
+    # Vertices Polygon 2
+    [start_position['x'] + (CUBOID_EXTENT_METERS / 2), start_position['y'] + (CUBOID_EXTENT_METERS / 2), start_position['z'] - METERS_BELOW_START],
+    [start_position['x'] - (CUBOID_EXTENT_METERS / 2), start_position['y'] + (CUBOID_EXTENT_METERS / 2), start_position['z'] - METERS_BELOW_START],
+    [start_position['x'] - (CUBOID_EXTENT_METERS / 2), start_position['y'] - (CUBOID_EXTENT_METERS / 2), start_position['z'] - METERS_BELOW_START],
+    [start_position['x'] + (CUBOID_EXTENT_METERS / 2), start_position['y'] - (CUBOID_EXTENT_METERS / 2), start_position['z'] - METERS_BELOW_START],
+  ]).astype("float64") 
+
+if __name__ == '__main__':
+  main()

crop-demo-2.py

@@ -0,0 +1,55 @@
+import open3d as o3d
+import numpy as np
+
+pcd = o3d.io.read_point_cloud("../laser1.ply")
+# ply_point_cloud = o3d.data.PLYPointCloud()
+# pcd = o3d.io.read_point_cloud(ply_point_cloud.path)
+
+corners = np.array([
+  [ 10, -3.21384387,  0.30217625,],
+  [ 10, -1.13804348,  0.29917539,],
+  [ 7.69983939, -1.16651864,  0.30329364,],
+  [ 7.67473496, -3.24231903,  0.3062945, ],
+  [ 10, -3.21276837,  1.03551451,],
+  [ 10, -1.13696798,  1.03251366,],
+  [ 7.69856999, -1.16544314,  1.03663191,],
+  [ 7.67346556, -3.24124353,  1.03963277,],
+])
+
+# Convert the corners array to have type float64
+bounding_polygon = corners.astype("float64")
+
+# Create a SelectionPolygonVolume
+vol = o3d.visualization.SelectionPolygonVolume()
+
+# You need to specify what axis to orient the polygon to.
+# I choose the "Y" axis. I made the max value the maximum Y of
+# the polygon vertices and the min value the minimum Y of the
+# polygon vertices.
+vol.orthogonal_axis = "Y"
+vol.axis_max = np.max(bounding_polygon[:, 1])
+vol.axis_min = np.min(bounding_polygon[:, 1])
+
+# Set all the Y values to 0 (they aren't needed since we specified what they
+# should be using just vol.axis_max and vol.axis_min).
+bounding_polygon[:, 1] = 0
+
+# Convert the np.array to a Vector3dVector
+vol.bounding_polygon = o3d.utility.Vector3dVector(bounding_polygon)
+
+# Crop the point cloud using the Vector3dVector
+cropped_pcd = vol.crop_point_cloud(pcd)
+
+# Get a nice looking bounding box to display around the newly cropped point cloud
+# (This part is optional and just for display purposes)
+bounding_box = cropped_pcd.get_axis_aligned_bounding_box()
+bounding_box.color = (1, 0, 0)
+
+# Draw the newly cropped PCD and bounding box
+# o3d.visualization.draw_geometries([cropped_pcd, bounding_box],
+#                                   zoom=2,
+#                                   front=[5, -2, 0.5],
+#                                   lookat=[7.67473496, -3.24231903,  0.3062945],
+#                                   up=[1.0, 0.0, 0.0])
+# o3d.visualization.draw_geometries([pcd, bounding_box])
+o3d.visualization.draw([pcd, bounding_box],)

doc-open3d.md

@@ -0,0 +1,107 @@
+# Geometry
+
+## point cloud!!!
+open3d.data.PointCloud内部数据：
+```
+ |  Data descriptors defined here:
+ |
+ |  colors
+ |      ``float64`` array of shape ``(num_points, 3)``, range ``[0, 1]`` , use ``numpy.asarray()`` to access data: RGB colors of points.
+ |
+ |  covariances
+ |      ``float64`` array of shape ``(num_points, 3, 3)``, use ``numpy.asarray()`` to access data: Points covariances.
+ |
+ |  normals
+ |      ``float64`` array of shape ``(num_points, 3)``, use ``numpy.asarray()`` to access data: Points normals.
+ |
+ |  points
+ |      ``float64`` array of shape ``(num_points, 3)``, use ``numpy.asarray()`` to access data: Points coordinates.
+ |
+```
+
+### Voxel downsampling!!!
+可用于减少计算量。
+
+### crop point cloud!!!
+用一个立体的多边形polygon（用一个json数据描述）对点云做裁剪。这个json数据的规则呢？
+
+### Convex hull!!!
+
+### DBSCAN clustering
+对点云聚类，给每个点加上一个标签属性，-1表示噪点。参数:同类别中相邻点的距离、一个类别最少包含多少个点。
+
+计算量大！
+
+### Plane segmentation
+找到潜在的平面，把点云中各个点分类到不同的平面上。系用最小二乘法近似。
+
+### planar patch detection
+类似上一条，但找到的不是无限的平面而是平的有边界的“片”。
+
+### Hidden point removal
+在某个视角查看点云时，把按理说会在点云所代表的物体的背面的那些点去掉。
+
+## mesh
+
+### mesh filtering
+对mesh表面做滤波，效果可以是让mesh表面变得光滑一些。
+
+### Sampling
+对mesh表面采样形成点云。
+
+### Mesh subdivision
+把各个三角形拆成更小的三角形，从而让mesh表面更圆滑。
+
+### Mesh simplification
+把mesh表面三角形变大，从而让mesh表面更粗糙更简单。
+
+### Connected components
+把各个三角形按连通与否分组。进而能够知道每个组的大小。
+
+## Point cloud outlier removal
+去除点云中的噪声
+
+## Voxelization
+点云和mesh都可以素格化。
+
+进而可以知道给定点是否在素格以内。
+
+### voxel carving
+用点云和mesh得到的素格都在物体表面。用深度图则可以carve出实心的一堆素格。
+
+## Surface reconstruction
+就是从点云得到mesh。
+
+### alpha shapes！！！！
+点云内部也可以挖空。
+
+得到的mesh可能不封闭。
+
+### ball pivoting！！！
+要求点云有法线（normals）信息。
+
+封闭？？？？
+
+### Poisson surface reconstruction
+要求点云有法线（normals）信息。
+
+得到的mesh可能不封闭。
+
+### Normal estimation
+用于给点云加上法线信息。可以用一个方法使各点的法线协调起来。
+
+## Ray Casting
+创建一个场景，里边放些几何物，设置一个光线的远点和方向，能得到光线与几何物的交汇点。
+
+## Distance Queries
+用于知道给定的一个点到mesh的最短距离，以及这个点是否在mesh内部。！！！！
+
+## python interface
+用help(o3d.xxx.yyy)获取某个API的文档。！！！
+
+## Working with Numpy
+
+# Visualization
+
+## interactive_visualization.py
+

mission-voxelization.py

@@ -0,0 +1,117 @@
+# 如此设定每个点占据的体积：确保比例为ratio_factor的点周围至少有neighbor_point_min_num个点和自己连通。
+
+import open3d as o3d
+import numpy as np
+import math
+import sys
+import random
+
+ratio_factor = 0.9
+neighbor_point_min_num = 8
+
+# 读取点云文件
+pcd = o3d.io.read_point_cloud("./data/PointClound_01.ply")
+print(pcd)
+# o3d.visualization.draw([pcd])
+
+# 原始点数量
+pcd_size = np.asarray(pcd.points).size / 3
+print('原始点数量：', pcd_size)
+
+# # 粗暴降采样
+# down_sampling_size = 1
+# pcd_down = pcd.voxel_down_sample(voxel_size=down_sampling_size)
+# pcd_down_size = np.asarray(pcd_down.points).size / 3
+
+# # 酌情降采样，直到点数量降到阈值以下
+# down_sampling_size = 0.01
+# pcd_down_size = float('inf')
+# while(pcd_down_size > 1000):
+#   down_sampling_size = down_sampling_size + 0.02
+#   pcd_down = pcd.voxel_down_sample(voxel_size=down_sampling_size)
+#   pcd_down_size = np.asarray(pcd_down.points).size / 3
+#   if down_sampling_size > 10:
+#     break
+
+# print('降采样后点数量：', pcd_down_size)
+# print('降采样尺度：', down_sampling_size)
+# o3d.visualization.draw([pcd])
+
+# bbox情况
+# obb = pcd.get_oriented_bounding_box()
+# print('外包盒子体积：', obb.volume())
+
+# # fit to unit cube
+# pcd.scale(1 / np.max(pcd.get_max_bound() - pcd.get_min_bound()),
+#           center=pcd.get_center())
+# o3d.visualization.draw_geometries([pcd])
+
+# 建立kd tree
+print('build kdtree...')
+pcd_tree = o3d.geometry.KDTreeFlann(pcd)
+print('build over')
+
+# 对每个采样点，记录其与周围点的距离
+sample_num = 0
+if pcd_size <= 10000:
+  sample_num = pcd_size
+else:
+  sample_num = 10000
+neighbor_dist_list = np.zeros([sample_num], dtype = float)
+for sampleIdx in range(0, math.floor(sample_num)):
+  sampleIdxInPcd = random.randrange(0, pcd_size - 1, 1)
+  neighbor_num = neighbor_point_min_num + 1 # 真正的最近的那个点是自己
+  [k, neighbor_idx_list, _] = pcd_tree.search_knn_vector_3d(pcd.points[sampleIdxInPcd], neighbor_num)
+  pointSample = np.asarray(pcd.points)[sampleIdxInPcd]
+  for neighborIdx in range(1, neighbor_num):
+    pointNeighbor = np.asarray(pcd.points)[neighbor_idx_list[neighborIdx]]
+    new_val = math.sqrt(math.pow(pointNeighbor[0] - pointSample[0], 2) + math.pow(pointNeighbor[1] - pointSample[1], 2) + math.pow(pointNeighbor[2] - pointSample[2], 2))
+    if new_val > neighbor_dist_list[sampleIdx]:
+      neighbor_dist_list[sampleIdx] = new_val
+  # print(neighbor_dist_list[sampleIdx])
+  # if i < 100:
+  #   print('-------------')
+  #   print('iterate', i)
+  #   print(pointNeighbor, pointSample)
+  #   print(pointNeighbor[0] - pointSample[0])
+  #   print(math.pow(pointNeighbor[0] - pointSample[0], 2))
+  #   print(math.pow(pointNeighbor[0] - pointSample[0], 2) + math.pow(pointNeighbor[1] - pointSample[1], 2) + math.pow(pointNeighbor[2] - pointSample[2], 2))
+  #   print(neighbor_dist_list[i])
+# exit(0)
+
+# 确定素格尺寸
+voxel_size = 0
+neighbor_dist_min = np.min(neighbor_dist_list)
+neighbor_dist_max = np.max(neighbor_dist_list)
+print('距离最大值：', np.max(neighbor_dist_list))
+print('距离最小值：', np.min(neighbor_dist_list))
+if neighbor_dist_min == neighbor_dist_max:
+  voxel_size = neighbor_dist_max
+else:
+  binNum = 1000
+  histogram = np.zeros((binNum), dtype = int) 
+  for sampleIdx in range(0, sample_num):
+    belongBinIdx = math.floor((neighbor_dist_list[sampleIdx] - neighbor_dist_min) / (neighbor_dist_max - neighbor_dist_min) * binNum)
+    if belongBinIdx == binNum:
+      belongBinIdx = binNum - 1
+    histogram[belongBinIdx] = histogram[belongBinIdx] + 1
+
+  countThreshold = sample_num * ratio_factor
+  countNow = 0
+  for sampleIdx in range(0, binNum):
+    print('-----------')
+    print(histogram[sampleIdx], 'bin content value')
+    countNow = countNow + histogram[sampleIdx]
+    print(countNow, 'current total')
+    if countNow >= countThreshold:
+      voxel_size = neighbor_dist_min + (neighbor_dist_max - neighbor_dist_min) / binNum * (sampleIdx + 1)
+      break
+print('voxel_size:', voxel_size)
+
+voxel_grid = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd, voxel_size=voxel_size)
+print('voxel_grid: ', voxel_grid)
+total_volume = np.asarray(voxel_grid.get_voxels()).size * (math.pow(voxel_size, 3))
+print('点云占据体积：', total_volume)
+o3d.visualization.draw_geometries([voxel_grid])
+
+sys.exit(0)

mission_alpha_shaped_mesh.py

@@ -0,0 +1,32 @@
+import open3d as o3d
+import numpy as np
+
+# --------------------------- 加载点云 ---------------------------
+print("->正在加载点云... ")
+# ply = o3d.io.read_point_cloud("../laser1.ply")
+pcd = o3d.io.read_point_cloud("./data/PointClound_01.ply")
+
+# 展示点云、打印点云
+print('pcd data: ~~~~~~~~~~~~~~~~~~~')
+print(pcd)
+# print(np.asarray(pcd.points))
+# print(np.asarray(pcd.normals))
+# print(np.asarray(pcd.covariances))
+# print(np.asarray(pcd.colors))
+print('end of pcd data ~~~~~~~~~~~~~~~~~~~')
+axes = o3d.geometry.TriangleMesh.create_coordinate_frame()
+o3d.visualization.draw([pcd, axes])
+
+# 降采样
+pcd = pcd.voxel_down_sample(voxel_size=0.1)
+o3d.visualization.draw([pcd])
+
+# ------------------------- Alpha shapes -----------------------
+alpha = 20
+print(f"alpha={alpha:.3f}")
+mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(pcd, alpha)
+mesh.compute_vertex_normals()
+o3d.visualization.draw_geometries([mesh], mesh_show_back_face=True)
+volume = mesh.get_volume()      # 计算体积
+print("体积为：", volume)
+# ==============================================================

mission_convex_hull.py

@@ -0,0 +1,58 @@
+import open3d as o3d
+import numpy as np
+import math
+from scipy.spatial import Delaunay
+import math
+from functools import reduce
+import sys
+import copy
+
+pcd = o3d.io.read_point_cloud("data/stockpile.ply")
+# pcd = o3d.io.read_point_cloud("../laser1.ply")
+# print(pcd)  # 输出点云点的个数
+o3d.visualization.draw([pcd])
+print("Paint pointcloud")
+
+# 找出地面所在平面
+plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,
+                                         ransac_n=3,
+                                         num_iterations=10000)
+
+# 把地面上的点和其他点分开
+[a, b, c, d] = plane_model
+plane_pcd = pcd.select_by_index(inliers)
+plane_pcd.paint_uniform_color([1.0, 0, 0])
+stockpile_pcd = pcd.select_by_index(inliers, invert=True)
+stockpile_pcd.paint_uniform_color([0, 0, 1.0])
+# o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])
+
+# 移动到原点
+plane_pcd = plane_pcd.translate((0,0,d/c))
+stockpile_pcd = stockpile_pcd.translate((0,0,d/c))
+
+# 旋转至与坐标轴对齐
+cos_theta = c / math.sqrt(a**2 + b**2 + c**2)
+sin_theta = math.sqrt((a**2+b**2)/(a**2 + b**2 + c**2))
+u_1 = b / math.sqrt(a**2 + b**2 )
+u_2 = -a / math.sqrt(a**2 + b**2)
+rotation_matrix = np.array([[cos_theta + u_1**2 * (1-cos_theta), u_1*u_2*(1-cos_theta), u_2*sin_theta],
+                            [u_1*u_2*(1-cos_theta), cos_theta + u_2**2*(1- cos_theta), -u_1*sin_theta],
+                            [-u_2*sin_theta, u_1*sin_theta, cos_theta]])
+plane_pcd.rotate(rotation_matrix)
+stockpile_pcd.rotate(rotation_matrix)
+# o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])
+
+# 去除噪点
+cl, ind = stockpile_pcd.remove_statistical_outlier(nb_neighbors=50,
+                                                    std_ratio=5.0)
+stockpile_pcd = stockpile_pcd.select_by_index(ind)
+
+o3d.visualization.draw([stockpile_pcd])
+ 
+hull, idx = stockpile_pcd.compute_convex_hull()
+hull.paint_uniform_color([1, 0.7, 0])  # 给凸包渲染颜色
+area = hull.get_surface_area()  # 计算表面积
+volume = hull.get_volume()      # 计算体积
+print("表面积为：", area)
+print("体积为：", volume)
+o3d.visualization.draw_geometries([stockpile_pcd, hull])

mission_stock_pile_triangularization.py

@@ -0,0 +1,83 @@
+import open3d as o3d
+import numpy as np
+import math
+from scipy.spatial import Delaunay
+import math
+from functools import reduce
+import sys
+
+pcd = o3d.io.read_point_cloud("data/stockpile.ply")
+print(pcd)
+
+axes = o3d.geometry.TriangleMesh.create_coordinate_frame()
+
+# 找出地面所在平面
+plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,
+                                         ransac_n=3,
+                                         num_iterations=10000)
+
+# 把地面上的点和其他点分开
+[a, b, c, d] = plane_model
+plane_pcd = pcd.select_by_index(inliers)
+plane_pcd.paint_uniform_color([1.0, 0, 0])
+stockpile_pcd = pcd.select_by_index(inliers, invert=True)
+stockpile_pcd.paint_uniform_color([0, 0, 1.0])
+o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])
+
+# 移动到原点
+plane_pcd = plane_pcd.translate((0,0,d/c))
+stockpile_pcd = stockpile_pcd.translate((0,0,d/c))
+
+# 旋转至与坐标轴对齐
+cos_theta = c / math.sqrt(a**2 + b**2 + c**2)
+sin_theta = math.sqrt((a**2+b**2)/(a**2 + b**2 + c**2))
+u_1 = b / math.sqrt(a**2 + b**2 )
+u_2 = -a / math.sqrt(a**2 + b**2)
+rotation_matrix = np.array([[cos_theta + u_1**2 * (1-cos_theta), u_1*u_2*(1-cos_theta), u_2*sin_theta],
+                            [u_1*u_2*(1-cos_theta), cos_theta + u_2**2*(1- cos_theta), -u_1*sin_theta],
+                            [-u_2*sin_theta, u_1*sin_theta, cos_theta]])
+plane_pcd.rotate(rotation_matrix)
+stockpile_pcd.rotate(rotation_matrix)
+# o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])
+
+# 去除噪点
+cl, ind = stockpile_pcd.remove_statistical_outlier(nb_neighbors=30,
+                                                    std_ratio=2.0)
+stockpile_pcd = stockpile_pcd.select_by_index(ind)
+# o3d.visualization.draw_geometries([stockpile_pcd])
+
+# # 降采样
+# downpcd = stockpile_pcd.voxel_down_sample(voxel_size=0.05)
+# print(downpcd)
+downpcd = stockpile_pcd
+
+# 三角剖分
+xyz = np.asarray(downpcd.points)
+xy_catalog = []
+for point in xyz:
+    xy_catalog.append([point[0], point[1]])
+tri = Delaunay(np.array(xy_catalog))
+
+surface = o3d.geometry.TriangleMesh()
+surface.vertices = o3d.utility.Vector3dVector(xyz)
+surface.triangles = o3d.utility.Vector3iVector(tri.simplices)
+surface.paint_uniform_color([0, 0, 1.0])
+# o3d.visualization.draw_geometries([surface], mesh_show_wireframe=True)
+
+def get_triangles_vertices(triangles, vertices):
+    triangles_vertices = []
+    for triangle in triangles:
+        new_triangles_vertices = [vertices[triangle[0]], vertices[triangle[1]], vertices[triangle[2]]]
+        triangles_vertices.append(new_triangles_vertices)
+    return np.array(triangles_vertices)
+
+def volume_under_triangle(triangle):
+    p1, p2, p3 = triangle
+    x1, y1, z1 = p1
+    x2, y2, z2 = p2
+    x3, y3, z3 = p3
+    return abs((z1+z2+z3)*(x1*y2-x2*y1+x2*y3-x3*y2+x3*y1-x1*y3)/6)
+
+print('sdfjsfksldfj')
+volume = reduce(lambda a, b:  a + volume_under_triangle(b), get_triangles_vertices(surface.triangles, surface.vertices), 0)
+print(f"The volume of the stockpile is: {round(volume, 4)} m3")

mission_stock_pile_triangularization_true_data.py

@@ -0,0 +1,87 @@
+import open3d as o3d
+import numpy as np
+import math
+from scipy.spatial import Delaunay
+import math
+from functools import reduce
+import sys
+
+# pcd = o3d.io.read_point_cloud("data/2box-SG-t-MnQ0zoEQIBX.ply")
+# pcd = o3d.io.read_point_cloud("data/4box-1-SG-t-AsSRO3iA0XT.ply")
+# pcd = o3d.io.read_point_cloud("data/4box-2-SG-t-cp3kSewxMKi.ply")
+pcd = o3d.io.read_point_cloud("data/4box-3-SG-t-kuuv0moEP0C.ply")
+print(pcd)
+
+axes = o3d.geometry.TriangleMesh.create_coordinate_frame()
+
+# 找出地面所在平面
+plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,
+                                         ransac_n=3,
+                                         num_iterations=10000)
+
+
+# # 把地面上的点和其他点分开
+# [a, b, c, d] = plane_model
+# plane_pcd = pcd.select_by_index(inliers)
+# plane_pcd.paint_uniform_color([1.0, 0, 0])
+# stockpile_pcd = pcd.select_by_index(inliers, invert=True)
+# stockpile_pcd.paint_uniform_color([0, 0, 1.0])
+# o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])
+
+# 移动到原点
+plane_pcd = plane_pcd.translate((0,0,d/c))
+stockpile_pcd = stockpile_pcd.translate((0,0,d/c))
+
+# 旋转至与坐标轴对齐
+cos_theta = c / math.sqrt(a**2 + b**2 + c**2)
+sin_theta = math.sqrt((a**2+b**2)/(a**2 + b**2 + c**2))
+u_1 = b / math.sqrt(a**2 + b**2 )
+u_2 = -a / math.sqrt(a**2 + b**2)
+rotation_matrix = np.array([[cos_theta + u_1**2 * (1-cos_theta), u_1*u_2*(1-cos_theta), u_2*sin_theta],
+                            [u_1*u_2*(1-cos_theta), cos_theta + u_2**2*(1- cos_theta), -u_1*sin_theta],
+                            [-u_2*sin_theta, u_1*sin_theta, cos_theta]])
+plane_pcd.rotate(rotation_matrix)
+stockpile_pcd.rotate(rotation_matrix)
+# o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])
+
+# 去除噪点
+cl, ind = stockpile_pcd.remove_statistical_outlier(nb_neighbors=30,
+                                                    std_ratio=2.0)
+stockpile_pcd = stockpile_pcd.select_by_index(ind)
+# o3d.visualization.draw_geometries([stockpile_pcd])
+
+# # 降采样
+# downpcd = stockpile_pcd.voxel_down_sample(voxel_size=0.05)
+# print(downpcd)
+downpcd = stockpile_pcd
+
+# 三角剖分
+xyz = np.asarray(downpcd.points)
+xy_catalog = []
+for point in xyz:
+    xy_catalog.append([point[0], point[1]])
+tri = Delaunay(np.array(xy_catalog))
+
+surface = o3d.geometry.TriangleMesh()
+surface.vertices = o3d.utility.Vector3dVector(xyz)
+surface.triangles = o3d.utility.Vector3iVector(tri.simplices)
+surface.paint_uniform_color([0, 0, 1.0])
+# o3d.visualization.draw_geometries([surface], mesh_show_wireframe=True)
+
+def get_triangles_vertices(triangles, vertices):
+    triangles_vertices = []
+    for triangle in triangles:
+        new_triangles_vertices = [vertices[triangle[0]], vertices[triangle[1]], vertices[triangle[2]]]
+        triangles_vertices.append(new_triangles_vertices)
+    return np.array(triangles_vertices)
+
+def volume_under_triangle(triangle):
+    p1, p2, p3 = triangle
+    x1, y1, z1 = p1
+    x2, y2, z2 = p2
+    x3, y3, z3 = p3
+    return abs((z1+z2+z3)*(x1*y2-x2*y1+x2*y3-x3*y2+x3*y1-x1*y3)/6)
+
+print('sdfjsfksldfj')
+volume = reduce(lambda a, b:  a + volume_under_triangle(b), get_triangles_vertices(surface.triangles, surface.vertices), 0)
+print(f"The volume of the stockpile is: {round(volume, 4)} m3")

mission_stock_pile_voxel.py

@@ -0,0 +1,138 @@
+import open3d as o3d
+import numpy as np
+import math
+import math
+from functools import reduce
+import sys
+import copy
+
+pcd = o3d.io.read_point_cloud("data/stockpile.ply")
+print(pcd)
+
+axes = o3d.geometry.TriangleMesh.create_coordinate_frame()
+
+# 找出地面所在平面
+plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,
+                                         ransac_n=3,
+                                         num_iterations=10000)
+
+# 把地面上的点和其他点分开
+[a, b, c, d] = plane_model
+plane_pcd = pcd.select_by_index(inliers)
+plane_pcd.paint_uniform_color([1.0, 0, 0])
+stockpile_pcd = pcd.select_by_index(inliers, invert=True)
+# stockpile_pcd = pcd
+stockpile_pcd.paint_uniform_color([0, 0, 1.0])
+# o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])
+
+# 移动到原点
+plane_pcd = plane_pcd.translate((0,0,d/c))
+stockpile_pcd = stockpile_pcd.translate((0,0,d/c))
+
+# 旋转至与坐标轴对齐
+cos_theta = c / math.sqrt(a**2 + b**2 + c**2)
+sin_theta = math.sqrt((a**2+b**2)/(a**2 + b**2 + c**2))
+u_1 = b / math.sqrt(a**2 + b**2 )
+u_2 = -a / math.sqrt(a**2 + b**2)
+rotation_matrix = np.array([[cos_theta + u_1**2 * (1-cos_theta), u_1*u_2*(1-cos_theta), u_2*sin_theta],
+                            [u_1*u_2*(1-cos_theta), cos_theta + u_2**2*(1- cos_theta), -u_1*sin_theta],
+                            [-u_2*sin_theta, u_1*sin_theta, cos_theta]])
+plane_pcd.rotate(rotation_matrix)
+stockpile_pcd.rotate(rotation_matrix)
+o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])
+
+# 去除噪点
+cl, ind = stockpile_pcd.remove_statistical_outlier(nb_neighbors=100,
+                                                    std_ratio=10.0)
+stockpile_pcd = stockpile_pcd.select_by_index(ind)
+print('有效点：', stockpile_pcd)
+o3d.visualization.draw_geometries([stockpile_pcd])
+
+# 降采样
+# downpcd = stockpile_pcd.voxel_down_sample(voxel_size=0.1)
+downpcd = stockpile_pcd
+print('有效点云：', np.asarray(downpcd.points))
+
+mean = 0
+meanCount = 0
+for point in np.asarray(downpcd.points):
+  mean += point[2]
+  meanCount += 1
+mean = mean / meanCount
+print('z均值：', mean)
+
+# 建立二维网格
+bbox = downpcd.get_axis_aligned_bounding_box()
+print('bbox: ', bbox)
+bin_num_1d = 100
+bin_size_x = (bbox.max_bound[0] - bbox.min_bound[0]) / bin_num_1d
+bin_size_y = (bbox.max_bound[1] - bbox.min_bound[1]) / bin_num_1d
+grid2d = np.zeros((bin_num_1d, bin_num_1d, 2))
+print('bin_size: ', bin_size_x, bin_size_y)
+
+points = np.asarray(downpcd.points)
+bin_idx_x = -1
+bin_idx_y = -1
+for point in points:
+  bin_idx_x = math.floor((point[0] - bbox.min_bound[0]) / (bbox.max_bound[0] - bbox.min_bound[0]) * bin_num_1d)
+  if bin_idx_x <= -1: # 浮点数溢出
+    bin_idx_x = 0
+  if bin_idx_x >= bin_num_1d:
+    bin_idx_x = bin_num_1d - 1
+
+  bin_idx_y = math.floor((point[1] - bbox.min_bound[1]) / (bbox.max_bound[1] - bbox.min_bound[1]) * bin_num_1d)
+  if bin_idx_y <= -1: # 浮点数溢出
+    bin_idx_y = 0
+  if bin_idx_y >= bin_num_1d:
+    bin_idx_y = bin_num_1d - 1
+
+  # print('bin before ', bin_idx_x, bin_idx_y, grid2d[bin_idx_x, bin_idx_y, 0], grid2d[bin_idx_x, bin_idx_y, 1])
+  # print('point', point)
+  grid2d[bin_idx_x, bin_idx_y, 0] += (point[2]) # 注意！认为地面高度为0！
+  grid2d[bin_idx_x, bin_idx_y, 1] += 1
+  # print('bin after  ', bin_idx_x, bin_idx_y, grid2d[bin_idx_x, bin_idx_y, 0], grid2d[bin_idx_x, bin_idx_y, 1])
+
+# 建立二维kdtree以便查询邻居，对二维网格中没有点的网格做插值
+downpcd_2d = copy.deepcopy(downpcd)
+for point in downpcd_2d.points:
+  point[2] = 0
+# print('build kdtree...')
+kd_tree_2d = o3d.geometry.KDTreeFlann(downpcd_2d)
+# print('build over')
+
+# 二维网格插值
+sum = 0
+interCount = 0
+interSum = 0
+for i in range(0, bin_num_1d):
+  for j in range(0, bin_num_1d):
+    # print('~~~~~~~~~~~~~')
+
+    if grid2d[i, j, 1] == 0:
+      # print('做插值……')
+      interCount += 1
+      [k, neighbor_idx_list, _] = kd_tree_2d.search_knn_vector_3d([bin_size_x * (i + 0.5) + bbox.min_bound[0], bin_size_y * (j + 0.5) + bbox.min_bound[1], 0], 1)
+      # print(neighbor_idx_list)
+      pointNeighbor = np.asarray(downpcd.points)[neighbor_idx_list[0]]
+      # if pointNeighbor[2] > 0.08:
+      #   print('最近点：', pointNeighbor)
+      grid2d[i, j, 0] = pointNeighbor[2]
+      grid2d[i, j, 1] = 1
+      interSum += pointNeighbor[2]
+
+      
+    if(grid2d[i, j, 1]):
+      sum += (grid2d[i, j, 0] / grid2d[i, j, 1] * bin_size_x * bin_size_y)
+    # print('本网格idx', i, j)
+    # print('本网格sum',grid2d[i, j, 0])
+    # print('本网格点数',grid2d[i, j, 1])
+    
+    # print('sum: ', sum)
+    # print('~~~~~~~~~~~~~')
+
+print('插值均值：', interSum / interCount)
+
+print('体积：', sum)
+print('插值次数：', interCount)
+
+exit(0)

mission_stock_pile_voxel_true_data.py

@@ -0,0 +1,147 @@
+import open3d as o3d
+import numpy as np
+import math
+import math
+from functools import reduce
+import sys
+import copy
+
+# pcd = o3d.io.read_point_cloud("data/2box-SG-t-MnQ0zoEQIBX.ply")
+# pcd = o3d.io.read_point_cloud("data/4box-1-SG-t-AsSRO3iA0XT.ply")
+# pcd = o3d.io.read_point_cloud("data/4box-2-SG-t-cp3kSewxMKi.ply")
+pcd = o3d.io.read_point_cloud("data/4box-3-SG-t-kuuv0moEP0C.ply")
+print(pcd)
+assert (pcd.has_normals())
+
+axes = o3d.geometry.TriangleMesh.create_coordinate_frame()
+
+# 找出地面所在平面
+plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,
+                                         ransac_n=3,
+                                         num_iterations=10000)
+
+# 把地面上的点和其他点分开
+[a, b, c, d] = plane_model
+plane_pcd = pcd.select_by_index(inliers)
+plane_pcd.paint_uniform_color([1.0, 0, 0])
+stockpile_pcd = pcd.select_by_index(inliers, invert=True)
+# stockpile_pcd = pcd
+stockpile_pcd.paint_uniform_color([0, 0, 1.0])
+o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])
+
+# 移动到原点
+plane_pcd = plane_pcd.translate((0,0,d/c))
+stockpile_pcd = stockpile_pcd.translate((0,0,d/c))
+
+# 旋转至与坐标轴对齐
+cos_theta = c / math.sqrt(a**2 + b**2 + c**2)
+sin_theta = math.sqrt((a**2+b**2)/(a**2 + b**2 + c**2))
+u_1 = b / math.sqrt(a**2 + b**2 )
+u_2 = -a / math.sqrt(a**2 + b**2)
+rotation_matrix = np.array([[cos_theta + u_1**2 * (1-cos_theta), u_1*u_2*(1-cos_theta), u_2*sin_theta],
+                            [u_1*u_2*(1-cos_theta), cos_theta + u_2**2*(1- cos_theta), -u_1*sin_theta],
+                            [-u_2*sin_theta, u_1*sin_theta, cos_theta]])
+plane_pcd.rotate(rotation_matrix)
+stockpile_pcd.rotate(rotation_matrix)
+o3d.visualization.draw_geometries([plane_pcd, stockpile_pcd, axes])
+
+# 去除噪点
+cl, ind = stockpile_pcd.remove_statistical_outlier(nb_neighbors=80,
+                                                    std_ratio=0.1)
+stockpile_pcd = stockpile_pcd.select_by_index(ind)
+print('有效点：', stockpile_pcd)
+o3d.visualization.draw_geometries([stockpile_pcd])
+
+# 降采样
+# downpcd = stockpile_pcd.voxel_down_sample(voxel_size=0.1)
+downpcd = stockpile_pcd
+print('有效点云：', np.asarray(downpcd.points))
+
+mean = 0
+meanCount = 0
+for point in np.asarray(downpcd.points):
+  mean += point[2]
+  meanCount += 1
+mean = mean / meanCount
+print('z均值：', mean)
+
+# 建立二维网格
+bbox = downpcd.get_axis_aligned_bounding_box()
+print('bbox: ', bbox)
+bin_num_1d = 200
+bin_size_x = (bbox.max_bound[0] - bbox.min_bound[0]) / bin_num_1d
+bin_size_y = (bbox.max_bound[1] - bbox.min_bound[1]) / bin_num_1d
+grid2d = np.zeros((bin_num_1d, bin_num_1d, 2))
+print('bin_size: ', bin_size_x, bin_size_y)
+
+points = np.asarray(downpcd.points)
+bin_idx_x = -1
+bin_idx_y = -1
+for point in points:
+  bin_idx_x = math.floor((point[0] - bbox.min_bound[0]) / (bbox.max_bound[0] - bbox.min_bound[0]) * bin_num_1d)
+  if bin_idx_x <= -1: # 浮点数溢出
+    bin_idx_x = 0
+  if bin_idx_x >= bin_num_1d:
+    bin_idx_x = bin_num_1d - 1
+
+  bin_idx_y = math.floor((point[1] - bbox.min_bound[1]) / (bbox.max_bound[1] - bbox.min_bound[1]) * bin_num_1d)
+  if bin_idx_y <= -1: # 浮点数溢出
+    bin_idx_y = 0
+  if bin_idx_y >= bin_num_1d:
+    bin_idx_y = bin_num_1d - 1
+
+  # print('bin before ', bin_idx_x, bin_idx_y, grid2d[bin_idx_x, bin_idx_y, 0], grid2d[bin_idx_x, bin_idx_y, 1])
+  # print('point', point)
+  grid2d[bin_idx_x, bin_idx_y, 0] += (point[2]) # 注意！认为地面高度为0！
+  grid2d[bin_idx_x, bin_idx_y, 1] += 1
+  # print('bin after  ', bin_idx_x, bin_idx_y, grid2d[bin_idx_x, bin_idx_y, 0], grid2d[bin_idx_x, bin_idx_y, 1])
+
+# 建立二维kdtree以便查询邻居，对二维网格中没有点的网格做插值
+downpcd_2d = copy.deepcopy(downpcd)
+for point in downpcd_2d.points:
+  point[2] = 0
+# print('build kdtree...')
+kd_tree_2d = o3d.geometry.KDTreeFlann(downpcd_2d)
+# print('build over')
+
+# 二维网格体积累积
+sum = 0
+interCount = 0
+interSum = 0
+
+for i in range(0, bin_num_1d):
+  for j in range(0, bin_num_1d):
+    # print('~~~~~~~~~~~~~')
+
+    if grid2d[i, j, 1] == 0:
+      # print('做插值……')
+      interCount += 1
+      [k, neighbor_idx_list, _] = kd_tree_2d.search_knn_vector_3d([bin_size_x * (i + 0.5) + bbox.min_bound[0], bin_size_y * (j + 0.5) + bbox.min_bound[1], 0], 100)
+      
+      # 针对物体有垂直的侧面的情况，拿k近邻中z值最小的
+      neighbor_pcd = downpcd.select_by_index(neighbor_idx_list)
+      neighbor_pcd_pos = np.asarray(neighbor_pcd.points)
+      z_min = 100
+      for pos in neighbor_pcd_pos:
+        if pos[2] < z_min:
+          z_min = pos[2]
+      
+      grid2d[i, j, 0] = z_min
+      grid2d[i, j, 1] = 1
+      interSum += z_min
+
+    if(grid2d[i, j, 1]):
+      sum += (grid2d[i, j, 0] / grid2d[i, j, 1] * bin_size_x * bin_size_y)
+    # print('本网格idx', i, j)
+    # print('本网格sum',grid2d[i, j, 0])
+    # print('本网格点数',grid2d[i, j, 1])
+    
+    # print('sum: ', sum)
+    # print('~~~~~~~~~~~~~')
+
+print('插值均值：', interSum / interCount)
+
+print('体积：', sum)
+print('插值次数：', interCount)
+
+exit(0)

point_cloud_paint.py

@@ -0,0 +1,23 @@
+# ----------------------------------------------------------------------------
+# -                        Open3D: www.open3d.org                            -
+# ----------------------------------------------------------------------------
+# Copyright (c) 2018-2023 www.open3d.org
+# SPDX-License-Identifier: MIT
+# ----------------------------------------------------------------------------
+
+import open3d as o3d
+import numpy as np
+
+if __name__ == "__main__":
+    # pcd = o3d.io.read_point_cloud("./data/laser1.ply")
+    pcd = o3d.io.read_point_cloud("./data/PointClound_01.ply")
+    print('pcd data: ~~~~~~~~~~~~~~~~~~~')
+    print(pcd)
+    # print(np.asarray(pcd.points))
+    # print(np.asarray(pcd.normals))
+    # print(np.asarray(pcd.covariances))
+    print(np.asarray(pcd.colors))
+    print('end of pcd data ~~~~~~~~~~~~~~~~~~~')
+    o3d.visualization.draw([pcd])
+    # pcd.paint_uniform_color([1, 0.706, 0])
+    # o3d.visualization.draw([pcd])

readme.md

@@ -0,0 +1,112 @@
+# .ply文件内部啥样、数字单位是多少
+是二进制文件。坐标单位是米。有坐标、颜色、法线信息，无协方差信息。
+
+# google上有吗？
+
+## ~~https://jose-llorens-ripolles.medium.com/stockpile-volume-with-open3d-fa9d32099b6f~~
+只适用于土堆等（stockpile）堆在平地上的、从上方向地面看去表面一览无余的物体
+
+## https://stackoverflow.com/questions/76752991/open3d-how-to-calculate-the-volume-of-a-mesh-created-by-point-cloud
+对alpha shaped mesh求体积时报错：invalid tetra in TetraMesh，是因为mesh内部而非表面也有点。
+
+可通过计算每个点到表面的距离来去除内部和外部的点。
+
+用tri_mesh.is_watertight()判断是否封闭。
+
+~~用这个算法同样不能确保封闭。create_from_point_cloud_ball_pivoting~~
+
+~~也不用考虑Poisson surface reconstruction，这个算法可以保证mesh封闭，但open3d没有实现。~~
+
+## Struggling to create watertight meshes out of point cloud data using Open3D in Python
+三个表面重构方法都不能保证封闭。但泊松方法洞最少，alpha shapes方法洞最多。
+
+用PyVista配合pymeshfix倒是可以修复mesh的漏洞、计算体积。但如何从点云得到PyVista mesh？？？这又是一个未知领域。
+```
+import pymeshfix
+import numpy as np
+import pyvista as pv
+
+pv.set_plot_theme('document')
+
+array = np.genfromtxt('ct_prostate_contour_data.csv', delimiter=',')
+
+point_cloud = pv.PolyData(array)
+surf = point_cloud.reconstruct_surface(nbr_sz=20, sample_spacing=2)
+
+mf = pymeshfix.MeshFix(surf)
+mf.repair()
+repaired = mf.mesh
+
+pl = pv.Plotter()
+pl.add_mesh(point_cloud, color='k', point_size=10)
+pl.add_mesh(repaired)
+pl.add_title('Reconstructed Surface')
+pl.show()
+```
+
+## ~~https://stackoverflow.com/questions/61269980/open3d-crop-pointcloud-with-polygon-volume~~
+讨论如何裁剪点云。
+
+提问者的代码已过时，而且肯定有问题他才会提问。所以不管了。
+
+第一个回答者的代码？？？
+
+第二个回答者的代码？？？
+
+## https://stackoverflow.com/questions/70449880/how-to-crop-a-point-cloud-in-open3d-by-using-a-polygon-volume
+裁剪相关
+
+## https://github.com/isl-org/Open3D/pull/1218
+裁剪相关
+
+# 我的流程
+
+## open3d
+前提：明确点云遗漏的漏洞最大直径r（即直径多少的球刚好能掉入漏洞）不超过多少。显然r完美情况下是0.01
+
+点云降采样（颗粒度g，10厘米？）
+
+(丢弃噪点？remove_statistical_outlier， https://jose-llorens-ripolles.medium.com/stockpile-volume-with-open3d-fa9d32099b6f)
+
+表面重构，采用的参数依据r和g中较大者。
+
+（mesh简化？）
+
+(mesh求连同区域，舍弃零碎的？)
+
+求体积
+
+### 切割后补点
+提供原始点和切割多面体
+
+
+# 测试环境管理帐号
+* 88888888888
+* 2222222222a
+
+# 下载点云文件
+https://laser-data.oss-cn-shenzhen.aliyuncs.com/testdata/SG-t-MnQ0zoEQIBX/data/SG-t-MnQ0zoEQIBX_ply.zip
+
+# 查看方框的四个点
+控制台输入：
+```
+viewer.modules.clipping.getAllBoxes().filter(e=>!e.isNew).map(e=>e.box.geometry.vertices.map(p=>p.clone().applyMatrix4(e.matrixWorld)))
+```
+
+# 结果
+
+## https://uat-laser.4dkankan.com/uat/index.html?m=SG-t-MnQ0zoEQIBX#/
+* （按尺寸标注计算：0.470立方米）
+* stockpile voxel: 0.43 立方米
+
+## https://uat-laser.4dkankan.com/uat/index.html?m=SG-t-AsSRO3iA0XT#/
+* （按尺寸标注计算：0.689立方米）
+* stockpile voxel: 0.75立方米
+
+## https://uat-laser.4dkankan.com/uat/index.html?m=SG-t-cp3kSewxMKi#/
+* （按尺寸标注计算：0.689立方米）
+* stockpile voxel: 0.70立方米
+
+## https://uat-laser.4dkankan.com/uat/index.html?m=SG-t-kuuv0moEP0C#/
+*（按尺寸标注计算：0.689立方米）
+* stockpile voxel: 0.84立方米

watch-help.py

@@ -0,0 +1,4 @@
+import open3d as o3d
+import numpy as np
+# get_oriented_bounding_box
+help(o3d.geometry.AxisAlignedBoundingBox)