# Determine rotation axis (matrix) based on two positions and an offset vector

I am working on a problem where some algorithm is creating a triangle mesh, and each triangle is neighbour to another and they always share an edge. Once I determined which edge is the shared edge with the next triangle, I have to build a local coordinate system, and then use input coordinates on the local origin (which is the mid-point of the shared edge) to get the third point of the new triangle.

This is generally working fine for me, however I am running into a problem where I get a "bad triangle" (where the third point is along the shared edge, making it a line rather than a triangle). The subsequent triangle of that "bad triangle" then cannot build the local coordinate system because the last triangle as a line and as such doesn't have a "proper" normal.

However, I have the source data, and as such I know what the next point should be. And I assume based on that data, I can determine what the local coordinate system needs to be in order to correctly arrive at that target point.

I assume that generally there could be an infinite amount of local coordinate systems by which I could arrive at the correct target location, however I have specific coordinates that have to be used in that local coordinate system in order to arrive at the destination - therefore I assume (hope!) that it should be possible to calculate a single result based on the input data that I have.

(Note: I do have a fallback algorithm which is supposed to be used in order to arrive at the correct location, if the previous triangle was invalid, however that does not work at all - therefore I want to see if I can find the correct algorithm by finding the local coordinate system for this case)

Background info:

This is how to normally (ie when you have a valid previous triangle) determine the local coordinate system used to get to the third point of the new triangle (right hand coordinate system):

• set as origin the mid-point of the shared edge (ie mit between v0 and v1)
• build the x-axis from the shared edge (ie v1-v0)
• build the z-axis from the normal of that previous triangle (use direction vector of third point of previous triangle and above calculated origin, and
cross it with the above calculated x-axis)
• build the y-axis as cross product of z- and x-axis

My concrete example data:

I am providing the concrete data as a reference to check the resulting calculations that somebody can hopefully help me to find.

Coordinates of the previous triangle:

v0 = { 1214.7296 -2190.7231 -16437.4473 }
v1 = { 1189.5955 -2157.4277 -16397.6074 }
v2 = { 1199.9772 -2171.1804 -16414.0625 }


Shared edge:

v0-v2


Origin (midpoint of shared edge)

{ 1207.3534, -2180.95175, -16425.7549 }


Target-point (ie third point of the new triangle)

{ 1216.7325 -2177.5632 -16393.6191 }


Coordinates applied to local coord system in order to arrive at this target point

{ -20.063538445448748486743693311438, -10.031769222724374743384233267162, 25.079423056810936358448196556462 }


(Note: I used the highest precision I got from the calculator for these local coords)

I would highly appreciate any input to help me find a formula/algorithm that allows me to find a "local coordinate system" based on an origin point, a "target point" and the x/y/z coord which needs to be applied to that coord-system in order to arrive at the target point. So I do know how to normally traverse the list of triangles, however I am running into the problem where the input coordinates for the local coordinate system are (x,0,0) thereby creating a line rather than a triangle. Then the following triangle will run into a problem as it cannot correctly build the y- and z-axis. By finding a formula/algorithm to "move backwards" (from the (in this case) already known target point and the transformation coordinates) I hope to be able to later deduce a general formula to apply for such (x,0,0) cases.

• So, you have a list of triangles but some of them might be degenerate, and you run into trouble when you try to build a local coordinate system on one of these bad triangles. Is that a fair assessment of the situation? If so, how about keeping the last “good” set of axes (basis vectors) and when you detect one of these bad cases, use the saved coordinate system, adjusting the origin as necessary to the midpoint of the shared edge common to the bad triangle and the next one? – amd Aug 9 '16 at 21:35
• Yes, that assessment is correct.The official procedure for cases when the y- or z-axis are 0 (or less than minimum floating point precision) is to build an ortho representation using the x-axis and y(0/1/0) as input to build the z-axis and then the final y-axis as cross product of x and z. However in the case of this data it doesn't bring me to the "target point" (I know how the model should look like, so I know where it went wrong). That's why I tried to find the coord system that was used to arrive there, and hope to deduce from that what the real procedure for degenerate data is – Marko T. Aug 10 '16 at 2:29

