トップ | puarts.com
ロゴ
「C言語」に関連する記事一覧
0  

4 x 4の逆行列の計算面倒だなぁと思って検索したら、下記ページに貼ってあるgluの逆行列計算関数がそのまま使えたので転載しておきます。適当な視点行列を逆変換してみて、正しく動作することも確認できました。

http://stackoverflow.com/questions/1148309/inverting-a-4x4-matrix


bool gluInvertMatrix(const double m[16], double invOut[16])
{
    double inv[16], det;
    int i;

    inv[0] = m[5]  * m[10] * m[15] - 
             m[5]  * m[11] * m[14] - 
             m[9]  * m[6]  * m[15] + 
             m[9]  * m[7]  * m[14] +
             m[13] * m[6]  * m[11] - 
             m[13] * m[7]  * m[10];

    inv[4] = -m[4]  * m[10] * m[15] + 
              m[4]  * m[11] * m[14] + 
              m[8]  * m[6]  * m[15] - 
              m[8]  * m[7]  * m[14] - 
              m[12] * m[6]  * m[11] + 
              m[12] * m[7]  * m[10];

    inv[8] = m[4]  * m[9] * m[15] - 
             m[4]  * m[11] * m[13] - 
             m[8]  * m[5] * m[15] + 
             m[8]  * m[7] * m[13] + 
             m[12] * m[5] * m[11] - 
             m[12] * m[7] * m[9];

    inv[12] = -m[4]  * m[9] * m[14] + 
               m[4]  * m[10] * m[13] +
               m[8]  * m[5] * m[14] - 
               m[8]  * m[6] * m[13] - 
               m[12] * m[5] * m[10] + 
               m[12] * m[6] * m[9];

    inv[1] = -m[1]  * m[10] * m[15] + 
              m[1]  * m[11] * m[14] + 
              m[9]  * m[2] * m[15] - 
              m[9]  * m[3] * m[14] - 
              m[13] * m[2] * m[11] + 
              m[13] * m[3] * m[10];

    inv[5] = m[0]  * m[10] * m[15] - 
             m[0]  * m[11] * m[14] - 
             m[8]  * m[2] * m[15] + 
             m[8]  * m[3] * m[14] + 
             m[12] * m[2] * m[11] - 
             m[12] * m[3] * m[10];

    inv[9] = -m[0]  * m[9] * m[15] + 
              m[0]  * m[11] * m[13] + 
              m[8]  * m[1] * m[15] - 
              m[8]  * m[3] * m[13] - 
              m[12] * m[1] * m[11] + 
              m[12] * m[3] * m[9];

    inv[13] = m[0]  * m[9] * m[14] - 
              m[0]  * m[10] * m[13] - 
              m[8]  * m[1] * m[14] + 
              m[8]  * m[2] * m[13] + 
              m[12] * m[1] * m[10] - 
              m[12] * m[2] * m[9];

    inv[2] = m[1]  * m[6] * m[15] - 
             m[1]  * m[7] * m[14] - 
             m[5]  * m[2] * m[15] + 
             m[5]  * m[3] * m[14] + 
             m[13] * m[2] * m[7] - 
             m[13] * m[3] * m[6];

    inv[6] = -m[0]  * m[6] * m[15] + 
              m[0]  * m[7] * m[14] + 
              m[4]  * m[2] * m[15] - 
              m[4]  * m[3] * m[14] - 
              m[12] * m[2] * m[7] + 
              m[12] * m[3] * m[6];

    inv[10] = m[0]  * m[5] * m[15] - 
              m[0]  * m[7] * m[13] - 
              m[4]  * m[1] * m[15] + 
              m[4]  * m[3] * m[13] + 
              m[12] * m[1] * m[7] - 
              m[12] * m[3] * m[5];

    inv[14] = -m[0]  * m[5] * m[14] + 
               m[0]  * m[6] * m[13] + 
               m[4]  * m[1] * m[14] - 
               m[4]  * m[2] * m[13] - 
               m[12] * m[1] * m[6] + 
               m[12] * m[2] * m[5];

    inv[3] = -m[1] * m[6] * m[11] + 
              m[1] * m[7] * m[10] + 
              m[5] * m[2] * m[11] - 
              m[5] * m[3] * m[10] - 
              m[9] * m[2] * m[7] + 
              m[9] * m[3] * m[6];

    inv[7] = m[0] * m[6] * m[11] - 
             m[0] * m[7] * m[10] - 
             m[4] * m[2] * m[11] + 
             m[4] * m[3] * m[10] + 
             m[8] * m[2] * m[7] - 
             m[8] * m[3] * m[6];

    inv[11] = -m[0] * m[5] * m[11] + 
               m[0] * m[7] * m[9] + 
               m[4] * m[1] * m[11] - 
               m[4] * m[3] * m[9] - 
               m[8] * m[1] * m[7] + 
               m[8] * m[3] * m[5];

    inv[15] = m[0] * m[5] * m[10] - 
              m[0] * m[6] * m[9] - 
              m[4] * m[1] * m[10] + 
              m[4] * m[2] * m[9] + 
              m[8] * m[1] * m[6] - 
              m[8] * m[2] * m[5];

    det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];

    if (det == 0)
        return false;

    det = 1.0 / det;

    for (i = 0; i < 16; i++)
        invOut[i] = inv[i] * det;

    return true;
}

