ODE

ODE(Open Dynamics Engine)とは†

ODEは剛体ダイナミクスをシミュレーションするためのオープンソースで高性能なC/C++ライブラリです．その特徴は次のとおり

• 安定
• プラットフォーム非依存
• 先進的なジョイントタイプの取り扱い
• 摩擦まで考慮した衝突検出・応答

インストール†

• ダウンロード
  ODE Download(SourceForge)から，Windows用の ode-(バージョン番号).zip をダウンロードして，適当な解凍ソフトで解凍する．解凍してできたフォルダ(ode-*)を適当なところに置く．ここでは置いたフォルダのパスを $ODE と表現する．

• ライブラリのビルド
  
  • 0.10.*より前のバージョンでは，$ODE\build に各環境に合わせたコンパイル・ビルド用ファイルがあるので，自分の環境に合わせてビルド．Windows/Visual Studio 2005 の場合，$ODE\build\vs2005\ode.sln を用いる．用意されているプロジェクトの構成は，DebugDLL,DebugLib,ReleaseDLL,ReleaseLib の四つ．
    
  • 0.10.*以降では$ODE\build\にあるpremakeを使ってビルド用ファイルを作成する． コマンドプロンプトからpremake --helpでコマンドの説明が出るので， それを参考にここではVisual Studio 2005用のビルドファイルを作ってみる．

    premake --with-demos --with-tests --target vs2005

  • 0.11からpremake.exeとpremake4.exe，0.11.1からはpremake4.exeのみが同梱されている．

    premake4 --help

    でコマンドの説明が出る．
    1. デモも含めて Visual Studio 2005用のビルドファイルを作る場合は，

       premake4 --with-demos vs2005

       をコマンドプロンプトから打ち込む．
       
    2. vs2005フォルダが作られるので，その中のode.slnを使ってライブラリをビルドする． 用意されているプロジェクトの構成は，

       DebugDoubleDLL,DebugDoubleLib,DebugSingleDLL,DebugSingleLib,
       ReleaseDoubleDLL,ReleaseDoubleLib,ReleaseSingleDLL,ReleaseSingleLib

       の八つ． Visual Studioのツールバーからビルドしたい構成を選択するか(下図左)， メニュー「ビルド」->「構成マネージャ」からアクティブソリューション構成で選択する(下図右)． 全構成を一気にビルドしたい場合は，「ビルド」->「バッチビルド」を使うと便利である．
       
       ここでは，*Libを用いるが，コード生成のランタイムライブラリがマルチスレッドDLLになっているので， スタティックリンクにしたい場合，プロジェクトのプロパティから「コード生成」->「ランタイムライブラリ」をマルチスレッドにする(下図参照)．(複数のプロジェクトの設定を一度に変更する場合は，ソリューションエクスプローラでShiftキーを押しながらクリックで複数プロジェクトを選択後，設定する)．
       
    3. ビルドされたファイルは，$ODE\lib\debuglib\,$ODE\lib\releaselib\ などに入っている．ode_*.lib が必須ファイルで，これを自分のプロジェクトにライブラリファイルとして追加する．また，ODEが用意したGUI環境を用いるなら drawstuff.lib も必要． 追加のインクルードディレクトリは，$ODE\include ．
    4. プロジェクトの設定で，ode_*.lib drawstuff.lib を追加のライブラリファイルにする(*はsingleもしくはdouble)． drawstuff を用いないなら，ode_*.lib だけでもよい．

• 倍精度 0.10.*以降では，上記のプロジェクト構成で*Double*のものを用いると倍精度でビルドされる．
  それ以前のバージョンでは，
  ODE のライブラリは標準では単精度(float)となっているが，ビルド時に倍精度(double)に変更可能．$ODE\include\ode\config.h 内の最初のほうにある

  #define dSINGLE 1

  を

  #define dDOUBLE 1

  で置き換えてから，ライブラリをビルドしなおす．

シミュレーション空間†

ODEのシミュレーション空間はWorldやSpace,Body,Geometryなどで構成される．

• World : ODE内の動力学ワールド．剛体(Body)を登録し，動力学に従い動作させる．
• Space : 形状データ(Geometry)を登録する空間．GeometryはBodyと関連づけられ，衝突判定に用いられる．
• Joint and Constraint : 2つのBody間の相対位置，方向のことで，単に接触している2物体の他に現実のジョイント(ドアのヒンジや車軸など)も含む．

