Ken Xiao
tech

饶手测量

November 21, 2025

使用3D模型测量手掌面积

摘要

本文介绍了一种通过3D扫描模型投影与数学公式结合的方法，计算手掌的投影面积。通过加载3D模型，生成深度图，提取轮廓并计算面积，可以精确地测量手掌的平面投影面积。本文详细介绍了实现的关键步骤，并辅以数学公式和代码说明。

方法概述

核心步骤：

  1. 加载3D模型： 加载手掌的3D扫描模型（.obj 文件），获取顶点和三角面片数据。这里说一下使用iPhone 13 pro 的lidar雷达做的3d识别

  2. 生成深度图： 将3D模型投影到指定平面（如 z 平面），生成2D深度图。

  3. 提取轮廓： 基于深度图提取手掌的边界轮廓。

  4. 计算面积： 使用数学方法根据轮廓点计算手掌的投影面积。

数学公式

1. 模型投影公式

将3D网格投影到二维平面时，顶点的坐标投影公式如下：

投影点 (x,y,z)={(x,y)投影到 z 平面(x,z)投影到 y 平面(y,z)投影到 x 平面\text{投影点} \, (x, y, z) = \begin{cases} (x, y) & \text{投影到 $z$ 平面} \\ (x, z) & \text{投影到 $y$ 平面} \\ (y, z) & \text{投影到 $x$ 平面} \end{cases}

其中，(x,y,z)(x, y, z) 是3D顶点的坐标，(x′,y′)(x', y') 是投影后的2D坐标。

2. 轮廓面积计算公式

对于提取的轮廓点 {(x1,y1),(x2,y2),…,(xn,yn)}\{(x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\}，使用多边形面积公式计算面积：

A=12∣∑i=1n(xiyi+1−yixi+1)∣A = \frac{1}{2} \left| \sum_{i=1}^{n} (x_i y_{i+1} - y_i x_{i+1}) \right|

其中，xn+1=x1x_{n+1} = x_1 和 yn+1=y1y_{n+1} = y_1，即轮廓点首尾相连。

3. 深度图插值与平滑

为提高深度图的质量，使用高斯滤波平滑深度值：（由傅里叶变换推出，了解即可）

G(x,y)=12πσ2e−x2+y22σ2G(x, y) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2 + y^2}{2\sigma^2}}

where x is the distance from the origin in the horizontal axis, y is the distance from the origin in the vertical axis, and σ is the standard deviation of the Gaussian distribution（正态分布）. 这可以减少噪声和不规则性，使轮廓提取更加精确。

方法1 插入手机录屏 屏幕截图 2024-11-26 093545.png

屏幕截图 2024-11-26 093507.png

import numpy as np  
import open3d as o3d  
import matplotlib.pyplot as plt  
from scipy.spatial import ConvexHull  
from skimage.measure import find_contours  
from scipy.ndimage import gaussian_filter  # 添加这行  
from scipy.interpolate import griddata    # 添加这行  
  
def load_obj_file(file_path):  
    """  
    加载 .obj 文件并返回三角网格和顶点数据  
    """    mesh = o3d.io.read_triangle_mesh(file_path)  
    if not mesh.has_triangles():  
        raise ValueError("The provided .obj file does not contain triangle mesh data.")  
    return mesh  
  
  
