三维模型(.obj)体素化 matlab实现

通过计算几何，将输入obj模型根据设置的体素边长进行体素化

为什么我不是赵天宇呢

2434人浏览 · 2023-09-01 15:17:35

为什么我不是赵天宇呢 · 2023-09-01 15:17:35 发布

基于计算几何体素化obj三维模型：

所包含函数已经更新，抱歉看漏了

写的时候函数嵌套太严重了，如果还有缺失可以再联系我

邮箱1318823905@qq.com

原理：

函数read_obj2：网上一找一大堆（已经更新）
- 读取之后归零化，后续计算会带来便利（大概吧，我是这么觉得的）
- 这个函数只支持纯三角面的obj模型，建议打上断点排查下是否这里有问题（后续再更新）
函数voxelization：主要是计算三角面和体素相交，判断哪个位置的体素逻辑值设置1
- 为减少运算时间，循环中计算三角面片的三维外包框所占据的体素，再进行相交判断
- 具体计算原理后续更新（懒）
其余函数有注释，可当黑箱子
matlab版本为2020a 稍微高级点的版本应该都可以

输入内容&参数设置：

obj文件
体素边长（正方体）
体素边界预留个数（不需要可设置0）

需要那个模型可以联系我，随缘登录看到了回复（大概吧）

输出内容：

voxels 结构体
- logical ：体素逻辑矩阵
- centerpoints ：体素中心坐标（x,y,z）

结果演示：

原始模型
体素边长设置为0.02
体素边长设置为0.05

函数read_obj2

function [vertices, faces, normals] = read_obj2(filename)
    fid = fopen(filename, 'r');
    if fid == -1
        error(['Cannot open file ' filename]);
    end
    
    vertices = [];
    normals = [];
    faces = [];

    while ~feof(fid)
        line = fgetl(fid);
        if isempty(line) || line(1) == '#' % 忽略空行和注释
            continue;
        end
        
        tokens = strsplit(line);
        
        if strcmp(tokens{1}, 'v')
            % 读取顶点坐标
            x = str2double(tokens{2});
            y = str2double(tokens{3});
            z = str2double(tokens{4});
            vertices(end+1,:) = [x y z];
        elseif strcmp(tokens{1}, 'vn')
            % 读取法向量
            x = str2double(tokens{2});
            y = str2double(tokens{3});
            z = str2double(tokens{4});
            normals(end+1,:) = [x y z];
        elseif strcmp(tokens{1}, 'f')
            % 读取面
            v1 = str2double(strsplit(tokens{2}, '/'));
            v2 = str2double(strsplit(tokens{3}, '/'));
            v3 = str2double(strsplit(tokens{4}, '/'));
            faces(end+1,:) = [v1(1) v2(1) v3(1)];
        end
    end
    
    fclose(fid);
end

函数voxelization

function [voxelGrid, voxelCentersX, voxelCentersY, voxelCentersZ] = voxelization(vertices, faces, voxelSize, borderSize)
minVertex = min(vertices);
maxVertex = max(vertices);
NN=voxelSize;
% Extend the bounding box by the border size
extendedMinVertex = minVertex - borderSize * voxelSize;
extendedMaxVertex = maxVertex + borderSize * voxelSize;

dimensions = ceil((extendedMaxVertex - extendedMinVertex) / voxelSize);
voxelGrid = zeros(dimensions);

% Calculate voxel centers' coordinates
voxelCentersX = extendedMinVertex(1) + voxelSize/2 + (0:dimensions(1)-1) * voxelSize;
voxelCentersY = extendedMinVertex(2) + voxelSize/2 + (0:dimensions(2)-1) * voxelSize;
voxelCentersZ = extendedMinVertex(3) + voxelSize/2 + (0:dimensions(3)-1) * voxelSize;



