AlxisDataBlog-プログラミング


MathJax = {
  chtml: {
    matchFontHeight: false
  },
  tex: {
    inlineMath: [['$', '$']]
  }
};






「サクッとコーディング」できて正しく動けばうれしいと思いますが、ちょっとだけ考えてコーディングして「サクッと動く」方がうれしくはないでしょうか。


ちゃんと考えてコーディングしていますか。考えないと頭脳は衰えます。試しに次の式を計算してみましょう。

\begin{align}
\sum_{m=0}^{65535}\Big\{\sum_{n=0}^{65535}\sin(2+3m+5n+7m^2+11n^2)\Big\}
\end{align}

有効桁の話をすると、一気に難しくなるので今回はとりあえず計算できれば良いとしておきましょう。


できましたか。


正解は「37745.31849...」です。答えはさて置き、問題は「どんなコードを書いたか」です。



int i,j;
double temp,sum=0.0;

for(i=0;i<65536;i++){
  for(j=0;j<65536;j++){
    temp=2.0+3.0*i+5.0*j+7.0*i*i+11.0*j*j;
    sum+=sin(temp);
  }
}


このようなコードは、サクッと動かずよくありません。(私のマシンでは 188869ms も掛かりました)


分り易いところだけ重複を回避してみましょう。



int i,j;
double offset,temp,sum=0.0;

for(i=0;i<65536;i++){
  offset=2.0+(3.0+7.0*i)*i;
  for(j=0;j<65536;j++){
    temp=offset+(5.0+11.0*j)*j;
    sum+=sin(temp);
  }
}


セオリー通りですが本質的な改善には程遠いです。もうひとつ、セオリー通りですが、テーブル化という方法もあります。



int i,j;
double temp,sum=0.0;
double *table1=NULL,*table2=NULL;

table1=malloc(sizeof(double)*65536);
table2=malloc(sizeof(double)*65536);
for(i=0;i<65536;i++){
  table1[i]=2.0+(3.0+7.0*i)*i;
  table2[i]=(5.0+11.0*i)*i;
}
for(i=0;i<65536;i++){
  for(j=0;j<65536;j++){
    temp=table1[i]+table2[j];
    sum+=sin(temp);
  }
}
free(table1);
free(table2);
table1=NULL;
table2=NULL;


179712ms。少し速くなりましたが、まだまだです。$\sin$ 関数の引数を再帰的に計算すればどうなるでしょうか。



int i,j;
double offset[4],temp,sum=0.0;

offset[0]=2.0;
offset[1]=10.0;
for(i=0;i<65536;i++){
  offset[2]=0.0;
  offset[3]=16.0;
  for(j=0;j<65536;j++){
    temp=offset[0]+offset[2];
    sum+=sin(temp);
    offset[2]+=offset[3];
    offset[3]+=22.0;
  }
  offset[0]+=offset[1];
  offset[1]+=14.0;
}


乗算がなくなりましたが、これでも私のマシンでは 175249ms も掛かってしまいます。結局のところ入れ子のループ文が処理速度のボトルネックになっていますね。ここで登場するのが「数学」です。

\begin{align}
\sum\sin(2+3m+5n+7m^2+11n^2)&=\sum\{\sin(2+3m+7m^2)\cos(5n+11n^2)+\cos(2+3m+7m^2)\sin(5n+11n^2)\}\\
&=\sum\sin(2+3m+7m^2)\cos(5n+11n^2)+\sum\cos(2+3m+7m^2)\sin(5n+11n^2)\\
&=\sum\sin(2+3m+7m^2)\sum\cos(5n+11n^2)+\sum\cos(2+3m+7m^2)\sum\sin(5n+11n^2)\\
\end{align}

このように、加法定理を使えば変数を分離して、変数ごとの総和計算で値を求めることができます。


これでループ文の入れ子が解消されます。



int i;
double offset1,offset2,temp[4]={0.0},sum;

for(i=0;i<65536;i++){
  offset1=2.0+(3.0+7.0*i)*i;
  offset2=(5.0+11.0*i)*i;
  temp[0]+=sin(offset1);
  temp[1]+=cos(offset1);
  temp[2]+=sin(offset2);
  temp[3]+=cos(offset2);
}
sum=temp[0]*temp[3]+temp[1]*temp[2];