def generate_depth_map(mesh, plane_axis='z', resolution=1024):  
    """  
    将 3D 模型投影到一个特定平面以生成深度图  
    """    # 将 open3d Mesh 转换为 numpy 格式  
    vertices = np.asarray(mesh.vertices)  
    face_indices = np.asarray(mesh.triangles)  
  
    # 打印模型顶点范围  
    print(f"模型顶点范围: min={vertices.min(axis=0)}, max={vertices.max(axis=0)}")  
  
    # 投影到给定平面  
    if plane_axis == 'z':  
        projection_axes = [0, 1]  # x, y  
        depth_axis = 2  # z  
    elif plane_axis == 'y':  
        projection_axes = [0, 2]  # x, z  
        depth_axis = 1  # y  
    elif plane_axis == 'x':  
        projection_axes = [1, 2]  # y, z  
        depth_axis = 0  # x  
    else:  
        raise ValueError("Invalid plane_axis. Must be one of ['x', 'y', 'z'].")  
  
    # 提取投影坐标  
    xy = vertices[:, projection_axes]  
    depth = vertices[:, depth_axis]  
  
    # 打印深度值统计信息  
    print(f"深度值统计: min={depth.min()}, max={depth.max()}, mean={depth.mean()}")  
  
    # 构建深度图网格  
    min_xy = xy.min(axis=0)  
    max_xy = xy.max(axis=0)  
    extent = max_xy - min_xy  
    pixel_size = extent / resolution  
    grid_shape = (resolution, resolution)  
  
    # 初始化网格  
    grid = np.full(grid_shape, np.nan, dtype=np.float32)  
  
    # 映射点到网格  
    for face in face_indices:  
        face_vertices = vertices[face]  
        face_xy = face_vertices[:, projection_axes]  
        face_depth = face_vertices[:, depth_axis]  
  
        # 计算面片的边界框  
        min_face_xy = face_xy.min(axis=0)  
        max_face_xy = face_xy.max(axis=0)  
  
        # 映射到网格坐标  
        grid_min = np.floor((min_face_xy - min_xy) / pixel_size).astype(int)  
        grid_max = np.ceil((max_face_xy - min_xy) / pixel_size).astype(int)  
  
        # 确保在有效范围内  
        grid_min = np.maximum(grid_min, 0)  
        grid_max = np.minimum(grid_max, resolution - 1)  
  
        # 填充该区域  
        for y in range(grid_min[1], grid_max[1] + 1):  
            for x in range(grid_min[0], grid_max[0] + 1):  
                if 0 <= x < resolution and 0 <= y < resolution:  
                    grid[y, x] = np.nanmax([grid[y, x], face_depth.mean()])  
  
    # 添加平滑处理  
    grid = gaussian_filter(grid, sigma=2)  
  
    # 填充空洞  
    xx, yy = np.meshgrid(np.arange(grid.shape[1]), np.arange(grid.shape[0]))  
    valid_mask = ~np.isnan(grid)  
    grid = griddata(  
        (xx[valid_mask], yy[valid_mask]),  
        grid[valid_mask],  
        (xx, yy),  
        method='nearest'  
    )  
  
    # 打印深度图统计信息  
    print(  
        f"深度图统计: min={np.nanmin(grid)}, max={np.nanmax(grid)}, mean={np.nanmean(grid)}, 非零值数量={np.count_nonzero(~np.isnan(grid))}")  
  
    return grid, extent  
  
def extract_contours_and_area(depth_map, threshold=None, pixel_size=(1.0, 1.0)):  
    """  
    提取轮廓并计算面积  
    """    if threshold is None:  
        # 使用百分位数而不是平均值  
        valid_depths = depth_map[~np.isnan(depth_map)]  
        threshold = np.percentile(valid_depths, 25)  # 使用第25百分位数作为阈值  
        print(f"使用的动态阈值: {threshold}")  
  
    # 创建二值图像  
    binary_image = np.zeros_like(depth_map)  
    binary_image[~np.isnan(depth_map) & (depth_map > threshold)] = 1  
  
    # 可视化深度图直方图  
    valid_depth_values = depth_map[~np.isnan(depth_map)]  
    plt.hist(valid_depth_values.flatten(), bins=50, color='blue', alpha=0.7)  
    plt.title("Depth Value Histogram")  
    plt.xlabel("Depth Value")  
    plt.ylabel("Frequency")  
    plt.show()  
  
    # 提取轮廓  
    contours = find_contours(binary_image, level=0.5)  
    if not contours:  
        raise ValueError("No contours found in the depth map.")  
  
    print(f"找到的轮廓数量: {len(contours)}")  
    for i, contour in enumerate(contours):  
        print(f"轮廓 {i} 点数量: {len(contour)}")  
  
    # 获取最大轮廓  
    contour = max(contours, key=len)  
    print(f"最大的轮廓点数: {len(contour)}")  
  
    # 如果轮廓点不足，直接抛出错误  
    if len(contour) < 3:  
        raise ValueError("轮廓点不足，无法计算凸包面积。")  
  
    # 计算面积  
    contour = np.array(contour)  
    x = contour[:, 1] * pixel_size[0]  
    y = contour[:, 0] * pixel_size[1]  
    area = 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1)))  
  
    return area, contour  
  
  
