Recently, I came across the beautiful magic tree on OpenProcessing (check out processing.org if you have no idea what it is). As a daily MATLAB user, I asked my self a question: is it possible to recreate something similar using MATLAB? After some trials, I managed to get MATLAB producing something with similar looking. Here is what it looks like:

Before jumping into the code, let us make a checklist about MATLAB's drawing capabilities:

- Can we draw lines (for the trunk and branches) in MATLAB? Yes.
Use the
`line`

function. - Can we draw circles (for the leaves) in MATLAB? Yes. Use the
`rectangle`

function (yes, in MATLAB you use the`rectangle`

function to draw circles). - Does MATLAB support anti-aliasing? Yes. Since 2014b, MATLAB introduced major updates in its graphics engine, and the resulting figures looks much nicer.
- Does MATLAB support transparency (so I can get nicer looking leaves)? Yes.
Since 2014b, alpha value can be specified in the
`color`

property.

Now we know we need to use MATLAB 2014b or newer to get some good looking results. Let us start by setting up the figure and defining the parameters:

```
function plot_tree
% Plots a beautiful looking tree.
% set up figure
figure;
canvas_width = 1000;
canvas_height = 800;
axis('equal');
axis([0 canvas_width 0 canvas_height]);
% thickness of the tree trunk
base_thickness = 16;
% length of the tree trunk
base_length = 200;
% determines how fast the length and thickness of the branches decay
branch_decay_factor = 0.680;
% determines the randomness of branch_decay_factor
branch_decay_randomness = 0.3;
% determines the direction of the new branch, 0 = same direction
branch_angle = deg2rad(25);
% determines the randomness of branch_angle
branch_angle_randomness = 0.4;
% color of the tree trunk and branches
branch_color = [51 25 0]/255;
% minimum length of the branches
min_branch_length = base_length * branch_decay_factor^7;
% determines how far the leaves can deviate from the branch tips
leaf_offset = 8;
% gives the leaves a random color
leaf_color = hsv2rgb([unifrnd(0, 1) 1 1]);
% drawing logic
% nested functions
end
```

Note that we added lots of randomness related parameters. Without randomness, the tree will look pretty dull. Next let us outline the drawing logic.

```
% drawing logic
leaf_locations = zeros(1024, 2);
leaf_count = 0;
% draw trunk, adding an extra y offset to make the junction look better
line([canvas_width/2 canvas_width/2], [0 100+base_thickness], ...
'LineWidth', base_thickness, 'Color', branch_color);
% draw branches
draw_branch_at(canvas_width/2, 100, ...
branch_decay_factor * base_thickness, branch_decay_factor * base_length, 1j);
% draw leaves
leaf_locations = leaf_locations(1:leaf_count, :);
% draw three layers of leaves to get a nicer look
draw_leaves(10, 0.2, 0.6, 0.2);
draw_leaves(40, 0.3, 0.1, 0.3);
draw_leaves(100, 0.6, 0.03, 0.1);
```

We first preallocate a matrix to store the leaf locations (actually we do not
know how many leaves will be generated because of randomness). We then draw the
tree trunk, the branches, and finally the leaves. The `draw_leaves`

function
is called three times with different parameters to get a more interesting
look. The whole procedure can be visualized by the following four figures.

Now we need to implement the nested functions `draw_branch_at`

and `draw_leaves`

. Let us first take a look
at `draw_branch_at`

, which recursively calls itself to draw all the branches:

