四元数和旋转轴

4
这个问题已经得到回答。以下是大部分代码,您需要使用它来使其工作!希望能帮助其他人。
感谢@Aki Suihkonen,@David Hammen和@MBo。
两个角函数都给出了正确的答案。
我有三个点:
A: 12 4 5
B: 6 8 -10
C: 5 6 7

我已经实现了四元数。我想旋转点C,使得角度(A,B,C)比之前高40度。
我的问题是:根据哪个轴线进行旋转?我想象因为A、B和C创建了一个平面,所以我必须根据向量BA和BC的垂直轴对点C进行旋转。我用它们单位向量的叉积得到它,但是当我尝试获取角度(A,B,C)时,它并没有给我正确的结果。
这是我获得角度的方法:(旧方法)
results:
Angle of ABC before rotation = 67.3895.
Angle of ABC after rotation = 107.389.

float Angle( float x1, float y1, float z1,
             float x2, float y2, float z2 )
{
  float x, y, z;
  CrossProduct( x1, y1, z1, x2, y2, z2, &x, &y, &z );

  float result = atan2 ( L2Norm( x, y, z ),
                         DotProduct( x1, y1, z1, x2, y2, z2 ) );

  return result;
}

当 x1、y1、z1、x2、y2 和 z2 分别为矢量 B-A 和 B-C 的单位向量时,求出的结果为:

更新后的角度函数:

results:
Angle of ABC before rotation = 67.3895.
Angle of ABC after rotation = 107.389.

Angle( Atom &atom1, Atom &atom2, Atom &atom3 )
{
float px1, py1, pz1;
    float px2, py2, pz2;

    UnitVector( atom1.X(), atom1.Y(), atom1.Z(),
            atom2.X(), atom2.Y(), atom2.Z(),
            &px1, &py1, &pz1 );


    UnitVector( atom3.X(), atom3.Y(), atom3.Z(),
            atom2.X(), atom2.Y(), atom2.Z(),
            &px2, &py2, &pz2 );

  float dot_product = DotProduct( px1, py1, pz1, px2, py2, pz2 );
  float length_BA = sqrt( px1*px1 + py1*py1 + pz1*pz1 );
  float length_BC = sqrt( px2*px2 + py2*py2 + pz2*pz2 );

  return acos( dot_product / ( length_BA * length_BC ) );
}

float DotProduct( float x1, float y1, float z1,
                  float x2, float y2, float z2 )
{
  return x1*x2 + y1*y2 + z1*z2;
}

void CrossProduct( float x1, float y1, float z1,
                   float x2, float y2, float z2,
                   float *ox, float *oy, float *oz )
{
  *ox = (y1*z2) -(z1*y2);
  *oy = -((x1*z2) -(z1*x2));
  *oz = (x1*y2) -(y1*x2);
}

所以我的问题是:我需要将点C绕哪个轴旋转,才能使角度(A,B,C)比之前大40度?
  // The 3 points.
  Atom A( "A", -4, 2 , 8 );
  Atom B( "B", -1, 3 , 4 );
  Atom C( "C", -2, -4 , 5 );

  float x1, y1, z1;
  float x2, y2, z2;
  float x, y, z;

  // Get the cross product. Create the perpendicular vector to the BA and BC vectors.
  PointVector( A.X(), A.Y(), A.Z(), B.X(), B.Y(), B.Z(), &x1, &y1, &z1 );
  PointVector( C.X(), C.Y(), C.Z(), B.X(), B.Y(), B.Z(), &x2, &y2, &z2 );

  CrossProduct( x1, y1, z1, x2, y2, z2, &x, &y, &z );

  // Normalize the coordinates.
  float length = sqrt( x*x + y*y + z*z );
  length = 1.0 / length;

  x *= length;
  y *= length;
  z *= length;

  // Create the 40 degrees angle. It is supposed to increment the current ABC angle by 40 degrees.
  float angle = 40*M_PI/180;
  float sinAngleOver2 = sin(angle/2);
  float w = cos(angle/2);

  // Create the axis quat.
  Quatd q(w, x * sinAngleOver2, y * sinAngleOver2, z * sinAngleOver2);

  // Create the point quaternal. The angle of it equals to the current ABC angle.
  angle = Angle( A, B, C ) *180/M_PI;
  angle *= M_PI/180;
  sinAngleOver2 = sin(angle/2);
  w = cos(angle/2);

  // Normalize the coordinates. The coordinates are the C point coordinates.


    x = C.X() - B.X();
    y = C.Y() - B.Y();
    z = C.Z() - B.Z();

    length = sqrt( x*x + y*y + z*z );
    length = 1.0 / length;

    x *= length;
    y *= length;
    z *= length;


  Quatd qpt(w, x * sinAngleOver2, y * sinAngleOver2, z * sinAngleOver2);

  // Rotate.
  qpt = q*qpt*q.unitInverse();

  Atom new_C;
  new_C.X( qpt.x + B.X() );
  new_C.Y( qpt.y + B.Y() );
  new_C.Z( qpt.z + B.Z() );

  cout << "Angle: " << Angle( A, B, new_C )*180/M_PI << '\n';

1
函数L2Norm是做什么用的?你可以通过计算结果=arccos(DotProduct(x1, y1, z1, x2, y2, z2)/(BA.Length*CA.Length))来计算角度。 - MBo
不确定。atan2公式是这样解释的。我相信我在Matlab论坛中看到过它。至于轴,你有解决方案吗? - user985611
是的,旋转轴将是平面法线,它是定义平面的两个向量的归一化叉乘。然而,正如MBo(以及您在其他主题中的我)所建议的那样,您的Angle函数可能需要修订。它给出89.926,而arccos给出80.841(度)。 - Aki Suihkonen
我指的是B-A和C-A,因为你使用了“向量B-A和C-A”,而我认为A是旋转点。如果使用单位向量,请忘记分母。 - MBo
我更新了以反映:其中x1、y1、z1、x2、y2、z2是单位向量B-A和B-C的结果。 - user985611
显示剩余10条评论
1个回答

3
这段代码是不正确的:
// Normalize the coordinates. The coordinates are the C point coordinates.
length = sqrt( C.X()*C.X() + C.Y()*C.Y() + C.Z()*C.Z() );
length = 1.0 / length;

x = C.X() / length;
y = C.Y() / length;
z = C.Z() / length;

您已经混淆了自己,通过重新使用您的“长度”变量来存储“1 /长度”。在此处,您需要将C的组件乘以长度来规范化向量。
我不熟悉四元数数学,无法理解您在末尾尝试做什么,但似乎您正在围绕原点旋转C。您想要围绕B旋转,这意味着从C中减去B,围绕原点进行旋转,然后将结果添加到B上以返回到原始空间。
个人而言,我会实现某种四元数-向量乘法,或将四元数转换为矩阵以执行变换。

胜利了!非常感谢。现在它可以工作了。我将使用正确的代码更新主贴,以便每个人都能受益。如果这是 Reddit,我会给每个人点赞! - user985611

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接