ラドン変換を用いた回転・拡大の算出（プログラミング言語はMatlabです。）

現在、画像I0に対してI1の画像がどれだけ回転・拡大しているのかというのを算出するプログラムを組もうとしていますが、プログラム途中でゆきずまっています。（ただし、I1はI0をアフィン変換した画像です。）
はじめに理論から説明して、次に自分が作っている途中のプログラムをお見せしようと思います。
以下理論の説明です。

画像I0(x,y)をアフィン変換した画像をI1(x',y')とすると
I0(x,y)＝I1(x',y')　　（１）
x'=s*(x*cosφ-y*sinφ)+Δx　　　（２）
y'=s*(x*sinφ+y*cosφ)+Δy　　　（３）
となります。ここで、s,φ,Δx,Δyはそれぞれ拡大率、回転角、平行移動量とします。

　画像それぞれに対して、ラドン変換します。
r0(θ,ρ)=∬I0(x,y)δ(ρ-xcosθ-ysinθ)dxdy　　　（４）
r1(θ,ρ)=∬I1(x',y')δ(ρ-x'cosθ-y'sinθ)dx'dy'　　　（５）

上の４式をとくと、
r1(θ,ρ)＝r0(θ-φ,{ρ-tcos(θ-k)}/s)　　　（６）
t=sqrt(Δx^2+Δy^2)　　　（７）
k=tan^-1(Δx/Δy)　　　（８）
となります。

次に、ラドン変換平面をρ軸方向に一次元フーリエ変換した後にパワースペクトルを計算すると、r0(θ,ρ) , r1(θ,ρ)はそれぞれ
　R0(θ,f)=|∫r0(θ,ρ)exp(-j2πfρ)dρ|　　　（９）
　R1(θ,f)=|∫r1(θ,ρ)exp(-j2πfρ)dρ|　　　（１０）
　 　　
となります。
(5),(9),(10)より、
r1(θ,f)＝r0(θ-φ,sf)　
の関係が成り立ち、平行移動の影響を取り除けるます。
　ここで、
r1(θ,f)＝r0(θ-φ,q-λ)　　　　(１１)
q=ln(f)　　　(１２)
λ=-ln(s)　　　(１３)
のようにf軸の対数変換を行なうと、拡大、回転の影響がθ、q軸のシフト量として表せます。


これらの理論を踏まえた上で質問が２点あります。
上の対数変換をどうプログラミングすればよいのか分からないのですがどうすればよいでしょうか？

あと、はじめのアフィン変換の(2)式ですが、これをプログラムしようとすると拡大率sやらΔｘ、Δｙなどを定義しなくてはならないのですが、これらは最後に求めるもの、あるいは計算途中で必要なくなるものであって、定義のしようがないのですがどうすればよいでしょうか？

プログラム言語はMatlabです。プログラミング初心者で至らないところがあると思いますが助言の方よろしくお願いします。


以下、プログラミング途中のものです。
＊ラドン変換をハフ変換にかえてあります。
A1=imread('tekubi-small.jpg');
A2=imread('tekubi-small.jpg');
A1=imresize(A1,[128 128]); A2=imresize(A2,[128 128]);
B1=rgb2gray(A1);B2=rgb2gray(A2);B2=rot90(B2);
画像Ｂ１とをＢ１を９０度回転させたＢ２とで回転角・拡大率を算出できるか実験しています。
B1=edge(B1,'canny');
B2=edge(B2,'canny');

%ハフ変換
[H, THETA, RHO] = hough(B1);
[H2, THETA2, RHO2] = hough(B2);
%ハフ変換平面をρ軸方向に一次元フーリエ変換する
C1=fft(H,[],1);
C2=fft(H2,[],1);
%パワースペクトルを計算する
C11=abs(C1);
C22=abs(C2);
というところまでやりました。

No.1 ベストアンサー 回答者： rabbit_cat 回答日時： 2008/07/21 20:02

仰っている理論通りなら、

明示的にに対数変換の計算をしないでも、そのまま、C11とC22を比べてx方向とy方向にどれだけシフトさせたときに一番ぴったり合うかを調べて、そのシフト量から、φとλを算出すればいいと思います。
ほんとは対数変換するってことを考えると、どこがぴったり合うか調べるときに、低周波成分を高周波成分よりも重みを高くして適合度を調べるとかしたほうがいいかな。

この回答へのお礼

書き込みありがとうございます！

実は上に書いた理論だけがすべてじゃなくて続きがあります。
　確かに上の理論だけ見るとシフト量から求められそうですが、(説明不足で申し訳ないです)それだけじゃダメみたいで、の後に、位相限定相関法という手法を用いてシフト量を計算するらしいですが、その前段階でつまづいてしまっている状況です。
　
上のプログラムからどうやったらできるか教えていただけるとうれしいです。

お礼日時：2008/07/24 18:47