C++言語で衝突判定のプログラムを作りたいのですが…

Question

初めまして．

現在，C++言語とopenGLを用いて『複数の円同士を衝突させ，水玉模様を得る』というプログラムを作成しています．

衝突の方法に関して，自然な水玉模様を得るために，お盆の上に円を配置し，揺らす　という方法を取っています．

プログラムの流れは，

①3つの異なる半径とゆらぎ係数を入力し，それぞれの円の個数を得る

②①で求めた個数を規則正しくお盆の上に配置させる（このとき，お盆のサイズは規則正しく配置させた円が収まる適当なサイズです）

③お盆を揺らす＋お盆のサイズを縮小させる

④円同士，円とお盆の壁との衝突

⑤運動量保存，エネルギー保存

⑥目標のお盆のサイズになれば終了（目標サイズは最初に入力します）

と．このようになっています．

研究の一環としてこの作業を行っており，このプログラム自体は過去に作成した方から頂いたものです．

②までは確認できているのですが，それ以降が入力するパラメータによって結果がうまく表示されるときとされないときがあり，その原因がわからず困っています．

わかりづらくてすみません．どなたかアドバイス頂けると幸いです．

よろしくお願い致します．

環境は，Windows XP, Microsoft Visual.Net2003, openGLです．

 

Answer 1

OpenGLは関係ないかもしれませんね。
どのような状態なのか不明なので問題点を確実に指摘はできませんが、
両者とも真円(楕円でない)の場合なら、
円C1の半径をr1、中心座標がP1、円C2の半径をr2、中心座標をP2とした場合
円C1と円C2が接している、又は重なっている状態とは、
  (P1、P2間の距離) <= (r1+r2)
です。簡単ですね。
ちなみに(P1、P2間の距離)は
double distP1P2 =
    sqrt(( double)(P2.x - P1.x)*(P2.x - P1.x)
　　　　+ ( double)(P2.y - P1.y)*(P2.y - P1.y));
ですが、比較元(r1+r2)を2乗しておけば高コストなsqrt()を使わずに済むことは
おわかりですね。似たような計算をしている所を疑ってみてはどうでしょう。
最初は2個の円でやってみるのが吉かも。

Answer 2

>研究の一環としてこの作業を行っており，このプログラム自体は過去に作成した方から頂いたものです．

>②までは確認できているのですが，それ以降が入力するパラメータによって結果がうまく表示されるときとされないときがあり，その原因がわからず困っています．

Algorithmは明確になっていますか？

明確ならば、まずは計算処理を検証しましょう。

期待される結果に対して、そうでない結果が得られたのならば、計算処理を具体的に確認しましょう。

計算単位ごとに、Logに計算途中の結果を書き出したり、Debuggerで値を確認すると良いでしょう。

計算処理が正しいことが確認できたのなら、次にOpenGLの処理を確認します。

OpenGLの仕様書がありますので、それを参考にされると良いでしょう。

そうすれば、もっと具体的に問題点や不明点を洗い出せると思います。

それを質問されたほうが、良いと思いますよ。

Answer 3

仲澤＠失業者さん

回答ありがとうございます．

OpenGLは，最終的な円の状態を描画するのに使っています．説明不足でした．すみません＞＜

衝突判定の式は理解しました．

この判定式に当てはまっている場合の式が合っているのかわかりません．

以下は，衝突判定に対する式のプログラムです．もしお時間に余裕がありましたらアドバイスお願い致します．

 

#include <iostream>
#include<stdio.h>
#include <glut.h>
#include <math.h>

using namespace std;

int num1=3;//半径の種類の個数
int n_all;//石畳の総数
double yoko;//お盆初期サイズ
double tate;
double YOKO;//お盆目標サイズ
double TATE;
double sokudo_heikoukabe;//お盆の長辺が縮小する速度
double sokudo_suichokukabe;//お盆の短辺が縮小する速度
const int time1=150000;
int watch_time1;
int n[400];//半径r[i]の円の個数
double r[100];//半径(値のみ)
double rr[600];//i番目の半径（それぞれの円の半径）
double S[100];
double SS[100];
static double X;//充填率
static double Sall;//円が占める総面積
double SSall=0.;
static double A;//ゆらぎ係数（-1なら1/fゆらぎを表す）
static double B;//
double Rx[time1+1];//リサージュ曲線
double Ry[time1+1];
double xx[600][time1+1];//円の座標
double yy[600][time1+1];
double a[600];
double b[600];
double dis[600][600];
double heikoukabe1_dis[600];
double heikoukabe2_dis[600];
double suichokukabe1_dis[600];
double suichokukabe2_dis[600];
double PI=3.1415926;

