3元連立方程式をガウスの消去法で解いてみた

タイトルまんまです（ガウスの消去法についてはwiki参照のこと）。逆行列を求めるのに似たような方法を使うのですが、微妙に忘れかけてたので復習を兼ねて解いてみた。本当は精度の問題とかもあって、もっとちゃんとしたアルゴリズムがあるらしいが、そんなのは知らん。あまりにひどいハードコーディングなので、明日はN元まで拡張するつもり。

#include <iostream>
#include <iomanip>

void PrintArray(const double *pArray, int col, int row)
{
    for(int i=0; i<col; ++i)
    {
        for(int j=0; j<row; ++j)
        {
            std::cout << std::setw(10) << pArray[i*row+j] << " ";
        }
        std::cout << std::endl;
    }
    std::cout << "------------------------------" << std::endl;

}

// 入力用データ
// x,y,z = 2, 3, 4
// 2x + 3y + 2z = 21 (1)
// 3x + 4y + 3z = 30 (2)
// 5x + 2y + 4z = 32 (3)
double mat[4*4] = 
{
    2, 3, 2, 21,
    3, 4, 3, 30,
    5, 2, 4, 32,
    1, 1, 1, 1,
};

int main(int argc, char **argv)
{
    //計算用配列
    double tmp[4*4] = {1};

    //まずは(1)式のxの係数を1にする
    for(int i=0; i<4; ++i)
    {
        tmp[i] = mat[i] / mat[0];
    }
    PrintArray(tmp, 4, 4);

    //(2),(3)式のxの係数を１に合わせる
    for(int i=0; i<4; ++i)
    {
        tmp[1*4+i] = mat[1*4+i] / mat[1*4];
        tmp[2*4+i] = mat[2*4+i] / mat[2*4];
    }
    PrintArray(tmp, 4, 4);

    //(1)式から(2),(3)式を引いたものをtmp[2,3][0-3]に格納
    for(int i=0; i<4; ++i)
    {
        tmp[1*4+i] = tmp[0*4+i] - tmp[1*4+i];
        tmp[2*4+i] = tmp[0*4+i] - tmp[2*4+i];
    }
    PrintArray(tmp, 4, 4);

    //次に(2)式のyの係数と(3)式のyの係数を合わせる
    double yk = tmp[1*4+1] / tmp[2*4+1];
    for(int i=1; i<4; ++i)
    {
        tmp[2*4+i] = tmp[2*4+i] * yk;
    }
    PrintArray(tmp, 4, 4);

    //(2)式から(3)式を引く
    for(int i=1; i<4; ++i)
    {
        tmp[2*4+i] = tmp[1*4+i] - tmp[2*4+i];
    }
    PrintArray(tmp, 4, 4);

    //z成分を求める
    double z = tmp[3*4+2] = tmp[2*4+3] / tmp[2*4+2];
    PrintArray(tmp, 4, 4);
    
    //y成分を求める
    double y = tmp[3*4+1] = (tmp[1*4+3] - (z*tmp[1*4+2])) / tmp[1*4+1];
    PrintArray(tmp, 4, 4);

    //xを求める
    double x = tmp[3*4+0] = (tmp[0*4+3] - (z*tmp[0*4+2] + y*tmp[0*4+1])) / tmp[0*4+0];
    PrintArray(tmp, 4, 4);

    std::cout <<"(x,y,z) = (" << x << ", " << y << ", " << z << ")" << std::endl;  

    return 0;
} 

どうでもいいがなんで2次元配列使わなかったんだろう俺…。あと所々でPrintArrayを呼び出してるのはデバッグのためです。あしからず。

linuxからこんにちは

知り合いからいただいたPCを修理することに成功した。電源まわりの調子が悪くて、電源とハードディスクを交換しました。1万弱で新しいPCが手に入ったと思えば安いもんだ。


OSにはlinux Debianを選択。CPUはAthronX64でグラフィックボードはnvidia Quadro4XLG980と6年ほど前の構成。グラボがDirectX8世代なので、そのうちこいつは交換しようと思いますが、それさえやってしまえばグラフィックプログラミングには困らないスペックになりそうです。

この投稿もlinuxからのものです。ブラウザさえ入っていれば日常の用途にはまったく困らない。いい時代ですね。

これから開発環境を整えようと思います。eclipseも考えたんですが、debianだとちょっと面倒らしいので止めました。なので、ガチのemacsで整えちゃおうかなと。基本的なキーバインドすら覚えていないですが（笑）

OpenGLでテクスチャ表示

楽勝だと思っていたが意外と苦労した。アラインメントや配列について学びなおすいい機会だったと思う。とりあえず、テクスチャを表示するだけコードを載せてみよう。

#include <gl/glut.h>

