Sampling from surfaces of a 3D cube

Question

I have a task to sample a fixed number of points from a 3D cube, so that these points can be a point cloud. But I do not know 1) how to describe the surfaces of an arbitrary 3D cubes and 2) how to sample from that surface.

I can think of an approach that I first describe a unit cube with x = 1, x = -1, y = 1, y = -1, z = 1, z = -1, then sample from that cube first. Then I can scale and rotate every point. But then, what's the formula for rotating and scaling 3D objects? Is rotation for any single point just multiplying a 3-by-3 matrix? How about scaling?

I will also need to do these in a computational way so I can implement it in Python.

Any pointers to reference would be very helpful! I have not learnt geometry except for a course on vector calculus, so I have difficulty even find the correct term to Google my question. Thanks!

Answer 1

I am assuming you wish to sample points on the surface of a cuboid with vertices at $$\bbox{\begin{array}{l} \vec{v}_0 \\ \vec{v}_0 + \vec{v}_x \\ \vec{v}_0 + \vec{v}_y \\ \vec{v}_0 + \vec{v}_x + \vec{v}_y \\ \vec{v}_0 + \vec{v}_z \\ \vec{v}_0 + \vec{v}_x + \vec{v}_z \\ \vec{v}_0 + \vec{v}_y + \vec{v}_z \\ \vec{v}_0 + \vec{v}_x + \vec{v}_y + \vec{v}_z \\ \end{array}}$$ i.e. $\vec{v}_0$ is one of the vertices, and $\vec{v}_x$, $\vec{v}_y$, and $\vec{v}_z$ are the edges from that vertex to the three other vertices that vertex shares an edge with.

For the unit cube, you'd have $\vec{v}_0 = (0, 0, 0)$, $\vec{v}_x = (1, 0, 0)$, $\vec{v}_y = (0, 1, 0)$, and $\vec{v}_z = (0, 0, 1)$.

The cuboid has six faces. Each pair of opposing faces has the same area: $$\bbox{\begin{aligned} A_{x y} &= \left\lVert\vec{v}_x \times \vec{v}_y \right\rVert \\ A_{x z} &= \left\lVert\vec{v}_x \times \vec{v}_z \right\rVert \\ A_{y z} &= \left\lVert\vec{v}_y \times \vec{v}_z \right\rVert \\ \end{aligned}}$$ Note that if you have $\vec{a} = ( a_x , a_y , a_z )$ and $\vec{b} = ( b_x , b_y , b_z )$, $$\bbox{\vec{a} \times \vec{b} = ( a_y b_z - a_z b_y , a_z b_x - a_x b_z , a_x b_y - a_y b_x )}$$ and $$\bbox{\left\lVert\vec{a}\right\rVert = \sqrt{\vec{a} \cdot \vec{a}} = \sqrt{ a_x^2 + a_y^2 + a_z^2 }}$$

Generate two random scalars $0 \le i \le 1$ and $0 \le j \le 1$, and a random scalar $s$, $0 \le s \lt 2 A_{x y} + 2 A_{x z} + 2 A_{y z}$: $$\bbox{\begin{array}{l|c} \vec{p}(i, j) & s \\ \hline \vec{v}_0 + i \vec{v}_x + j \vec{v}_y & 0 \le s \lt A_{x y} \\ \vec{v}_0 + \vec{v}_z + i \vec{v}_x + j \vec{v}_y & A_{x y} \le s \lt 2 A_{x y} \\ \vec{v}_0 + i \vec{v}_x + v \vec{v}_z & 2 A_{x y} \le s \lt A_{x z} + 2 A_{x y} \\ \vec{v}_0 + \vec{v}_y + i \vec{v}_x + j \vec{v}_z & 2 A_{x y} + A_{x z} \le s \lt 2 A_{x z} + 2 A_{x y} \\ \vec{v}_0 + i \vec{v}_y + j \vec{v}_z & 2 A_{x y} + 2 A_{x z} \le s \lt A_{y z} + 2 A_{x z} + 2 A_{x y} \\ \vec{v}_0 + \vec{v}_x + i \vec{v}_y + j \vec{v}_z & 2 A_{x y} + 2 A_{x z} + A_{y z} \le s \lt 2 A_{y z} + 2 A_{x z} + 2 A_{x y} \\ \end{array}}$$ The uniform random point on the surface of the cuboid is then $\vec{p}(i, j)$, as shown in the left column corresponding to the range on the right side $s$ falls into.

