%{

Impossible Shape Generator

based on "Anomalous Mirror Symmetry Generated by Optical Illusion"
by Prof. Kokichi Sugihara
https://www.mdpi.com/2073-8994/8/4/21/htm

code by Dr. Ben Finio, Science Buddies (www.sciencebuddies.org)
thanks to Dr. Mark Veillette for his contributions

This code takes a two-dimensional shape defined in the XY plane and
generates an "impossible" three-dimensional version of the shape, e.g.
an arrow that always appears to point to the right even when rotated
180 degrees, like that shown in Prof. Sugihara's paper. 
It then outputs STL files for 3D printing, which define
the surface of a 3D shape using triangles.

To edit the original 2D shape, change the y values of the
points in the a1y and a2y vectors in Section 1 of the code. Alternatively,
uncomment the section of code to load the points from a CSV file. 

The code is broken up into sections to make it easier to understand the
output. In the Matlab editor, use the "Run and Advance" button in the
SECTION part of the top menu to run the code one section at a time. 

%}

% clear variables, figures, and output from previous program
close all
clear all
clc
disp('Use the "Run and Advance" button to run the code one section at a time.')
disp('Read the comments in Section 1 to learn how to define your own 2D shape.')
disp('Then, run that section to plot the shape.')

%% SECTION 1: Setup and 2D shape

% define variables for generating the 3D shape (you can change these)
height = 0.5;        % height of the 3D shape
viewing_angle = 45;  % desired viewing angle for the 3D shape in degrees

% define some appearance variables for the plots (you can change these)
linewidth = 2;       % line width for plots
linecolor = 'blue';
facecolor = 'cyan';  % color for plots of 3D objects

% INPUT: define your 2D shape using parametric curves a1 and a2 in the x-y
% plane. Both curves must meet at their endpoints on the x-axis. Change the
% y values to generate different shapes. 

x = [-1:0.1:1]';  % x coordinates as column vector. Make sure your x values are evenly spaced and symmetric
                  % about the y axis.
a1x = x;          % curve a1 x values 
a1y = [0	      % curve a1 y values, change these to change the top half of the 2D shape
       0.5	      % you can also enter functions of t (lines, parabolas, sine/cosine etc.) 
       0.5        % to generate shapes. Functions must be monotone in x (only one y-value
       0.5	      % per x-value, i.e. no vertical lines).
       0.5	 
       0.5  
       0.5  
       0.5  
       0.5	 
       0.5  
       1  
       0.9  
       0.8  
       0.7	
       0.6	 
       0.5   
       0.4	
       0.3	 
       0.2  
       0.1 
       0];
a2x = x;        % curve a2 x values
a2y = [0	    % curve a2 y values, change these to change the bottom half of the 2D shape
    -0.5 
    -0.5 
    -0.5 
    -0.5	
    -0.5 
    -0.5 
    -0.5 
    -0.5	
    -0.5 
    -1	
    -0.9 
    -0.8 
    -0.7	
    -0.6 
    -0.5 
    -0.4 
    -0.3	
    -0.2 
    -0.1 
    0];

% Load data from a CSV file instead. The CSV file should consist of three
% columns, each with a header row. One column for the x coordinate, one
% column for the y coordinate of a1, and one column for the y coordinate of
% a2. Save your CSV in the same directory as this m-file, then uncomment the
% following lines and change the filename to match your file. This will
% overwrite the values for t, a1x, a1y, a2x, and a2y entered above.

% points = readmatrix('car.csv') ;       % load data from the CSV file
% x = points(:,1);                            % get t from first column
% a1x = x;                                    % curve a1 x values 
% a1y = points(:,2);                          % curve a1 y values
% a2x = x;                                    % curve a2 x values
% a2y = points(:,3);                          % curve a2 y values

% plot the 2D shape

figure                                        % open a new figure
subplot(2,3,1)                                % plot on the first graph of a 2x3 grid
plot(a1x,a1y,linecolor,'LineWidth',linewidth) % line plot of curve a1
hold on                                       % turn on hold so we don't overwrite the first plot
plot(a2x,a2y,linecolor,'LineWidth',linewidth) % line plot of curve a2
grid on                                       % turn the grid on
axis equal                                    % set the x and y axes to equal scaling
axis tight                                    % set axis limits to the limits of the shape
xlabel('x')                                   % label the x axis
ylabel('y')                                   % label the y axis
title('1: Original 2D Shape')                 % give the plot a title