//頂点座標
float vertices[][2] =
{
   {-0.5f,  0.5f},
   {-0.5f, -0.5f},
   {0.5f, -0.5f},
   {0.5f,  0.5f},
};

//uv(st)座標
float uvs[][2] =
{
   {0.f, 1.f},
   {0.f, 0.f},
   {1.f, 0.f},
   {1.f, 1.f},
};

//テクスチャハンドル
GLuint hTextures[1];
const unsigned TEXSIZE = 64;

//画素バッファ
GLubyte texImage[TEXSIZE][TEXSIZE][4];

//アイドル時にも画面は更新する
void Idle()
{
   glutPostRedisplay();
}

//描画関数
void Display()
{
   glClear(GL_COLOR_BUFFER_BIT);
   glEnable(GL_TEXTURE_2D);
   glBindTexture(GL_TEXTURE_2D, hTextures[0]);

   //反時計回りに板ポリを描画
   glBegin(GL_QUADS);
   {
      glVertex2fv(vertices[0]); glTexCoord2fv(uvs[0]);
      glVertex2fv(vertices[1]); glTexCoord2fv(uvs[1]);
      glVertex2fv(vertices[2]); glTexCoord2fv(uvs[2]);
      glVertex2fv(vertices[3]); glTexCoord2fv(uvs[3]);
   }

   glEnd();
   glBindTexture(GL_TEXTURE_2D, 0);

   glDisable(GL_TEXTURE_2D);
   glutSwapBuffers();
}

//テクスチャ初期化
void InitTexture()
{
   //テクスチャを１枚生成する
   glGenTextures(1, hTextures);

   //今回は計算でテクスチャを作る
   for(unsigned y=0; y<TEXSIZE; ++y)
   {
      int yy = y / 8;
      for(unsigned x=0; x<TEXSIZE; ++x)
      {
         int xx = x / 8;
         int col = ((yy+xx)&1)*255;

         texImage[y][x][0] = col; //r
         texImage[y][x][1] = 0;   //g
         texImage[y][x][2] = ~col;//b
         texImage[y][x][3] = 0xff;//a
      }
   }

   glPixelStorei(GL_UNPACK_ALIGNMENT, 4);

   //hTexures[0]の設定をここから行うという宣言
   glBindTexture(GL_TEXTURE_2D, hTextures[0]);

   //VRAMにイメージを転送
   glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, TEXSIZE, TEXSIZE, 0, GL_RGBA, GL_UNSIGNED_BYTE, texImage);

   //テクスチャマッピング方式。今回は繰り返しなし
   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);

   //補完フィルタ設定。市松模様なのでnearestでおｋ
   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);

   //ポリゴンカラーとの混色設定
   glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);

   //hTexures[0]の設定終了
   glBindTexture(GL_TEXTURE_2D, 0);
}

//エントリポイント
int main(int argc, char **argv)
{
   //お決まりの処理
   glutInit(&argc, argv);
   glutInitDisplayMode(GLUT_DOUBLE| GLUT_RGBA);
   glutCreateWindow("DrawTexture");
   glutDisplayFunc(Display);
   glutIdleFunc(Idle);

   glClearColor(0,0,0,1);
   InitTexture();

   //メインループに入る
   glutMainLoop();

   //終了処理
   glDeleteTextures(1, hTextures);
   return 0;
}

テクスチャを表示するだけでこの面倒さ。実際にはＢＭＰなどの画像を読み込む処理も要りますから、これからが大変ですな・・・。あと、3次元配列を最初から用意しておくのもどうかと思うので、こいつをクラス化したりもしました。この辺はまたあとで。

VC++でアラインメントを指定する

ビットマップ読み込み関連で失敗したのでメモ。

struct Hoge
{
    char c;
    int i;
};

上記の構造体サイズは理屈上は5バイトであるが、パディングされてるので実際には8バイトになっている。こいつを確実に5バイトにしたいときは#pragma pack命令を使う。

#pragma pack(push, 1)
struct Hoge
{
    char c;
    int i;
};
#pragma pack(pop)

1バイト単位でアラインメントを調整するので、結果としてパディングされないということになる。

サンプルコード

#include <stdio.h>

struct A
{
    char c;
    int i;
};

#pragma pack(push,1)
struct B
{
    char c;
    int i;
};
#pragma pack(pop)

int main()
{
    printf("sizeof A: %d\n", sizeof(A));// >> sizeof A: 8
    printf("sizeof B: %d\n", sizeof(B));// >> sizeof B: 5
    return 0;
}

バイナリでファイルを読み込むときに、使えそうなテクだ。
参考サイト

• Visual C++でアラインメントを扱う3つの機能、#pragma pack, __declspec(align(#)), __alignof演算子についてのメモ

OpenGL+GLSLでポイントライトと距離減衰

距離によって強さが変わるライトの実装に成功した。

