エクセルで作る地形解析図（6）
６．二次元フーリエ変換
画像処理の本を見ていると、たいていフーリエ変換の項目がある。「フーリエ変換は、周期的に変化する波形に対してその特徴をつかむことを目的に利用される手法の一つである」 というような説明がついている。ここまではわかる。「画像解析の分野においては、フーリエ変換は画像のフィルタリング、復元・再構成、符号化などに利用されている」 となってくると、なんだかさっぱりイメージがわかない。にもかかわらず、なぜかやってみた。
地形形状を自然界に見られる複雑な関数とみなし、これにフーリエ変換を適用して正弦波の周波数成分などがどの程度含まれるかを解析することで、その関数の特徴をつかむことができないかなぁと期待してみたんだけど。 理屈やら解説はプロに任せる。ここでは素人なりにエクセルで遊んでみた結果を載せる。 （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において、中心部の値が大きいと周波数の低いパターンが多く存在することを意味し、逆に遠い場所は周波数の高いパターンが存在することを意味する。これを地形に当てはめると、中心部の値が大きいと沢密度の小さい平滑な大斜面の分布が想定され、遠い場所が大きいと細かい沢が発達する地形が想定される。 （と思っていたが、間違いであった。2019.2.15追記）  教科書的には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セル以内のセル値 （右図下の黄緑部分） は無視することにした。当初は半径5セル以内を無視するとしたが、3セル以上であればあまり影響がなさそうなので、方向幅の定義を見直し、2019更新版から無視するのは2セル以内に変更した。 理屈は分からないまま、とにかく以上により変換ボタンをクリックすれば、入力地形画像・パワースペクトル・動径方向分布・角度方向分布が出力され、さらに逆変換ボタンでフーリエ変換出力に対してローパスフィルタ・ハイパスフィルタ・バンドパスフィルタ・ガウシアンフィルタなどを通して逆変換もできるマクロファイルを作成した（6.3）。 6.2 計算例 128×128で計算した例。画像をクリックすると拡大する。ブラウザの 「戻る」 ボタンでこの画面に戻る。
自然対数パワースペクトルは、パワースペクトルを自然対数で表示したもの。 正規化スペクトルは、それを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分かかったがエラーも出ず動いた。
「2D_Fourier_19.zip（4.26MB）」 ： 2019.2.13 ファイル更新
・角度方向スペクトルのレーダーチャートに卓越方向を示すラインを追加。 ※ファイル履歴 2019.2.7 ： 2D_Fourier_19.zip（4.24MB） ・角度方向スペクトル和を計算する際、半径N以内の該当セルを合計すべきを、誤って半径N-1以内となっていたのを修正。 ・方向幅定義を修正。 ・合わせて、「512_seki.zip（7.42MB）」も更新。 ・ほか、個人的興味で微修正。 2014.12.16 ： 2D_Fourier_f.xlsm ・ 任意の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 計算本体の変更はない。半径設定など、私の使用目的に応じてマイナーチェンジ。あと、骨粗鬆症に関する論文にスペクトルの規格化について述べられていたので、それを取り込んでみた。 6.4 地形の二次元フーリエ変換 6.4.1 地形形成に寄与する成分6.4.2 尾根・谷の方向性 6.4.3 尾根と谷とで方向性が異なる事例 6.4.1 地形形成に寄与する成分 地形のフーリエ変換の論文を探しているのだが、最近のものはなく、50mDEMが世の中に出回った頃より以前のものしか見つからない。やはり使えないからなんだろうか。あれこれ、時間が空いた時に思いつくまま試してみるのだが、どうもイメージと合わない。なぜ合わないのか探ってみた。 地形の凹凸を波の合成と考え、二次元フーリエ変換によりその地形を形成するのに最も寄与する成分（波長）を取り出せないかと考えた。フーリエスペクトルを半径ごとに逆変換すると、その成分（波長）が復元される。 128×128セル（ほぼ1.5km四方）で起伏が明瞭な地形の例を示す。図は、エクセルのグラフで等高線表示させたものである。下図は入力した元の地形、右図はフーリエスペクトルの半径1（中心4セル）からの逆変換、半径3以内の逆変換、同様に半径6以内、半径10以内、半径30以内、半径60以内の逆変換結果を示す。 
 二次元だと起伏の差がわかりにくいので、S列とCA列測線の断面を示す。  次に、半径ごと、すなわち波長ごとの逆変換結果を示す。半径1は、直流成分である標高値が含まれていることがわかる。半径6以上は、標高0m付近に細かく振幅するだけなので、ここでは省略する。細かく振幅するということは、元地形の大きな起伏には寄与せず、微地形に反映するということになる。 （ちなみに、直流成分は、標高値の合計 2,878,105.2mをセル数128×128で割った 175.6656mになる。）  どの成分が寄与しているか着目すると、S列ではr3（半径3）が大きいように見える。そこで地形断面との相関係数を求めてみた。S列では確かに半径3が大きくなっているが、CA列のほか、計算してみたBA列、DA列でも、半径2が一番大きい。二次元地形で求めると、同様である。すなわち、局部的には半径3成分に合致する地形もあるということらしい。     尾根筋と断面測線とが一致していると、半径1と地形の相関が最も高くなり、斜交していくと半径2の方が高くなっていくことが想像される。また、尾根や谷の幅が卓越する波長と一致するのかと想像していたが、地形に影響するのは波長ではなく振幅の大きさであり、半径6以上、すなわち1/6波長以上になると振幅は数m以下となり、微地形にしか現れてこないことがわかった。 ということは、地形を形成するのに寄与するのは、半径2〜4の成分の振幅の大きさなんだろうか？ 6.4.2 尾根・谷の方向性 尾根・谷の方向性を示せないものか、試してみた。データは 3×3 のメディアンフィルタによって、尾根を +1 、谷を -1 とした尾根谷図である。ほぼ同じ地質の範囲を取り出し、128×128 セル （ほぼ1.5km四方） でフーリエ変換し、角度方向スペクトル分布を並べてみた。 図の説明をすると、左上は 128×128セルの等高線図（エクセルではなぜこれを等高線図というのだろう？）。右上は条件付き書式で +1 を赤茶、 -1 を水色にセルを塗りつぶしたexcelのシートを縮小表示したもの（作図したものではない）。右下は角度方向スペクトルのレーダーチャートに、卓越すると思われる方向を書き加えたもの。左下はその範囲全体の地形で、左上はそれを分割したものとなっている。 新第三紀堆積岩については、この地域はほぼ南北走向の地層であるが、東西走向の尾根・谷が卓越するらしい。細粒岩の方は極端に谷密度が低いこともあって、地形図を見ても東西走向の尾根・谷は明らかである。粗粒岩の方は、尾根・谷図を見るだけだとあまり方向性がないような気がするが、レーダーチャートと見比べてみると、確かにその方向が卓越しているような気がしてくる。   花崗岩・付加コンプレックス地域は、産総研のシームレス地質図では同一色に塗られているのだが、花崗岩地域は地形を見るだけでも岩相変化が想像される。いずれにせよ、このスケールでは方向性がないことが角度方向スペクトルでも言える。 ※ 花崗岩の右上区域はほぼゴルフ場が占め、人工改変地形。  
6.4.3 尾根と谷の方向性の違い  図は、128×128 セル （ほぼ1.5km四方） の範囲に3×3 のメディアンフィルタを使い、上から 尾根を +1 、谷を -1 とした尾根谷図、尾根だけの図、谷だけの図（水系図）である。 地形を見ると、南北に主尾根が走り、東西に支尾根が延びている。 角度方向スペクトルを見ると、中段の尾根の方は南北と東西の二方向が明瞭に示されている。一方、谷の方は南北方向のピークはあるが大きくない。 左の尾根図では、南北に長く連続する線は3筋（左寄りに2筋、右寄りに1筋）見えるが、水系図の方は2筋である。この違いが出ているように思う。これはメディアンフィルタによる抽出の特徴ではないかと思う。すなわち、地形としては南北方向が優勢であるが、南北方向で抽出されたセルの数よりも、東西方向の沢地形で抽出されたセルの数の方が多いことによるものと思われる。 このことは、他所でも感じる。 ということは、中央に南北の大きな尾根が1筋あるだけの場合、水系図には南北方向は出てこないし、同様に河川が南北にあれば、隣接範囲に尾根が同方向にあっても反映されない。計算する範囲の取り方の問題にもなってくるが、留意しなければならないということだろうか。 |
