logo标签列表

在C或C++中实现在三维空间中通过3个点构建圆的方法

c++geometry
6
我们有3个(xyz)点在三维空间中定义了一个圆，需要将这个圆转换成折线(以便进一步渲染)。我正在寻找一个现成的C或C++函数或免费库来完成这项工作。
不理解为什么this被关闭。我甚至无法在那里回答自己的问题。你们真可耻。但你们不会阻止知识的传播！
我不明白为什么你不理解为什么你的问题被关闭了。你已经收到了两位发帖人的评论来解释原因。你可能在这里待了足够长的时间，看到了数十个“请给我代码”的问题，所以你应该明白为什么社区会想关闭任何类似的问题。请不要假设关闭投票者是出于恶意行事，或者为了“阻止知识传播”。 - Kevin
1
不，我没有看到这样的问题。相反，我看到了大量要求代码并且有带有代码或链接的答案的问题。所以我不明白为什么这个问题与众不同？如果我表述问题不好，请帮我重新表述。 - Dmitriy
3个回答
12

寻找真实3D圆形参数的更简单解决方案，只需查看http://en.wikipedia.org/wiki/Circumscribed_circle中的“重心坐标”部分即可。您可以从中提取以下优化代码：

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

你可以使用OpenGL中的GL_LINE_STRIP将真正的三维圆上的点计算出来并绘制它们。 这比使用二维sin / cos方式要快得多。

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

创建变换矩阵（即makeRotMatrix4d()）请参见http://paulbourke.net/geometry/rotate/。

请注意，我没有测试上述代码是否真正编译通过，但它应该给您足够的提示。

使用更直接的重心坐标实现，我得到了正确的答案，但是我无法让您的“优化”解决方案正常工作。您能否更详细地展示它们的推导过程？我假设（u*v）表示点积，但这又如何从方程中推出呢？ - Quantum7
@Mark wsl和iwsl2是什么？ - David Doria
@Mark，const double tt = t*t; 中的符号 t*t 是指 t.dot(t) 吗？ - David Doria
1
哈，这是一个很酷的解决方案。谢谢！ :) 我已经将原始代码移植到C#（这导致了一些对t*t和类似计算的混淆），并将其发布为附加答案。提示：确实意味着点积。 - j00hi
6
以下是Mark的answer的C#/Unity移植版。它使用了Unity脚本API中的类型和实用函数。
// triangle "edges"
var t = p2 - p1;
var u = p3 - p1;
var v = p3 - p2;

// triangle normal
var w = Vector3.Cross(t, u);
var wsl = Vector3.Dot(w, w);
// TODO: 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
var iwsl2 = 1f / (2f * wsl);
var tt = Vector3.Dot(t, t);
var uu = Vector3.Dot(u, u);

// result circle
Vector3 circCenter = p1 + (u * tt * (Vector3.Dot(u, v)) - t * uu * (Vector3.Dot(t, v))) * iwsl2;
var     circRadius = Mathf.Sqrt(tt * uu * (Vector3.Dot(v, v)) * iwsl2 * 0.5f);
Vector3 circAxis   = w / Mathf.Sqrt(wsl);

使用 Unity的Gizmos，可以如下绘制圆形（在这种情况下使用30条线段来近似它）：
// Draw the circle:
Gizmos.color = Color.white;
for (int i = 0; i < 30; ++i) 
{
    Gizmos.DrawLine(
        circCenter + Quaternion.AngleAxis(360f / 30f *  i     , circAxis) * (p1 - circCenter),
        circCenter + Quaternion.AngleAxis(360f / 30f * (i + 1), circAxis) * (p1 - circCenter)
    );
}

以下是顶点位置的结果：var p1 = new Vector3(0f, 1.44f, 0f); var p2 = new Vector3(0f, 0.73f, 0.65f); var p3 = new Vector3(0f, -1.04f, 0f);

Circle through three points

这就是我一直在寻找的。我最初尝试了这个答案，但它没有起作用。 - John Alexiou
3

这篇文章介绍了如何在二维平面中通过三个点构建圆的优美代码示例。

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

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

要在三维空间中构建一个圆，我们需要：

  • 将三个点旋转到二维平面中
  • 计算圆心
  • 使用文章中的代码在二维平面上构建圆
  • 将圆旋转回原来的平面

对于旋转，最好使用四元数。

为了找到正确的四元数，我查看了Ogre3d源代码： void Quaternion::FromAngleAxis (const Radian& rfAngle, const Vector3& rkAxis)

那里还有一个更有用的函数： Quaternion getRotationTo(const Vector3& dest, const Vector3& fallbackAxis = Vector3::ZERO) const 但我没有使用它。

对于四元数和向量，我使用了我们自己的类。以下是执行此操作的函数的完整源代码：

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);
    }
}
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接