Transform a direction from world space to local space?

Question

I've been reading this pdf about vector transformations and I don't quite understand how to implement it in a computer program. On page 10, it shows you how to transform a vector from one coordinate system to another, exactly how would this look on paper? Say I had two 3d vectors that represented directions and I wanted to get the relative direction of $vector A$ to $vector B$. So if $A = (0,0,0)$ and $B = (0,1,0)$, $C$ would also be $(0,1,0)$. If $A = (0,1,0)$, $C$ would then be $(0,0,-1)$. Would I make matrices of $A$ and $B$ and just multiply them together? Other places have talked about getting the dot product? What would the $i$ and $i'$ vectors be in this case? Thank you.

I'm pretty sure this is all the info I need but in case there is an obviously better way of doing it I will explain my project. It's a simple voxel 3D raycaster. I'm actually using Unity and rendering to a texture using a compute shader.

The way I am calculating the rays is I am centering the pixel coordinates as if they were 3D space coordinates, moving them a little forward and getting the directions from 0,0,0 to each of those points (and shrinking the dimensions for field-of-view). Then, I transform those directions relative to the player's transform every frame to get the ray directions for the GPU. Obviously, this also needs to be calculated on the GPU so I can't use the handy Unity method I used just to test if it would work (Transform.TransformDirection). So I guess I could also use a relative point transform function too and just send those points and the player's transform to the GPU.

Answer 1

So this page ended up helping me. I simply send the positions of the pixels as if they were a screen slightly in front of the player at 0,0,0 and center them and whatnot, then in the compute shader I transform those points to world space relative to the player's rotation for each ray like this:

float3 transformDirectionFromPoint(float3 p) {

  float3 u1 = p.x * playerWorldRight;
  float3 u2 = p.y * playerWorldUp;
  float3 u3 = p.z * playerWorldForward;

  return u1 + u2 + u3; // the direction to that point
}

And then I use that direction to cast a ray out from the player's position. I found with higher resolutions I needed to decrease the field of view.

Answer 2

Voxel tracing differs from ordinary ray tracing in that you want to use a coordinate system where the voxels are defined for each unit cell, i.e. for unit (integer) coordinates.

You can use one vector, say $\vec{o}$ for observer (eye) location (in the voxel coordinates), and an unit quaternion $\mathbf{q} = q_r + q_x \mathbf{i} + q_y \mathbf{j} + q_z \mathbf{k} = (q_x, q_y, q_z; q_r)$ to describe the observer orientation; it corresponds to rotation matrix $\mathbf{R}$ via $$\mathbf{R} = \left[ \begin{matrix} 1-2( q_y^2 + q_z^2 ) & 2( q_x q_y- q_z q_r) & 2( q_x q_z+ q_y q_r) \\ 2( q_x q_y+ q_z q_r) & 1-2( q_x^2 + q_z^2 ) & 2( q_y q_z- q_x q_r) \\ 2( q_x q_z- q_y q_r) & 2( q_y q_z+ q_x q_r) & 1-2( q_x^2 + q_y^2 ) \\ \end{matrix} \right] = \left[ \begin{matrix} u_x & v_x & w_x \\ u_y & v_y & w_y \\ u_z & v_z & w_z \\ \end{matrix} \right] = \left[ \begin{matrix} \hat{u} & \hat{v} & \hat{w} \end{matrix} \right]$$ where $\hat{w}$ is the "forward" unit vector in voxel coordinates, $\hat{v}$ is "up", and $\hat{u}$ is "right". For the identity unit quaternion $\mathbf{q} = 1 = (0, 0, 0; 1)$, $\hat{u} = (1, 0, 0)$, $\hat{v} = (0, 1, 0)$, and $\hat{w} = (0, 0, 1)$.

