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

为什么我不是赵天宇呢

已于 2023-11-17 00:32:10 修改

阅读量882

点赞数 2

文章标签： matlab 拓扑学几何学

于 2023-09-01 15:17:35 首次发布

本文链接：https://blog.csdn.net/weixin_52969081/article/details/132624197

版权

基于计算几何体素化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);