Let $Q$ be a convex quadrilateral in $\mathbb{R}^2$ with vertices $a_1,a_2,a_3,a_4 \in \mathbb{R}^2$. Consider the bilinear map $f \colon [0,1]^2 \to Q$, $$f(x,y)=a_1+(a_2−a_1)x+(a_4−a_1)y+(a_1+a_3−a_2−a_4)x y$$
Note that $f$ is defined in such a way that it maps the vertices of $[0,1]^2$ to the vertices of $Q$. To prove: $$f([0,1]^2)=Q.$$
Remark: the implication $\subseteq$ is clear since for any $(x,y)\in [0,1]^2$ $$f(x,y)=(1-x-y+xy)a_1+(x-xy)a_2+xya_3+(y-xy)a_4 \in \text{conv}\{a_1,a_2,a_3,a_4\}=Q.$$ How to prove the other direction?
Short summary: you have rewritten $f$ as convex combination of the vertices. You just have to solve for $x$ and $y$ given the coefficients.
We have to show that given $q \in Q$ there exist $x, y \in [0, 1]$ such that $f(x, y) = q$.
As $q \in Q = \text{conv}(a_1, \ldots, a_4)$, there exist $\lambda_1, \ldots, \lambda_4 \ge 0$ with $\sum_{k = 1}^{4} \lambda_k = 1$ such that $$ q = \sum_{k = 1}^{4} \lambda_k a_k \overset{!}{=} (1 - x - y + x y) a_1 + x(1 - y) a_2 + x y a_3 + y(1 - x) a_4. $$ Hence we have to show that there exists a solution $(x, y) \in [0, 1]^2$ to $$ \begin{cases} 1 - x - y + x y &= \lambda_1, \\ x(1 - y) &= \lambda_2, \\ x y &= \lambda_3, \\ y(1 - x) &= \lambda_4. \end{cases} $$ Adding the first two equations gives $\lambda_1 + \lambda_2 = 1 - y$, that is, $y = 1 - \lambda_1 - \lambda_2$. Adding the first and last equation gives $\lambda_1 + \lambda_4 = 1 - x$ and thus $x = 1 - \lambda_1 - \lambda_4$. (Alternatively, adding the second and third equation gives $\lambda_2 + \lambda_3 = x$ and adding the third and last equation gives $\lambda_3 + \lambda_4 = y$, which, together with $\sum_{k = 1}^{4} \lambda_k = 1$ gives the same result.)
Since $0 \le \lambda_1 + \lambda_2 \le \sum_{k = 1}^{4} \lambda_k = 1$ and $0 \le \lambda_1 + \lambda_4 \le \sum_{k = 1}^{4} \lambda_k = 1$, we have $x, y \in [0, 1]$.