In principle you need two rotations and a matrix-multiplication. Call the recentered matrices Q ad S and P as T.
Then find the rotation s, which rotates S to "triangular" form of its coordinates, so that the coordinates of the first point becomes $[x_{s,1},0,0]$, of the second becomes $[x_{s,2},y_{s,2},0]$. This can be done with your ansatz of using unknowns, if you assume three rotations (for each pair of planes 1 rotation, and gives a 3x3-matrix whose 3'rd column is zero and describes thus a triangle in a plane)
Then do the same with another rotation t which rotes T in the same manner.
Then $T = S \cdot s \cdot t^{-1} = S \cdot r$ where r is the matrix for the complete rotation.
A bit of pseudocode (code in my proprietary MatMate-tool):
S = Q - Meanzl(Q) // do the translation to the origin
T = P - Meanzl(P)
// get the rotation-matrices which rotate some matrix to triangular shape
ts = gettrans(S,"tri")
tt = gettrans(T,"tri")
// make one rotation-metrix. the quote-symbol means transposition
tr = ts * tt'
// difference should be zero
CHK = T - S*tr
Example:
We assume
S and
T being centered. Let
$ \small \text{ S =} \begin{bmatrix}
14.469944&22.964690&-7.581631\\
-15.275348&5.923432&23.720255\\
0.805404&-28.888122&-16.138624
\end{bmatrix}
$
and
$ \small \text{ T =} \begin{bmatrix}
22.808501&2.515200&16.361035\\
8.393637&-5.071089&-27.109127\\
-31.202138&2.555889&10.748092
\end{bmatrix}
$
Now you can solve for a rotation in
y/z-plane, which makes the entry in $S_{1,3}=0$. The rotation-parameters are some cos/sin-values. Apply this and you get
$\small \text{ S}^{(1)} = \begin{bmatrix}
14.469944&24.183840&0.000000\\
-15.275348&-1.811476&24.381471\\
0.805404&-22.372364&-24.381471
\end{bmatrix}$
Now you can solve for a rotation in
x/y-plane, which makes the entry in $S_{1,2}=0$. The rotation-parameters are some other cos/sin-values. Apply this and you get
$\small \text{S}^{(2)} = \begin{bmatrix}
28.182218&0.000000&0.000000\\
-9.397482&12.178056&24.381471\\
-18.784736&-12.178056&-24.381471
\end{bmatrix}$
After that a third set of rotation-parameters cos/sin-values make
S triangular. It looks then like this
$\small \text{ S}^{(3)} = \begin{bmatrix}
28.182218&0.000000&0.000000\\
-9.397482&27.253645&0.000000\\
-18.784736&-27.253645&0.000000
\end{bmatrix}$
Because the center of your original triangle was moved to the origin, the last column (the z-coordinates) are zero/not needed, since 3 points can always be placed in a plane.
From the three rotation with their cos/sin-values you can construct a 3x3-rotation-matrix, say s.
The same can be done using the matrix T leading to a rotation-matrix t. If S and T describe in fact the same triangle, only rotated, the results are equal: $T^{(3)}=S^{(3)}$.
Then you can use the nice fact, that the inverse of a rotation-matrix is just its transpose, such that with
$\small \text{ tr =} \begin{bmatrix}
0.205215&0.860645&0.466022\\
0.946329&-0.295966&0.129867\\
0.249696&0.414359&-0.875190
\end{bmatrix}$
we get $ T = S \cdot tr $
In principle this can also be solved using the concept of
pseudoinverses: we demand
$ tr = S^{-1} \cdot T $ .
But because
S has reduced rank the inverse means to divide by zero. Using SVD-decomposition (see wikipedia) the pseudoinverse can be computed if the inverse of the diagonal matrix of the SVD-factors is used (where zeros are simply left zero). This should lead to the same solution.
Your code is almost completely correct (and avoids the use of atan2, which is nice). The one problem you have is that when the vector from the disk-center to your point happens to point perpendicular to the disk, i.e., when the test point happens to be exactly above the disk's center, it falls apart: projecting that vector to the $xy$-plane gives the zero vector, whose tip does not lie on the disk's boundary circle.
Here's replacement code (mostly yours):
vec1 = vector between disk center and particle
vec1 = [(5 - 0), (5 - 0), (5 - 0)]
vec1 = [5,5,5]
pvec1 = vec1 with its z-coordinate set to 0.
if (norm (pvec1) = 0 ):
answer = sqrt( norm(vec1) * norm(vec1) + radius * radius )
return
unitvec1 = unit vector in direction of pvec1
unitvec1 = pvec1/norm(pvec1)
unitvec1 = [0.7071, 0.7071, 0]
vec2 = vector between disk center and point on the perimeter closest to the particle
vec2 = disk radius * unitvec1
vec2 = 3 * [0.7071, 0.7071, 0]
vec2 = [2.1..., 2.1..., 0]
vec3 = vector between particle and point on the perimeter closest to the particle
vec3 = vec1 - vec2
So the min distance is
norm(vec3) = 6.8087
That'll handle the special case on which your current code is failing. Note that you might want to replace the test "norm (pvec1) = 0" with "norm (pvec1) < .001" to make the code a bit more numerically stable, and introduce almost no error.
Best Answer
I don't know whether this is a correct solution. But I will make a try.
Let $P_1$, $P_2$, $P_3$ be three closest points to $X$. Consider plane $\pi$ containing $P_1$, $P_2$, $P_3$.
Let $P$ be the orthogonal projection of $X$ on $\pi$. If $P$ is inside triangle $P_1P_2P_3$ then the distance is $PX$, otherwise $\min(\mathrm{dist}(P_1,X),\mathrm{dist}(P_2,X),\mathrm{dist}(P_3,X))$. Here is a picture, for clarification
Well, how to determine that $P$ is inside $P_1P_2P_3$? Just check the equality $$ \mathrm{area}(P_1P_2P_3)= \mathrm{area}(PP_2P_3)+\mathrm{area}(P_1PP_3)+\mathrm{area}(P_1P_2P) $$
In order to determine projectoin $P$ you can use the following formula $$ P=X+tN $$ where $$ N=[P_3P_1,P_2P_1],\qquad t=\frac{\langle P_1-X,N\rangle}{\langle N,N\rangle} $$