Question

Draw a Koch Snowflake fractal image using MATLAB programming.

I'm pretty new to MATLAB and im not too familar with how all the plotting and drawing functions work.

Answer 1

Matlab CODE:

%KOCH: Plots 'Koch Curve' Fractal
% koch(n) plots the 'Koch Curve' Fractal after n iterations
% (be patient for n>8, depending on Computer speed)
%--------------------------------------------------------------------
function []=Koch(n)
if (n==0)
   x=[0;1];
   y=[0;0];
   line(x,y,'Color','r');
   axis equal
   set(gca,'Visible','off')
else
   levelcontrol=10^n;
   L=levelcontrol/(3^n);
   l=ceil(L);
   kline(0,0,levelcontrol,0,l);
   axis equal
   set(gca,'Visible','off')
   set(gcf,'Name','Koch Curve')
end
%--------------------------------------------------------------------
function kline(x1,y1,x5,y5,limit)
length=sqrt((x5-x1)^2+(y5-y1)^2);
if(length>limit)
   x2=(2*x1+x5)/3;
   y2=(2*y1+y5)/3;
   x3=(x1+x5)/2-(y5-y1)/(2.0*sqrt(3.0));
   y3=(y1+y5)/2+(x5-x1)/(2.0*sqrt(3.0));
   x4=(2*x5+x1)/3;
   y4=(2*y5+y1)/3;
   % recursive calls
   kline(x1,y1,x2,y2,limit);
   kline(x2,y2,x3,y3,limit);
   kline(x3,y3,x4,y4,limit);
   kline(x4,y4,x5,y5,limit);
else
   plotline(x1,y1,x5,y5);
end
%--------------------------------------------------------------------
function plotline(a1,b1,a2,b2)
x=[a1;a2];
y=[b1;b2];
line(x,y);
%--------------------------------------------------------------------

Result:

Command: Koch(3)

if you want full curve, use the below code:

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
nn= input('Enter number of iterations needed : ');
for iteration = 1:nn
    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);
%------------------------------------------------------------------

Save this function in your Matlab location,

function p_vec = perp(vec)
% Returns the unit vector p_vec, which is perpendicular to the input
% 2D vector vec
p_vec(1) = -vec(2);
p_vec(2) = vec(1);

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

end

________________________________________

Result: