clear
clc
% you need to have this script, the stlread function (by Eric Johnson
% https://www.mathworks.com/matlabcentral/fileexchange/22409-stl-file-reader),
% and the file to be analyzed all in the same working directory.
% The output rotations should be applied one at a time, in the order of
% Z then X.
% Choose your file.
filename='your_model_here.stl';
% Acceptable overhang angle (in degrees). Must be between 0-90 with 45
% generally considered the standard.
acceptable_angle=45;
% Choose how many iterations you want to test. This will be automatically
% rounded if it does not produce equaly sized steps. 4050 will produce
% orientations in four degree increments; 16200 for two degrees; 64800
% for one degree.
iterations=64800;
% Choose the radius (in degrees) used to exclude nearby local minimums when
% culling the results. This prevents a noisy plot from producing many
% local minimums all within a small range to the exclusion of candiates
% in other orientations.
search_radius=30;
% Choose maximum number of local minimums to consider for refinement after
% culling is completed.
minimum_limit=10;
% Choose refinement iterations. Set to zero to not run refinement. This
% should be at least 50 to be worth the time. Like with iterations this
% will be rounded if necessary to produce even spacing.
refinements=225;
% Show graphs? 1 or 0.
show_graphs=1;
% % % From here on is the functional code. DO NOT CHANGE.
% Import the file
[F,V,N]=stlread(filename);
F=F';
V=V';
N=N';
[whatever,faces]=size(F);
% define steps
Rx_steps=round(sqrt(iterations/2));
Rz_steps=Rx_steps*2;
[Rx,Rz]=meshgrid(0:180/Rx_steps:180, 0:360/Rz_steps:360-360/Rz_steps);
% Define storage matrix
N_temp=N;
overhang_factor=zeros(Rx_steps+1,Rz_steps)';
% Find each face area
vector1=zeros(3,faces);
vector2=zeros(3,faces);
for x=1:faces
vector1(:,x)=(V(:,x*3)-V(:,x*3-1))';
vector2(:,x)=(V(:,x*3)-V(:,x*3-2))';
end
temp=cross(vector1,vector2);
Area=(sqrt(temp(1,:).^2+temp(2,:).^2+temp(3,:).^2)/2);
% Set plane vectors for finding angles
plane=zeros(size(N));
plane(3,:)=-1;
%do top orientation, all z rotations will be the same
N_temp=N;
% find angles, setting anything above the threshold to 90 so they
% don't count towards the area calculation
plate_angle=acosd(dot(N_temp,plane));
angle_modified=plate_angle;
angle_modified(angle_modified(:,:)>acceptable_angle)=90;
projected_factor=sum(abs(cosd(angle_modified).*Area));
overhang_factor(:,1)=projected_factor;
%do bottom orientation, all z rotations will be the same
% rotate to bottom
Ty=[cosd(180),0,sind(180);0,1,0;-sind(180),0,cosd(180)];
N_temp=Ty*N;
% angles
plate_angle=acosd(dot(N_temp,plane));
angle_modified=plate_angle;
angle_modified(angle_modified(:,:)>acceptable_angle)=90;
projected_factor=sum(abs(cosd(angle_modified).*Area));
overhang_factor(:,Rx_steps+1)=projected_factor;
%do the rest of the orientations
for X=2:Rx_steps
% show the progam is doing work
clc
Initial_Iteration_Progress=X/Rx_steps*100
for Z=1:Rz_steps
% do the rotations in sequence, z then x
Tz=[cosd(Rz(Z,X)),-sind(Rz(Z,X)),0; sind(Rz(Z,X)),cosd(Rz(Z,X)),0; 0,0,1];
Tx=[1,0,0; 0, cosd(Rx(Z,X)),-sind(Rx(Z,X)); 0,sind(Rx(Z,X)),cosd(Rx(Z,X))];
N_temp=Tz*N;
N_temp=Tx*N_temp;
% angles, again
plate_angle=acosd(dot(N_temp,plane));
angle_modified=plate_angle;
angle_modified(angle_modified(:,:)>acceptable_angle)=90;
projected_factor=sum(abs(cosd(angle_modified).*Area));
overhang_factor(Z,X)=projected_factor;
end
end
% % % Reduce Data
% vectorize
Rx_vector=Rx(:);
Rz_vector=Rz(:);
overhang_factor_vector=overhang_factor(:);
% sort
[overhang_factor_vector,order]=sort(overhang_factor_vector);
Rx_vector=Rx_vector(order);
Rz_vector=Rz_vector(order);
% initialize 'while' count data
[reduced_size,whatever]=size(overhang_factor_vector);
reduce_index=1;
% loop through the whole thing, need while loop because data it is
% dynamically resized
while reduce_index<=reduced_size
deleteme=zeros(size(overhang_factor_vector));
% set all points within range of current point to be deleted
deleteme(Rx_vector<=Rx_vector(reduce_index)+search_radius ...
& Rx_vector>=Rx_vector(reduce_index)-search_radius ...
& Rz_vector<=Rz_vector(reduce_index)+search_radius ...
& Rz_vector>=Rz_vector(reduce_index)-search_radius)...
=1;
deleteme(reduce_index)=0; %don't delete yourself, dummy
% delete candidates
overhang_factor_vector(deleteme==1)=[];
Rx_vector(deleteme==1)=[];
Rz_vector(deleteme==1)=[];
% move on to next step
reduce_index=reduce_index+1;
[reduced_size,whatever]=size(overhang_factor_vector);
% show the program is doing work
end
% keep only the best candidates
if reduced_size>minimum_limit
overhang_factor_vector=overhang_factor_vector(1:minimum_limit);
Rx_vector=Rx_vector(1:minimum_limit);
Rz_vector=Rz_vector(1:minimum_limit);
end
if refinements~=0
% find the array of sub-angles to test at each location
sides=round((sqrt(refinements)-1)/2);
[basic_square_x,basic_square_z]=meshgrid(-sides:sides,-sides:sides);
basic_square_x=basic_square_x*(180/Rx_steps/(sides+1));
basic_square_z=basic_square_z*(180/Rx_steps/(sides+1));
% do the refining
[points_to_refine,whatever]=size(overhang_factor_vector);
for refine=1:points_to_refine
refine_square_x=basic_square_x+Rx_vector(refine);
refine_square_z=basic_square_z+Rz_vector(refine);
refined_overhang_factor=zeros(size(refine_square_x));
for X=1:sides*2+1
for Z=1:sides*2+1
% do the rotations in sequence, z then x
Tz=[cosd(refine_square_z(Z,X)),-sind(refine_square_z(Z,X)),0; sind(refine_square_z(Z,X)),cosd(refine_square_z(Z,X)),0; 0,0,1];
Tx=[1,0,0; 0, cosd(refine_square_x(Z,X)),-sind(refine_square_x(Z,X)); 0,sind(refine_square_x(Z,X)),cosd(refine_square_x(Z,X))];
N_temp=Tz*N;
N_temp=Tx*N_temp;
% angles, again
plate_angle=acosd(dot(N_temp,plane));
angle_modified=plate_angle;
angle_modified(angle_modified(:,:)>acceptable_angle)=90;
projected_factor=sum(abs(cosd(angle_modified).*Area));
refined_overhang_factor(Z,X)=projected_factor;
end
end
% vectorize
refine_square_x=refine_square_x(:);
refine_square_z=refine_square_z(:);
refined_overhang_factor=refined_overhang_factor(:);
% sort
[refined_overhang_factor,order]=sort(refined_overhang_factor);
refine_square_x=refine_square_x(order);
refine_square_z=refine_square_z(order);
% save
Rx_vector(refine)=refine_square_x(1);
Rz_vector(refine)=refine_square_z(1);
overhang_factor_vector(refine)=refined_overhang_factor(1);
clc
refining_progress=refine/points_to_refine*100
end
end
% Give output of candidates
clc
ordered_canditates_Rx_Rz=[Rx_vector,Rz_vector]
if show_graphs==1
% Plot the Candidates
figure('name','Results and Candidates','units','normalized','outerposition',[0 0.05 0.6 0.9])
hold on
surf(Rx,Rz,overhang_factor,'linestyle','none');
title(['Projected Surface Area for Faces Below ',num2str(acceptable_angle),' Degrees'])
xlabel('X Rotation (Degrees)')
ylabel('Z Rotation (Degrees)')
axis equal
axis([0 180 0 360])
view(0,-90)
h1=scatter3(Rx_vector,Rz_vector,overhang_factor_vector,'black','filled');
h2=scatter3(Rx_vector(1),Rz_vector(1),overhang_factor_vector(1),'red','filled');
legend([h2,h1],'Best Candidate','Other Minimum Candidates','Location','eastoutside')
% Show the best orientations
orientations_to_plot=6;
if minimum_limit<6
orientations_to_plot=minimum_limit;
end
figure('name',['Top ',num2str(orientations_to_plot),' Candidates'],'units','normalized','outerposition',[0.05 0.05 0.9 0.9])
for x=1:orientations_to_plot
Tz=[cosd(Rz_vector(x)),-sind(Rz_vector(x)),0; sind(Rz_vector(x)),cosd(Rz_vector(x)),0; 0,0,1];
Tx=[1,0,0; 0, cosd(Rx_vector(x)),-sind(Rx_vector(x)); 0,sind(Rx_vector(x)),cosd(Rx_vector(x))];
V_temp=Tz*V;
V_temp=Tx*V_temp;
subplot(2,3,x)
patch('Faces',F','Vertices',V_temp','FaceColor', [0.8 0.8 1.0], ...
'EdgeColor', 'none', ...
'FaceLighting', 'gouraud', ...
'AmbientStrength', 0.15);
camlight('headlight');
material('dull');
axis('image');
title(['Overhang Area (',num2str(round(overhang_factor_vector(x))),') Orientation (Z ',num2str(round(Rz_vector(x))),' Degrees, X ',num2str(round(Rx_vector(x))),' Degrees)'])
view([45 45]);
end
end