This is the fourth article in my "MATLAB Plots Can Be Beautiful" series. To view the previous three articles, check the following links:

- MATLAB Plots Can Be Beautiful I - Drawing a Tree
- MATLAB Plots Can Be Beautiful II - Mountains
- MATLAB Plots Can Be Beautiful III - Light Dots

In this article, we will again use MATLAB to generate wallpaper quality plots. More specifically, we will generate plots similar to the following:

At first glimpse, the above picture looks more like a Photoshop product than a MATLAB plot. Upon further inspection, we can observe that the above picture consists of only triangles with different colors. Recall that we can indeed plot filled triangles in MATLAB. Hence, it is completely plausible to generate the above picture with MATLAB's plot functions. Actually, the generation procedure is quite simple, and can be summarized in three steps:

- Sample some random points in a 2D canvas.
- Triangulate the canvas region with these points.
- Color each triangle.

Next, we will cover the implementation of each step in detail.

## Poisson-Disc Sampling and Random Points Generation

The first step is to sample random points in a 2D canvas. The most straightforward sampling method is uniform sampling. Suppose we need to sample $n$ points, and the canvas width and height are stored in `canvas_width`

and `canvas_height`

, respectively. Basically, we can use `unifrnd(0, canvas_width, n, 1)`

to sample the x-coordinates and `unifrnd(0, canvas_height, n, 1)`

to sample the y-coordinates. However, this approach has a problem: it may generate two points that are too close to each other, and we cannot enforce the minimal distance between any too points. If we triangulate the canvas using these points, the result does not look very natural.

To tackle this issue, we will use *Poisson-disc sampling*, which will guarantee the minimal distance between any two points, leading to a more nature pattern. We will use Bridson's algorithm to generate the sample points. More details can be found in this one-page paper. Below is a MATLAB implementation for the 2D case:

```
function points = poisson_disc_2d(w, h, n, r, k)
%POISSON_DISC_2D Generate 2D points using the Bridson's Poisson-disc
%sampling algorithm.
%Inputs:
% w - Width of the region.
% h - Height of the region.
% n - Maximum number of points.
% r - Radius of the disc.
% k - (Optional) Maximum number of trials. Default value is 30.
%Output:
% points - n x 2 matrix where each row represents a point. n may be less
% than the specified n.
if nargin < 5
k = 30;
end
% Background grid for fast occupancy search.
grid_size = r / sqrt(2);
n_w = ceil(w / grid_size);
n_h = ceil(h / grid_size);
grid = zeros(n_w, n_h);
% Generate the initial point.
points = zeros(n, 2);
points(1, :) = [unifrnd(0, w) unifrnd(0, h)];
active_points = [1];
n_points_generated = 1;
grid(floor(points(1,1) / grid_size) + 1, floor(points(1,2) / grid_size) + 1) = 1;
% Generate the remaining points
while ~isempty(active_points) && n_points_generated < n
% Pick a random point from the active set.
idx = randi(size(active_points, 1));
p = points(active_points(idx), :);
% Try generate a new point from this point within the annulus [r, 2r].
success = false;
for kk = 1:k
% Uniformly samples a point within an annulus.
new_r = sqrt(r * r + 3 * r * r * rand());
new_theta = unifrnd(0, pi*2);
new_p = [new_r * cos(new_theta) + p(1)...
new_r * sin(new_theta) + p(2)];
% Check validity.
if new_p(1) < 0 || new_p(2) < 0 || new_p(1) > w || new_p(2) > h
continue;
end
grid_x = floor(new_p(1) / grid_size) + 1;
grid_y = floor(new_p(2) / grid_size) + 1;
if grid_x <= 0 || grid_y <= 0 || grid_x > n_w || grid_y > n_h
continue;
end
if grid(grid_x, grid_y) > 0
continue;
end
% Check eight adjacent grids.
dirs = [-1 -1;-1 0;-1 1;0 -1;0 1;1 -1;1 0;1 1];
valid = true;
for ii = 1:8
new_grid_x = grid_x + dirs(ii, 1);
new_grid_y = grid_y + dirs(ii, 2);
if new_grid_x == 0 || new_grid_y == 0 || new_grid_x > n_w || new_grid_y > n_h
continue;
end
if grid(new_grid_x, new_grid_y) == 0
continue;
end
p = points(grid(new_grid_x, new_grid_y), :);
if sqrt((new_p(1) - p(1)).^2 + (new_p(2) - p(2)).^2) < r
valid = false;
break;
end
end
if ~valid
continue;
end
% Add this point.
n_points_generated = n_points_generated + 1;
grid(grid_x, grid_y) = n_points_generated;
active_points = [active_points; n_points_generated];
points(n_points_generated, :) = new_p;
success = true;
break;
end
% Remove this point if the generation fails.
if ~success
active_points(idx,:) = [];
end
end
if n_points_generated < n
points = points(1:n_points_generated, :);
end
end
```

Note:the above implementation is not optimized. The line`active_points(idx,:) = []`

takes linear time and can be further optimized. It is also possible to use MEX to develop an even faster implementation.

The two figures below show the difference between uniform sampling and Poisson-disc sampling. We can observe that Poisson-disc sampling produces a more nature pattern for triangulation.

With our implementation of Poisson-disc sampling, we can write down the MATLAB code for generating the points as follows:

```
x_min = -0.05 * canvas_width;
x_max = 1.05 * canvas_width;
y_min = -0.05 * canvas_height;
y_max = 1.05 * canvas_height;
points = poisson_disc_2d(canvas_width * 1.1, canvas_height * 1.1, max_n_points, r);
% Shift to center.
x = points(:, 1);
y = points(:, 2);
x = x - 0.05 * canvas_width;
y = y - 0.05 * canvas_height;
% Fix four corners.
x = [x; x_min; x_min; x_max; x_max];
y = [y; y_min; y_max; y_min; y_max];
n_points = size(x, 1);
```

