Crea cerchio da 3 punti nell'implementazione dello spazio 3D in C o C++

Question

Abbiamo 3 (tre) punti xyz che definiscono un cerchio nello spazio 3D, questo cerchio deve essere convertito in una polilinea (per ulteriori rendering). Sto cercando una funzione C o C++ pronta o una libreria libera che possa fare il lavoro.

Non capisco perché this era chiuso. E non posso nemmeno rispondere alla mia domanda lì. Vergognatevi ragazzi. Ma non fermerai la diffusione della conoscenza!

Answer 1

C'è un bell'articolo e un esempio di codice su come costruire un cerchio di 3 punti in 2D, piano XY.

http://paulbourke.net/geometry/circlesphere/

http://paulbourke.net/geometry/circlesphere/Circle.cpp

Per costruire un cerchio 3D dovremo:

• ruotare i nostri 3 punti in XY
• Calcolare cerchio centrale
• costruire un cerchio nel piano XY utilizzando il codice nell'articolo
• ruotarlo di nuovo nel suo piano originale

Per le rotazioni è preferibile utilizzare i quaternioni.

Per trovare un corretto quaternione ho guardato il codice sorgente Ogre3D: void Quaternion::FromAngleAxis (const Radian& rfAngle, const Vector3& rkAxis)

C'è un'altra funzione utile lì: Quaternion getRotationTo(const Vector3& dest, const Vector3& fallbackAxis = Vector3::ZERO) const ma non l'ho usato.

Per quaterioni e vettori ho usato le nostre classi. Ecco il codice sorgente completo della funzione che fa il lavoro:

bool IsPerpendicular(Point3d *pt1, Point3d *pt2, Point3d *pt3); 
double CalcCircleCenter(Point3d *pt1, Point3d *pt2, Point3d *pt3, Point3d *center); 

void FindCircleCenter(const Point3d *V1, const Point3d *V2, const Point3d *V3, Point3d *center) 
{ 
    Point3d *pt1=new Point3d(*V1); 
    Point3d *pt2=new Point3d(*V2); 
    Point3d *pt3=new Point3d(*V3); 


    if (!IsPerpendicular(pt1, pt2, pt3))  CalcCircleCenter(pt1, pt2, pt3, center); 
    else if (!IsPerpendicular(pt1, pt3, pt2)) CalcCircleCenter(pt1, pt3, pt2, center); 
    else if (!IsPerpendicular(pt2, pt1, pt3)) CalcCircleCenter(pt2, pt1, pt3, center); 
    else if (!IsPerpendicular(pt2, pt3, pt1)) CalcCircleCenter(pt2, pt3, pt1, center); 
    else if (!IsPerpendicular(pt3, pt2, pt1)) CalcCircleCenter(pt3, pt2, pt1, center); 
    else if (!IsPerpendicular(pt3, pt1, pt2)) CalcCircleCenter(pt3, pt1, pt2, center); 
    else { 
     delete pt1; 
     delete pt2; 
     delete pt3; 
     return; 
    } 
    delete pt1; 
    delete pt2; 
    delete pt3; 

} 

bool IsPerpendicular(Point3d *pt1, Point3d *pt2, Point3d *pt3) 
// Check the given point are perpendicular to x or y axis 
{ 
    double yDelta_a= pt2->y - pt1->y; 
    double xDelta_a= pt2->x - pt1->x; 
    double yDelta_b= pt3->y - pt2->y; 
    double xDelta_b= pt3->x - pt2->x; 

    // checking whether the line of the two pts are vertical 
    if (fabs(xDelta_a) <= 0.000000001 && fabs(yDelta_b) <= 0.000000001){ 
     return false; 
    } 

    if (fabs(yDelta_a) <= 0.0000001){ 
     return true; 
    } 
    else if (fabs(yDelta_b) <= 0.0000001){ 
     return true; 
    } 
    else if (fabs(xDelta_a)<= 0.000000001){ 
     return true; 
    } 
    else if (fabs(xDelta_b)<= 0.000000001){ 
     return true; 
    } 
    else 
     return false ; 
} 

double CalcCircleCenter(Point3d *pt1, Point3d *pt2, Point3d *pt3, Point3d *center) 
{ 
    double yDelta_a = pt2->y - pt1->y; 
    double xDelta_a = pt2->x - pt1->x; 
    double yDelta_b = pt3->y - pt2->y; 
    double xDelta_b = pt3->x - pt2->x; 

    if (fabs(xDelta_a) <= 0.000000001 && fabs(yDelta_b) <= 0.000000001){ 
     center->x= 0.5*(pt2->x + pt3->x); 
     center->y= 0.5*(pt1->y + pt2->y); 
     center->z= pt1->z; 

     return 1; 
    } 

    // IsPerpendicular() assure that xDelta(s) are not zero 
    double aSlope=yDelta_a/xDelta_a; // 
    double bSlope=yDelta_b/xDelta_b; 
    if (fabs(aSlope-bSlope) <= 0.000000001){ // checking whether the given points are colinear. 
     return -1; 
    } 

    // calc center 
    center->x= (aSlope*bSlope*(pt1->y - pt3->y) + bSlope*(pt1->x + pt2 ->x) 
         - aSlope*(pt2->x+pt3->x))/(2* (bSlope-aSlope)); 
    center->y = -1*(center->x - (pt1->x+pt2->x)/2)/aSlope + (pt1->y+pt2->y)/2; 

    return 1; 
} 