```
function can_branch = draw_branch_at(x, y, thickness, len, direction)
% check length
if len < min_branch_length
% add leaf here
leaf_count = leaf_count + 1;
leaf_locations(leaf_count, :) = [x y];
can_branch = false;
return;
end
% generate branch angles
theta_1 = branch_angle * (1.0 + branch_angle_randomness * max(min(randn(), 1), -1));
theta_2 = -branch_angle * (1.0 + branch_angle_randomness * max(min(randn(), 1), -1));
% we use a unit length complex number to represent the directional
% vector here, and rotation is done via complex number multiplication
new_dir_1 = direction * (cos(theta_1) + 1j*sin(theta_1));
new_dir_2 = direction * (cos(theta_2) + 1j*sin(theta_2));
% compute end locations of the new branches
new_x_1 = x + real(new_dir_1) * len;
new_y_1 = y + imag(new_dir_1) * len;
new_x_2 = x + real(new_dir_2) * len;
new_y_2 = y + imag(new_dir_2) * len;
% when drawing, draw a litte longer to make the junctions look better
extra_x_1 = real(new_dir_1) * thickness;
extra_y_1 = imag(new_dir_1) * thickness;
extra_x_2 = real(new_dir_2) * thickness;
extra_y_2 = imag(new_dir_2) * thickness;
% draw the two branches
line([x new_x_1+extra_x_1], [y new_y_1+extra_y_1], ...
'LineWidth', thickness, 'Color', branch_color);
line([x new_x_2+extra_x_2], [y new_y_2+extra_y_2], ...
'LineWidth', thickness, 'Color', branch_color);
% make new branches
decay = branch_decay_factor * (1.0 + branch_decay_randomness * max(min(randn(), 1), -1));
new_length = decay * len;
new_thickness = decay * thickness;
if ~draw_branch_at(new_x_1, new_y_1, new_thickness, new_length, new_dir_1) || ...
~draw_branch_at(new_x_2, new_y_2, new_thickness, new_length, new_dir_2)
% we have already branched here, but cannot further branch
% add a leaf at the starting location of the current branch
leaf_count = leaf_count + 1;
leaf_locations(leaf_count, :) = [x y];
end
can_branch = true;
end
```

We first check if the current branch length is less than `min_branch_length`

.
If so we add a new leave and return. Otherwise we draw the current branch and
add two new branches by recursively calling `draw_branch_at`

.

Compared with `draw_branch_at`

, the implementation of `draw_leaves`

is pretty
straightforward:

```
function draw_leaves(leave_size, leave_size_randomness, alpha, alpha_randomness)
for ii = 1:leaf_count
% add randomness to size and transparency
cur_size = leave_size * (1.0 + leave_size_randomness * max(min(randn(), 1), -1));
cur_alpha = alpha * (1.0 + alpha_randomness * randn());
cur_alpha = max(min(cur_alpha, 1), 0);
cur_color = [leaf_color cur_alpha];
% compute and add variance to the location of the current leaf
loc = leaf_locations(ii,:) - cur_size/2;
theta = unifrnd(0, 2*pi);
offset = unifrnd(0, leaf_offset);
loc = loc + [cos(theta) sin(theta)] * offset;
% draw the current leaf
rect = [loc cur_size cur_size];
rectangle('Position', rect, 'FaceColor', cur_color, ...
'EdgeColor', 'none', 'LineWidth', 1, 'Curvature', [1 1]);
end
end
```

Putting them together:

