Many textbooks describe GA as being ‘exception-less’, i.e., exception cases are not always needed. E.g., the meet of two planes is always a line—although the line may be ideal if the planes are parallel. How far does this concept extend? I’m using PGA-3D to solve the problem of bisecting a triangular-mesh volume with a cutting plane. Namely, given an implicit 3D volume defined by a closed, coherently oriented triangular mesh and an arbitrary 3D plane which bisects the mesh volume, find the volume of the mesh ‘above’ and ‘below’ the plane. So my question is this: can this bisection of the triangular mesh be done without resource to conditional, branching code?
My current strategy is as follows. First, I compute the distances of all mesh vertices to the cutting plane. Next I disposition each triangle based on the sign of these vertex distances: positive, negative, and zero. For example, if the signs are all positive, the triangle is ‘above’ the cutting plane. Likewise, if the signs are all negative, the triangle is ‘below’ the cutting plane. If the signs are mixed or zero, then the triangle is either bisected by the plane, touches the plane, or lies entirely in the plane. If the triangle is bisected, it is fragmented along the meet of the plane and the triangle. This process allows me to update the mesh and split it into exactly two sets, those ‘above’ and those ‘below’ the plane. Due to the multiple topologies that need to be handled (namely, 3^3 =27 combinations, as each vertex distance can be negative, zero, or positive), I see no way to avoid coding multiple conditional branches—the disposition of each topological condition needs to be uniquely handled. Or does it?
@DrQ There are limits to the seamlessness of GA expressions. It can’t handle all the possible ways a plane can cut a triangle.
May I suggest couple of improvements to the proposed algorithm? First, what about considering just two cases: on or above the plane (d >= 0) and below the plane (d<0). This may result in degenerate triangles but the area formulas will I think work correctly. Also, if the plane is not completely on one side, then you know that exactly one vertex will be above/on and 2 vertices below, or vice-versa. The cutting line will cut the triangle into a smaller triangle (that includes this singleton vertex) and a quadrilateral (including the two vertices on the other side of the plane). By examining the array of distances you can figure out a set of indices for performing this cut.
As I’m in the process of learning to write JavaScript I took this as an opportunity to write a short ganja demo (in 2D) that illustrates the main ideas with the example of 1 triangle and 1 cutting line. Perhaps it can be of use to you.
Suggestions always welcome. Looking through the code you highlight the kind of things I was looking for. E.g., being able to exploit degenerate cases (i.e., quadrilaterals that degenerate to triangles when the cutting plane cuts through a vertex, triangles that degenerate to zero area, etc.). I’ve also thought of first computing the meet points of all mesh edges with the cutting plane, then sorting out the results afterwards. Or maybe working in geometric duality (e.g., a point is the intersection of all lines through it). By the way, I’ve implemented both 2D and 3D versions of PGA in Matlab. I’m working on single-instruction, multiple-data hardware, so over-computing to exploit parallelism, even if I throw away some of those results later, is okay. Thanks for the response.
The implementation of PGA in MatLab sounds interesting! Do you anticipate making this implementation available to others? That would be a great resource for the MatLab community, I think.
I’m using Matlab to: understand PGA math, understand the issues in implementing PGA code, and to solve a sufficiently hard problem. Hence, my code is being written in the most readable form, not for efficiency in speed or memory usage. The motivation for my question is because matlab becomes very efficient if the underlying code is ‘vectorized’, meaning one avoids for-loops and relies on the fact that Matlab inherently understands arrays and matrix operations. This in turn implies that exception-less algorithms should be very fast.
I am mostly interest in Euclidean 3D-space problems, so my work will be expanded to include: Vector-GA, Projective-GA, and Conformal-GA models. I started in the middle with PGA. I wish to understand the trade-offs between computational & memory costs verse the expressive power of a particular model. Solving the same problem in all three models will help in that understanding.
Normally I try make my efforts available, but this go around is not very bullet-proof. The visualization part is very problematic, but gets the job done. Most Matlab GA/Clifford packages are very broad, meaning they handle arbitrary signatures, and basis sizes. So, if the Matlab community is interested in this rather narrow slice, they should find it useful.
@tionkje thanks for the examples. I’m note very proficient with this language, but I think I can decipher enough to understand what is going on. I saw the same inconsistency problems in my code which account for the majority of the 27 topological cases. I’m trying to arrange my code so that these problems are down in the computational noise floor or go ‘ideal’. Again, thanks.
If the task is just to find the volume above and below, I think there is no need for an explicit triangulization, or for calculating ‘caps’. I haven’t tested it yet but I suspect the following algorithm should do the trick:
Start with two running sums, lets call them A_\text{tot} and B_\text{tot} for Above and Below respectively. Call M the motor that takes one of the basis planes P_\text{basis} to the cutting plane P_\text{cut}. (M = \sqrt{P_\text{basis}P_\text{cut}}). Now for each triangle of the mesh:
transform its vertices v_i using w_i = \tilde M v_i M.
calculate the triangle plane in this frame as P_\text{triangle} = w_1 \vee w_2 \vee w_3
check for each v_i if it is above or below the basis plane P_\text{basis}. (a single coordinate check)
if all three w_i are on or above P_\text{basis} add P_\text{triangle} to A_\text{tot}
if all three w_i are on or below P_\text{basis} add P_\text{triangle} to B_\text{tot}
else find the single vertex, lets say w_1 that is on the other side as the other vertices. w_2,w_3.
create two new vertices u_2 = (w_1 \vee w_2) \wedge P_\text{basis} and u_3 = (w_1 \vee w_3) \wedge P_\text{basis}
find the plane of this small triangle P_\text{small} = w_1 \vee \hat u_2 \vee \hat u_3.
if w_1 was above the basis plane, update the sums as A_\text{tot} = A_\text{tot} + P_\text{small} and B_\text{tot} = B_\text{tot} + P_\text{triangle} - P_\text{small}
if w_1 was below the basis plane, update the sums the other way around.
When done, the volume above the cutting plane will be \frac 1 6 \lVert A_\text{tot} \rVert_\infty and the volume below will be \frac 1 6 \lVert B_\text{tot} \rVert_\infty.
Transforming all the vertices with the motor, then checking if a single coordinate of each is above or below the basis plane seems more expensive than simply computing the oriented distance s_i of each v_i to the cutting plane P_\text{cut} directly as s_i = v_i \vee P_\text{cut}. (Unless, I suppose, if the underlying code can exploit sparse multi-vectors.) Otherwise your algorithm is the same as mine, except I work with \{v_i, P_\text{cut}\} not \{w_i, P_\text{basis}\}—with one exception: the ‘caps’.
I’m confused why you do not need to compute the ‘caps’? Typically volume formulae for a mesh requires the mesh be closed, i.e., the mesh divides inside and outside. Don’t you need to compute properly weighted ‘cap’ planes by a scaling a normalized cutting plane P_\text{cut} (or basis plane P_\text{basis} in your case) by the area of the plane within the mesh in order to ‘close’ the mesh of either above or below volumes? Are you saying the volume formula works if the mesh is not closed, provided the edges of a mesh’s hole are co-planar?
Glad you ask - I should’ve been more clear on the tricks in that algorithm. The volume formula (\frac 1 6 \lVert \sum P_i \rVert_\infty) works by adding signed volumes of the tetrahedrons (or more accurately, parallelepipeds, hence the \frac 1 6) from the three vertices of your triangle to the origin. This is why you can avoid the explicit construction of caps (as any triangle in a plane that contains the origin will create a zero volume tetrahedron). It is also why I’m moving the vertices so that the cutting plane becomes a basis plane (you could make a similar version that only translates the points so that the origin is in the cutting plane). This is a simple, non-conditional, linear operation which should easily parallelize. Now the only conditions needed are to check if vertices are above or below that basis plane so that you can update the correct running sums with the volume of the small tetrahedron and the difference of the volumes of the large and small one respectively (this difference in volume trick also only works when the cutting line and the origin are both in the cutting plane).
Combined it removes the complicated step of constructing cap geometry and guaranteeing mesh closedness completely and additionally simplifies the remaining processing.
Below is an illustration of the situation in 2D - for areas. Blue triangles have a positive area, red ones a negative area. The left geometry is capped at the bottom, so it generates one big ‘negative area’ triangle that takes out the excess area from the blueish triangles. On the top right you can see that leaving the ‘cap’ out is not an option, but for the one below at the cutting line, the cap is not needed. (as it effectively would create a zero-area triangle anyway).
Let us know if you find a better method
Cheers!
(disclaimer this is a method I describe in an upcomming paper called ‘clean up your mesh!’)
@enki, Now I see (pun intended). That formula now makes perfect sense. Brilliant!
Further questions:
Should 3rd bullet say \quad “check \color{red}w_i if it is above or below the basis plane P_\text{basis}…”
or should 4th and 5th bullet items say \quad “if all three \color{red}v_i are on or…”? One of these appears wrong.
I think the motor is backwards, either M = \sqrt{P_\text{cut} P_\text{basis}} or w_i = M v_i \tilde{M}. The later appears to place the volume mesh correctly.
Correct on both remarks. The third bullet should be w_i and the conjugation should be w_i = M v_i \tilde M.
(you could alternatively transform just the plane equation using the motor instead of transforming the three vertices, and to optimize even further only the \mathbf e_0 coefficient of the transformed plane needs to be calculated and summed - in this case you’ll need v_i \vee P_\text{cut} to determine the side the v_i are on.)
@enki, Actually you need only transform the final sums A_\text{tot} and B_\text{tot} to determine the volumes as \frac 1 6 \| M A_\text{tot} \tilde M \|_\infty and \frac 1 6 \| M B_\text{tot} \tilde M \|_\infty. The oriented distance equation d_i = \hat v_i \vee \hat P_\text{cut} is considerably less expense to compute than a full sandwich product (M v_i \tilde M) across all the vertices. (Also, this leaves the v_i untransformed with the added benefit of having to compute P_\text{triangle} = v_1 \vee v_2 \vee v_3 only once for use by multiple cutting planes P_{\text{cut}, k}.)
I look forward to your paper ‘Clean up your mesh!’