The idea here is that we generate each point uniformly distributed within their face, and choose the face randomly weighted by their areas. This yields an uniform distribution on the face of the cuboid.

In Python:

from math import sqrt
"""
SPDX-License-Identifier: CC0-1.0
"""

def UniformCuboidSurface(generator, v0, vx, vy, vz):
    """Return uniform random points on the surface of cuboid,
       with vertices at (v0), (v0+vx), (v0+vy), (v0+vx+vy),
       (v0+vz), (v0+vx+vz), (v0+vy+vz), and (v0+vx+vy+vz).
       Uses generator.uniform(0.0, max) for uniform random reals."""
    v0, vx, vy, vz = Vec3d(v0), Vec3d(vx), Vec3d(vy), Vec3d(vz)
    Axy = vx.crossnorm(vy)
    Axz = vx.crossnorm(vz)
    Ayz = vy.crossnorm(vz)
    s = float(generator.uniform(0.0, 2.0*(Axy + Axz + Ayz)))
    i = float(generator.uniform(0.0, 1.0))
    j = float(generator.uniform(0.0, 1.0))
    limit = Axy
    if s < limit:
        return v0 + vx.scale(i) + vy.scale(j)
    limit += Axy
    if s < limit:
        return v0 + vz + vx.scale(i) + vy.scale(j)
    limit += Axz
    if s < limit:
        return v0 + vx.scale(i) + vz.scale(j)
    limit += Axz
    if s < limit:
        return v0 + vy + vx.scale(i) + vz.scale(j)
    limit += Ayz
    if s < limit:
        return v0 + vy.scale(i) + vz.scale(j)
    else:
        return v0 + vx + vy.scale(i) + vz.scale(j)