double center_dis(double a1, double a2, double b1, double b2)
//(a1,b1)(a2,b2)中心間の距離を求める
{    double dis;
     dis = sqrt((a2-a1)*(a2-a1)+(b2-b1)*(b2-b1));
     return dis;
}

double x_shoutotsu(double a1_t, double a2_t, double a1_t1, double b1_t, double b2_t, double b1_t1)
//石畳同士で衝突したときのベクトルのｘ成分を求める
//(自分のｘ成分,相手のｘ成分,自分のｘ成分の一段階前,自分のｙ成分,相手のｙ成分,自分のｙ成分の一段階前)
{    double a;
     a=(a1_t-a2_t)/sqrt((a1_t-a2_t)*(a1_t-a2_t)+(b1_t-b2_t)*(b1_t-b2_t))
                  +(a1_t-a1_t1)/sqrt((a1_t-a1_t1)*(a1_t-a1_t1)+(b1_t-b1_t1)*(b1_t-b1_t1));
     return a;
}

double y_shoutotsu(double a1_t, double a2_t, double a1_t1, double b1_t, double b2_t, double b1_t1)
//石畳同士で衝突したときのベクトルのｙ成分を求める
//(自分のｘ成分,相手のｘ成分,自分のｘ成分の一段階前,自分のｙ成分,相手のｙ成分,自分のｙ成分の一段階前)
{    double b;
     b=(b1_t-b2_t)/sqrt((a1_t-a2_t)*(a1_t-a2_t)+(b1_t-b2_t)*(b1_t-b2_t))
                  +(b1_t-b1_t1)/sqrt((a1_t-a1_t1)*(a1_t-a1_t1)+(b1_t-b1_t1)*(b1_t-b1_t1));
     return b;
}

double center_heikoukabe_dis(double a, double b, double c)
//石畳の中心(a,b)と平行壁(y=c)との距離を求める
{    double dis;
     dis = fabs(c-b);
     return dis;
}

double center_suichokukabe_dis(double a, double b, double c)
//石畳の中心(a,b)と垂直壁(x=c)との距離を求める
{    double dis;
     dis = fabs(c-a);
     return dis;
}

double heikoukabe_x_shoutotsu(double a_t, double a_t1, double b_t, double b_t1)
//平行壁に衝突したときの単位ベクトルのｘ成分を求める(ｘ成分,ｘ成分の一段階前,ｙ成分,ｙ成分の一段階前)
{    double a;
     a=(a_t-a_t1)/sqrt((a_t-a_t1)*(a_t-a_t1)+(b_t1-b_t)*(b_t1-b_t));
     return a;
}

double heikoukabe_y_shoutotsu(double a_t, double a_t1, double b_t, double b_t1)
//平行壁に衝突したときの単位ベクトルのｙ成分を求める(ｘ成分,ｘ成分の一段階前,ｙ成分,ｙ成分の一段階前)
{    double b;
     b=(b_t1-b_t)/sqrt((a_t-a_t1)*(a_t-a_t1)+(b_t1-b_t)*(b_t1-b_t));
     return b;
}

double suichokukabe_x_shoutotsu(double a_t, double a_t1, double b_t, double b_t1)
//垂直壁に衝突したときの単位ベクトルのｘ成分を求める(ｘ成分,ｘ成分の一段階前,ｙ成分,ｙ成分の一段階前)
{    double a;
     a=(a_t1-a_t)/sqrt((a_t1-a_t)*(a_t1-a_t)+(b_t-b_t1)*(b_t-b_t1));
     return a;
}

double suichokukabe_y_shoutotsu(double a_t, double a_t1, double b_t, double b_t1)
//垂直壁に衝突したときの単位ベクトルのｙ成分を求める(ｘ成分,ｘ成分の一段階前,ｙ成分,ｙ成分の一段階前)
{    double b;
     b=(b_t-b_t1)/sqrt((a_t1-a_t)*(a_t1-a_t)+(b_t-b_t1)*(b_t-b_t1));
     return b;
}

double hozonsoku(double r1, double r2, double a1_t, double a1_t1, double a2_t, double a2_t1,
              double b1_t, double b1_t1, double b2_t, double b2_t1)