大体検索すればコードが出てくるので、大学時代に行列演算をやったのが全く生きてこないですね。笑

自分でチューニングをするというのは逆にプログラムを遅くしてしまう場合があるということを検証してみたかったので、コンパイラで最適化しずらそうなコードと最適化しやすそうなコードで最適化前と後で速度の変化はどのようになるかを実験してみました。

実験環境は以下です。

OS: Windows 7
CPU: core i7 860
メモリ: 8GB

次のようなコードでコンパイラによる最適化の前と後でfunc1とfunc2それぞれの速度を比較しました。

#include <stdlib.h>

void func1(int *data, int len)
{
    int i, j;
    for(i=0; i<len; i++){
        data[i]=0;
        for(j=0; j<i; j++){
            data[i]+=j;
        }
    }
}
 
void func2(int *data, int len)
{
    int i, j;
    int *ptr;
    for(i=0; i<len; i++){
        ptr = &data[i];
        *ptr = 0;
        for(j=0; j<i; j++){
            *ptr += j;
        }
    }
}

void main(){
    int len = 100000;
    int *data = (int*)malloc(len*sizeof(int));
 
    func1(data, len);
    func2(data, len);
    free(data);
}

このコードのままであれば、func2の方が余計なポインタ演算が減るので高速になるはずです。下が実験の結果で、最適化を行わないとfunc2の方が高速になるという予想通りの結果になりました。

最適化なしの結果
func1: 18.88 secs
func2: 15.65 secs

では次に最適化オプションをつけて速度を比較してみます。下が実験結果で、最適化前とうってかわってfunc1の方が圧倒的に速くなりました。

Visual C++の実行速度の最大化 (/O2)オプションで最適化したときの結果
func1: 3.13 secs
func2: 9.26 secs

実験結果から、コンパイラの挙動をよく知らない上でチューニングを行うことは逆にコードを複雑にしてしまい、コンパイラの自動的な最適化を妨げる原因になり得て、プログラムの動作速度低下を招く危険をはらんでいるということがわかりました。

どこでどのようなチューニングを行えば速度が向上するのかは当然コンパイラに依存するところなので、そこまで考慮するのはとても高度な知識が必要になるかと思います。私にはそこまで考えられる技量がないので、試しに速くなるかやってみて速くなれば採用するという風にしています。

ただ、コンパイラを変えたら遅くなるということも当然あり得ることなので、シンプルなコードとチューニングしたコードは両方とも残しておくようにして、いつでも切り替えられるようにしておくのがベストなのかなと思います。
OpenMPで並列化によりずいぶんと処理の重い画像処理系のC言語で書いたプログラムの速度を高めることができたのですが、逆にループ回数が多くても一回の処理が軽いときはOpenMPを使うと逆に遅くなってしまうような場合もありました。
 
OpenMPを使っていてなんとなく分かってきたのが、OpenMPによる並列化は出来るだけループの外側に持ってくるのが良くて、さらに高速化をするためには内側にある細かいループを出来る限り自分でチューニングする必要があるということです。
 
ということで、今使っている処理速度をほとんど考えていないコードに簡単なチューニングをしてみたところ結構速くなったので、どんなチューニングを行ったのかを記しておこうと思います。
 