The reason for using a quaternion (or a bivector, which yields the same math) is that the orientation can be safely normalized by dividing each of the four components by $\sqrt{q_r^2 + q_x^2 + q_y^2 + q_z^2}$, because for unit quaternions, $q_r^2 + q_x^2 + q_y^2 + q_z^2 = 1$. To apply a rotation $\mathbf{t}$ to the orientation, you use Hamilton product: $$\left\lbrace\begin{aligned} q_r^\prime &= t_r q_r - t_x q_x - t_y q_y - t_z q_z \\ q_x^\prime &= t_r q_x + t_x q_r + t_y q_z - t_z q_y \\ q_y^\prime &= t_r q_y - t_x q_z + t_y q_r + t_z q_x \\ q_z^\prime &= t_r q_z + t_x q_y - t_y q_x + t_z q_r \\ \end{aligned} \right.$$ If the projection plane (view in voxel coordinates) is $2 \chi$ wide and $2 \gamma$ tall, at distance $d$ from the observer, each ray $(i, j)$ ($0 \le i \le i_\max$, $0 \le j \le j_\max$) has unnormalized ray direction vector $\vec{s}$, $$\vec{s} = d \hat{w} + \left(\frac{i}{i_\max} - \frac{1}{2}\right) \chi \hat{u} + \left(\frac{j}{j_\max} - \frac{1}{2}\right) \gamma \hat{v}$$ and starts at $\vec{p}$, $$\vec{p} = \vec{o} + \vec{s}$$ The interesting thing is that if you calculate the unit ray vector $\hat{r}$, $$\hat{r} = \frac{\vec{s}}{\lVert\vec{s}\rVert} = \left[\begin{matrix} r_x \\ r_y \\ r_z \end{matrix} \right]$$ its reciprocal components tell you the interval at which the ray steps one voxel unit along each axis, $$\left\lbrace \begin{aligned} L_x &= \frac{1}{\lvert r_x \rvert}, \\ L_y &= \frac{1}{\lvert r_y \rvert}, \\ L_z &= \frac{1}{\lvert r_z \rvert}, \\ \end{aligned} \right. \quad \begin{aligned} S_x &= \operatorname{sign}(r_x) = \frac{r_x}{\lvert r_x \rvert} \\ S_y &= \operatorname{sign}(r_y) = \frac{r_y}{\lvert r_y \rvert} \\ S_z &= \operatorname{sign}(r_z) = \frac{r_z}{\lvert r_z \rvert} \\ \end{aligned}$$ where $S_x$, $S_y$, $S_z$ are $-1$, $0$, or $+1$ depending on the direction of the ray, and division by zero yields either $\infty$ or a value larger than the voxel map size.

The initial distance $d_x$ where the ray intersects a voxel $x$ face, is $$d_x = \begin{cases} ( p_x - \lfloor p_x \rfloor ) L_x, & S_x = -1 \\ \infty, & S_x = 0 \\ ( \lceil p_x \rceil - p_x ) L_x, & S_x = +1 \\ \end{cases}$$ $d_x \ge 0$; and similarly for the $y$ and $z$ axes.

Instead of progressing the ray at fixed length steps, you examine where the ray next intersects a voxel face. You do this by choosing the smallest of $d_x$, $d_y$, and $d_z$.

If $d_x$ is the smallest of the three, the ray length $d$ is set to $d_x$, then $d_x$ increased by $L_x$ (corresponding to the distance where the ray intersects the next $x$ face), and the ray intersects a voxel $x$ face at $\vec{p} + d \hat{r}$. Similarly for the other two coordinates.

You may need to handle the rare cases $d_x = d_y$, $d_x = d_z$, etc. where the ray intersects a voxel edge, and $d_x = d_y = d_z$ where the ray intersects a voxel vertex, separately.

The voxel is identified by the integer part of coordinates of $\vec{p} + d \hat{r}$. The fractional parts correspond to the fractional coordinates within the voxel, remembering that the above method only examines voxel faces; at least one of the coordinates is always an integer. This allows you to use textures for each voxel face, if you want.

---

Consider the following 2D scenario for a single ray starting at lower left, progressing up right: The distance along the ray between the blue tabs is $L_x$, and between the red tabs is $L_y$. We walk in increasing distance.

If we call vertical/blue faces $x$ (because they occur at integer $x$ coordinates), and horizontal/red faces $y$, the order in which the ray encounters these faces is $x$, $y$, $x$, $y$, $x$, $x$, $y$, $x$, $x$, $y$, $x$, $y$, $x$, $x$, $y$, and so on.