class Vec3d(tuple):
    """3D vector type with real components."""

    def __new__(class_, *args):
        try:
            if len(args) == 1 and isinstance(args[0], Vec3d):
                return args[0]
            if len(args) == 1 and len(args[0]) == 3:
                return tuple.__new__(class_, ( float(args[0][0]),
                                               float(args[0][1]),
                                               float(args[0][2]) ))
            if len(args) == 3:
                return tuple.__new__(class_, ( float(args[0]),
                                               float(args[1]),
                                               float(args[2]) ))
        except:
            pass
        raise ValueError("Invalid components for a 3D vector.")

    def __repr__(self):
        """String representation"""
        return "(%.6f, %.6f, %.6f)" % self

    def __eq__(self, other):
        """Exactly equal to"""
        return ( self[0] == float(other[0]) and
                 self[1] == float(other[1]) and
                 self[2] == float(other[2]) )

    def __ne__(self, other):
        """Not exactly equal to"""
        return ( self[0] != float(other[0]) or
                 self[1] != float(other[1]) or
                 self[2] != float(other[2]) )

    def __mul__(self, other):
        """Vector scaling or cross product"""
        if isinstance(other, (int, float)):
            other = float(other)
            return Vec3d(self[0]*other, self[1]*other, self[2]*other)
        if len(other) == 3:
            return self.cross(other)
        return NotImplemented

    def __rmul__(self, other):
        """Vector scaling or cross product"""
        if isinstance(other, (int, float)):
            other = float(other)
            return Vec3d(other*self[0], other*self[1], other*self[2] )
        if len(other) == 3:
            other = ( float(other[0]), float(other[1]), float(other[2]) )
            return Vec3d(other[1]*self[2] - other[2]*self[1],
                         other[2]*self[0] - other[0]*self[2],
                         other[0]*self[1] - other[1]*self[0])
        return NotImplemented

    def __div__(self, other):
        """Vector scaling or dot product"""
        if isinstance(other, (int, float)):
            other = float(other)
            return Vec3d(self[0]/other, self[1]/other, self[2]/other)
        if len(other) == 3:
            return ( self[0]*float(other[0]) + 
                     self[1]*float(other[1]) + 
                     self[2]*float(other[2]) )
        return NotImplemented

    def __add__(self, other):
        """Vector addition"""
        if len(other) == 3:
            return Vec3d( self[0] + float(other[0]),
                          self[1] + float(other[1]),
                          self[2] + float(other[2]) )
        return NotImplemented

    def __sub__(self, other):
        """Vector subtraction"""
        if len(other) == 3:
            return Vec3d( self[0] - float(other[0]),
                          self[1] - float(other[1]),
                          self[2] - float(other[2]) )
        return NotImplemented

    def cross(self, other):
        """Vector cross product"""
        other = ( float(other[0]), float(other[1]), float(other[2]) )
        return Vec3d( self[1]*other[2] - self[2]*other[1],
                      self[2]*other[0] - self[0]*other[2],
                      self[0]*other[1] - self[1]*other[0] )

    def crossnorm(self, other):
        """Euclidean norm of vector cross product"""
        other = ( float(other[0]), float(other[1]), float(other[2]) )
        return sqrt( (self[1]*other[2] - self[2]*other[1])**2 +
                     (self[2]*other[0] - self[0]*other[2])**2 +
                     (self[0]*other[1] - self[1]*other[0])**2 )

    def scale(self, factor):
        """Scale vector"""
        factor = float(factor)
        return Vec3d(self[0]*factor, self[1]*factor, self[2]*factor)

    def dot(self, other):
        """Vector dot product"""
        return ( self[0]*float(other[0]) + 
                 self[1]*float(other[1]) + 
                 self[2]*float(other[2]) )

    @property
    def x(self):
        """X coordinate"""
        return self[0]

    @property
    def y(self):
        """Y coordinate"""
        return self[1]

    @property
    def z(self):
        """Z coordinate"""
        return self[2]

Most of the code is the 3D real vector class, Vec3d, implementation. As an example, you could do e.g.

import random

rng = random.Random()
v0 = Vec3d(1.00, 2.00, 3.00)
vx = Vec3d(2.00, 4.00, 1.00)
vy = Vec3d(0.00,-1.00, 4.00)
vz = Vec3d(4.25,-2.00,-0.50)

for n in range(0, 1000):
    print("%.6f %.6f %.6f" % UniformCuboidSurface(rng, v0, vx, vy, vz))

to print 1000 uniform random points on the surface of a cuboid with vertices at $(1.00, 2.00, 3.00)$, $(3.00, 6.00, 4.00)$, $(1.00, 1.00, 7.00)$, $(3.00, 5.00, 8.00)$, $(5.25, 0.00, 2.50)$, $(7.25, 4.00, 3.50)$, $(5.25, -1.00, 6.50)$, and $(7.25, 3.00, 7.50)$.

---

Sampling points using a regular grid on each face is simple, if you don't care whether the grid is orthogonal (or if the sides are all perpendicular to each other, $\vec{v}_x \cdot \vec{v}_y = \vec{v}_x \cdot \vec{v}_z = \vec{v}_y \cdot \vec{v}_z = 0$). In that case, just use a regular grid in $i$, $j$ on each face (separately), noting that the true length of the $i$ axis is $\left\lVert\vec{v}_x\right\rVert$ on the first four faces, and $\left\lVert\vec{v}_y\right\rVert$ on the last two faces; and the true length of the $j$ axis is $\left\lVert\vec{v}_y\right\rVert$ on the first two faces, and $\left\lVert\vec{v}_z\right\rVert$ on the last four faces.

---

Sampling points on a regular grid in Cartesian coordinates $x$, $y$, $z$ is much more interesting.

