Cartesian formula for incenter of tangential quadrilateral from vertices only (no angles) I'm looking for the Cartesian formula to express the coordinates of the incenter of a tangential quadrilateral based on the coordinates of the vertices (without angles). I couldn't find such a formula anywhere. I was only able to find the general formula for the incenter of a triangle. The algebra also seems to be very tedious and I don't have the patience to derive it in pen and paper. If that helps, I need the formula for programming purposes. So, if anyone knows where I can find it, that would be great. Thank you!
Edit: Thank you all for coming up with these great solutions and describing them in such depth. I really appreciate it!
 A: Edit : in fact the solution I presented at first is not simple ; I keep it as "Solution n°2" in order to understand in particular the very interesting comment of Blue.
Solution n°1 (the simplest):

(this figure is adapted from this site).
Idea: the center is plainly the incenter (center of incircle) of triangle $BDE$ where $E$ is the intersection point of $AB$ and $CD$ (assumed non parallel : otherwise take the other pair of opposite sides, of course if $ABCD$ isn't a square).
Details of computation : The coordinates of $E$ are given by considering that there exists $\lambda, \ \mu$ such that:
$$E=B+\lambda(A-B)=D+\mu(C-D)\tag{a}$$
Implying $$\lambda(A-B)+\mu(D-C)=D-B\tag{b}$$
giving rise to the linear system:
$$\begin{cases}\lambda(x_A-x_B)&+&\mu(x_D-x_C)&=&x_D-x_B\\\lambda(y_A-y_B)&+&\mu(y_D-y_C)&=&y_D-y_B\end{cases}$$
Solving it, we get in particular a formula for $\lambda$:
$$\lambda=\dfrac{(x_D-x_B)(y_D-y_C)-(y_D-y_B)(x_D-x_C)}{(x_A-x_B)(y_D-y_C)-(y_A-y_B)(x_D-x_C)}$$
that we just have to "plug" into (a), giving:
$$\begin{cases}x_E&=&x_B+\lambda (x_A-x_B)\\y_E&=&y_B+\lambda (y_A-y_B)\end{cases}$$
We are now able to deduce the coordinates of incenter $I=k(aA + bB +eE) \ \text{with} \  k:=\dfrac{1}{a+b+c}$:
$$\begin{cases}x_I&=&k(ax_A+bx_B+ex_E)\\y_I&=&k(ay_A+by_B+ey_E)\end{cases}$$
i.e., the weighted average of the coordinates of $A,B,E$ with weights the lengths $a=BE=\sqrt{(x_E-x_B)^2+(y_E-y_B)^2}$, $b=AE$, e=$AB$.
Solution n°2 (less simple):
It is based on Newton's theorem as described and proved in the reference given above.
Let:
$$\begin{cases}x_E&=&\tfrac12(x_A+x_C)\\ y_E&=&\tfrac12(y_A+y_C)\end{cases} \ \ 
   \ \begin{cases}x_F&=&\tfrac12(x_B+x_D)\\ y_F&=&\tfrac12(y_B+y_D)\end{cases}$$
the coordinates of center $O$ are :
$$\begin{cases}x_O&=&\lambda x_E  + (1-\lambda)x_F\\y_O&=&\lambda y_E  + (1-\lambda)y_F \end{cases}\tag{1}$$
for a certain $\lambda, 0 \le \lambda \le 1$.
Besides, consider the following double formula that can be found in this article
$$K=r.s=\frac12 \sqrt{p^2q^2-(ac-bd)^2}\tag{2}$$
where $K$ is the area of the quadrilateral, $a,b,c,d$ its side lengths, $p,q$ its diagonal lengths, $r$ the radius of the inscribed circle, $s=\frac12(a+b+c+d)$ the semi-perimeter.
From (2), radius $r$ can be expressed as a function of $a,b,c,d,p,q$.
It remains to constrain the distance from $O$ to line $AB$ to be equal to $r$. For this we need the equation of straight line $AB$, which is :
$$\begin{vmatrix}x&x_A&x_B\\y&y_A&y_B\\1&1&1\end{vmatrix}=0$$
or, in an expanded form:
$$ux+vy+w=0 \ \text{with} \ \begin{cases}u&=&y_A-y_B\\v&=&x_B-x_A\\w&=&x_Ay_B-x_By_A\end{cases}$$
then the distance from $O$ to line $AB$ is given by a classical formula:
$$r=\dfrac{1}{\sqrt{u^2+v^2}}|ux_0+vy_0+w|\tag{3}$$
(3) is a first degree equation in $\lambda$ from which one can deduce its value (in fact one has to choose between a plus and a minus sign for the absolute value).
A: Let us derive a formula for incenter valid for all tangential polygons.
We will identify Euclidean plane $\mathbb{R}^2$ with complex plane $\mathbb{C}$ and use same symbol to refer to a geometric point and corresponding complex number.
Let $p_1,p_2,\ldots,p_n$ be vertices of a tangential $n$-gon with inradius $r$ and incenter $c$. We will assume $p_k$ are ordered
in a counterclockwise manner along the circumference of the $n$-gon.
Extend labeling of indices by periodicity (i.e $p_0 = p_n, p_1 = p_{n+1},\ldots$, etc.)
Let $\psi_k$ be the interior angle $\angle p_{k-1}p_k p_{k+1}$.
Let $q_k = c + re^{i\theta_k}$ be the contact point of edge $p_kp_{k+1}$ with the incircle.
In terms of $c$ and $\theta_k$, we have
$$\psi_k = \pi - (\theta_k - \theta_{k-1})
\quad\text{ and }\quad
p_k = c + re^{\frac{i}{2}(\theta_k + \theta_{k-1})}\sec\left(\frac{\theta_k - \theta_{k-1}}{2}\right)
$$
This leads to
$$(p_k - c)\sin\psi_k = 2re^{\frac{i}{2}(\theta_k + \theta_{k-1})}
\sin\left(\frac{\theta_k - \theta_{k-1}}{2}\right)
= -i(re^{i\theta_k} - re^{i\theta_{k-1}}) = -i(q_k - q_{k-1})
$$
Summing over $k$ gives us $\sum\limits_{k=1}^n(p_k - c)\sin\psi_k = 0$. This leads to

The incenter $c$ of any tangential $n$-gon is a weighted sum of it vertices with weight proportional to $\sin\psi_k$: $$c = \frac{\sum\limits_{k=1}^n p_k \sin\psi_k}{\sum\limits_{k=1}^n \sin\psi_k}$$

Let $\Delta_k$ be the area of triangle $p_{k-1}p_kp_{k+1}$ and $\ell_k$ be the length of edge $p_kp_{k+1}$. If you know the coordinates of $p_k = (x_k,y_k)$, you can compute these by the formula:
$$2\Delta_k = \left|\begin{matrix}
x_{k-1} & y_{k-1}  & 1\\
x_{k} & y_k & 1 \\
x_{k+1} & y_{k+1} & 1\end{matrix}\right|\quad\text{ and }\quad
\ell_k = \sqrt{(x_k - x_{k+1})^2+(y_k - y_{k+1})^2}
$$
This will allow you to compute the weight $\sin\psi_k$ by the formula:
$$\sin\psi_k = \frac{2\Delta_k}{\ell_{k-1}\ell_k}$$
