v0.14.0
tricircumcenter.c
Go to the documentation of this file.
1 // Jonathan R Shewchuk <jrs+@cs.cmu.edu>
2 /*****************************************************************************/
3 /* */
4 /* tricircumcenter() Find the circumcenter of a triangle. */
5 /* */
6 /* The result is returned both in terms of x-y coordinates and xi-eta */
7 /* coordinates, relative to the triangle's point `a' (that is, `a' is */
8 /* the origin of both coordinate systems). Hence, the x-y coordinates */
9 /* returned are NOT absolute; one must add the coordinates of `a' to */
10 /* find the absolute coordinates of the circumcircle. However, this means */
11 /* that the result is frequently more accurate than would be possible if */
12 /* absolute coordinates were returned, due to limited floating-point */
13 /* precision. In general, the circumradius can be computed much more */
14 /* accurately. */
15 /* */
16 /* The xi-eta coordinate system is defined in terms of the triangle. */
17 /* Point `a' is the origin of the coordinate system. The edge `ab' extends */
18 /* one unit along the xi axis. The edge `ac' extends one unit along the */
19 /* eta axis. These coordinate values are useful for linear interpolation. */
20 /* */
21 /* If `xi' is NULL on input, the xi-eta coordinates will not be computed. */
22 /* */
23 /*****************************************************************************/
24 
25 #include <stdlib.h>
26 
27 void tricircumcenter_tp(a, b, c, circumcenter, xi, eta)
28 double a[2];
29 double b[2];
30 double c[2];
31 double circumcenter[2];
32 double *xi;
33 double *eta;
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 }
76 
77 /*****************************************************************************/
78 /* */
79 /* tricircumcenter3d() Find the circumcenter of a triangle in 3D. */
80 /* */
81 /* The result is returned both in terms of xyz coordinates and xi-eta */
82 /* coordinates, relative to the triangle's point `a' (that is, `a' is */
83 /* the origin of both coordinate systems). Hence, the xyz coordinates */
84 /* returned are NOT absolute; one must add the coordinates of `a' to */
85 /* find the absolute coordinates of the circumcircle. However, this means */
86 /* that the result is frequently more accurate than would be possible if */
87 /* absolute coordinates were returned, due to limited floating-point */
88 /* precision. In general, the circumradius can be computed much more */
89 /* accurately. */
90 /* */
91 /* The xi-eta coordinate system is defined in terms of the triangle. */
92 /* Point `a' is the origin of the coordinate system. The edge `ab' extends */
93 /* one unit along the xi axis. The edge `ac' extends one unit along the */
94 /* eta axis. These coordinate values are useful for linear interpolation. */
95 /* */
96 /* If `xi' is NULL on input, the xi-eta coordinates will not be computed. */
97 /* */
98 /*****************************************************************************/
99 
100 void tricircumcenter3d_tp(a, b, c, circumcenter, xi, eta)
101 double a[3];
102 double b[3];
103 double c[3];
104 double circumcenter[3];
105 double *xi;
106 double *eta;
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 }
tricircumcenter3d_tp
void tricircumcenter3d_tp(a, b, c, circumcenter, double *xi, double *eta)
Definition: tricircumcenter.c:100
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
a
constexpr double a
Definition: approx_sphere.cpp:30
eta
double eta
Definition: free_surface.cpp:170
tricircumcenter_tp
void tricircumcenter_tp(a, b, c, circumcenter, double *xi, double *eta)
Definition: tricircumcenter.c:27
Generated by Doxygen 1.8.17 and hosted at University of Glasgow