Cris’ Image Analysis Blog

Several algorithms require a set of uniformly distributed points across the image. For example, superpixel algorithms typically start with a regular grid. A rectangular grid of points is trivial to draw into an image. One simply generates the set of points of a grid that covers the image, rounds those points to integer coordinates, and sets those pixels. A rotated grid is a bit more challenging–one needs to compute bounds of the rotated grid such that the full image is covered. With other than rectangular grids (for example a hexagonal grid) the math gets a bit more complicated, but this can still be computed. But how to generalize such an algorithm to three dimensions? And to an arbitrary number of dimensions?

I played around with this task for a bit and came up with a very simple, general solution.

First, let’s review how points on a grid can be generated in general. A grid is defined by a set of basis vectors. These vectors do not need to be orthogonal. For example, for an hexagonal grid, the basis vectors are v1 = [1,0] and v2 = [cos(60°),sin(60°)] (using degrees, not radian, for simplicity of notation). Given these two vectors, a*v1 + b*v2, with a and b integers, yields the grid points.

This notion generalizes trivially to any number of dimensions. For example, in 3D the FCC (face-centered cubic) grid is given by vectors v1 = [1,1,0], v2 = [1,0,1] and v3 = [0,1,1]. There are always as many vectors as number of dimensions. These vectors form the columns of a matrix M, such that M*s is a grid point, with s a column vector of integer values. This s can be thought of as the grid index.

So the difficulty is to determine which integer-valued vectors p lead to grid points that intersect the image. For a rectangular grid that is aligned with the axes, it is also a rectangular domain of s that needs to be considered. When rotating a 2D grid, it is still quite simple to determine the larger rectangular domain that completely covers the image, then discard the points that do not intersect the image. The same can be said for a non-rectangular grid in 2D, it is always possible to find the rectangular domain of s that covers the whole image. But these computations are more complex in the 3D case, and quite difficult to generalize to an arbitrary number of dimensions. For higher-dimensional spaces, these domains of s become increasingly large, making the point discarding solution increasingly inefficient.

Instead, we can start with the image domain itself, and for each pixel find out if it is the pixel nearest a grid point. For this, we first need to find the inverse of the matrix M, let’s call it inv_M. For a given pixel p = [x,y,z,...] the expression s = round(inv_M*p) is an integer grid index. As we saw before, q = M*s is the coordinates of the grid point for grid index s. round(q) is the nearest pixel. If round(q)==p, then the pixel at p should be marked as a grid point.

Using DIPimage 3, we can write as follows (this code will not work in DIPimage 2.9 or earlier without some changes):

% Define image domain
sz = [256,140];
% Define grid
d = 10.2;
M = [[1;0],[cosd(60);sind(60)]] * d;
% Let's rotate the grid too
phi = pi/12;
M = [cos(phi),sin(phi);-sin(phi),cos(phi)] * M;
% The inverse of the grid matrix
inv_M = inv(M);
% Convert M and inv_M to dip_image objects so we can do the multiplication
M = dip_image(M(:),[2,2]);          % Create single-pixel tensor image
inv_M = dip_image(inv_M(:),[2,2]);  % ...yes, this syntax is a bit awkward still, hopefully this will be easier soon
% Do the computation
p = coordinates(sz);                % An image where each pixel's value is it's own coordinate vector
s = round(inv_M * p);
q = round(M * s);
grid = all(q == p,'tensor');

Note that in the code above, we determined the image domain sz as 2D, and defined the grid matrix M as 2D. Other than that, there is no mention of the dimensionality. Changing these two matrices we can instead generate a grid in a 3D image, or a 5D image.

In case you don’t want to bother installing DIPimage, below is a version of the code that uses only standard MATLAB functions. This version, however, is explicitly for 2D.

% Define image domain
sz = [256,140];
% Define grid
d = 10.2;
M = [[1;0],[cosd(60);sind(60)]] * d;
% Let's rotate the grid too
phi = pi/12;
M = [cos(phi),sin(phi);-sin(phi),cos(phi)] * M;
% The inverse of the grid matrix
inv_M = inv(M);
% Do the computation
[x,y] = meshgrid(1:sz(1),1:sz(2));
p = [x(:).'; y(:).'];
s = round(inv_M * p);
q = round(M * s);
grid = all(q == p, 1);
grid = reshape(grid, size(x));