v0.9.1
tricircumcenter.c File Reference
#include <stdlib.h>

Go to the source code of this file.

## Functions

void tricircumcenter_tp (a, b, c, circumcenter, double *xi, double *eta)

void tricircumcenter3d_tp (a, b, c, circumcenter, double *xi, double *eta)

## ◆ tricircumcenter3d_tp()

 void tricircumcenter3d_tp ( a , b , c , circumcenter , double * xi, double * eta )

Definition at line 100 of file tricircumcenter.c.

107 {
108  double xba, yba, zba, xca, yca, zca;
109  double balength, calength;
110  double xcrossbc, ycrossbc, zcrossbc;
111  double denominator;
112  double xcirca, ycirca, zcirca;
113
114  /* Use coordinates relative to point a' of the triangle. */
115  xba = b[0] - a[0];
116  yba = b[1] - a[1];
117  zba = b[2] - a[2];
118  xca = c[0] - a[0];
119  yca = c[1] - a[1];
120  zca = c[2] - a[2];
121  /* Squares of lengths of the edges incident to a'. */
122  balength = xba * xba + yba * yba + zba * zba;
123  calength = xca * xca + yca * yca + zca * zca;
124
125  /* Cross product of these edges. */
126 #ifdef EXACT
127  /* Use orient2d() from http://www.cs.cmu.edu/~quake/robust.html */
128  /* to ensure a correctly signed (and reasonably accurate) result, */
129  /* avoiding any possibility of division by zero. */
130  xcrossbc = orient2d(b[1], b[2], c[1], c[2], a[1], a[2]);
131  ycrossbc = orient2d(b[2], b[0], c[2], c[0], a[2], a[0]);
132  zcrossbc = orient2d(b[0], b[1], c[0], c[1], a[0], a[1]);
133 #else
134  /* Take your chances with floating-point roundoff. */
135  xcrossbc = yba * zca - yca * zba;
136  ycrossbc = zba * xca - zca * xba;
137  zcrossbc = xba * yca - xca * yba;
138 #endif
139
140  /* Calculate the denominator of the formulae. */
141  denominator = 0.5 / (xcrossbc * xcrossbc + ycrossbc * ycrossbc +
142  zcrossbc * zcrossbc);
143
144  /* Calculate offset (from a') of circumcenter. */
145  xcirca = ((balength * yca - calength * yba) * zcrossbc -
146  (balength * zca - calength * zba) * ycrossbc) * denominator;
147  ycirca = ((balength * zca - calength * zba) * xcrossbc -
148  (balength * xca - calength * xba) * zcrossbc) * denominator;
149  zcirca = ((balength * xca - calength * xba) * ycrossbc -
150  (balength * yca - calength * yba) * xcrossbc) * denominator;
151  circumcenter[0] = xcirca;
152  circumcenter[1] = ycirca;
153  circumcenter[2] = zcirca;
154
155  if (xi != (double *) NULL) {
156  /* To interpolate a linear function at the circumcenter, define a */
157  /* coordinate system with a xi-axis directed from a' to b' and */
158  /* an eta-axis directed from a' to c'. The values for xi and eta */
159  /* are computed by Cramer's Rule for solving systems of linear */
160  /* equations. */
161
162  /* There are three ways to do this calculation - using xcrossbc, */
163  /* ycrossbc, or zcrossbc. Choose whichever has the largest */
164  /* magnitude, to improve stability and avoid division by zero. */
165  if (((xcrossbc >= ycrossbc) ^ (-xcrossbc > ycrossbc)) &&
166  ((xcrossbc >= zcrossbc) ^ (-xcrossbc > zcrossbc))) {
167  *xi = (ycirca * zca - zcirca * yca) / xcrossbc;
168  *eta = (zcirca * yba - ycirca * zba) / xcrossbc;
169  } else if ((ycrossbc >= zcrossbc) ^ (-ycrossbc > zcrossbc)) {
170  *xi = (zcirca * xca - xcirca * zca) / ycrossbc;
171  *eta = (xcirca * zba - zcirca * xba) / ycrossbc;
172  } else {
173  *xi = (xcirca * yca - ycirca * xca) / zcrossbc;
174  *eta = (ycirca * xba - xcirca * yba) / zcrossbc;
175  }
176  }
177 }

## ◆ tricircumcenter_tp()

 void tricircumcenter_tp ( a , b , c , circumcenter , double * xi, double * eta )

Definition at line 27 of file tricircumcenter.c.

34 {
35  double xba, yba, xca, yca;
36  double balength, calength;
37  double denominator;
38  double xcirca, ycirca;
39
40  /* Use coordinates relative to point a' of the triangle. */
41  xba = b[0] - a[0];
42  yba = b[1] - a[1];
43  xca = c[0] - a[0];
44  yca = c[1] - a[1];
45  /* Squares of lengths of the edges incident to a'. */
46  balength = xba * xba + yba * yba;
47  calength = xca * xca + yca * yca;
48
49  /* Calculate the denominator of the formulae. */
50 #ifdef EXACT
51  /* Use orient2d() from http://www.cs.cmu.edu/~quake/robust.html */
52  /* to ensure a correctly signed (and reasonably accurate) result, */
53  /* avoiding any possibility of division by zero. */
54  denominator = 0.5 / orient2d(b, c, a);
55 #else
56  /* Take your chances with floating-point roundoff. */
57  denominator = 0.5 / (xba * yca - yba * xca);
58 #endif
59
60  /* Calculate offset (from a') of circumcenter. */
61  xcirca = (yca * balength - yba * calength) * denominator;
62  ycirca = (xba * calength - xca * balength) * denominator;
63  circumcenter[0] = xcirca;
64  circumcenter[1] = ycirca;
65
66  if (xi != (double *) NULL) {
67  /* To interpolate a linear function at the circumcenter, define a */
68  /* coordinate system with a xi-axis directed from a' to b' and */
69  /* an eta-axis directed from a' to c'. The values for xi and eta */
70  /* are computed by Cramer's Rule for solving systems of linear */
71  /* equations. */
72  *xi = (xcirca * yca - ycirca * xca) * (2.0 * denominator);
73  *eta = (ycirca * xba - xcirca * yba) * (2.0 * denominator);
74  }
75 }