Generate camera trajectory with cameras' lookat following shape

Question

I am trying with no avail to create a camera path that follows a lemniscate shape as shown in the figure.

I do it in 3D but just set the vertical dimension (y) to a fixed number so really is 2D. the x and z coordinates are generated with a closed form function. I calculated the derivative of the function and that should give me the normal at each point. With the normal and a world UP axis I estimate the lookat direction and build a lookat matrix. However, this is not working for me. The camera does not turn properly at the lemniscate curves. Is there any other way to do this?

I have the 3D points of the lemniscate and I know world coordinates are X-Right, Y-Down and Z-Forward.

This is the code:

    # Parameters
    a = 0.2  # Scale of the lemniscate

    # Create theta values
    theta1 = np.linspace(0.5*np.pi, 0.75*np.pi, 25)
    theta2 = np.linspace(0.75*np.pi, 1.25*np.pi, 80)
    theta3 = np.linspace(1.25*np.pi, 1.75*np.pi, 45)
    theta4 = np.linspace(1.75*np.pi, 2.25*np.pi, 80)
    theta5 = np.linspace(2.25*np.pi, 2.5*np.pi, 25)
    theta = np.concatenate((theta1, theta2, theta3, theta4, theta5))

    # Convert to Cartesian coordinates
    x = a * np.cos(theta) / (np.sin(theta)**2 + 1)
    z = a * np.cos(theta) * np.sin(theta) / (np.sin(theta)**2 + 1)
    y = a * 0.2 * np.cos(4*theta)
    Cs = np.stack((x, y, z)).T

    dx_dtheta = -a*(np.sin(theta) * (np.sin(theta) ** 2 + 2*np.cos(theta)**2 + 1)) / (np.sin(theta)**2 + 1)**2
    dz_dtheta = -a*(np.sin(theta)**4 + np.sin(theta)**2 + (np.sin(theta)**2-1)*np.cos(theta)**2) / (np.sin(theta)**2 + 1)**2
    dy_dtheta = - 4 * a * 0.2 * np.sin(4*theta)  # If I use this, it gets worse

    normal = np.stack((dx_dtheta, np.zeros_like(x), dz_dtheta)).T
    normal = normal / np.linalg.norm(normal, axis=1, keepdims=True)

    UP = np.array([0, -1, 0])[None, :]
    lookat = np.cross(normal, UP)
    lookat = lookat / np.linalg.norm(lookat, axis=1, keepdims=True)
    right = normal
    up = np.cross(right, lookat)
    up /= np.linalg.norm(up, axis=1, keepdims=True)

    Rs = np.concatenate((right, up, lookat), axis=1).reshape(-1, 3, 3).transpose(0, 2, 1)

Answer 1

In practice if you have a continuously twice differentiable regular curve $p(t) = (x(t), y(t), z(t))$ you can compute its velocity as $v(t) = \dot{p}(t) = \frac{dp}{dt}(t)$ and its acceleration as $a(t) = \dot{v}(t) = \ddot{p}(t)$ and use those to form an orthonormal basis (this wouod fail if $v(t)=0$ or $a(t)=0$). The velocity $v$ is tangent to the curve, but $a$ is not necessarily orthogonal to it, so if we wish to construct an orthonormal basis we need to orthonormalize those, e.g., by using Gram-Schmidt: \begin{equation} f_1(t) = \frac{v(t)}{\|v(t)\|_2}, \quad a_{\perp v}(t) = a(t) - (f_1(t)\cdot a(t))f_1(t), \quad f_2(t) = \frac{a(t)}{\|a(t)\|_2}. \end{equation}

We can compute the third vector by using the cross product $f_3(t) = f_1(t)\times f_2(t)$. Now $f_1(t), f_2(t), f_3(t)$ form an orthonormal basis, and you can use the following for your camera matrix $C$ and view matrix $V$: \begin{align} C &= \begin{bmatrix} | & | & | & | \\ f_1(t) & f_2(t) & f_3(t) & p(t) \\ | & | & | & | \\ 0 & 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} F(t) & p(t) \\ 0^T & 1 \end{bmatrix}, \\ V &= C^{-1} = \begin{bmatrix} (F(t))^T & -(F(t))^Tp(t) \\ 0^T & 1 \end{bmatrix}. \end{align}

Looking at your code the specific curve you have is parwmetrized as follows: \begin{align} x(t) = \frac{a \cos(t)}{1+\sin^2(t)}, \quad y(t) = \frac{a \cos(t)\sin(t)}{1+\sin^2(t)},\quad z(t) = b \cos(c t). \end{align}

The first derivatives (the components of the velocity) can be computed in the following way (I derived this without pen and paper so you should double check it): \begin{align} \dot{x}(t) &= -\frac{a\sin(t)}{1+\sin^2(t)}-\frac{2a\cos^2(t)\sin(t)}{(1+\sin^2(t))^2}, \\ \dot{y}(t) &= \frac{a \cos(2t)}{1+\sin^2(t)} - \frac{a\sin^2(2t)}{2(1+\sin^2(t))^2}, \\ \dot{z}(t) &= -bc\sin(ct). \end{align}

The second derivatives (components of the acceleration) are too tedious to compute so I suggest using wolfram alpha or some symbolic solver to do that for you.