for i=1:size(faces,1)  %对每一个面，都进行计算

    %找出该次循环需要计算的三角形面片的三个顶点坐标
    tri_point1=vertices(faces(i,1),:);
    tri_point2=vertices(faces(i,2),:);
    tri_point3=vertices(faces(i,3),:);
    tri_points=[tri_point1;tri_point2;tri_point3];
    %计算外包框最大最小值
    tri_outbox_min=min([tri_point1;tri_point2;tri_point3], [], 1);
    tri_outbox_max=max([tri_point1;tri_point2;tri_point3], [], 1);

    %该次循环中三角形面片外包框所占据的体素范围
    this_min=floor(tri_outbox_min/NN)+1+borderSize;
    this_max=floor(tri_outbox_max/NN)+1+borderSize;

    for ii=this_min(1):this_max(1)
        for jj=this_min(2):this_max(2)
            for kk=this_min(3):this_max(3)
                this_center_point=[voxelCentersX(ii) voxelCentersY(jj) voxelCentersZ(kk)];%当前判断的立方体中心坐标


                for ll=1:12  %对于每一个体素（12条边），每个边对该次三角形面片进行射线检测
                    [segementStart,segementEnd]=center_point_to_line_point(NN,this_center_point,ll);
                    if(checkIntersection(segementStart,segementEnd,tri_points))
                        voxelGrid(ii,jj,kk)=1;
                        break;
                    end
                end


                %三角形三个边对每个体素进行计算
                for mm=1:3
                    switch mm
                        case 1
                            segementStart=tri_point1;
                            segementEnd=tri_point2;
                        case 2
                            segementStart=tri_point2;
                            segementEnd=tri_point3;
                        case 3
                            segementStart=tri_point3;
                            segementEnd=tri_point1;
                    end

                    for nn=1:6
                        squarepoints=get_square_faces_points(NN,this_center_point,nn);

                        if(checkIntersection(segementStart,segementEnd,squarepoints(1:3,:)))
                            voxelGrid(ii,jj,kk)=1;
                            break;
                        end
                        if(checkIntersection(segementStart,segementEnd,squarepoints([1 3 4],:)))
                            voxelGrid(ii,jj,kk)=1;
                            break;
                        end
                    end
                end
            end
        end
    end
end




end

函数 center_point_to_line_point

function [line_start line_end] = center_point_to_line_point(NN,centerpoint,index)
%输入正方体中心坐标和边长和对应编号，获取对应的线段起始结束坐标
    squarepoints=center_point_to_squarepoints(NN,centerpoint);
    p1=squarepoints(1,:);
    p2=squarepoints(2,:);
    p3=squarepoints(3,:);
    p4=squarepoints(4,:);
    p5=squarepoints(5,:);
    p6=squarepoints(6,:);
    p7=squarepoints(7,:);
    p8=squarepoints(8,:);
    switch index
        case 1
            line_start=p1;
            line_end=p2;
        case 2
            line_start=p2;
            line_end=p3;
        case 3
            line_start=p3;
            line_end=p4;
        case 4
            line_start=p4;
            line_end=p1;
        case 5
            line_start=p5;
            line_end=p6;
        case 6
            line_start=p6;
            line_end=p7;
        case 7
            line_start=p7;
            line_end=p8;
        case 8
            line_start=p8;
            line_end=p5;
        case 9
            line_start=p6;
            line_end=p2;
        case 10
            line_start=p5;
            line_end=p1;
        case 11
            line_start=p8;
            line_end=p4;
        case 12
            line_start=p7;
            line_end=p3;
    end

end

函数checkIntersection

function isIntersect = checkIntersection(segmentStart, segmentEnd, triangleVertices)
    % 计算线段的方向向量
    segmentVector = segmentEnd - segmentStart;
    
    % 计算三角形的法向量
    triangleNormal = cross(triangleVertices(2,:) - triangleVertices(1,:), triangleVertices(3,:) - triangleVertices(1,:));
    
    % 判断线段与三角形是否平行
    if norm(segmentVector) < eps || norm(triangleNormal) < eps
        isIntersect = false; % 线段与三角形平行，不相交
        return;
    end
    
    % 计算线段参数方程的分母
    segmentDenominator = dot(triangleNormal, segmentVector);
    
    % 判断线段与三角形是否共面
    if abs(segmentDenominator) < eps
        isIntersect = false; % 线段与三角形共面，不相交
        return;
    end
    
    % 计算线段的起点到三角形顶点1的向量
    startToVertex1 = triangleVertices(1,:) - segmentStart;
    
    % 计算线段参数方程的参数t
    t = dot(triangleNormal, startToVertex1) / segmentDenominator;
    
    % 判断t是否在[0,1]的范围内
    if t < 0 || t > 1
        isIntersect = false; % t不在[0,1]范围内，线段不与三角形相交
        return;
    end
    
    % 计算交点
    intersectionPoint = segmentStart + t * segmentVector;
    
    % 判断交点是否在三角形内部
    isIntersect = isInsideTriangle(intersectionPoint, triangleVertices);