def analyze_model_dimensions(vertices):  
    """  
    分析模型在xyz三轴上的尺寸  
    """    # 获取xyz轴的最大最小值  
    x_min, y_min, z_min = vertices.min(axis=0)  
    x_max, y_max, z_max = vertices.max(axis=0)  
  
    # 计算各轴的范围  
    x_range = x_max - x_min  
    y_range = y_max - y_min  
    z_range = z_max - z_min  
  
    # 计算各轴的均值和中位数  
    x_mean, y_mean, z_mean = vertices.mean(axis=0)  
    x_median = np.median(vertices[:, 0])  
    y_median = np.median(vertices[:, 1])  
    z_median = np.median(vertices[:, 2])  
  
    # 打印详细信息  
    print("\n=== 模型尺寸分析 ===")  
    print(f"\nX轴 (宽度):")  
    print(f"  最小值: {x_min:.2f}mm")  
    print(f"  最大值: {x_max:.2f}mm")  
    print(f"  范围: {x_range:.2f}mm")  
    print(f"  均值: {x_mean:.2f}mm")  
    print(f"  中位数: {x_median:.2f}mm")  
  
    print(f"\nY轴 (高度):")  
    print(f"  最小值: {y_min:.2f}mm")  
    print(f"  最大值: {y_max:.2f}mm")  
    print(f"  范围: {y_range:.2f}mm")  
    print(f"  均值: {y_mean:.2f}mm")  
    print(f"  中位数: {y_median:.2f}mm")  
  
    print(f"\nZ轴 (深度):")  
    print(f"  最小值: {z_min:.2f}mm")  
    print(f"  最大值: {z_max:.2f}mm")  
    print(f"  范围: {z_range:.2f}mm")  
    print(f"  均值: {z_mean:.2f}mm")  
    print(f"  中位数: {z_median:.2f}mm")  
  
    print("\n模型总体尺寸:")  
    print(f"  宽 × 高 × 深: {x_range:.2f}mm × {y_range:.2f}mm × {z_range:.2f}mm")  
    print(f"  表面积估计: {x_range * z_range:.2f}平方毫米")  
  
    return {  
        'x': {'min': x_min, 'max': x_max, 'range': x_range, 'mean': x_mean, 'median': x_median},  
        'y': {'min': y_min, 'max': y_max, 'range': y_range, 'mean': y_mean, 'median': y_median},  
        'z': {'min': z_min, 'max': z_max, 'range': z_range, 'mean': z_mean, 'median': z_median}  
    }  
  
  
  
  
def main(obj_file_path):  
    # 加载 .obj 文件  
    mesh = load_obj_file(obj_file_path)  
  
    # 分析模型尺寸  
    vertices = np.asarray(mesh.vertices)  
    dimensions = analyze_model_dimensions(vertices)  
  
    # 生成深度图（投影到 Z 平面）  
    depth_map, extent = generate_depth_map(mesh, plane_axis='z', resolution=1024)  
  
    # 调试：可视化深度图  
    plt.imshow(depth_map, cmap='viridis')  
    plt.colorbar(label="Depth")  
    plt.title("Generated Depth Map")  
    plt.show()  
  
    # 输出像素大小  
    pixel_size = extent / depth_map.shape  
    print(f"像素大小: x={pixel_size[0]:.6f}, y={pixel_size[1]:.6f}")  
  
    # 提取轮廓并计算面积  
    try:  
        area, contour = extract_contours_and_area(depth_map, threshold=None, pixel_size=pixel_size)  
        print(f"投影平面的面积: {area:.6f}")  
  
        # 可视化深度图和轮廓  
        plt.figure(figsize=(10, 5))  
        plt.subplot(1, 2, 1)  
        plt.title("Depth Map")  
        plt.imshow(depth_map, cmap='viridis')  
        plt.colorbar(label="Depth")  
        plt.subplot(1, 2, 2)  
        plt.title("Contours")  
        plt.imshow(depth_map, cmap='gray')  
        plt.plot(contour[:, 1], contour[:, 0], color='red')  # 绘制轮廓  
        plt.show()  
  
    except ValueError as e:  
        print(f"错误: {e}")  
  
  
if __name__ == "__main__":  
    obj_file = "C://Users//xfc05//Downloads//Untitled_Scan_19_39_18//textured_output.obj"  # 替换为您的 .obj 文件路径  
    main(obj_file)

如何计算动态阈值（25百分位数）

动态阈值的计算基于 百分位数 (Percentile) 概念。以下将详细解释 25百分位数 的含义、计算方法以及代码实现。

百分位数的定义

百分位数 是一个统计学术语，用于描述数据分布的特定位置。第 p 百分位数 表示一组数据中有 p% 的数据小于或等于该值。

  • 25百分位数 是第一四分位数，表示数据中 最小的25%的值的上限。
  • 百分位数公式的数学表达： Qp=排序后的数据中第 p% 的位置对应的值Q_p = \text{排序后的数据中第 } p\% \text{ 的位置对应的值} 其中，( Q_p ) 是第 ( p ) 百分位数。

计算25百分位数的步骤

1. 数据准备

确定需要计算百分位数的一组数据（如深度图中的非 NaN 有效深度值）。

2. 数据排序

将数据按从小到大的顺序排列。

3. 确定位置索引

对于第 ( p ) 百分位数，位置索引 ( i ) 的计算公式为：

i=(N−1)⋅p100i = (N - 1) \cdot \frac{p}{100}

