Given 4 random points in 2D, how do I compute the area of the quadrilateral formed by the points? I'm aware of formulae giving the area when I know the sides a,b,c,d and the diagonals p & q. But how do I decide algorithmically which of the 6 connecting lines between the 4 points are sides and which are diagonals?

Here's a literal answer to your question: Sort the points, so that $a$ is has the leftmost $x$coordinate, $ba$ has the smallest argument (choosing the branch where arguments take values in $(\pi,\pi]$), and $ca$ has the largest argument. Take half the ($z$coordinate of the) crossproduct of $ba$ with $ca$ to get the area of the triangle $abc$. Add either $0$ or half the ($z$coordinate of the) crossproduct of $dc$ and $db$, depending on which is bigger. You might want to start by checking triples for collinearity and throwing out the internal points, to avoid ambiguity. 


Since the value depends on the order of the points, why not calculating the area on all the possible permutations (4!=24) ? or on some subset (the "necklaces"): The largest area value seems to corresponds to nonself crossing quadrilaterals: some python code here 


I guess that you want something like "invariants of ordertypes". An order type is a configuration of a finite number number of points which is generic (no triplets of aligned points) and up to isotopy. An easy invariant of the order type is the number of vertices in its convex hull. This invariant is complete for four points since there are only two possible configurations. Let me mention another invariant which is much less know and also very easy to compute: Given an order type $\mathcal C$, there is a unique equivalence relation of $\mathcal C$ having at most two classes such that the following holds: If $\mathcal C$ is in convex position and consists of an odd number of points, the equivalence relation is trivial. If $\mathcal C$ is in convex position and contains an even number of points, then two consecutive vertices of the convex hull are nonequivalent. If two configurations $\mathcal C_1,\mathcal C_2$ are related by a move which goes through a unique (unavoidable) alignement involving three points $\lbrace P,Q,R\rbrace$, then the equivalence relation $\mathcal C_1$ and $\mathcal C_2$ is the same for $\mathcal C_i\setminus \lbrace P,Q,R\rbrace$. The three points $P,Q,R$ change classes (with respect to a fixed class represented by a point in $\mathcal C_i\setminus \lbrace P,Q,R\rbrace$). This determines the equivalence relation completely since any two generic configurations containing involving $n$ points can be related by isotopies and moves involving $3$point alignements. Concretely, this relation can be computed as follows: Given two points $P,Q\in \mathcal C$ of a generic configuration $\mathcal C$ with $n$ points let $s(P,Q)$ denote the number of lines defined by two points of $\mathcal C\setminus\lbrace P,Q\rbrace$ which separate $P$ and $Q$. Then $P\sim Q$ if and only if $s(P,Q)\equiv n+1\pmod 2$. (This relation has in fact a grouptheoretic interpretation which gives generalizations to arbitrary dimensions. There is even a "fermionic" version which I did not publish.) Application: This gives a rational rule for planting plum and cherrytrees in an orchard (up to transposition of plum and cherry), if the positions of the trees are generic. For four points, you get points with alternating classes in the case of convex position. In the case of a triangle with an interior point, the three vertices of the triangle form one class, the remaining interior point forms the other class. References: MR2038321 (the general case) or MR2300644 and MR2300648 (only the two dimensional babycase). 


A simple way of looking at the problem is the following. You can compute the Area by computing the integral: $A = \int_C d^2x$ where $C$ is yor quadrilateral domain bounded by the 4 corners with coordinates $P_i = (x_i,y_i)$ with $i=1,...,4$. It is easy to check that the transformation $(x,y) \rightarrow (u,v)$ defined by $x = \sum_i x_i Q_i(u,v); y = \sum_i y_i Q_i(u,v)$ with $Q_1(u,v) = \frac{(1u)(1v)}{4}$ $Q_2(u,v) = \frac{(1+u)(1v)}{4}$ $Q_3(u,v) = \frac{(1+u)(1+v)}{4}$ $Q_4(u,v) = \frac{(1u)(1+v)}{4}$ transform the interior of your quadrilateral in the interior of the "standard quadrilateral" $SQ = \{(u,v) / 1\le u \le 1, 1\le v \le 1\}$ Now your original integral $A = \int_c d^2x = \int_{1}^{1}\int_{1}^1 J du dv$ Where $J$ is the Jacobian of the transformation, that can be easily computed to be $J = \alpha_0 + \alpha_1 u + \alpha_2 v$ with $\alpha_0 = \frac{1}{8}\left[ (x_4x_2)(y_1  y_3) + (x_3x_1)(y_4y_2)\right]$ (the values of $\alpha_{1,2}$ are irrelevant for your qustion because do not contribute to the value of teh Area.). Now the value of the area is simply $A = 4\alpha_0$ Hope this is what you want! 


Under the assumption, that the quadrilateral is convex, the edges constituting to the matching with maximal edgeweight sum, are the diagonals. There is also a wiki page dedicated to quadrialterals (http://en.wikipedia.org/wiki/Quadrilateral), which also discusses, how to calculate the are of convex quadrilaterals in various ways; among it the remarkable Bretschneider's formula. 