## 2 Answers

It’s possible to recover the local frame in the non-degenerate case. The key pieces of information that enable this are that the local basis is right-handed orthonormal and you have one of the axes.

Let $\mathbf m$ be the midpoint of the shared edge, with associated unit vector $\mathbf u_1$ (the local $x$-axis). We only need to find one other unit basis vector, say $\mathbf u_2$ (local $y$-axis), because $\mathbf u_3=\mathbf u_1\times\mathbf u_2$. For the target vertex $\mathbf v$, we have its coordinates $[\mathbf v]_G$ in the global frame and $[\mathbf v]_L$ in the local frame. Let $U$ be the matrix with these local basis vectors as its columns. This matrix is orthogonal (in fact, it’s a rotation), so in particular $U^{-1}=U^T$. We then have $$[\mathbf v]_L=U^T[\mathbf v-\mathbf m]_G,\tag{1}$$ which yields two linear equations for the components of the unknown $\mathbf u_2$. The orthogonality condition $\mathbf u_1\cdot\mathbf u_2=0$ gives us a third. Solving this system is a straightforward computation. You can, of course, instead take $\mathbf u_3$ as the unknown and set $\mathbf u_2=\mathbf u_3\times\mathbf u_1$.

Note that this procedure works whether or not the “previous” triangle was degenerate. It only uses the shared edge, not the entire triangle. On the other hand, if $\mathbf v$ is colinear with the shared edge, the procedure fails and the linear system only has the trivial solution. That’s OK, though, since we can choose any other basis vectors for $\mathbf u_2$ and $\mathbf u_3$ without changing $[\mathbf v]_L$.

The reason you only get two equations from (1) above is that the top row of the matrix is the known vector $\mathbf u_1$, so you end up with $\text{const}=\text{const}$ for the $x$-coordinates. Comparing these values is a good sanity check. In practice, they’re not going to be exactly equal because of the loss of precision that comes with all of the multiplication and root extraction, but if they don’t match, then you’ve got the wrong local $x$-axis. In that case, you’ll need to dig deeper.

Here’s a sketch of a procedure to find the local $x$-axis. First, translate so that the origin is at $\mathbf m$ and let $\mathbf v'=\mathbf v-\mathbf m$, that is, the coordinates of $\mathbf v$ in this translated frame. Since the components of a vector are its projections onto the basis vectors, each of these components lies on circle centered on the line defined by $\mathbf v'$ and in a plane orthogonal to it. Wlog assume that $\mathbf v'$ is in the direction of the positive $z$-axis. Then each circle can be parametrized as $$\mathbf r_k(t_k)=w_k(\sin{\alpha_k}\cos{t_k},\sin{\alpha_k}\sin{t_k},\cos{\alpha_k}),$$ where $\cos{\alpha_k}=w_k/\|\mathbf v'\|$, which then gives $$\mathbf u_k(t_k)=(\sin{\alpha_k}\cos{t_k},\sin{\alpha_k}\sin{t_k},\cos{\alpha_k}).\tag{2}$$ That this looks like spherical coordinates is no coincidence: the $\mathbf u_k$ all lie on a unit sphere. The mutual orthogonality requirement $\mathbf u_j\cdot\mathbf u_k=0\,(j\ne k)$ works out to be $$\tan{\alpha_j}\tan{\alpha_k}\cos{(t_j-t_k)}=-1$$ so once you fix a value for one of these parameters, the other two are determined up to sign. You will need to use other information to find the appropriate value of $t_1$ and resolve the remaining sign ambiguities. Another way to find the other possible basis vectors once you’ve chosen $\mathbf u_1$ is to find the intersection of its perpendicular plane with the corresponding circles.

All of the above depends on your being able to determine the local frame’s origin $\mathbf m$. If you also have to reverse-engineer that, another two degrees of freedom are introduced since $\mathbf m$ must lie on a sphere of radius $\|\mathbf w\|$ centered at $\mathbf v$.