処理速度は、なんと 15ms でした。数学の偉大さを感じるとともに、少し考えて工夫することの重要性を理解いただければ良いです。(written by J.Nishigata)


Copy Right 2020-現在 ALXIS DATA, Inc. All Right Reserved.

コーディングの工数削減よりも重要なこと

AlxisDataBlog-プログラミング


MathJax = {
  chtml: {
    matchFontHeight: false
  },
  tex: {
    inlineMath: [['$', '$']]
  }
};






FIR フィルタリングは画像処理に頻繁に登場します。画像処理系のエンジニアにとっては、FIR フィルタリングの高速化は必須のテクニックと考えても良いでしょう。まずは FIR フィルタを定義しましょう。何でも良いのですが、$20\times 20$ の平滑化フィルタにしましょう。



double BLUR[400];
void defBlur(){
  int i;
  
  for(i=0;i<400;i++){
    BLUR[i]=1.0/400.0;
  }
}


あと、画像ファイルの R/W も定義しておきます。openCV で実装すると簡単です。



IplImage *IMG;
void readImage(
char *filename){
  IMG=cvLoadImage(filename,CV_LOAD_IMAGE_ANYDEPTH | CV_LOAD_IMAGE_ANYCOLOR);
}

void writeImage(
char *filename){
  cvSaveImage(filename,IMG,0);
  cvReleaseImage(&IMG);
}


IplImage 構造体から RGB への R/W も定義します。



double *R,*G,*B;
void splitRGB(){
  int i,j;
  unsigned char value[3];
  
  R=malloc(sizeof(double)*IMG->width*IMG->height);
  G=malloc(sizeof(double)*IMG->width*IMG->height);
  B=malloc(sizeof(double)*IMG->width*IMG->height);
  for(j=0;j<IMG->height;j++){
    for(i=0;i<IMG->width;i++){
      value[0]=(unsigned char)IMG->imageData[3*i+j*IMG->widthStep];
      value[1]=(unsigned char)IMG->imageData[1+3*i+j*IMG->widthStep];
      value[2]=(unsigned char)IMG->imageData[2+3*i+j*IMG->widthStep];
      B[i+j*IMG->width]=(double)value[0];
      G[i+j*IMG->width]=(double)value[1];
      R[i+j*IMG->width]=(double)value[2];
    }
  }
}

void mergeRGB(){
  int i,j;
  
  for(j=0;j<IMG->height;j++){
    for(i=0;i<IMG->width;i++){
      IMG->imageData[3*i+j*IMG->widthStep]=(unsigned char)B[i+j*IMG->width];
      IMG->imageData[1+3*i+j*IMG->widthStep]=(unsigned char)G[i+j*IMG->width];
      IMG->imageData[2+3*i+j*IMG->widthStep]=(unsigned char)R[i+j*IMG->width];
    }
  }
  free(R);
  free(G);
  free(B);
}


高速化の比較対象として、実空間でフィルタリングを実行すると以下のようになります。



void exeFilterOnReal(
int width,
int height,
double *io){
  int i,j,k,l;
  double sum;
  
  for(j=0;j<height-20;j++){
    for(i=0;i<width-20;i++){
      sum=0.0;
      for(l=0;l<20;l++){
        for(k=0;k<20;k++){
          sum+=BLUR[k+20*l]*io[i+k+(j+l)*width];
        }
      }
      io[i+j*width]=sum;
    }
  }
}

int main(){
  clock_t start,stop;
  
  defBlur();
  readImage("input.bmp");
  splitRGB();
  start=clock();
  exeFilterOnReal(IMG->width,IMG->height,R);
  exeFilterOnReal(IMG->width,IMG->height,G);
  exeFilterOnReal(IMG->width,IMG->height,B);
  stop=clock();
  stop-=start;
  mergeRGB();
  writeImage("output.bmp");
  printf("%d\n",stop);
  getch();
  return 0;
}


この実装を FFT で書き直しましょう。FFT ライブラリは FFTW を使ってみます。

高速フーリエ変換、高速化の秘密を知ろう