disp(' ')
disp('Run the next section to plot the 3D curve.')
disp('Rotate the graph to see the illusion. The outline of the 2D')
disp('shape is only visible from a certain angle.')

%% SECTION 2: 3D Curve

% calculate the 3D parametric curves c1(t) and c2(t) using equations 5 and
% 6 from Prof. Sugihara's paper https://www.mdpi.com/2073-8994/8/4/21/htm

alpha = tan(deg2rad(viewing_angle)); % parameter used in the equations for c1z and c2z

c1x = x;                          % x component of c1
c1y = (a1y-flip(a2y))/2;          % y component of c1
c1z = alpha*(a1y+flip(a2y))/2;    % z component of c1

c2x = x;                          % x component of c2
c2y = (a2y-flip(a1y))/2;          % y component of c2
c2z = alpha*(a2y+flip(a1y))/2;    % z component of c2

% plot the 3D curve. Use the 'rotate 3D' feature in the figure to rotate
% the graph and see the illusion. You should only be able to see the
% equivalent of the 2D shape when viewed from the correct angle.

subplot(2,3,2)                                      % plot on second graph of 2x3 grid
plot3(c1x,c1y,c1z,linecolor,'LineWidth',linewidth)  % plot curve c1
hold on                                             % turn hold on so we don't overwrite the first plot
plot3(c2x,c2y,c2z,linecolor,'LineWidth',linewidth)  % plot curve of c2
label_plot('2: 3D Curve (rotate to see illusion)',viewing_angle) % label the plot and set 3D viewing properties
% see label_plot function definition at the end of this file

disp(' ')
disp('Run the next section to plot an offset of the 3D curve.')
disp('This defines the top and bottom edges of the 3D shape.')
                                       

%% SECTION 3: Offset 3D Curve

% offset the curves c1 and c2 in the z direction to create the top edge of
% the 3D shape using curves c3 and c4

c3z = c1z + height;     % increase z values of c1 curve by 'height' and store in new variable
c4z = c2z + height;     % increase z values of c1 curve by 'height' and store in new variable

% plot the top and bottom edges of the 3D shape

subplot(2,3,3)                                      % plot on third graph of 2x3 grid
plot3(c1x,c1y,c1z,linecolor,'LineWidth',linewidth)  % plot curve c1
hold on                                             % turn hold on so we don't overwrite the first plot
plot3(c2x,c2y,c2z,linecolor,'LineWidth',linewidth)  % plot curve of c2
plot3(c1x,c1y,c3z,linecolor,'LineWidth',linewidth)  % plot curve of c3
plot3(c2x,c2y,c4z,linecolor,'LineWidth',linewidth)  % plot curve of c4
label_plot('3: Top and Bottom Edges of 3D Shape',viewing_angle) % label the plot and set 3D viewing properties

disp(' ')
disp('Run the next section to plot all the points on the edges of the 3D shape.')
disp('We need to connect these points using triangles to make an STL file.')
disp('There are two different ways to do this.')

%% SECTION 4: Plot 3D Points

% store each of the four curves in a matrix with 3 columns for the x, y,
% and z coordinates. note that C3 is directly above C1 so has the same x
% and y coordinates, the same for C4 and C2. We will need the data in this
% format later when generating the triangularization for the STL files. 

C1 = [c1x,c1y,c1z];     % store the points for c1 in a single matrix
C2 = [c2x,c2y,c2z];     % store the points for c2 in a single matrix
C3 = [c1x,c1y,c3z];     % store the points for c3 in a single matrix
C4 = [c2x,c2y,c4z];     % store the points for c4 in a single matrix

% plot all of the points (vertices of triangles on the surface of the 3D
% shape)
subplot(2,3,4)          % plot on fourth graph of 2x3 grid
plot3(C1(:,1),C1(:,2),C1(:,3),'b.'); % plot the points stored in C1
hold on
plot3(C2(:,1),C2(:,2),C2(:,3),'b.');
plot3(C3(:,1),C3(:,2),C3(:,3),'b.');
plot3(C4(:,1),C4(:,2),C4(:,3),'b.');
label_plot('4: Vertices',viewing_angle)

% store all of the points in a single matrix, which we will need later for
% the triangularization function. Note that instead of simply concatinating
% the C matrices, i.e. P1 = [C1;C2;C3;C4], we are re-ordering the points to
% loop around vertical "slices" of the shape (e.g. the first slice is
% points 1, 2, 3, 4; the second slice is points 5, 6, 7, 8...). This makes
% defining the vertices of the triangles easier later. 