重力や何も設定しないときのBodyの空気抵抗などはWorldに対して設定する． また，Jointに関するパラメータとしてERPやCFMなどを設定する．これらは，

• ERP(error reduction parmeter) : 接続されているジョイントがユーザ操作やシミュレーションにより離れてしまった場合， 2接続点が近づくような力を負荷する．このときの力をコントロールするためのパラメータ．[0,1]．
• CFM(constraint force mixing) : 剛体間のめり込みを許す "ソフトな" ジョイントを許すためのパラメータ． 0なら完全に "ハードな" ジョイントとなり，正の0でない値なら "ソフトな" ジョイントとなる．

ODEのシミュレーション空間の構築例は以下

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35

dWorldID g_dWorld;                dSpaceID g_dSpace;                dJointGroupID g_dContactGroupID;    bool g_bWorldCreated = false;
 
void InitODE(void)
{
    dInitODE2(0);
 
        if(g_bWorldCreated){
        dJointGroupDestroy(g_dContactGroupID);
        dSpaceDestroy(g_dSpace);
        dWorldDestroy(g_dWorld);
    }
 
        g_dWorld = dWorldCreate();
    g_dSpace = dHashSpaceCreate(0);
    g_dContactGroupID = dJointGroupCreate(0);
 
        dWorldSetGravity(g_dWorld, 0, -9.8, 0);                dWorldSetCFM(g_dWorld, 1e-5);                        dWorldSetContactMaxCorrectingVel(g_dWorld, 0.1);
    dWorldSetERP(g_dWorld, 0.8);                     
    dWorldSetLinearDamping(g_dWorld, 0.00001);            dWorldSetAngularDamping(g_dWorld, 0.005);            dWorldSetMaxAngularSpeed(g_dWorld, 200);            dWorldSetContactSurfaceLayer(g_dWorld, 0.05);        dWorldSetAutoDisableFlag(g_dWorld, 1);             
    g_bWorldCreated = true;
}

ちなみに，デフォルトのパラメータ値はそれぞれ，

Gravity = (0, 0, 0)
CFM = 10^-5 (Typical values 0.1〜0.8)
ERP = 0.2 (Typical values 10^-9〜1)
LinearDamping = 0
AngularDamping = 0
ContactSurfaceLayer = 0
AutoDisableFlag = 0

また，空間を破棄するときは，

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15

void CleanODE(void)
{
    if(g_bWorldCreated){
        for(int i = 0; i < (int)g_vRObjs.size(); ++i){
            dBodyDestroy(g_vRObjs[i]);
        }
        dGeomDestroy(g_dGround);
 
        dJointGroupDestroy(g_dContactGroupID);
        dSpaceDestroy(g_dSpace);
        dWorldDestroy(g_dWorld);
 
        dCloseODE();
    }
}

剛体の追加†

BodyをODEのシミュレーション空間に登録する手順は以下．

1. dMassSet* でBodyの質量をその大きさから算出(計算はODE側で実装されているので密度と大きさを与える)
2. dCreate* で形状(Geometry)を作成
3. dBodyCreate で剛体オブジェクトを作成し，質量を設定(dBodySetMass)，形状を登録(dGeomSetBody)
4. dBodySetPosition，dBodySetRotation，dBodySetLinearVelなどで位置，姿勢，速度などを設定

立方体形状を登録するサンプルコードを以下に示す．

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20

dBodyID SetODECube(Vec3 pos, double len)
{
        dMass mass;
    dMassSetBox(&mass, 1.0, (dReal)len, (dReal)len, (dReal)len);
    dMassAdjust(&mass, 1.0);
 
        dGeomID geom = dCreateBox(g_dSpace, (dReal)len, (dReal)len, (dReal)len);
 
        dBodyID body = dBodyCreate(g_dWorld);
    dBodySetMass(body, &mass);
    dGeomSetBody(geom, body);
 
        dBodySetPosition(body, pos[0], pos[1], pos[2]);
 
    return body;
}

ODEでは剛体の特性として以下の3つを指定できる．