First, we consider the internal unit coordinate system $(i, j, k)$ within the cube, $0 \le i, j, k \le 1$: $$\bbox{\vec{p}(i, j, k) = \vec{v}_0 + i \vec{v}_x + j \vec{v}_y + k \vec{v}_z}$$ This means that point $\vec{p}(i, j, k) = (x, y, z)$ corresponds to coordinates $$\bbox{\begin{cases} d = \vec{v}_x \cdot ( \vec{v}_y \times \vec{v}_z ) \\ i = \frac{(\vec{p} - \vec{v}_0) \cdot ( \vec{v}_y \times \vec{v}_z )}{d} \\ j = \frac{(\vec{p} - \vec{v}_0) \cdot ( \vec{v}_z \times \vec{v}_x )}{d} \\ k = \frac{(\vec{p} - \vec{v}_0) \cdot ( \vec{v}_x \times \vec{v}_z )}{d} \\ \end{cases}}$$ Because the cuboid is a linear construct, we can use the above to solve the internal coordinates where any plane of our choosing (in the real 3D coordinates) intersects the cuboid.

If we use $\vec{n}$ for the normal of the plane we are interested in, and $d$ is the signed distance from origin to that plane in units of $\left\lVert\vec{n}\right\rVert$ (so point $\vec{p}$ is on the plane if $\vec{p} \cdot \vec{n} - d = 0$), we can write the equation for that plane in the cuboid internal coordinates $(i, j, k)$ as $$\bbox{ (\vec{n}\cdot\vec{v}_x) i + (\vec{n}\cdot\vec{v}_y) j + (\vec{n}\cdot\vec{v}_z) k + (\vec{n} \cdot \vec{v}_0 - d) = 0}$$

This plane intersects the twelve edges of our cuboid at $$\begin{array}{lll} \displaystyle i = \frac{d - \vec{n} \cdot \vec{v}_0 }{\vec{n} \cdot \vec{v}_x } , & j = 0 , & k = 0 \\ \displaystyle i = \frac{d - \vec{n} \cdot (\vec{v}_0 + \vec{v}_z) }{\vec{n} \cdot \vec{v}_x } , & j = 0 , & k = 1 \\ \displaystyle i = \frac{d - \vec{n} \cdot (\vec{v}_0 + \vec{v}_y) }{\vec{n} \cdot \vec{v}_x } , & j = 1 , & k = 0 \\ \displaystyle i = \frac{d - \vec{n} \cdot (\vec{v}_0 + \vec{v}_y + \vec{v}_z) }{\vec{n} \cdot \vec{v}_x } , & j = 1 , & k = 1 \\ \displaystyle j = \frac{d - \vec{n} \cdot \vec{v}_0 }{\vec{n} \cdot \vec{v}_y } , & i = 0 , & k = 0 \\ \displaystyle j = \frac{d - \vec{n} \cdot (\vec{v}_0 + \vec{v}_z) }{\vec{n} \cdot \vec{v}_y } , & i = 0 , & k = 1 \\ \displaystyle j = \frac{d - \vec{n} \cdot (\vec{v}_0 + \vec{v}_x) }{\vec{n} \cdot \vec{v}_y } , & i = 1 , & k = 0 \\ \displaystyle j = \frac{d - \vec{n} \cdot (\vec{v}_0 + \vec{v}_x + \vec{v}_z) }{\vec{n} \cdot \vec{v}_y } , & i = 1 , & k = 1 \\ \displaystyle k = \frac{d - \vec{n} \cdot \vec{v}_0}{\vec{n} \cdot \vec{v}_z } , & i = 0 , & j = 0 \\ \displaystyle k = \frac{d - \vec{n} \cdot (\vec{v}_0 + \vec{v}_x)}{\vec{n} \cdot \vec{v}_z } , & i = 1 , & j = 1\\ \displaystyle k = \frac{d - \vec{n} \cdot (\vec{v}_0 + \vec{v}_y)}{\vec{n} \cdot \vec{v}_z } , & i = 0 , & j = 0\\ \displaystyle k = \frac{d - \vec{n} \cdot (\vec{v}_0 + \vec{v}_x + \vec{v}_y)}{\vec{n} \cdot \vec{v}_z } , & i = 1 , & j = 1\\ \end{array}$$ If you intend to draw a 2D projection of the 3D cube, the above allows you to map any plane perpendicular to your view (i.e., a line in 2D) to one to six line segments in internal cuboid coordinates on the surface of the cuboid. (Intersection between a cube an a plane is a point, line, triangle, quadrilateral (or tetragon), pentagon, or hexagon.)

