Find a tetrahedron circumcenter

C++ code to compute tetrahdron circumcenter - 05/2021 - #Geometry

Tetrahdron's circumcenter

Find the center (a.k.a circumcenter) of the circumscribed sphere of a Tetrahedron.

Definition

For all tetrahedra, there exists a sphere called the circumsphere which completely encloses the tetrahedron. The tetrahedron's vertices all lie on the surface of its circumsphere. The point at the centre of the circumsphere is called the circumcentre.

WebGL demo

"use strict";

import * as THREE from "https://cdn.jsdelivr.net/npm/three@0.114/build/three.module.js";
import { OrbitControls } from "https://cdn.jsdelivr.net/npm/three@0.114/examples/jsm/controls/OrbitControls.js";
import { DragControls } from "https://cdn.jsdelivr.net/npm/three@0.114/examples/jsm/controls/DragControls.js";




/*
    Up vector: 'y' axis
    grid's square len: 100units

*/

var g_camera, g_scene, g_renderer;
var g_controls;

var g_canvasName = document.querySelector("script[type=module]").attributes['canvas_name'].value;
var g_canvaContainer = document.getElementById(g_canvasName);

var g_canvasRootDivId = "threeJsWrap_" + g_canvasName;
var g_canvasRoot = document.getElementById(g_canvasRootDivId );

var g_datGuiInnerDivId = "datGuiInner_" + g_canvasName;
var g_datGuiBottomDivId = "datGuiBottom_" + g_canvasName;
var g_datGuiBottomDiv = document.getElementById(g_datGuiBottomDivId);


/* CodePen: 
function windowWidth() { window.innerWidth;  }
function windowHeight(){ window.innerHeight; }
*/

// Standalone:
function windowWidth() { 
    return g_canvasRoot.clientWidth;  
}

function windowHeight(){ 
    return 450;
    //return g_datGuiBottomDiv.clientHeight; 
}

var g_tetrahedron = null;

// The mesh representing our "center of gravity"
// (i.e. barycenter, circumcenter or incenter according to the selected option)
var g_cogPointMesh;

// The transparent sphere that either:
// - fit the 4 vertices of our tetrahedron. (circumcenter)
// - fit inside the faces of our tetrahedron (incenter)
// - somewhat fits inside the tetrahedron (barycenter/ radius == nearest vertex)
var g_fittedSphereMesh;


var DatGuiContext = function () {
    this.lightColor1 = 0xe3e3e3;
    this.lightColor2 = 0x6b6b6b;
    this.algoType = "Circumcenter";
};


var g_datGuiContext = new DatGuiContext();

// callback when the html page finish loading:
window.onload = function () {
    setValue();
    

    var gui = new dat.GUI({ autoPlace: false, width:"100%"  });  // width: windowWidth()
    var div = document.getElementById(g_datGuiBottomDivId);
    div.appendChild(gui.domElement);

    //gui.addColor(g_datGuiContext, 'lightColor1').onChange(setValue);
    //gui.addColor(g_datGuiContext, 'lightColor2').onChange(setValue);    
    gui.add(g_datGuiContext, 'algoType', ["Circumcenter", "Barycenter", "Incenter"]).onChange(setValue);

    
    

};


function setValue() {
    //d_light.color.setHex( g_datGuiContext.lightColor1 );  
    //h_light.color.setHex( g_datGuiContext.lightColor2 );
}



function mat(col) {
    return new THREE.MeshLambertMaterial({ color: col });
}



function cross( a, b ) {
    let res = new THREE.Vector3();
    const ax = a.x, ay = a.y, az = a.z;
    const bx = b.x, by = b.y, bz = b.z;
    res.x = ay * bz - az * by;
    res.y = az * bx - ax * bz;
    res.z = ax * by - ay * bx;
    return res;
}




class Tetrahedron {

    constructor() {
        let sphereGeometry = new THREE.SphereBufferGeometry(1.0, 50, 50);
        let material = new THREE.MeshPhongMaterial({ color: 0x3399dd });
        this.vertexMesh = [];
        // Initial values
        var vert_positions = {
            x1: new THREE.Vector3(-50, 20, 50),
            x2: new THREE.Vector3(50, 20, -50),
            x3: new THREE.Vector3(-50, 20, -50),
            x4: new THREE.Vector3(0, 70, 0)
        };

        // Add dragable spheres
        for (var v in vert_positions) {
            var sphereMesh = new THREE.Mesh(sphereGeometry, material);
            sphereMesh.scale.set(7, 7, 7);
            var vert = vert_positions[v];
            sphereMesh.position.set(vert.x, vert.y, vert.z);
            sphereMesh.castShadow = true;
            this.vertexMesh.push(sphereMesh);
        }

        this.buildFaces()
    }