• 動力学(動作)剛体(Dynamics (moving) rigidbodies)
  • BodyとGeometry
  • 毎シミュレーションステップで力学に基づきその位置と姿勢が更新される
• 静的剛体(Static rigidbodies)
  • Geometryのみ
  • 衝突が起こっても動かない
• 運動学剛体(Kinematic rigidbodies)
  • BodyとGeometry，dBodySetKinematicで設定．BodyがKinematicかどうかはdBodyIsKinematicで確認できる．
  • ユーザ操作に従って動くだけで，上記二つとの衝突では影響を与えるが，影響は受けない．

タイムステップ†

ODEのシミュレーションステップを進めるには，dWorldQuickStep関数などを用いる．

  1
  2
  3
  4
  5
  6
  7
  8
  9

void StepODE(void)
{
    NearCallData data;
    data.world = &g_dWorld;
    data.contact = &g_dContactGroupID;
    dSpaceCollide(g_dSpace, &data, &ODENearCallback);
    dWorldQuickStep(g_dWorld, g_fDt);
    dJointGroupEmpty(g_dContactGroupID);
}

ODEのステップを進める関数としては，dWorldStepの方がより精度が高いが， 大きな線形システム(サイズmとする)を解くため，最悪mの3乗のオーダの計算時間がかかる． これに対して，dWorldQuickStepは反復的な手法を用いているので，m*Nのオーダで済む(Nは反復回数)． ただし，精度は低くなることに注意が必要である．反復数Nは dWorldSetQuickStepNumIterationsで設定，dWorldGetQuickStepNumIterationsで取得できる(デフォルトは20)． さらに，dWorldQuickStepでは緩和法の一種であるSORを用いており， その緩和加速係数を，dWorldSetQuickStepW/dWorldGetQuickStepWで設定/取得できる(デフォルトは1.3)．

ステップ関数では第二引数でタイムステップ幅を指定できるが， ODEはタイムステップ幅が固定されていることを想定して設計されているため， 可変タイムステップは使用しない方がよい．

また，上記コード例では，NearCallback(剛体同士が近づいた時に呼び出される関数)も設定している． NearCallback関数の例は以下．

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47

struct NearCallData
{
    dWorldID* world;
    dJointGroupID* contact;
};
 
 

static void ODENearCallback(void *data, dGeomID o1, dGeomID o2)
{
        dBodyID b1 = dGeomGetBody(o1);
    dBodyID b2 = dGeomGetBody(o2);
    if(b1 && b2 && dAreConnected(b1, b2)){
        return;
    }
    
        dContact contact[RX_MAX_CONTACTS];
    int n = dCollide(o1, o2, RX_MAX_CONTACTS, &contact[0].geom, sizeof(dContact));
 
    if(n > 0){            for(int i = 0; i < n; ++i){
                        contact[i].surface.mode = dContactBounce | dContactSlip1 | dContactSlip2/* | dContactSoftERP | dContactSoftCFM | dContactApprox1*/;
 
            contact[i].surface.mu = 1.0;                    contact[i].surface.bounce = 0.1;                contact[i].surface.bounce_vel = 0.001;
            contact[i].surface.slip1 = 0.0;
            contact[i].surface.slip2 = 0.0;
//            contact[i].surface.soft_erp = 0.1;
//            contact[i].surface.soft_cfm = 0.01;
 
            dJointID c = dJointCreateContact(*(((NearCallData*)data)->world), *(((NearCallData*)data)->contact), contact+i);
            dJointAttach(c, dGeomGetBody(o1), dGeomGetBody(o2));
        }
 
    }
}

OpenGLでの描画†

ODEは独自の描画フレームワークdrawstuffを含んでいるが， より詳細な描画コントロールを行い時はOpenGLなどを使うこともできる． 基本的には，

1. dBodyGetPosition,dBodyGetRotation でBody の位置，姿勢を取得し，glTranslatef, glMultMatrixf でOpenGL に設定
2. Geometry をdBodyGetFirstGeom 関数でBody から取得
3. dGeomGetClass 関数で形状の種類をGeometry から取得
4. それぞれの形状に合わせてパラメータ(例えばBox ならdGeomBoxGetLengths で辺の長さ) を取得し，OpenGL で描画 の手順で行う．

Box,Sphere,Capsule,Cylinder形状に対してOpenGLで描画を行うコード例は以下である．

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46

void DrawODEBody(dBodyID body)
{
    glPushMatrix();
 
    Vec3 pos;
    GLfloat m[16];
 
    pos = ODE_GetPosition(body);
    GetRotationMatrix(body, m);
    
    glTranslatef(pos[0], pos[1], pos[2]);
    glMultMatrixf(m);
 
        dGeomID geom = dBodyGetFirstGeom(body);
 
        int type = dGeomGetClass(geom);
 
        if(type == dBoxClass){
        dReal sl[4];
        dGeomBoxGetLengths(geom, sl);
 
        glScalef(sl[0], sl[1], sl[2]);
        glutSolidCube(1.0);
    }
    else if(type == dSphereClass){
        dReal rad = dGeomSphereGetRadius(geom);
 
        glutSolidSphere(rad, 32, 16);
    }
    else if(type == dCapsuleClass){
        dReal rad, len;
        dGeomCapsuleGetParams(geom, &rad, &len);
 
        DrawCapsule(rad, len, 20);
    }
    else if(type == dCylinderClass){
        dReal rad, len;
        dGeomCylinderGetParams(geom, &rad, &len);
 
        DrawCylinder(rad, len, 20);
    }
    glPopMatrix();
}

ここで，GetRotationMatrixは下で説明しているODEからOpenGLの変換行列を得る関数である． また，DrawCapsule,DrawCylinderはそれぞれ，Capsule形状，Cylinder形状をOpenGLで描画する関数である (OpenGL - サンプルコードを参照)．

ODEとOpenGLの変換行列について†

OpenGLでプリミティブを描画するときは，dGeomGetPositionで現在の位置，dBodyGetRotationで回転行列などを 取得し，glTranslatefやglMultiMatrix()を使って描画する． このとき，ODEの回転行列は4x3であり，以下のような順序で配列に格納されている．

| m[0] m[1] m[2]  m[3] | 
| m[4] m[5] m[6]  m[7] |
| m[8] m[9] m[10] m[11] |

回転行列が4x3になっているのは，SIMD friendlyにするためらしい(ODE WikiのFAQ参照)．また，配列内の並びは row major である(これもODE WikiのFAQ参照)． 一方，OpenGLのモデルビュー変換行列は4x4で，

| m[0] m[4] m[8]  m[12] | 
| m[1] m[5] m[9]  m[13] |
| m[2] m[6] m[10] m[14] |
| m[3] m[7] m[11] m[15] |

のように column major であるので，変換が必要となる． つまり，以下のような関数を用いればよい．

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19

 void GetRotationMatrix(dBodyID body, GLfloat m[16])
 {
     m[0]  = dBodyGetRotation(body)[0];
     m[1]  = dBodyGetRotation(body)[4]; 
     m[2]  = dBodyGetRotation(body)[8]; 
     m[3]  = 0.0f; 
     m[4]  = dBodyGetRotation(body)[1]; 
     m[5]  = dBodyGetRotation(body)[5]; 
     m[6]  = dBodyGetRotation(body)[9]; 
     m[7]  = 0.0f; 
     m[8]  = dBodyGetRotation(body)[2]; 
     m[9]  = dBodyGetRotation(body)[6]; 
     m[10] = dBodyGetRotation(body)[10]; 
     m[11] = 0.0f; 
     m[12] = dBodyGetRotation(body)[3]; 
     m[13] = dBodyGetRotation(body)[7]; 
     m[14] = dBodyGetRotation(body)[11]; 
     m[15] = 1.0f;
 }

プリミティブ†

球体(Sphere)†

球体ジオメトリを作成する．

dGeomID dCreateSphere(dSpaceID space, dReal radius);

球体の半径を指定する．これ以降の形状でも同様だが，これらの命令はジオメトリを作成するものなので，その位置などはジオメトリが関連付けされたボディ(dBodyID)で

void dBodySetPosition(dBodyID, dReal x, dReal y, dReal z);
void dBodySetRotation(dBodyID, const dMatrix3 R);
void dBodySetQuaternion dBodyID, const dQuaternion q);
void dBodySetLinearVel(dBodyID, dReal x, dReal y, dReal z);
void dBodySetAngularVel(dBodyID, dReal x, dReal y, dReal z);