//石畳同士で衝突したときの運動量保存・エネルギー保存を考慮した速度ベクトルの大きさ
//（自分の半径,相手の半径,自分のｘ成分,自分のｘ成分の一段階前,相手のｘ成分,相手のｘ成分の一段階前,
//自分のｙ成分,自分のｙ成分の一段階前,相手のｙ成分,相手のｙ成分の一段階前)
{    double v;
     v=0.5*((r1*r1*PI-r2*r2*PI)*sqrt((a1_t-a1_t1)*(a1_t-a1_t1)+(b1_t-b1_t1)*(b1_t-b1_t1))/(r1*r1*PI+r2*r2*PI)
         +2*r2*r2*PI*sqrt((a2_t-a2_t1)*(a2_t-a2_t1)+(b2_t-b2_t1)*(b2_t-b2_t1))/(r1*r1*PI+r2*r2*PI));
     return v;
}

void move()//円を動かす
{    
    for( int i=0 ; i<n_all ; i+=6 ){//初期配置（規則正しく配置）
        xx[i][0]=2.3*rr[0]*(-(n_all/12)+i/6);
        yy[i][0]=-5*rr[0];
        xx[i+1][0]=2.3*rr[0]*(-(n_all/12)+i/6);
        yy[i+1][0]=-2.5*rr[0];
        xx[i+2][0]=2.3*rr[0]*(-(n_all/12)+i/6);
        yy[i+2][0]=0;
        xx[i+3][0]=2.3*rr[0]*(-(n_all/12)+i/6);
        yy[i+3][0]=2.5*rr[0];
        xx[i+4][0]=2.3*rr[0]*(-(n_all/12)+i/6);
        yy[i+4][0]=5*rr[0];
        xx[i+5][0]=2.3*rr[0]*(-(n_all/12)+i/6);
        yy[i+5][0]=7.5*rr[0];
    }

    //お盆の初期サイズ（規則正しく配置した円が収まるような大きさ）
    if(YOKO < 2*rr[0]*(n_all/10+10))yoko=-2*(2.3*rr[0]*(-(n_all/12))-4*rr[0]);
    else yoko=YOKO;
    if(TATE < 2*10*rr[0]) tate=25*rr[0];
    else tate=TATE;
    
    watch_time1=(yoko-YOKO)*time1/(yoko-YOKO+30);

    if(time1!=0){
        sokudo_heikoukabe=(tate/2-TATE/2)/watch_time1;
        sokudo_suichokukabe=(yoko/2-YOKO/2+15)/time1;
    }
    
    int T,check;

    for(int i=0 ; i<n_all-1 ; i+=3){//初期動作
        a[0+i]=-0.2/sqrt(5.0);
        b[0+i]=0.1/sqrt(5.0);

        a[1+i]=0.1/sqrt(2.0);
        b[1+i]=0.1/sqrt(2.0);

        a[2+i]=-0.1/sqrt(2.0);
        b[2+i]=-0.1/sqrt(2.0);
    }

    
    T=1;

    for(int hit=0 ; hit<50000001 ; hit++){
        for(int t=T ; t<=time1 ; t++){
            check = T;
            
            for(int i=0 ; i<n_all ; i++){
                xx[i][t]=xx[i][t-1]+a[i];
                yy[i][t]=yy[i][t-1]+b[i];
            }
    
            
            for(int i=0 ; i<n_all ; i++){
            
                        heikoukabe1_dis[i]=center_heikoukabe_dis(xx[i][t],yy[i][t],tate/2-sokudo_heikoukabe*t);
                        //平行壁1(初期y=tate/2)とi番目の円板の中心との距離
                        heikoukabe2_dis[i]=center_heikoukabe_dis(xx[i][t],yy[i][t],-tate/2+sokudo_heikoukabe*t);
                        //平行壁2(初期y=-tate/2)とi番目の円板の中心との距離
                        suichokukabe1_dis[i]=center_suichokukabe_dis(xx[i][t],yy[i][t],yoko/2-sokudo_suichokukabe*t);
                        //垂直壁1(初期x=yoko/2)とi番目の円板の中心との距離
                        suichokukabe2_dis[i]=center_suichokukabe_dis(xx[i][t],yy[i][t],-yoko/2+sokudo_suichokukabe*t);
                        //垂直壁2(初期x=-yoko/2)とi番目の円板の中心との距離
                                    
                        if(heikoukabe1_dis[i]<=rr[i]+0.001 || heikoukabe2_dis[i]<=rr[i]+0.001){//平行壁1(初期y=tate/2), 平行壁2(初期y=-tate/2)と円との衝突判定
                            a[i]=heikoukabe_x_shoutotsu(xx[i][t], xx[i][t-1], yy[i][t], yy[i][t-1])
                                *pow(0.999999999999,hit);//衝突毎に少しずつ速度を落とす

                            b[i]=heikoukabe_y_shoutotsu(xx[i][t], xx[i][t-1], yy[i][t], yy[i][t-1])
                                *pow(0.999999999999,hit);
                                
                            T=t;
                            goto END_OF_OUTER_LOOP;
                        }

                        if(suichokukabe1_dis[i]<=rr[i]+0.001 || suichokukabe2_dis[i]<=rr[i]+0.001){//垂直壁1(初期y=tate/2), 垂直壁2(初期y=-tate/2)と円との衝突判定
                            a[i]=suichokukabe_x_shoutotsu(xx[i][t], xx[i][t-1], yy[i][t], yy[i][t-1])
                                *pow(0.999999999999,hit);

                                
                            b[i]=suichokukabe_y_shoutotsu(xx[i][t], xx[i][t-1], yy[i][t], yy[i][t-1])
                                *pow(0.999999999999,hit);
                                
                            T=t;
                            goto END_OF_OUTER_LOOP;
                        }
                                        
            }
                   
            
            for(int i=1 ; i<n_all ; i++){
                for(int j=0 ; j<i ; j++){
                        dis[i][j]=center_dis(xx[i][t],xx[j][t],yy[i][t],yy[j][t]);//i番目の石畳とｊ番目の石畳の距離
                                    
                        if(dis[i][j]<=(rr[i]+rr[j]+0.001)){//円同士の衝突判定
                            a[i]=x_shoutotsu(xx[i][t],xx[j][t],xx[i][t-1],yy[i][t],yy[j][t],yy[i][t-1])
                                *hozonsoku(rr[i],rr[j],xx[i][t],xx[i][t-1],xx[j][t],xx[j][t-1],yy[i][t],yy[i][t-1],yy[j][t],yy[j][t-1])
                                *pow(0.999999999999,hit);

                            b[i]=y_shoutotsu(xx[i][t],xx[j][t],xx[i][t-1],yy[i][t],yy[j][t],yy[i][t-1])
                                *hozonsoku(rr[i],rr[j],xx[i][t],xx[i][t-1],xx[j][t],xx[j][t-1],yy[i][t],yy[i][t-1],yy[j][t],yy[j][t-1])
                                *pow(0.999999999999,hit);
                                
                            a[j]=x_shoutotsu(xx[j][t],xx[i][t],xx[j][t-1],yy[j][t],yy[i][t],yy[j][t-1])
                                *hozonsoku(rr[j],rr[i],xx[j][t],xx[j][t-1],xx[i][t],xx[i][t-1],yy[j][t],yy[j][t-1],yy[i][t],yy[i][t-1])
                                *pow(0.999999999999,hit);

                            b[j]=y_shoutotsu(xx[j][t],xx[i][t],xx[j][t-1],yy[j][t],yy[i][t],yy[j][t-1])
                                *hozonsoku(rr[j],rr[i],xx[j][t],xx[j][t-1],xx[i][t],xx[i][t-1],yy[j][t],yy[j][t-1],yy[i][t],yy[i][t-1])
                                *pow(0.999999999999,hit);
                            
                            T=t;
                            goto END_OF_OUTER_LOOP;
                        }
                }
            }

        
        }
        
        if (T==check)  break;
        END_OF_OUTER_LOOP:
        ;

           cout << T << endl;

    }
    
}