Having said all of that, I’d recommend an interactive experimental approach. Set one or two of the bad cases up in something like GeoGebra. Build the circles that I described above and the two possible basis sets that they define, then move them around to see if they line up with anything plausible.

• Thanks for your reply! Gonna have to think through that (am a bit rusty lol) to make sure I understand it well (and before I could mark it as answer), and will try it with known points first. Just a quick question: Is it possible to determine all 3 axis as well, in case it turns out the local x-axis is not along the shared edge? (I strongly hope it is still along the shared edge but just wonder if there is still just one solution if it is not) – Marko T. Aug 10 '16 at 13:45
• @MarkoT. I tried it out with the data you provided and did crank out a solution. It was a good exercise, because when I followed the procedure for building the local frame in your question, I came up with rather different local coordinates for the target, so I was able to find the frame you used. – amd Aug 10 '16 at 18:21
• @MarkoT. If you can’t nail down one of the axes in some way, then the problem is likely underconstrained. There are effectively 5 degrees of freedom instead of 3. You might be able to make use of the requirement that the basis consist of unit vectors, but I’m not sure that’ll be enough. – amd Aug 10 '16 at 18:24
• The whole thing is built upon a spec (which sadly has many holes, though, most of which I have already reverse engineered by now). Contrary to the spec it actually turned out that I only get the correct result (for target point) when I negate the x- and y-axis. As I mentioned, the original fallback for degenerate data was to use a function that made an ortho representation from an input x-axis, and it seemed like they take this v1-v0 x-axis, but it didn't work (they wrote "unit x-axis" actually, but didn't work either). So I am trying what it really is... – Marko T. Aug 11 '16 at 2:27
• @MarkoT. Unfortunately, there’s always the spec and what actually gets implemented. Back in my Xerox days, there was one guy who kept implementing to the spec and would end up with incorrect behavior as a result. I’ll take another look at this tomorrow and see what can be done if the assumed $x$-axis doesn’t match the data. BTW, the difference that I saw when I tried the forward direction wasn’t a simple sign difference, which can be fixed by flipping axes. The coordinate values were very different from yours to boot. – amd Aug 11 '16 at 5:28

I didn't quite understand what is given and what you are trying to find. Is your input data a triangle with vertices $v_0, v_1, v_2$ and three coordinates $(x, y, z)$ if the new (target) point in local coordinates? Is it guaranteed that the new triangle should always be a true triangle? That is, is it true that your local coordinates never have the form $(x,0,0)$?

Anyway, if I understand correctly, then the local coordinate system makes sense. Basically you take the midpoint $m=(v_1+v_0)/2$ as as the origin, and the three orthonormal vectors $$e_1 =\frac{v_1-v_0}{\|v_1-v_0\|}$$ $$e_3 =-\frac{(v_0-v_2)\times (v_1-v_2)}{\|v_0-v_2\| \|v_1-v_2\|}$$ $$e_2=e_3\times e_1$$ Then the target point is calculated as $$v_{target}=m + x e_1 + y e_2 + z e_3$$ To have an actual triangle you need at least one of the coordinates $y$ or $z$ to be nonzero.

Is this what you have been trying to achieve?

Edit

OK let's see if this is what you are looking for.

Assume you have a preceding triangle $u v_1 v_2$ which is a true triangle, non degenerate, and we focus on the current triangle $v_0 v_1 v_2$, which shares with $u v_1 v_2$ the common edge $v_1 v_2$ and is allowed to be degenerate ( that is point $v_0$ lies on the line determined by $v_1 v_2$ ). You know the next point $w$ which forms the subsequent triangle $w v_0 v_1$ and it's coordinates $(x,y,z)$ in an unknown orthonormal coordinate system centered at the midpoint $$m =\frac{v_0+v_1}{2}$$ of the segment $v_0v_1$ and aligned with the edge $v_0v_1$.

