I have a bunch of points in 3D all on the same plane, having the general equation $$ \hat{\mathbf{n}}\cdot\left(\mathbf{r}-\mathbf{r}_0\right)=0 $$ The unit normal is computed by using 3 points on the surface.
The best choice of system to investigate these points is the tangential frame of the plane. Therefore, I use the unit normal $\hat{\mathbf{n}}$ for the first axes and compute the other two base vectors the following way:
For starters, an arbitrary vector not parallel to $\hat{\mathbf{n}}$ is needed. For the purpose I choose one of the usual 3D unit vectors $\hat{\mathbf{i}}=(1,0,0)$, $\hat{\mathbf{j}}=(0,1,0)$ or $\hat{\mathbf{k}}=(0,0,1)$. The choice is done based on the size of the angle $\hat{\mathbf{i}}\cdot\hat{\mathbf{n}}$,$\hat{\mathbf{j}}\cdot\hat{\mathbf{n}}$ or $\hat{\mathbf{k}}\cdot\hat{\mathbf{n}}$. Using this, the other two base vectors of the tangential coordinate system are given by $$ \begin{aligned} \hat{\mathbf{t}}_1 &= \hat{\mathbf{n}}\times\hat{\mathbf{i}}\\ \hat{\mathbf{t}}_2 &= \hat{\mathbf{n}}\times\hat{\mathbf{t}}_1 \end{aligned} $$ Choosing one of the $\hat{\mathbf{i}},\hat{\mathbf{j}},\hat{\mathbf{k}}$ is supposed to prevent the random rotation around $\hat{\mathbf{n}}$ of the vectors $\hat{\mathbf{t}}_1$ and $\hat{\mathbf{t}}_2$(actually it is just a guess).
To transform between the frames following transformation matrix can be used.
$$
\mathbf{R}=\left (
\begin{matrix}
\hat{\mathbf{i}}\cdot \hat{\mathbf{t}}_1 & \hat{\mathbf{j}}\cdot \hat{\mathbf{t}}_1 & \hat{\mathbf{k}}\cdot \hat{\mathbf{t}}_1 \\
\hat{\mathbf{i}}\cdot \hat{\mathbf{t}}_2 & \hat{\mathbf{j}}\cdot \hat{\mathbf{t}}_2 & \hat{\mathbf{k}}\cdot \hat{\mathbf{t}}_2 \\
\hat{\mathbf{i}}\cdot \hat{\mathbf{n}} & \hat{\mathbf{j}}\cdot \hat{\mathbf{n}} & \hat{\mathbf{k}}\cdot \hat{\mathbf{n}}
\end{matrix}
\right )
$$
This transformation matrix between the base vectors of the old and the new frame can be derived by noting that the vectors in the new frame $\hat{\mathbf{n}}$, $\hat{\mathbf{t}}_1$, $\hat{\mathbf{t}}_2$ can be expressed in terms of the vectors in the old frame $\hat{\mathbf{i}},\hat{\mathbf{j}},\hat{\mathbf{k}}$
$$
\begin{matrix}
\hat{\mathbf{t}}_1 = R_{11}\hat{\mathbf{i}} &+& R_{21}\hat{\mathbf{j}} &+& R_{31}\hat{\mathbf{k}} \\
\hat{\mathbf{t}}_2 = R_{12}\hat{\mathbf{i}} &+& R_{22}\hat{\mathbf{j}} &+& R_{32}\hat{\mathbf{k}} \\
\hat{\mathbf{n}} = R_{13}\hat{\mathbf{i}} &+& R_{23}\hat{\mathbf{j}} &+& R_{33}\hat{\mathbf{k}}
\end{matrix}
$$
From the above the entries of the matrix can be found, for instance $R_{11}=\hat{\mathbf{i}}\cdot\hat{\mathbf{t}}_1$
Subsequently, the coordinates of all points $P_i$ in the new frame can be computed by using the formula
$$
(\tilde{x}_i,\tilde{y}_i,\tilde{z}_i)=\mathbf{R}^{-1}\cdot(x_i,y_i,z_i)
$$
This can be derived by looking at the expansion of the vector $\mathbf{x}$ in the old base and the new base and comparing (Einstein summation is used). To keep the notation terse, lets use $\hat{\mathbf{e'}}_i$ for the new base and $\hat{\mathbf{e}}_i$ for the old. $$ \mathbf{x}=x'^j\hat{\mathbf{e'}}_j=x'^iR^j_i\hat{\mathbf{e}}_j=x^j\hat{\mathbf{e}}_j $$
Therefore, by comparing $$ x'^iR^j_i=x^j. $$
Where $\mathbf{R}^{-1}$ is the inverse matrix. Since all points are on the same plane, the resulting set $(\tilde{x}_i,\tilde{y}_i,\tilde{z}_i)$ is in 2D.
Now I would like to do the following. Start with a rectangle located in the plane $z=0$ with $x\in(-1,1)$ and $y\in(-1,1)$ (not roated) and project it on the original plane defined by the unit vectors $\hat{\mathbf{t}}_1$ and $\hat{\mathbf{t}}_2$ and $\hat{\mathbf{n}}$.
This should be done by computing $\mathbf{R}$ using the base vectors case of a non rotted rectangle $$ \begin{aligned} \hat{\mathbf{i}} &= (1,0,0)\\ \hat{\mathbf{j}} &= (0,1,0)\\ \hat{\mathbf{k}} &= (0,0,1)\\ \end{aligned} $$
Unfortunately, the rectangle gets rotated when on the surface. Where I expect something like this:
I get:
Question
Therefore, I would like to ask the experts how to do this rotation. I basically need to figure out a way to "glue/ align" the rectangle with $\hat{\mathbf{i}},\hat{\mathbf{j}},\hat{\mathbf{k}}$ with the surface defined by $\hat{\mathbf{n}}$, $\hat{\mathbf{t}}_1$ and $\hat{\mathbf{t}}_2$.
Running some trials, I noted that to get the desired transformation, I should use the inverse $\mathbf{R}^{-1}$ to rotate $\hat{\mathbf{i}},\hat{\mathbf{j}},\hat{\mathbf{k}}$, which comes a bit surprising.
First you have to project the rectangle onto the plane, and then write the projection in the new coordinates.
If $\mathbf{p}$ is a point in $\Bbb{R}^3$, then the projection of $\mathbf{p}$ onto your plane $\hat{\mathbf{n}}\cdot (\mathbf{r}-\mathbf{r}_0)$ is given by $$\mathbf{q} = \mathbf{p} - (\hat{\mathbf{n}}\cdot (\mathbf{p}-\mathbf{r}_0))\hat{\mathbf{n}}.$$ So for for your rectangle $P= [-1,1]^2\times \{0\}$ the projection is $$Q = \{\mathbf{p} - (\hat{\mathbf{n}}\cdot (\mathbf{p}-\mathbf{r}_0))\hat{\mathbf{n}} : \mathbf{p} \in [-1,1]^2\times \{0\}\}.$$ Now we have to write $\mathbf{q} \in Q$ in your new coordinates $(\hat{x},\hat{y},\hat{z})$ w.r.t. axes $\hat{\mathbf{n}}, \hat{\mathbf{t}}_1, \hat{\mathbf{t}}_2$. Since $\hat{\mathbf{n}}$ corresponds to $(1,0,0)$, we have $$(\hat{x},\hat{y},\hat{z}) = R^{-1}\mathbf{p} - (\hat{\mathbf{n}}\cdot (\mathbf{p}-\mathbf{r}_0),0,0).$$