ある地点における光の強さは、光源までの距離を d 、一定減衰率を kc 、1時減衰率を k1 、2時減衰率を k2 とすると以下の式で表すことができます。

減衰率kc, k1, k2はそれぞれGL側（ホストプログラム）からシェーダ側に送ることができる。

以下はホスト側のプログラム

//距離減衰
//kc: 一定減衰率 = 0
//k1: 1次減衰率 = 0
//k2: 2次減衰率 = 1.0 / (lightPos.y * lightPos.y)
glLightf(GL_LIGHT0, GL_CONSTANT_ATTENUATION, kc);
glLightf(GL_LIGHT0, GL_LINEAR_ATTENUATION, k1);
glLightf(GL_LIGHT0, GL_QUADRATIC_ATTENUATION, k2);

んで、フラグメントシェーダ側の記述は以下のとおり

varying vec3 P;
varying vec3 N;

void main()
{
    //法線やライティング情報を受け取る
    vec3 L = gl_LightSource[0].position.xyz - P;

    //ライトまでの距離
    float d = length(L);

    //ライトへの向きベクトルを正規化する
    L = normalize(L);

    //減衰係数
    float attenuation = 1.0 / (gl_LightSource[0].constantAttenuation
                       + gl_LightSource[0].linearAttenuation * d
                       + gl_LightSource[0].quadraticAttenuation * d * d);

    vec3 N = normalize(N);

    //アンビエント確定
    vec4 ambient = gl_FrontLightProduct[0].ambient;
    float dotNL = dot(N, L);

    //ディフューズ確定
    vec4 diffuse = gl_FrontLightProduct[0].diffuse * max(0.0, dotNL);

    vec3 V = normalize(-P);
    vec3 H = normalize(L + V);
    float powNH = pow(max(dot(N, H), 0.0), gl_FrontMaterial.shininess);

    if(dotNL <= 0.0)
    {
        powNH = 0.0;
    }
    //スペキュラ確定
    vec4 specular = gl_FrontLightProduct[0].specular * powNH;

    //カラー値統合
    gl_FragColor = ambient + diffuse + specular;

    //減衰係数を乗じる
    gl_FragColor *= attenuation;
}

まあほとんど教本のまんまなんですが・・・。

gl_LightSource[0].constantAttenuation
gl_LightSource[0].linearAttenuation
gl_LightSource[0].quadraticAttenuation

がそれぞれホスト側で設定した減衰係数です。

画像の式に当たるのがコメントで減衰係数と書いてあるところ。んで、最後にカラー値を統合したgl_FragColorに減衰係数をかけて明るさを確定してるわけですね。このあたりはレイトレースでやったことがかなり活きてるかんじ。

OpenGLで簡易シャドウ

OpenGLでのお話。今まで敬遠してきたシャドウイングにチャレンジ。簡易シャドウというらしい。原理はまだよくわからないが、使用している行列や平面のパラメトリック表現から察するに、ライトの位置から平面にびたーっと張り付くように頂点変形をさせるのが影の変換行列ってやつらしい。

影の変換行列
平面の方程式を ax + by + cz + d = 0 、光源ベクトルをe(e_x, e_y, e_z) とすると、影の変換行列(Ms)は

となります。行列の意味についていろいろ考えたいが、とりあえず後回し。

しかし、よくよく見るとポリゴンが重なってる部分は影が濃く出てますね。自分の影を見ればわかりますが、透明なガラスならともかく、この影は不自然ですね。おそらくこれが「簡易」シャドウといわれる所以でしょう。

上記のマトリクスを最初はHTMLに書き込もうと思ったが、画像化してくれる超便利なツールがオンラインで使用できるのでそちらを使った。LATEX形式だが、知らなくても適当にボタンを押してるだけでそれっぽく仕上がるナイスな出来！

VC++のキーワードハイライト機能を拡張する

VC++2008のお話。

"C:\Program Files\Microsoft Visual Studio 9.0\Common7\IDE"フォルダ下にある"usertype.dat"を開く。なければ作成。

登録したいキーワードを追加する。記述の仕方はいたって簡単で、ハイライトしたいキーワードを１行ずつ追加するだけ。以下の例はGLSLの組み込み型。

次にVCのメニューバーから"ツール"->"オプション"->"テキストエディタ"->"ファイル拡張子"を選択。追加したい拡張子と、編集用のエディタを選択する。

私の場合はバーテクスシェーダファイルは"vert"、ピクセルシェーダファイルは"frag"で登録しました。

最後にVCを再起動して終了。

これだけだとほかのブログにもあるので、さらに小ネタを追加。先ほど追加した"usertype.dat"ですが、新しい型が出てくるたびにフォルダから探して編集するのが面倒なので、ショートカットを作成しちゃいましょう。

最後にショートカットの"プロパティ"->"リンク先"をテキストエディタに変更して、もともとのリンク先はそのまま引数として使う。

画像の1番が追加する使用したいエディタへのパス。2番がもともと記述されていた"usertype.dat"へのパスです。ショートカットのアイコンがエディタのものに変わっていれば成功。このショートカットをデスクトップや、プロジェクトのディレクトリに放り込んでおけばいつでも1クリックで追加したい型を編集可能になります。

OpenSSH と Squidでプロキシサーバー

学校のwebフィルターがあまりにも鬱陶しかったので家に鯖を立てて串を通すことにした。
Linuxはdebian, opensshとsquidを入れて、学校のwinからputtyというソフトを使ってキャッシュをコピーしてくる。
ほとんど担任の入れ知恵だが、これで調べ物が楽にできる。


　　　　　　　 ／￣￣￣＼
　 　 　 　 ／　　　　　　　 ＼
　　　　　/　　 ─　　　─　　ヽ
　　　 　 |　　 （●） 　（●）　　|
　　　　　＼　　 （__人__）　__,／
　　　　　／　　　｀ ⌒´　　　＼
　　　＿/((┃))＿＿＿＿_＿i　|　ｷｭｯｷｭｯ
..　／ /ヽ,,⌒) ￣￣￣￣￣（,,ノ ＼
／　 /＿＿＿＿＿＿＿＿＿ヽ.. 　＼
.￣￣￣￣￣￣￣￣￣￣￣￣￣￣
　　　　 　　　　＿＿＿_
　　　　 　 　／⌒ 　 ー ＼
　 　 　 　／ （●） 　（●） ＼　　＋
　　　 　/ :::::⌒（__人__）⌒:::::ヽ
　 　 　 |　　　　 |r┬-|　　　　|　　+
.　 　 　 ＼_　　　`ー'´ 　 _,／
　 　 　　／　 　 　 　 　 　 ＼　　　　　+
　 　 　 |￣￣￣￣￣￣￣￣ |　　ﾄﾝ
　　　_（,,) 　　 自　由　 　　 (,,）_
　 ／ 　|　　 　　　　　　　 　 　|　 ＼
／　 　 |＿＿＿＿＿＿＿＿_|　　　＼

それにしてもクリエイティブ系の学科なのにゲーム会社のサイトすら見れないってどういうことよ・・・。

PythonでFizzBuzz

4月→新入社員→研修→FizzBuzz!!

ということでやってみた。社会人じゃないけど。
一番苦労したのがprintだと改行が上手くいかないので最後にカンマをつけるってところ。

※下記のプログラムには誤りがあります！

for i in range(1, 100):
    if(i % 3 == 0) or (i % 5 == 0):
        print 'fizz',
        if(i % 5 == 0):
            print 'buzz',
    else:
        print i,
    print '\n',

この手のアルゴリズムネタはスクリプト言語でやるのが楽しいね。

（追記）
とか格好つけて嘘っぱちなコードをネット上に晒すとかアホもいいとこだ。改めて実行すると上記のアルゴリズムじゃダメダメですね。即効で直したものを上げておこう。あー恥ずかしい！

for i in range(1, 100):
    if(i % 3 == 0) or (i % 5 == 0):
        if(i % 3 == 0):
            print 'fizz',
        if(i % 5 == 0):
            print 'buzz',
    else:
        print i,
    print '\n',

参考サイト

• とりあ えず: [Python][お勉強] Python入門(16) - printステートメント

VC++のビルド時に作られた中間ファイルを削除する

久々にプログラマらしいお仕事をしたので書き残そう。
詳細はあとで編集。とりあえず体裁整える前のコードだけ上げておく。

# -*- coding: mbcs -*-
from fnmatch import fnmatch
import os

root = os.getcwd()
print root

for dirs in os.walk(root):
    path = dirs[0] + '\\'
    for f in dirs[2]:
        if fnmatch(f, '*.obj') or fnmatch(f, '*.pch'):
            print f
            os.remove(path+f)

FBX SDKのメモリリークに対応する

Autodesk FBX SDK 2010.2を使っていると単純なプログラムでもメモリリークを起こすことがわかった。
メモリリークを起こす例

//マネージャの生成
mySdkManager = KFbxSdkManager::Create();

//マネージャの破棄（多分ここでリークしている）
mySdkManager->Destroy();

こいつに対応するにはIOSREF.FreeIOSettings()という命令をマネージャを破棄する直前に追加すればよい。

//マネージャの生成
mySdkManager = KFbxSdkManager::Create();

//マネージャの破棄
IOSREF.FreeIOSettings();//<-追加
mySdkManager->Destroy();

参考リンク

• Memory Leaks in 2010.0
• FBX SDKでのメモリリーク