The solution: shear matrices
A shear matrix is better explained visually.
In this case we applied a vertical shear (parallel to the y-axis). Our
x-values are the same, but our y-values have been transformed in a way
that they match the line \(y'=mx+y\) (see the matrix representation).
It has almost all of the properties required for our 2.5d rotation -
we're only distorting the matrix on one axis and preserving the other.
The only issue is that our transformed object will appear longer.
In my case, I do not need to know the angle - I can use the gradient
between the silo \(\mathbf{S}\) and target \(\mathbf{T}\).
\[m=\frac{\mathbf{S}_y-\mathbf{T}_y}{\mathbf{S}_x-\mathbf{T}_x}\] \(m\)
can also be obtained with a trig ratio, yielding the following shear
transformation that involve angles rather than gradients. \begin{align}
\text{Parallel to the y-axis, Vertical shear}\quad \mathbf{X'} &=
\begin{bmatrix} 1 & 0\\ \sin(\theta) & 1 \end{bmatrix}\mathbf{X}\\\\
\text{Parallel to the x-axis, Horizontal shear}\quad \mathbf{X'} &=
\begin{bmatrix} 1 & \sin(\theta)\\ 0 & 1 \end{bmatrix}\mathbf{X}\\\\
\end{align} To resolve our issue with the transformed object being
'lengthened', we need to scale it back.
Let's choose two points, \(\mathbf{P_1}\) and \(\mathbf{P_2}\).
We obtain the original distance, \(d\) and the transformed distance,
\(d'\) \begin{align} d &= \|\mathbf{P_1}-\mathbf{P_2}\| \\ &=
\sqrt{(\mathbf{P_1}_x-\mathbf{P_2}_x)^2+(\mathbf{P_1}_y-\mathbf{P_2}_y)^2}
\\ d' &= \|\mathbf{P'_1}-\mathbf{P'_2}\| \\ &=
\sqrt{(\mathbf{P'_1}_x-\mathbf{P'_2}_x)^2+(\mathbf{P'_1}_y-\mathbf{P'_2}_y)^2}\\\\
\text{ where }\quad\mathbf{P'_1},\mathbf{P'_2} &= \begin{bmatrix} 1 &
0\\ \frac{\mathbf{S}_y-\mathbf{T}_y}{\mathbf{S}_x-\mathbf{T}_x} & 1
\end{bmatrix}\mathbf{P_1},\mathbf{P_2}\\ &= \begin{bmatrix} 1 & 0\\
\sin(\theta) & 1 \end{bmatrix}\mathbf{P_1},\mathbf{P_2}\\ \end{align}
Note that \(d\) is also the distance from the silo \(\mathbf{S}\) and
target \(\mathbf{T}\)
\(\mathbf{P_1}\) and \(\mathbf{P_2}\) are not the same as \(\mathbf{S}\)
and \(\mathbf{T}\), because we've created these points with
\(\begin{bmatrix}0\\0\end{bmatrix}\) as the origin,
not \(\mathbf{S}\). We will translate all of the points later so
that our starting point is \(\mathbf{S}\) after rotating.
However, the distances are the same, so we might as well use
\(\mathbf{S}\) and \(\mathbf{T}\) to calculate \(d\). \[d =
\|\mathbf{P_1}-\mathbf{P_2}\| = \|\mathbf{T}-\mathbf{S}\|\] Our scale
factor is simply the ratio between the two. \begin{align} s &=
\frac{d}{d'} \\ &=
\frac{\|\mathbf{T}-\mathbf{S}\|}{\|\mathbf{P'_1}-\mathbf{P'_2}\|}
\end{align} Putting it all together, our 2.5d rotation is now:
\begin{align} \mathbf{X'} &= s\begin{bmatrix} 1 & 0\\ m & 0
\end{bmatrix}\mathbf{X}\\ &=
\frac{\|\mathbf{T}-\mathbf{S}\|}{\|\mathbf{P'_1}-\mathbf{P'_2}\|}\cdot\begin{bmatrix}
1 & 0\\ \frac{\mathbf{S}_y-\mathbf{T}_y}{\mathbf{S}_x-\mathbf{T}_x} & 1
\end{bmatrix}\mathbf{X}\\ \text{ or }\quad &=
\frac{\|\mathbf{T}-\mathbf{S}\|}{\|\mathbf{P'_1}-\mathbf{P'_2}\|}\cdot\begin{bmatrix}
1 & 0\\ \sin(\theta) & 1 \end{bmatrix}\mathbf{X}\\\\ \end{align} So we
need a 'pre transformation' on two points, \(\mathbf{P_1}\) and
\(\mathbf{P_2}\) in order to obtain our scale factor
The limitation of this technique is that we need to perform the rotation
about \(\begin{bmatrix}0\\0\end{bmatrix}\), the origin. After this we
can translate the rotated matrix to match the silo and target
location.
You can use homogeneous coordinates to achieve this (see the yellow box
earlier), or programatically do it.
We also have to solve one last issue with this method: it only works
within the domain \(\left(\frac{\pi}{2}, -\frac{\pi}{2}\right)\).
This is because the gradient can only describe so much information.
For instance, the line \(f(x)=1\cdot x\).
We cannot tell if the line is going from the top-right to bottom-left,
or bottom-left to top-right because it describes both situations.
We need to adjust the signs on the shear matrix according to the
quadrant our target is in.
import math
import numpy as np
def quad(A,B): #A is our origin/fix point, B is our other point.
if B[0] > A[0] and B[1] > A[1]:
return 1 #First quadrant
elif B[0] < A[0] and B[1] > A[1]:
return 2 #Second quadrant
elif B[0] < A[0] and B[1] < A[1]:
return 3 #Third quadrant
else:
return 4 #Fourth quadrant
#By elimination. We could do another if but it's unnecessary
# elif B[0] < A[0] and B[1] > A[1]:
def trajectory(distance, n):
# Here we create a list of points. n is the amount we create
# (higher n results in greater precision)
# Not going to include it here - uses bezier interpolation and stuff to create
# the illusion of a ballistic trajectory with gravity (not a simple quadratic function)
#
return points
S = np.array([10,4]) #Silo position vector
T = np.array([5204,954]) #Target position vector
d = np.linalg.norm(T - S) #Get the distance: d = ((Tx-Sx)^2+(Ty-Sy)^2)^(0.5)
m = (T[1]-S[1])/(T[0]-S[0]) #Gradient of line between target and silo.
X = trajectory(d, 50) #50 is a constant for n. only affects trajectory precision
#X now contains a 'flat trajectory', a list of points which have the correct
#distance but the wrong angle and offset.
#It does not intersect with the target T.
def shear(m, X):
q = quad(S,T) #Silo and target is passed to the quadrant finding function.
#q is the quadrant
if q == 1 or q == 4:
return np.array([[1,0],
[m,1]]) #shear matrix.
else:
return np.array([[-1,0],
[m,-1]]) #shear matrix, but for quadrants 2 and 3 (negative x).
# AN easier way to do this would be to just compare S[0] > T[0]. Forget about the quadrant stuff.
# because we only need to know if it's in the negative or positive x.
SHEAR = shear(m, X) #we get the shear matrix, which the function returns
P_ = np.matmul(SHEAR, X[(0,-1),]) #first, take a slice of the array X to obtain the first (P1)
# and last (P2) points. Then multiply this new 2x2 matrix
# containing P1 and P2 by the shear matrix
# to get P_, which contains P'1 and P'2
d_ = np.linalg.norm(P_[1] - P_[0]) #adjusted distance
s = d/d_ #scale factor
X_ = s * np.matmul(SHEAR, X) + S # shear, then scale every point by s.
# Then add S, the silo position, to each position to yield the
# correct positions
# X_ (X') is now a numpy array containing the shifted position vectors
# that have been transformed through the 2.5d rotation.
# The trajectory will still look 'realistic', not curving to the side but still
# 'visually correct' and will be functionally correct in that it
# starts at silo S and ends at target T.