Tuesday, October 18, 2016

A Note on Projective Splines

In this post, I'll make a seemingly in-obvious remark on the great work of Projective Splines by Thomas Lewiner and Marcos Craizer. Special thanks to Marcos Craizer for the discussion. The paper of interest is Projective splines and estimators for planar curves. The goal of the paper is an invariant spiral construction in the canonical projective frame. The setting looks like:
The idea is to use \(3\) points augmented by the tangents to generate new points and exploit the cross ratio in estimation of two projective invariants: projective length and projective curvature. The spline construction is basically done by solving a linear system in Section 3c. There, the last pair of equations require \((P_{x,3},P_{y,3})\), which are not explicitly given in the paper. Here I'll provide the expression for that: Following Section 3a, knowing \(P\), the standard spiral, the spiral corresponding to the newly generated point \(x_{i,3}\) could be written as \(\mathbf{P}_3(\sigma_i) = \mathbf{P}(\sigma_i) + t\mathbf{W}\). Here $$W(\sigma)=a e^{a\sigma}$$ and \(t\) is given in Section 3a. Also, \(a=\frac{3}{2}+\frac{i\sqrt{3}}{2}\). Here the complex notation is used to jointly represent coordinates \((x,y)\). If we plug \(a\) into \(W\): $$ \begin{align*} W(\sigma)&=a e^{a\sigma}\\ &=(\frac{3}{2}+\frac{i\sqrt{3}}{2}) e^{3\sigma/2} e^{\sqrt{3}i\sigma/2} \\ &=(\frac{3}{2}e^{3\sigma/2}+\frac{i\sqrt{3}}{2}e^{3\sigma/2}) e^{\sqrt{3}i\sigma/2} \\ &=\bigg(\frac{3}{2}e^{3\sigma/2}+\frac{i\sqrt{3}}{2}e^{3\sigma/2}\bigg) \bigg(cos(\frac{\sqrt{3}\sigma}{2})+i sin(\frac{\sqrt{3}\sigma}{2})\bigg) \text{ (via Euler's Formula)}\\ \text{fast forwarding...} & \\ &=\frac{\sqrt{3}}{2} e^{3\sigma/2} \bigg( \frac{\sqrt{3}}{2} cos(\frac{\sqrt{3}\sigma}{2}) - sin(\frac{\sqrt{3}\sigma}{2}) \bigg) + i \frac{\sqrt{3}}{2} e^{3\sigma/2} \bigg( \frac{\sqrt{3}}{2} sin(\frac{\sqrt{3}\sigma}{2}) + cos(\frac{\sqrt{3}\sigma}{2}) \bigg) \end{align*} $$ We could now define the vector \(\mathbf{W}\) to be: $$ \mathbf{W} = (W_x, W_y, W_z)^T = \begin{bmatrix} \frac{\sqrt{3}}{2} e^{3\sigma/2} \bigg( \frac{\sqrt{3}}{2} cos(\frac{\sqrt{3}\sigma}{2}) - sin(\frac{\sqrt{3}\sigma}{2}) \bigg) \\ \frac{\sqrt{3}}{2} e^{3\sigma/2} \bigg( \frac{\sqrt{3}}{2} sin(\frac{\sqrt{3}\sigma}{2}) + cos(\frac{\sqrt{3}\sigma}{2}) \bigg) \\ 1 \end{bmatrix} $$ Note that the linear system uses only two components of \(\mathbf{P}_3\), letting us to set \(W_z=w=1\). Now, let's get back to Section 2b to define: $$ \mathbf{P} = (P_x, P_y, P_z) = \begin{bmatrix} e^{\sigma/2}cos(\frac{\sqrt{3}\sigma}{2})\\ e^{\sigma/2}sin(\frac{\sqrt{3}\sigma}{2})\\ e^{-\sigma} \end{bmatrix} $$ being the components of standard spiral. To obtain the projective curve at w=1, we have to multiply all components by \(e^{\sigma}\). Plugging in \(\mathbf{P}_3(\sigma_i) = \mathbf{P}(\sigma_i) + t\mathbf{W}\) results in: $$ \mathbf{P}_3 = (P_{x,3}, P_{y,3}, P_{z,3})^T = \begin{bmatrix} e^{3\sigma/2}cos(\frac{\sqrt{3}\sigma}{2}) + t \frac{\sqrt{3}}{2} e^{3\sigma/2} \bigg( \frac{\sqrt{3}}{2} cos(\frac{\sqrt{3}\sigma}{2}) - sin(\frac{\sqrt{3}\sigma}{2}) \bigg) \\ e^{3\sigma/2}sin(\frac{\sqrt{3}\sigma}{2}) + t\frac{\sqrt{3}}{2} e^{3\sigma/2} \bigg( \frac{\sqrt{3}}{2} sin(\frac{\sqrt{3}\sigma}{2}) + cos(\frac{\sqrt{3}\sigma}{2}) \bigg) \\ t \end{bmatrix} $$ Here is the MATLAB code:
% requires P (the standard spiral), t (step), sigma (projective length)
% computes the canonical coordinate P_3
function [P3] = compute_p3(P, t, sigma)
s32 = sqrt(3)/2;
coeff = s32 * exp(1.5*sigma);
W = [
    coeff .* (s32.*cos(s32.*sigma) - sin(s32.*sigma));
    coeff .* (s32.*sin(s32.*sigma) + cos(s32.*sigma));
P3 = P + t.*W;

Sunday, October 16, 2016

A Better Approach to Plane Intersection

One interesting problem, arises in many geometric algorithms is the intersection of two planes. For sure, this problem has a well studied textbook answer, which I don't find particularly elegant. This is because, the solution is depends on picking an arbitrary point, which lacks geometric intuition. Ideally, we'd like this point to have some meaning, such as being close to the planes, or the line or etc. For that, in this post, I'd like to remind you of a solution by John Krumm, which remains to be unnoticed by many. While the plane-plane intersection problem is well studied, for our purposes, we would like to summarize a different, more effective approach. The textbook solution to the problem depends on picking an arbitrary point, which might lack geometric intuition. Ideally, we'd like this point to have some meaning, such as being close to the points, planes, or the line or etc. Let \(\mathbf{p}=\{p_x,p_y,p_z\}\) and \(\mathbf{n}=\{n_x,n_y,n_z\}\) compose the plane \(\mathbf{P}=\{\mathbf{n},\mathbf{p}\}\). Let there be two planes \(P_1\) and \(P_2\), for which we'd like to compute the intersecting line \(\mathbf{l}\). Let's first review the standard solution. It's trivial to compute the direction as the cross-product: \begin{equation} \mathbf{l}_d=\mathbf{n}_1 \times \mathbf{n}_2 \end{equation} We then pick an arbitrary point (origin, or the midpoint) \(\mathbf{p}_0\), plug it into both the equations and solve for the unknowns. At this stage there is no meaning of \(\mathbf{p}_0\). In the following, we'll pick this point to be the midpoint of the input points. Basically, we desire that the resulting point \(\mathbf{p}\) is as close to the chosen point \(\mathbf{p}_0\) as possible. Thus, we could write a distance : \begin{equation} \lVert \mathbf{p}-\mathbf{p}_0 \rVert = (p_x-p_{0x})^2 + (p_y-p_{0y})^2 + (p_z-p_{0z})^2 \end{equation} Incorporating the other points in the similar fashion, and writing this constraint using Lagrange multipliers into an objective function results in: \begin{equation} J=\lVert \mathbf{p}-\mathbf{p}_0 \rVert+\lambda(\mathbf{p}-\mathbf{p}_1)^2 + \mu(\mathbf{p}-\mathbf{p}_2)^2 \end{equation} Using the standard Lagrange framework (omitting the details), one establishes a nice matrix, in the form: \begin{equation} M= \left[ {\begin{array}{ccccc} 2 & 0 & 0 & n_{1x} & n_{2x}\\ 0 & 2 & 0 & n_{1y} & n_{2y}\\ 0 & 0 & 2 & n_{1z} & n_{2z} \\ n_{1x} & n_{1y} & n_{1z} & 0 & 0\\n_{2x} & n_{2y} & n_{2z} & 0 & 0 \end{array} } \right] \end{equation} This matrix can now be used in a system of linear equations: \begin{equation} \mathbf{M}\left[ {\begin{array}{c} p_x \\ p_y \\ p_z \\ \lambda \\ \mu \end{array} } \right] = \left[ {\begin{array}{c} 2p_{0x} \\ 2p_{0y} \\ 2p_{0z} \\ \mathbf{p}_1 \cdot \mathbf{n}_1 \\ \mathbf{p}_2 \cdot \mathbf{n}_2 \end{array} } \right] \end{equation} to solve for the unknown point, \(\mathbf{p}\), as well as the Lagrange multipliers, \(\{\lambda, \mu\}\). While the multipliers are not of particular interest they would be interesting for understanding configuration of points, or for different parametrizations. Here is the MATLAB script:
% Robust plane to plane intersection
% Plane is parameterized by a point and a normal.
function [p, n] = intersect_planes(p1, n1, p2, n2, p0)

M = [2 0 0 n1(1) n2(1)
     0 2 0 n1(2) n2(2)
     0 0 2 n1(3) n2(3)
     n1(1) n1(2) n1(3) 0 0
     n2(1) n2(2) n2(3) 0 0

b4 = p1(1).*n1(1) + p1(2).*n1(2) + p1(3).*n1(3);
b5 = p2(1).*n2(1) + p2(2).*n2(2) + p2(3).*n2(3);
b = [2*p0(1) ; 2*p0(2) ; 2*p0(3); b4 ; b5];

x = M\b;
p = x(1:3);
n = cross(n1, n2);