Or, in other words, the above allows you to map any 2D line in your orthographic 2D projection of the 3D cube to zero to six line segments on the surface of the cube.

(If you use a perspective projection, the sampling along each line does not occur at fixed intervals, if you wish the sampling to be regular in the 2D coordinate system.)

---

If you wish to sample the surface using lines parallel to the $z$ axis, i.e. find the point on the cuboid surface that corresponds to 3D coordinates $(x, y, \text{any})$, using $\vec{v}_0 = ( x_0 , y_0 , z_0 )$, $\vec{v}_x = ( x_1 , y_1 , z_1 )$, $\vec{v}_y = ( x_2 , y_2 , z_2 )$, and $\vec{v}_z = ( x_3 , y_3 , z_3 )$, you need to find the internal coordinates on each face of the cuboid: $$\bbox{\begin{array}{lll} i = \frac{(x - x_0) y_2 - (y - y_0) x_2}{x_1 y_2 - x_2 y_1} , & j = \frac{(y - y_0) x_1 - (x - x_0) y_1}{x_1 y_2 - x_2 y_1} , & k = 0 \\ i = \frac{(x - x_0 - x_3) y_2 - (y - y_0 - y_3) x_2}{x_1 y_2 - x_2 y_1 } , & j = \frac{(y - y_0 - y_3) x_1 - (x - x_0 - x_3) y_1}{x_1 y_2 - x_2 y_1 } , & k = 1 \\ i = \frac{(x - x_0) y_3 - (y - y_0) x_3}{x_1 y_3 - x_3 y_1} , & j = 0 , & k = \frac{(y - y_0) x_1 - (x - x_0) y_1}{x_1 y_3 - x_3 y_1} \\ i = \frac{(x - x_0 - x_2) y_3 - (y - y_0 - y_2) x_3}{x_1 y_3 - x_3 y_1 } , & j = 1 , & k = \frac{(y - y_0 - y_2) x_1 - (x - x_0 - x_2) y_1}{x_1 y_3 - x_3 y_1 } \\ i = 0 , & j = \frac{(x - x_0) y_3 - (y - y_0) x_3}{x_2 y_3 - x_3 y_2} , & k = \frac{(y - y_0) x_2 - (x - x_0) y_2}{x_2 y_3 - x_3 y_2} \\ i = 1 , & j = \frac{(x - x_0 - x_1) y_3 - (y - y_0 - y_1) x_3}{x_2 y_3 - x_3 y_2} , & k = \frac{(y - y_0 - y_1) x_2 - (x - x_0 - x_1) y_2}{x_2 y_3 - x_3 y_2} \\ \end{array}}$$ Of those that have $0 \le i \le 1$, $0 \le j \le 1$, $0 \le k \le 1$, compute their $$\bbox{ z = z_0 + i z_1 + j z_2 + k z_3 }$$ and pick the one that has the smallest $z$ if smaller $z$ coordinates are in front of larger $z$ coordinates; the largest $z$ otherwise.

This is especially nice if you want to do sampling with regular real $x$ and $y$ coordinates. The above equations simplify to form $$\bbox{\begin{cases} i = C_0 + C_1 x + C_2 y \\ j = C_3 + C_4 x + C_5 y \\ k = 0 \end{cases}} , \quad \bbox{\begin{cases} i = C_6 + C_7 x + C_8 y \\ j = C_9 + C_{10} x + C_{11} y \\ k = 1 \end{cases}}$$ et cetera for each cuboid (36 scalar constants for each cuboid, in addition to the four vectors defining the cuboid), making this a viable option in a computer program.