Voxelized obj three-dimensional model based on computational geometry:
Principle:
- Function read_obj2: I can find a lot of them online.
- Resetting to zero after reading will bring convenience to subsequent calculations (probably, this is what I think)
- Function voxelization: mainly calculates the intersection of the triangle surface and the voxel, and determines which position of the voxel logical value is set to 1
- In order to reduce the calculation time, the voxels occupied by the three-dimensional outer box of the triangular patch are calculated in the loop, and then the intersection judgment is performed.
- The specific calculation principle will be updated later (lazy)
Input content & parameter settings:
- obj file
- Voxel side length (cube)
- Number of voxel boundaries reserved (can be set to 0 if not required)
- If you need that model, you can contact me. I logged in casually and saw the reply (probably)
Output content:
- voxels structure
- logical: voxel logical matrix
- centerpoints: voxel center coordinates (x, y, z)
Result demonstration:
- original model
- Voxel side length is set to 0.02
- Voxel side length is set to 0.05
Function read_obj2
function [vertices, faces, normals] = read_obj2(filename)
fid = fopen(filename, 'r');
iffid==-1
error(['Cannot open file ' filename]);
end
vertices = [];
normals = [];
faces = [];
while ~feof(fid)
line = fgetl(fid);
if isempty(line) || line(1) == '#' % ignore empty lines and comments
continue;
end
tokens = strsplit(line);
if strcmp(tokens{<!-- -->1}, 'v')
% Read the vertex coordinates
x = str2double(tokens{<!-- -->2});
y = str2double(tokens{<!-- -->3});
z = str2double(tokens{<!-- -->4});
vertices(end + 1,:) = [x y z];
elseif strcmp(tokens{<!-- -->1}, 'vn')
% Read normal vector
x = str2double(tokens{<!-- -->2});
y = str2double(tokens{<!-- -->3});
z = str2double(tokens{<!-- -->4});
normals(end + 1,:) = [x y z];
elseif strcmp(tokens{<!-- -->1}, 'f')
% reading side
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
Function 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) % Calculate each face
%Find the three vertex coordinates of the triangle patch that need to be calculated for this loop
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];
% Calculate the maximum and minimum values of the outer box
tri_outbox_min=min([tri_point1;tri_point2;tri_point3], [], 1);
tri_outbox_max=max([tri_point1;tri_point2;tri_point3], [], 1);
%The voxel range occupied by the outer frame of the triangle patch in this cycle
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)];%Currently judged cube center coordinates
for ll=1:12 %For each voxel (12 edges), each edge performs ray detection on the sub-triangular patch
[segementStart,segementEnd]=center_point_to_line_point(NN,this_center_point,ll);
if(checkIntersection(segmentStart,segmentEnd,tri_points))
voxelGrid(ii,jj,kk)=1;
break;
end
end
%The three sides of the triangle are calculated for each voxel
for mm=1:3
switch mm
case 1
segmentStart=tri_point1;
segmentEnd=tri_point2;
case 2
segmentStart=tri_point2;
segmentEnd=tri_point3;
case 3
segmentStart=tri_point3;
segmentEnd=tri_point1;
end
fornn=1:6
squarepoints=get_square_faces_points(NN,this_center_point,nn);
if(checkIntersection(segmentStart,segmentEnd,squarepoints(1:3,:)))
voxelGrid(ii,jj,kk)=1;
break;
end
if(checkIntersection(segmentStart,segmentEnd,squarepoints([1 3 4],:)))
voxelGrid(ii,jj,kk)=1;
break;
end
end
end
end
end
end
end
end
Main script:
#Running script
clear
clc
NN = 0.05; % side length of cube
borderSize=5; % reserved number of border voxels
[vertices, faces, normals] = read_obj2('haha.obj'); %Select the obj file you need to read
min_vertex = min(vertices, [], 1); %Find the minimum xyz value of the model
vertices = vertices - min_vertex; % normalize to 0 to avoid negative values affecting subsequent calculations
[voxels.logical, x, y, z] = voxelization(vertices, faces, NN,borderSize);% voxelization
for i =1:size(x,2) % calculates the voxel center coordinates to facilitate subsequent calculations
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];
%draw
figure
patch('Vertices', vertices, 'Faces', faces, 'FaceColor', 'green', 'FaceAlpha', 0.7, 'EdgeColor', 'black') ;
axis equal;
view(-30, 30);
hold on
% Draw the minimum bounding box in three-dimensional space
% 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);
% For each voxel, if its value is 1, a small cube is generated at the corresponding position
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; ];
% Create a face index of 6 faces
faces = [ 1, 2, 3, 4; 2, 6, 7, 3; 6, 5, 8, 7; 5, 1, 4, 8; 1, 2, 6, 5; 4, 3, 7, 8 ; ];
% Add cube to graphics object
patch('Vertices', verticess, 'Faces', facess, 'FaceColor', 'none','EdgeColor', 'black' ,'LineWidth', 1) ;
end
end
end
end
axis equal;
view(-30, 30);