---

Nothing beats an example, I guess. Here is a horrible-crude-simple tracer in Python 3:

# SPDX-License-Identifier: CC0-1.0
# -*- coding: utf-8 -*-
from math import inf as _inf, sqrt as _sqrt, copysign as _copysign, floor as _floor

class Vector(tuple):
    """Immutable, compact 3D vector"""

    def __new__(cls, *args):
        """Create a new 3D vector

           Usage:
             v = Vector(x, y, z)
             v = Vector((x, y, z))
             v = Vector([x, y, z])
             v = Vector(v)
        """
        if len(args) == 1:
            if isinstance(args[0], (tuple, list)):
                args = args[0]
            else:
                raise TypeError("Cannot construct a Vector from a %s." % str(type(args[0])))
        if len(args) != 3:
            raise ValueError("Vector() needs three components; %d given." % len(args))

        return tuple.__new__(cls, (float(args[0]), float(args[1]), float(args[2])))

    def __init__(self, *args):
        """Vectors are immutable"""
        pass

    def __str__(self):
        return "(%.6f, %.6f, %.6f)" % self

    @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]

    @property
    def normsqr(self):
        """Squared norm, i.e. Euclidean length squared"""
        return self[0]*self[0] + self[1]*self[1] + self[2]*self[2]

    @property
    def length(self):
        """Euclidean length"""
        return _sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])

    @property
    def normalized(self):
        """Vector normalized to unit length"""
        n = _sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])
        if n == 1:
            return self
        elif n > 0:
            return tuple.__new__(self.__class__, (self[0]/n, self[1]/n, self[2]/n))
        else:
            raise ValueError("Cannot normalize a zero vector to unit length!")

    def dot(self, other):
        """Vector dot product"""
        if not isinstance(other, Vector):
            other = Vector(other)
        return self[0]*other[0] + self[1]*other[1] + self[2]*other[2]

    def cross(self, other):
        """Vector cross product"""
        if not isinstance(other, Vector):
            other = Vector(other)
        return tuple.__new__(self.__class__, ( 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 proj(self, other):
        """Vector projection to another vector"""
        if not isinstance(other, Vector):
            other = Vector(other)
        n = ( self[0]*other[0] +
              self[1]*other[1] +
              self[2]*other[2] ) / ( other[0]*other[0] +
                                     other[1]*other[1] +
                                     other[2]*other[2] )
        return tuple.__new__(self.__class__, ( n*other[0], n*other[1], n*other[2] ))

    def perp(self, other):
        """The part of current vector perpendicular to the other vector"""
        unit = Vector(other).normalized
        return self - unit*self.dot(unit)

    def __add__(self, other):
        """Vector addition"""
        if isinstance(other, (tuple, list)) and len(other) == 3:
            return tuple.__new__(self.__class__, ( self[0] + float(other[0]),
                                                   self[1] + float(other[1]),
                                                   self[2] + float(other[2]) ))
        else:
            return NotImplemented

    def __iadd__(self, other):
        """Vector addition"""
        if isinstance(other, (tuple, list)) and len(other) == 3:
            return tuple.__new__(self.__class__, ( self[0] + float(other[0]),
                                                   self[1] + float(other[1]),
                                                   self[2] + float(other[2]) ))
        else:
            return NotImplemented

    def __radd__(self, other):
        """Vector addition"""
        if isinstance(other, (tuple, list)) and len(other) == 3:
            return tuple.__new__(self.__class__, ( float(other[0]) + self[0],
                                                   float(other[1]) + self[1],
                                                   float(other[2]) + self[2] ))
        else:
            return NotImplemented

    def __sub__(self, other):
        """Vector subtraction"""
        if isinstance(other, (tuple, list)) and len(other) == 3:
            return tuple.__new__(self.__class__, ( self[0] - float(other[0]),
                                                   self[1] - float(other[1]),
                                                   self[2] - float(other[2]) ))
        else:
            return NotImplemented

    def __isub__(self, other):
        """Vector subtraction"""
        if isinstance(other, (tuple, list)) and len(other) == 3:
            return tuple.__new__(self.__class__, ( self[0] - float(other[0]),
                                                   self[1] - float(other[1]),
                                                   self[2] - float(other[2]) ))
        else:
            return NotImplemented

    def __rsub__(self, other):
        """Vector subtraction"""
        if isinstance(other, (tuple, list)) and len(other) == 3:
            return tuple.__new__(self.__class__, ( float(other[0]) - self[0],
                                                   float(other[1]) - self[1],
                                                   float(other[2]) - self[2] ))
        else:
            return NotImplemented

    def __mul__(self, other):
        """Vector-scalar multiplication"""
        if isinstance(other, (int, float)):
            return tuple.__new__(self.__class__, ( self[0]*other,
                                                   self[1]*other,
                                                   self[2]*other ))
        else:
            return NotImplemented

    def __imul__(self, other):
        """Vector-scalar multiplication"""
        if isinstance(other, (int, float)):
            return tuple.__new__(self.__class__, ( self[0]*other,
                                                   self[1]*other,
                                                   self[2]*other ))
        else:
            return NotImplemented

    def __rmul__(self, other):
        """Scalar-vector multiplication"""
        if isinstance(other, (int, float)):
            return tuple.__new__(self.__class__, ( other*self[0],
                                                   other*self[1],
                                                   other*self[2] ))
        else:
            return NotImplemented

    def __truediv__(self, other):
        """Vector-scalar division"""
        if isinstance(other, (int, float)):
            return tuple.__new__(self.__class__, ( self[0]/other,
                                                   self[1]/other,
                                                   self[2]/other ))
        else:
            return NotImplemented

    def __itruediv__(self, other):
        """Vector-scalar division"""
        if isinstance(other, (int, float)):
            return tuple.__new__(self.__class__, ( self[0]/other,
                                                   self[1]/other,
                                                   self[2]/other ))
        else:
            return NotImplemented

    def __rtruediv__(self, other):
        """Division by vector not implemented"""
        return NotImplemented

    def __neg__(self):
        """Unary - (Vector negation)"""
        return tuple.__new__(self.__class__, ( -self[0], -self[1], -self[2] ))

    def __pos__(self):
        """Unary + (No effect)"""
        return self

    def __abs__(self):
        """abs() (Euclidean length)"""
        return _sqrt(self[0]*self[0] + self[1]*self[1] + self[2]*self[2])

    def __invert__(self):
        """Unary ~ (Vector scaled to reciprocal length)"""
        n = self[0]*self[0] + self[1]*self[1] + self[2]*self[2]
        return tuple.__new__(self.__class__, ( self[0]/n, self[1]/n, self[2]/n ))

    def __bool__(self):
        """Zero vectors are False, nonzero vectors True"""
        return (self[0]*self[0] + self[1]*self[1] + self[2]*self[2]) > 0.0

def _step_start_length(direction, start):
    one = abs(direction)

    try:
        len_x = one / abs(direction.x)
        sgn_x = int(_copysign(1.0, direction.x))
    except ZeroDivisionError:
        len_x = _inf
        sgn_x = 0

    try:
        len_y = one / abs(direction.y)
        sgn_y = int(_copysign(1.0, direction.y))
    except ZeroDivisionError:
        len_y = _inf
        sgn_y = 0

    try:
        len_z = one / abs(direction.z)
        sgn_z = int(_copysign(1.0, direction.z))
    except ZeroDivisionError:
        len_z = _inf
        sgn_z = 0

    if sgn_x > 0:
        dist_x = (_floor(start.x + 1) - start.x) * len_x
    elif sgn_x < 0:
        dist_x = (start.x - _floor(start.x)) * len_x
    else:
        dist_x = len_x

    if sgn_y > 0:
        dist_y = (_floor(start.y + 1) - start.y) * len_y
    elif sgn_y < 0:
        dist_y = (start.y - _floor(start.y)) * len_y
    else:
        dist_y = len_y

    if sgn_z > 0:
        dist_z = (_floor(start.z + 1) - start.z) * len_z
    elif sgn_z < 0:
        dist_z = (start.z - _floor(start.z)) * len_z
    else:
        dist_z = len_z

    return (sgn_x, sgn_y, sgn_z), Vector(dist_x, dist_y, dist_z), Vector(len_x, len_y, len_z)


class View:
    """2D view into a voxel lattice"""

    def __init__(self, **kwargs):
        self._voxel = None
        self._position = Vector(0, 0, -50)
        self._basis = (Vector(1,0,0), Vector(0,1,0), Vector(0,0,1))
        self._picture_distance = 10
        self._picture_halfsize = 8
        self._xsize = 320
        self._ysize = 160
        self._colors = [ b'\x00\x00\x00',   # Background #000000
                         b'\x33\x99\xff',   #    X faces #3399ff
                         b'\x33\xff\x99',   #    Y faces #33ff99
                         b'\x33\xcc\xcc',   #   XY edges #33cccc
                         b'\xff\x99\x33',   #    Z faces #ff9933
                         b'\xcc\x66\xcc',   #   XZ edges #cc66cc
                         b'\xcc\xcc\x66',   #   YZ edges #cccc66
                         b'\xcc\xcc\xcc',   #  XYZ vertx #cccccc
                       ]

    def renderPPM(self, output, voxel):
        output.write(b'P6\n%d %d 255\n' % (self._xsize, self._ysize))
        picx = (self._picture_halfsize / self._xsize) * self._basis[0]
        picy = (-self._picture_halfsize / self._xsize) * self._basis[1]
        pic0 = self._position + self._picture_distance * self._basis[2] - picx * (0.5 * self._xsize) - picy * (0.5 * self._ysize)

        for y in range(0, self._ysize):
            for x in range(0, self._xsize):
                start = pic0 + x*picx + y*picy
                delta = (start - self._position).normalized
                step, nextdist, length = _step_start_length(delta, start)
                dist = 0

                while True:

                    # Blank wall at distance 200
                    if dist >= 200:
                        face = 0
                        break

                    elif nextdist.x < nextdist.y and nextdist.x < nextdist.z:
                        dist = nextdist.x
                        nextdist += (length.x, 0, 0)
                        face = 1

                    elif nextdist.y < nextdist.x and nextdist.y < nextdist.z:
                        dist = nextdist.y
                        nextdist += (0, length.y, 0)
                        face = 2

                    elif nextdist.z < nextdist.x and nextdist.z < nextdist.y:
                        dist = nextdist.z
                        nextdist += (0, 0, length.z)
                        face = 4

                    elif nextdist.x == nextdist.y and nextdist.x < nextdist.z:
                        dist = nextdist.x
                        nextdist += (length.x, length.y, 0)
                        face = 1 + 2

                    elif nextdist.x == nextdist.z and nextdist.x < nextdist.y:
                        dist = nextdist.x
                        nextdist += (length.x, 0, length.z)
                        face = 1 + 4

                    elif nextdist.y == nextdist.z and nextdist.y < nextdist.x:
                        dist = nextdist.y
                        nextdist += (0, length.y, length.z)
                        face = 2 + 4

                    else:
                        dist = nextdist.x
                        nextdist += (length.x, length.y, length.z)
                        face = 1 + 2 + 4

                    if voxel(start + dist * delta, face):
                        break

                output.write(self._colors[face])

            print("Row %d of %d rendered" % (y + 1, self._ysize))


def voxel(point, face):

    # Voxel map has a spheroid-thingy of radius 6 at (10,4,3)
    if (point - (10,4,3)).normsqr <= 36:
        return True

    # and a bigger one at (-20,0,10)
    if (point - (-20,0,10)).normsqr <= 144:
        return True

    return False


if __name__ == '__main__':

    view = View()

    with open('out.ppm', 'wb') as out:
        view.renderPPM(out, voxel)

    print("Saved 'out.ppm'.")

Near the end, the voxel() function determines what is visible. It saves a NetPBM P6 image (Portable Pixmap format) out.ppm, which looks like: