画像の周波数データ変換(DFT)と逆変換(IDFT)
    離散フーリエ変換: 周波数空間データの更なる応用


高速DFTによる画像の周波数データ形式への変換と、その応用の考察 (2017/02/01)

目次
0、はじめに
1、基本的な演算の考え方 (離散フーリエ変換)
2、プログラムのソースコードとWindows用実行形式の提供
3、64bit(マルチスレッド対応)版 dft・idft ソフトを使用する時の注意事項
4、変換処理サンプル
5、データファイルとヘダー情報
6、コマンドの機能と使い方
8、ライセンス及び免責事項
9、終わりに


0、はじめに

古くから画像処理の技術の1つとして利用されている離散フーリエ変換ですが、
計算量が非常に多いので高速化が難しく利用が限られてしまいます。
周波数領域のデータの特徴性を利用する処理が再認識されて来ている様ですので
2次元の画像(空間領域)データに対するDFT・IDFT処理の実装を考えてみます。
DFTの高速処理として、FFT が有りますが入力データの一辺が2のn乗の正方形に
限られてしまいますので、任意の縦横サイズを処理出来る DFT を検討しますが、
FFTに類する高速化アルゴリズムが無いのが難点です。( 画像サイズ:2n x 2m )


1、基本的な演算の考え方 (離散フーリエ変換)

基本原理は、詳しく説明している他のサイトでの確認をお願いします。
 キーワード:フーリエ変換、離散フーリエ変換、オイラーの公式、三角関数の対称性

DFT/IDFT
1次元の DFT・IDFT の考え方とデータ操作方法

DFT/IDFT Y_axis
2次元の画像データに対する DFT・IDFT の考え方(Step1、横1行分の処理を縦方向分繰り返す)

DFT/IDFT X_axis
2次元の画像データに対する DFT・IDFT の考え方(Step2、縦1行分の処理を横方向分繰り返す)


2次元の画像データに対する DFT・IDFT 操作を行う C 言語の ソースコードを掲載。
 画像データは既に、複素数型(COMPLEX)の実部に登録されている事が前提(虚部は,0で初期化済み)
typedef struct COMPLEX { double r, i; } COMPLEX;
#define RotorCOMPLEX(a)   (COMPLEX){ cos(a), sin(a) }
#define AddCOMPLEX(a,b)   (COMPLEX){ (a).r + (b).r, (a).i + (b).i }
#define MulCOMPLEX(a,b)   (COMPLEX){ (a).r*(b).r - (a).i*(b).i, (a).r*(b).i + (a).i*(b).r }
#define ImeDivCOMPLEX(a,b)(COMPLEX){ (a).r / (b), (a).i / (b) }

void dft_TwoDimensional( COMPLEX **d, int lx, int ly, int flag ) {
int x, y, u, v, elements;
double unit;
COMPLEX t[(lx>ly?lx:ly)];

elements = lx; unit = (flag==IDFT) ? 2*M_PI/elements : -2*M_PI/elements;
for( y=0; y<ly; y++ ) { // Σdft_oneDimensional( for Y axis ...);
for( u=0; u<elements; u++ ) t[u]=(COMPLEX){0,0}; // From 0 To LX -1
for( u=0; u<elements; u++ ) {
for( v=0; v<elements; v++ ) t[u] = AddCOMPLEX( t[u], MulCOMPLEX(d[y][v],RotorCOMPLEX(unit*u*v)) );
}
if(flag==IDFT) for( u=0; u<elements; u++ ) d[y][u] = ImeDivCOMPLEX(t[u],elements);
else for( u=0; u<elements; u++ ) d[y][u] = t[u];
} //////////////////////////////////////////////////////////////////////////
elements = ly; unit = (flag==IDFT) ? 2*M_PI/elements : -2*M_PI/elements;
for( x=0; x<lx; x++ ) { // Σdft_oneDimensional( for X axis ...);
for( u=0; u<elements; u++ ) t[u]=(COMPLEX){0,0}; // From 0 To LY -1
for( u=0; u<elements; u++ ) {
for( v=0; v<elements; v++ ) t[u] = AddCOMPLEX( t[u], MulCOMPLEX(d[v][x],RotorCOMPLEX(unit*u*v)) );
}
if(flag==IDFT) for( u=0; u<elements; u++ ) d[u][x] = ImeDivCOMPLEX(t[u],elements);
else for( u=0; u<elements; u++ ) d[u][x] = t[u];
} //////////////////////////////////////////////////////////////////////////
}
補足事項:正規化は逆変換時にまとめて行われるので、この関数で 順変換 逆変換を行う事。
     順変換 逆変換の間で、周波数空間でのフィルター処理を行う場合は、本点を考慮する事


2、プログラムのソースコードとWindows用実行形式の提供

