геометрия — Построить круг из 3-х точек в трехмерной космической реализации в C или Stack Overflow

Есть хорошая статья и пример кода о том, как построить круг на 3 точки в плоскости 2D, XY.

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

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

Чтобы построить трехмерный круг, нам нужно:

• поверните наши 3 точки в плоскости XY
• Рассчитать центр круга
• построить круг в плоскости XY, используя код в статье
• поверните его обратно в исходную плоскость

Для вращений лучше всего использовать кватернионы.

Чтобы найти правильный кватернион, я посмотрел исходный код 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);
}
}