    updateGeom() {
        if (false) {
            // The deprecated way
            // and I don't know how to handle flat shading and update normals with this:
            for (let v = 0; v < 4; v++) {
                this.tetMesh.geometry.vertices[v].copy(this.get_vert(v));
            }
            this.tetMesh.geometry.verticesNeedUpdate = true;
        } else {
            for (let v = 0; v < 4; v++) {
                let p = this.get_vert(v);
                this.tetMesh.geometry.attributes.position.setXYZ(v, p.x, p.y, p.z);
            }
            this.tetMesh.geometry.attributes.position.needsUpdate = true;
            this.tetMesh.geometry.computeVertexNormals();            
        }

       
    }

    /*
    buildFaces()
    {
        // Deprecated way
        let material = new THREE.MeshLambertMaterial({ 
                    color: 0xffa94b, 
                    transparent: true, 
                    opacity: 0.30, 
                    side: THREE.DoubleSide, 
                    depthWrite: false,
                    flatShading : true
            });
        
        let geom = new THREE.Geometry();     
        geom.dynamic = true;            
        geom.vertices.push(this.get_vert(0));
        geom.vertices.push(this.get_vert(1));
        geom.vertices.push(this.get_vert(2));
        geom.vertices.push(this.get_vert(3));
        geom.faces.push( new THREE.Face3(0, 1, 2) );
        geom.faces.push( new THREE.Face3(0, 2, 3) );
        geom.faces.push( new THREE.Face3(0, 3, 1) );
        geom.faces.push( new THREE.Face3(1, 3, 2) );       
        this.tetMesh = new THREE.Mesh(geom, material);
    }
    */

    buildFaces() {
        let material = new THREE.MeshPhongMaterial ({
            color: 0xffa94b,
            transparent: true,
            opacity: 0.30,
            side: THREE.DoubleSide,
            depthWrite: true,
            flatShading: true
        });

        let geom = new THREE.BufferGeometry();
        geom.dynamic = true;
        let v0 = this.get_vert(0);
        let v1 = this.get_vert(1);
        let v2 = this.get_vert(2);
        let v3 = this.get_vert(3);
        const vertices = new Float32Array([
            v0.x, v0.y, v0.z,
            v1.x, v1.y, v1.z,
            v2.x, v2.y, v2.z,
            v3.x, v3.y, v3.z
        ]);

        geom.addAttribute('position', new THREE.BufferAttribute(vertices, 3));
        geom.setIndex([0, 1, 2,
            0, 2, 3,
            0, 3, 1,
            1, 3, 2]);
        
        this.tetMesh = new THREE.Mesh(geom, material);
        this.tetMesh.renderOrder = 1;
    }

    addToScene(scene) {
        for (let mesh of this.vertexMesh) {
            scene.add(mesh);
        }
        scene.add(this.tetMesh)
    }

    addDraggableObjects(list) {
        for (let mesh of this.vertexMesh) {
            list.push(mesh);
        }
    }

    minDistanceFrom(cog) {
        var minDist = cog.distanceTo(this.get_vert(0));
        for (let v = 1; v < 4; v++) {
            var dist = cog.distanceTo(this.get_vert(v));
            if (minDist > dist)
                minDist = dist;
        }
        return minDist;
    }



    circumcenter() {
        var e1 = new THREE.Vector3();
        e1.copy(this.get_vert(1));

        var e2 = new THREE.Vector3();
        e2.copy(this.get_vert(2));

        var e3 = new THREE.Vector3();
        e3.copy(this.get_vert(3));

        e1.sub(this.get_vert(0));
        e2.sub(this.get_vert(0));
        e3.sub(this.get_vert(0));

        var A = new THREE.Matrix3();
        A.set(
            e1.x, e1.y, e1.z,
            e2.x, e2.y, e2.z,
            e3.x, e3.y, e3.z);

        if (true) {

            var B = new THREE.Vector3();
            B.set(this.get_vert(1).lengthSq() - this.get_vert(0).lengthSq(),
                this.get_vert(2).lengthSq() - this.get_vert(0).lengthSq(),
                this.get_vert(3).lengthSq() - this.get_vert(0).lengthSq());

            B.multiplyScalar(0.5);           
            A.getInverse(A);
            return B.applyMatrix3(A);
        } 
        else 
        {
           // alternative:
           let a = cross(e2, e3).multiplyScalar( e1.lengthSq() );
           let b = cross(e3, e1).multiplyScalar( e2.lengthSq() );           
           let c = cross(e1, e2).multiplyScalar( e3.lengthSq() );
           
           const alpha = A.determinant() * 2.0;
           let sum = new THREE.Vector3(0,0,0);
           sum.add(a);
           sum.add(b);
           sum.add(c);
           sum.divideScalar( alpha );
           sum.add( this.get_vert(0) );
           return sum;
        }
    }

