# Computing circumcenter of triangle in 2D with MATLAB

I'm writing a finite volume program over a 2D triangular mesh, and at one point I need to calculate the circumcenters of the triangles. The equation given in class and that on Wikipedia give different results.

In class:

$$\vec{r}_c = \vec{r}_i + \frac{( {\lVert \vec{r}_{ik} \rVert}^2 \vec{r}_{ij} - {\lVert \vec{r}_{ij} \rVert}^2 \vec{r}_{ik}) \times \vec{r}_{ik} \times \vec{r}_{ij}} { 2 {\lVert \vec{r}_{ik} \times \vec{r}_{ij} \rVert}^2}$$

$$x_c = \frac{{\lVert A \rVert}^2(B_y-C_y)+{\lVert B \rVert}^2(C_y-A_y)+{\lVert C \rVert}^2(A_y-B_y)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}$$

$$y_c = \frac{{\lVert A \rVert}^2(C_x-B_x)+{\lVert B \rVert}^2(A_x-C_x)+{\lVert C \rVert}^2(B_x-A_x)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}$$

Where $\vec{r}_{ab} = \vec{r}_b - \vec{r}_a$, and $(\vec{r}_i,\vec{r}_j,\vec{r}_k)$ and $(A,B,C)$ are the triangle vertices.

I implented both after running into problems with the final solution (cross section of fluid velocity in a strangely-shaped tube), and found they gave different solutions.

Plotting the triangles and calculated circumcenters gives the following.

Using in-class formula:

Using Wikipedia formula:

Circumcenters were calculated with the following lines of code. In-class formula:

C = [I,0]+(cross(cross(((norm(K-I))^2*[J-I,0]-(norm(J-I))^2*[K-I,0]),[K-I,0]),[J-I,0]))/(2*(norm(cross([K-I,0],[J-I,0])))^2);
C = C(1,1:2);


Wikipedia:

D = 2*(I(1)*(J(2)-K(2))+J(1)*(K(2)-I(2))+K(1)*(I(2)-J(2)));
xc = ((I(2)^2+I(1)^2)*(J(2)-K(2))+(J(2)^2+J(1)^2)*(K(2)-I(2))+(K(2)^2+K(1)^2)*(I(2)-J(2)))/D;
yc = ((I(1)^2+I(2)^2)*(J(1)-K(1))+(J(1)^2+J(2)^2)*(K(1)-I(1))+(K(1)^2+K(2)^2)*(I(1)-J(1)))/D;
C = [xc -yc];


I figure that there must be something wrong with my implementation of the in-class formula, as many acute triangles don't have a calculated circumcenter inside of them, but I can't find the error. Can you?

PS: Yeah, I see the more elegant methods of implementing the Wikipedia formula. Meh.

-
I don't see where any serious roundoff error could be generated here. I can't find an error in the code, either. It's noticeable that some of the circumcentres come out exactly right and others don't. Are you sure you've initialized the $z$ components to zero? The in-class version of the code assumes that (in forming the cross products) whereas the other one doesn't. – joriki Nov 18 '11 at 8:24
Your 'in-class' formula is ambiguous because the cross product isn't associative, that is $a\times (b\times c) \neq (a\times b)\times c$. Which one is intended? – Chris Taylor Nov 18 '11 at 8:26
Yes, @joriki, I set the z-components to zero by passing [P,0] to the cross(...) function (where P==[Px,Py]). – tyblu Nov 18 '11 at 8:53
@Chris Taylor, I'm not really sure -- that is a serious omission. It gives me an idea as to how to answer my own question... – tyblu Nov 18 '11 at 8:56
@tyblu: I see, sorry, I don't know the language you're using; should have guessed that that's what that's for... – joriki Nov 18 '11 at 9:03

Substituting variables $(\vec{r}_i,\vec{r}_j,\vec{r}_k) = (A,B,C)$, with all z-coordinates equal to zero, into the 'in-class' formula gives the following:

$$\vec{r}_c = A+\frac{({\lVert C-A \rVert}^2(B-A)-{\lVert B-A \rVert}^2(C-A)) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2}$$

$$\vec{r}_c = A+\frac{(({\lVert A \rVert}^2+{\lVert C \rVert}^2-2A \cdot C)(B-A)-({\lVert A \rVert}^2+{\lVert B \rVert}^2-2A \cdot B)(C-A)) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2}$$

$$\vec{r}_c = A+\frac{\left({\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A) -2\left[(A\cdot C)(B-A)-(A\cdot B)(C-A)\right]\right) \times (C-A) \times (B-A)}{2{\lVert (C-A) \times (B-A) \rVert}^2}$$

Taking a hint from the denominator, I assume the cross product follows the following order. (The terms to the right combine to form a unit vector.)

$$\vec{r}_c = A+\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A) -2\left[(A\cdot C)(B-A)-(A\cdot B)(C-A)\right]}{2{\lVert (C-A)\times(B-A) \rVert}} \times \frac{(C-A)\times(B-A)}{\lVert (C-A)\times(B-A)\rVert}$$

$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{(A\cdot C)(B-A)-(A\cdot B)(C-A)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$

Let $Y=B-A$, $Z=C-A$, $A\cdot B=X\cdot Y$, and $A\cdot C=X\cdot Z$, giving the following:

$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{(X\cdot Z)Y-(X\cdot Y)Z}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$

$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}-\frac{X\times(Y\times Z)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$

$$\vec{r}_c = A+\left(\frac{{\lVert A \rVert}^2 (B-C) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}+X\times\frac{(C-A)\times(B-A)}{\lVert (C-A)\times(B-A)\rVert}\right)\times \hat z$$

$$\vec{r}_c = A+(X\times\hat z)\times\hat z +\frac{{\lVert A \rVert}^2 (C-B) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}\times\hat z$$

Now, if I knew how to solve for $-X$, I'd show the following:

$$\vec{r}_c = \frac{{\lVert A \rVert}^2 (C-B) + {\lVert B \rVert}^2 (A-C) + {\lVert C \rVert}^2 (B-A)}{2(A_x(B_y-C_y)+B_x(C_y-A_y)+C_x(A_y-B_y))}\times\hat z$$

In other words, the two formulas are off by a factor of -1, and it appears to be the Wikipedia one that maps improperly, as I had been 'fixing' it by setting yc=-yc -- the x-axis is symmetric, so I didn't notice that it was also flipped.

The 'in-class' formula needs to be implemented like this: (@Chris Taylor hit it on the head!)

C = [I,0]+(cross((norm(K-I))^2*[J-I,0]-(norm(J-I))^2*[K-I,0],cross([K-I,0],[J-I,0])))/(2*(norm(cross([K-I,0],[J-I,0])))^2);
C = C(1,1:2);


I'll have to try it out later, as it's time to head to class!

-
Would be great if someone could double-check this, so I or another can correct the Wikipedia page (perhaps with this or similar derivation...). – tyblu Nov 18 '11 at 16:57