8 replies to this topic

[Gasim]

Hello, I want to ask a question about Physics and Math together: 1. Can i use verlet integration with Orientation. Something Like that:

orientation += (orientation - prev_orientation) + angular_acceleration * dt * dt;

Orientation is Quaternion Angular Acceleration is Vector3 2. If i can do the 1st: Can i convert Angular acceleration to Quaternion using Euler to Quaternion? Like That:

double dSY = sin(nz * 0.5f);
double dSP = sin(ny * 0.5f);
double dSR = sin(nx * 0.5f);
double dCY = cos(nz * 0.5f);
double dCP = cos(ny * 0.5f);
double dCR = cos(nx * 0.5f);

x = dSR * dCP * dCY - dCR * dSP * dSY;
y = dCR * dSP * dCY + dSR * dCP * dSY;
z = dCR * dCP * dSY - dSR * dSP * dCY;
w = dCR * dCP * dCY + dSR * dSP * dSY;

nx, nz, ny are coordinates of Vector3 Thanks, Kasya [Edited by - Kasya on December 13, 2008 7:52:52 AM]

[alvaro]

I would probably stay away from verlet integration for rotations. When your variable is a position, the previous position captures essentially the same information as the velocity. However, when dealing with rotations, if the previous rotation and the current rotation are the same, you can't know whether the object is not moving or is rotating at one turn per frame.

Other than that, your formulas seem quite wrong. Your proposed formula for verlet integration wouldn't work for positions either.

Quaternions are a good representation for orientations, but not for angular velocities and angular accelerations. Those are properly described by pseudo-vectors.

[Gasim]

Quote:

I would probably stay away from verlet integration for rotations. When your variable is a position, the previous position captures essentially the same information as the velocity. However, when dealing with rotations, if the previous rotation and the current rotation are the same, you can't know whether the object is not moving or is rotating at one turn per frame.

Then its better to use Euler Integration For Orientations or there is something else? If there is, which?

Quote:

Other than that, your formulas seem quite wrong. Your proposed formula for verlet integration wouldn't work for positions either.

dt = 0.001;
temp = position;
position += (position - prev_position) + acceleration * dt * dt
prev_position = temp;

Is that wrong?

Quote:

Quaternions are a good representation for orientations, but not for angular velocities and angular accelerations. Those are properly described by pseudo-vectors.

I used angular velocities and accelerations as vectors here, but the orientations as Quaternion.

Thanks,
Kasya

P.S i put = instead of += in my last post.

[Gasim]

I just wrote code like that:

	void RigidBody::Integrate(real time) {

		real dt = 0.01f;

		Vector3 temp = position;
		Quat accel, qTemp = orientation;
		Math::QuatFromEuler(accel, angaccel);
		position += (position - prev_position) + acceleration * dt * dt;
		prev_position = temp;

		orientation += (orientation - prev_orientation) + accel * dt * dt;
		prev_orientation = qTemp;

	}

position, prev_position: Vector3
orientation, prev_orientation: Quat
acceleration, angaccel: Vector3
accel: Quat

is that right?

Thanks,
Kasya

[alvaro]

You sort of cheated me when you edited the code to make it correct. :)

I've thought about it a little, and I think the basic problem with your code is that you should be multiplying and dividing instead of adding and subtracting, since multiplication is the operation that represents composition in quaternions.

I am trying to think of the correct mathematical abstraction that would cover both the case of points in R^3 and the case of orientations, and I am not exactly sure (Differentiable manifolds with a notion of parallel transport? Lie groups?). My Geometry is too rusty.

In any case, the important thing is that at the very least you need the inertial part to work well: If no acceleration is applied, an object should rotate forever. In order to achieve that, you need to do

orientation *= orientation/old_orientation;

Then think of how to apply an angular acceleration to the result. It should be something like, after multiplying by dt*dt, expressing the rotation as a quaternion, and then multiplying that in.

You may want to renormalize the result after each update, to make sure that your quaternion's length remains 1 (otherwise it won't represent a rotation).

[alvaro]

If instead of using verlet integration you are OK with using Euler integration, I came up with this code, and it seems to work:

#include <iostream>
#include <cmath>

struct Vector3D {
  double x, y, z;
  
  Vector3D(double x, double y, double z)
    : x(x), y(y), z(z) {
  }
  
  Vector3D &operator+=(Vector3D const &v) {
    x+=v.x;
    y+=v.y;
    z+=v.z;
    
    return *this;
  }
};

Vector3D operator*(double scale, Vector3D const &v) {
  return Vector3D(scale*v.x, scale*v.y, scale*v.z);
}

double length(Vector3D const &v) {
  return std::sqrt(v.x*v.x+v.y*v.y+v.z*v.z);
}