    static faceArea(a, b, c) {
        let edge0 = b.sub(a);
        let edge1 = c.sub(a);
        return (edge0.cross(edge1)).length() / 2.0;
    }

    area1() { return Tetrahedron.faceArea(this.get_vert(0), this.get_vert(1), this.get_vert(2)); }
    area2() { return Tetrahedron.faceArea(this.get_vert(0), this.get_vert(2), this.get_vert(3)); }
    area3() { return Tetrahedron.faceArea(this.get_vert(0), this.get_vert(3), this.get_vert(1)); }
    area4() { return Tetrahedron.faceArea(this.get_vert(1), this.get_vert(3), this.get_vert(2)); }
   

    sumArea(){
        return this.area1() + this.area2() + this.area3() + this.area4();
    }

    incenter() {
        let sum = this.sumArea();
        let incenter = new THREE.Vector3(0, 0, 0);
        incenter.add(this.get_vert(0).multiplyScalar(this.area4()));
        incenter.add(this.get_vert(1).multiplyScalar(this.area2()));
        incenter.add(this.get_vert(2).multiplyScalar(this.area3()));
        incenter.add(this.get_vert(3).multiplyScalar(this.area1()));
        incenter.divideScalar(sum);
        return incenter;
    }

    volume(){
        // to be checked
        // from wikipedia
        let a = this.get_vert(0);
        let b = this.get_vert(1);
        let c = this.get_vert(2);
        let d = this.get_vert(3);

        return Math.abs( (a.sub(d)).dot( (b.sub(d)).cross(c.sub(d)) ) ) / 6.0;
    }

    inradius() {

        if (false) {
            // https://people.sc.fsu.edu/~jburkardt/presentations/cg_lab_tetrahedrons.pdf sec. 14
            
            let a = this.get_vert(0);
            let b = this.get_vert(1);
            let c = this.get_vert(2);
            let d = this.get_vert(3);

            var A = new THREE.Matrix4();
            A.set(
                a.x, a.y, a.z, 1.0,
                b.x, a.y, b.z, 1.0,
                c.x, c.y, c.z, 1.0,
                d.x, d.y, d.z, 1.0);

            let alpha = A.determinant();
            // Notice that in the paper they compute the cross product of edges 
            // but don't divide by 2, effectively computing : area*2.0 
            return Math.abs(alpha) / (this.sumArea()*2.0);
        }
        else 
        {
            // from wikipedia
            return (this.volume() * 3.0) / this.sumArea();
        }
    }

    barycenter() {
        var cog = new THREE.Vector3(0, 0, 0);
        for (let v = 0; v < 4; v++) {
            cog.add(this.get_vert(v));
        }
        cog.divideScalar(4);
        return cog;
    }

    get_vert(i) {
        var vert = new THREE.Vector3();
        vert.copy(this.vertexMesh[i].position);
        return vert;
    }

    set_vert(i, pos) {
        this.vertexMesh[i].position = pos;
    }
}


