I add our explanation and the origin of the problem.
To obtain the formula we just need to verify the following lemma (by a straightforward computation).
Magic Lemma. Let $(a,b),(c,d)$ be as in our summation in the main posting. Draw tangents to the unit circle which are orthogonal to the directions $(a,b),(c,d),(a+c,b+d)$. Then the triangle in the intersection has the area $\frac{f(a,b,c,d)^2}{2}$.
Now consider the circle given by $x^2+y^2=1$. At every rational point on it we draw the tangent line. Clearly, our disk is the intersection of all the half-planes given by these lines.
So we can approximate the area of the disk by the following sequential cutting. We start with the square $[-1,1]\times [-1,1]$. It has four points of tangency with our circle. We will construct a sequence of polygons converging to the circle. Initially we have sides of directions orthogonal to $(1,0),(0,1),(-1,0),(0,1)$. At every step, we take two vectors $v_1,v_2$(in the same quadrant) which give the basis of $\mathbb Z^2$ and add to our polygon a new side, which is given by the tangent line of the direction orthogonal to $v_1+v_2$.
For example, the first step will be cutting the square by the line $x+y=\sqrt 2$, if we take the vectors $(1,0),(0,1)$. The area of the eliminated triangle is $\frac{(\sqrt2-1)^2}{2}=\frac{f(1,0,0,1)^2}{2}$.
It seems that this procedure can be generalized to higher dimensions — at least a preliminary computation shows that something like the above magic lemma should happen.