Xây dựng vòng tròn từ 3 điểm trong triển khai không gian 3D trong C hoặc C++

Question

Chúng tôi có 3 (ba) điểm xyz xác định vòng tròn trong không gian 3D, vòng tròn này cần phải được chuyển đổi thành polyline (để hiển thị thêm). Tôi đang tìm một hàm C hoặc C++ đã sẵn sàng hoặc thư viện miễn phí có thể thực hiện công việc.

Không hiểu tại sao this bị đóng. Và tôi thậm chí không thể trả lời câu hỏi của riêng tôi ở đó. Xấu hổ với các bạn. Nhưng bạn sẽ không ngăn chặn sự lan truyền tri thức!

Answer 1

Có một bài viết hay và một mẫu mã về cách xây dựng một vòng tròn bằng 3 điểm trong mặt phẳng 2D, XY.

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

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

Để xây dựng một vòng tròn 3D chúng tôi sẽ phải:

• xoay 3 điểm của chúng tôi vào XY máy bay
trung tâm vòng tròn
• Tính
• xây dựng một vòng tròn trên mặt phẳng XY sử dụng mã trong bài viết
• xoay nó trở lại vào mặt phẳng ban đầu của nó

Để quay tốt nhất là sử dụng quaternions.

Để tìm một quaternion đúng Tôi nhìn mã nguồn Ogre3d: void Quaternion::FromAngleAxis (const Radian& rfAngle, const Vector3& rkAxis)

Có thêm một chức năng hữu ích có: Quaternion getRotationTo(const Vector3& dest, const Vector3& fallbackAxis = Vector3::ZERO) const Nhưng tôi đã không sử dụng nó.

Đối với quaterions và vectơ tôi đã sử dụng các lớp của riêng mình. Đây là mã nguồn đầy đủ các chức năng mà không được công việc:

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ó một giải pháp đơn giản hơn nhiều để tìm các thông số vòng tròn trong không gian 3D thực tế, chỉ có một cái nhìn tại "tọa barycentric" trong http://en.wikipedia.org/wiki/Circumscribed_circle. Bạn có thể trích xuất mã tối ưu hóa sau đây từ đó:

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

Sau đó bạn có thể tính toán các điểm trên vòng tròn trong không gian 3D thực quá và ví dụ vẽ chúng bằng GL_LINE_STRIP trong OpenGL. Điều này sẽ nhanh hơn nhiều so với cách sử dụng phương pháp tiếp cận 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 
} 

Để tạo ma trận chuyển đổi (ví dụ: makeRotMatrix4d()) xem http://paulbourke.net/geometry/rotate/ chẳng hạn.

Xin lưu ý rằng tôi đã không kiểm tra xem mã trên có thực sự biên dịch hay không, nhưng nó sẽ cung cấp cho bạn đủ gợi ý.