Answer 4

kozzさん

回答ありがとうございます．

アルゴリズム自体が本当に合っているのか不安がある状態です．

仲澤＠失業者さんへの返信に，衝突判定の部分のプログラムを記述させて頂いたのですが，．

衝突判定が出た後の処理はこれでよいのでしょうか？

Answer 5

これは、C++ ではなく、C ですね :p

どの様な計算規則で目的のものを求めようとしているのか、わかりません。最初の質問には「水玉模様」とか「円」と書いてありましたが、コード中には「石畳」という言葉があります。これらの関係は、どうなっているのでしょうか。

たとえば、「heikoukabe1_dis[i]」変数に、平行壁（水平壁、では？）との距離が入るようです。これが、円の半径をわずか（0.001）に大きくした値よりも小さいなら、「衝突した」と判断しているようです。この、わずかな量の根拠は、何でしょうね。

「それ以降が入力するパラメータによって結果がうまく表示されるときとされないときがあり」と言うことですが、「うまく表示されないとき」というのは、どの様な期待値に対してどの様な実測値なのでしょうか。

「複数の円同士を衝突させ，水玉模様を得る」と言うことですが、最終的に、どの様になるのでしょうか。円が衝突して分裂し、それが盆を埋め尽くすのかと思いきや、コードには盆が小さくなると書かれている。すると、最初に与えられた円がギュウギュウ詰めになるようにしたい？

