\input LetterFormat \input GeoMacros \input KineMacros

\begfile{ForrestReply.tex of December 9, 1983 7:00 am --- Stolfi}

\xeroxheader

\date{November 29, 1983}

\sendto{Prof. Robin Forrest\lb
School of Computer Studies and Accountancy\lb
University of East Anglia\lb
??? ??? (England)}

\oh{Dear Prof. Forrest:}

A few months ago, you left us some comments on our Voronoi paper,\foot\dag{{\it Primitives for the manipulation of general subdivisions and the computation of Voronoi diagrams}, by Leo Guibas and Jorge Stolfi. Proc. 15th. Ann. ACM Symposium on Theory of Computing (STOC), April 1983, 221--234} due (if I am not mistaken) to a student of yours.

We were pleased to know about his independent implementation of our ideas. We found his comments quite valuable; the two errors he found in the paper had indeed escaped our attention. May I ask you to send him a copy of the the attached reply? Please transmit him also our thanks, and my apologies for the delay.

Thank you very much,

\signed{Jorge Stolfi}

\fillit\eject

I will try to answer your coments one by one, but let me give first a general disclaimer. The neat typesetting of our paper, which is largely due to the wonders of Don Knuth's TeX and to the skill of our secretary, may hide in some measure the haste in which it was written. Most of it was put together in three frantic weeks, and some of the ideas (particularly in the subdivision section) were developed in detail only a few days before the submission deadline. That you found errors in it is not so surprising; what is surprising is that you didn't find more of them. We had to omit many important details, for lack of both time and space; indeed, had we got more time, we would probably have shortened the paper even more.

{\bf The sign of the InCircle determinant.} You are quite right on that one: InCircle is true if the determinant is positive, rather than negative. That was truly an error in our algebra, not just in out typing, and we thank you for pointing it out.

You will perhaps be interested in the story behind that error. We implemented the two Voronoi algorithms (recursive and incremental) in the summer of '82, about six months before writing the final version of the paper. The quad-edge data structure and the $\incircle$ predicate were ``invented'' at that time, but were still missing many of the details and the theory given in the paper. The quad edge structure, for example, was restricted to embeddings in orientable manifolds, and we were still unaware of the determinant form of the $\incircle$ predicate and of its symmetry properties.

We implemented $\incircle(a,b,c,d)$ in a very geometrical way, by computing the bisectors of $ac$ and $bd$, finding their intersection point $o$, and testing whether $dist(a,o)>\dist(b,o)$. For this computation we assumed $abcd$ to be a non-degenerate simple quadrilateral (not necessarily convex), and $abc$ a counterclockwise triangle; these assumptions were true in our Voronoi algorithms. Except for input and output, we did all geometric computations in homogeneous coordinates (mostly to handle Voronoi vertices at infinity), but this was not particularly important for the $\incircle$ predicate, in view of the above assumption.

It was only after the implementation was completed, while we were trying to clean up the theory for the paper, that we fully realized the symmetry properties of $\incircle$, namely $\incircle(a,b,c,d)\eqv\incircle(b,c,d,a)$ and $\incircle(a,b,c,d)\eqv\outcircle(b,a,c,d)$. The latter required us to redefine the ``inside'' of a {\it clockwise} oriented circle as being what we would normally call the outside. All these symmetries hinted at a simple algebraic formula for $\incircle$.

We thought at first that formula would come out of Ptolemy's theorem, which says that a simple quadrilateral $abcd$ is inscribed in a circle if and only if it satisfies $$\bar{ab}\times\bar{cd}+\bar{bc}\times\bar{da}=\bar{ac}\times\bar{bd}.\eqno(1)$$
We guessed that the equal sign would be replaced by $>$ or $<$ depending on whether we had $\incircle(a,b,c,d)$ or $\outcircle(a,b,c,d)$ (or vice-versa). So, we wrote down (1) with $>$ instead of $=$, and removed the square roots implicit in $bar{ab}, \bar{bc}$, etc., by the textbook method. The result was the inequality
$$4\times{\bar{ab}}^2\times{\bar{cd}}^2\times{\bar{bc}}^2\times{\bar{da}}^2 - ({\bar{ac}}^2\times{\bar{bd}}^2-{\bar{ab}}^2\times{\bar{cd}}^2-{\bar{bc}}^2\times{\bar{da}}^2)^2 > 0\eqno(2)$$
which is an eight-degree polynomial inequality on the coordinates of the points.