などを用いて設定する． 球体のパラメータの設定と取得は，

void dGeomSphereSetRadius(dGeomID sphere, dReal radius);
dReal dGeomSphereGetRadius(dGeomID sphere);

で行う．また，ある点(x,y,z)が球体内部にあるかどうかなどの判定に用いることができるDepth値の取得には，

dReal dGeomSpherePointDepth(dGeomID sphere, dReal x, dReal y, dReal z);

を用いることができる．ここで，Depth値は球体内部で正，外部で負，球体表面でゼロとなる． !ボックス(Box)

dGeomID dCreateBox(dSpaceID space, dReal lx, dReal ly, dReal lz);

ボックスのx,y,z軸長さを指定

平面(Plane)†

dGeomID dCreatePlane(dSpaceID space, dReal a, dReal b, dReal c, dReal d);

平面の方程式 ax+by+cz+d=0 の各係数を指定．n = (a, b, c)が平面の法線，d は原点からの符号付距離．

円筒(Cylinder)†

両端面がフラットな円筒形状．

dGeomID dCreateCylinder(dSpaceID space, dReal radius, dReal length);

パラメータの設定と取得は

void dGeomCylinderSetParams(dGeomID cylinder, dReal radius, dReal length);
void dGeomCylinderGetParams(dGeomID cylinder, dReal *radius, dReal *length);

円筒(Capped Cylinder, Capsule)†

dGeomID dCreateCapsule(dSpaceID space, dReal radius, dReal length);

球体を一方向にのばしたような形状の円筒(カプセル形状)で半径と長さを指定．dCreateCCylinder でも可

#define dCreateCCylinder dCreateCapsule

初期状態ではジオメトリのローカル座標のz軸方向に平行．質量設定の計算は，

void dMassSetCappedCylinder(dMass *, dReal density, int direction, dReal radius, dReal length);
void dMassSetCappedCylinderTotal(dMass *, dReal total_mass, int direction, dReal radius, dReal length);

の二つで行う．上が密度指定，下が全体質量指定．direction は円筒の length 方向がどの軸に平行かを設定(x=1, y=2, z=3)． !光線(Ray)

dGeomID dCreateRay(dSpaceID space, dReal length);

他のジオメトリと異なり体積を持たない線を記述する．衝突時の接触点は1ヶ所になる．

三角形メッシュ(Triangle Mesh)†

dGeomID dCreateTriMesh(dSpaceID space, dTriMeshDataID Data,
                       dTriCallback *Callback,
                       dTriArrayCallback * ArrayCallback,
                       dTriRayCallback* RayCallback);

三角形メッシュは，メッシュデータ(dTriMeshDataID)を生成してから登録する．メッシュデータの生成は，

dTriMeshDataID dGeomTriMeshDataCreate();
void dGeomTriMeshDataBuild(dTriMeshDataID g, const void* Vertices,
                           int VertexStride, int VertexCount,
			    const void* Indices, int IndexCount,
                           int TriStride, const void* Normals);

で行う．

関数など†

位置と回転の取得†

const dReal* dBodyGetPosition(dBodyID);
const dReal* dBodyGetRotation(dBodyID);
const dReal* dBodyGetQuaternion(dBodyID);

の三つ． dBodyGetPositionは位置の取得で，返値は dReal[3] ．(x,y,z)座標値． dBodyGetRotationは回転行列(4x3)の取得で，返値は dReal[4*3] ．配置はソース中では R[(i)*4+(j)] (i=0,1,2, j=0,1,2,3)と定義されている． dBodyGetQuaternionは四元数の取得で，返値は dReal[4] ．q = (cos(q0/2), sin(q0/2)*q1, sin(q0/2)*q2, sin(q0/2)*q3) = q0+q1*i+q2*j+q3*k．

サンプル†

立方体落下(OpenGL描画)†

• box.zip
  • Visual Studio 2005用プロジェクトファイル同梱
  • ODEに登録したBodyからShapeを参照して，OpenGLで描画する．

リンク†

• Open Dynamics Engine(本家)
• Open Dynamics Engine v0.5 User Guide
• ODE Wiki