# Plotting the Koch Snowflake

The Koch snowflake is a well known fractal. Beginning with an equilateral triangle, a smaller equilateral triangle is placed halfway along each edge of the shape. This process then repeats on each edge of the new shape. Here is my attempt to plot this in MATLAB.

clc, clearvars, close all;

length = 1; % Original side length of triangle
theta = 60*pi/180; % Angle of the equilateral triangle in radians

% Points of the original equilateral triangle
p1 = [0,0];
p2 = [1/2*length,length*sin(theta)];
p3 = [length,0];

flake_points = [p1; p2; p3]; % Vertices of the snowflake

for iteration = 1:7
n_points = size(flake_points,1); % No. points currently in data
new_triangle_points = nan(n_points*3,2); % New points to be added
for i = 1:n_points
p = flake_points(i,:);
if i == n_points
vector = flake_points(1,:)-flake_points(i,:);
else
% Vector between next point and current point
vector = flake_points(i+1,:)-flake_points(i,:);
end
b1 = p + 1/3*vector; % Base point 1 of the new triangle
b2 = p + 2/3*vector; % Base point 2 of the new triangle
top = p + 1/2*vector + perp(vector)*length/3*sin(theta);
new_triangle_points(3*i-2:3*i,:) = [b1; top; b2];
end
flake_points(end+1:end+9,:) = nan(9,2); % Add rows to data
new_flake_points = nan(n_points*3+n_points,2);
n = 1;
for i = 1:n_points
new_flake_points(n,:) = flake_points(i,:);
index = 3*i-2;
new_flake_points(n+1:n+3,:) = new_triangle_points(index:index+2,:);
n = n + 4; % Three new points have been inserted,
% so we now skip to the next empty spaces
end
flake_points = new_flake_points;
length = 1/3*length; % Revise length
end

plot(flake_points(:,1),flake_points(:,2),'w')
hold on; % Plot last data point with 1st (ie, close the loop)
plot([flake_points(1,1),flake_points(end,1)],...
[flake_points(1,2),flake_points(end,2)],'w')
set(gca,'color',[0.1 0.1 0.1])
axis equal;
str = sprintf('Iteration: %d',iteration);
title(str, 'FontSize', 20);


Output: Edit - This function is also required:

function p_vec = perp(vec)
% Returns the unit vector p_vec, which is perpendicular to the input
% 2D vector vec
%
% Inputs
%
%   vec : double
%       n_dimensions x 1
%       Input vector
%
% Outputs
%
%   p_vec : double
%       n_dimensions x 1
%       Perpendicular unit vector

p_vec(1) = -vec(2);
p_vec(2) = vec(1);

% Convert to unit vector
p_vec = p_vec / norm(p_vec);

end

• top = p + 1/2*vector + perp(vector)*length/3*sin(theta); What does this perp(vector) ? – lola Oct 27 '17 at 8:11
• @Iola my mistake, I must have written a function without including it. It returns the unit vector p_vec, which is perpendicular to the input 2D vector – Vermillion Oct 29 '17 at 22:26
• i am trying to understand your code but your code references a 'perp' function I am assuming it finds the perpendicular vector which must satisfy the dot product of such a vector which when under dot product operation will yield zero for the vector that is being inputted into the equation, vector * perpendicularvector=o by any chance do you have a copy of this function? – Adrian Mar 4 '18 at 15:48
• @Adrian I've edited the post now to include this function. – Vermillion Mar 6 '18 at 2:19

Very nice script!

I suggest you put this into a function, that takes n as an input argument instead of hardcoding it. It requires very little alterations, and will make the function a lot easier to use. You can skip all the clearing in the beginning and substitute it with:

function [] = koch_snowflake(iterations)  % or


and add an end at the bottom of your script.

length = 1; % Original side length of triangle


This is a no no. I strongly recommend you not to use length as a variable name as it's a very commonly used builtin function. This can cause a lot of errors that are hard to find. Note that you should always use the function numel instead of length, but it's still not a good idea to overload it.

theta = 60*pi/180; % Angle of the equilateral triangle in radians


I guess you think 60*pi/180 is easier to read than pi/3, but I disagree. If you want to use the number of degrees, and convert it to radians, you might want to use deg2rad.

I suggest substituting the iterator i to for instance ii if you ever work with complex numbers, to avoid potential errors and confusion. i is by default the imaginary unit in MATLAB (so is j, 1i, 1j). This is a matter of preference though.

Also, personal preference again, but I like it when comments are aligned (if possible):

length = 1;                   % Original side length of triangle
theta = 60*pi/180;            % Angle of the equilateral triangle in radians

% Points of the original equilateral triangle
p1 = [0,0];
p2 = [1/2*length,length*sin(theta)];
p3 = [length,0];

flake_points = [p1; p2; p3];  % Vertices of the snowflake

for iteration = 1:7
n_points = size(flake_points,1);          % No. points currently in data
new_triangle_points = nan(n_points*3,2);  % New points to be added
...
...