Welcome to M.S.E.
This is an amazing problem! You're going right on creating the incircle, tangents, and the altitudes. As a hint, let the four angles at the vertices be equal to $2\alpha,2\beta,2\gamma,2\delta,$ respectively. Let $O$ be the incenter. Then the lines $AO,BO, CO,$ and $DO$ would all bisect the angles of the quadrilateral, and the divided angles will have magnitudes of $\alpha, \beta, \gamma,$ and $\delta.$
Then try to apply trigonometry on the individual triangles formed by the tangency. You'll get;
$$\tan{\alpha}=\frac{r}{a},\tan{\beta}=\frac{r}{b},\text{ and } \tan{\gamma}=\frac{r}{c}.$$
But then what is $\tan{(\alpha+\beta+\gamma+\delta)}$ equal to? You can use the tangent addition formula to make a cool observation;
$$\dfrac{\tan(\alpha+\beta)+\tan(\gamma+\delta)}{1-\tan(\alpha+\beta)\tan(\alpha+\beta)}=\tan{(\alpha+\beta+\gamma+\delta)}.$$
It will be easy after this. I'll leave it to you to fill in the blanks. Best of Luck!
Let $I$ denote the center of the circle inscribed in quadrilateral $ABCD,$ (i.e. the incenter of quadrilateral $ABCD$). Also let $H_{AB},$ $H_{BC},$ $H_{CD},$ and $H_{DA}$ denote the foot of perpendiculars from $I$ to sides $AB,$ $BC,$ $CD,$ and $DA,$ respectively. We know that the inradius of quadrilateral $ABCD$ is$$IH_{AB}=IH_{BC}=IH_{CD}=IH_{DA}.$$Since $AH_{AB}=AH_{DA}=a,$ $IH_{AB}=IH_{BC},$ and $\angle AH_{AB}I=\angle AH_{DA}I=90^\circ,$ we know that $\triangle AH_{AB}I\cong\triangle AH_{DA}I$ by the SAS congruence theorem. Hence, we have;$$\angle H_{DA}AI=\angle H_{AB}AI.$$Similarly, we can show that $IA,$ $IB,$ $IC,$ and $ID$ bisect $\angle DAB,$ $\angle ABC,$ $\angle BCD,$ and $\angle CDA,$ respectively. We let the angles of quadrilateral $ABCD$ be $2\alpha,$ $2\beta,$ $2\gamma,$ and $2\delta;$ $$\angle DAB=2\alpha$$ $$\angle ABC=2\beta$$ $$\angle BCD=2\gamma$$ $$\angle CDA=2\delta.$$ 
Since $IA,$ $IB,$ $IC,$ and $ID$ bisect $\angle DAB,$ $\angle ABC,$ $\angle BCD,$ and $\angle CDA,$ respectively, we know that;$$\angle DAI=\angle IAB=\alpha$$$$\angle ABI=\angle IBC=\beta$$$$\angle BCI=\angle ICD=\gamma$$$$\angle CDI=\angle IDA=\delta.$$We know that $\angle DAB+\angle ABC+\angle BCD+\angle CDA=360^\circ,$ hence, we have;$$2\alpha+2\beta+2\gamma+2\delta=360^\circ$$$$\implies \alpha+\beta+\gamma+\delta=180^\circ.$$$$\implies \tan(\alpha+\beta+\gamma+\delta)=\tan{180^\circ}=0.$$Now we use the Tangent Sum Identity, i.e. for any two angles $x$ and $y$;$$\tan(x+y) = \frac{\tan x + \tan y}{1-\tan x \tan y}.$$We have;$$\tan(\alpha+\beta+\gamma+\delta)=\frac{\tan(\alpha+\beta)+\tan(\gamma+\delta)}{1-\tan(\alpha+\beta)\tan(\gamma+\delta)}=0.$$Hence, we know that;$$\tan(\alpha+\beta)+\tan(\gamma+\delta)=0$$Applying the Tangent Sum Identity again gives us;$$\frac{\tan \alpha + \tan \beta}{1-\tan \alpha \tan \beta}+\frac{\tan \gamma + \tan \delta}{1-\tan \gamma \tan \delta}=0\tag{Equation 1}$$Analogously, we also have;$$\frac{\tan \alpha + \tan \delta}{1-\tan \alpha \tan \delta}+\frac{\tan \beta + \tan \gamma}{1-\tan \beta\tan \gamma}=0,\tag{Equation 2}$$and$$\frac{\tan \alpha + \tan \gamma}{1-\tan \alpha \tan \gamma}+\frac{\tan \beta + \tan \delta}{1-\tan \beta\tan \delta}=0.\tag{Equation 3}$$We know that;$$\tan\alpha=\frac{r}{a}$$$$\tan\beta=\frac{r}{b}$$$$\tan\gamma=\frac{r}{c}$$$$\tan\delta=\frac{r}{d}.$$Hence, we have;$$r=a\cdot\tan\alpha=b\cdot\tan\beta=c\cdot\tan\gamma=d\cdot\tan\delta.$$We know that;$$\tan\alpha+\tan\beta=\frac{r}{a}+\frac{r}{b}=\frac{r(a+b)}{ab}.$$and,$$\tan\alpha\tan\beta=\frac{r}{a}\cdot\frac{r}{b}=\frac{r^2}{ab}.$$Hence, we have;\begin{align*}\dfrac{\tan \alpha + \tan \beta}{1-\tan \alpha \tan \beta}&=\dfrac{\frac{r(a+b)}{ab}}{1-\frac{r^2}{ab}}\\&=\dfrac{\frac{r(a+b)}{ab}}{\frac{ab-r^2}{ab}}\\&=\dfrac{r(a+b)}{ab-r^2}.\end{align*}Similarly, we also have;$$\frac{\tan \gamma + \tan \delta}{1-\tan \gamma \tan \delta}=\dfrac{r(c+d)}{cd-r^2},$$$$\frac{\tan \alpha + \tan \delta}{1-\tan \alpha \tan \delta}=\dfrac{r(a+d)}{ad-r^2},$$and$$\frac{\tan \beta + \tan \gamma}{1-\tan \beta\tan \gamma}=\dfrac{r(b+c)}{bc-r^2}.$$Now, we can substitute these values into Equation 1, and Equation 2, we have;$$\frac{\tan \alpha + \tan \beta}{1-\tan \alpha \tan \beta}+\frac{\tan \gamma + \tan \delta}{1-\tan \gamma \tan \delta}=\dfrac{r(a+b)}{ab-r^2}+\dfrac{r(c+d)}{cd-r^2}=0$$and$$\frac{\tan \alpha + \tan \delta}{1-\tan \alpha \tan \delta}+\frac{\tan \beta + \tan \gamma}{1-\tan \beta\tan \gamma}=\dfrac{r(a+d)}{ad-r^2}+\dfrac{r(b+c)}{bc-r^2}=0.$$Since in these two Equations, the two values add up to 0, we can reciprocate each value and yet keep their sum the same, we have;$$\dfrac{ab-r^2}{r(a+b)}+\dfrac{cd-r^2}{r(c+d)}=0$$and$$\dfrac{ad-r^2}{r(a+d)}+\dfrac{bc-r^2)}{r(b+c)}=0.$$Simplifying these Equations gives;$$\begin{align*} & \dfrac{(ab-r^2)(c+d)}{r(a+b)(c+d)}+\dfrac{(cd-r^2)(a+b)}{r(a+b)(c+d)}=0\\ \implies & \dfrac{(ab-r^2)(c+d)+(cd-r^2)(a+b)}{r(a+b)(c+d)}=0\\ \implies & (ab-r^2)(c+d)+(cd-r^2)(a+b)=0\\ \implies & abc+abd-(c+d)r^2+acd+bcd-(a+b)r^2=0 \end{align*}$$$$\begin{align*} \implies abc+abd+acd+bcd & =(a+b)r^2+(c+d)r^2 \\ & =(a+b+c+d)r^2 \end{align*}$$$$\implies r^2=\frac{abc+abd+acd+bcd}{a+b+c+d}.\tag{Equation 4}$$We know that;$$\begin{align*}[ABCD]&=[AIB]+[BIC]+[CID]+[DIA]\\&=\frac{1}{2}AB\cdot r+\frac{1}{2}BC\cdot r+\frac{1}{2}CD\cdot r+\frac{1}{2}DA\cdot r\\&=\frac{1}{2}(AB+BC+CD+DA)\cdot r\\&=\frac{AB+BC+CD+DA}{2}\cdot r\\&=\frac{a+b+b+c+c+d+d+a}{2}\cdot r\\&=(a+b+c+d)\cdot r.\end{align*}$$Hence, we have;$$[ABCD]^2=r^2(a+b+c+d)^2.$$By Equation 4, we know that;$$\implies r^2=\frac{abc+abd+acd+bcd}{a+b+c+d}.$$Substituting this value in Equation 4 gives us;$$\begin{align*}[ABCD]^2&=r^2(a+b+c+d)^2\\&=(a+b+c+d)^2\cdot\frac{abc+abd+acd+bcd}{a+b+c+d}\\&=(a+b+c+d)(abc+abd+acd+bcd),\end{align*}$$as desired.
You may be better off generating a parallelogram within a single quadrilateral. Simply connect the midpoints of the sides in rotational order. The parallelogram thus obtained has half the area of the original quadrilateral, regardless of the shape of the latter figure.
Divide and Conquer
How does that work? Let $Q$ be any quadrilateral and $P$ be the parallelogram defined by the midpoints of the sides of $Q$. Divide the plane in half along a diagonal of $Q$, choosing the diagonal through the concave vertex if such exists.
Then the area of $Q$ is divided into two triangular areas. Taking the base to lie along the dividing line, the area of $Q$ is then measured by half the length of the dividing diagonal tines the sum of the triangle altitudes.
Meanwhile the area of $P$ is divided into two parallelogram areas because the sides of $P$ are constructed parallel to the diagonals of $Q$, and the base of $P$ may also be rendered along the dividing line. Then this base measures half the dividing diagonal inside $Q$ and the corresponding altitude measures half the sum of the triangle altitudes, from which the area of $P$ works out to half the previously rendered area of $Q$.
Best Answer
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).