function init() {

    const width = windowWidth();
    const height = windowHeight();

    g_camera = new THREE.PerspectiveCamera(45, width / height, 1, 2000);
    g_camera.position.set(50, 100, 150);

    g_scene = new THREE.Scene();
    g_scene.background = new THREE.Color(0xffffff);//0xa0a0a0
    g_scene.fog = new THREE.Fog(0xffffff, 200, 600);//0xa0a0a0

    var h_light = new THREE.HemisphereLight(0x555555, 0x555555);
    h_light.position.set(0, 200, 0);
    g_scene.add(h_light);

    var d_light = new THREE.DirectionalLight(0xffffff);
    d_light.position.set(0, 200, 100);
    d_light.castShadow = true;
    d_light.shadow.camera.top = 180;
    d_light.shadow.camera.bottom = - 100;
    d_light.shadow.camera.left = - 120;
    d_light.shadow.camera.right = 120;

    g_scene.add(d_light);

    //g_scene.add(new THREE.CameraHelper(light.shadow.camera));
    // ground
    var mesh = new THREE.Mesh(new THREE.PlaneBufferGeometry(2000, 2000), new THREE.MeshPhongMaterial({ color: 0xffffff, depthWrite: false }));
    mesh.rotation.x = - Math.PI / 2;
    mesh.receiveShadow = true;

    g_scene.add(mesh);

    var grid = new THREE.GridHelper(2000, 20, 0x000000, 0x000000);
    grid.material.opacity = 0.2;
    grid.material.transparent = true;

    g_scene.add(grid);

    var canvas = document.querySelector('#' + g_canvasName);
    g_renderer = new THREE.WebGLRenderer({ antialias: true, canvas });

    g_renderer.setSize(width, height);
    g_renderer.shadowMap.enabled = true;

    g_controls = new OrbitControls(g_camera, g_renderer.domElement);
    g_controls.target.set(0, 40, 0);
    g_controls.update();

    // base, position, axis, limits, size, graphicsOffset, material
    const w = 5.0;
    const h = 15.0;

    var draggableObjects = [];

    var sphereGeometry = new THREE.SphereBufferGeometry(1.0, 50, 50);
    var red = new THREE.MeshPhongMaterial({ color: 0xdd9933 });



    var invertEllipsoid = false;
    var material = new THREE.MeshLambertMaterial({ 
        color: 0x777777, 
        transparent: true, 
        opacity: 0.45, 
        side: invertEllipsoid ? THREE.BackSide : THREE.FrontSide, 
        depthWrite: false });
        

    g_tetrahedron = new Tetrahedron()
    g_tetrahedron.addToScene(g_scene)
    g_tetrahedron.addDraggableObjects(draggableObjects)

    // Set up transparent sphere
    g_fittedSphereMesh = new THREE.Mesh(sphereGeometry, material);
    g_fittedSphereMesh.scale.set(150, 150, 150);
    g_fittedSphereMesh.castShadow = false;
    g_fittedSphereMesh.renderOrder = 0;
    g_scene.add(g_fittedSphereMesh);


    g_cogPointMesh = new THREE.Mesh(sphereGeometry, red);
    g_scene.add(g_cogPointMesh);
    g_cogPointMesh.scale.set(2, 2, 2);
    g_cogPointMesh.position.set(0, 0, 0);
    g_cogPointMesh.castShadow = true;



    var dragControls = new DragControls(draggableObjects, g_camera, g_renderer.domElement);
    dragControls.addEventListener('dragstart', function () {
        g_controls.enabled = false;
    });
    dragControls.addEventListener('dragend', function () {
        g_controls.enabled = true;
    });

    window.addEventListener('resize', onWindowResize, false);
}

function onWindowResize() {
    const width = windowWidth();
    const height = windowHeight();
    g_camera.aspect = width / height;
    g_camera.updateProjectionMatrix();
    g_renderer.setSize(width, height);
}


function animate() {
    var cog;
    var radius = 1.0;
    g_fittedSphereMesh.visible = true;
    switch (g_datGuiContext.algoType) {
        case "Circumcenter": 
            cog = g_tetrahedron.circumcenter(); 
            radius = g_tetrahedron.minDistanceFrom(cog);
        break;

        case "Barycenter": 
            cog = g_tetrahedron.barycenter(); 
            radius = g_tetrahedron.minDistanceFrom(cog);
            g_fittedSphereMesh.visible = false;
        break;

        case "Incenter": 
            cog = g_tetrahedron.incenter(); 
            radius = g_tetrahedron.inradius();
        break;

        default: 
            cog = new THREE.Vector3(0, 0, 0); 
        break;
    }

    // update objects:
    g_cogPointMesh.position.copy(cog);
    g_fittedSphereMesh.position.copy(cog);   
    g_fittedSphereMesh.scale.set(radius, radius, radius);
    // Must be last:
    g_tetrahedron.updateGeom()

    requestAnimationFrame(animate);
    g_renderer.render(g_scene, g_camera);

}

init();
animate();

Javascript source code (Editable version)

Derivation

The circumcenter of a tetrahedron can be computed as intersection of three bisector planes. A bisector plane is defined as the plane centered on, and orthogonal to an edge of the tetrahedron. Alternatively, a perpendicular edge bisector can be defined as the set of points equidistant from two vertices of the tetrahedron.

Formula

