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 を使ってみます。

FIR フィルタリングを FFT で高速化してみよう

AlxisDataBlog-画像処理


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






最も簡単な矩形 PSF による 1 次元 Blur の話をします。

\begin{align}
g(y)=\int_{y-1/2}^{y+1/2}f(x)\,dx
\end{align}

という劣化モデルで $g$ から $f$ を求める問題を考えましょう。よくあるアイデアはスペクトル空間で回復させるものです。しかし、この方法には 2 つの欠点が存在します。


欠点1. 計算量の問題
欠点2: アーチファクトの問題


スペクトル空間を用いたデコンボリューションは、PassMark 700 換算で 9M ピクセル毎秒のスループットがあります。FA 分野では生産性向上のため 1G 以上のスループットが要求されます。これを達成するためには、単純計算で PassMark 78000 相当のマシンが必要になります。少し非現実的ですね。


アーチファクトの問題は PSF のスペクトルを計算すれば分ります。

\begin{align}
\int_{-1/2}^{1/2}e^{-2\pi ixy}\,dx&=-\dfrac{e^{-\pi iy}-e^{\pi iy}}{2\pi iy}\\
&=\dfrac{\sin(\pi y)}{\pi y}
\end{align}

零を除くすべての格子点でスペクトルが消えます。これらのスペクトルの零点がアーチファクトになるのですが、実際は最も増幅率の大きい低周波の零点に対応するアーチファクトが支配的に発生します。これは理論と現実の差が最も顕著な部分で、理論的にはアーチファクトは問題になりませんが、現実は上手くいきません。


スペクトルへ変換せずに実空間でデコンボリューションを考えることにしましょう。$g$ を微分すると

\begin{align}
\partial g(x)=f(x+1/2)-f(x-1/2)
\end{align}

したがって、区間 $[-1/2,1/2]$ で初期値 $f_0$ が与えられていれば帰納的に $f$ が計算できます。

\begin{align}
f(x+n+1/2)=f_0(x-1/2)+\sum_{i=0}^n\partial g(x+i)
\end{align}

求めた関数は $f_0$ の積分値が $g(0)$ ならば、必ずデコンボリューションの解になります。デコンボリューションの解の族の中で、どの解が妥当であるのか重みをつけることは $g$ の情報のみからではできません。


初期値 $f_0$ に形をもたせるためには、一旦論理的思考を緩めざるを得ません。例えば $f(x-1/2),f(x+1/2)$ が張る台形で $g(x)$ が近似できると仮定してみましょう。すると

\begin{align}
f(x+1/2)\approx g(x)+\dfrac{\partial g(x)}{2}
\end{align}

という近似が成立します。とりあえず、これを用いて初期値を与えておきましょう。さて、$f$ を計算するときに初期値へ加算する $\partial g$ の級数項は、理論通りであれば $f(x+n+1/2)-f(x-1/2)$ となります。したがって $f(x-1/2)$ を初期値が上手く相殺していなければ、その差分が $f(x+n+1/2)$ にもそのまま乗ってしまう、すなわちアーチファクトになります。


この問題も仮説で回避してみましょう。$\partial g(x)=f(x+1/2)-f(x-1/2)$ の左辺は実験値です。実験値に含まれるノイズにより、この方程式の等号は成立しません。


初期値に含まれる高周波成分はノイズを含む級数項により均されて、充分大きな $n$ に対する $f(x+n+1/2)$ には反映されません。


これで一旦初期値の問題は仮説的に解決しました。アーチファクトが充分減衰した後の校正方法については、また別の機会に紹介したいと思います。(written by J.Nishigata)


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

バーコードの安定符号化技術