月火水木金土日
C++でI×J行列 × J×K行列
4x4行列の積を求めるプログラムをほんの少し改良しただけなのだが、C++でIxJ行列とJxK行列の積を求めるプログラムを書く。与える行列の表現は一次元配列で、OpenGLと同じ表現方法。
コード
//! @brief IxJ行列の i,j の要素にアクセスする //! @param [in] i 行番号 //! @param [in] j 列番号 //! @param [in] I 行数 inline int matijI(const int i, const int j, const int I) { return i + I * j; } //! @brief IxJ行列 × JxK行列 の積 //! @param [out] C Cik 計算結果格納先 //! @param [in] A Aijの行列 //! @param [in] B Bjkの行列 //! @param [in] I Aの行列の行数 //! @param [in] J Aの行列の列数 //! @param [in] K Bの行列の列数 //! @return C template< typename real_tC, typename real_tA, typename real_tB> real_tC* mult_matrixIJK(real_tC* C, const real_tA* A, const real_tB* B, const int I, const int J, const int K) { for (int i = 0; i < I; i++) for (int k = 0; k < K; k++) { C[matijI(i, k, I)] = 0.0; for (int j = 0; j < J; j++) { C[matijI(i, k, I)] += A[matijI(i, j, I)] * B[matijI(j, k, J)]; } } return C; }
動作検証
Eigenとの結果比較
#include <Eigen/Dense>
void mult_matrix_check() {
std::cout << "---------- Eigen ------------" << std::endl; Eigen::MatrixXf mat53(5,3); mat53(0, 0) = 00; mat53(0, 1) = 01; mat53(0, 2) = 02; mat53(1, 0) = 10; mat53(1, 1) = 11; mat53(1, 2) = 12; mat53(2, 0) = 20; mat53(2, 1) = 21; mat53(2, 2) = 22; mat53(3, 0) = 30; mat53(3, 1) = 31; mat53(3, 2) = 32; mat53(4, 0) = 40; mat53(4, 1) = 41; mat53(4, 2) = 42; Eigen::MatrixXf mat34(3, 4); mat34(0, 0) = 100; mat34(0, 1) = 101; mat34(0, 2) = 102; mat34(0, 3) = 103; mat34(1, 0) = 110; mat34(1, 1) = 111; mat34(1, 2) = 112; mat34(1, 3) = 113; mat34(2, 0) = 120; mat34(2, 1) = 121; mat34(2, 2) = 122; mat34(2, 3) = 123; std::cout << "---------- Eigen mat 1" << std::endl; std::cout << mat53 << std::endl; std::cout << "---------- Eigen mat 2" << std::endl; std::cout << mat34 << std::endl; std::cout << "---------- Eigen mat1*mat2" << std::endl; std::cout << mat53 * mat34 << std::endl;
std::cout << "---------- 自前実装 ------------" << std::endl; int I; float fmat53[5 * 3]; I = 5; fmat53[matijI(0, 0, I)] = 00; fmat53[matijI(0, 1, I)] = 01; fmat53[matijI(0, 2, I)] = 02; fmat53[matijI(1, 0, I)] = 10; fmat53[matijI(1, 1, I)] = 11; fmat53[matijI(1, 2, I)] = 12; fmat53[matijI(2, 0, I)] = 20; fmat53[matijI(2, 1, I)] = 21; fmat53[matijI(2, 2, I)] = 22; fmat53[matijI(3, 0, I)] = 30; fmat53[matijI(3, 1, I)] = 31; fmat53[matijI(3, 2, I)] = 32; fmat53[matijI(4, 0, I)] = 40; fmat53[matijI(4, 1, I)] = 41; fmat53[matijI(4, 2, I)] = 42; float fmat34[3 * 4]; I = 3; fmat34[matijI(0, 0, I)] = 100; fmat34[matijI(0, 1, I)] = 101; fmat34[matijI(0, 2, I)] = 102; fmat34[matijI(0, 3, I)] = 103; fmat34[matijI(1, 0, I)] = 110; fmat34[matijI(1, 1, I)] = 111; fmat34[matijI(1, 2, I)] = 112; fmat34[matijI(1, 3, I)] = 113; fmat34[matijI(2, 0, I)] = 120; fmat34[matijI(2, 1, I)] = 121; fmat34[matijI(2, 2, I)] = 122; fmat34[matijI(2, 3, I)] = 123; std::cout << "---------- 自前実装 mat 1" << std::endl; std::cout << mat_to_string(fmat53, "%3.0lf", 5, 3); std::cout << "---------- 自前実装 mat 2" << std::endl; std::cout << mat_to_string(fmat34, "%3.0lf", 3, 4); std::cout << "---------- 自前実装 mat1*mat2" << std::endl; float fmat54[5 * 4]; mult_matrixIJK(fmat54, fmat53, fmat34, 5, 3, 4); std::cout << mat_to_string(fmat54, "%5.0lf", 5, 4);
}
出力結果
---------- Eigen ------------ ---------- Eigen mat 1 0 1 2 10 11 12 20 21 22 30 31 32 40 41 42 ---------- Eigen mat 2 100 101 102 103 110 111 112 113 120 121 122 123 ---------- Eigen mat1*mat2 350 353 356 359 3650 3683 3716 3749 6950 7013 7076 7139 10250 10343 10436 10529 13550 13673 13796 13919
---------- 自前実装 ------------ ---------- 自前実装 mat 1 0 1 2 10 11 12 20 21 22 30 31 32 40 41 42 ---------- 自前実装 mat 2 100 101 102 103 110 111 112 113 120 121 122 123 ---------- 自前実装 mat1*mat2 350 353 356 359 3650 3683 3716 3749 6950 7013 7076 7139 10250 10343 10436 10529 13550 13673 13796 13919
OpenGLの行列の表記とC/C++での扱い
やっていると時々混乱し、致命的なミスを犯すので確認しておきたい。
まず、以下はOpenGLの移動行列。一番右側の列に移動パラメータが入る。
そしてこれは揺るがない。数学の世界、概念の世界の話なので。
次に、C/C++でこれがどう表現されているのか。以下のプログラムで確認する。
//OpenGLの機能で変換行列を作成する
glMatrixMode(GL_MODELVIEW); glTranslated(-1, 0, -2);
//作成した変換行列を取得する GLfloat modelf[16]; glGetFloatv(GL_MODELVIEW_MATRIX, modelf);
//取得した変換行列を表示する for (int i = 0; i < 16; i++) { printf("[%2d] = %3.0lf\n",i, modelf[i]); }
結果
[ 0] = 1
[ 1] = 0
[ 2] = 0
[ 3] = 0
[ 4] = 0
[ 5] = 1
[ 6] = 0
[ 7] = 0
[ 8] = 0
[ 9] = 0
[10] = 1
[11] = 0
[12] = -1
[13] = 0
[14] = -2
[15] = 1
[ 1] = 0
[ 2] = 0
[ 3] = 0
[ 4] = 0
[ 5] = 1
[ 6] = 0
[ 7] = 0
[ 8] = 0
[ 9] = 0
[10] = 1
[11] = 0
[12] = -1
[13] = 0
[14] = -2
[15] = 1
これをそのまま４要素区切りで表示すると、以下のように転置でもしたかのように見えてしまう。これはバグでも何でも無く、「現実世界の実物がプログラム中でどう表現されているか」という問題でしかない。
//OpenGLの機能で変換行列を作成する glMatrixMode(GL_MODELVIEW); glTranslated(-1, 0, -2); //作成した変換行列を取得する GLfloat modelf[16]; glGetFloatv(GL_MODELVIEW_MATRIX, modelf); //取得した変換行列を表示する for (int i = 0; i < 16; i++) { if (i % 4 == 0) puts(""); printf("%3.0lf", modelf[i]); }
結果
1 0 0 0 0 1 0 0 0 0 1 0 -1 0 -2 1
そして、OpenGLでは”行列”を、配列を用いて以下の形で表現している。
この配列に対して、数学的にアクセスしたい場合、以下のような計算で要素番号を求める。
//! @brief 4x4行列の i,j の要素にアクセスする //! @param [in] i 行番号 //! @param [in] j 列番号 inline int matij4(const int i, const int j) { return i + 4 * j; }
より一般的なアクセス
OpenGLの行列は4x4や3x3だけ考えればいいが、より一般的な行列演算関数の中でi,jから行列の添え字に変換するにはどうすればいいか。以下のような計算を行う。
//! @brief IxJ行列の i,j の要素にアクセスする //! @param [in] i 行番号 //! @param [in] j 列番号 //! @param [in] I 行数 inline int matijI(const int i, const int j, const int I) { return i + I * j; }
この考え方ではまず行を調整し、現在位置から一列の要素数×列番号だけ移動する。（※一列の要素数＝＝行数、一行の要素数＝＝列数）
OpenGL形式の表現の行列を表示する関数
関数本体
//! @brief 指定した行列を表示する関数
//! @param [in] m 行列を表す一次元配列
//! @param [in] format printf形式の書式
//! @param [in] I 行数
//! @param [in] J 列数
template<typename real_t> inline std::string mat_to_string(const real_t* m, const char* format,int I,int J) { std::stringstream s; char tmp[100]; for (int i = 0; i < I; i++) { for (int j = 0; j < J; j++) { sprintf(tmp, format, m[matijI(i, j,I)] ); s << tmp << " "; } s << "\n"; } return s.str(); }
使用例
//OpenGLの機能で変換行列を作成する glMatrixMode(GL_MODELVIEW); glTranslated(-1, 0, -2); //作成した変換行列を取得する GLfloat modelf[16]; glGetFloatv(GL_MODELVIEW_MATRIX, modelf); std::cout << mat_to_string(modelf, "%02.3f", 4, 4);
結果
数学的に表示できていることを確認する。
1.000 0.000 0.000 -1.000 0.000 1.000 0.000 0.000 0.000 0.000 1.000 -2.000 0.000 0.000 0.000 1.000
4×4の逆行列を求める関数
3DCGには欠かせない4×4行列の逆行列を求める関数を自前実装する。
動作チェックはEigenとの差を見て判断する。
逆行列を求める関数
namespace matcalc {
//誤差の定義 template<typename T> struct real_error {static constexpr float value = (float)1e-7;}; template<> struct real_error<double> { static constexpr double value = 1e-15; }; //! @brief 誤差の閾値未満かをチェック //! @param [in] value 検査する値 //! @retval true 値の絶対値は閾値より小さい //! @retval false 値の絶対値は閾値より大きい template<typename T> inline bool is_error(const T value) { return abs(value) < real_error<T>::value; }
//! @brief 4x4行列の行列式を取得 //! @param [in] dm16 4x4行列 //! @return 行列式 template<typename real_t> inline auto determinant44(const real_t* dm16) { using scalar_t = decltype(dm16[0]); const scalar_t a11 = dm16[0]; const scalar_t a12 = dm16[1]; const scalar_t a13 = dm16[2]; const scalar_t a14 = dm16[3]; const scalar_t a21 = dm16[4]; const scalar_t a22 = dm16[5]; const scalar_t a23 = dm16[6]; const scalar_t a24 = dm16[7]; const scalar_t a31 = dm16[8]; const scalar_t a32 = dm16[9]; const scalar_t a33 = dm16[10]; const scalar_t a34 = dm16[11]; const scalar_t a41 = dm16[12]; const scalar_t a42 = dm16[13]; const scalar_t a43 = dm16[14]; const scalar_t a44 = dm16[15]; return a11 * a22* a33* a44 + a11 * a23 * a34 * a42 + a11 * a24 * a32 * a43 + a12 * a21 * a34 * a43 + a12 * a23 * a31 * a44 + a12 * a24 * a33 * a41 + a13 * a21 * a32 * a44 + a13 * a22 * a34 * a41 + a13 * a24 * a31 * a42 + a14 * a21 * a33 * a42 + a14 * a22 * a31 * a43 + a14 * a23 * a32 * a41 - a11 * a22 * a34 * a43 - a11 * a23 * a32 * a44 - a11 * a24 * a33 * a42 - a12 * a21 * a33 * a44 - a12 * a23 * a34 * a41 - a12 * a24 * a31 * a43 - a13 * a21 * a34 * a42 - a13 * a22 * a31 * a44 - a13 * a24 * a32 * a41 - a14 * a21 * a32 * a43 - a14 * a22 * a33 * a41 - a14 * a23 * a31 * a42; }
//! @brief 4x4行列の逆行列を取得 //! @param [out] dst_dm16 結果を格納する一次元配列 //! @param [in] src_dm16 元の行列 //! @retval true 成功した //! @retval false 行列式の大きさが小さすぎて失敗した template<typename real_t> inline bool inverse_matrix44(real_t* dst_dm16, const real_t* src_dm16) { using s_scalar_t = decltype(src_dm16[0]); using d_scalar_t = decltype(dst_dm16[0]); real_t a11 = src_dm16[0]; real_t a12 = src_dm16[1]; real_t a13 = src_dm16[2]; real_t a14 = src_dm16[3]; real_t a21 = src_dm16[4]; real_t a22 = src_dm16[5]; real_t a23 = src_dm16[6]; real_t a24 = src_dm16[7]; real_t a31 = src_dm16[8]; real_t a32 = src_dm16[9]; real_t a33 = src_dm16[10]; real_t a34 = src_dm16[11]; real_t a41 = src_dm16[12]; real_t a42 = src_dm16[13]; real_t a43 = src_dm16[14]; real_t a44 = src_dm16[15]; real_t& b11 = dst_dm16[0]; real_t& b12 = dst_dm16[1]; real_t& b13 = dst_dm16[2]; real_t& b14 = dst_dm16[3]; real_t& b21 = dst_dm16[4]; real_t& b22 = dst_dm16[5]; real_t& b23 = dst_dm16[6]; real_t& b24 = dst_dm16[7]; real_t& b31 = dst_dm16[8]; real_t& b32 = dst_dm16[9]; real_t& b33 = dst_dm16[10]; real_t& b34 = dst_dm16[11]; real_t& b41 = dst_dm16[12]; real_t& b42 = dst_dm16[13]; real_t& b43 = dst_dm16[14]; real_t& b44 = dst_dm16[15]; b11 = a22 * a33 * a44 + a23 * a34 * a42 + a24 * a32 * a43 - a22 * a34 * a43 - a23 * a32 * a44 - a24 * a33 * a42; b12 = a12 * a34 * a43 + a13 * a32 * a44 + a14 * a33 * a42 - a12 * a33 * a44 - a13 * a34 * a42 - a14 * a32 * a43; b13 = a12 * a23 * a44 + a13 * a24 * a42 + a14 * a22 * a43 - a12 * a24 * a43 - a13 * a22 * a44 - a14 * a23 * a42; b14 = a12 * a24 * a33 + a13 * a22 * a34 + a14 * a23 * a32 - a12 * a23 * a34 - a13 * a24 * a32 - a14 * a22 * a33; b21 = a21 * a34 * a43 + a23 * a31 * a44 + a24 * a33 * a41 - a21 * a33 * a44 - a23 * a34 * a41 - a24 * a31 * a43; b22 = a11 * a33 * a44 + a13 * a34 * a41 + a14 * a31 * a43 - a11 * a34 * a43 - a13 * a31 * a44 - a14 * a33 * a41; b23 = a11 * a24 * a43 + a13 * a21 * a44 + a14 * a23 * a41 - a11 * a23 * a44 - a13 * a24 * a41 - a14 * a21 * a43; b24 = a11 * a23 * a34 + a13 * a24 * a31 + a14 * a21 * a33 - a11 * a24 * a33 - a13 * a21 * a34 - a14 * a23 * a31; b31 = a21 * a32 * a44 + a22 * a34 * a41 + a24 * a31 * a42 - a21 * a34 * a42 - a22 * a31 * a44 - a24 * a32 * a41; b32 = a11 * a34 * a42 + a12 * a31 * a44 + a14 * a32 * a41 - a11 * a32 * a44 - a12 * a34 * a41 - a14 * a31 * a42; b33 = a11 * a22 * a44 + a12 * a24 * a41 + a14 * a21 * a42 - a11 * a24 * a42 - a12 * a21 * a44 - a14 * a22 * a41; b34 = a11 * a24 * a32 + a12 * a21 * a34 + a14 * a22 * a31 - a11 * a22 * a34 - a12 * a24 * a31 - a14 * a21 * a32; b41 = a21 * a33 * a42 + a22 * a31 * a43 + a23 * a32 * a41 - a21 * a32 * a43 - a22 * a33 * a41 - a23 * a31 * a42; b42 = a11 * a32 * a43 + a12 * a33 * a41 + a13 * a31 * a42 - a11 * a33 * a42 - a12 * a31 * a43 - a13 * a32 * a41; b43 = a11 * a23 * a42 + a12 * a21 * a43 + a13 * a22 * a41 - a11 * a22 * a43 - a12 * a23 * a41 - a13 * a21 * a42; b44 = a11 * a22 * a33 + a12 * a23 * a31 + a13 * a21 * a32 - a11 * a23 * a32 - a12 * a21 * a33 - a13 * a22 * a31; auto detA = determinant44(src_dm16); if (is_error(detA) == true) return false; detA = static_cast<decltype(detA)>(1.0) / detA; for (int i = 0; i < 16; i++) { dst_dm16[i] *= detA; } return true; }
}
動作チェック
#include<Eigen/Dense> #include "inversematrix.hpp" //自前実装版 void inversematrix_check() { ///////////////////////////////////////////// // 元データ double base[16]; for (size_t i = 0; i < 16; i++) { base[i] = rand() % 1000 / 1000.0; }
///////////////////////////////////////////// // Eigenで求める Eigen::MatrixXd eg(4, 4); eg << base[0], base[1], base[2], base[3], base[4], base[5], base[6], base[7], base[8], base[9], base[10], base[11], base[12], base[13], base[14], base[15]; auto inv = eg.inverse(); std::cout << "eigen----------" << std::endl; std::cout << inv << std::endl; std::cout << "----------" << std::endl;
///////////////////////////////////////////// //自前実装版 double* A = new double[4 * 4]{ base[0], base[1], base[2], base[3] , base[4], base[5], base[6], base[7] , base[8], base[9], base[10], base[11] , base[12], base[13], base[14], base[15] }; double* B = new double[4 * 4]; matcalc::inverse_matrix44(B, A); #pragma warning(push) #pragma warning(disable:4996) std::cout << "自前----------" << std::endl; for (int i = 0; i < 4; i++) { for (int j = 0; j < 4; j++) { printf("%0.5lf", B[i*4+ j]); printf(" "); } puts(""); }
std::cout << "差--------" << std::endl; std::cout << "----------" << std::endl; for (int i = 0; i < 4; i++) { for (int j = 0; j < 4; j++) { printf("%+0.5lf", B[i * 4 + j]-inv(i,j)); printf(" "); } puts(""); }
std::cout << std::endl << std::endl << std::endl; #pragma warning(pop) }
テスト結果
DIBの４バイト境界について。あとDIBの保存。
BMPファイルは１スキャンラインを４バイト境界に合わせなければならないという所を前回確認した。
実はDIB Sectionでアクセスできる画素データは、この４バイト境界が考慮されている。つまり２４ビットDIBで幅が４バイトの倍数にない場合、画像内に４バイト境界に合わせるためのダミーデータが点在している。
これはこれでいいのだが、画素データを直接編集したいだけの場合は注意が必要で、普通に画像幅×y座標＋x座標の式で目的の画素のアドレスを求められないことになる。といっても画素幅が画素バイトになるだけなのだが。
４の倍数ということは素因数分解したときに２が二つ以上出てくるということで、この条件を満たす数はかなり多い。例えば１００は2^2×5^2なので、１００の倍数は全て４バイト境界の条件を満たす。だから実験として300x300とかの画像で確認などすると正常に動作してしまい、考慮漏れとなる可能性がある。言うまでもないことだが256や512など2のn乗系の数字は0,1,2以外全部当てはまる。
実験コード
このプログラムでは24bit DIBを作り画面に表示して保存している。画素のアドレスはCalcPosで行っているが、移動方法は
１スキャンラインバイト数×y座標 +１画素バイト数×x座標
としている。
なお１スキャンラインバイト数の計算とBMPファイルへの書き出しは前回作成したBMPwriter.hppの関数を使っている。
#include<windows.h> #include "BMPwriter.hpp"
BITMAPINFO bmpInfo; //!< CreateDIBSectionに渡す構造体 LPDWORD lpPixel; //!< 画素へのポインタ HBITMAP hBitmap; //!< 作成したMDCのビットマップハンドル HDC hMemDC; //!< 作成したMDCのハンドル //24bit DIB 作成 void createMemDC_24(HDC hdc,int width,int height) { //DIBの情報を設定する bmpInfo.bmiHeader.biSize = sizeof(BITMAPINFOHEADER); bmpInfo.bmiHeader.biWidth = width; bmpInfo.bmiHeader.biHeight = -height; //-を指定しないと上下逆になる bmpInfo.bmiHeader.biPlanes = 1; bmpInfo.bmiHeader.biBitCount = 24; bmpInfo.bmiHeader.biCompression = BI_RGB; hBitmap = CreateDIBSection(hdc, &bmpInfo, DIB_RGB_COLORS, (void**)& lpPixel, nullptr, 0); hMemDC = CreateCompatibleDC(hdc); SelectObject(hMemDC, hBitmap); }
//１スキャンラインのバイト数（何バイト移動すれば次の行へいけるのか）を計算 int get_scanline_bytes(const int width,const int bitcount) { return width * (bitcount / 8) + getbmp4line_err(width, bitcount); //getbmp4line_errはBMPwriter.hpp内 } //! @brief 与えられた座標のアドレスを返す //! @param [in] head 画像の先頭アドレス //! @param [in] width 画像幅（画素数） //! @param [in] height 画像高さ（画素数） //! @param [in] x x座標 //! @param [in] y y座標 //! @param [in] bitcount １画素のビット数 //! @return x,yの要素のアドレス unsigned char* calcPos( unsigned char* head, int width, int height, int x, int y, int bitcount) { int pixelbyte = bitcount / 8; int scanlinebyte = get_scanline_bytes(width, bitcount); unsigned char* yhead = head + scanlinebyte * y; return yhead + x * pixelbyte; }
//! @brief 画像を灰色で塗りつぶす //! @param [in] 対象となる画像の先頭アドレス //! @param [in] width 画像幅（画素数） //! @param [in] height 画像高さ（画素数） //! @param [in] bitcount DIBの１画素のビット数 //! @return なし void clear_data(unsigned char* img, int width, int height, int bitcount) { int scanlinebyte = get_scanline_bytes(width, bitcount); int size = scanlinebyte * height; for (size_t i = 0; i < size; i++) { img[i] = 200; } } //! @brief 画像の真ん中に線を引く //! @param [in] 対象となる画像の先頭アドレス //! @param [in] width 画像幅（画素数） //! @param [in] height 画像高さ（画素数） //! @param [in] bitcount DIBの１画素のビット数 //! @return なし void drawline_data(unsigned char* img, int width, int height,int bitcount) { for (int y = 0; y < height; y++) { for (int x = 0; x < width; x++) { unsigned char* pos = calcPos(img, width, height, x, y, bitcount); if (x == width / 2) { pos[0] = 255; pos[1] = 0; pos[2] = 0; } } } }
LRESULT CALLBACK WndProc(HWND hwnd, UINT msg, WPARAM wp, LPARAM lp) { PAINTSTRUCT ps; static HPEN hpen1; static HPEN hpen2; static HPEN open; int iwidth = 4*30+1; //4の倍数にならない値を作成 int iheight = 80; switch (msg) { case WM_CREATE: { HDC hdc = GetDC(hwnd); createMemDC_24(hdc, iwidth, iheight);//24bit DIB 作成 ReleaseDC(hwnd, hdc); //DIBを灰色でクリア clear_data((unsigned char*)lpPixel, iwidth, iheight, 24); hpen1 = CreatePen(PS_SOLID, 5, RGB(255, 0, 0));// 何かしら描画する hpen2 = CreatePen(PS_SOLID, 3, RGB(0, 0, 255)); open = (HPEN)SelectObject(hMemDC, hpen1); MoveToEx(hMemDC, 0, 0, nullptr); LineTo(hMemDC, 20, 50); SelectObject(hMemDC, hpen2); MoveToEx(hMemDC, 0, 0, nullptr); LineTo(hMemDC, 50, 20); SelectObject(hMemDC, open); drawline_data((unsigned char*)lpPixel, iwidth, iheight,24);//縦に一本線を入れる。座標計算が間違っていたら斜めの線になる
//bmpファイルに保存。lpPixelは4byte境界が考慮されているのでそのまま出力できる。bmp_writeはBMPwriter.hpp内 bmp_write(L"c:\\data\\a.bmp", lpPixel, iwidth, iheight, 24);
} return 0; case WM_PAINT: { HDC hdc = BeginPaint(hwnd, &ps); BitBlt(hdc, 0, 0, iwidth, iheight, hMemDC, 0, 0, SRCCOPY);// 画像を画面に表示 EndPaint(hwnd, &ps); } return 0; case WM_DESTROY: DeleteObject(hpen1); DeleteObject(hpen2); PostQuitMessage(0); return 0; } return DefWindowProc(hwnd, msg, wp, lp); }
//! @brief エントリポイント //! @param [in] hInstance 現在のインスタンスハンドル //! @param [in] hPrevInstance 以前のインスタンスハンドル //! @param [in] lpCmdLine コマンドライン引数の文字列 //! @param [in] ウィンドウの表示状態 int WINAPI WinMain( HINSTANCE hInstance, HINSTANCE hPrevInstance, PSTR lpCmdLine, int nCmdShow) { WNDCLASSEX wc; wc.cbSize = sizeof(WNDCLASSEX); wc.style = CS_HREDRAW | CS_VREDRAW; wc.lpfnWndProc = WndProc; wc.cbClsExtra = 0; wc.cbWndExtra = 0; wc.hInstance = hInstance; wc.hIcon = LoadIcon(NULL, IDI_APPLICATION); wc.hCursor = LoadCursor(NULL, IDC_ARROW); wc.hbrBackground = (HBRUSH)GetStockObject(WHITE_BRUSH); wc.lpszMenuName = NULL; wc.lpszClassName = TEXT("SZLCODE"); wc.hIconSm = LoadIcon(NULL, IDI_APPLICATION); if (!RegisterClassEx(&wc)) return 0; const int top = CW_USEDEFAULT; const int left = CW_USEDEFAULT; const int width = 300; const int height = 200; HWND hwnd; hwnd = CreateWindowEx( WS_EX_COMPOSITED, TEXT("SZLCODE"), TEXT("win32api example"), WS_OVERLAPPEDWINDOW | WS_VISIBLE, top, left, width, height, NULL, NULL, hInstance, NULL ); if (hwnd == NULL) return 0; // MSG msg; while (GetMessage(&msg, NULL, 0, 0)) { TranslateMessage(&msg); DispatchMessage(&msg); } return (int)msg.wParam; }
WindowsからBMPファイルを出力
BMPファイル形式について
BMPはMicrosoftとIBMが開発したファイル形式で、Windows版とOS/2版がある。ここではWindows版の32bit,24bitしか考えない。
以下Wikipedia:
https://ja.wikipedia.org/wiki/Windows_bitmap
BMPファイルの４バイト境界について
BMPファイルは、横一行のデータのバイト数が4の倍数でなければいけないという制約がある。32bit BMPの場合は１画素が4バイトなので必ずこの条件を満たす。従って考える必要はない（考える必要がないだけで制約が取り払われているわけではない）。
問題は24bitの時で、例えば3×2の画像の一行の画素数は3なので、１行3×(24/8)=9バイトになる。9より大きいもっとも小さな4の倍数は12なので、10,11,12バイト目はダミーのデータで埋めてやる必要がある。
しかも、32bit BMPは24bit BMPに比べて扱えないソフトが多いとの理由で、.bmp形式なら24bit BMPが世の中では推奨されている。だから24bit BMPを無視できない。
BMPファイルの出力手順
１．データが4byte境界の条件を満たしていないなら、満たしているデータを作成する。
２．BITMAPFILEHEADER構造体に値を設定する
３．BITMAPINFO構造体に値を設定する
４．上記２，３をバイナリでファイルに書き出す
５．画像データをファイルに書き出す
BMPファイル形式はMicrosoftが考案した形式だが、Windowsに用意されているのはヘッダの構造体だけだ。データの4byte境界を修正する関数はおろかそのエラーチェックをする関数も用意されていない。圧縮もサポートしていると言ってもヘッダに圧縮状態のフラグを指定できるだけで、使いたかったら自分で圧縮しなければならない。ヘッダの出力もやってることはただのfwriteで、sizeof(BITMAPFILEHEADER)バイトだけ、sizeof(BITMAPINFO)バイトだけバイナリ出力しているに過ぎない。
サンプルコード
BMPwriter.hpp
#pragma once //! @brief 与えられた画素数がビットマップの要求する4byte境界に何バイト足りないかを返す //! @param [in] width 画素数 //! @param [in] pixelsize １ピクセルのビット数 //! @return 足して４の倍数になる値 int getbmp4line_err(const unsigned int width, const unsigned int pixelbits); //! @brief 幅と高さと画素のビット数からビットマップデータにしたときに何バイトになるかを算出する //! @param [in] width 画像のピクセル数（幅） //! @param [in] height 画像のピクセル数（高さ） //! @param [in] pixelbit １ピクセルのビット数 //! @return 4バイト境界に合わせたときのバイト数 unsigned int getbmpdata_memsize(const unsigned int width, const unsigned int height, const unsigned int pixelbits); //! @brief データをスキャンライン4倍バイトデータにしたものを作成 //! @param [out] dst 調整後の画素データの格納先アドレス //! @param [in] src 元データの先頭アドレス //! @param [in] width 画像幅（ピクセル） //! @param [in] height 画像高さ（ピクセル） //! @param [in] pixelbit １ピクセルのビット数 //! @return そのままBMPにできるデータ配列 //! @note dstは確保済みの必要がある。このバイト数はgetbmpdatasizeで取得できる void data_to_bmp4lines(unsigned char* dst,const void* src, const unsigned int width, const unsigned int height, const unsigned int pixelbits); //! @brief ビットマップファイルを保存する //! @param [in] pathname ファイル名 //! @param [in] pixeldata 画像データの先頭ポインタ。４バイト境界に調整済みの必要がある //! @param [in] width 画像のピクセル幅 //! @param [in] height 画像のピクセル高さ //! @param [in] １画素のビット数 //! @return なし void bmp_write(const wchar_t* pathname, const void* pixeldata, const unsigned int width, const unsigned int height, const unsigned int pixelbits); //! @brief ビットマップファイルを保存する関数のインタフェース //! @param [in] pathname ファイル名 //! @param [in] pixeldata 画像データの先頭ポインタ。４バイト境界の調整前のもの //! @param [in] width 画像のピクセル幅 //! @param [in] height 画像のピクセル高さ //! @param [in] １画素のビット数 //! @return なし void call_bmp_write(const wchar_t* pathname, const void* pixeldata, const unsigned int width, const unsigned int height, const unsigned int pixelbits);
BMPWriter.cpp
#include "BMPwriter.hpp" #include <Windows.h> #include <vector>
int getbmp4line_err(const unsigned int width, const unsigned int pixelbits) { return (4 - (width * pixelbits /8) % 4) % 4; }
unsigned int getbmpdata_memsize(const unsigned int width, const unsigned int height, const unsigned int pixelbits) { unsigned int widthbytes = width * (pixelbits / 8) + getbmp4line_err(width, pixelbits); return widthbytes * height; }
void data_to_bmp4lines(unsigned char* dst, const void* src, const unsigned int width, const unsigned int height, const unsigned int pixelbits) { unsigned int pixelbytes = (pixelbits / 8); std::uint8_t * p = (std::uint8_t*)src; int err = getbmp4line_err(width, pixelbits);//4バイト境界に何バイト足りないかを算出 size_t pos=0; for (size_t h = 0; h < height; h++) { //一行コピー for (size_t x = 0; x < width * pixelbytes; x++) { dst[pos++] = *p; p++; } if (err != 0) { // width * pixelsizeが4になるようにデータを挿入 for (size_t k = 0; k < err; k++) { dst[pos++] = 0; } } } }
void bmp_write( const wchar_t* pathname, const void* pixeldata, const unsigned int width, const unsigned int height, const unsigned int pixelbits) { const unsigned long bytecount = pixelbits / 8; const unsigned long headSize = (sizeof(BITMAPFILEHEADER) + sizeof(BITMAPINFO)); //４バイト境界を考慮した１スキャンラインのバイト数 const unsigned long widthBytes = width * bytecount + getbmp4line_err(width, pixelbits); const unsigned long imageBytes = (widthBytes * height); // BITMAPFILEHEADERの初期化 BITMAPFILEHEADER bmpHead = { 0 }; bmpHead.bfType = 0x4D42; // "BM" bmpHead.bfSize = headSize + imageBytes; bmpHead.bfReserved1 = 0; bmpHead.bfReserved2 = 0; bmpHead.bfOffBits = headSize; // BITMAPINFOの初期化 BITMAPINFO bmpInfo = { 0 }; bmpInfo.bmiHeader.biSize = sizeof(BITMAPINFOHEADER); bmpInfo.bmiHeader.biWidth = width; bmpInfo.bmiHeader.biHeight = height; bmpInfo.bmiHeader.biPlanes = 1; bmpInfo.bmiHeader.biBitCount = pixelbits; bmpInfo.bmiHeader.biCompression = BI_RGB; bmpInfo.bmiHeader.biSizeImage = 0; bmpInfo.bmiHeader.biXPelsPerMeter = 0; bmpInfo.bmiHeader.biYPelsPerMeter = 0; bmpInfo.bmiHeader.biClrUsed = 0; bmpInfo.bmiHeader.biClrImportant = 0; FILE * fp = _wfopen(pathname, L"wb" ); fwrite(&bmpHead, sizeof(BITMAPFILEHEADER), 1, fp); fwrite(&bmpInfo, sizeof(BITMAPINFO), 1, fp); fwrite(pixeldata, 1 , imageBytes, fp); fclose(fp); }
void call_bmp_write(const wchar_t* pathname, const void* pixeldata, const unsigned int width, const unsigned int height, const unsigned int pixelbits) { //BMP保存時の画像のメモリサイズを求める int retsize = getbmpdata_memsize(width, height, pixelbits); std::vector<std::uint8_t> save; save.resize(retsize); //４バイト境界の条件を満たしたデータを作成する data_to_bmp4lines(&save[0], pixeldata, width, height, pixelbits); //BMPファイルに保存する bmp_write(pathname, &save[0], width, height, pixelbits); }
使用例
#include <iostream> #include <vector> #include "BMPwriter.hpp" struct data32bit { unsigned char p[4]; data32bit() {} data32bit( unsigned char b, unsigned char g, unsigned char r, unsigned char u ) { p[0] = b; p[1] = g; p[2] = r; p[3] = u; } data32bit( unsigned char b, unsigned char g, unsigned char r ) { p[0] = b; p[1] = g; p[2] = r; p[3] = 0; } }; struct data24bit { unsigned char p[3]; data24bit() {} data24bit( unsigned char b, unsigned char g, unsigned char r ) { p[0] = b; p[1] = g; p[2] = r; } }; template<typename dtype> void create_data(std::vector<dtype>& img,int width,int height) { for (int y = 0; y < height; y++) { for (int x = 0; x < width; x++) { int pos = y * width + x; if (y == x) { img[pos] = dtype(255, 0, 0); } else { img[pos] = dtype(0, 0, 255); } if (x == width - 1) img[pos] = dtype(0, 255, 0); } } }
int main() { //画像データ std::vector<data24bit> img; unsigned int bits = 24; unsigned int width = 7; unsigned int height = 9; img.resize(width*height); create_data(img,width,height); std::cout << "scanline bytes % 4 = " << (width * bits / 8) % 4; call_bmp_write(L"C:\\data\\a.bmp",&img[0], width, height, bits); int k; std::cin >> k; }
Blender 2.8 How to Create the World in 1 Minute 抜粋 2-シェーダの設定 + 雲
前回の続き。まずざっくり以下のようにノードを設定する。
次に、夜景のテクスチャを追加し、Shader to RGBとMultiplyしてEmissionへつなげる
雲
[Shift+D]で地球を複製し、僅かに拡大する。このままだとマテリアルが共通なので、数字をクリックして新しいマテリアルを作成する。
複製したほうの球体に対して雲のテクスチャを設定する。雲のテクスチャは恐らく雲＝白、その他＝黒になっているので、そのままAlphaに接続すると雲だけがレンダリングできる。
結果(アニメーションPng)
Blender 2.8 How to Create the World in 1 Minute 抜粋 1-モデル作成
使用テクスチャ
テクスチャ
http://www.shadedrelief.com/natural3/pages/textures.html
使用アドオン
Planeを配置
１．Image as Planeで画像を貼り付けたPlaneを作成
２．Editモードへ移行し、[R][Y][90]でY軸に90度回転
３．Ctr+RでX軸に水平に一回ループカット
４．Subdivision Surfaceモディファイアを追加し6,simpleに設定
５．Simple Deformモディファイアを追加し、Bend,Axis = X,Angle=180
６．もう一つSimple Deformモディファイアを追加し、Bend,Axis=Z,Angle=360
７．モディファイアを適用し、Editモードで[Mesh]→[Clean Up]→[Decimate Geometry]で重複点をmergeする。
続く
Blender で地面を作る（テクスチャ使用）
以下はHow to make Realistic Ground の内容。＋α。
1.基本
1-1.テクスチャのダウンロード
まずテクスチャを準備する。わかりやすいように石が多いものがいい。また、displacementマップを持っている必要がある。
今回は以下を使用：
https://cc0textures.com/view?id=Ground022
1-2.Planeの設置
[Shift + A]→[Mesh]→[Plane]で平面を追加
[Tab]でエディットモードに入り、右クリック→[Subdivide]を6回ほど繰り返しPlaneを細かく分割する
1-3.マテリアルの設定
[Shift+A]→[Texture]→[Image Texture]でテクスチャノードを追加し、カラー画像を設定する
1-4.Displaceモディファイアの設定
1-4-１．テクスチャの設定
Displaceモディファイアを追加し、テクスチャに、マテリアルに設定したテクスチャのDisplacementテクスチャを追加する
1-4-2.調整、細分化
Displaceモディファイアの設定に戻り、Strengthを0.04程度にする。
さらにSubdivisionモディファイア追加し、Subdivisionsの設定を4等の大きい値にする。
2.応用 拡大・縮小
2-1.displacementモディファイアの設定
この方法では、マテリアルに設定したカラーのテクスチャと、displaceモディファイアに設定したDisplacementテクスチャの座標系が１：１に対応していなければならない。マテリアル側とDisplace側で座標系がやや違うのでそれを考慮しながら作業する
2-1-1.Planeの追加
[Shift+A]→[Mesh]→[Plane]でplaneを追加する（Plane.001）。X,Y方向の位置が重要になるので、適当に移動させずに、真下に移動させるか別コレクションとしておく
2-1-2.DisplaceモディファイアにPlane.001を設定
Displaceモディファイアの[Texture Coordinates]に追加したPlane.001を指定
このままPlane.001をスケーリングするとDisplaceも拡大縮小するが、スケーリングの中心がPlaneの中央になっている。テクスチャ側のスケーリング中心が左端なので合わせる必要がある。
2-1-3.Plane.001を(-1,-1)へ移動
Plane.001を-1,-1へ移動する。z座標は邪魔にならなければいくつでも良い。
2-1-4.テクスチャの座標系を変更
[Shift+A]→[Vector]→[Mapping]を追加し、Locationを0.5,0.5に設定
Texture CoordinateはUVから繋いでおく。
これで、再びDisplaceモディファイアとカラーのテクスチャが一致するようになる
2-1-5.テクスチャとdisplaceをスケーリング
Plane.001側を0.5,0.5倍、テクスチャ側を2.0,2.0倍にする。
3.その他
Displaceモディファイアをもう一つ追加し、Cloudsなどを設定するとよりリアルな地形になる。
主成分分析で三次元点群から平面を求める – 2 – 固有ベクトルを求める
前回はデータから分散共分散行列を求めた。
分散共分散行列（3×3）に対してヤコビ法などを用いると、行列の固有値と固有ベクトルがそれぞれ３つずつ求まる。固有値が最も小さい固有ベクトルが面法線となる。
プログラム
リンク先、先頭のヤコビ法のプログラムをそのまま使用する
固有ベクトルと固有値
該当部分プログラム
printf("分散共分散行列\n"); printf("% 1.5lf % 1.5lf % 1.5lf\n", Sxx, Sxy, Sxz); printf("% 1.5lf % 1.5lf % 1.5lf\n", Sxy, Syy, Syz); printf("% 1.5lf % 1.5lf % 1.5lf\n", Sxz, Syz, Szz); puts(""); float a[3 * 3] = { Sxx,Sxy,Sxz, Sxy,Syy,Syz, Sxz,Syz,Szz }; float v[3 * 3]; EigenJacobiMethod(a, v, 3); printf("対角要素固有値\n"); printf("% 1.5lf % 1.5lf % 1.5lf\n", a[0], a[1], a[2]); printf("% 1.5lf % 1.5lf % 1.5lf\n", a[3], a[4], a[5]); printf("% 1.5lf % 1.5lf % 1.5lf\n", a[6], a[7], a[8]); puts(""); printf("固有ベクトル\n"); printf("% 1.5lf % 1.5lf % 1.5lf\n", v[0], v[1], v[2]); printf("% 1.5lf % 1.5lf % 1.5lf\n", v[3], v[4], v[5]); printf("% 1.5lf % 1.5lf % 1.5lf\n", v[6], v[7], v[8]);
面法線の選択
主成分分析で得られるベクトルは、データがどの方向にばらついているかを表していて、一番ばらつきが少ない方向が面法線という事になる
参考：
主成分分析の考え方
ので、一番小さい固有値に対応した固有ベクトルが面法線になる。
参考：
”主成分分析における固有値は、その固有ベクトルの方向に沿ったデータの分散の大きさに対応します。固有値が大きい固有ベクトルほど、データの分散をよく説明している、つまり、データの重要な特徴を捉えていると考えられます。なので、各固有値を固有値の合計で割ったものを寄与率と呼んだり、固有値（寄与率）が大きい主成分から順に、第１主成分（PC1）、第２主成分（PC2）……と呼んだりします。”
上記結果では固有値が0.00072のベクトル( 0.00957 , -0.00150 , 0.99995 )が答えになる。
面の特定
面はその面に含まれる一点と法線ベクトルがわかれば求められる。
面に含まれる一点・・・頂点群の重心
面法線・・・頂点群の主成分分析の結果の第三主成分
で面が決定できる。
プログラム全体
#include <iostream> #include <vector> #include <array> //! @struct point //! @brief 座標を表すデータ型 struct point { double p[3]; point(double x, double y, double z) { p[0] = x; p[1] = y; p[2] = z; } point() {} }; //! @brief 固有値・固有ベクトル struct eigen { double val; //!< 固有値 double vec[3];//!< 固有ベクトル }; int EigenJacobiMethod(float* a, float* v, int n, float eps = 1e-8, int iter_max = 100); void data(std::vector<point>& pt);
//! @brief データの平均を求める //! @param [in] vec データの配列 //! @return 各要素の平均 point Covariance_Ave(const std::vector<point>& vec) { // 初期化 point ave; ave.p[0] = 0; ave.p[1] = 0; ave.p[2] = 0; // 各要素平均 for (size_t i = 0; i < vec.size(); i++) { ave.p[0] += vec[i].p[0]; ave.p[1] += vec[i].p[1]; ave.p[2] += vec[i].p[2]; } ave.p[0] /= vec.size(); ave.p[1] /= vec.size(); ave.p[2] /= vec.size(); return ave; }
//! @brief 共分散を求める //! @param [in] average 各要素の平均 //! @param [in] vec データ //! @param [in] param どの要素に対して求めるか。例えばxyzの時、x,yに対する共分散なら{0,1}を与える。 //! @return 偏差の積の和の要素数分の一 double Covariance(const point& average,const std::vector<point>& vec,const std::array<int,2>& param) { double sum = 0.0; for (size_t i = 0; i < vec.size(); i++) { //指定したパラメータの偏差を求める point deviation; for (size_t j = 0; j < param.size(); j++) { int target = param[j]; deviation.p[target] = (vec[i].p[target] - average.p[target]); } //偏差の積 double product = 1.0; for (size_t j = 0; j < param.size(); j++) { int target = param[j]; product *= deviation.p[target]; } //偏差の積の和を更新 sum += product; } //偏差の積の和のN分の一 return 1.0 / vec.size() * sum; }
//! @brief データを主成分分析する //! @param [out] ev 固有値と固有ベクトル //! @param [in] average 全データの各要素の平均 //! @param [in] pt 三次元の座標値一覧 //! @return なし void eigen_vector(eigen* ev, point* average, const std::vector<point>& pt) { //各要素の平均を求める。 //これは共分散を求めるときに (x[i] - xの平均)×(y[i] - yの平均) 等の計算が必要なため *average = Covariance_Ave(pt); //共分散を求める //第三引数の{0,0}はxxを表す。xyなら{0,1}。これはデータがxyzの順に並んでいる事が前提。 double Sxx = Covariance(*average, pt, { 0, 0 }); double Sxy = Covariance(*average, pt, { 0, 1 }); double Sxz = Covariance(*average, pt, { 0, 2 }); double Syy = Covariance(*average, pt, { 1, 1 }); double Syz = Covariance(*average, pt, { 1, 2 }); double Szz = Covariance(*average, pt, { 2, 2 }); // 分散共分散行列 float a[3 * 3] = { Sxx,Sxy,Sxz, Sxy,Syy,Syz, Sxz,Syz,Szz }; float v[3 * 3]; EigenJacobiMethod(a, v, 3); ev[0].val = a[0]; //固有値 ev[0].vec[0] = v[0];//固有値に対応した固有ベクトル ev[0].vec[1] = v[3]; ev[0].vec[2] = v[6]; ev[1].val = a[4]; //固有値 ev[1].vec[0] = v[1];//固有値に対応した固有ベクトル ev[1].vec[1] = v[4]; ev[1].vec[2] = v[7]; ev[2].val = a[8]; //固有値 ev[2].vec[0] = v[2];//固有値に対応した固有ベクトル ev[2].vec[1] = v[5]; ev[2].vec[2] = v[8]; }
int main() { std::vector<point> pt; data( pt ); //データの準備 eigen ev[3]; point average; eigen_vector(ev, &average,pt); //固有値と固有ベクトルを表示 printf("[%lf] %lf %lf %lf\n", ev[0].val, ev[0].vec[0], ev[0].vec[1], ev[0].vec[2]); printf("[%lf] %lf %lf %lf\n", ev[1].val, ev[1].vec[0], ev[1].vec[1], ev[1].vec[2]); printf("[%lf] %lf %lf %lf\n", ev[2].val, ev[2].vec[0], ev[2].vec[1], ev[2].vec[2]); printf("\n重心\n"); printf("%.5lf %.5lf %.5lf\n", average.p[0], average.p[1], average.p[2]); int i; std::cin >> i; }
void data(std::vector<point>& pt) { pt.emplace_back(-0.9282366037368774, -4.60659646987915, 4.179419040679932); pt.emplace_back(-0.6398676633834839, -4.345795631408691, 3.8638782501220703); pt.emplace_back(-0.3469808101654053, -4.115365028381348, 3.5088584423065186); pt.emplace_back(-0.05463826656341553, -3.881274938583374, 3.1585941314697266); pt.emplace_back(0.2328449785709381, -3.614520311355591, 2.850792646408081); pt.emplace_back(-1.0035791397094727, -3.9636778831481934, 3.8052620887756348); pt.emplace_back(-0.7170133590698242, -3.690755844116211, 3.5054771900177); pt.emplace_back(-0.4283676743507385, -3.4318137168884277, 3.1875195503234863); pt.emplace_back(-0.13146033883094788, -3.2284107208251953, 2.7973649501800537); pt.emplace_back(0.15122388303279877, -2.929394483566284, 2.531501293182373); pt.emplace_back(-1.0883535146713257, -3.2573564052581787, 3.513524055480957); pt.emplace_back(-0.7981395721435547, -3.0089566707611084, 3.181861639022827); pt.emplace_back(-0.5119532346725464, -2.733482599258423, 2.88539457321167); pt.emplace_back(-0.21177172660827637, -2.5520896911621094, 2.466628313064575); pt.emplace_back(0.08094996213912964, -2.3205485343933105, 2.113051414489746); pt.emplace_back(-1.1616857051849365, -2.6279513835906982, 3.1217992305755615); pt.emplace_back(-0.8680680990219116, -2.402432918548584, 2.7603931427001953); pt.emplace_back(-0.5798090696334839, -2.140892267227173, 2.4458136558532715); pt.emplace_back(-0.2915937602519989, -1.87905752658844, 2.1316158771514893); pt.emplace_back(0.002897977828979492, -1.6594164371490479, 1.7625699043273926); pt.emplace_back(-1.2396777868270874, -1.967220664024353, 2.7707955837249756); pt.emplace_back(-0.9490185976028442, -1.7218148708343506, 2.43524169921875); pt.emplace_back(-0.6570032835006714, -1.4855259656906128, 2.087836503982544); pt.emplace_back(-0.36362549662590027, -1.2583959102630615, 1.7285257577896118); pt.emplace_back(-0.0764693021774292, -0.9894416332244873, 1.4235832691192627); }
/*! * Jacobi法による固有値の算出 * @param[inout] a 実対称行列．計算後，対角要素に固有値が入る * @param[out] v 固有ベクトル(aと同じサイズ) * @param[in] n 行列のサイズ(n×n) * @param[in] eps 収束誤差 * @param[in] iter_max 最大反復回数 * @return 反復回数 * @see http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?%B8%C7%CD%AD%C3%CD/%B8%C7%CD%AD%A5%D9%A5%AF%A5%C8%A5%EB */ int EigenJacobiMethod(float* a, float* v, int n, float eps, int iter_max) { float* bim, * bjm; float bii, bij, bjj, bji; bim = new float[n]; bjm = new float[n]; for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { v[i * n + j] = (i == j) ? 1.0 : 0.0; } } int cnt = 0; for (;;) { int i, j; float x = 0.0; for (int ia = 0; ia < n; ++ia) { for (int ja = 0; ja < n; ++ja) { int idx = ia * n + ja; if (ia != ja && fabs(a[idx]) > x) { i = ia; j = ja; x = fabs(a[idx]); } } } float aii = a[i * n + i]; float ajj = a[j * n + j]; float aij = a[i * n + j]; float alpha, beta; alpha = (aii - ajj) / 2.0; beta = sqrt(alpha * alpha + aij * aij); float st, ct; ct = sqrt((1.0 + fabs(alpha) / beta) / 2.0); // sinθ st = (((aii - ajj) >= 0.0) ? 1.0 : -1.0) * aij / (2.0 * beta * ct); // cosθ // A = PAPの計算 for (int m = 0; m < n; ++m) { if (m == i || m == j) continue; float aim = a[i * n + m]; float ajm = a[j * n + m]; bim[m] = aim * ct + ajm * st; bjm[m] = -aim * st + ajm * ct; } bii = aii * ct * ct + 2.0 * aij * ct * st + ajj * st * st; bij = 0.0; bjj = aii * st * st - 2.0 * aij * ct * st + ajj * ct * ct; bji = 0.0; for (int m = 0; m < n; ++m) { a[i * n + m] = a[m * n + i] = bim[m]; a[j * n + m] = a[m * n + j] = bjm[m]; } a[i * n + i] = bii; a[i * n + j] = bij; a[j * n + j] = bjj; a[j * n + i] = bji; // V = PVの計算 for (int m = 0; m < n; ++m) { float vmi = v[m * n + i]; float vmj = v[m * n + j]; bim[m] = vmi * ct + vmj * st; bjm[m] = -vmi * st + vmj * ct; } for (int m = 0; m < n; ++m) { v[m * n + i] = bim[m]; v[m * n + j] = bjm[m]; } float e = 0.0; for (int ja = 0; ja < n; ++ja) { for (int ia = 0; ia < n; ++ia) { if (ia != ja) { e += fabs(a[ja * n + ia]); } } } if (e < eps) break; cnt++; if (cnt > iter_max) break; } delete[] bim; delete[] bjm; return cnt; }
Blenderで可視化
もっと良いツールがあるのかもしれないがＢｌｅｎｄｅｒで重心から各主成分に線を伸ばしたものを作成して結果を確認する
主成分分析で三次元点群から平面を求める – 1 – 分散共分散行列を求める
まず、手順としては下の図のようになる。
頂点データを行列として扱う所は特に問題は無いが、行列が苦手な人（自分もだが）のためにコメントすると3×頂点数の行列ができる。つまり例えば3×50000の行列とかにする。コード的な視点から見るとわざわざ「行列型」にする必要はない。構造体の配列で問題ない。
その行列に対して処理をして、最終的に3×3の分散共分散行列を作成する。その行列の固有ベクトル（ xyzのベクトルが三つ ）を求めればいい。
データ
例えば以下のようなデータの場合、分散共分散行列は3x3の行列になる。
| x | y | z |
|---|---|---|
| -0.9282366037368774 | -4.60659646987915 | 4.179419040679932 |
| -0.6398676633834839 | -4.345795631408691 | 3.8638782501220703 |
| -0.3469808101654053 | -4.115365028381348 | 3.5088584423065186 |
| -0.05463826656341553 | -3.881274938583374 | 3.1585941314697266 |
| 0.2328449785709381 | -3.614520311355591 | 2.850792646408081 |
| -1.0035791397094727 | -3.9636778831481934 | 3.8052620887756348 |
| -0.7170133590698242 | -3.690755844116211 | 3.5054771900177 |
| -0.4283676743507385 | -3.4318137168884277 | 3.1875195503234863 |
| -0.13146033883094788 | -3.2284107208251953 | 2.7973649501800537 |
| 0.15122388303279877 | -2.929394483566284 | 2.531501293182373 |
| -1.0883535146713257 | -3.2573564052581787 | 3.513524055480957 |
| -0.7981395721435547 | -3.0089566707611084 | 3.181861639022827 |
| -0.5119532346725464 | -2.733482599258423 | 2.88539457321167 |
| -0.21177172660827637 | -2.5520896911621094 | 2.466628313064575 |
| 0.08094996213912964 | -2.3205485343933105 | 2.113051414489746 |
| -1.1616857051849365 | -2.6279513835906982 | 3.1217992305755615 |
| -0.8680680990219116 | -2.402432918548584 | 2.7603931427001953 |
| -0.5798090696334839 | -2.140892267227173 | 2.4458136558532715 |
| -0.2915937602519989 | -1.87905752658844 | 2.1316158771514893 |
| 0.002897977828979492 | -1.6594164371490479 | 1.7625699043273926 |
| -1.2396777868270874 | -1.967220664024353 | 2.7707955837249756 |
| -0.9490185976028442 | -1.7218148708343506 | 2.43524169921875 |
| -0.6570032835006714 | -1.4855259656906128 | 2.087836503982544 |
| -0.36362549662590027 | -1.2583959102630615 | 1.7285257577896118 |
| -0.0764693021774292 | -0.9894416332244873 | 1.4235832691192627 |
分散共分散行列
分散共分散行列の形
データがx,y,zの場合、以下のような形をしている
問題はこのSxx,Sxy,Sxz,...をどうやって求めるかなのだが、このSxyは「xとyの共分散」である。
だから「xとxの共分散（=xの分散）」、「xとyの共分散」、...を求める。では共分散はどうやって求めるのか。
共分散の求め方
当然だが、Sxxの場合、「Sxx = (x[i]の偏差)の二乗の合計×(1/要素数)」になる。
偏差の求め方がわかれば共分散が求められ、共分散を求めることができれば分散共分散行列を組み立てることができる。
Excelで動作確認
分析ツールの追加
プログラムに起こしたので、Excelに同じデータを入れて正しい答えが出ているかを確認する。
Excelを開き、ファイル→アドオン→分析ツール→設定→分析ツール→ＯＫ
で分析ツールを追加し、データ→データ分析→共分散と操作する。
Excelによる計算結果
自作コードでの計算結果
ソースコード
#include <iostream> #include <vector> #include <array> //! @brief 座標を表すデータ型 struct point { double p[3]; point(double x, double y, double z) { p[0] = x; p[1] = y; p[2] = z; } point() {} };
//! @brief データの平均を求める //! @param [in] vec データの配列 //! @return 各要素の平均 point Covariance_Ave(const std::vector<point>& vec) { // 初期化 point ave; ave.p[0] = 0; ave.p[1] = 0; ave.p[2] = 0; // 各要素平均 for (size_t i = 0; i < vec.size(); i++) { ave.p[0] += vec[i].p[0]; ave.p[1] += vec[i].p[1]; ave.p[2] += vec[i].p[2]; } ave.p[0] /= vec.size(); ave.p[1] /= vec.size(); ave.p[2] /= vec.size(); return ave; }
//! @brief 共分散を求める //! @param [in] average 各要素の平均 //! @param [in] vec データ //! @param [in] param どの要素に対して求めるか。例えばxyzの時、x,yに対する共分散なら{0,1}を与える。 //! @return 偏差の積の和の要素数分の一 double Covariance(const point& average,const std::vector<point>& vec,const std::array<int,2>& param) { double sum = 0.0; for (size_t i = 0; i < vec.size(); i++) { //指定したパラメータの偏差を求める point deviation; for (size_t j = 0; j < param.size(); j++) { int target = param[j]; deviation.p[target] = (vec[i].p[target] - average.p[target]); } //偏差の積 double product = 1.0; for (size_t j = 0; j < param.size(); j++) { int target = param[j]; product *= deviation.p[target]; } //偏差の積の和を更新 sum += product; } //偏差の積の和のN分の一 return 1.0 / vec.size() * sum; }
void data(std::vector<point>& pt);
int main() { std::vector<point> pt; data( pt ); //データの準備 //各要素の平均を求める。 //これは共分散を求めるときに (x[i] - xの平均)×(y[i] - yの平均) 等の計算が必要なため point average = Covariance_Ave(pt); //共分散を求める //第三引数の{0,0}はxxを表す。xyなら{0,1}。これはデータがxyzの順に並んでいる事が前提。 double Sxx = Covariance(average, pt, { 0, 0 }); double Sxy = Covariance(average, pt, { 0, 1 }); double Sxz = Covariance(average, pt, { 0, 2 }); double Syy = Covariance(average, pt, { 1, 1 }); double Syz = Covariance(average, pt, { 1, 2 }); double Szz = Covariance(average, pt, { 2, 2 }); // 分散共分散行列を表示 printf("分散共分散行列\n"); printf("%.5lf %.5lf %.5lf\n", Sxx, Sxy, Sxz); printf("%.5lf %.5lf %.5lf\n", Sxy, Syy, Syz); printf("%.5lf %.5lf %.5lf\n", Sxz, Syz, Szz); int i; std::cin >> i; }
void data(std::vector<point>& pt) { pt.emplace_back(-0.9282366037368774, -4.60659646987915, 4.179419040679932); pt.emplace_back(-0.6398676633834839, -4.345795631408691, 3.8638782501220703); pt.emplace_back(-0.3469808101654053, -4.115365028381348, 3.5088584423065186); pt.emplace_back(-0.05463826656341553, -3.881274938583374, 3.1585941314697266); pt.emplace_back(0.2328449785709381, -3.614520311355591, 2.850792646408081); pt.emplace_back(-1.0035791397094727, -3.9636778831481934, 3.8052620887756348); pt.emplace_back(-0.7170133590698242, -3.690755844116211, 3.5054771900177); pt.emplace_back(-0.4283676743507385, -3.4318137168884277, 3.1875195503234863); pt.emplace_back(-0.13146033883094788, -3.2284107208251953, 2.7973649501800537); pt.emplace_back(0.15122388303279877, -2.929394483566284, 2.531501293182373); pt.emplace_back(-1.0883535146713257, -3.2573564052581787, 3.513524055480957); pt.emplace_back(-0.7981395721435547, -3.0089566707611084, 3.181861639022827); pt.emplace_back(-0.5119532346725464, -2.733482599258423, 2.88539457321167); pt.emplace_back(-0.21177172660827637, -2.5520896911621094, 2.466628313064575); pt.emplace_back(0.08094996213912964, -2.3205485343933105, 2.113051414489746); pt.emplace_back(-1.1616857051849365, -2.6279513835906982, 3.1217992305755615); pt.emplace_back(-0.8680680990219116, -2.402432918548584, 2.7603931427001953); pt.emplace_back(-0.5798090696334839, -2.140892267227173, 2.4458136558532715); pt.emplace_back(-0.2915937602519989, -1.87905752658844, 2.1316158771514893); pt.emplace_back(0.002897977828979492, -1.6594164371490479, 1.7625699043273926); pt.emplace_back(-1.2396777868270874, -1.967220664024353, 2.7707955837249756); pt.emplace_back(-0.9490185976028442, -1.7218148708343506, 2.43524169921875); pt.emplace_back(-0.6570032835006714, -1.4855259656906128, 2.087836503982544); pt.emplace_back(-0.36362549662590027, -1.2583959102630615, 1.7285257577896118); pt.emplace_back(-0.0764693021774292, -0.9894416332244873, 1.4235832691192627); }
参考
https://sci-pursuit.com/math/statistics/covariance.html
http://www.geisya.or.jp/~mwm48961/statistics/syuseibun1.htm
http://www.it.is.tohoku.ac.jp/lecture/pattern/pdf/2012/slide3.pdf
http://www.nomuraz.com/denpa/prog005.htm
http://sysplan.nams.kyushu-u.ac.jp/gen/edu/Algorithms/PlaneFitting/index.html
次回
固有ベクトルを求めるのだけれど、いかんせんヤコビ法が全く分からないので、その部分はコピペで済ます。