That was not precisely the ``simple algebraic formula'' we were expecting. Instead of thinking a bit more, we went straight for the high-tech, high-power solution: we connected across country to the MACLISP symbolic manipulation system at MIT, and commanded it to expand and then simplify the left-hand side of (2), substituting. MACLISP chewed at the problem for a few minutes, and then asked if it was OK for it to request more memory space. Yes, we replied. This happened another four or five times, until MACLISP ran completely out of memory and had to give up. Only then did we realize that the expanded form of (2) has $6^4+(3\cdot6^2)^2=12960$ terms...

A bit discouraged, but still unwilling to concede defeat, we decided to let $d$ be the point $(0,0)$ and try again (after all, $\incircle$ should be independent of a translation). This substitution reduced the number of terms to a mere $6^2\cdot2
^2 + (3\cdot6\cdot2)^2=1450,$ which was within MACLISP limits. The simplified output was not much shorter than (2), but had a totally different and more regular structure:
$$\biglp (x
a^2 +y
a^2)(x
by
c - x
cy
b) - (x
b^2 +y
b^2)(x
ay
c + x
cy
a) + (x
c^2 +y
c^2)(x
ay
b + x
by
a)\bigrp^2,\eqno(3)$$
which we immediately recognized as the {\it square} of the determinant
$$\detmat{\matrix3{
x
ay
a(x
a^2 +y
a^2)\cr
x
by
b(x
b^2 +y
b^2)\cr
x
cy
c(x
c^2 +y
c^2)\cr}}\eqno(4)$$
This showed us we were wrong in assuming (2) to be $\incircle$ (or its negation): the left-hand side is never negative, and is zero only when the four points are cocircular. It follows that in (1) the left-hand side is never smaller than the right-hand side. Worse than that, we should have noticed somehing was wrong just by looking at (1) a bit more carefully: a cyclic permutation of $a,b,c,d$, which should negate $\incircle$, does not affect either side of (1).

Anyway, it was natural to guess that $\incircle$ would be determined by the sign of the determinant in (4), before squaring. Also, considering the alternating behavior of $\incircle$ under permutations, it was natural to guess that the general case of (4) (with arbitrary $d$) would be the $4\times4$ determinant
$$\detmat{\matrix4{
x
ay
a(x
a^2 +y
a^2)1\cr
x
by
b(x
b^2 +y
b^2)1\cr
x
cy
c(x
c^2 +y
c^2)1\cr
x
dy
d(x
d^2 +y
d^2)1\cr
}}\eqno(5)$$
It was easy to prove this determinant vanishes if and only if the four points are cocircular, and that it changes sign when one of them crosses the circle determined by the other three. Indeed, formula (5) is so obvious that we should have derived it easily, had we approached the problem the right way. We knew that a point $d$ is inside or outside a circle $C$ depending on the sign of a quadratic form
$$P(d) = U(x
d^2 + y
d^2) + V x
d + W y
d - 1.\eqno(6)$$
If we had bothered to compute the coefficients $U,V,W$ for the circle passing though $a,b,c$, we woud have got
{\bf STOPPED HERE}
$$U={\detmat{\matrix3{
x
ay
a1\cr
x
by
b1\cr
x
cy
c1\cr}}\over
\detmat{\matrix3{
x
ay
a(x
a^2 +y
a^2)\cr
x
by
b(x
b^2 +y
b^2)\cr
x
cy
c(x
c^2 +y
c^2)\cr}}\eqno(4)$$

\fillit\end