Sutherland-Hodgman

We will implement the Sutherland-Hodgman polygon clipping algorithm ^[1], as it is pretty intuitive and not too hard to implement.

First, we will define, what a clipping plane is. This is basically the same concept as with the lines in the 2D case. There are different ways to define a plane. The first is with a normal $\mathbf{n}$ and a point $\mathbf{p}$. We want to check, if a point $\mathbf{x}$ is part of the plane. To do that, we compute the vector between both points $\mathbf{x} - \mathbf{p}$ (Note: With this and following definitions, you can switch the order of vectors, which will change some signs along the way. Just be careful to be consistent, otherwise it doesn't matter). The difference vector is part of the plane, so it has to be perpendicular to the plane normal.

$$\mathbf{n}\cdot(\mathbf{x} - \mathbf{p}) = 0$$

Rearranging a bit we will go from this vector equation to the scalar form:

$$\begin{align*} \mathbf{n}\cdot(\mathbf{x} - \mathbf{p}) &= 0\\ \mathbf{n}\cdot\mathbf{x} \underbrace{-\mathbf{n}\cdot\mathbf{p}}_{d} &= 0\\ \mathbf{n}\cdot\mathbf{x} +d &= 0 \end{align*}$$

If $\mathbf{n}$ is normalized, $d$ has the nice geometric meaning of being the signed distance of the origin to the plane.

So the function $\operatorname{d}(\mathbf{x}) = \mathbf{n}\cdot\mathbf{x} +d$ gives us the signed distance (scaled by the length of the normal) to the plane, with positive values meaning "$\mathbf{x}$ is in front" and negative values meaning "$\mathbf{x}$ is in the back".

The function is also linear (there are only multiplications, additions and now powers or anything like that), which means, we can write the following:

$$\operatorname{d}(\mathbf{x} + s\mathbf{y}) = \operatorname{d}(\mathbf{x}) + s\operatorname{d}(\mathbf{y})$$

For clipping polygons, we only need two operations:

1. Check, if a point lies in front or to the back of a plane
2. Intersect a line with a plane

The first part is already covered with the sign of $\operatorname{d}$. We will now use the linearity of $\operatorname{d}$ to solve the second point.

The line between two points $\mathbf{A}$ and $\mathbf{B}$ is defined by:

$$\begin{align*} \mathbf{p}(t) &= \mathbf{A} + t (\mathbf{B} - \mathbf{A}) \\ &= \mathbf{A} + t\mathbf{v} \\ t &= \in[0,1] \end{align*}$$

The intersection of this line with the plane will be a point, that has a distance of zero to the plane. We will just plug in the line and solve for the parameter.

$$\begin{align*} \operatorname{d}(\mathbf{p}(t)) &= 0 \\ \operatorname{d}(\mathbf{A} + t (\mathbf{B} - \mathbf{A})) &= 0\\ \operatorname{d}(\mathbf{A}) + t\operatorname{d}(\mathbf{B} - \mathbf{A}) &= 0\\ \operatorname{d}(\mathbf{A}) + t(\operatorname{d}(\mathbf{B}) - \operatorname{d}(\mathbf{A})) &= 0\\ \operatorname{d}(\mathbf{A}) + t(\operatorname{d}(\mathbf{B}) - \operatorname{d}(\mathbf{A})) &= 0\\ t(\operatorname{d}(\mathbf{B}) - \operatorname{d}(\mathbf{A})) &= -\operatorname{d}(\mathbf{A})\\ t &= -\frac{\operatorname{d}(\mathbf{A})}{\operatorname{d}(\mathbf{B}) - \operatorname{d}(\mathbf{A})}\\ t &= \frac{\operatorname{d}(\mathbf{A})}{\operatorname{d}(\mathbf{A}) - \operatorname{d}(\mathbf{B})}\\ \end{align*}$$

One last part that we will do, to also be prepared for later, is writing our vectors in a specific way. For now, you can just think about it as being a different notation. But doing it this way will actually save us rewriting thing in a later step. We will also use 3D vectors in general now, not only 2D. As you have seen from the formulas above, they don't depend on the dimension. And we can just set the z coordinate to 0 for 2D.

We will write a point $\mathbf{a} = \begin{pmatrix}a_x \\ a_y \\ a_z\end{pmatrix}$ as $\overline{\mathbf{a}} = \begin{pmatrix}a_x \\ a_y \\ a_z \\ 1\end{pmatrix}$. So we will just add a fourth coordinate with a $1$ at the end. Calculating a difference vector between to points $\overline{\mathbf{a}} - \overline{\mathbf{b}}$ then obviously results in a $0$ in the last place.

This makes writing $\operatorname{d}$ and defining the plane pretty easy:

$$\begin{align*} \operatorname{d}(\mathbf{x}) &= \mathbf{n}\cdot\mathbf{x} +d \\ &= \begin{pmatrix}x_x \\ x_y \\ x_z \\ 1\end{pmatrix} \cdot \begin{pmatrix}n_x \\ n_y \\ n_z\\ d \end{pmatrix} \\ &= \overline{\mathbf{x}} \cdot \mathbf{l} \end{align*}$$

So we can just write our plane values in a 4D vector and do a dot product to find $\operatorname{d}$. The definition of the line and calculating the intersection point with the computed $t$ value does not change when we augment our vectors. We just treat our vectors as normal 4D vectors when doing any computations, such as addition or scaling. This will come in very handy later.

Now on to finally compute the clipping (we already have the most important stuff down now).

We define a polygon by its indexed vertices $\mathbf{P}_0, \dots, \mathbf{P}_n$. The edge $i$ is just the line from point $\mathbf{P}_i$ to $\mathbf{P}_{i+1}$, where we wrap around the last index. So the last edge is $\mathbf{P}_n$ to $\mathbf{P}_0$.

The reason we don't just work with triangles here, is that if we use multiple clipping planes, a triangle might become a quadriliteral or just some general polygon after multiple clippings. It will stay a convex shape though!

The Sutherland-Hodgman algorithm works by producing a sequence of vertices that correspond to the clipped polygon. If nothing needs to be clipped, the original polygon is produced, if it is fully clipped, you will get an empty list.

It works the following way:

1. Define an empty output array points

2. Iterate over all polygon edges, starting from the last one (from point $n$ to point $0$)

   1. Compute distances for both edge points

   2. Handle edge (4 cases)

      1. Both points outside? -> Nothing to do

      2. Both points inside? -> Append endpoint of edge to points

      3. Startpoint inside, endpoint outside? -> Append intersection to points

      4. Startpoint outside, endpoint inside? -> Append intersection and then the endpoint to points

3. Return points as result

As we are starting with the last edge, we only ever output the end of a line, as that ensures, that with no clipping, the first vertex of the result will be the first vertex of the input. Notes for the 4 cases:

1. When both points are outside, that edge can be discarded. Cases 3 and 4 handle the case that the polygon goes out of the plane and comes back in at some point
2. When both points are inside, we do not need to clip. As the last edge in the point definitions is implicit, we output only the endpoint
3. The polygon edge will stop at the intersection, so we replace the actual endpoint (it is clipped) with the intersection
4. The edge comes back into the plane. Normally we would just put the endpoint there, but since we were outside before, the startpoint was clipped and couldn't have been added to the list before, so we need to add the intersection before the endpoint

Determining, where a point lies with respect to the clip plane is simply done by checking the signs of the distances, positive for inside, negative for outside. Those distances can then also be used for the intersection calculation.

This took some effort, but we can now implement the steps of the algorithm ourselves! For now, we will use the inbuilt JavaScript rendering to visualize it, but in the next step, we will integrate it into our framework. This will once again be just a bit of bookkeeping, as we need to do something to handle drawing polygons with our current triangle mechanism.

This implementation will take an array of clip planes and just walk through them to clip the results one after another.

Exercise:

• Implement the clipping algorithm described above in clip_polygon in the clip.js file

  • Follow the instruction in the code

Solution: