GLSLで画面描画に使う行列でZ軸の値がおかしくなる原因が知りたい

Question

シェーダーで行う計算をCPUで普通に計算するテストコードを実装したいのですが予測通り最後で3次元目がよくわからない値に変わっているため画面に表示されない等の不具合が発生しているということを突き止めたのですがこれはなぜこんな値になっているのでしょうか？

追記　ロドリゲス回転行列がおかしいとの修正依頼を多数受けまして自分なりに修正しましたがもうどうしてもわかりません何をしたのか教えてくれますでしょうか？printf();でデバッグするとどうやらrotateの部分で値がおかしいためすべて狂ってるみたいです。

追記　
参考サイト: http://marupeke296.sakura.ne.jp/DXG_No58_RotQuaternionTrans.html
② クォータニオン→回転行列変換より回転行列を使いかました。

下記サイトの163ページの透視射形変換行列の数式を移しました。
参考サイト: https://tokoik.github.io/GLFWdraft.pdf

#version 400 
//頂点シェーダー

in vec3  position;//頂点座標

uniform mat4 scale;//スケール行列
uniform mat4 rotate;//回転行列
uniform mat4 move;//平行移動行列

uniform mat4 MP;

out vec4 mt;

void main()
{
	vec4 t = vec4(position,1.0);
	mat4 M = mat4(scale * rotate * move);




	mt = vec4(MP * M * t);




	gl_Position = MP * M * t;
}

#include "stdio.h"
#include <iostream>

//透視射形行列
float mp[16];

//座標構造体
struct position
{
public:
    float x;
    float y;
    float z;

    float w;
};


//頂点バッファー
float Vertex[6][4] =
{
    {-0.5, 0.5, 1.0,1.0},
    {-0.5, -0.5, 1.0,1.0},
    {0.5, -0.5, 1.0,1.0},

    {-0.5, 0.5, 1.0,1.0},
    {0.5, 0.5, 1.0,1.0 },
    {0.5, -0.5, 1.0,1.0}
};



//透視投影変換行列を作る
void create_matri_mp(float top, float bottom, float left, float right,
    float near, float far, float result[16])
{
    result[0] = (2 * near) / (right - left);
    result[5] = (2 * near) / (top - bottom);

    result[8] = (right + left) / (right - left);

    result[9] = (top + bottom) / (top - bottom);

    result[10] = (far + near) / (far - near);

    result[11] = -1;
    result[14] = (2 * far * near) / (far - near);


}

//行列を表示
void print_matrix(const char* str, float mp[4][4])
{
    printf("name: %s\n", str);
    for (int i = 0; i < 4; i++)
    {
        for (int j = 0; j < 4; j++)
        {
            printf("[  %.2f   ] , ", mp[i][j]);
        }
        printf("\n");
    }
}


//行列掛け算
//void mul(float a[4][4], float b[4][4])
void mul(float a[16], float b[16], float result[16])
{
    result[0] = (a[0] * b[0]) + (a[4] * b[1]) + (a[8] * b[2]) + (a[12] * b[3]);

    result[1] = (a[1] * b[0]) + (a[5] * b[1]) + (a[9] * b[2]) + (a[13] * b[3]);

    result[2] = (a[2] * b[0]) + (a[6] * b[1]) + (a[10] * b[2]) + (a[14] * b[3]);

    result[3] = (a[3] * b[0]) + (a[7] * b[1]) + (a[11] * b[2]) + (a[15] * b[3]);

    /////////

    result[4] = (a[0] * b[4]) + (a[4] * b[5]) + (a[8] * b[6]) + (a[12] * b[7]);

    result[5] = (a[1] * b[4]) + (a[5] * b[5]) + (a[9] * b[6]) + (a[13] * b[7]);

    result[6] = (a[2] * b[4]) + (a[6] * b[5]) + (a[10] * b[6]) + (a[14] * b[7]);

    result[7] = (a[3] * b[4]) + (a[7] * b[5]) + (a[11] * b[6]) + (a[15] * b[7]);

    /////////
    result[8] = (a[0] * b[8]) + (a[4] * b[9]) + (a[8] * b[10]) + (a[12] * b[11]);

    result[9] = (a[1] * b[8]) + (a[5] * b[9]) + (a[9] * b[10]) + (a[13] * b[11]);

    result[10] = (a[2] * b[8]) + (a[6] * b[9]) + (a[10] * b[10]) + (a[14] * b[11]);

    result[11] = (a[3] * b[8]) + (a[7] * b[9]) + (a[11] * b[10]) + (a[15] * b[11]);

    /////////

    result[12] = (a[0] * b[12]) + (a[4] * b[13]) + (a[8] * b[14]) + (a[12] * b[15]);

    result[13] = (a[1] * b[12]) + (a[5] * b[13]) + (a[9] * b[14]) + (a[13] * b[15]);

    result[14] = (a[2] * b[12]) + (a[6] * b[13]) + (a[10] * b[14]) + (a[14] * b[15]);

    result[15] = (a[3] * b[12]) + (a[7] * b[13]) + (a[11] * b[14]) + (a[15] * b[15]);


}

//行列計算 a x b
void mul_A_B(float a[16], float b[16], float result[16])
{
    result[0] = (a[0] * b[0]) + (a[1] * b[4]) + (a[2] * b[8] ) + (a[3] * b[12]);
    result[1] = (a[0] * b[1]) + (a[1] * b[5]) + (a[2] * b[9] ) + (a[3] * b[13]);
    result[2] = (a[0] * b[2]) + (a[1] * b[6]) + (a[2] * b[10]) + (a[3] * b[14]);
    result[3] = (a[0] * b[3]) + (a[1] * b[7]) + (a[2] * b[11]) + (a[3] * b[15]);

    result[4] = (a[4] * b[0]) + (a[5] * b[4]) + (a[6] * b[8]) + (a[7] * b[12]);
    result[5] = (a[4] * b[1]) + (a[5] * b[5]) + (a[6] * b[9]) + (a[7] * b[13]);
    result[6] = (a[4] * b[2]) + (a[5] * b[6]) + (a[6] * b[10]) + (a[7] * b[14]);
    result[7] = (a[4] * b[3]) + (a[5] * b[7]) + (a[6] * b[11]) + (a[7] * b[15]);


    result[8]  = (a[8] * b[0]) + (a[9] * b[4]) + (a[10] * b[8]) + (a[11] * b[12]);
    result[9]  = (a[8] * b[1]) + (a[9] * b[5]) + (a[10] * b[9]) + (a[11] * b[13]);
    result[10] = (a[8] * b[2]) + (a[9] * b[6]) + (a[10] * b[10]) + (a[11] * b[14]);
    result[11] = (a[8] * b[3]) + (a[9] * b[7]) + (a[10] * b[11]) + (a[11] * b[15]);

    result[12] = (a[12] * b[0]) + (a[13] * b[4]) + (a[14] * b[8]) + (a[15] * b[12]);
    result[13] = (a[12] * b[1]) + (a[13] * b[5]) + (a[14] * b[9]) + (a[15] * b[13]);
    result[14] = (a[12] * b[2]) + (a[13] * b[6]) + (a[14] * b[10]) + (a[15] * b[14]);
    result[15] = (a[12] * b[3]) + (a[13] * b[7]) + (a[14] * b[11]) + (a[15] * b[15]);
}



//ベクトルと行列の掛け算
void mul_vec3_matrix(float m[16], float v[3], float r[3])
{
    printf("\n\n");
    for (int j = 0; j < 3; j++)
    {
        printf("[ %.2f ] ", r[j]);
    }
    printf("\n");
}


//画面表示
void print_(float tmp[16])
{
    for (int i = 0; i < 16; i++)
    {
        if (i % 4 == 0)
        {
            printf("\n");
        }

        printf(" [ %.2f ] ", tmp[i]);
    }

    printf("\n\n");

}


int main()
{
    create_matri_mp(1.0f, -1.0f, -1.0f, 1.0f, 1.0, 10.f, mp);//透視射形変換行列mp


    //回転行列
    struct position pos;

    pos.x = 1.0f;
    pos.y = 0;
    pos.z = 0;

    pos.w = 1.0f;
    float r = 1.0f;

    float rotate[16] = {
        (pos.x * pos.x * (1 - cos(r)) + cos(r)),
        (pos.x * pos.y * (1 - cos(r)) - (pos.z * sin(r))),
        (pos.x * pos.z * (1 - cos(r)) + (pos.y * sin(r))),
        0,

        (pos.x * pos.y * (1 - cos(r)) + (pos.z * sin(r))),
        (pos.y * pos.y * (1 - cos(r)) + cos(r)),
        (pos.y * pos.z * (1 - cos(r)) - (pos.x * sin(r))),
        0,

        (pos.x * pos.z * (1 - cos(r)) - (pos.y * sin(r))),
        (pos.y * pos.z * (1 - cos(r)) + (pos.x * sin(r))),
        (pos.z * pos.z * (1 - cos(r)) + cos(r)),
        0,

        0,0,0,1
    };

  //  print_(rotate);

    //平行移動
    struct position p;
    p.x = 0;
    p.y = 0;
    p.z = 0;
    p.w = 1;

    float move[16] =
    {
        1,0,0,0,
        0,1,0,0,
        0,0,1,0,
        p.x,p.y,p.z,1,
    };

 //   print_(move);

    float t[16];
    mul(rotate, move, t);
 //   print_(t);

    //スケール行列
    float scale[16] = {
        1,0,0,0,
        0,1,0,0,
        0,0,1,0,
        0,0,0,1
    };

    print_(scale);

    float t2[16];
    mul(scale, t, t2);
    print_(t);





    for (int i = 0; i < 6; i++)
    {
        float t3[16];
        mul(Vertex[i], t2, t3);
//         print_(t2);


        float t4[16];
        mul(mp, t2, t4);
 //       print_(t4);
    }


    float a[16] =
    {
        0,3,2,5,
        0,2,2,1,

        0,7,4,3,
        0,3,3,7,

    };


    float b[16] =
    {
        0,3,5,5,
        0,5,8,6,

        8,7,3,3,
        5,8,3,8,

    };


    float c[16] = { 0 };

    //mul_A_B(a,b,c);

//    print_(c);



    return 0;
}

Comments

[fana]

デバッグのための行いなのですから，そこは頑張りましょうよ．

・行列の乗算処理が間違っていないか？
→　これが疑わしいのならば，乗算処理を単体でテストして，まともなことを確認しましょう．真っ先に．
「適当に要素値を決めた行列２個と，手計算なりで求めたその積」をテストパターンとして何個か用意すれば確認できるでしょう．

・create_matri_mp()がおかしいのか？
　→　これで計算された行列には「（これこれこういう作用だ，という）意味」があるハズ．
つまり，「ある座標に乗じたならばその結果はこうなるハズだ」というINとOUTのセットを，その意味から，前述の行列の積と同じようにテストパターンとして用意して確認することができるハズ．

[退会済みユーザー]

なるほど。しかしなにか間違えているというのは事実でしょうか？

[fana]

それを調べるのですよね？

やっている計算は
最終結果 = MP * scale * rotate * move * v
なわけですが，
例えば，vに何か適当な値を定めたときに，
move * v
は期待通りの値になっているか？　を確認する．
ここが大丈夫そうならば，次は，
rotate * (move * v )
までが期待通りの値になっているか？　を同様に確認する．
ここまでが大丈夫ならば……

として，見ていけばよいのでは．
どこかで「期待と異なる結果」に遭遇したならば，怪しいのはそこです．

各段階での「期待される結果」をあなたが把握していることが必要です．（それが無いと確認のしようがないので）

[1T2R3M4]

SOJに似たような質問があったから定期的にチェックしてみては。
https://ja.stackoverflow.com/questions/70444

[退会済みユーザー]

わかりました。質問ですがこの手のバグはどうやったらいいのでしょうか？

[stdio]

create_matri_mpの result[3][3]の値が0になっているのか気になります。記憶は曖昧ですが、1.0だったような気がします。

[退会済みユーザー]

result[3][3] = 1.0;としましたがやはり値が同じくマイナスですｗ

[退会済みユーザー]

質問内容を編集しましたので再度御覧ください

[1T2R3M4]

デバッグについてはfanaさんが書いているように
入力に対して各関数の出力が期待値通りか確認して
違ったら関数内で引数は正しくわたっているか、
計算結果は正しいかとかを確認していけばいいのではないでしょうか。

[退会済みユーザー]

質問内容を大幅に変更しましたので御覧ください。

[Soei]

検証などはしていないので確かなことは言えませんが、
回転行列がおかしくないですか?
12個しか初期化要素がないため(4,4)成分にあたるrot[15]が0になってしまっています。
また、ロドリゲスの回転公式を表現したいのでしょうが、そうであれば回転軸posは単位ベクトル(大きさが1のベクトル)でなければいけません。
※追記
ロドリゲスの回転公式としてもおかしいですね。
(1 - cos(r))は足すのではなく掛けないといけませんし、ほかにもいくつか間違いがあるみたいです。

[stdio]

それともう一つ。
rotateの初期化部分で、16個初期化しないといけないにも関わらず、最初の12個しか初期化されていないのは何故でしょうか？

[退会済みユーザー]

提示コードを編集しましたので御覧ください

[退会済みユーザー]

以前にも貴方に言いましたが、ここは他人にデバッグをしてもらうサイトではないですよ？

[退会済みユーザー]

そう見えてしまうのですが。デバッグ依頼をするつもりは一切ないのですが自分の知らない範囲で何かやっているものがあるのかと思い質問しました。削除したほうがよろしいのでしょうか？

[退会済みユーザー]

別に後ろめたいことが無いなら消さなくてもいいと思いますが、Soeiさん、stdioさんにコメントで指摘されてる事くらいまともに返信したらどうでしょう。

[退会済みユーザー]

わかりました。　Soeiさん、stdioさん　修正したのですが治りませんでした。また追記の文章を御覧ください。実装を変えましたがだめでした。 rotateの最後の４つの数値に0,0,0,1というコードを忘れていたので追加しました。が結果が変わりません

[Soei]

指摘した部分が直っていませんよ。
・回転軸posは単位ベクトル(大きさが1のベクトル)でなければいけません。
・(1 - cos(r))は足すのではなく掛けないといけません。
・ほかにもいくつか間違いがあります。細かく指摘しないので自分で探してください。

[退会済みユーザー]

様々な場所のコードを編集しました。演算子のミスなどの間違えを修正しましたので御覧ください。

[退会済みユーザー]

何故私が他人にデバッグをして貰ってるように見えるかというと、コンピュータは外部から入力されるランダムなデータなど不確定要素が入らない限りは数式通りの答えしか出しません。
デバッグでステップ実行したり、変数値を出力しながら少しづつ入力値と計算結果、数式を照合しながら調べていけば、通常は自分で間違いに気づく筈なのです。

なので、きちんとデバッグをしている人が、
> 3次元目がよくわからない値に変わっているため
などという事を言うのがおかしいのです。これは、つまり自分で計算結果を何も検証していないということです。
fana さんが指摘されている、
> 各段階での「期待される結果」をあなたが把握していることが必要です．（それが無いと確認のしようがないので）
これが全てだと思います。

Answer 1

**自己解決しました。
**

まずは構造体のstruct poisition pos;を初期化していないことで変な値が入っているということも原因のうちでした。提示コードは修正してありますが

シェーダーの掛け算の順番もおかしいためそれを修正

mat4 M = mat4(move * rotate * scale);

各行列は数学上では左上から右上に進みますがOpenGLの場合は左上から左下に進みむのでまずその数式とプログラム数式の違いを修正しました。

また編集と修正にも多数コメントが有る通り値の設定の忘れていつ要素への数値の設定と演算子や符号のつけ忘れや間違いなどを修正。rotate行列

自分はなぜ行うのか不明ですが最後にある

        float r2[4];
        mul_vec4_matrix(t4, Vertex[i], r2);

コードはツイッターで回答をもらいましたて「　透視投影行列を掛けた結果のvec4はwが1ではない値になると思います。その場合はxyzをwで割った値が求めたい頂点座標　」というのがあるみたいでその数式になりますがシェーダーで本番する場合はこれは必要ないですがCPU計算による数式を作る場合は必要みたいです。