盆を揺らすとは書いてありますが、盆が揺れるベクトル、円が何時なぜ動くのかについては書かれていません。コードにも表現されていないように思います。いきなり move 関数が現れています。そういうのがわからないと、アルゴリズムがあっているかどうかもわからないと思うのですが、どうでしょう？そして、アルゴリズムは、研究室で検証するものではないでしょうか？

追加：
あれれ？なんか、おかしくないですか？仲澤さんが示されたのは、「衝突したかどうかを判定する式」であって、「衝突した後にどう動くかという式」ではないですよね？たとえば、摩擦係数が0として、止まっている壁にぶつかると、…垂直壁に 横方向の移動量10mm/sec、縦方向の移動量5mm/sec でぶつかると、縦方向の移動量は変わらず、横方向は移動量が＋－反転することになるでしょう。円同士の場合は、接線にぶつかったとして、新しいベクトルを計算することになるでしょう。一度に全部動かそうとせず、まず1つの円をバウンドさせる、2つの円をぶつける、というように、段階を踏んだ方がいいと思いますよ。

---

Jitta＠わんくま同盟

Answer 6

jittaさん

回答ありがとうございます．

まず，「水玉模様」，「円」，「石畳」は全て同じものを指しています．書き方が全く統一できてなかったです．すみません．

 

「heikoukabe_dis[i]」のところで，円の半径に0.001を足しているのは，最終的なデザインを得るときに円同士がくっついてしまう状態を避けたかったからです．

円が衝突してその後位置が移動すれば，円同士がくっつくという状態にはならないですが，このプログラムでは，衝突の度に少しずつ速度を落とすようにしたかったので円同士がくっついてしまう可能性があるかと思い，0.001という値を足したものを衝突判定の式に使用しました．

 

デザインの最終的な形は，うまくいったときは以下のようになります．（ちゃんと表示できているでしょうか・・・；）

ぎゅうぎゅう詰めになる場合もありますが，円の個数によってはまばらになることもあります．（move関数の前に円の個数を決めるプログラムを組みます）

しかし，うまくいかないときは以下のように真ん中にポツンと一つだけ表示されるだけになってしまいます．

この原因がわかりません．

衝突後の位置の計算方法もイマイチよくわかりません．

 

どちらにせよ，一度にいろいろやろうと焦りすぎですよね；

jittaさんのおっしゃるように，一つ一つ段階を踏んで行こうと思います．ありがとうございました．

 

Answer 7

スクリーンショットは貼れないようですね；すみません．

Answer 8

>アルゴリズム自体が本当に合っているのか不安がある状態です．
まずはAlgorithmを確定する事から始める必要がありますね。
例えば、小数点の誤差や精度をどこまで求めるのか、なども考察する必要がありますね。 
そうしないと、何が正しくて、何が間違いか、誰も分からないですよ。

「円 衝突 アルゴリズム」で検索すると、いくつか参考情報が見つかりますが、私からは確実なことが言えません。

>スクリーンショットは貼れないようですね；すみません．
SkyDriveなどに画像をUploadして、その画像を右クリックで表示されるメニューからコピーします。

この画面（返信等の編集画面）で右クリックで表示されるメニューから貼り付けます。

うまくいかなければ、以下を参考にしてください。
  http://social.answers.microsoft.com/Forums/ja-JP/answersfeedbackja/thread/2907d2bd-d218-422c-ac40-284df488845f

テスト投稿フォーラムで試してみてください。

Answer 9

　とりあえず、「編集」して、編集ウインドウの上にある[HTML]というボタンをクリック、出てきた HTML 編集ウインドウで img 要素を探して削除しておきましょう。

　「ボール 衝突モデル」で検索して、「ボールの衝突運動シミュレーション」というページを見つけました。衝突後の運動については、「2ページ」目にかかれています。「4ページ」目に、コードがあります。円の運動は、この様にされたいのではないでしょうか。

---

Jitta＠わんくま同盟