```
function plot_tree
% Plots a beautiful looking tree.
% set up figure
figure;
canvas_width = 1000;
canvas_height = 800;
axis('equal');
axis([0 canvas_width 0 canvas_height]);
% thickness of the tree trunk
base_thickness = 16;
% length of the tree trunk
base_length = 200;
% determines how fast the length and thickness of the branches decay
branch_decay_factor = 0.680;
% determines the randomness of branch_decay_factor
branch_decay_randomness = 0.3;
% determines the direction of the new branch, 0 = same direction
branch_angle = deg2rad(25);
% determines the randomness of branch_angle
branch_angle_randomness = 0.4;
% color of the tree trunk and branches
branch_color = [51 25 0]/255;
% minimum length of the branches
min_branch_length = base_length * branch_decay_factor^7;
% determines how far the leaves can deviate from the branch tips
leaf_offset = 8;
% gives the leaves a random color
leaf_color = hsv2rgb([unifrnd(0, 1) 1 1]);
% drawing logic
leaf_locations = zeros(1024, 2);
leaf_count = 0;
% draw trunk, adding an extra y offset to make the junction look better
line([canvas_width/2 canvas_width/2], [0 100+base_thickness], ...
'LineWidth', base_thickness, 'Color', branch_color);
% draw branches
draw_branch_at(canvas_width/2, 100, ...
branch_decay_factor * base_thickness, branch_decay_factor * base_length, 1j);
% draw leaves
leaf_locations = leaf_locations(1:leaf_count, :);
% draw three layers of leaves to get a nicer look
draw_leaves(10, 0.2, 0.6, 0.2);
draw_leaves(40, 0.3, 0.1, 0.3);
draw_leaves(100, 0.6, 0.03, 0.1);
% nested functions
function can_branch = draw_branch_at(x, y, thickness, len, direction)
% check length
if len < min_branch_length
% add leaf here
leaf_count = leaf_count + 1;
leaf_locations(leaf_count, :) = [x y];
can_branch = false;
return;
end
% generate branch angles
theta_1 = branch_angle * (1.0 + branch_angle_randomness * max(min(randn(), 1), -1));
theta_2 = -branch_angle * (1.0 + branch_angle_randomness * max(min(randn(), 1), -1));
% we use a unit length complex number to represent the directional
% vector here, and rotation is done via complex number multiplication
new_dir_1 = direction * (cos(theta_1) + 1j*sin(theta_1));
new_dir_2 = direction * (cos(theta_2) + 1j*sin(theta_2));
% compute end locations of the new branches
new_x_1 = x + real(new_dir_1) * len;
new_y_1 = y + imag(new_dir_1) * len;
new_x_2 = x + real(new_dir_2) * len;
new_y_2 = y + imag(new_dir_2) * len;
% when drawing, draw a litte longer to make the junctions look better
extra_x_1 = real(new_dir_1) * thickness;
extra_y_1 = imag(new_dir_1) * thickness;
extra_x_2 = real(new_dir_2) * thickness;
extra_y_2 = imag(new_dir_2) * thickness;
% draw the two branches
line([x new_x_1+extra_x_1], [y new_y_1+extra_y_1], ...
'LineWidth', thickness, 'Color', branch_color);
line([x new_x_2+extra_x_2], [y new_y_2+extra_y_2], ...
'LineWidth', thickness, 'Color', branch_color);
% make new branches
decay = branch_decay_factor * (1.0 + branch_decay_randomness * max(min(randn(), 1), -1));
new_length = decay * len;
new_thickness = decay * thickness;
if ~draw_branch_at(new_x_1, new_y_1, new_thickness, new_length, new_dir_1) || ...
~draw_branch_at(new_x_2, new_y_2, new_thickness, new_length, new_dir_2)
% we have already branched here, but cannot further branch
% add a leaf at the starting location of the current branch
leaf_count = leaf_count + 1;
leaf_locations(leaf_count, :) = [x y];
end
can_branch = true;
end
function draw_leaves(leave_size, leave_size_randomness, alpha, alpha_randomness)
for ii = 1:leaf_count
% add randomness to size and transparency
cur_size = leave_size * (1.0 + leave_size_randomness * max(min(randn(), 1), -1));
cur_alpha = alpha * (1.0 + alpha_randomness * randn());
cur_alpha = max(min(cur_alpha, 1), 0);
cur_color = [leaf_color cur_alpha];
% compute and add variance to the location of the current leaf
loc = leaf_locations(ii,:) - cur_size/2;
theta = unifrnd(0, 2*pi);
offset = unifrnd(0, leaf_offset);
loc = loc + [cos(theta) sin(theta)] * offset;
% draw the current leaf
rect = [loc cur_size cur_size];
rectangle('Position', rect, 'FaceColor', cur_color, ...
'EdgeColor', 'none', 'LineWidth', 1, 'Curvature', [1 1]);
end
end
end
```

We can further improve the results by using gradient fill instead of solid color fill for the leaves. However, getting gradient fill working can be more tricky. Interested readers are directed to here.