//! Builds a circle in 3D space by 3 points on it and an optional center 
void buildCircleBy3Pt(const float *pt1, 
         const float *pt2, 
         const float *pt3, 
         const float *c,  // center, can be NULL 
         std::vector<float> *circle) 
{ 
    /* Get the normal vector to the triangle formed by 3 points 
     Calc a rotation quaternion from that normal to the 0,0,1 axis 
     Rotate 3 points using quaternion. Points will be in XY plane 
     Build a circle by 3 points on XY plane 
     Rotate a circle back into original plane using quaternion 
    */ 
    Point3d p1(pt1[0], pt1[1], pt1[2]); 
    Point3d p2(pt2[0], pt2[1], pt2[2]); 
    Point3d p3(pt3[0], pt3[1], pt3[2]); 
    Point3d center; 
    if (c) 
    { 
     center.set(c[0], c[1], c[2]); 
    } 

    const Vector3d p2top1 = p1 - p2; 
    const Vector3d p2top3 = p3 - p2; 

    const Vector3d circle_normal = p2top1.crossProduct(p2top3).normalize(); 
    const Vector3d xy_normal(0, 0, 1); 


    Quaternion rot_quat; 
    // building rotation quaternion 
    { 
     // Rotation axis around which we will rotate our circle into XY plane 
     Vector3d rot_axis = xy_normal.crossProduct(circle_normal).normalize(); 
     const double rot_angle = xy_normal.angleTo(circle_normal); // radians 

     const double w = cos(rot_angle * 0.5); 
     rot_axis *= sin(rot_angle * 0.5); 

     rot_quat.set(w, rot_axis.x, rot_axis.y, rot_axis.z); 
    } 

    Quaternion rot_back_quat; 
    // building backward rotation quaternion, same as prev. but -angle 
    { 
     const double rot_angle = -(xy_normal.angleTo(circle_normal)); // radians 
     const double w_back = cos(rot_angle * 0.5); 
     Vector3d rot_back_axis = xy_normal.crossProduct(circle_normal).normalize(); 
     rot_back_axis *= sin(rot_angle * 0.5); 
     rot_back_quat.set(w_back, rot_back_axis.x, rot_back_axis.y, rot_back_axis.z); 
    } 

    rot_quat.rotate(p1); 
    rot_quat.rotate(p2); 
    rot_quat.rotate(p3); 
    rot_quat.rotate(center); 

    if (!c) 
    { 
     // calculate 2D center 
     FindCircleCenter(&p1, &p2, &p3, &center); 
    } 

    // calc radius 
    const double radius = center.distanceTo(p1); 

    const float DEG2RAD = 3.14159f/180.0f; 
    // build circle 
    for (int i = 0; i < 360; ++i) 
    { 
     float degInRad = i * DEG2RAD; 
     Point3d pt(cos(degInRad) * radius + center.x, sin(degInRad) * radius + center.y, 0); 

     // rotate the point back into original plane 
     rot_back_quat.rotate(pt); 

     circle->push_back(pt.x); 
     circle->push_back(pt.y); 
     circle->push_back(pt.z); 
    } 
} 

Answer 2

C'è una soluzione molto più semplice per trovare i parametri cerchio in 3D reale, basta dare un'occhiata alla sezione "coordinate baricentriche" in http://en.wikipedia.org/wiki/Circumscribed_circle. È possibile estrarre il seguente codice ottimizzato da questo:

// triangle "edges" 
const Vector3d t = p2-p1; 
const Vector3d u = p3-p1; 
const Vector3d v = p3-p2; 

// triangle normal 
const Vector3d w = t.crossProduct(u); 
const double wsl = w.getSqrLength(); 
if (wsl<10e-14) return false; // area of the triangle is too small (you may additionally check the points for colinearity if you are paranoid) 

// helpers 
const double iwsl2 = 1.0/(2.0*wsl); 
const double tt = t*t; 
const double uu = u*u; 

// result circle 
Vector3d circCenter = p1 + (u*tt*(u*v) - t*uu*(t*v)) * iwsl2; 
double circRadius = sqrt(tt * uu * (v*v) * iwsl2*0.5); 
Vector3d circAxis = w/sqrt(wsl); 

È quindi possibile calcolare i punti sul cerchio in 3D reale troppo e per esempio disegnali usando GL_LINE_STRIP in OpenGL. Questo dovrebbe essere molto più veloce rispetto all'utilizzo dell'approccio sin/cos 2D.

// find orthogonal vector to the circle axis 
const Vector3d an = circAxis.getNormalized(); 
const Vector3d ao = Vector3d(4.0+an[0], 4.0+an[0]+an[1], 4.0+an[0]+an[1]+an[2]).crossProduct(an).getNormalized(); 

// 4x4 rotation matrix around the circle axis 
const int steps = 360; // maybe adjust according to circle size on screen 
Matrix4d R = makeRotMatrix4d(circCenter, circAxis, 2.0*M_PI/double(steps)); 

// one point on the circle 
Vector3d cp = circCenter + ao*circRadius; 

// rotate point on the circle 
for (int i=0; i<steps; ++i) 
{ 
    circlePoints.push_back(cp); 
    cp = transformPoint(cp, R); // apply the matrix 
} 

Per la realizzazione della matrice di trasformazione (cioè makeRotMatrix4d()) vedi http://paulbourke.net/geometry/rotate/ per esempio.

Si prega di notare che non ho provato se il codice sopra compila davvero, ma dovrebbe darti abbastanza suggerimenti.