MATLAB Plots Can Be Beautiful IV - Triangles

Date: 6-8 2018

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

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

example_thumbnail

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))

https://en.wikipedia.org/wiki/Delaunay_triangulation ↩ ↩

Mianzhi Wang