まずはコードのクリティカルな部分を探さなければならないのですが、画像処理でよく出てくる処理は何かと探すと、そのひとつに画像データのコピーがありました。なので、画像のコピーを例にチューニング例を書きたいと思います。
 
具体的にチューニングする処理内容なのですが、画像を読み込むときにpng画像をunsigned charバイト列で読み込み、それを少数演算用にfloatのデータ配列にコピーするという処理が結構あったので、そのコピーの部分に着目します。
 
まず元となる何のチューニングもないコード例は次のようなコードとします。dst_dataがコピー先float型1次元配列、src_dataがコピー元unsigned char型1次元配列で、widthが横解像度、heightが縦解像度です。

 
int x, y, i, j;
int dst_index, src_index;
int num_channels_src = 4; //コピー元のチャンネル数
int num_channels_dst = 3; //コピー先のチャンネル数
for(y=0; y<height; y++){ //heightは縦画素数
    for(x=0; x<width; y++){ //widthは横画素数
        i = y*width + x;
        src_index = i * num_channels_src;
        dst_index = i * num_channels_dst;
        for(j=0; j<num_channels; j++){
            dst_data[dst_index+j] = (float)src_data[src_index+j];
        }
    }
}


画像処理では縦横の解像度の他にRGBのようなチャンネル数を考えなければならないので、処理速度に関して何も考えていない場合は3次元的なループを書く人もいると思います。ちなみにコード中のnum_channels_srcはコピー元データのチャンネル数で、コピー元unsigned char型配列のデータsrc_dataはアルファチャンネルを含んだRGBA構造なので4チャンネル、そしてnum_channels_dstはコピー先のチャンネル数で、コピー先のfloat型配列dst_dataはアルファチャンネルを省いたRGBの3チャンネルのみを格納するという場合を考えています。
 
ループ回数が多いほど少しの無駄の省きで速度が上がるので、結果がわかりやすいように10000x10000の画像を実験に使いました。30回処理速度を計測して平均を取りました。チューニング前の平均速度は下記のようになりました。
 
30回の平均速度: 3.454667 sec
 
ということでまず最初のチューニングを行います。まずはかなり初歩的なところですが、画像の縦横の2重ループを1重ループに直します。
 
int i, j;
int dst_index, src_index;
int num_channels_src = 4;
int num_channels_dst = 3;
int size = width*height;
for(i=0; i<size; i++){ //sizeは画素数
    src_index = i * num_channels_src;
    dst_index = i * num_channels_dst;
    for(j=0; j<num_channels; j++){
        dst_data[dst_index+j] = (float)src_data[src_index+j];
    }
}

30回の平均速度: 3.269134 sec
 
若干処理速度が改善されました。
さらに、このコードではsrc_index, dst_indexの演算も無駄です。続いてそれも省きます。
 
int i, j;
int num_channels_src = 4;
int num_channels_dst = 3;
int nc_diff = num_channels_src - num_channels_dst;
int size = width*height;
float *dst_ptr = dst_data;
uchar *src_ptr = src_data;
for(i=0; i<size; i++){
    for(j=0; j<num_channels_dst; j++){
        *dst_ptr = (float)(*src_ptr);
        dst_ptr++;
        src_ptr++;
    }
    src_ptr+=nc_diff;
}


30回の平均速度: 2.635134 sec
 
これだけでずいぶん処理速度が向上しました。次にポインタをカウンタ代わりに使えばカウンタiが不要になるので、カウンタiを省きます。
 
int num_channels_src = 4;
int num_channels_dst = 3;
int nc_diff = num_channels_src - num_channels_dst;
int size = width*height;
float *dst_ptr = dst_data;
uchar *src_ptr = src_data;
float *end_ptr = dst_data + size*num_channels_dst;
int cnt=0;
for(; dst_ptr!=end_ptr; dst_ptr++, src_ptr++){
    *dst_ptr = (float)(*src_ptr);
    cnt++;
    if(cnt==o_num_channels){
        src_ptr+=nc_diff;
        cnt=0;
    }
}

 
30回の平均速度: 2.558067 sec
 
これだけでも目に見えて処理時間が縮まりました。私の知識ではチューニングはここまでが限界ですので、画像データコピーのチューニングは以上になります。
最終的なコードはずいぶんみづらくなってしまいましたが、ちゃんとわかりやすいコメントを入れておけば問題はないと思います。
今回は10000x10000の画像で0.9秒高速になりましたが、これをパーセンテージで言えば処理時間は74パーセントほどに縮まったことになります。これはかなり大きいことだと思います。
 
