Computing a projective transformation
A projective transformation of the (projective) plane is uniquely defined by four projected points, unless three of them are collinear. Here is how you can obtain the $3\times 3$ transformation matrix of the projective transformation.
Step 1: Starting with the 4 positions in the source image, named $(x_1,y_1)$ through $(x_4,y_4)$, you solve the following system of linear equations:
$$\begin{pmatrix}
x_1 & x_2 & x_3 \\
y_1 & y_2 & y_3 \\
1 & 1 & 1
\end{pmatrix}\cdot
\begin{pmatrix}\lambda\\\mu\\\tau\end{pmatrix}=
\begin{pmatrix}x_4\\y_4\\1\end{pmatrix}$$
The colums form homogenous coordinates: one dimension more, created by adding a $1$ as the last entry. In subsequent steps, multiples of these vectors will be used to denote the same points. See the last step for an example of how to turn these back into two-dimensional coordinates.
Step 2: Scale the columns by the coefficients you just computed:
$$A=\left(\begin{array}{lll}
\lambda\cdot x_1 & \mu\cdot x_2 & \tau\cdot x_3 \\
\lambda\cdot y_1 & \mu\cdot y_2 & \tau\cdot y_3 \\
\lambda & \mu & \tau
\end{array}\right)$$
This matrix will map $(1,0,0)$ to a multiple of $(x_1,y_1,1)$, $(0,1,0)$ to a multiple of $(x_2,y_2,1)$, $(0,0,1)$ to a multiple of $(x_3,y_3,1)$ and $(1,1,1)$ to $(x_4,y_4,1)$. So it will map these four special vectors (called basis vectors in subsequent explanations) to the specified positions in the image.
Step 3: Repeat steps 1 and 2 for the corresponding positions in the destination image, in order to obtain a second matrix called $B$.
This is a map from basis vectors to destination positions.
Step 4: Invert $A$ to obtain $A^{-1}$ (or use the adjugate as discussed below).
$A$ maps from basis vectors to the source positions, so the inverse matrix maps in the reverse direction.
Step 5: Compute the combined Matrix $C = B\cdot A^{-1}$.
$A^{-1}$ maps from source positions to basis vectors, while $B$ maps from there to destination positions. So the combination maps source positions to destination positions. This is the matrix of the transformation you were requesting.
Step 6: To map a location $(x,y)$ from the source image to its corresponding location in the destination image, compute the product
$$\begin{pmatrix}x'\\y'\\z'\end{pmatrix} =
C\cdot\begin{pmatrix}x\\y\\1\end{pmatrix}$$
These are the homogenous coordinates of your transformed point.
Step 7: Compute the position in the destination image like this:
\begin{align*}
x'' &= \frac{x'}{z'} \\
y'' &= \frac{y'}{z'}
\end{align*}
This is called dehomogenization of the coordinate vector.
How to use this projective transformation with CSS
In general such a transformation will not be an affine transformation, so you cannot express this in terms of affine transformations like scaling, rotating and shearing, since these cannot express perspectivity. You might however try to simply set the first two entries of the last row to zero, so you get an affine transformation which might be close enough to your desired transformation.
If on the other hand you can use a matrix3d transformation, then you can take the 2D projective transformation matrix $C$ computed as described above, and use its entries to build a 3D projective transformation matrix like this:
$$\begin{pmatrix}
C_{1,1} & C_{1,2} & 0 & C_{1,3} \\
C_{2,1} & C_{2,2} & 0 & C_{2,3} \\
0 & 0 & 1 & 0 \\
C_{3,1} & C_{3,2} & 0 & C_{3,3}
\end{pmatrix}$$
This transformation will transform $x$ and $y$ coordinate as above, but leave the $z$ coordinates of the homogenous coordinate vectors alone. Dehomogenization might still change the value of the $z$ coordinate in space, but as you don't really care about these, this should be good enough.
I've written a proof-of-concept implementation. The user interface is pretty crude, but the math works well enough. The implementation there uses the adjugate matrix instead of the inverse, both in solving the linear equations in step 1 and for the reverse transform in step 4. The result differs only by a scalar factor, which is irrelevant for homogenous coordinates. The benefit is that this avoids computing a bunch of determinants and performing a bunch of divisions.
If you wanted to, you could play the same game for five points in 3D space, in order to compute the full spatial projective transformation matrix. But that only makes sense if you actually have depth, sice no four of the five points may be coplanar.