end

函数 get_square_faces_points

function squareVertices = get_square_faces_points(NN,this_center_point,i)

    squarepoints=center_point_to_squarepoints(NN,this_center_point);
    p1=squarepoints(1,:);
    p2=squarepoints(2,:);
    p3=squarepoints(3,:);
    p4=squarepoints(4,:);
    p5=squarepoints(5,:);
    p6=squarepoints(6,:);
    p7=squarepoints(7,:);
    p8=squarepoints(8,:);
    
    switch i
        case 1
            squareVertices=[p1;p2;p3;p4];
        case 2
            squareVertices=[p5;p6;p7;p8];
        case 3
            squareVertices=[p1;p2;p6;p5];
        case 4
            squareVertices=[p3;p4;p8;p7];
        case 5
            squareVertices=[p2;p3;p7;p6];
        case 6
            squareVertices=[p1;p4;p8;p5];
    end
end

函数 center_point_to_squarepoints

function square_points = center_point_to_squarepoints(NN,centerpoint)
%输入正方体中心坐标和边长,获取所有顶点坐标
    square_points=zeros(8,3);
    square_points(1,:)=[centerpoint(1)-NN/2 centerpoint(2)-NN/2 centerpoint(3)-NN/2];
    square_points(2,:)=[centerpoint(1)+NN/2 centerpoint(2)-NN/2 centerpoint(3)-NN/2];
    square_points(3,:)=[centerpoint(1)+NN/2 centerpoint(2)+NN/2 centerpoint(3)-NN/2];
    square_points(4,:)=[centerpoint(1)-NN/2 centerpoint(2)+NN/2 centerpoint(3)-NN/2];
    square_points(5,:)=[centerpoint(1)-NN/2 centerpoint(2)-NN/2 centerpoint(3)+NN/2];
    square_points(6,:)=[centerpoint(1)+NN/2 centerpoint(2)-NN/2 centerpoint(3)+NN/2];
    square_points(7,:)=[centerpoint(1)+NN/2 centerpoint(2)+NN/2 centerpoint(3)+NN/2];
    square_points(8,:)=[centerpoint(1)-NN/2 centerpoint(2)+NN/2 centerpoint(3)+NN/2];
end

函数 isInsideTriangle

% 辅助函数，判断点是否在三角形内部
function inside = isInsideTriangle(point, triangleVertices)
    % 计算三角形的法向量和边向量
    triangleNormal = cross(triangleVertices(2,:) - triangleVertices(1,:), triangleVertices(3,:) - triangleVertices(1,:));
    edge1 = triangleVertices(2,:) - triangleVertices(1,:);
    edge2 = triangleVertices(3,:) - triangleVertices(2,:);
    edge3 = triangleVertices(1,:) - triangleVertices(3,:);
    
    % 计算点到三个边的向量
    toEdge1 = point - triangleVertices(1,:);
    toEdge2 = point - triangleVertices(2,:);
    toEdge3 = point - triangleVertices(3,:);
    
    % 判断点是否在三个边的同侧
    cross1 = dot(triangleNormal, cross(edge1, toEdge1));
    cross2 = dot(triangleNormal, cross(edge2, toEdge2));
    cross3 = dot(triangleNormal, cross(edge3, toEdge3));
    
    % 判断点是否在三角形内部
    inside = cross1 >= 0 && cross2 >= 0 && cross3 >= 0;
end

主脚本：