n = size(C1,1);     % number of points in a single c curve
N = n*4;            % total number of points in the entire shape

for i=1:n
    P1(4*i-3,:) = C1(i,:);
    P1(4*i-2,:) = C2(i,:);
    P1(4*i-1,:) = C3(i,:);
    P1(4*i,:) = C4(i,:);
end

disp(' ')
disp('Run the next section draw the solid 3D shape using the first method.')
disp('This method draws triangles across the top and bottom surfaces of the shape,')
disp('connecting the opposite edges.')


%% SECTION 5: Plot Triangularization (Method 1)

% break the surface up into triangles for export as an STL

% method 1: vertical triangles on the side walls, triangles going across
% the top and bottom surfaces to opposite edges

% move along each "slice" of the shape and define the triangles connecting
% this slice to the next slice. Each section will have 4 sides (top,
% bottom, left, and right) and each side is broken into 2 triangles, for 8
% total triangles per section. Think of this like doing "connect the dots,"
% where each triangle is defined by its 3 vertices, and your goal is to
% fill in the entire surface of the shape.

for i=1:n-1
     
    % define the indices of the points that make up the corners of this
    % slice of the shape. the index is the row number of the P1 matrix that 
    % contains the x,y,z coordinates for that point.
    % note: "left" and "right" here are defined when viewing the shape with
    % the positive x-axis pointing away from you, the positive y-axis to
    % your left, and the positive z-axis pointing up
    BL1 = 4*i-3;  % bottom left
    BR1 = 4*i-2;  % bottom right
    TL1 = 4*i-1;  % top left
    TR1 = 4*i;    % top right

    % define corners of next slice
    BL2 = 4*(i+1)-3;  % bottom left
    BR2 = 4*(i+1)-2;  % bottom right
    TL2 = 4*(i+1)-1;  % top left
    TR2 = 4*(i+1);    % top right

    % define the "triangulation connectivity matrix" T1. T1 is a matrix
    % where each row represents one triangle. Each triangle is defined by
    % the indices of its three vertices. The vertices are defined in order
    % using the right-hand rule so the surface normal of the triangle
    % points outward (not into the interior of the shape).

    T1(8*i-7,:) = [BL1, TL1, BL2]; % triangle 1 (left face) 
    T1(8*i-6,:) = [BL2, TL1, TL2]; % triangle 2 (left face)
    T1(8*i-5,:) = [TL1, TR1, TL2]; % triangle 1 (top face)
    T1(8*i-4,:) = [TL2, TR1, TR2]; % triangle 2 (top face)
    T1(8*i-3,:) = [TR2, BR1, BR2]; % triangle 1 (right face)
    T1(8*i-2,:) = [TR1, BR1, TR2]; % triangle 2 (right face)
    T1(8*i-1,:) = [BR1, BL1, BR2]; % triangle 1 (bottom face)
    T1(8*i,:)   = [BR2, BL1, BL2]; % triangle 2 (bottom face)

end

% add "end caps" at the beginning and end of the shape. Note that if the
% original 2D curves a1 and a2 have the same endpoints, this will result in
% triangles of zero width, so they are technically not needed, but
% keeping this code will make sure you still get a solid shape without any
% gaps if a1 and a2 do not have the same endpoints.

i = 8*i+1;
T1(i,:) = [1,4,3];
i = i+1;
T1(i,:) = [1,2,4];
i = i+1;
T1(i,:) = [N-2,N-1,N];
i = i+1;
T1(i,:) = [N-2,N-3,N-1];

% use matlab's built-in triangularization function with our connectivity
% matrix T1 and our matrix of points P1
TRI1 = triangulation(T1,P1);

% make a surface plot of the triangularization
subplot(2,3,5)
trisurf(TRI1,'FaceColor',facecolor)
label_plot('5: STL Method 1 (triangles across)',viewing_angle)

disp(' ')
disp('Run the next section draw the solid 3D shape using the second method.')
disp('This method connects triangles on the top and bottom surfaces to')
disp('a central point in the middle of the surface.')

%% SECTION 6: Plot Triangularization (Method 2)

% STL method 2: instead of making triangles go across to opposite edges of
% the top and bottom of the shape, we will add two more points in the
% center of the top and bottom, and make triangles that connect to these
% points.

