質問: フーリエ変換における問題点がわかりません。教えてもらえれば幸いですが、この問題点について説明しているホームページを教えてほしいです。よろしくお願いします。
回答: 2次元画像に関する、平行投影からの画像再構成法の話かと思います。投影からの再構成が二次元のフーリエ変換を使ってきれいに定式化できるということを最初に発見したのが誰だかは知らないですが、理屈はとても簡単ですんで、解説してみましょう。 [1] 画像を、直交座標系x-y上の2次元密度分布f(x,y)で表します。さて、f(x,y)を二次元フーリエ変換したものをF(α,β)としましょう。Fを作るには、まずf(x,y)をyだけについて一次元フーリエ変換したものg(x,β)を作り、 g(x,β)=∫f(x,y) exp(-iyβ)dy (積分はy=-∞~∞の定積分) さらにg(x,β)をxだけについて一次元フーリエ変換したものF(α,β)を作る F(α,β)=∫g(x,β) exp(-ixα)dx (積分はx=-∞~∞の定積分) というやりかたで出来ます。 exp(-ixα)exp(-iyβ)=exp(-i(xα+yβ)) を利用して、まとめて書けば F(α,β)=∫∫f(x,y) exp(-i(xα+yβ))dydx (積分はx=-∞~∞, y=-∞~∞の定積分) です。 [2] ここで、Fを原点を中心にしてθだけ回転した座標系で見たものをFθと書くことにしましょう。すなわち A=αcosθ-βsinθ B=αsinθ+βcosθ という回転変換を考える訳です。逆変換は(θ→-θ, α→A, β→B, A→α,B→βと置き換えればよく、) α=Acosθ+Bsinθ β=-Asinθ+Bcosθ だから、fを二次元フーリエ変換してから、回転して見たものFθ(A,B)は Fθ(A,B)=F(Acosθ+Bsinθ,-Asinθ+Bcosθ)=F(α,β) と表せます。代入してみれば、 Fθ(A,B)=∫∫f(x,y) exp(-i((Acosθ+Bsinθ)x+(-Asinθ+Bcosθ)y))dydx (積分はx=-∞~∞, y=-∞~∞の定積分) ですね。 [3]画像f(x,y)を原点を中心に角度θだけ回転して見たものをfθ(X,Y)と書けば、同様に X=xcosθ-ysinθ Y=xsinθ+ycosθ ですから、 fθ(X,Y)=f(Xcosθ+Ysinθ,-Xsinθ+Ycosθ)=f(x,y) です。 さて、fθを二次元フーリエ変換したものをUθとしましょう。 Uθ(A,B)=∫∫fθ(X,Y) exp(-i(XA+YB))dYdX (積分はX=-∞~∞, Y=-∞~∞の定積分) だから、 Uθ(A,B)=∫∫fθ(xcosθ-ysinθ,xsinθ+ycosθ) exp(-i((xcosθ-ysinθ)A+(xsinθ+ycosθ)B))dYdX =∫∫f(x,y) exp(-i((xcosθ-ysinθ)A+(xsinθ+ycosθ)B))dYdX =∫∫f(x,y) exp(-i(xAcosθ-yAsinθ+xBsinθ+yBcosθ))dYdX =∫∫f(x,y) exp(-i((Acosθ+Bsinθ)x+(-Asinθ+Bcosθ)y))dYdX =∫∫f(x,y) exp(-i((Acosθ+Bsinθ)x+(-Asinθ+Bcosθ)y))dydx =Fθ(A,B) つまり、「二次元フーリエ変換してから回転しても、回転してから二次元フーリエ変換しても同じ」ということです。 なおここで、dYdX = dxdyは、積分変数を変える置換積分の操作(ヤコビアンを使う)によって示されます。 [4] 一方、投影とは、画像f(x,y)を原点を中心に角度θだけ回転して見たものをfθ(X,Y)を、Yだけについて積分したものです。すなわち角度θ方向への投影とは pθ(X)=∫fθ(X,Y)dY (積分はY=-∞~∞の定積分) です。さらにこれをXについてフーリエ変換したものをPθとすると、 Pθ(A)=∫pθ(X) exp(-iXA) dX (積分はX=-∞~∞の定積分) =∫∫fθ(X,Y) exp(-iXA) dYdX です。この式と Uθ(A,B)=Fθ(A,B)=∫∫fθ(X,Y) exp(-i(XA+YB))dYdX とを見比べれば直ちに、 Pθ(A)=Fθ(A,0) と分かりますね。投影とは、Yに沿ってB=0(直流)の成分を取り出す、ということである。言われてみれば当たり前、という訳です。 そしてこれこそが、画像再構成の基本原理である「中央断面定理(中央切断定理、central slice theorem)」なんです。 角度θ方向の投影pθ(X)をフーリエ変換したものPθ(A)を作って、これをF(α,β)の二次元平面上で、角度θ方向の直径に沿って並べてやる。これがFθ(A,0)です。θ=0~πについてこの操作をやればF(α,β)が至るところ決定できるから、あとはこいつを二次元逆フーリエ変換すればf(x,y)が再現できることになる。 平行投影(parallel beam projection)の場合のconvolution-backprojection法は、この定理からさらにごにょごにょ計算すると出てきます。 また、convergent beam projectionの場合のconvolution-backprojection法は、平行投影の場合のconvolution-backprojection法から導かれます。 [5] 画像再構成理論には、この他にもフーリエ変換が出てくる場面があります。それは、平行投影の場合のconvolution-backprojection法において、畳み込み積分(convolution)をフーリエ変換で計算する、という部分です。また、convergent beam projectionの場合のconvolution-backprojection法において、畳み込み積分(convolution)をフーリエ変換で計算する、という部分にも出てきますが、こちらは、厳密に言えば「周期的畳み込み積分を離散化した積和計算を速く行うのに、離散フーリエ変換を用いる」と言った方が正確です。どちらについても、フーリエ変換は畳み込み積分の計算手段に過ぎないのであって、上記の中央断面定理におけるような本質的な役割を持っている訳ではありません。 [6] さて、中央断面定理に戻り、この定理をそのまま数値計算でやるとどうか、を考えます。 すると、実際にpθ(X)のデータが取れるθは有限個である。離散的にしかデータが得られません。さらにXも離散的です。従って、F(α,β)の至る所の値が分かる訳ではない。分かるのは、F(α,β)の平面上に描いた極座標の格子点に於ける値だけである。F(α,β)を二次元離散フーリエ変換するには、直交座標の格子点上の値が必要なんだけど、それは得られません。 だから、直交座標の格子点での値を、最寄りの分かっている値から内挿して推定する必要があります。「理想的にやれば正確な推定が可能である」ということはサンプリング定理を使ったいささか込み入った議論によって保証できますが、実際の計算ではそうは行かず、推定には誤差が生じます。 Pθ(A)のAの絶対値が小さい部分(つまりpθ(X)の低周波成分)については、(α,β)平面上にデータを置くと、すんごくたくさんのデータがぎっしりあることになり、一方、Pθ(A)のAの絶対値が大きい部分(つまりpθ(X)の高周波成分)については、半径Aが大きい分だけ(α,β)の平面でサンプル点の間隔が大きくなる。このような不均衡が生じるわけです。このために、高周波成分に関しては直交座標の格子点上の値を推定する際の誤差がどうしても大きくなり、従って、F(α,β)を二次元離散フーリエ変換して得られる画像には、本来存在しないはずの嘘の像(偽像 artifact)が生じることになります。(多分、ご質問はこの偽像のことじゃないか、と思うんですが、どうでしょうかね。) 推定の方法として最も簡単なのは、値を求めたい直交座標の格子点 c に対して、それを囲む、値が分かっている極座標の格子点4つを取って、双一次式で内挿することです。でも、これではartifactが強く出てしまう。cに近い格子点をもっとたくさん使って、たとえば双三次式で内挿してやるとましになります。(そして、上記の理想的な推定では、極座標の格子点全部を使う訳です。)