A compact expression of the circumcenter \( \mathbf c \) of a tetrahedron with vertices \(\mathbf v_0, \mathbf v_1, \mathbf v_2, \mathbf v_3\) can be formulated as a matrix product:

\[ \mathbf c = \mathbf A^{-1} \mathbf B  \text{  with  }
\mathbf A =
\left(\begin{matrix}
\left[\mathbf v_1 - \mathbf v_0\right]^T \\
\left[\mathbf v_2 - \mathbf v_0\right]^T \\
\left[\mathbf v_3 - \mathbf v_0\right]^T
\end{matrix}\right)   \text{  and  }
\mathbf B =
\frac{1}{2}
\left( \begin{matrix}
\| \mathbf v_1 \|^2 - \| \mathbf v_0 \|^2 \\
\| \mathbf v_2 \|^2 - \| \mathbf v_0\|^2 \\
\| \mathbf v_3 \|^2 - \| \mathbf v_0\|^2
\end{matrix}\right) 
\]

Where:
\[ \vec e_{1} = \left[\mathbf v_1 - \mathbf v_0\right]^T = \{ v_{1x} - v_{0x}, v_{1y} - v_{0y}, v_{1z} - v_{0z} \} \\ \vec e_{2} = \left[\mathbf v_2 - \mathbf v_0\right]^T = \{ v_{2x} - v_{0x}, v_{2y} - v_{0y}, v_{2z} - v_{0z} \} \\ \vec e_{3} = \left[\mathbf v_3 - \mathbf v_0\right]^T = \{ v_{3x} - v_{0x}, v_{3y} - v_{0y}, v_{3z} - v_{0z} \}
\]

are the row vectors representing an edge of the tetrahedron and:

\[ \| \mathbf p \|^2 = p_x^2 + p_y^2 + p_z^2 \]

is the squared Euclidean norm.

In contrast to the centroid, the circumcenter may not always lay on the inside of a tetrahedron.
Analogously to an obtuse triangle, the circumcenter is outside of the object for an obtuse tetrahedron.


Expanded expression

Let's expand the above matrix formula in more details, here we use \( \vec e_{i} = \mathbf v_i - \mathbf v_0 \) for brievity and \( \times \) is the usual vector cross-product:
\[ c = \mathbf v_0 + \frac{ 1 } {2 ~ det \left| \mathbf A \right |} \left ( \| \vec e_{3} \|^2 \big ( \vec e_{1} \times \vec e_{2} \big ) + \| \vec e_{2} \|^2 \big ( \vec e_{3} \times \vec e_{1} \big ) + \| \vec e_{1} \|^2 \big ( \vec e_{2} \times \vec e_{3} \big ) \right ) \]

The determinant of the \( 3 \times 3 \) \( \mathbf A \) matrix can be expanded as well recursively:

\[ \begin{align*} det \left| \mathbf A \right | &= \left| \begin{matrix} v_{1x} - v_{0x} & v_{1y} - v_{0y} & v_{1z} - v_{0z} \\ v_{2x} - v_{0x} & v_{2y} - v_{0y} & v_{2z} - v_{0z} \\ v_{3x} - v_{0x} & v_{3y} - v_{0y} & v_{3z} - v_{0z} \end{matrix} \right | \\ ~\\ &= \begin{vmatrix} e_{1x} & e_{1y} & e_{1z} \\ e_{2x} & e_{2y} & e_{2z} \\ e_{3x} & e_{3y} & e_{3z} \end{vmatrix} \\~\\ &= e_{1x} \begin{vmatrix} e_{2y} & e_{2z} \\ e_{3y} & e_{3z} \end{vmatrix} - e_{1y} \begin{vmatrix} e_{2x} & e_{2z} \\ e_{3x} & e_{3z} \end{vmatrix} + e_{1z} \begin{vmatrix} e_{2x} & e_{2y} \\ e_{3x} & e_{3y} \end{vmatrix}\\ ~\\ &= e_{1x}e_{2y}e_{3z} + e_{1y}e_{2z}e_{3x} + e_{1z}e_{2x}e_{3y} - e_{1z}e_{2y}e_{3x} - e_{1y}e_{2x}e_{3z} - e_{1x}e_{2z}e_{3y} \end{align*} \]

C/C++ implementation

Notes on stability: these expressions are unstable only in one case: if the denominator is close to zero. This instability, which arises if the tetrahedron is nearly degenerate, is unavoidable. Depending on your application, you may want to use exact arithmetic to compute the value of the determinant.