Step 1 Form the unit vector $$e_1 =\frac{v_1-v_0}{\|v_1-v_0\|}$$

Step 2 Check that the dot product $((w-m)\cdot e_1)=x.$

If not, something is not right with your data.

Step 3 Check if $\| w-m-xe_1\|^2=y^2+z^2$.

If not, something is not right with your data.

Step 4 Form the unit vector $$w_3=\frac{(v_1-u)\times (v_0-u)}{\|(v_1-u)\times (v_0-u)\|}.$$ It is orthogonal to $e_1$ by construction, where in the case of $v_0v_1v_2$ being non degenerate, we can replace $u$ by $v_2$ everywhere in the steps to follow.

Step 5 Form the unit vector $w_2=w_3\times e_1$.

Remark the two remaining unit vectors of the coordinate system we are looking for can be written down as $$e_2=\lambda w_2 + \mu w_3$$ $$e_3=-\mu w_2 + \lambda w_3$$ $$\lambda^2+\mu^2=1$$ Therefore they should satisfy the conditions $$((w-m)\cdot e_2)=y$$ $$((w-m)\cdot e_3)=z$$

Step 6 Solve the linear system $$((w-m)\cdot w_2)\lambda+((w-m)\cdot w_3)\mu=y$$ $$-((w-m)\cdot w_2)\mu+((w-m)\cdot w_3)\lambda = z$$ in the non-degenerate case, when $w$ does not lie on $v_0v_1$, and check whether $$\lambda^2+\mu^2=1$$ The system is easy to solve as the inverse matrix can be very easily computed (something like minus the transpose times it's determinant).

or alternatively, in the degenerate case, simply take $e_2=w_2$ and $e_3=w_3$.

Step 7 Form the orthonormal vectors (in the non degenerate case) $$e_2=\lambda w_2 + \mu w_3$$ $$e_3=-\mu w_2 + \lambda w_3$$

Remark The point $u$ does not have to be a vertex of the triangulation, it can be any point available that does not lie on the line $v_0v_1$.

• What is given: Based on the "previous triangle" I have the "origin" for the coord system. I know where I should "arrive" (target point). And I know the x/y/z coords that are to be applied to that coord system. What I need: The actual (concrete) coordinate system used to go from "origin" to "target" using those x/y/z coordinates (or a formula to find it) - in the hope of using that concrete coord system to deduce a general formula to use if (x/0/0) is the case. Sorry if my question was not clear enough (am not a native speaker, and lack some of the mathematical knowledge required too). Thanks! – Marko T. Aug 9 '16 at 13:18
• Is it always the case that the x axis of the local coordinate is aligned with the edge $v_0v_1$? In other words, is it always true for your data that $(v_{target}-m). e_1=x$ (this is dot product)? – Futurologist Aug 9 '16 at 15:31
• In my general cases it is always true. However, I don't know what's happening in the case where one triangle gets (x,0,0) as input and the subsequent triangle has to use an "alternative local coord system". That's why I was hoping to find (based on the fact that I happen to know what the next point (target point) should be) a way to determine what coordinate system was used to arrive at that point (based on x,y,z being {-20.0635, -10.0317, 25.0794}). I was hoping there was a way to find one single specific local coord system, that using those input coords would get from origin to target – Marko T. Aug 9 '16 at 16:13
• So that "background info" part I added to my question was really just that: background info on how things usually work. The key problem is figuring out what to do if they don't work. And by having data of the complete set of triangles available, I was hoping there was a way to deduce the local coord system based, and once I have that I was hoping to deduce a general algorithm for those cases (ie (x,0,0) and how to calculate subsequent coord systems). Sorry, I hope I could explain myself well enough. Thanks for taking the time to help me analyze. – Marko T. Aug 9 '16 at 16:16
• Hi again, thanks for the update on your answer, have to look at it a bit tonight (and probably better tmw morning when I am fresh), and will first go through it with known axis to confirm I understand it well. – Marko T. Aug 10 '16 at 13:55