Edited 06-29: updated how to compute the length of the direction vector. See the end of this post.

Intuitively, the Chinese stroke “㇐”U+31D0 has a horizontal direction, while “㇓”U+31D3 has a mostly top-right to bottom-left direction. How do we determine the “direction” of a stroke?

First we need to figure out what is a stroke and what is a “direction”. Since we are working in a 2D plane, we can instead use a method for any closed loop without self-intersection. Most frequently, these shapes are represented in piecewise parametric curves[1], which is also the case for font-shapped characters from free-type fonts. One possible definition of “direction” is the main axis of the smallest bounding ellipsis. Computing the smallest bounding ellipsis is hard to do analytically, but we can deal with another form of ellipsis much more easily: statistics has given us a tool for dealing with fitted ellipsis: principal component analysis (PCA).

In practice, PCA is always used with samples, but instead we can abuse it to compute the ellipsis fitted to a multivariate distribution with known PDF. In our case, we can define the “direction” of a shape LL to be the main axis of the ellipsis fitted to the uniform distribution over the internal of LL.

What do we need

First, the PDF of the distribution over the internal of some shape closed LL with internal int(L)int(L) and area A(L)A(L) is:

f(p)={0pint(L),1A(L)pint(L), f(\mathbf{p}) = \begin{cases} 0& \mathbf{p} \notin int(L),\\ 1 \over A(L) & \mathbf{p} \in int(L), \end{cases}

So first, we need to compute the area of LL. Another subtle thing that needs to be determine is the orientation of LL, because it’s not guaranteed that it has conter-clockwise orientation. Both of these can be computed in a single step by first assuming LL has the correct orientation, compute the area A(L)A(L) as-is, and then judge the orientation by the sign of A(L)A(L).

A(L)=int(L)dA=Lxdy A(L) = \iint_{int(L)} dA = \oint_L x dy

Now with the PDF fully defined, we need to compute the covariance matrix of this distribution. This in turn requires the following values:

  • Var(x)Var(x) and Var(y)Var(y)
  • Cov(x,y)=Cov(y,x)Cov(x, y) = Cov(y, x)

Since Var(x)=E(x2)E(x)2Var(x) = E(x^2) - E(x)^2 and Var(y)=E(y2)E(y)2Var(y) = E(y^2) - E(y)^2, so we need the first and second central moment of x and y, which are:

  • E(x)=1A(L)int(L)xdA=12A(L)Lx2dyE(x) = \frac{1}{A(L)} \iint_{int(L)} x dA = \frac{1}{2A(L)} \oint_L x^2 dy
  • E(y)=1A(L)int(L)ydA=12A(L)Ly2dxE(y) = \frac{1}{A(L)} \iint_{int(L)} y dA = - \frac{1}{2A(L)} \oint_L y^2 dx
  • E(x2)=1A(L)int(L)x2dA=13A(L)Lx3dyE(x^2) = \frac{1}{A(L)} \iint_{int(L)} x^2 dA = \frac{1}{3A(L)} \oint_L x^3 dy
  • E(y2)=1A(L)int(L)y2dA=13A(L)Ly3dxE(y^2) = \frac{1}{A(L)} \iint_{int(L)} y^2 dA = - \frac{1}{3A(L)} \oint_L y^3 dx

And for the covariance:

Cov(x,y)=1A(L)int(L)xydA=14A(L)Lx2ydyy2xdx Cov(x, y) = \frac{1}{A(L)} \iint_{int(L)} x y dA = \frac{1}{4A(L)} \oint_L x^2 y dy - y^2 x dx

Finally, compute the largest eigenvector & eigenvalue for the covariance matrix:

[Var(x)Cov(x,y)Cov(x,y)Var(y)] \begin{bmatrix} Var(x) & Cov(x, y) \\ Cov(x, y) & Var(y) \end{bmatrix}

Since this is a symmetric matrix, a real eigenvalue & eigenvector is guaranteed.

Computation

In practice, all curves are SVG paths. We first preprocess the paths to split them into non-overlapping components that do not self-intersect.

Then we can use the lyon_algorithms crate for sampling along a SVG path. Extra caution should be taken for the left-over segment at the end of the integration.

For computing the eigenvalue and eigenvector, we used two tricks for 2x2 matrices, one from 3b1b (YouTube), and another one from math stackexchange. For the following matrix:

A=[abcd] A = \begin{bmatrix} a & b \\ c & d \end{bmatrix}

Recall that the sum of eigenvalues λ1+λ2=tr(A)\lambda_1 + \lambda_2 = tr(A), and the product of eigenvalues λ1λ2=det(A)\lambda_1 \lambda_2 = det(A). So we can directly solve for the larger eigenvalue λ1=12tr(A)+14tr(A)2det(A)\lambda_1 = \frac{1}{2} tr(A) + \sqrt{\frac{1}{4} tr(A)^2 - det(A)}.

Then, we can compute the (unnormalized) eigenvector corresponding to λ1\lambda_1 to be:

(bλ1a) \big( \begin{matrix} b \\ \lambda_1 - a \end{matrix} \big)

BTW, since we are dealing with a symmetric real matrix (b=cb = c), it’s guaranteed that:

(ad)20    (a+d)24ad4ad4bc    14tr(A)2=14(a+d)2adbc=det(A) \begin{aligned} & (a - d)^2 \geq 0 \\ \implies & (a + d)^2 \geq 4 a d \geq 4 a d - 4 b c \\ \implies & \frac{1}{4} tr(A)^2 = \frac{1}{4} (a + d)^2 \geq ad - bc = det(A) \end{aligned}

Finally, we scale the eigenvector to have the length λ1\sqrt{\lambda_1}, because now it repersent the “standard derivation” along that axis.

Result

Title from the last post

The above figure shows the computed direction vector scaled by 10. Because the length of the original vector is the standard derivation along the principal axis, it’s almost always being covered by the stroke itself.

Also because of this, the length of the vector does not directly reflect the “outer dimension” of the stroke. The distribution of thickness along the axis will affect the length of the vector. One may argue this is not a good indicator of the dimension of a shape. In that case, a distribution along the boundary of the shape may be used instead of the internal of the shape.

But for me, this is good enough™, and does reflect the optical heaviness of the various part of a stroke, so I’m settling with this.


Edit 06/29: We noticed a problem with the aforementioned approach to compute the length of direction vectors. When the primary and secondary component has similar derivation, the direction of the primary component alone may not actually suggest the shape of the original distribution. An example is the dot in the figure above.

Instead, we choose to scale the direction vector based on the ratio between the derivations along primary and secondary component (square root of the two eigenvalues):

L=1λ2λ1 L = 1 - \sqrt{\lvert \frac{\lambda_2}{\lambda_1} \rvert}

This mapping has the property that if λ1λ2\lambda_1 \approx \lambda_2, then L=0L = 0, effectively suggests that the shape has no direction. The more eccentric the fitted ellipsis is, the closer LL is to 1. One can also just the use ellipsis eccentricity.

Numerically, λ10\lambda_1 \approx 0 indicates a small vector shape. A good fallback is to just return L=0L = 0.


- [1]: more concretely, SVG paths, or bezier curve segments with various orders.