エクセルで作る地形解析図(6)

目次参考

6. 2次元フーリエ変換

 画像処理の本を見ていると、たいていフーリエ変換の項目がある。「フーリエ変換は、周期的に変化する波形に対してその特徴をつかむことを目的に利用される手法の一つである」 というような説明がついている。ここまではわかる。「画像解析の分野においては、フーリエ変換は画像のフィルタリング、復元・再構成、符号化などに利用されている」 となってくると、なんだかさっぱりイメージがわかない。にもかかわらず、なぜかやってみた。
 地形形状を自然界に見られる複雑な関数とみなし、これにフーリエ変換を適用して正弦波の周波数成分などがどの程度含まれるかを解析することで、その関数の特徴をつかむことができないかなぁと期待してみたんだけど。
 理屈やら解説はプロに任せる。ここでは素人なりにエクセルで遊んでみた結果を載せる。 (2014.8.7)

6.1 なぜかフーリエ変換

 エクセルには、1次元の高速フーリエ変換機能および逆変換機能がある。これを利用してマクロを組み、2次元高速フーリエ変換および逆変換をできるようにした。エクセルのバージョンは、最初は2010だったが途中から2013になった。まぁ、影響はないだろう。
(2014.10.15追記) Windows版を前提に書いていたが、Mac版エクセルの2008、2011には分析ツールがないそうだ。 6.3項の【注意】参照
 画像処理の分野において、画像の輝度パターンの特徴を数値に変換する方法として2次元フーリエ変換を用いる。ちょっと教科書的になるが、2次元フーリエ変換による出力は空間周波数 (u,v)、およびフーリエスペクトル |F(u,v)| の3次元情報によって表現される。 |F(u,v)|^2 はパワースペクトルと呼ばれ、その3次元情報はパワースペクトル画像 (以下、PSP:Power Spectrum Picture) で表現される。これは一般にはモノクロの濃淡画像としてあらわされるが、せっかくエクセルを使っているのだから、ここではこのグラフ機能を使い、3次元等高線画像として出力した。
 変換後のPSPにおいて、中心部の値が大きいと周波数の低いパターンが多く存在することを意味し、逆に遠い場所は周波数の高いパターンが存在することを意味する。これを地形に当てはめると、中心部の値が大きいと沢密度の小さい平滑な大斜面の分布が想定され、遠い場所が大きいと細かい沢が発達する地形が想定される。
 教科書的にはPSPの特徴ベクトルとして、動径方向分布 p(r) と、角度方向分布 q(θ) を使用して解析するらしい。
 動径方向分布 p(r) は中心から r=sqrt(u^2+v^2) の距離に存在する環状領域内のパワースペクトルの和を表す。右図上は、エクセルで環状領域を 1 、それ以外を 0 とし、条件付き書式で 1 のセルに網かけしたものを縮小して画面コピーし、それに文字等を加筆した図である。つまり、 r を指定することで 1 と 0 にセルを分け、これとパワースペクトルから@SUMPRODUCT関数で和を算出している。
 角度方向分布 q(θ) は水平軸から角度θの線形領域内のパワースペクトルの和を表す。同様に右図下のように角度を指定することでその角度方向和が求められる。なお、角度方向分布 q(θ) において、中心の4セルは45度と135度の2方向となる。4セルのうち右下は周波数成分0 (直流成分) の値であり、地形データの場合は標高値の合計となるため、常に135度方向にピークが出ることになる。ここでは10度ごとの角度方向分布 q(θ) を求めることとした。半径5セル以上であると隣接セルと10度程度以内の角度となるため、半径5セル以内のセル値 (右図下の黄緑部分) は無視することにした。

 理屈は分からないまま、とにかく以上により変換ボタンをクリックすれば、入力地形画像・パワースペクトル・動径方向分布・角度方向分布が出力され、さらに逆変換ボタンでフーリエ変換出力に対してローパスフィルタ・ハイパスフィルタ・バンドパスフィルタ・ガウシアンフィルタなどを通して逆変換もできるマクロファイルを作成した(6.3)。


6.2 計算例

 128×128で計算した例。画像をクリックすると拡大する。ブラウザの 「戻る」 ボタンでこの画面に戻る。