其中：

  • ( N )：数据点的个数。
  • ( p )：百分位数（例如，25百分位数对应 ( p = 25 )）。

4. 线性插值

如果 ( i ) 是整数，则直接取第 ( i ) 个数据点的值作为百分位数。如果 ( i ) 是小数，则通过线性插值计算百分位数：

Qp=(1−f)⋅x⌊i⌋+f⋅x⌊i⌋+1Q_p = (1 - f) \cdot x_{\lfloor i \rfloor} + f \cdot x_{\lfloor i \rfloor + 1}

其中：

  • x⌊i⌋x_{\lfloor i \rfloor}：位置索引的整数部分对应的数据点。
  • ff：位置索引的小数部分。

用代码实现动态阈值

以下代码展示如何计算25百分位数：

代码实现

import numpy as np

# 示例数据：深度图中的有效深度值（非 NaN）
depth_values = np.array([0.2, 0.5, 0.4, 0.1, 0.3, 0.6, 0.7, 0.8, 0.4, 0.2, 0.5])

# 使用 NumPy 计算 25 百分位数
percentile_25 = np.percentile(depth_values, 25)
print(f"25百分位数: {percentile_25}")

解释

  1. np.percentile 函数：

    • np.percentile(data, p) 计算数据集中第 ( p ) 百分位数的值。
    • 在代码中，p=25 表示计算第25百分位数。

  2. 计算过程：

    • 将 depth_values 排序为 [0.1, 0.2, 0.2, 0.3, 0.4, 0.4, 0.5, 0.5, 0.6, 0.7, 0.8]。
    • 根据公式： i=(N−1)⋅25100=(11−1)⋅0.25=2.5i = (N - 1) \cdot \frac{25}{100} = (11 - 1) \cdot 0.25 = 2.5
    • ( i = 2.5 ) 表示第25百分位数位于第2和第3个数据点之间（索引从0开始计数）。
    • 使用线性插值计算： Q25=(1−0.5)⋅0.2+0.5⋅0.3=0.25Q_{25} = (1 - 0.5) \cdot 0.2 + 0.5 \cdot 0.3 = 0.25

  3. 输出结果：

    25百分位数: 0.25

动态阈值的优点

  1. 自适应性：

    • 数据分布可能在不同深度图中有所不同，动态阈值根据当前数据的分布自动调整，避免使用固定阈值带来的误差。

  2. 鲁棒性：

    • 相比平均值或固定值，百分位数更能减少极值（如噪声或异常点）的影响。例如，25百分位数可以屏蔽掉较小的异常值。

  3. 适应性强：

    • 通过调整百分位数（如25%、50%、75%），可以灵活选择感兴趣的区域。例如：
      • 第25百分位数：选择较深的区域。
      • 第75百分位数：选择较浅的区域。

总结

  • 25百分位数 是一种动态阈值计算方法，能够有效过滤掉无关或背景区域。
  • 通过 np.percentile 函数可以快速计算百分位数，适用于各种数据分布。
  • 在深度图处理中，使用25百分位数可以提取较深的感兴趣区域（如手掌部分），为后续的轮廓提取和面积计算提供基础。

详细解释：轮廓提取与轮廓面积计算

什么是轮廓（Contour）？

1. 基本概念

  • 轮廓 是一个闭合或开放的曲线，表示图像中强度（像素值）相同的像素点的边界。
  • 在二值图像中（像素值仅为0或1），轮廓是值为1的像素点构成的边缘。

2. 轮廓的意义

  • 轮廓是物体的边界，可以表示感兴趣区域（如手掌）的形状。
  • 提取轮廓后，可以进一步计算面积、周长、形状特征等。

5. 确保轮廓点数足够

if len(contour) < 3:
    raise ValueError("轮廓点不足，无法计算凸包面积。")
  • 为什么至少需要3个点？
    • 面积计算需要闭合的多边形，而多边形至少需要3个点。
    • 如果轮廓点不足，无法形成一个有效的多边形，直接抛出错误。

总结

1. 轮廓的作用

  • 轮廓是描述二值图像中感兴趣区域形状的基础。
  • 提取轮廓后，可以进行面积计算、周长计算、形状分析等操作。

2. 面积公式的优点

  • 面积公式适用于任意闭合多边形。
  • 通过 pixel_size 转换，可以将像素面积精确转换为实际单位面积。

3. 可视化轮廓

  • 可以使用 matplotlib 可视化轮廓，验证提取结果： 
    plt.imshow(binary_image, cmap='gray')
plt.plot(contour[:, 1], contour[:, 0], color='red')
plt.title("Contour Visualization")
plt.show()

通过这段代码，我们可以从二值图像中提取手掌区域的轮廓，并精确计算其面积，为后续分析提供基础。

Comments