Build-Kreis von 3 Punkten in 3D-Raum Implementierung in C oder C++

Question

Wir haben 3 (drei) xyz Punkte, die einen Kreis im 3D-Raum definieren, dieser Kreis muss in eine Polylinie umgewandelt werden (für weitere Rendering). Ich suche nach einer fertigen C- oder C++ - Funktion oder einer freien Bibliothek, die diese Aufgabe erfüllen kann.

Verstehe nicht, warum this geschlossen war. Und ich kann dort meine eigene Frage nicht einmal beantworten. Schande über euch Leute. Aber Sie werden nicht aufhören, das Wissen zu verbreiten!

Answer 1

Es gibt einen schönen Artikel und ein Codebeispiel zum Erstellen eines Kreises um 3 Punkte in 2D, XY Ebene.

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

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

Um einen 3D-Kreis bauen wir müssen:

• unsere 3 Punkte in XY-Ebene drehen
• berechnen Kreismittel
• einen Kreis bauen in XY-Ebene unter Verwendung des Codes in dem Artikel
• drehen Sie es zurück in seine ursprüngliche Ebene

Für Rotationen ist es am besten, Quaternionen zu verwenden.

Um eine korrekte quaternion finde ich bei Ogre3D Quellcode sah: void Quaternion::FromAngleAxis (const Radian& rfAngle, const Vector3& rkAxis)

Es gibt eine weitere nützliche Funktion gibt: Quaternion getRotationTo(const Vector3& dest, const Vector3& fallbackAxis = Vector3::ZERO) const Aber ich habe es nicht verwenden.

Für Quaterionen und Vektoren habe ich unsere eigenen Klassen verwendet. Hier ist der vollständige Quellcode der Funktion, die die Arbeit erledigt:

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

Es gibt eine viel einfachere Lösung, die die Kreisparameter in Echtzeit 3D zu finden, nehmen Sie nur einen Blick auf die „Schwerpunktkoordinaten“ in http://en.wikipedia.org/wiki/Circumscribed_circle. Sie können den folgenden optimierten Code aus, dass extrahieren:

// 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); 

Sie können dann die Punkte auf dem Kreis in echten 3D-Berechnung zu und z.B. zeichne sie mit GL_LINE_STRIP in OpenGL. Dies sollte viel schneller sein als der 2D-sin/cos-Ansatz.

// 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 
} 

Für die Erstellung der Transformations-Matrix (d.h. makeRotMatrix4d()), siehe beispielsweise http://paulbourke.net/geometry/rotate/.

Bitte beachten Sie, dass ich nicht getestet habe, ob der obige Code wirklich kompiliert, aber es sollte Ihnen genug Hinweise geben.