Note that to avoid artifacts cases by boundary points, we sample the points in a slightly larger area. We also manually put four corner points just in case the sampling algorithm fails to generate points near the corners. We need to adjust `max_n_points`

according to our choice of `r`

, the minimal distance between any two points.

## Triangulation

Given the generated points, we can now triangulate the 2D canvas. Luckily, MATLAB has a built-in function named `delaunay`

that implements the Delaunay triangulation^{[1]}. The documentation can be found here. Therefore, the triangulation step can be done with a single function call:

```
% Triangulation
% =============
tri = delaunay(x, y);
```

Each row of the variable `tri`

represents a triangle defined by indices with respect to the points defined by the coordinate vectors `x`

and `y`

. For instance, if `tri(1,:) = [1 5 3]`

, it means the first triangle's three vertices are given by the first, the fifth, and the third sample points. The figure below shows what these triangles look like:

One interesting characteristic of the Delaunay triangulation is that it maximizes the minimum angle of all the angles of the triangles in the triangulation^{[1]}. Therefore, the resulting triangles usually do not have extremely acute angles.

## Drawing

We now arrive at the final part: we need to draw and color these triangles. Because the indices of the vertices are already stored in `tri`

, we can simply use `patch`

to draw all the triangles in one batch. More details about the `patch`

function can be found here. Before calling `patch`

, we still have one more thing to do: determine the color of each triangle. Here, to make our implementation more flexible, we do not directly compute the colors. Instead, we compute a scalar value for each triangle, and then use `colormap`

to map these scalar values to the desired colors (similar to `imagesc`

). These scalar values are calculated as follows:

- Compute the centroids of these triangles.
- Evaluate the pattern function which outputs a scalar for each centroid.
- Normalize the outputs of the pattern function to $[0, 1]$.
- Add some Gaussian noise for a little bit of variations.

In MATLAB, the above steps are implemented as follows:

```
% Compute the centroids of the triangles.
c_x = mean(x(tri), 2);
c_y = mean(y(tri), 2);
% Evaluate the value for each face.
cc = pattern_func(c_x, c_y);
cc = cc - min(cc);
cc = cc / max(cc);
cc = cc + variation * randn(size(cc, 1), 1);
```

The simplest `pattern_func`

is the linear one: `x + y`

. You can generate more interesting color patterns by tweaking this function. Check the last part of this article for some interesting pattern functions.

We can now call `patch`

to draw these triangles:

```
% Use patch to draw the triangles.
patch('Faces', tri, 'Vertices', [x y], 'FaceVertexCData', cc, ...
'FaceColor', 'flat', 'EdgeColor', 'none');
```

Finally, we color then using a custom color map:

```
h = unifrnd(0, 1) * ones(n_colors, 1);
s = linspace(1, 0.2, n_colors)';
v = linspace(0, 1, n_colors)';
colormap(hsv2rgb([h s v]));
```

To save the plot as a high-resolution picture, simply use:

```
set(gcf, 'PaperPositionMode', 'auto');
print(gcf, 'triangles.png', '-dpng', '-r300');
```

## Putting Everything Together

Now we can put everything together and finalize our MATLAB script:

```
clear();
% Configurations
% ==============
canvas_width = 1080;
canvas_height = 720;
r = 40;
n_colors = 100;
variation = 0.02;
max_n_points = floor(canvas_height * canvas_width / r / r) + 1;
pattern_func = @(x, y) x + y;
% Points Generation
% =================
x_min = -0.05 * canvas_width;
x_max = 1.05 * canvas_width;
y_min = -0.05 * canvas_height;
y_max = 1.05 * canvas_height;
points = poisson_disc_2d(canvas_width * 1.1, canvas_height * 1.1, max_n_points, r);
% Shift to center.
x = points(:, 1);
y = points(:, 2);
x = x - 0.05 * canvas_width;
y = y - 0.05 * canvas_height;
% Fix four corners.
x = [x; x_min; x_min; x_max; x_max];
y = [y; y_min; y_max; y_min; y_max];
n_points = size(x, 1);
% Triangulation
% =============
tri = delaunay(x, y);
% Drawing
% =======
hf = figure('Position', [50 50 canvas_width canvas_height]);
axis('equal');
set(gca, 'Position', [0 0 1 1]);
axis([0 canvas_width 0 canvas_height]);
axis('off');
% Compute the centroids of the triangles.
c_x = mean(x(tri), 2);
c_y = mean(y(tri), 2);
% Evaluate the value for each face.
cc = pattern_func(c_x, c_y);
cc = cc - min(cc);
cc = cc / max(cc);
cc = cc + variation * randn(size(cc, 1), 1);
% Use patch to draw the triangles.
patch('Faces', tri, 'Vertices', [x y], 'FaceVertexCData', cc, ...
'FaceColor', 'flat', 'EdgeColor', 'none');
% Apply the color map.
h = unifrnd(0, 1) * ones(n_colors, 1);
s = linspace(1, 0.2, n_colors)';
v = linspace(0, 1, n_colors)';
colormap(hsv2rgb([h s v]));
```

## Bonus: Changing Pattern Functions

We can change the pattern function `pattern_func`

to generate more interesting color patterns:

`(x - canvas_width / 2).^2 - (y - canvas_height / 2).^2`

`sin(1.5 * pi * x / canvas_width) .* cos(1.1 * pi * y / canvas_height)`

`abs(x - canvas_width / 2 + canvas_width / canvas_height * (y - canvas_height / 2))`