Let $C(c_x,c_y)$. Then, we have two equations about $c_x,c_y$ :
$$\vec{BA}\cdot\vec{BC}=0\iff (a_x-b_x)(c_x-b_x)+(a_y-b_y)(c_y-b_y)=0$$$$BC=L\iff (c_x-b_x)^2+(c_y-b_y)^2=L^2$$
If you have $c_x\lt b_x$, then you can get $c_x,c_y$ by solving these.
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
Let the coordinates of $p_n$ be $(x_n,y_n)$. Then the slope of $A$ is $m_A=\frac{y_2-y_1}{x_2-x_1}$. The slope of $B$ is $m_B=\frac{-1}{m_A}=\frac{x_1-x_2}{y_2-y_1}$. Then $p_3=p_1\pm B(\frac{1}{\sqrt{1+m_B^2}},\frac{m_B}{\sqrt{1+m_B^2}})$ where the sign ambiguity corresponds to two orientations of the triangle. I have ignored issues when the sides are vertical or horizontal, which can lead to division by zero