source.zip dft_2D.c
idft_2D.c
complex_array_info.c
showPowerSpectrum.c
dft・idft ソフト及び、複素数型データのファイルヘッダー部確認用ツール、
パワースペクトル表示ツールの各ソースコード
MD5SUM:e661dafdae9633c3aec4f51e6cea25dc
x32.zip dft_2D.exe
idft_2D.exe
complex_array_info.exe
showPowerSpectrum.exe
32bit Windows用 バイナリファイル(実行形式のファイル)
離散フーリエ変換ソフトの構造は非常にシンプルだが、処理時間は長い
MD5SUM:2828fcadbd0761eb28019a6b1ae6ccc8



4、変換処理サンプル

lena 処理対象のBMPファイル
color2gray dft_2D.exe(dft_2Dx64.exe)は、
入力画像を濃淡画像に変換し、
変換後の画像データに対して
DFT 変換を行います。

濃淡画像に変換された画像は、
BMPファイルとして保存されます。
PowerSpectrum dft_2D.exe(dft_2Dx64.exe)は同時に、
フーリエ変換後の周波数データに対し、
そのパワースペクトラム像を作成し、
BMPファイルとして、保存します。
(波形の状態を確認し易いように、
 輝度情報は、log10の対数表示)

また、複素データ形式(2つの連続した
DOUBLE型の領域に実部、虚部を格納)
に変換された情報は、dft_complex.dat
と言う名前で、保存されます。
idft_2D 逆変換により生成された画像
idft_2D.exe(idft_2Dx64.exe)は、
dft_2D.exe(dft_2Dx64.exe)が作成した
ヘダー情報付き複素データ形式の
データを読み込み、逆フーリエ変換
して、結果をBMPファイルとして保存
します。
ヘダーに記載されている情報に
関しては、別項目の”データファイルと
ヘダー情報”を参照してください。
 
 
   
filter dft_2D.exe(dft_2Dx64.exe)が作成した、
周波数変換情報(複素データ形式)に、
フィルター処理を施した、複素数型
データを、idft_2D.exe(idft_2Dx64.exe)
で、逆変換した画像。

"ネコ" ではなくて、"レナ" 画像。
周波数領域での特徴量を更に符号化
すると......
新たな有益な利用方法が、見つかる
かも知れない   記憶、認識とは?
idft_2D
周波数空間で、異なる画像に
属している、周波数情報間の
加算処理をした結果に対して
逆変換を行い再生成した画像





idft_2D
波の独立性より、加算データに
対して、片方の画像に属している
周波数のデータ分を削除処理した
結果に対して、逆変換を行う事で
再生成した画像 説明用(意図的な残部有)

恰も、記憶の一部が削除された
様な結果を生む


1領域(本説明に於いては画面が変換された周波数領域に相当)に記録可能な情報量は2つに制限されない。
総情報量は2以上(公知の事実)であって良い。 lena, Traffic sign, etc
これは、一つの境域に異なる情報を多重化して、追加、削除等の操作が出来る事を意味する。

また、画像による説明は、単なるメタファーに過ぎない。(画像処理に限定されない)



5、データファイルとヘダー情報


typedef struct COMPLEX { double r, i; } COMPLEX;
複素数型のデータ構造は、二つの double 型変数(r,i)から構成されており、
r 部には実部情報が格納され、 i 部には虚部情報が格納される。

struct HEADERCX { // COMPLEX => bfType = 0x585043 (0x43,0x50,0x58)
unsigned int cxType, zero; // identification 0x43='C',0x50='P',0x58='X'
unsigned int cxSize; // ?? cxOffByts+cxSizeData +sizeof(message)
unsigned int cxDataByte; // 16 sizeof(COMPLEX) 固定値
unsigned int cxOffByts; // ?? 32+α sizeof(HEADER) +sizeof(message)
unsigned int cxSizeData; // ?? lx*ly*sizeof(COMPLEX)
signed int cxWidth; // ?? lx
signed int cxHeight; // ?? ly
///////////////////////////////////////////////////////////////////////////
char message[]; // オプション:使用バイト数は cxSize に含む
} HEADERCX ; // cxOffByts の +α の部分
#define COMPLEXID 0x585043 // マジックナンバー(識別子) 0x585043 'CPX'

本ツール群が使用する複素数型のデータは、ファイルの先頭に上記ヘッダー情報を持ち、この後に COMPLEX 型のデータが lx * ly 個続く単純な構造をしている。
ヘッダー部(64バイト)に続けてユーザーの任意の情報を付加する事が可能であるが、追加情報
を使用する場合は、ファイル全体のバイト数を表す cxSize とファイルの先頭から COMPLEX型の実
データが開始される場所までのオフセットを表す cxOffByts が正しく設定されなければならない。
HEADERCX h = { COMPLEXID, // unsigned int 識別子 0x585043 'CPX'
0, // unsigned int 予約領域 0 固定値
HEADERSIZE +message +
lx*abs(ly)*sizeof(COMPLEX), // unsigned int ファイル全体の大きさ
sizeof(COMPLEX), // unsigned int 複素数型のバイト数
HEADERSIZE+message, // unsigned int ファイルの先頭から、
// 複素数型データ列までのバイト数
lx*ly*sizeof(COMPLEX), // unsigned int 複素データ総バイト数
lx, // signed int x軸方向の要素数
ly // signed int y軸方向の要素数
}; ※ 各変数は4バイト unsigned int, signed int
この後に、複素数型のデータが続く。 COMPLEX data[lx*ly]

