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 points, and the canvas width and height are stored in
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 = ; 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.
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. 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
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. Therefore, the resulting triangles usually do not have extremely acute angles.
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 .
- 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);
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))