A common way to do this is to find a planar perspective transformation that “warps” the rectangle into the image quadrilateral and then use its inverse to map points on the image back to the rectangle. There are software libraries that include this as standard functionality.
If you want to code this up for yourself, it’s not terribly difficult to work out the necessary mapping. Without going into the detailed derivation, a general planar perspective transformation can be represented (in homogeneous coordinates) by a matrix of the form $$
M=\pmatrix{m_{00} & m_{01} & m_{02} \\ m_{10} & m_{11} & m_{12} \\ m_{20} & m_{21} & 1},
$$ which corresponds to the mapping $$\begin{align}
x' &= {m_{00}x+m_{01}y+m_{02} \over m_{20}+m_{21}+1} \\
y' &= {m_{10}x+m_{11}y+m_{12} \over m_{20}+m_{21}+1}.
\end{align}$$ Given a set of corner-to-corner correspondences between a pair of quadrilaterals, the eight coefficients can by found by solving a system of linear equations.
If one of the quads is a rectangle aligned with the coordinate axes, this transformation can be built up in stages, which can be more convenient for implementation in software. Assuming that we have the matrix $A$ which maps from the unit square to the quadrilateral, the rectangle-to-quadrilateral transformation can be derived by composing it with suitable translation and scaling transformations: $$
M = A\cdot S\left(\frac1w,\frac1h\right)\cdot T(-x_{LL},-y_{LL}).
$$ The inverse map is then $$
M^{-1}=T(x_{LL},y_{LL})\cdot S(w,h)\cdot A^{-1} = \pmatrix{w&0&x_{LL} \\ 0&h&y_{LL} \\ 0&0&1 }\cdot A^{-1},
$$ where $(x_{LL},y_{LL})$ are the coordinates of the rectangle’s lower-left corner. Left-handed coordinate systems can be accommodated by throwing in the appropriate reflections, and skew rectangles can be dealt with by adding a rotation to the transformation cascade.
Now, it’s just a matter of finding the matrix $A$. Let the corners of the quadrilateral be given by the points $q_i'$, and map the corners of the unit square to them as follows: $$\begin{align}
(0,0)&\mapsto q_0' \\
(1,0)&\mapsto q_1' \\
(0,1)&\mapsto q_2' \\
(1,1)&\mapsto q_3'.
\end{align}$$ Solving the resulting system of equations is tedious, but straightforward. One form of the solution is as follows: $$\begin{align}
a_{00} &= a_{20}x_1'+\Delta x_{10} \\
a_{10} &= a_{20}y_1'+\Delta y_{10} \\
a_{20} &= {\omega(\Delta_{10},\Delta_{32})+\omega(\Delta_{20},\Delta_{32})-\omega(\Delta_{30},\Delta_{32}) \over \omega(\Delta_{31},\Delta_{32})} \\
a_{01} &= a_{21}x_2'+\Delta x_{20} \\
a_{11} &= a_{21}y_2'+\Delta y_{20} \\
a_{21} &= -{\omega(\Delta_{10},\Delta_{31})+\omega(\Delta_{20},\Delta_{31})-\omega(\Delta_{30},\Delta_{31}) \over \omega(\Delta_{31},\Delta_{32})} \\
a_{02} &= x_0' \\
a_{12} &= y_0', \\
\end{align}$$ where $\Delta_{ij}=q_i'-q_j'$, $\Delta x_{ij}$ and $\Delta y_{ij}$ are the corresponding coordinate deltas, and $\omega$ is the symplectic form $\omega(\mathbf u, \mathbf v) = u_xv_y-u_yv_x$. (That this looks like a lot of cross products is no coincidence—there’s a geometric interpretation of the solution as the intersection of various planes.)
Since you’re mainly interested in mapping from a quadrilateral to a rectangle, it might be more efficient, and probably more computationally stable, to compute $A^{-1}$ directly from the corner coordinates. It’s similar to the solution for the other direction, but I’ll leave that calculation to you.
Best Answer
For a regular simplex in $d-$space, look at the centroid $G$ of an arbitrary $d-2$ face $s_{d-2}$. Furthermore, there is a pair of vertices, call them $V_i$ and $V_j$, which are from the total $d$ simplex, but not from the $d-2$ face $s_{d-2}$ (they are what's left of the vertices of the $d$ simplex when you remove all the vertices of the $d-2$ face $s_{n-2}$). This centroid is the center of the $d-2$ face $s_{d-2}$, which itself is a $d-2$ regular simplex obtained as the intersection of the pair of adjacent $d-1$ regular simplicial faces $V_i\,s_{d-2}$ and $V_j\,s_{d-2}$ of the $d$ simplex. Because the $d$ simplex is regular, all of its faces are also regular simplices of various dimensions. Hence, each segment $A_iG$ and $A_jG$ is perpendicular to the $d-2$ affine subspace spanned by the $d-2$ simplex $s_{d-2}$. Moreover, these two segments are equal in length, i.e. $V_iG = V_jG = h_{d-1}$ because they are the heights of two congruent regular simplices of dimension $d-1$, namely $V_is_{d-2}$ and $V_js_{d-2}$. Therefore, the triangle $V_iV_jG$ is isosceles, where $V_iG = V_jG = h_{d-1}$ and $V_iV_j=1$. The dihedral angle you are interested in is then $$\angle \, V_iGV_j \, = \, \theta$$
If you denote by $G_i$ and $G_j$ the centroids (centers) of the $d-1$ simplices
$V_i \, s_{d-2}$ and $V_j\,s_{d-2}$ respectively, then $G_i \in V_iG$ and $G_j \in V_jG$, as well as $V_iG_j \, \perp \, V_jG$ and $V_jG_i \, \perp \, V_iG$. The segments $V_iG_j$ and $V_jG_i$ are congruent heights in the isosceles triangle $V_iV_jG$ and the triangle $V_iG_jG$ is right-angled, so $$\cos(\theta) \, =\, \cos(\angle \, V_iGV_j) \, =\, \cos(\angle \, V_iGG_j) \, =\, \frac{G_jG}{V_iG}$$ and by the reflection symmetry of the isosceles triangle $V_iV_jG$, the two segments $G_iG = G_jG$ so $$\cos(\theta) \, =\, \frac{G_iG}{V_iG}$$ Observe that in the latter equality, $V_iG$ is the height of the regular $d-1$ simplex $V_i\,s_{d-2}$ that connects the vertex $V_i$ with the center $G$ of the opposing regular $d-2$ simplex $s_{d-2}$. Moreover, the point $G_i$ is the center of the $d-1$ simplex $V_i \, s_{d-2}$.
Now everything boils down to calculating the ratio in which the center $G_i$ of a regular $d-1$ simplex $V_i\,s_{d-2}$ splits the height $V_iG$ connecting $V_i$ with the center of the regular $d-2$ simplex $s_{d-2}$:
We will do this by showing how it works in one dimension higher ($d-1$ regular face of a $d$ regular simplex, instead of the current $d-2$ regular face of a $d-1$ regular simplex) and use it as an induction step in $d$, calculating the ratio in which the center $G_d$ of the total regular simplex $s_{d}$ splits the height $V_iG_j$. In the current notations, $G_d = V_iG_j \, \cap \, V_jG_i$. The triangles $V_iG_iG_d$ and $V_iG_jG$ are right-angled and similar, so $$\frac{G_jG}{V_iG} = \frac{G_iG_d}{V_iG_d} = \cos(\theta)$$ But by the fact that the traigle $V_iV_jG$ is isosceles, $$G_jG = G_iG \,\,\,\, \text{ and } \,\,\,\, G_iG_d = G_jG_d$$ which yields the ratio $$\frac{G_iG}{V_iG} = \frac{G_jG_d}{V_iG_d} = \cos(\theta)$$
Let me flip it
$$\frac{V_iG_d}{G_dG_j} = \frac{V_iG}{GG_i} = \frac{1}{\cos(\theta)}$$
I claim that $$\frac{V_iG_d}{G_dG_j} = d$$
The ratio $\frac{V_iG}{GG_i}$ is a ratio that comes from the analogous situation one dimension lower, i.e. the position of the canter $G_i$ of the $d-1$ regular simplex $V_i \, s_{d-2}$ on the height $V_iG$ of this same simplex. If, by induction, we assume that it has been shown that $$\frac{V_iG_i}{G_iG} = d-1$$ then $$\frac{V_iG}{GG_i} = \frac{V_iG_i + G_iG}{G_iG} = \frac{V_iG_i}{G_iG} + 1 = d - 1 + 1 = d$$ Hence
$$\frac{V_iG_d}{G_dG_j} = \frac{V_iG}{GG_i} = d = \frac{1}{\cos(\theta)}$$ and there you have that
$$\cos(\theta) = \frac{1}{d}$$ and consequently the dihedral angle you are after is
$$\theta = \arccos\left(\frac{1}{d}\right)$$
The beginning of the inductive step can start with $d=2$, in which case we have that we have an equilateral triangle, where the ratio is known to be $\frac{1}{2}$.