6、コマンドの機能と使い方

dft_2D.exe, dft_2Dx64.exe
入力されたBMP画像ファイルを読み込み、画像データに対し離散フーリエ変換を行います。
 ・入力されたBMPファイルを濃淡画像に変換(変換画像は、BMPファイルとして保存)
 ・濃淡画像を周波数データに変換(変換データは、datファイルとして保存)
 ・周波数データをパワースペクトラム画像に変換(変換データは、BMPファイルとして保存)

usage: dft_2D.exe [-t] [-h] [input_bmp [complex_out.dat]]

        -t: Show run time.
        -h: This option shows only this message.
 

idft_2D.exe, idft_2Dx64.exe
入力された周波数データファイルを読み込み、波形データに対し逆離散フーリエ変換を行います。
 ・逆離散フーリエ変換された波形データは、画像データとして復元されます。
  復元された画像データは、BMPファイルとして保存されます。
usage: idft_2D.exe [-t] [-h] [complex_in.dat [output_bmp]]

        -t: Show run time.
        -h: This option shows only this message.
 

complex_array_info.exe
入力された周波数データファイルを読み込み、ヘッダーに記録されている情報を表示します
 ・データの縦、横の大きさ等が確認出来ます。

usage: complex_array_info.exe [-h] [complex.dat]

        -h: This option shows only this message.
 

showPowerSpectrum.exe
入力された周波数データファイルを読み込み、パワースペクトラム画像作成します
 ・d オプションを指定すると、低周波成分が四隅に、高周波成分が中心に表示されます
 ・n オプションを指定すると、パワースペクトラムの輝度を log10 変換しません
 ・生成されたパワースペクトラム画像は、BMPファイルとして保存されます。
 ・ソースコードは、独自フィルタープログラムを作成する時の雛形

usage: showPowerSpectrum.exe [-d] [-n] [-h] [complex_in.dat [output_bmp]]

        -d: Not replace the power spectrum.
        -n: normal power spectrum. default display log10.
        -h: This option shows only this message.
 


7、実際の処理は速度的に優れた FFTW 等を利用して下さい。

マルチスレッド化した、64ビットアプリケーションも用意しましたが、FFTW がはるかに優れ
且つ、任意の画像サイズも扱えるので、このアプリケーションの提供は中止しました。
32ビット版のアプリは残してありますが、あまりに遅いので実用性は期待できません。

// fftw を使う場合の参考イメージ
#include "fftw3.h"

void dft_TwoDimensional( COMPLEX **d, int lx, int ly, int flag ) {
fftw_complex *in=(fftw_complex *)&d[0][0], *out=(fftw_complex *)&d[0][0];
fftw_plan    plan;

    if( flag==DFT ) plan=fftw_plan_dft_2d(ly,lx,in,out,FFTW_FORWARD, FFTW_ESTIMATE);
    else            plan=fftw_plan_dft_2d(ly,lx,in,out,FFTW_BACKWARD,FFTW_ESTIMATE);
    fftw_execute( plan );
    fftw_destroy_plan( plan );
}



8、ライセンス及び免責事項

ライセンスは、二条項BSDライセンスに準拠する。

本ソフトウェアは、著作権者およびコントリビューターによって「現状のまま」提供
されており、明示黙示を問わず、商業的な使用可能性、および 特定の目的に対する
適合性に関する暗黙の保証も含め、またそれに限定されない、いかなる保証もない。
著作権者もコントリビューターも、事由のいかんを問わず、損害発生の原因いかんを
問わず、かつ責任の根拠が契約であるか厳格責任であるか(過失その他の)不法行為
であるかを問わず、仮にそのような損害が発生する可能性を知らされていたとしても、
本ソフトウェアの使用によって発生した(代替品または代用サービスの調達、使用の
喪失、データの喪失、利益の喪失、業務の中断も含め、またそれに限定されない)
直接損害、間接損害、偶発的な損害、特別損害、懲罰的損害、または結果損害について、
一切責任を負わないものとする。


9、終わりに

コンピュータの性能が著しく向上しているので、今までは処理時間の問題で難しかった手法も、
利用できる様になった事の意義は大きい。手段の変更はまた、アルゴリズムをより柔軟に考え
る事を可能にすることも意義深い事実である。


  参考資料 (追加情報)

1、慶應大学講義 物理情報数学 C   (strongly recommend)
第一回 信号とシステム
第二回 三角関数と基本周期
信号の表現,重ね合わせの理とたたみこみ積分
第四回 たたみこみ積分の計算法,内積と直交
第五回 フーリエ級数
第六回 複素フーリエ級数からフーリエ変換
第七回 フーリエ変換
第八回 フーリエ変換の性質,ラプラス変換




go to TopPage go to CategoryTop