フーリエ変換例 (3例)フーリエ逆変換例


 自然対数パワースペクトルは、パワースペクトルを自然対数で表示したもの。
 正規化スペクトルは、それを256階諧調にしたもの。聞くところによると、フリーウェアのImageJの出力はこれらしい。
 動径方向分布図の半径は、距離(セル数)を自然対数で表示している。
 スペクトル和は、環状領域のスペクトルの和。環状領域は半径が増すとセル数も増えるので、その和を比較するのはおかしくないかと思いセル数で割って1セルあたりの平均を算出してみたが、差はないようだ。(2014.12.16追記: 当初はあまり差はないように思っていたが、いろいろやってみると実はこれがかなり重要なことがわかってきた。理由の一つは半径の取り方。ここではセルの座標をそのまま距離に換算しているが、これだと半径とセル数が比例せず、例えば半径5は36セルだが半径6は32セル、半径10は64セルだが半径11は60セルなど、逆転するところがある。計算方法はこちら。何を目的としているか次第だが、要注意。 )

 なお、角度方向分布図において、32×32サイズだと0度と90度方向がスペクトル和だと小さく出てしまう。これは、マクロファイルの方に記載しておいたが、半径11セル以上じゃないとカウントされないことによる。厳密に考えるなら、これは今後の課題になるのかもしれない。


6.3 マクロファイル

 マクロファイルを作成するに当たり、2次元フーリエ変換の勉強も含めて 友光ほか(2011) と、古堅(2008) を大いに参考にさせていただいた。
 ファイルは自分が使うために作成したので、他人が見てすぐに使い方がわからないかもしれない。一応、説明は入れてみた。
 画像のサイズは、任意の2のn乗セルの正方格子 (8x8、16x16、〜、128x128、〜) に対応させているが、私の環境では普通に動作するのは512x512までで、1024でも動くがちょっと厳しくなる。メモリあるいはクロックの問題と思うが、わからない。とりあえずサンプルとして入れている32x32サイズのデータで動作を試してみるのがいいと思う。128x128でも厳しいという人がいた。私の7年位前のVistaをWin7にアップデートしたExcel 2007では、128x128は10分かかったがエラーも出ず動いた。

 【注意】
このファイルを使ってあなたのPCがおかしくなっても一切関知しない。
マクロが動かないと連絡くれた人のほとんどは、エクセルのアドイン (データ分析 > 分析ツール および 分析ツール VBA) が有効になっていないのが原因だった。ファイルタブ ⇒ オプション でご確認のこと。
このアドインは Mac版のエクセルには入っていないと聞いた。検索してみたら、「Mac OSXのExcelには分析ツールは含まれていないので、 AnalystSoft のWebサイトから StatPlus:mac LE をダウンロードして使用する。」 という記事が出てきた。Macの人はお試しを。
2D_Fourier_f.xlsm(2.02MB)」 : 2014.12.16 ファイル更新
 ・ 任意の2のn乗サイズのデータに対するフーリエ変換と逆変換。
 ・ 動径方向分布と角度方向分布のグラフは 32〜512セルサイズに対応。
   ただし、512サイズの場合は別ファイル「512_seki.xlsm(7.24MB)」を開いておくこと。
 ・ バンドパスフィルタやガウシアンフィルタは 16〜256セルサイズに対応。
   LPF、HPF、BPF、ガウシアン(LPF、HPF、高域強調、LOG)
 ・ BPFオプションを加えた(これは私の個人的興味)。
 ・ このファイルで計算後、変換結果シートだけを別ファイルに移動・保存できる。

※ファイル履歴
2014.10.6 : 2D_Fourier.xlsm
 2種類作成していたファイルを統合。
2014.10.7 : 2D_Fourier_b.xlsm、512_seki.xlsm
 128以下のサイズだと問題なかったが、256以上になるとファイルサイズと処理速度が気になった。このため、不要と思いつつ残していた途中の計算式を削除するようにした。
2014.10.14 : 2D_Fourier_c.xlsm
 動径方向分布図の半径設定を少し改良。
2014.12.16 : 2D_Fourier_f.xlsm
 計算本体の変更はない。半径設定など、私の使用目的に応じてマイナーチェンジ。あと、骨粗鬆症に関する論文にスペクトルの規格化について述べられていたので、それを取り込んでみた。


目次参考