std::ostream &operator<<(std::ostream &os, Vector3D const &v) {
  return os << '('
	    << v.x << ','
	    << v.y << ','
	    << v.z << ')';
}

struct Quaternion {
  double r;
  double x,y,z;
  
  Quaternion(double r, double x, double y, double z)
    : r®, x(x), y(y), z(z) {
  }
  
  // Rotation of alpha radians around unit vector u
  Quaternion(double alpha, Vector3D const &u) {
    double half=alpha*.5;
    double s=std::sin(half);
    double c=std::cos(half);
    
    r=c;
    x=s*u.x;
    y=s*u.y;
    z=s*u.z;
  }
  
  // Convert pseudovector to quaternion
  explicit Quaternion(Vector3D const &rot) {
    double l=length(rot);
    if(l==0.0) {
      r=1.0;
      x=y=z=0.0;
      return;
    }
    double half=l*.5;
    double s_over_l=std::sin(half)/l;
    
    r=std::cos(half);
    x=s_over_l*rot.x;
    y=s_over_l*rot.y;
    z=s_over_l*rot.z;
  }
  
  Quaternion &operator*=(Quaternion const &R) {
    Quaternion L=*this;
    
    r = L.r*R.r - L.x*R.x - L.y*R.y - L.z*R.z;
    x = L.r*R.x + L.x*R.r + L.y*R.z - L.z*R.y;
    y = L.r*R.y + L.y*R.r + L.z*R.x - L.x*R.z;
    z = L.r*R.z + L.z*R.r + L.x*R.y - L.y*R.x;
    
    return *this;
  }
  
  Vector3D apply(Vector3D const &v) const {
    Quaternion acc = *this;
    acc *= Quaternion(0, v.x, v.y, v.z);
    acc *= Quaternion(r, -x, -y, -z);
    return Vector3D(acc.x, acc.y, acc.z);
  }
};

std::ostream &operator<<(std::ostream &os, Quaternion const &q) {
  return os << '('
	    << q.r << "+"
	    << q.x << "*i+"
	    << q.y << "*j+"
	    << q.z << "*k)";
}

struct Orientation_with_angular_velocity {
  Quaternion orientation;
  Vector3D angular_velocity;
  
  Orientation_with_angular_velocity(Quaternion orientation,
				    Vector3D angular_velocity)
    : orientation(orientation), angular_velocity(angular_velocity) {
  }
  
  void apply_timestep(Vector3D const &angular_acceleration, double dt) {
    angular_velocity += dt*angular_acceleration;
    Quaternion angular_velocity_rotation(dt*angular_velocity);
    orientation *= angular_velocity_rotation;
  }
};

int main() {
  Orientation_with_angular_velocity o(Quaternion(1.0,0.0,0.0,0.0),
				      Vector3D(0.0,0.0,0.0));
  for(int i=0;i<20;++i) {
    std::cout << o.orientation << '\n';
    o.apply_timestep(Vector3D(1.0,0.0,0.0),0.1);
  }
}

[Gasim]

Thank you very very much! :) I'll change my Verlet Orientation Integration into Euler and will test it.

I have some another questions?

1. The Equivalent of Inertia Tensor in position integration is mass. What is equivalent of Moment of Inertia in position integration?

2. How can i calculate the Inertia Tensor of an object if it's geometry figure is not one of These

Inertia Tensor is Matrix in mine


P.S Please don't give me resource about reading it. I read the whole Google and found nothing interesting. Then opened a source code and seen a comment about Inertia Tensor that it is like mass but not real mass.

[alvaro]

Quote:

Original post by Kasya Thank you very very much! :) I'll change my Verlet Orientation Integration into Euler and will test it. I have some another questions?

I'll try to help you, but I am not an expert in mechanics so I might say something wrong.

Quote:

1. The Equivalent of Inertia Tensor in position integration is mass. What is equivalent of Moment of Inertia in position integration?

My understanding is that the moment of inertia is the proper equivalent of mass in rotation computations. The moment of inertia depends on which axis you are rotating around. The tensor of inertia is a compact object from which you can compute the moment of inertia with respect to any axis. You probably want to use the tensor form, although I would have to look up the details myself, since I haven't done this in the last 14 years.

Quote:

2. How can i calculate the Inertia Tensor of an object if it's geometry figure is not one of These

This link contains the definition of the moment of inertia for a rigid object of N point masses. If your object is continuous, you can try using the triple-integral expression at the end of that section, or you can view your continuous model as a limit of discrete models.

[strimpa]

Hi,

This thread is from quite a while ago, but I thought I share my insights on here.
I struggled to get this to work, too, but figured it out eventually.
The main factors you have to consider when implementing integration on Quaternions are
- the linear integration issues on quats: Search for SLERP.
- non-unit Quats being added together creating axis problems.

Regards