Inputs are the vertices coordinates 'a[]' 'b[]' 'c[]' and 'd[]'. The result is returned both in terms of x y z coordinates in circumcenter[] and xi-eta-zeta coordinates. Both outputs are relative to the tetrahedron's vertex 'a' (that is, 'a' is the origin of both coordinate systems). In other words, the x y z coordinates returned are NOT absolute! One must add the coordinates of 'a' to find the absolute coordinates of the circumcircle.

The xi-eta-zeta are coordinate along the tetrahedron's edges. Vertex 'a' being the origin, xi axis represent one unit along the edge 'ab'. One unit of the eta axis is the edge 'ac'  and logically zeta is the coordinate along the edge 'ad'. These coordinate values are useful for linear interpolation.

Note: if 'xi' is NULL on input, the xi-eta-zeta coordinates will not be computed.

void tetrahedron_circumcenter(
        // In:
        const double a[3],
        const double b[3],
        const double c[3],
        const double d[3],
        // Out:
        double circumcenter[3],
        double* xi,
        double* eta,
        double* zeta)
{
    double denominator;


    // Use coordinates relative to point 'a' of the tetrahedron.

    // ba = b - a
    double ba_x = b[0] - a[0];
    double ba_y = b[1] - a[1];
    double ba_z = b[2] - a[2];

    // ca = c - a
    double ca_x = c[0] - a[0];
    double ca_y = c[1] - a[1];
    double ca_z = c[2] - a[2];

    // da = d - a
    double da_x = d[0] - a[0];
    double da_y = d[1] - a[1];
    double da_z = d[2] - a[2];

    // Squares of lengths of the edges incident to 'a'.
    double len_ba = ba_x * ba_x + ba_y * ba_y + ba_z * ba_z;
    double len_ca = ca_x * ca_x + ca_y * ca_y + ca_z * ca_z;
    double len_da = da_x * da_x + da_y * da_y + da_z * da_z;

    // Cross products of these edges.

    // c cross d
    double cross_cd_x = ca_y * da_z - da_y * ca_z;
    double cross_cd_y = ca_z * da_x - da_z * ca_x;
    double cross_cd_z = ca_x * da_y - da_x * ca_y;

    // d cross b
    double cross_db_x = da_y * ba_z - ba_y * da_z;
    double cross_db_y = da_z * ba_x - ba_z * da_x;
    double cross_db_z = da_x * ba_y - ba_x * da_y;

    // b cross c
    double cross_bc_x = ba_y * ca_z - ca_y * ba_z;
    double cross_bc_y = ba_z * ca_x - ca_z * ba_x;
    double cross_bc_z = ba_x * ca_y - ca_x * ba_y;

    // Calculate the denominator of the formula.
    denominator = 0.5 / (ba_x * cross_cd_x + ba_y * cross_cd_y + ba_z * cross_cd_z);

    // Calculate offset (from 'a') of circumcenter.
    double circ_x = (len_ba * cross_cd_x + len_ca * cross_db_x + len_da * cross_bc_x) * denominator;
    double circ_y = (len_ba * cross_cd_y + len_ca * cross_db_y + len_da * cross_bc_y) * denominator;
    double circ_z = (len_ba * cross_cd_z + len_ca * cross_db_z + len_da * cross_bc_z) * denominator;

    circumcenter[0] = circ_x;
    circumcenter[1] = circ_y;
    circumcenter[2] = circ_z;

    if (xi != (double*) nullptr) {
        // To interpolate a linear function at the circumcenter, define a
        // coordinate system with a xi-axis directed from 'a' to 'b',
        // an eta-axis directed from 'a' to 'c', and a zeta-axis directed
        // from 'a' to 'd'.  The values for xi, eta, and zeta are computed
        // by Cramer's Rule for solving systems of linear equations.
        denominator *= 2.0;
        *xi   = (circ_x * cross_cd_x + circ_y * cross_cd_y + circ_z * cross_cd_z) * denominator;
        *eta  = (circ_x * cross_db_x + circ_y * cross_db_y + circ_z * cross_db_z) * denominator;
        *zeta = (circ_x * cross_bc_x + circ_y * cross_bc_y + circ_z * cross_bc_z) * denominator;
    }
}

References

Find C/C++ code for the circumcenter and circumradius of tetrahedrons, 3D and 2D triangles here:

Other:

No comments

(optional field, I won't disclose or spam but it's necessary to notify you if I respond to your comment)
All html tags except <b> and <i> will be removed from your comment. You can make links by just typing the url or mail-address.
Anti-spam question: