AlxisDataBlog-数学・物理学


MathJax = {
  chtml: {
    matchFontHeight: false
  },
  tex: {
    inlineMath: [['$', '$']]
  }
};






今回は、ただひたすらに計算していくだけの回です。回転積分した反転公式

\begin{align}
8\pi^3f(z)=\int_0^{2\pi}\bigg\{\mathrm{pv.}\int_0^{2\pi}\dfrac{1+\cos(\theta-\tau)}{\sin(\theta-\tau)}\mathscr{H}f(\theta+\pi,z-re(\tau)+re(\theta))\,d\theta\bigg\}\,d\tau+2\pi\int_0^{2\pi}\psi(r,z-re(\tau))\,d\tau
\end{align}

からスタートしましょう。右辺の第 $2$ 項は既に半径 $r$ のフィルタになっているので、第 $1$ 項(以下 $A$ とおきます)だけ注目していきます。

\begin{align}
A&=\int_0^{2\pi}\bigg\{\mathrm{pv.}\int_{-\pi}^{\pi}\dfrac{1+\cos\theta}{\sin\theta}\mathscr{H}f(\theta+\tau+\pi,z-re(\tau)+re(\theta+\tau))\,d\theta\bigg\}\,d\tau\\
&=\int_0^{2\pi}\bigg\{\int_0^{\pi}\dfrac{1+\cos\theta}{\sin\theta}\{\mathscr{H}f(\theta+\tau+\pi,z-re(\tau)+re(\theta+\tau))\\
&\quad-\mathscr{H}f(-\theta+\tau+\pi,z-re(\tau)+re(-\theta+\tau))\}\,d\theta\bigg\}\,d\tau\\
&=\int_0^{\pi}\dfrac{1+\cos\theta}{\sin\theta}\bigg\{\int_0^{2\pi}\{\mathscr{H}f(\theta+\tau+\pi,z-re(\tau)+re(\theta+\tau))\\
&\quad-\mathscr{H}f(-\theta+\tau+\pi,z-re(\tau)+re(-\theta+\tau))\}\,d\tau\bigg\}\,d\theta
\end{align}

最後の項を $\tau=\tau+\theta-\pi$ に置き換えて

\begin{align}
A&=\int_0^{\pi}\dfrac{1+\cos\theta}{\sin\theta}\bigg\{\int_0^{2\pi}\{\mathscr{H}f(\theta+\tau+\pi,z-re(\tau)+re(\theta+\tau))\\
&\quad-\mathscr{H}f(\tau,z-re(\tau+\theta-\pi)+re(\tau-\pi))\}\,d\tau\bigg\}\,d\theta\\
&=\int_0^{\pi}\dfrac{1+\cos\theta}{\sin\theta}\bigg\{\int_0^{2\pi}\{\mathscr{H}f(\theta+\tau+\pi,z-re(\tau)+re(\theta+\tau))\\
&\quad-\mathscr{H}f(\tau,z+re(\tau+\theta)-re(\tau))\}\,d\tau\bigg\}\,d\theta
\end{align}

$\mathscr{H}f$ の引数を書き換えます。

\begin{align}
-e(\tau)+e(\theta+\tau)&=(-\cos\tau+\cos(\theta+\tau),-\sin\tau+\sin(\theta+\tau))\\
&=(-2\sin(\theta/2)\sin(\tau+\theta/2),2\sin(\theta/2)\cos(\tau+\theta/2))\\
&=2\sin(\theta/2)e(\tau+\theta/2+\pi/2)
\end{align}

を利用して

\begin{align}
A&=\int_0^{\pi}\dfrac{1+\cos\theta}{\sin\theta}\bigg\{\int_0^{2\pi}\{\mathscr{H}f(\theta+\tau+\pi,z+2r\sin(\theta/2)e(\tau+\theta/2+\pi/2))\\
&\quad-\mathscr{H}f(\tau,z+2r\sin(\theta/2)e(\tau+\theta/2+\pi/2))\}\,d\tau\bigg\}\,d\theta
\end{align}

$\tau=\tau-\theta/2-\pi/2$ に置き換えて

\begin{align}
A&=\int_0^{\pi}\dfrac{1+\cos\theta}{\sin\theta}\bigg\{\int_0^{2\pi}\{\mathscr{H}f(\theta/2+\tau+\pi/2,z+2r\sin(\theta/2)e(\tau))\\
&\quad-\mathscr{H}f(\tau-\theta/2-\pi/2,z+2r\sin(\theta/2)e(\tau))\}\,d\tau\bigg\}\,d\theta
\end{align}

$\mathscr{H}f$ の差分を積分で書くと

\begin{align}
A&=-\int_0^{\pi}\dfrac{1+\cos\theta}{\sin\theta}\bigg\{\int_0^{2\pi}\bigg\{\int_{-(\pi-\theta)/2}^{(\pi-\theta)/2}\partial_{\theta}\mathscr{H}f(\varphi+\tau-\pi,z+2r\sin(\theta/2)e(\tau))\,d\varphi\bigg\}\,d\tau\bigg\}\,d\theta
\end{align}

$\partial_{\theta}\mathscr{H}f$ を

\begin{align}
\dfrac{1}{r}\partial_{\theta}\mathscr{H}f(\theta,z)=(\partial_r-\partial_{e(\theta)})\psi(r,z+re(\theta))
\end{align}

を使って書き換えると

\begin{align}
A&=-r\int_0^{\pi}\dfrac{1+\cos\theta}{\sin\theta}\bigg\{\int_0^{2\pi}\bigg\{\int_{-(\pi-\theta)/2}^{(\pi-\theta)/2}(\partial_r+\partial_{e(\varphi+\tau)})\psi(r,z-re(\varphi+\tau)+2r\sin(\theta/2)e(\tau))\,d\varphi\bigg\}\,d\tau\bigg\}\,d\theta
\end{align}

この時点でフィルタサイズは半径 $r$ になっていますが、もう少し式を簡単にしていきます。$\tau$ だけ回転させる行列を $R_{\tau}$ として

\begin{align}
&\nu_0=\dfrac{-e(\varphi)+2\sin(\theta/2)e(0)}{\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}},\\
&\nu=R_{\tau}\nu_0,\ \tilde{\nu}=R_{\tau+\pi/2}\nu_0
\end{align}

とおきましょう。

\begin{align}
&\int_0^{2\pi}(\partial_r+\partial_{e(\varphi+\tau)})\psi(r,z-re(\varphi+\tau)+2r\sin(\theta/2)e(\tau))\,d\tau\\
&=\int_0^{2\pi}\partial_r\psi(r,z+r\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}R_{\tau}\nu_0)\,d\tau\\
&\quad+\int_0^{2\pi}\langle e(\varphi+\tau),\nu\rangle\partial_{\nu}\psi(r,z+r\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}R_{\tau}\nu_0)\,d\tau\\
&\quad+\int_0^{2\pi}\langle e(\varphi+\tau),\tilde{\nu}\rangle\partial_{\tilde{\nu}}\psi(r,z+r\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}R_{\tau}\nu_0)\,d\tau\\
&=\int_0^{2\pi}\partial_r\psi(r,z+r\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}R_{\tau}\nu_0)\,d\tau\\
&\quad+\int_0^{2\pi}\langle e(\varphi),\nu_0\rangle\partial_{\nu}\psi(r,z+r\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}R_{\tau}\nu_0)\,d\tau\\
&\quad+\int_0^{2\pi}\langle e(\varphi),R_{\pi/2}\nu_0\rangle\partial_{\tilde{\nu}}\psi(r,z+r\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}R_{\tau}\nu_0)\,d\tau\\
&=\int_0^{2\pi}\partial_r\psi(r,z+r\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}R_{\tau}\nu_0)\,d\tau\\
&\quad+\int_0^{2\pi}\dfrac{-1+2\sin(\theta/2)\cos\varphi}{\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}}\partial_{\nu}\psi(r,z+r\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}R_{\tau}\nu_0)\,d\tau\\
&\quad+\int_0^{2\pi}\dfrac{2\sin(\theta/2)\sin\varphi}{\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}}\partial_{\tilde{\nu}}\psi(r,z+r\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}R_{\tau}\nu_0)\,d\tau
\end{align}

第 $1,2$ 項は $\varphi$ に関して偶関数、第 $3$ 項は奇関数となることに注意して

\begin{align}
A&=-2r\int_0^{\pi}\dfrac{1+\cos\theta}{\sin\theta}\bigg\{\int_0^{2\pi}\bigg\{\int_0^{(\pi-\theta)/2}\partial_r\psi(r,z-re(\varphi+\tau)+2r\sin(\theta/2)e(\tau))\\
&\quad +\dfrac{-1+2\sin(\theta/2)\cos\varphi}{\sqrt{1-4\sin(\theta/2)\cos\varphi+4\sin^2(\theta/2)}}\partial_{\nu}\psi(r,z-re(\varphi+\tau)+2r\sin(\theta/2)e(\tau))\,d\varphi\bigg\}\,d\tau\bigg\}\,d\theta
\end{align}

$R=2\sin(\theta/2),\rho=\cos\varphi$ と変数変換しましょう。

\begin{align}
A&=-4r\int_0^2\dfrac{1}{R}\bigg\{\int_0^{2\pi}\bigg\{\int_{R/2}^1\dfrac{1}{\sqrt{1-\rho^2}}\partial_r\psi(r,z+r\sqrt{1-2R\rho+R^2}R_{\tau}\nu_0)\\
&\quad +\dfrac{1}{\sqrt{1-\rho^2}}\dfrac{-1+R\rho}{\sqrt{1-2R\rho+R^2}}\partial_{\nu}\psi(r,z+r\sqrt{1-2R\rho+R^2}R_{\tau}\nu_0)\,d\rho\bigg\}\,d\tau\bigg\}\,dR\\
&=-4r\int_0^2\dfrac{1}{R}\bigg\{\int_0^{2\pi}\bigg\{\int_{R/2}^1\dfrac{1}{\sqrt{1-\rho^2}}\partial_r\psi(r,z+r\sqrt{1-2R\rho+R^2}e(\tau))\\
&\quad +\dfrac{1}{\sqrt{1-\rho^2}}\dfrac{-1+R\rho}{\sqrt{1-2R\rho+R^2}}\partial_{e(\tau)}\psi(r,z+r\sqrt{1-2R\rho+R^2}e(\tau))\,d\rho\bigg\}\,d\tau\bigg\}\,dR
\end{align}

$X=\sqrt{1-2R\rho+R^2}$ と変数変換しましょう。

\begin{align}
A&=-4r\int_0^1\bigg\{\int_0^{2\pi}\bigg\{\int_{1-X}^{1+X}\dfrac{2X}{R\sqrt{(2X)^2-(R^2-1-X^2)^2}}\partial_r\psi(r,z+rXe(\tau))\\
&\quad +\dfrac{2X}{R\sqrt{(2X)^2-(R^2-1-X^2)^2}}\dfrac{R^2-1-X^2}{2X}\partial_{e(\tau)}\psi(r,z+rXe(\tau))\,dR\bigg\}\,d\tau\bigg\}\,dX
\end{align}

$\sin\phi=(R^2-1-X^2)/(2X)$ 変数変換しましょう。

\begin{align}
A&=-2r\int_0^1\bigg\{\int_0^{2\pi}\bigg\{\int_{-\pi/2}^{\pi/2}\dfrac{1}{\sin\phi+(1+X^2)/(2X)}\partial_r\psi(r,z+rXe(\tau))\\
&\quad +\dfrac{\sin\phi}{\sin\phi+(1+X^2)/(2X)}\partial_{e(\tau)}\psi(r,z+rXe(\tau))\,d\phi\bigg\}\,d\tau\bigg\}\,dX
\end{align}

$\phi$ に関する積分は実行できそうですね。今回はここまで。次回も引き続きこの式を簡潔にしていきましょう。(written by J.Nishigata)


Copy Right 2020-現在 ALXIS DATA, Inc. All Right Reserved.

最小な反転公式フィルタの導出

AlxisDataBlog-数学・物理学


MathJax = {
  chtml: {
    matchFontHeight: false
  },
  tex: {
    inlineMath: [['$', '$']]
  }
};






今回から、ヒルベルト変換の応用例を紹介していきます。第 $1$ 回目は X-線変換の話をします。


X-線変換は、簡単に説明すると X-線 CT の原理となっている観測方法です。弧長表示されている直線 $l$ に対して

\begin{align}
\mathscr{X}f(l)=\int_{\mathbb{R}}f(l(s))\,ds
\end{align}

を $f$ の X-線変換と呼びます。直線に沿って $f$ を積分するということは、CT の世界では、どのくらいの X-線を吸収したのかという意味になるわけです。X-線 CT では、観測された X-線変換の像から $f$ を逆算する、すなわち反転公式が使われています。


良く知られた反転公式は投影定理を利用したものです。平面上の直線全体を

\begin{align}
l(\theta,h)(s)=se(\theta+\pi/2)+he(\theta)
\end{align}

のようにパラメトライズしておきます。X-線変換を

\begin{align}
\mathscr{X}f(\theta,h)=\int_{\mathbb{R}}f(se(\theta+\pi/2)+he(\theta))\,ds
\end{align}

と定義して $h$ に関するフーリエ変換を実行すると

\begin{align}
\mathscr{F}_h\mathscr{X}f(\theta,u)=\int_{\mathbb{R}^2}f(se(\theta+\pi/2)+he(\theta))e^{-2\pi ihu}\,ds\,dh
\end{align}

となります。ここで $(x,y)=se(\theta+\pi/2)+he(\theta)$ で変数変換すると

\begin{align}
\mathscr{F}_h\mathscr{X}f(\theta,u)=\int_{\mathbb{R}^2}f(x,y)e^{-2\pi i(x\cos\theta+y\sin\theta)u}\,ds\,dh
\end{align}

右辺は $f$ の $2$ 次元フーリエ変換なので

\begin{align}
\mathscr{F}_h\mathscr{X}f(\theta,u)=\mathscr{F}f(u\cos\theta,u\sin\theta)
\end{align}

と書くことができます。$f$ のフーリエ像が計算できたので、逆フーリエ変換すれば $f$ が求まります。以上が、良く知られた方法です。続いてヒルベルト変換を用いた別解を考えましょう。まず直線のパラメトライズを

\begin{align}
l(t,u,v)(s)=(s-t)e(u-\pi/2)+ve(u)+(x_0,y_0)
\end{align}

としておきます。これに合わせて X-線変換を

\begin{align}
\mathscr{X}f(t,u,v)=\int_{\mathbb{R}}f((s-t)e(u-\pi/2)+ve(u)+(x_0,y_0))\,ds
\end{align}

と定義しましょう。$v$ 方向への微分を考えると

\begin{align}
\partial_v\mathscr{X}f(t,u,0)&=\int_{\mathbb{R}}\partial_{e(u)}f((s-t)e(u-\pi/2)+(x_0,y_0))\,ds\\
&=\int_{\mathbb{R}}\dfrac{1}{s-t}\partial_uf((s-t)e(u-\pi/2)+(x_0,y_0))\,ds
\end{align}

右辺にヒルベルト変換が出てきました。後はヒルベルト像に関して解き、逆ヒルベルト変換を掛ければ $f$ が求まります。この計算は以前何回も取り上げたので省略します。今回はここまで。次回も引き続き応用例を紹介します。(written by J.Nishigata)


Copy Right 2020-現在 ALXIS DATA, Inc. All Right Reserved.

最小な反転公式フィルタの導出の続き