#运行的脚本
clear
clc
NN = 0.05; % 正方体边长 
borderSize=5; % 边界体素预留个数
[vertices, faces, normals] = read_obj2('haha.obj'); %选取你需要读取的obj文件
min_vertex = min(vertices, [], 1); %找到模型xyz最小值
vertices = vertices - min_vertex;  %归0化，避免负数值影响后续计算
[voxels.logical, x, y, z] = voxelization(vertices, faces, NN,borderSize);%体素化
for i =1:size(x,2)   %计算体素中心坐标，为后续计算提供便利
    for j =1:size(y,2)
        for k =1:size(z,2)
            voxels.centerpoint_x(i,j,k)=x(i);
            voxels.centerpoint_y(i,j,k)=y(j);
            voxels.centerpoint_z(i,j,k)=z(k);
        end
    end
end
xx=x;
yy=y;
zz=z;
% [x,y,z] = meshgrid(x, y, z);  
voxel_num=[ size(voxels.centerpoint_x,1)  size(voxels.centerpoint_x,2)  size(voxels.centerpoint_x,3) ];
voxel_size = [NN,NN,NN]; 

%绘制
figure
patch('Vertices', vertices, 'Faces', faces, 'FaceColor', 'green', 'FaceAlpha', 0.7, 'EdgeColor', 'black');
axis equal;
view(-30, 30);

hold on
% 在三维空间中绘制最小包围盒
% patch('Vertices', bbox_vertices, 'Faces', [1 2 3 4; 5 6 7 8; 1 2 6 5; 2 3 7 6; 3 4 8 7; 4 1 5 8], 'FaceColor', 'none', 'EdgeColor', 'red', 'LineWidth', 2);
% 对于每个体素，如果其值为1，则在相应的位置生成一个小正方体

for i = 1:voxel_num(1)
    for j = 1:voxel_num(2)
        for k = 1:voxel_num(3)
%             if voxels.logical(i, j, k) == 1
            if voxels.logical(i, j, k) == 1
                x = xx(i);
                y = yy(j);
                z = zz(k);
                verticess = [x-NN/2, y-NN/2, z-NN/2;x+NN/2, y-NN/2, z-NN/2;x+NN/2, y+NN/2, z-NN/2;x-NN/2, y+NN/2, z-NN/2;x-NN/2, y-NN/2, z+NN/2;                    x+NN/2, y-NN/2, z+NN/2;                    x+NN/2, y+NN/2, z+NN/2;                    x-NN/2, y+NN/2, z+NN/2;                ];
                % 创建6个面的面索引
                facess = [ 1, 2, 3, 4; 2, 6, 7, 3;6, 5, 8, 7;  5, 1, 4, 8;  1, 2, 6, 5;  4, 3, 7, 8;  ];
                % 添加立方体到图形对象
                patch('Vertices', verticess, 'Faces', facess, 'FaceColor', 'none','EdgeColor', 'black' ,'LineWidth', 1);            
            end
        end
    end
end
axis equal;
view(-30, 30);

技术共进，成长同行——讯飞AI开发者社区

更多推荐

大数据领域数据架构的人工智能融合应用

随着企业数据规模以年均40%的速度爆炸式增长（Gartner, 2023），传统数据架构在处理多模态数据、支撑实时智能决策时面临效率瓶颈。如何通过AI实现数据治理的自动化与智能化机器学习如何优化数据存储与查询效率深度学习模型在实时数据流处理中的应用范式智能决策系统与数据中台的架构耦合机制背景部分定义核心概念与技术演进路径核心章节解析融合架构的技术组件与算法实现实战篇提供完整的端到端解决方案应用篇呈

讯飞AI开发者社区

JWT原理

例如：随着人工智能的不断发展，机器学习这门技术也越来越重要，很多人都开启了学习机器学习，本文就介绍了机器学习的基础内容。提示：以下是本篇文章正文内容，下面案例可供参考提示：这里对文章进行总结：例如：以上就是今天要讲的内容，本文仅仅简单介绍了pandas的使用，而pandas提供了大量能使我们快速便捷地处理数据的函数和方法。