Deconvolution - 天文はかせ幕下

Teams勉強会
そもそもDeconvolutionとは
Pixinsightでの操作法、パラメターの意味
NGC2903の再処理

Teams勉強会

昨年末、ほしぞloveログのSamさんの呼びかけで、Teams上の画像処理座談会が行われたことがありました。テーマは「デジカメのノイズ」。もともとSamさんはノイズの専門家だそうですし参加者にはあぷらなーと氏、天リフの編集長、タカsiさん、Kagayaさん、kさんなど有名な方々だけでなく、有名会社でデジカメのセンサーを開発されている方とかが集まって相当にインパクトが強い「学会」でした。

この「学会」は一回限り。それではもったいないということで、ダボさん（関東で大学院生していて、天体写真をテーマにした研究をされている）の呼びかけて小さな画像処理研究会が組織されました。月に一回、満月の頃にTeamsのミーティングをやっています。

今月、その3回目のミーティングがありました。テーマは"deconvolution"。プレゼンターはブログ

を管理されているkさんです。お分かりの方なら「おお！」と思われたのではないでしょうか？21cmのニュートンと決して最新型ではない冷却CCDの組み合わせで、けっこうマイナーな銀河をしかも光害地から撮影されています。その画像はすさまじく、天文ガイドの入選もほぼ毎月という方です。最優秀賞も複数回受賞されているのではないでしょうか（確認してませんが）。

そのkさんが、銀河の"deconvolution"について話してくださるというのですから、顧問は3日くらい前から興奮していました。しかも私とそーなのかー氏が合作で撮影したNGC2903を例に処理の仕方やマスクの掛け方を話していただけました。

ちなみに下の写真が、私が一生懸命に処理したNGC2903です。

f:id:snct-astro:20210313100628j:plain

同じソースを使って、kさんが処理をおこなうと次のようになります。

f:id:snct-astro:20210313100632j:plain

・・・（無言

今回のエントリは、ミーティング得た内容をもとにして、上のイマイチなNGC2903を再処理した過程について、紹介したいと思います。

そもそもDeconvolutionとは

もし究極の望遠鏡とカメラがあったとしたら、それで撮影した星は1ピクセルの点として記録されるはずです。そのとき星雲や銀河はどんな姿なのか想像するだけでヨダレがでそうです。実際の天体写真では、星はぼんやりした楕円に写っています

f:id:snct-astro:20210313113406p:plain

こんな感じです。しかし，この星像には

「理想的には１ピクセルの点であるはずの星が、光学系と大気の揺らぎでどれだけぼやけるか？」

という貴重な情報が含まれているわけです。画像上での座標を $(x,y)$ とすれば，ボヤけた星を平均した輝度の分布は２変数の関数

$p(x-x_0,y-y_0)$

で表すことができます。これはPSF（Point Spread Function）と呼ばれていて、 $(x_0,y_0)$ にある点光源の星が $(x,y)$ の位置にどれだけの輝度を与えるか？つまり画像のボケ具合を表しています。

さて、究極のレンズとカメラがないと撮影できない「真の」天体像の輝度の分布を仮に $F(x,y)$ とし、顧問の普通のカメラとレンズで撮影した「ぼやけた」天体像の輝度の分布を $g(x,y)$ とします。それらの間には余計なノイズなどを度外視すれば

$g(x,y)=\displaystyle \int_0^{L_x}\int_0^{L_y} p(x-x_0,y-y_0)F(x_0,y_0)~dx_0dy_0~~~~~(1)$

という関係があります（この右辺は畳み込み(convolution)積分と名前がついています）。 $L_x,~L_y$ は画像の大きさを表しています。

いま $g(x,y)$ と $p(x,y)$ がわかっているので、原理的にはこの式を逆に辿れば $F(x,y)$ をもとめることができるはず、というのが"deconvolution（逆たたみこみ）"の考え方です。PSF関数を星の画像そのものから得られるこの方法は，天体写真にとって最適なシャープ化の手法といえるかもしれません。

(1)式から $F(x,y)$ を解く方法はいくつか提案されていて、中でも"Regularrized Richardson-Lucy" というのがその標準的な手法のようです。たとえば

の解説が比較的わかりやすいです。顧問も完全には理解できていないのですが、要は反復法とよばれる逐次近似計算によって、適当に仮定した近似解を，漸近的に本当の解へ収束させることをおこなっているようです。

Pixinsightでの操作法、パラメターの意味

Pixinsightには"deconvolution"の機能が搭載されています。

まず"Dynamic PSF"という機能で星の平均値からPSF関数を作った後，それを"deconvolution"に渡して計算を行います。その手続きはおなじみNiwaさんのブログ

にわかりやすく説明されていますので，ご覧ください。

いかではそのときに設定するパラメータの意味を考えてみます。下はDeconvolutionのウインドウです。右上の小さなウインドウは"Dynamic PSF"で得られた星像。ウインドウの最上段"PSF"のタブのところで，この星像を指定しています：

f:id:snct-astro:20210313123427p:plain

次に２段目の"Algorithm"です。

f:id:snct-astro:20210313123827p:plain

１行目では，(1)式を解くための複数の手法のうちいづれか（ここでは標準的なRichardson-Lucyを選んでます）を指定しています。つぎの"Iteration"は、上記の反復計算の回数ですね。基本的にはたくさん繰り返せばより真の $F(x,y)$ へ画像が近づいていくはずですが、画像にノイズなどが含まれていると、繰り返しを増やしすぎることによって、そのノイズが増幅されてしまったり、または数値振動を起こしたりと不具合も起こります（ちなみにkさんは100回を基本にしているとのことでした）。具体的な値は元画像との相談になるようですが，20~30回くらいがよいようです。

次に重要なのが、３段目の"Deringing"の設定です。どうしても生じてしまう星の周りに位落ち込みが生じるのを防ぐのが目的です。

f:id:snct-astro:20210313125307p:plain

local supportに星マスクを与えて、Local amountでその強度を調整するというのはNiwaさんのブログにあるとおりです。

Global BrightとGlobal Darkは、Deconvolution処理による輝度の修正の上限値と下限値を与えているもと思われます（これはちゃんとしたソースを見つけることができておらず、顧問の推測です）。つまりGlobal Dark=0.01とした場合、Deconvolutionの解の輝度が0.01以下になるとしたら、そこで強制カットオフして輝度を0.01以下にならないようにしていると予想しています。

kさんは、このGlobal Darkを基本的には０とし、最大限暗部を強調する手法をとられているようです。代償として生じるリンギングはマスク処理によって徹底的に排除されているようでした。書籍Inside Pixinsightや、Light Vortex Astronomy をみてもGlobal Darkを有限の値で調整してリンギングを防ぐと書いてあるので、これを0にしてしまう発想はありませんでした。ここが重要です。