-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmake_STL_of_Array.m
199 lines (171 loc) · 9.35 KB
/
make_STL_of_Array.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
function [Vertices, Triangle, Quads] = make_STL_of_Array(FileName,Data,scaleX,scaleY,scaleZ)
% make_STL_of_Array Convert a voxelised object contained within a 3D logical array into an STL surface mesh
%==========================================================================
% AUTHOR Amir-Hosein Safari
% CONTACT [email protected]
% INSTITUTION Max-Planck institute for informatic
% DATE 25th Aug 2018
%
% EXAMPLE make_STL_of_Array(FileName,Data,scaleX,scaleY,scaleZ)
% ..or.. [Vertices, Triangle, Quads] = make_STL_of_Array(FileName,Data,scaleX,scaleY,scaleZ)
%
% INPUTS FileName - string - Filename of the STL file.
%
% Data - 3D logical array - Voxelised data
% 1 => Inside the object
% 0 => Outside the object
% (FOR PRINTING WITH TWO MATERIALS YOU SHOULD INVERSE
% Data (Data = ~Data) AND RUN CODE AGAIN TO HAVE THE SECOND TYPE)
%
% scaleX - A number which means the X size of every
% voxel in mm
% scaleY - A number which means the Y size of every
% voxel in mm
% scaleZ - A number which means the Z size of every
% voxel in mm
%
%
% OUTPUTS vertices - Nx3 array - A list of the x,y,z coordinates of
% each vertex in the mesh.
%
% Triangle - Nx3 array - A list of the vertices used in
% each Triangle of the mesh, identified using the row
% number in the array vertices.
%
% Quads - Nx4 array - A list of the vertices used in
% each Quads of the mesh, identified using the row
% number in the array vertices.
%==========================================================================
%======================================================
% ALLOCATING THE MEMORY FOR MATRIXES
%======================================================
Vertices = zeros(8*size(Data,1)* size(Data,2)* size(Data,3), 3);
Quads = zeros(6*size(Data,1)* size(Data,2)* size(Data,3), 4);
Triangle = zeros(12*size(Data,1)* size(Data,2)* size(Data,3), 3);
Vertices_Pointer = 1;
Quads_Pointer = 1;
Triangle_Pointer = 1;
%======================================================
% LOOPING THROUGH EVERY ELEMENTS OF 3D LOGICAL ARRAY(Data)
%======================================================
for i3 = 1 : size(Data,3)
for i2 = 1 : size(Data,2)
for i1 = 1 : size(Data,1)
% i1 and i2 will change inside a for loop so we use i11 and
% i22 to save the value before changing and we will use them
i11 = i1;
i22 = i2;
% check that Data(i1,i2,i3) is not surronded by 1s!
if i1 ~= 1 && i2 ~= 1 && i3 ~= 1 && i1 ~= size(Data,1) && i2 ~= size(Data,2) && i3 ~= size(Data,3) && Data(i1 - 1,i2,i3) == 2 && Data(i1,i2 - 1,i3) == 2 && Data(i1,i2,i3 - 1) == 2 && (Data(i1 + 1,i2,i3) == 1 || Data(i1 + 1,i2,i3) == 2 ) && (Data(i1,i2 + 1,i3) == 1 || Data(i1,i2 + 1,i3) == 2 ) && (Data(i1,i2,i3 + 1) == 1 || Data(i1,i2,i3 + 1) == 2 ) && Data(i1,i2,i3) ~= 0
Data(i1,i2,i3) = 2;
continue;
end
if Data(i1,i2,i3) == 1
%======================================================
% LOOKING FOR THE BIGGEST RECTANGLE (BY ITS AREA) FROM THIS PONT
% IN 2D!
% FOR THIS WE LOOPING THROUGH X AND Y TILL FINDING A ZERO IN THAT DIRECTION!
% WE SAVE THE INDEX OF X AND Y OF THE BIGGEST AREA IN
% maxX, maxY
%======================================================
maxX = 0;
maxY = 0;
loop2Len = size(Data,2);
for i1 = i1 : size(Data,1)
i2 = i22;
% IF WE REACH A VOXEL FROM OUTSIDE WE BREAK FROM THE
% LOOP
if Data(i1,i2,i3) ~= 1
i1 = i1 - 1;
break;
end
for i2 = i2 : loop2Len
% IF WE REACH A VOXEL FROM OUTSIDE WE BREAK FROM THE
% LOOP
if Data(i1,i2,i3) ~= 1
loop2Len = i2 -1;
i2 = i2 - 1;
break;
else
if maxX*maxY < (i1 - i11 + 1) * (i2 - i22 + 1)
maxX = i1 - i11;
maxY = i2 - i22;
end
end
end
end
% WE CHANGE THE VALUE OF THE BIGGEST AREA TO 2, IT
% MEANS THAT WE'LL NOT CONSIDER THIS AREA AGAIN
Data(i11:i11 + maxX, i22:i22 + maxY, i3) = 2;
% CREATE THE VERTICES AND ADD THEM TO Vertices MATRIX
% AND INCREAS THE Vertices_Pointer
V1 = [ (i11 - 1) * scaleX , (i22 - 1) * scaleY, (i3 - 1) * scaleZ];
V2 = [ (i11 + maxX) * scaleX , (i22 - 1) * scaleY, (i3 - 1) * scaleZ];
V3 = [ (i11 - 1) * scaleX , (i22 + maxY) * scaleY, (i3 - 1) * scaleZ];
V4 = [ (i11 + maxX) * scaleX , (i22 + maxY) * scaleY, (i3 - 1) * scaleZ];
V5 = [ (i11 - 1) * scaleX , (i22 - 1) * scaleY, (i3) * scaleZ];
V6 = [ (i11 + maxX) * scaleX , (i22 - 1) * scaleY, (i3) * scaleZ];
V7 = [ (i11 - 1) * scaleX , (i22 + maxY) * scaleY, (i3) * scaleZ];
V8 = [ (i11 + maxX) * scaleX , (i22 + maxY) * scaleY, (i3) * scaleZ];
if Vertices_Pointer == 0
Vertices = cat(1,Vertices,V1);
end
Vertices(Vertices_Pointer,:) = V1;
V_pointer1 = Vertices_Pointer;
Vertices_Pointer = Vertices_Pointer + 1;
Vertices(Vertices_Pointer,:) = V2;
V_pointer2 = Vertices_Pointer;
Vertices_Pointer = Vertices_Pointer + 1;
Vertices(Vertices_Pointer,:) = V3;
V_pointer3 = Vertices_Pointer;
Vertices_Pointer = Vertices_Pointer + 1;
Vertices(Vertices_Pointer,:) = V4;
V_pointer4 = Vertices_Pointer;
Vertices_Pointer = Vertices_Pointer + 1;
Vertices(Vertices_Pointer,:) = V5;
V_pointer5 = Vertices_Pointer;
Vertices_Pointer = Vertices_Pointer + 1;
Vertices(Vertices_Pointer,:) = V6;
V_pointer6 = Vertices_Pointer;
Vertices_Pointer = Vertices_Pointer + 1;
Vertices(Vertices_Pointer,:) = V7;
V_pointer7 = Vertices_Pointer;
Vertices_Pointer = Vertices_Pointer + 1;
Vertices(Vertices_Pointer,:) = V8;
V_pointer8 = Vertices_Pointer;
Vertices_Pointer = Vertices_Pointer + 1;
% CREATE THE QUADS AND ADD THEM TO Quads MATRIX
% AND INCREAS THE Quads_Pointer
Quads(Quads_Pointer,:) = [V_pointer1, V_pointer2, V_pointer3, V_pointer4];
Quads_Pointer = Quads_Pointer + 1;
Quads(Quads_Pointer,:) = [V_pointer5, V_pointer6, V_pointer7, V_pointer8];
Quads_Pointer = Quads_Pointer + 1;
Quads(Quads_Pointer,:) = [V_pointer5, V_pointer6, V_pointer1, V_pointer2];
Quads_Pointer = Quads_Pointer + 1;
Quads(Quads_Pointer,:) = [V_pointer3, V_pointer4, V_pointer7, V_pointer8];
Quads_Pointer = Quads_Pointer + 1;
Quads(Quads_Pointer,:) = [V_pointer5, V_pointer1, V_pointer7, V_pointer3];
Quads_Pointer = Quads_Pointer + 1;
Quads(Quads_Pointer,:) = [V_pointer2, V_pointer6, V_pointer4, V_pointer8];
Quads_Pointer = Quads_Pointer + 1;
end
end
end
end
% CUT THE UNUSED ELEMENTS OF MATRIXES
Vertices = Vertices(1:Vertices_Pointer - 1,:);
Quads = Quads(1:Quads_Pointer - 1,:);
% REMOVE REPETITIOUS QUADS
Quads = unique(Quads,'rows','stable');
% CREATE THE TRIAGNLES FROM QUADS AND ADD THEM TO Triangle MATRIX
% AND INCREAS THE Triangle_Pointer
for i = 1 : size(Quads,1)
Triangle(Triangle_Pointer,:) = [Quads(i,1) Quads(i,2) Quads(i,3)];
Triangle_Pointer = Triangle_Pointer + 1;
Triangle(Triangle_Pointer,:) = [Quads(i,4) Quads(i,2) Quads(i,3)];
Triangle_Pointer = Triangle_Pointer + 1;
end
% CUT THE UNUSED ELEMENTS OF MATRIXES
Triangle = Triangle(1:Triangle_Pointer - 1,:);
% WRITE THE STL FILE
stlwrite(FileName, Triangle, Vertices);