しかし、何でもかんでもチューニングすれば良いという話ではなく、コードの可読性と効率の良さのバランスを考えることがとても大切だといろいろなプログラミングの参考書に書いてあったのを覚えています。この両者のバランスというのは多くの人と一緒に開発するような経験がないと判断できないものだと思うので、私にはまだよくわかりませんが、いつかは最適なコードが書けるようになりたいです。
私がC言語でエラーメッセージを吐き出すときにいつも使っている関数を紹介します。
エラーが起こったコードのファイル名、行番号と任意のメッセージを出力可能なので一回作っておくと結構便利です。さらにエラーメッセージに決まりがある場合はこの関数を自分用にカスタマイズすれば、もっと便利な関数になります。

使用例を含めたソースコードを以下に記しておきますので、エラー出力をfprintfなどのみで吐き出していた人は是非コピペして活用して下さい。

Sample Source Code

#include <stdio.h>
#include <stdarg.h>

inline void MyError(char *file, int line, char *fmt, ...){
    va_list ap;
    va_start(ap, fmt);
    fprintf(stderr, "%s: %d ", file, line);
    vfprintf(stderr, fmt, ap);
    va_end(ap);
}
#define MYERROR(fmt, ...) MyError(__FILE__, __LINE__, fmt, __VA_ARGS__)

void main(){
    int code=-1;
    MYERROR("error code is %d\n", code);
}


Result
main.cpp: 15 error code is -1
BVHの階層構造をファイルに書き込むときにC言語でインデントの出力が必要でした。

そのとき使った方法をメモしておきます。


#include <stdio.h>
int main(){ char fmt[128]; sprintf(fmt, "%% %ds",5); printf(fmt,""); printf("INDENT1\n"); sprintf(fmt, "%% %ds",10); printf(fmt,""); printf("INDENT2"); return 0; } // Result // INDENT1 // INDENT2
C言語でフーリエ変換した画像を周波数成分の振幅値で除算して位相部分のみ抜き出して画像に出力するプログラムを作りました。FFTにはFFTW、画像の入出力にはOpenCVを使っています。
下の画像が入力と出力です。

I made a C program for extration of phase channel from gray scale image. In this progaram, I used FFTW for FFT, and OpenCV for image I/O.

Here is an input image.
inputImg.jpg

Here is a result image.
4a375fbe.jpeg

ソースコードも置いておきます。

Here is source code
phaseExtraction.c

デバッグ用、リリース用のライブラリをコード上から指定する方法を昔書いた勉強用ノートを読み返していてそういえばこんなのあったなと思い出しました。

Visual Studioでは環境が変わるごとにライブラリの指定をしなおさなくて済むのがとても便利です。また忘れぬようブログにもメモしておきます。
 

OpenCVのライブラリを指定する例

#ifdef _DEBUG
#pragma comment(lib, "C:/OpenCV210/lib/cxcore210d.lib")
#pragma comment(lib, "C:/OpenCV210/lib/cvaux210d.lib")
#pragma comment(lib, "C:/OpenCV210/lib/cv210d.lib")
#pragma comment(lib, "C:/OpenCV210/lib/highgui210d.lib")
#else
#pragma comment(lib, "C:/OpenCV210/lib/cxcore210.lib")
#pragma comment(lib, "C:/OpenCV210/lib/cvaux210.lib")
#pragma comment(lib, "C:/OpenCV210/lib/cv210.lib")
#pragma comment(lib, "C:/OpenCV210/lib/highgui210.lib")
#endif

intを0詰めの任意の桁数文字列に変換するという何気ないことですが、一瞬迷ってしまったのでメモを残しておきます。


#include <string.h>
#include <stdio.h>
int main(){ char num[64],format[64]; sprintf(format, "%%0%dd", 5); sprintf(num, format, 12); printf(num); return 0; } // Output: 00012

C言語で DFT のプログラムを書いてみました。

0  

0.0234 sec
にほんブログ村 ゲームブログ ファイアーエムブレムへ にほんブログ村 デザインブログ コンピュータグラフィックスへ

Copyright(C)2006-2018 wsp All Rights Reserved