時系列画像の相互相関を利用した 移動体の軌道

平成 16 年度
卒業論文
時系列画像の相互相関を利用した
移動体の軌道ベクトルの計測：基礎理論と
気象衛星データへの応用
Trajectory vector estimation of moving objects using
cross-correlation of multi-temporal images: principle
and application to meteorological satellite data
指導教員
大内和夫教授
高知工科大学
工学部
物質環境システム工学科
1050026
柏井聡
1050042
新名宏規
1
要約
本研究の目的は、時系列画像のフーリエ変換をベースとした相互相関係数から移動体の軌道ベ
クトルを計測する基礎研究を行い、開発したアルゴリズムを使用し、静止気象衛星の画像から雲
の動きのトレースに応用することである。一つの移動体画像を moving window で切り出し、もう
一方の画像では moving window を 1 ピクセルずつずらしながら、その 2 つの画像の相互相関関数
より相関値を求め、その相関値をピクセル値とした画像（相関画像）を作成する。切り出した画
像が違っていれば相関値は低く、逆に同一あるいは同様の画像の場合は高い相関値となることか
ら、移動体の検出を行う。移動体を鮮明に検出するために閾値を設け 2 値化した相関画像を利用
する。
本研究では、まず基礎理論として 5 パターンのシミュレーションを行った。その結果 moving
window は移動体と同じ大きさが最適であること、移動体の画像強度が低くなるにつれて相関値も
低くなるという結果が得られた。そして、開発した手法を気象衛星画像に適用させたところ、台
風の目のトレースに成功した。しかし、一般的な雲では時間的に形状が変化するため、本手法は
適用できないことが判明した。
2
目次
第1章
第2章
序論
1.1
研究の背景・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
（3）
1.2
研究の目的・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
（3）
1.3
論文の構成・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
（3）
基礎理論
2.1
フーリエ変換・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
（4）
2.2
相互相関関数・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
（4）
2.3
第3章
2.2.1
相互相関関数を求める手順・・・・・・・・・・・・・・・・・・・・
（5）
2.2.2
相互相関画像の導出方法・・・・・・・・・・・・・・・・・・・・・
（5）
シミュレーション・・・・・・・・・・・・・・・・・・・・・・・・・・・
（7）
2.3.1
シミュレーション 1（Basic）
・・・・・・・・・・・・・・・・・・・
（7）
2.3.2
シミュレーション 2（明度の変化）
・・・・・・・・・・・・・・・・・
（8）
2.3.3
シミュレーション 3（移動体の形状）
・・・・・・・・・・・・・・・・
（9）
2.3.4
シミュレーション 4（window size の変化）
・・・・・・・・・・・・・
（10）
2.3.5
シミュレーション 5（正しい移動体を発見できるか）
・・・・・・・・
（11）
気象データへの応用
3.1
気象衛星データについて・・・・・・・・・・・・・・・・・・・・・・・
（12）
3.2
気象衛星データを使用して CCF を抽出（気象衛星データ 1）
・・・・・・・
（13）
3.3
二値化について・・・・・・・・・・・・・・・・・・・・・・・・・・・（14）
3.4
気象衛星データ 2（二値化画像の重ね合わせ）
・・・・・・・・・・・・・・
（15）
3.5
気象衛星データ 3（二値化画像の重ね合わせ、失敗例）
・・・・・・・・・・
（17）
第４章 結論・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
（20）
謝辞・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
（20）
参考文献・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
（20）
付録・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
（21）
3
第１章
序論
1.1 研究の背景
近年、気象衛星ひまわり 5 号が引退し後継機として気象衛星ひまわり 6 号が打ち上げられ、さ
まざまな衛星画像を入手できるようになり、特に雲画像は毎日の生活での天気予報に欠かせない
ものとなってきた。しかし、その画像は最終的に人間の目で雲の移動を予想している。将来的に
は全ての作業をコンピュータによって自動化する必要がある。
1.2
研究の目的
本研究の目的は、時系列画像のフーリエ変換をベースとした相互相関係数から移動体の軌道ベ
クトルを計測する基礎研究を行い、開発したアルゴリズムを使用し静止気象衛星の画像から雲の
動きのトレースに応用することである。
まず数値解析ソフトを使用して、プログラムを作成する。そして、そのプログラムを使用して
シミュレーションを 5 つパターンに分けて作成し、相互相関画像を作成する。最後に応用として
気象衛星データで相互相関画像を作成し、雲のトレースを確認し軌道ベクトルを生成する。
1.3
論文の構成
第 2 章では本研究を行うにあたり重要と思われる基礎理論（フーリエ変換、相互相関関数）と
シミュレーションを解説する。第 3 章では、実際に気象衛星データで解析し、第 4 章で結論を算
出する。
また、付録としてシミュレーションと実データで使用した MATLAB のプログラムを載ている。
4
第２章 基礎理論
2.1 フーリエ変換（Fourier Transform）
フーリエ変換とは空間的に変化する入力信号の中に、どの空間周波数成分がどれだけ含まれて
いるかを抽出する処理をフーリエ変換（Fourier Transform）という。
∞
G( f ) =
∫ g ( x) exp(−2πxf )dx
( 2.1)
−∞
式（2.1）がフーリエ変換の式であり、 f は空間周波数、 x は距離である。いま、 g (x ) を解析
される信号とすると、 G ( f ) はその「フーリエスペクトル」と呼ばれる。信号 g (x) は実数である
ことが多いが、フーリエスペクトル G ( f ) は一般的に複素数であり、多くの場合その絶対値の２
2
乗 G ( f ) が評価される。これを「パワースペクトル」という。
図 2-1 は、フーリエ変換の例である。sin 関数をフーリエ変換すると二本のスペクトルを算出
することが出来る。
G( f )
図 2-1：フーリエ変換の例
2.2 相互相関関数
相互相関関数とは、異なる 2 つの信号波形に対して相関処理を行い、2 つの類似度を調べるこ
とで、式（2.2）のようにフーリエ変換を使うと非常に高速で相互相関関数を算出できる。
φ12 (τ ) = ∫ φ12 ( f ) exp(2πifτ )df = ∫ g1 ( x) g 2 ( x − τ )dx
( 2 .2 )
ここで、 φ12 と φ12 はそれぞれクロスパワースペクトルと逆フーリエ変換したクロスパワースペ
クトルで、 τ は 2 つの信号のうち一方の波形をずらした量の関数である。
5
2.2.1
相互相関関数を求める手順
相互相関関数を求める手順は二つの信号から直接相互相関関数を求めることも出来るが、処理
時間が長くなってしまう。本研究では演算が早いフーリエ変換を使って相互相関関数を算出した。
図 2-2 にあるように二つの信号をフーリエ変換しフーリエスペクトルを求める。そして、二乗す
ることでパワースペクトルを求め、最後に逆フーリエ変換をして、相互相関関数を得ることが出
来る
図 2-2：相互相関関数を求める手順
2.2.2
相互相関関数の導出方法
CCF 画像の導出方法は、まず二つのバックグラウンド画像を作成し移動体を埋め込む。次に各
画像を moving window で画像を切り出す。そのとき Look1 の moving window は固定し Look2
の moving window は移動できるようにしておく。そして、それぞれの画像でフーリエ変換
（FFT）を行い、フーリエ変換によって生成された各画像の画像スペクトルＳ1 とＳ2 を掛け
合わせ逆高速フーリエ変換（FFT-）を行い、相互相関関数（CCF）を算出する。そして、相互
相関関数の座標中心を相関画像のピクセル値とし、最後に moving window を 1 ピクセルずつ
移動させながら画像全域で同様の手順を行い、相関画像を作成する。
（図 2-3 参照）
6
図 2-3：相互相関画像を求める手順
前述にあるように Look2 では図 2-4 のように moving window を 1 ピクセルずつ横に移動させ
ながら相互相関関数を算出し相関画像を生成する。
図 2-4：Look2 の moving window について
7
2.3 シミュレーション
今回の研究では、数値解析ソフトを使用してシミュレーションデータを 5 パターン作成した。
Look1 を現在、Look2 を 1 時間後の画像と仮定する。Look1 画像に探してほしい移動体を moving
window で切り出し固定し、Look1、Look2 それぞれの画像で高速フーリエ変換（FFT）を行い
相互相関画像を作成した。
2.3.1
シミュレーション 1（Basic）
まず 1 つ目は基本的なパターンである。図 2-5 で示すように、128 ピクセル×128 ピクセルの
各 Look 画像に 4×4 の正方形の移動体を埋め込み、Look1 の移動体を moving window で固定し
Look2 の移動体を探し出すことができるかシミュレーションを行った。×印は Look1 の移動体の
位置を表している。シミュレーションをおこなった結果、相互相関画像が作成され、二つの移動
体の類似度を表す相互相関係数は 0.51 と算出された。この相互相関係数（0.51）という値を基準
としてその他のシミュレーションを見ていくことにする。
図 2-5：シミュレーション１（Basic）
8
2.3.2
シミュレーション 2（明度の変化）
シミュレーション 2 では図 2-6 で示すように、移動体の明度を変化させて相互相関画像を生成
した。明度を低くしたために相関係数はシミュレーション 1 よりも低くなり、0.26 と出力された。
図 2-6：シミュレーション 2（明度の変化）
9
2.3.3
シミュレーション 3（移動体の形状）
シミュレーション 3 では図 2-7 で示すように、移動体の形状を変化させて相互相関画像を生成
した。このシミュレーションでは相関係数はシミュレーション 1 よりも低くなり、0.44 と出力さ
れた。
図 2-7：シミュレーション 3（移動体の形状）
10
2.3.4
シミュレーション 4（window size の変化）
シミュレーション 4 では図 2-8 で示すように、シミュレーション 3 の移動体の形状は変化なく、
window size を変化（大きく）させて相互相関画像を生成した。このシミュレーションでは相関
係数はシミュレーション 3 よりも低くなり、0.33 と出力された。
図 2-8：シミュレーション 4（window size の変化）
11
2.3.5
シミュレーション 5（正しい移動体を発見できるか）
シミュレーション 5 では図 2-9 で示すように、シミュレーション 3 の移動体にもう一つ違った
移動体を付け加えて、相互相関画像を生成した。このシミュレーションでは正しい移動体、違っ
た移動体共に相関値の高い部分が出力されたが、正しい移動体のほうが相関値が高く相関係数も
0.44 と出力された。
図 2-9：シミュレーション 5（正しい移動体を発見できるか）
12
第3章
気象衛星データへの応用
3.1 気象衛星データについて
図 3-1：本研究で使用した画像
今回の研究で使用した画像は 2002 年 7 月 12 日の静止気象衛星ひまわりによる日本付近の画像の
一部を使用している。
画像中の台風（7 号）
7 月 12 日 16 時現在
進路
中心気圧
北北西
945ｈＰａ
13
3.2 気象衛星データを使用した、CCF の抽出
静止気象衛星データ 1
図 3-2：気象衛星データ 1 の CCF 画像と二値化画像
まず最初に、静止気象衛星データの画像を二枚使用し、トレースを行った。本研究では、赤で
囲んである台風の目の部分を移動体として相関をとった。なお Look1 を原画像とし、12 時間後の
画像を Look2 としている。
図 3-2 の左下の画像は Look1 と Look2 より相関をとった相互相関画像です。中下の画像がこれ
を鮮明にした二値化した相互相関画像である。この画像の白い部分が Look1 の移動体と一番相関
値が高いことを示し、そのときのピークのグラフを右に示した。なお画像中の赤い×印は Look1
の移動体の位置を示し、ベクトルでトレースを表示している。
14
3.3 二値化
図 3-3：二値化画像の生成
二値化とは図 3-3 にあるように、フーリエ変換によって作成された相互相関画像の相関値の高
い部分とそうでない部分を区別するために閾値を設け、それより高い部分を 1 低い部分を 0 とお
くことによって、移動体をより鮮明に検出する方法である。
15
3.4 静止気象衛星データ 2
二値化画像の重ね合わせ
図 3-4：気象衛星データ 2
16
次により正確にトレースするために 4 時間ごとの画像データを用いて相関をとった。
図 3-4 上段の画像は左から原画像を Look1、4 時間後を Look2、8 時間後を Look3、12 時間後
を Look4 としている。図 3-4 の二段目の画像がそれぞれ相互相関画像です。3 枚の相互相関画像
からわかるように、移動体である台風の目が時間とともに移動していることがわかる。
図 3-4 の三段目の画像は左から、Look1 と Look2、Look1 と Look3、Look1 と Look4 のピー
クのグラフで、四段目が二値化した相互相関画像である。
上図はこの二値化画像を重ね合わせたものである。これを見ると始めのトレース画像よりも細
かく正確にトレース出来ていることがわかる。
17
3.5 静止気象衛星データ 3
二値化画像の重ね合わせ（失敗例）
図 3-5：気象衛星データ 3
図 3-4 にある画像データの Look4 から更に 4 時間後のデータを用いて相関をとった。なお
Look4 から 4 時間後の画像を Look5 としている。
18
図 3-6：ピーク値の時間的変化
図 3-6 はピーク値のグラフをに示している。図のグラフでは Look5 のピークにばらつきがある
ことがわかる。
図 3-7：二値化した CCF 画像
図 3-7 は二値化した相互相関画像を示している。このとき Look5 は 3 つの相関が出ていること
が確認できる。
19
1
2
3
図 3-8
重ね合わせた二値化相互相関画像軌道ベクトル（トレース）
2 の囲みのところに検出されている Look5 の移動体である台風の目の相関が他のものと比べる
と低いことがわかる。これは時間の経過とともに台風の目が崩れてきたためといえる。また 1 と
3 で囲んだ所に相関が出ている。これは Look1 の移動体に類似した雲が Look5 にあったために相
関したものと思われる。
20
第4章
結論
移動体検出アルゴリズムの開発におけるシミュレーションによる検証の結論としては、移動体
の画像強度が時間と共に低下するにしたがって、相互相関係数のピーク値は低下することがあげ
られる。Moving Window を大きくするとピーク値は低下し、最適な Moving Window のサイズ
は移動体のサイズに相当する。形状の異なる移動体でも相関があるが、正しい移動体の相関係数
のピーク値が大きいことから移動体の形状が変化すると、本手法の適用が困難であることが判明
した。
実データへの応用の結論として、台風の目の形があまり変化しないある一定の時間内では、鮮
明なトレースの自動検出に成功した。しかし一定の時間を越えるとシミュレーションの結果と同
様に台風の目が崩れてくるために正確にトレースできなかった。時間的変化の少ない台風の目の
トレースに限らず、形状が一定な全ての移動体（氷山、船舶など）の軌道ベクトルの自動算出に
利用可能である。
今後の課題として、相互相関画像を二値化する際にどこに閾値を設けるかの正確な決め方や特
異データ（方向、速度などから）の相関を出さないアルゴリズムの開発が挙げられる。
謝辞
本研究は、高知工科大学物質環境システム工学科の大内和夫教授のご指導のもとで行いました。
本研究を進めるにあたり、データの解析、卒業研究発表及び論文作成に際し、ご指導ならびに
ご助言をいただいた大内和夫教授に厚く御礼申し上げ、ここに感謝の意を表します。
本研究における、データの解析、卒業研究発表及び論文作成に際し、高知工科大学物質環境シ
ステム工学科の王海鵬氏、原利允氏、泉宏樹氏、竹中絵里子氏に多大なご協力を賜り、厚く御礼
申し上げます。
本研究における画像の引用には財団法人気象業務支援センター発行、気象庁提供の雲解析事例
集を高知工科大学で画像化したものです。
参考文献
[1] フーリエの冒険、言語交流研究所
ピッポファミリークラブ、1988.
[2] 理工学生のための MATLAB ビギナーズガイド、デローレス・M・エッター、2001.
[3] Matlab の総合応用、高谷邦夫、森北出版、2002.
[4] 雲解析事例集 2001・2002 の顕著事例、気象衛星センター
（財）気象業務支援センター、2003.
21
付録
シミュレーション 1 のプログラム
%シミュレーション 1（Basic）
win=7;
%window size の大きさ
m=128;
n=128;
win1=zeros(128,128);
win1(2:6,2:6)=1;
%移動体の大きさ
win2=zeros(128,128);
win2(62:66,62:66)=1;
%移動体の大きさ
ws1=win1(1:7,1:7);
%moving window の位置
ws2=zeros(win,win);
f1=fft2(ws1);
output=zeros((m-win+1),(win-n+1));
for M=1:(m-win+1)
%相互相関関数
for N=1:(n-win+1)
xstart=M;
xend=M+win-1;
ystart=N;
yend=N+win-1;
ws2=win2(ystart:yend,xstart:xend);
f2=fft2(conj(ws2));
xc=abs((1/win)^2*ifft2(f1.*conj(f2)));
output(N,M)=xc(1,1);
end
end
22
figure
axis off
axis square
subplot(2,2,1);
%Look1 の表示
imshow(win1);
subplot(2,2,2);
%Look2 の表示
imshow(win2);
subplot(2,2,3);
%相互相関画像の表示
imshow(output);
subplot(2,2,4);
%相互相関画像の表示
mesh(output);
axis([0,128,0,128,0,1]);
max(max(output))
%ピーク値の表示
シミュレーション 2 のプログラム
%シミュレーション 2（明度の変化）
win=7;
%window size の大きさ
m=128;
n=128;
win1=zeros(128,128);
win1(2:6,2:6)=1;
%移動体の大きさ
win2=zeros(128,128);
win2(62:66,62:66)=0.5;
%移動体の大きさ
ws1=win1(1:7,1:7);
%moving window の位置
ws2=zeros(win,win);
23
f1=fft2(ws1);
output=zeros((m-win+1),(win-n+1));
for M=1:(m-win+1)
%相互相関関数
for N=1:(n-win+1)
xstart=M;
xend=M+win-1;
ystart=N;
yend=N+win-1;
ws2=win2(ystart:yend,xstart:xend);
f2=fft2(conj(ws2));
xc=abs((1/win)^2*ifft2(f1.*conj(f2)));
output(N,M)=xc(1,1);
end
end
figure
axis off
axis square
subplot(2,2,1);
%Look1 の表示
imshow(win1);
subplot(2,2,2);
%Look2 の表示
imshow(win2);
subplot(2,2,3);
%相互相関画像の表示
imshow(output);
subplot(2,2,4);
%相互相関画像の表示
mesh(output);
axis([0,128,0,128,0,1]);
24
max(max(output))
%ピーク値の表示
シミュレーション 3 のプログラム
%シミュレーション 3（移動体の形状の変化）
win=13;
%window size の大きさ
m=128;
n=128;
win1=zeros(128,128);
win1(14:18,8:12)=1;
%移動体の大きさ
win1(14:18,14:18)=1;
win1(20:24,14:18)=1;
win2=zeros(128,128);
win2(107:111,26:30)=1;
%移動体の大きさ
win2(101:105,26:30)=1;
win2(101:105,20:24)=1;
ws1=win1(13:25,7:19);
%moving window の位置
ws2=zeros(win,win);
f1=fft2(ws1);
output=zeros((m-win+1),(win-n+1));
for M=1:(m-win+1)
%相互相関関数
for N=1:(n-win+1)
xstart=M;
xend=M+win-1;
ystart=N;
yend=N+win-1;
ws2=win2(ystart:yend,xstart:xend);
25
f2=fft2(conj(ws2));
xc=abs((1/win)^2*ifft2(f1.*conj(f2)));
output(N,M)=xc(1,1);
end
end
figure
axis off
axis square
subplot(2,2,1);
%Look1 の表示
imshow(win1);
subplot(2,2,2);
%Look2 の表示
imshow(win2);
subplot(2,2,3);
%相互相関画像の表示
imshow(output);
subplot(2,2,4);
%相互相関画像の表示
mesh(output);
axis([0,128,0,128,0,1]);
max(max(output))
%ピーク値の表示
シミュレーション 4 のプログラム
%シミュレーション 4（window size の変化）
win=15;
%window size の大きさ
m=128;
n=128;
win1=zeros(128,128);
26
win1(14:18,8:12)=1;
%移動体の大きさ
win1(14:18,14:18)=1;
win1(20:24,14:18)=1;
win2=zeros(128,128);
win2(107:111,26:30)=1;
%移動体の大きさ
win2(101:105,26:30)=1;
win2(101:105,20:24)=1;
ws1=win1(12:26,6:20);
%moving window の位置
ws2=zeros(win,win);
f1=fft2(ws1);
output=zeros((m-win+1),(win-n+1));
for M=1:(m-win+1)
%相互相関関数
for N=1:(n-win+1)
xstart=M;
xend=M+win-1;
ystart=N;
yend=N+win-1;
ws2=win2(ystart:yend,xstart:xend);
f2=fft2(conj(ws2));
xc=abs((1/win)^2*ifft2(f1.*conj(f2)));
output(N,M)=xc(1,1);
end
end
figure
axis off
axis square
subplot(2,2,1);
%Look1 の表示
27
imshow(win1);
subplot(2,2,2);
%Look2 の表示
imshow(win2);
subplot(2,2,3);
%相互相関画像の表示
imshow(output);
subplot(2,2,4);
%相互相関画像の表示
mesh(output);
axis([0,128,0,128,0,1]);
max(max(output))
%ピーク値の表示
シミュレーション 5 のプログラム
%シミュレーション 5（正しい移動体を見つけられるか）
win=13;
%window size の大きさ
m=128;
n=128;
win1=zeros(128,128);
win1(14:18,8:12)=1;
%移動体の大きさ
win1(14:18,14:18)=1;
win1(20:24,14:18)=1;
win2=zeros(128,128);
win2(90:94,90:94)=1;
%移動体の大きさ
win2(84:88,84:88)=1;
win2(84:88,96:100)=1;
win2(107:111,26:30)=1;
%移動体の大きさ
win2(101:105,26:30)=1;
win2(101:105,20:24)=1;
28
ws1=win1(13:25,7:19);
%moving window の位置
ws2=zeros(win,win);
f1=fft2(ws1);
output=zeros((m-win+1),(win-n+1));
for M=1:(m-win+1)
%相互相関関数
for N=1:(n-win+1)
xstart=M;
xend=M+win-1;
ystart=N;
yend=N+win-1;
ws2=win2(ystart:yend,xstart:xend);
f2=fft2(conj(ws2));
xc=abs((1/win)^2*ifft2(f1.*conj(f2)));
output(N,M)=xc(1,1);
end
end
figure
axis off
axis square
subplot(2,2,1);
%Look1 の表示
imshow(win1);
subplot(2,2,2);
%Look2 の表示
imshow(win2);
subplot(2,2,3);
%相互相関画像の表示
imshow(output);
subplot(2,2,4);
%相互相関画像の表示
29
mesh(output);
axis([0,128,0,128,0,1]);
max(max(output))
%ピーク値の表示
%静止気象衛星データ 1
a=imread('kumo20.bmp','bmp');
%Look1 画像
b=imread('kumo19.bmp','bmp');
%Look2 画像
win=6;
%window size の大きさ
n=128;
m=128;
img1=double(a)/255;
img10=imresize(img1,[128 128]);
aa=rgb2gray(img10);
img2=double(b)/255;
img20=imresize(img2,[128 128]);
bb=rgb2gray(img20);
win1=aa(116:121,51:56);
%Look1 の moving window の位置
win2=zeros(win,win);
%Look2 の moving window の位置
f1=fft2(win1);
output=zeros((m-win+1),(win-n+1));
%相互相関関数
for M=1:(n-win+1)
for N=1:(n-win+1)
xstart=M;
xend=M+win-1;
ystart=N;
yend=N+win-1;
30
win2=bb(ystart:yend,xstart:xend);
f2=fft2(conj(win2));
xc=abs((1/win)^2*ifft2(f1.*conj(f2)));
i1=mean2(win1);
i2=mean2(win2);
ccf=xc/(i1*i2);
output(N,M)=ccf(1,1);
end
end
figure
subplot(2,2,1);
%Look1 の表示
imshow(aa);
subplot(2,2,2);
%Look2 の表示
imshow(bb);
subplot(2,2,3);
%相互相関画像の表示
imagesc(output);
axis square
text(51,116,'x','color','r');
subplot(2,2,4);
%相互相関画像の表示
output2=output/1.0641;
meshz(output2);
axis([0,128,0,128,0.91,1.0]);
axis square
31
%静止気象衛星データ 2
a=imread('kumo21.bmp','bmp');
%Look1 画像
b=imread('kumo22.bmp','bmp');
%Look2 画像
c=imread('kumo23.bmp','bmp');
%Look3 画像
d=imread('kumo24.bmp','bmp');
%Look4 画像
win=6;
%window size の大きさ
n=128;
m=128;
img1=double(a)/255;
img10=imresize(img1,[128 128]);
aa=rgb2gray(img10);
img2=double(b)/255;
img20=imresize(img2,[128 128]);
bb=rgb2gray(img20);
img3=double(c)/255;
img30=imresize(img3,[128 128]);
cc=rgb2gray(img30);
img4=double(d)/255;
img40=imresize(img4,[128 128]);
dd=rgb2gray(img40);
win1=aa(116:121,51:56);
%Look1 の moving window の位置
win2=zeros(win,win);
%Look2 の moving window の位置
win3=zeros(win,win);
%Look3 の moving window の位置
win4=zeros(win,win);
%Look4 の moving window の位置
f1=fft2(win1);
output=zeros((m-win+1),(win-n+1));
%相互相関関数
32
for M=1:(n-win+1)
for N=1:(n-win+1)
xstart=M;
xend=M+win-1;
ystart=N;
yend=N+win-1;
win2=bb(ystart:yend,xstart:xend);
win3=cc(ystart:yend,xstart:xend);
win4=dd(ystart:yend,xstart:xend);
f2=fft2(conj(win2));
f3=fft2(conj(win3));
f4=fft2(conj(win4));
xc1=abs((1/win)^2*ifft2(f1.*conj(f2)));
xc2=abs((1/win)^2*ifft2(f1.*conj(f3)));
xc3=abs((1/win)^2*ifft2(f1.*conj(f4)));
i1=mean2(win1);
i2=mean2(win2);
i3=mean2(win3);
i4=mean2(win4);
ccf1=xc1/(i1*i2);
ccf2=xc2/(i1*i3);
ccf3=xc3/(i1*i4);
output1(N,M)=ccf1(1,1);
output2(N,M)=ccf2(1,1);
output3(N,M)=ccf3(1,1);
end
end
figure
33
subplot(2,2,1);
%Look1 の表示
imshow(aa);
subplot(2,2,2);
%Look2 の表示
imshow(bb);
subplot(2,2,3);
%Look3 の表示
imshow(cc);
subplot(2,2,4);
%Look4 の表示
imshow(dd);
%相互相関画像の表示
figure
subplot(2,3,1);
imagesc(output1);
axis square
colormap(gray);
subplot(2,3,2);
imagesc(output2);
axis square
subplot(2,3,3);
imagesc(output3);
axis square
subplot(2,3,4);
output4=output1/1.0637;
meshz(output4);
axis([0,123,0,123,0.9,1.01]);
axis square
subplot(2,3,5);
output5=output2/1.0613;
meshz(output5);
axis([0,123,0,123,0.9,1.01]);
axis square
34
subplot(2,3,6);
output6=output3/1.0629;
meshz(output6);
axis([0,123,0,123,0.9,1.01]);
axis square
%静止気象衛星データ 3
a=imread('kumo21.bmp','bmp');
%Look1 画像
b=imread('kumo22.bmp','bmp');
%Look2 画像
c=imread('kumo23.bmp','bmp');
%Look3 画像
d=imread('kumo24.bmp','bmp');
%Look4 画像
e=imread('kumo25.bmp','bmp');
%Look5 画像
win=6;
%window size の大きさ
n=128;
m=128;
img1=double(a)/255;
img10=imresize(img1,[128 128]);
aa=rgb2gray(img10);
img2=double(b)/255;
img20=imresize(img2,[128 128]);
bb=rgb2gray(img20);
img3=double(c)/255;
img30=imresize(img3,[128 128]);
cc=rgb2gray(img30);
img4=double(d)/255;
img40=imresize(img4,[128 128]);
dd=rgb2gray(img40);
img5=double(e)/255;
img50=imresize(img5,[128 128]);
35
ee=rgb2gray(img50);
win1=aa(116:121,51:56);
%Look1 の moving window の位置
win2=zeros(win,win);
%Look2 の moving window の位置
win3=zeros(win,win);
%Look3 の moving window の位置
win4=zeros(win,win);
%Look4 の moving window の位置
win5=zeros(win,win);
%Look5 の moving window の位置
f1=fft2(win1);
output=zeros((m-win+1),(win-n+1));
%相互相関関数
for M=1:(n-win+1)
for N=1:(n-win+1)
xstart=M;
xend=M+win-1;
ystart=N;
yend=N+win-1;
win2=bb(ystart:yend,xstart:xend);
win3=cc(ystart:yend,xstart:xend);
win4=dd(ystart:yend,xstart:xend);
win5=ee(ystart:yend,xstart:xend);
f2=fft2(conj(win2));
f3=fft2(conj(win3));
f4=fft2(conj(win4));
f5=fft2(conj(win5));
xc1=abs((1/win)^2*ifft2(f1.*conj(f2)));
xc2=abs((1/win)^2*ifft2(f1.*conj(f3)));
xc3=abs((1/win)^2*ifft2(f1.*conj(f4)));
xc4=abs((1/win)^2*ifft2(f1.*conj(f5)));
i1=mean2(win1);
i2=mean2(win2);
i3=mean2(win3);
36
i4=mean2(win4);
i5=mean2(win5);
ccf1=xc1/(i1*i2);
ccf2=xc2/(i1*i3);
ccf3=xc3/(i1*i4);
ccf4=xc4/(i1*i5);
output1(N,M)=ccf1(1,1);
output2(N,M)=ccf2(1,1);
output3(N,M)=ccf3(1,1);
output4(N,M)=ccf4(1,1);
end
end
figure
subplot(3,2,1);
%Look1 の表示
imshow(aa);
subplot(3,2,2);
%Look2 の表示
imshow(bb);
subplot(3,2,3);
%Look3 の表示
imshow(cc);
subplot(3,2,4);
%Look4 の表示
imshow(dd);
subplot(3,2,5);
%Look5 の表示
imshow(ee);
%相互相関画像の表示
figure
subplot(2,2,1);
imagesc(output1);
37
colormap(gray);
axis square
subplot(2,2,2);
imagesc(output2);
axis square
subplot(2,2,3);
imagesc(output3);
axis square
subplot(2,2,4);
imagesc(output4);
axis square
figure
subplot(2,2,1);
output6=output1/1.0637;
meshz(output6);
axis([0,128,0,128,0.92,1.01]);
axis square
colormap(gray);
subplot(2,2,2);
output7=output2/1.0613;
meshz(output7);
axis([0,128,0,128,0.92,1.01]);
axis square
subplot(2,2,3);
output8=output3/1.0629;
meshz(output8);
axis([0,128,0,128,0.92,1.01]);
axis square
subplot(2,2,4);
output9=output4/1.0465;
meshz(output9);
axis([0,128,0,128,0.92,1.01]);
axis square
figure
output5=output1+output2+output3+output4;
38
subplot(1,1,1);
output10=output5/4.0783;
colormap(gray);
meshz(output10);
axis([0,128,0,128,0.983,1.01]);
axis square
39