TC = N+1; % index of top center point
BC = N+2; % index of bottom center point
P2 = P1;  % create copy of points matrix P1
zavg = 0.5*(mean(c1z)+mean(c2z)); % calculate average z value of lower curve, this might not be zero
P2 = [P2;0 0 zavg+height;0 0 zavg]; % add two more rows to store xyz coordinates of top and bottom center points

% move along each slice of the shape, this time defining the top and bottom
% triangles from the edges to the center point

for i=1:n-1
    
    % define corners of this slice
    BL1 = 4*i-3;  % bottom left
    BR1 = 4*i-2;  % bottom right
    TL1 = 4*i-1;  % top left
    TR1 = 4*i;    % top right

    % define corners of next slice
    BL2 = 4*(i+1)-3;  % bottom left
    BR2 = 4*(i+1)-2;  % bottom right
    TL2 = 4*(i+1)-1;  % top left
    TR2 = 4*(i+1);    % top right

    % do left and right faces like before

    T2(8*i-7,:) = [BL1, TL1, BL2]; % triangle 1 (left face) 
    T2(8*i-6,:) = [BL2, TL1, TL2]; % triangle 2 (left face)
    T2(8*i-5,:) = [TR2, BR1, BR2]; % triangle 1 (right face)
    T2(8*i-4,:) = [TR1, BR1, TR2]; % triangle 2 (right face)

    % do top and bottom faces
    T2(8*i-3,:) = [TC,TL2,TL1]; % triangle 1 (top face)
    T2(8*i-2,:) = [TC,TR1,TR2]; % triangle 2 (top face)
    T2(8*i-1,:) = [BC,BL1,BL2]; % triangle 1 (bottom face)
    T2(8*i,:) = [BC,BR2,BR1]; % triangle 2 (bottom face)

end

% this section will result in triangles of zero thickness if the original
% curves a1 and a2 have the same endpoints

% add vertical end caps at the front and back of the shape
i = 8*i+1;
T2(i,:) = [1,4,3];
i = i+1;
T2(i,:) = [1,2,4];
i = i+1;
T2(i,:) = [N-2,N-1,N];
i = i+1;
T2(i,:) = [N-2,N-3,N-1];
i = i+1;

% add the first and last triangles on the top and bottom of the shape
% (these triangles go across the x-axis and are not included by the for
% loop above
T2(i,:) = [BC,2,1];
i = i+1;
T2(i,:) = [BC,N-3,N-2];
i = i+1;
T2(i,:) = [TC,3,4];
i = i+1;
T2(i,:) = [TC,N,N-1];

% use matlab's built-in triangularization function with our connectivity
% matrix T2 and our matrix of points P2
TRI2 = triangulation(T2,P2);

% make a surface plot of the triangularization
subplot(2,3,6)
trisurf(TRI2,'FaceColor',facecolor)
label_plot('6: STL Method 2 (triangles to center)',viewing_angle)

disp(' ')
disp('Run the final section of code to save STL files for both methods.')

%% SECTION 7: Export STL files

% export your STL files for viewing in a CAD program or 3D printing

STL_filename = 'myshape';           % name for the files - change this if you want a different name
ext = '.stl';                       % file extension - don't change this
file_num = 0;                       % file number to append to the file name (don't change this), it will
                                    % automatically increment using the code below

% check for existing files with file_num in the folder, if the file number
% exists, increment the number so we don't overwrite the existing files
while isfile(strcat(STL_filename,'_method1_',num2str(file_num),ext))
    file_num = file_num + 1;
end

% save one STL file for each triangularization method
stlwrite(TRI1,strcat(STL_filename,'_method1_',num2str(file_num),ext))
stlwrite(TRI2,strcat(STL_filename,'_method2_',num2str(file_num),ext))

disp(' ')
disp('The STL files should be saved in the same directory as this m-file.')
disp('You can import these STL files into a CAD program or 3D print them.')
disp('If you run the program again, it will automatically increment the number')
disp('at the end of the STL filename, so you do not need to worry about')
disp('overwriting the files.')

%% SECTION 8: Functions

% function to set up and label 3D plots
function [] = label_plot(plot_title,angle)
    axis vis3d                 % lock axis scaling for 3D rotation
    axis equal                 % set axes to equal scale
    grid on                    % turn grid on
    xlabel('x')                % label x axis
    ylabel('y')                % label y axis
    zlabel('z')                % label z axis
    title(plot_title)          % give the graph a title
    view(0,angle)              % set the view angle
end