MeshKit  1.0
circumcenter.cpp
Go to the documentation of this file.
```00001
00002 #include "meshkit/circumcenter.hpp"
00003
00004 #include <stdlib.h>
00005 #include <stdio.h>
00006
00007 #include <iostream>
00008 using namespace std;
00009
00010 //
00011 //Let a, b, c, d, and m be vectors in R^3.  Let ax, ay, and az be the components
00012 //of a, and likewise for b, c, and d.  Let |a| denote the Euclidean norm of a,
00013 //and let a x b denote the cross product of a and b.  Then
00014 //
00015 //    |                                                                       |
00016 //    | |d-a|^2 [(b-a)x(c-a)] + |c-a|^2 [(d-a)x(b-a)] + |b-a|^2 [(c-a)x(d-a)] |
00017 //    |                                                                       |
00018 //r = -------------------------------------------------------------------------,
00019 //                              | bx-ax  by-ay  bz-az |
00020 //                            2 | cx-ax  cy-ay  cz-az |
00021 //                              | dx-ax  dy-ay  dz-az |
00022 //
00023 //and
00024 //
00025 //        |d-a|^2 [(b-a)x(c-a)] + |c-a|^2 [(d-a)x(b-a)] + |b-a|^2 [(c-a)x(d-a)]
00026 //m = a + ---------------------------------------------------------------------.
00027 //                                | bx-ax  by-ay  bz-az |
00028 //                              2 | cx-ax  cy-ay  cz-az |
00029 //                                | dx-ax  dy-ay  dz-az |
00030 //
00031 //Some notes on stability:
00032 //
00033 //- Note that the expression for r is purely a function of differences between
00034 //  coordinates.  The advantage is that the relative error incurred in the
00035 //  computation of r is also a function of the _differences_ between the
00036 //  vertices, and is not influenced by the _absolute_ coordinates of the
00037 //  vertices.  In most applications, vertices are usually nearer to each other
00038 //  than to the origin, so this property is advantageous.  If someone gives you
00039 //  a formula that doesn't have this property, be wary.
00040 //
00041 //  Similarly, the formula for m incurs roundoff error proportional to the
00042 //  differences between vertices, but does not incur roundoff error proportional
00043 //  to the absolute coordinates of the vertices until the final addition.
00044
00045 //- These expressions are unstable in only one case:  if the denominator is
00046 //  close to zero.  This instability, which arises if the tetrahedron is nearly
00047 //  degenerate, is unavoidable.  Depending on your application, you may want
00048 //  to use exact arithmetic to compute the value of the determinant.
00049 //  Fortunately, this determinant is the basis of the well-studied 3D orientation
00050 //  test, and several fast algorithms for providing accurate approximations to
00051 //  the determinant are available.  Some resources are available from the
00052 //  "Numerical and algebraic computation" page of Nina Amenta's Directory of
00053 //  Computational Geometry Software:
00054
00055 //  http://www.geom.umn.edu/software/cglist/alg.html
00056
00057 //  If you're using floating-point inputs (as opposed to integers), one
00058 //  package that can estimate this determinant somewhat accurately is my own:
00059
00060 //  http://www.cs.cmu.edu/~quake/robust.html
00061
00062 //- If you want to be even more aggressive about stability, you might reorder
00063 //  the vertices of the tetrahedron so that the vertex `a' (which is subtracted
00064 //  from the other vertices) is, roughly speaking, the vertex at the center of
00065 //  the minimum spanning tree of the tetrahedron's four vertices.  For more
00067 //  Algorithms for 2D Delaunay Triangulations," International Journal of
00068 //  Computational Geometry & Applications 5(1-2):193-213, March-June 1995.
00069
00070 //For completeness, here are stable expressions for the circumradius and
00071 //circumcenter of a triangle, in R^2 and in R^3.  Incidentally, the expressions
00072 //given here for R^2 are better behaved in terms of relative error than the
00073 //suggestions currently given in the Geometry Junkyard
00074 //(http://www.ics.uci.edu/~eppstein/junkyard/triangulation.html).
00075
00076 //Triangle in R^2:
00077 //
00078 //     |b-a| |c-a| |b-c|            < Note: You only want to compute one sqrt, so
00079 //r = ------------------,             use sqrt{ |b-a|^2 |c-a|^2 |b-c|^2 }
00080 //      | bx-ax  by-ay |
00081 //    2 | cx-ax  cy-ay |
00082 //
00083 //          | by-ay  |b-a|^2 |
00084 //          | cy-ay  |c-a|^2 |
00085 //mx = ax - ------------------,
00086 //            | bx-ax  by-ay |
00087 //          2 | cx-ax  cy-ay |
00088 //
00089 //          | bx-ax  |b-a|^2 |
00090 //          | cx-ax  |c-a|^2 |
00091 //my = ay + ------------------.
00092 //            | bx-ax  by-ay |
00093 //          2 | cx-ax  cy-ay |
00094 //
00095 //Triangle in R^3:
00096 //
00097 //    |                                                           |
00098 //    | |c-a|^2 [(b-a)x(c-a)]x(b-a) + |b-a|^2 (c-a)x[(b-a)x(c-a)] |
00099 //    |                                                           |
00100 //r = -------------------------------------------------------------,
00101 //                         2 | (b-a)x(c-a) |^2
00102 //
00103 //        |c-a|^2 [(b-a)x(c-a)]x(b-a) + |b-a|^2 (c-a)x[(b-a)x(c-a)]
00104 //m = a + ---------------------------------------------------------.
00105 //                           2 | (b-a)x(c-a) |^2
00106 //
00107 //Finally, here's some C code that computes the vector m-a (whose norm is r) in
00108 //each of these three cases.  Notice the #ifdef statements, which allow you to
00109 //choose whether or not to use my aforementioned package to approximate the
00110 //determinant.  (No attempt is made here to reorder the vertices to improve
00111 //stability.)
00112
00113 /*****************************************************************************/
00114 /*                                                                           */
00115 /*  tetcircumcenter()   Find the circumcenter of a tetrahedron.              */
00116 /*                                                                           */
00117 /*  The result is returned both in terms of xyz coordinates and xi-eta-zeta  */
00118 /*  coordinates, relative to the tetrahedron's point `a' (that is, `a' is    */
00119 /*  the origin of both coordinate systems).  Hence, the xyz coordinates      */
00120 /*  returned are NOT absolute; one must add the coordinates of `a' to        */
00121 /*  find the absolute coordinates of the circumcircle.  However, this means  */
00122 /*  that the result is frequently more accurate than would be possible if    */
00123 /*  absolute coordinates were returned, due to limited floating-point        */
00124 /*  precision.  In general, the circumradius can be computed much more       */
00125 /*  accurately.                                                              */
00126 /*                                                                           */
00127 /*  The xi-eta-zeta coordinate system is defined in terms of the             */
00128 /*  tetrahedron.  Point `a' is the origin of the coordinate system.          */
00129 /*  The edge `ab' extends one unit along the xi axis.  The edge `ac'         */
00130 /*  extends one unit along the eta axis.  The edge `ad' extends one unit     */
00131 /*  along the zeta axis.  These coordinate values are useful for linear      */
00132 /*  interpolation.                                                           */
00133 /*                                                                           */
00134 /*  If `xi' is NULL on input, the xi-eta-zeta coordinates will not be        */
00135 /*  computed.                                                                */
00136 /*                                                                           */
00137 /*****************************************************************************/
00138
00139 /*****************************************************************************/
00140
00141 void tetcircumcenter(double a[3], double b[3], double c[3], double d[3],
00142                      double circumcenter[3], double *xi, double *eta, double *zeta)
00143 {
00144   double xba, yba, zba, xca, yca, zca, xda, yda, zda;
00145   double balength, calength, dalength;
00146   double xcrosscd, ycrosscd, zcrosscd;
00147   double xcrossdb, ycrossdb, zcrossdb;
00148   double xcrossbc, ycrossbc, zcrossbc;
00149   double denominator;
00150   double xcirca, ycirca, zcirca;
00151
00152   /* Use coordinates relative to point `a' of the tetrahedron. */
00153   xba = b[0] - a[0];
00154   yba = b[1] - a[1];
00155   zba = b[2] - a[2];
00156   xca = c[0] - a[0];
00157   yca = c[1] - a[1];
00158   zca = c[2] - a[2];
00159   xda = d[0] - a[0];
00160   yda = d[1] - a[1];
00161   zda = d[2] - a[2];
00162   /* Squares of lengths of the edges incident to `a'. */
00163   balength = xba * xba + yba * yba + zba * zba;
00164   calength = xca * xca + yca * yca + zca * zca;
00165   dalength = xda * xda + yda * yda + zda * zda;
00166   /* Cross products of these edges. */
00167   xcrosscd = yca * zda - yda * zca;
00168   ycrosscd = zca * xda - zda * xca;
00169   zcrosscd = xca * yda - xda * yca;
00170   xcrossdb = yda * zba - yba * zda;
00171   ycrossdb = zda * xba - zba * xda;
00172   zcrossdb = xda * yba - xba * yda;
00173   xcrossbc = yba * zca - yca * zba;
00174   ycrossbc = zba * xca - zca * xba;
00175   zcrossbc = xba * yca - xca * yba;
00176
00177   /* Calculate the denominator of the formulae. */
00178 #ifdef EXACT
00179   /* Use orient3d() from http://www.cs.cmu.edu/~quake/robust.html     */
00180   /*   to ensure a correctly signed (and reasonably accurate) result, */
00181   /*   avoiding any possibility of division by zero.                  */
00182   denominator = 0.5 / orient3d(b, c, d, a);
00183 #else
00184   /* Take your chances with floating-point roundoff. */
00185   printf( " Warning: IEEE floating points used: Define -DEXACT in makefile \n");
00186   denominator = 0.5 / (xba * xcrosscd + yba * ycrosscd + zba * zcrosscd);
00187 #endif
00188
00189   /* Calculate offset (from `a') of circumcenter. */
00190   xcirca = (balength * xcrosscd + calength * xcrossdb + dalength * xcrossbc) *
00191            denominator;
00192   ycirca = (balength * ycrosscd + calength * ycrossdb + dalength * ycrossbc) *
00193            denominator;
00194   zcirca = (balength * zcrosscd + calength * zcrossdb + dalength * zcrossbc) *
00195            denominator;
00196   circumcenter[0] = xcirca;
00197   circumcenter[1] = ycirca;
00198   circumcenter[2] = zcirca;
00199
00200   if (xi != (double *) NULL) {
00201     /* To interpolate a linear function at the circumcenter, define a    */
00202     /*   coordinate system with a xi-axis directed from `a' to `b',      */
00203     /*   an eta-axis directed from `a' to `c', and a zeta-axis directed  */
00204     /*   from `a' to `d'.  The values for xi, eta, and zeta are computed */
00205     /*   by Cramer's Rule for solving systems of linear equations.       */
00206     *xi = (xcirca * xcrosscd + ycirca * ycrosscd + zcirca * zcrosscd) *
00207           (2.0 * denominator);
00208     *eta = (xcirca * xcrossdb + ycirca * ycrossdb + zcirca * zcrossdb) *
00209            (2.0 * denominator);
00210     *zeta = (xcirca * xcrossbc + ycirca * ycrossbc + zcirca * zcrossbc) *
00211             (2.0 * denominator);
00212   }
00213 }
00214
00215 /*****************************************************************************/
00216 /*****************************************************************************/
00217 /*                                                                           */
00218 /*  tricircumcenter()   Find the circumcenter of a triangle.                 */
00219 /*                                                                           */
00220 /*  The result is returned both in terms of x-y coordinates and xi-eta       */
00221 /*  coordinates, relative to the triangle's point `a' (that is, `a' is       */
00222 /*  the origin of both coordinate systems).  Hence, the x-y coordinates      */
00223 /*  returned are NOT absolute; one must add the coordinates of `a' to        */
00224 /*  find the absolute coordinates of the circumcircle.  However, this means  */
00225 /*  that the result is frequently more accurate than would be possible if    */
00226 /*  absolute coordinates were returned, due to limited floating-point        */
00227 /*  precision.  In general, the circumradius can be computed much more       */
00228 /*  accurately.                                                              */
00229 /*                                                                           */
00230 /*  The xi-eta coordinate system is defined in terms of the triangle.        */
00231 /*  Point `a' is the origin of the coordinate system.  The edge `ab' extends */
00232 /*  one unit along the xi axis.  The edge `ac' extends one unit along the    */
00233 /*  eta axis.  These coordinate values are useful for linear interpolation.  */
00234 /*                                                                           */
00235 /*  If `xi' is NULL on input, the xi-eta coordinates will not be computed.   */
00236 /*                                                                           */
00237 /*****************************************************************************/
00238
00239
00240 /*****************************************************************************/
00241 void tricircumcenter(double a[2], double b[2], double c[2], double circumcenter[2],
00242                      double *xi, double *eta)
00243 {
00244   double xba, yba, xca, yca;
00245   double balength, calength;
00246   double denominator;
00247   double xcirca, ycirca;
00248
00249   /* Use coordinates relative to point `a' of the triangle. */
00250   xba = b[0] - a[0];
00251   yba = b[1] - a[1];
00252   xca = c[0] - a[0];
00253   yca = c[1] - a[1];
00254   /* Squares of lengths of the edges incident to `a'. */
00255   balength = xba * xba + yba * yba;
00256   calength = xca * xca + yca * yca;
00257
00258   /* Calculate the denominator of the formulae. */
00259 #ifdef EXACT
00260   /* Use orient2d() from http://www.cs.cmu.edu/~quake/robust.html     */
00261   /*   to ensure a correctly signed (and reasonably accurate) result, */
00262   /*   avoiding any possibility of division by zero.                  */
00263   denominator = 0.5 / orient2d(b, c, a);
00264 #else
00265   /* Take your chances with floating-point roundoff. */
00266   denominator = 0.5 / (xba * yca - yba * xca);
00267 #endif
00268
00269   /* Calculate offset (from `a') of circumcenter. */
00270   xcirca = (yca * balength - yba * calength) * denominator;
00271   ycirca = (xba * calength - xca * balength) * denominator;
00272   circumcenter[0] = xcirca;
00273   circumcenter[1] = ycirca;
00274
00275   if (xi != (double *) NULL) {
00276     /* To interpolate a linear function at the circumcenter, define a     */
00277     /*   coordinate system with a xi-axis directed from `a' to `b' and    */
00278     /*   an eta-axis directed from `a' to `c'.  The values for xi and eta */
00279     /*   are computed by Cramer's Rule for solving systems of linear      */
00280     /*   equations.                                                       */
00281     *xi = (xcirca * yca - ycirca * xca) * (2.0 * denominator);
00282     *eta = (ycirca * xba - xcirca * yba) * (2.0 * denominator);
00283   }
00284 }
00285
00286 /****************************************************************************/
00287
00288 /*****************************************************************************/
00289 /*                                                                           */
00290 /*  tricircumcenter3d()   Find the circumcenter of a triangle in 3D.         */
00291 /*                                                                           */
00292 /*  The result is returned both in terms of xyz coordinates and xi-eta       */
00293 /*  coordinates, relative to the triangle's point `a' (that is, `a' is       */
00294 /*  the origin of both coordinate systems).  Hence, the xyz coordinates      */
00295 /*  returned are NOT absolute; one must add the coordinates of `a' to        */
00296 /*  find the absolute coordinates of the circumcircle.  However, this means  */
00297 /*  that the result is frequently more accurate than would be possible if    */
00298 /*  absolute coordinates were returned, due to limited floating-point        */
00299 /*  precision.  In general, the circumradius can be computed much more       */
00300 /*  accurately.                                                              */
00301 /*                                                                           */
00302 /*  The xi-eta coordinate system is defined in terms of the triangle.        */
00303 /*  Point `a' is the origin of the coordinate system.  The edge `ab' extends */
00304 /*  one unit along the xi axis.  The edge `ac' extends one unit along the    */
00305 /*  eta axis.  These coordinate values are useful for linear interpolation.  */
00306 /*                                                                           */
00307 /*  If `xi' is NULL on input, the xi-eta coordinates will not be computed.   */
00308 /*                                                                           */
00309 /*****************************************************************************/
00310 /*****************************************************************************/
00311 void tricircumcenter3d(double a[3], double b[3], double c[3], double circumcenter[3],
00312                        double *xi, double *eta)
00313 {
00314   double xba, yba, zba, xca, yca, zca;
00315   double balength, calength;
00316   double xcrossbc, ycrossbc, zcrossbc;
00317   double denominator;
00318   double xcirca, ycirca, zcirca;
00319
00320   /* Use coordinates relative to point `a' of the triangle. */
00321   xba = b[0] - a[0];
00322   yba = b[1] - a[1];
00323   zba = b[2] - a[2];
00324   xca = c[0] - a[0];
00325   yca = c[1] - a[1];
00326   zca = c[2] - a[2];
00327   /* Squares of lengths of the edges incident to `a'. */
00328   balength = xba * xba + yba * yba + zba * zba;
00329   calength = xca * xca + yca * yca + zca * zca;
00330
00331   /* Cross product of these edges. */
00332 #ifdef EXACT
00333   /* Use orient2d() from http://www.cs.cmu.edu/~quake/robust.html     */
00334   /*   to ensure a correctly signed (and reasonably accurate) result, */
00335   /*   avoiding any possibility of division by zero.                  */
00336
00337   A[0] = b[1]; A[1] = b[2];
00338   B[0] = c[1]; B[1] = c[2];
00339   C[0] = a[1]; C[1] = a[2];
00340   xcrossbc = orient2d(A, B, C);
00341
00342   A[0] = c[0]; A[1] = c[2];
00343   B[0] = b[0]; B[1] = b[2];
00344   C[0] = a[0]; C[1] = a[2];
00345   ycrossbc = orient2d(A, B, C);
00346
00347   A[0] = b[0]; A[1] = b[1];
00348   B[0] = c[0]; B[1] = c[1];
00349   C[0] = a[0]; C[1] = a[1];
00350   zcrossbc = orient2d(A, B, C);
00351
00352   /*
00353   xcrossbc = orient2d(b[1], b[2], c[1], c[2], a[1], a[2]);
00354   ycrossbc = orient2d(b[2], b[0], c[2], c[0], a[2], a[0]);
00355   zcrossbc = orient2d(b[0], b[1], c[0], c[1], a[0], a[1]);
00356   */
00357 #else
00358   printf( " Warning: IEEE floating points used: Define -DEXACT in makefile \n");
00359   /* Take your chances with floating-point roundoff. */
00360   xcrossbc = yba * zca - yca * zba;
00361   ycrossbc = zba * xca - zca * xba;
00362   zcrossbc = xba * yca - xca * yba;
00363 #endif
00364
00365   /* Calculate the denominator of the formulae. */
00366   denominator = 0.5 / (xcrossbc * xcrossbc + ycrossbc * ycrossbc +
00367                        zcrossbc * zcrossbc);
00368
00369   /* Calculate offset (from `a') of circumcenter. */
00370   xcirca = ((balength * yca - calength * yba) * zcrossbc -
00371             (balength * zca - calength * zba) * ycrossbc) * denominator;
00372   ycirca = ((balength * zca - calength * zba) * xcrossbc -
00373             (balength * xca - calength * xba) * zcrossbc) * denominator;
00374   zcirca = ((balength * xca - calength * xba) * ycrossbc -
00375             (balength * yca - calength * yba) * xcrossbc) * denominator;
00376   circumcenter[0] = xcirca;
00377   circumcenter[1] = ycirca;
00378   circumcenter[2] = zcirca;
00379
00380   if (xi != (double *) NULL) {
00381     /* To interpolate a linear function at the circumcenter, define a     */
00382     /*   coordinate system with a xi-axis directed from `a' to `b' and    */
00383     /*   an eta-axis directed from `a' to `c'.  The values for xi and eta */
00384     /*   are computed by Cramer's Rule for solving systems of linear      */
00385     /*   equations.                                                       */
00386
00387     /* There are three ways to do this calculation - using xcrossbc, */
00388     /*   ycrossbc, or zcrossbc.  Choose whichever has the largest    */
00389     /*   magnitude, to improve stability and avoid division by zero. */
00390     if (((xcrossbc >= ycrossbc) ^ (-xcrossbc > ycrossbc)) &&
00391         ((xcrossbc >= zcrossbc) ^ (-xcrossbc > zcrossbc))) {
00392       *xi = (ycirca * zca - zcirca * yca) / xcrossbc;
00393       *eta = (zcirca * yba - ycirca * zba) / xcrossbc;
00394     } else if ((ycrossbc >= zcrossbc) ^ (-ycrossbc > zcrossbc)) {
00395       *xi = (zcirca * xca - xcirca * zca) / ycrossbc;
00396       *eta = (xcirca * zba - zcirca * xba) / ycrossbc;
00397     } else {
00398       *xi = (xcirca * yca - ycirca * xca) / zcrossbc;
00399       *eta = (ycirca * xba - xcirca * yba) / zcrossbc;
00400     }
00401   }
00402 }
00403 /****************************************************************************/
00404 void TriCircumCenter2D( double *a, double *b, double *c, double *result,
00405                         double *param)
00406 {
00407    tricircumcenter(a, b, c, result, &param[0], &param[1]);
00408
00409    result[0] += a[0];
00410    result[1] += a[1];
00411 }
00412 /****************************************************************************/
00413 void TriCircumCenter3D( double *a, double *b, double *c, double *result,
00414                          double *param)
00415 {
00416    tricircumcenter3d(a, b, c, result, &param[0], &param[1]);
00417    result[0] += a[0];
00418    result[1] += a[1];
00419    result[2] += a[2];
00420 }
00421
00422 /****************************************************************************/
00423 void TriCircumCenter3D( double *a, double *b, double *c, double *result)
00424 {
00425    double xi, eta;
00426    tricircumcenter3d(a, b, c, result, &xi, & eta);
00427    result[0] += a[0];
00428    result[1] += a[1];
00429    result[2] += a[2];
00430 }
00431 /****************************************************************************/
00432 void TetCircumCenter( double *a, double *b, double *c, double *d, double *result,
00433                       double *param)
00434 {
00435    double orient = orient3d(a, b, c, d);
00436
00437    if(orient < 0.0)
00438       tetcircumcenter(a, c, b, d, result, &param[0], &param[1], &param[2]);
00439    else
00440       tetcircumcenter(a, b, c, d, result, &param[0], &param[1], &param[2]);
00441
00442    result[0] += a[0];
00443    result[1] += a[1];
00444    result[2] += a[2];
00445 }
00446 /****************************************************************************/
00447 int UnitTest:: test_tet_circumcenter()
00448 {
00449       Point3D  a, b, c, d;
00450       Point3D  result, param;
00451       Array4D  dist;
00452
00453       exactinit();
00454       for( int i = 0; i < 100; i++) {
00455            a[0] = -5.0 + 10*drand48();
00456            a[1] = -5.0 + 10*drand48();
00457            a[2] = -5.0 + 10*drand48();
00458
00459            b[0] = -5.0 + 10*drand48();
00460            b[1] = -5.0 + 10*drand48();
00461            b[2] = -5.0 + 10*drand48();
00462
00463            c[0] = -5.0 + 10*drand48();
00464            c[1] = -5.0 + 10*drand48();
00465            c[2] = -5.0 + 10*drand48();
00466
00467            d[0] = -5.0 + 10*drand48();
00468            d[1] = -5.0 + 10*drand48();
00469            d[2] = -5.0 + 10*drand48();
00470
00471            TetCircumCenter( &a[0], &b[0], &c[0], &d[0], &result[0], &param[0]);
00472
00473            dist[0]  = Math::length(a, result);
00474            dist[1]  = Math::length(b, result);
00475            dist[2]  = Math::length(c, result);
00476            dist[3]  = Math::length(d, result);
00477
00478            if( fabs(dist[1]-dist[0]) > 1.0E-05) {
00479                cout << "Info: Tet CircumCenter failed " << endl;
00480                return 1;
00481            }
00482            if( fabs(dist[2]-dist[0]) > 1.0E-05) {
00483                cout << "Info: Tet CircumCenter failed " << endl;
00484                return 1;
00485            }
00486            if( fabs(dist[3]-dist[0]) > 1.0E-05) {
00487                cout << "Info: Tet CircumCenter failed " << endl;
00488                return 1;
00489            }
00490       }
00491
00492       cout << "Info: Tetrahedra Circumcenter tests passed " << endl;
00493
00494       return 0;
00495 }
00496
00497 int UnitTest:: test_tri3d_circumcenter()
00498 {
00499       Point3D  a, b, c;
00500       Point3D  result, param;
00501       Array4D  dist;
00502
00503       exactinit();
00504       for( int i = 0; i < 100; i++) {
00505            a[0] = -5.0 + 10*drand48();
00506            a[1] = -5.0 + 10*drand48();
00507            a[2] = -5.0 + 10*drand48();
00508
00509            b[0] = -5.0 + 10*drand48();
00510            b[1] = -5.0 + 10*drand48();
00511            b[2] = -5.0 + 10*drand48();
00512
00513            c[0] = -5.0 + 10*drand48();
00514            c[1] = -5.0 + 10*drand48();
00515            c[2] = -5.0 + 10*drand48();
00516
00517            TriCircumCenter3D( &a[0], &b[0], &c[0], &result[0], &param[0]);
00518
00519            dist[0]  = Math::length(a, result);
00520            dist[1]  = Math::length(b, result);
00521            dist[2]  = Math::length(c, result);
00522
00523            if( fabs(dist[1]-dist[0]) > 1.0E-05) {
00524                cout << "Info: Tet CircumCenter failed " << endl;
00525                return 1;
00526            }
00527
00528            if( fabs(dist[2]-dist[0]) > 1.0E-05) {
00529                cout << "Info: Tet CircumCenter failed " << endl;
00530                return 1;
00531            }
00532
00533       }
00534
00535       cout << "Info: 3D Triangle Circumcenter tests passed " << endl;
00536
00537       return 0;
00538 }
```