お正月は満月期。２０２１年は画像処理のマニアックな話題からスタートしたい所存です。
はじめに
テーマはタイトルの通り。天体写真の画像強調（ストレッチ）において、ノイズ処理をどのタイミングで行うのが正解か？という問題です。
これについてはおそらく３つの「派閥」が存在すると思われます。
- ストレッチの前にやるよ派
ストレッチによってノイズもまた強調されてしまう。だから強調後にノイズを低減させるよりも、ノイズが小さいうちにこれを叩いておいて、そのあと思う存分にストレッチしたほうが良いはず！ - ストレッチの後にやるよ派
ノイズ処理というのはある種の平均化処理であり、その実行によって必ず何らかの構造が失われてしまう。だから、ストレッチによってまずは星雲・分子雲の構造を引き出し、それをはっきりと把握したのちに、注意を払いながらノイズ処理を行うのが肝要である！ - ストレッチの前と後、両方やるよ派
ノイズと星雲・分子雲の構造、それぞれの画像上の長さ（ピクセル量）に隔たりがあるなら、（巧くやれば）ノイズ処理によって構造が失われる心配は小さい。よって前後に一回づつノイズ処理を行うことで、「前にやるよ派」と「後にやるよ派」のいいとこ取りができる！
なるほど。三者三様に理があるように思えます。
しかしながら、これらを全て同じまな板の上に載せて比較・検証を行うことは容易ではありません。なぜなら一言にノイズ処理と言っても、単純なボカシからwaveletなどの構造解析を取り入れたノイズ処理、さらにはAIを使ったものまで実に様々だからです。同じことは強調処理についても言えます。さらにはノイズにも、輝度ノイズ・カラーノイズ・アンプノイズなどあり、何に着目するかで結果が変わってくるでしょう。
そこでこのエントリでは，もっとも典型的でかつ単純化した状況において上の三つを比較して結果を検証してみます。次にその具体的な方法を説明します。
検証の方法
基本的な考え方は，上の３つの派閥の画像について，ストレッチの強度を共通にしておいて，背景のノイズレベルがほぼ同一になるようにノイズ処理を加えた画像について，明部の構造の様子を比較する，という考え方です。
ストレッチとしてはレベル補正の中間・暗部スライダを動かすだけの最も単純な方法と，wavelet を使った構造強調の組み合わせに限定します。またノイズ処理は，これもwaveletを基にして，長さスケールごとにボカシをかけるおそらく最も典型的な方法だけを想定します（Photoshop 内 Camera Rawの輝度ノイズ低減など）。
以下では，その内容をもう少し詳しく書いておきます。
まず検証の対象にするのは次の画像です
これはデジカメ(EOS6D)で撮影し、Pixinsightで前処理（フラット・バイアス補正、スタック）を行なった直後の画像の一部を切り出したものです（Pixinsightは、もともと14bitのデジカメ画像の上側にデータを追加して32bitに伸張するので、出力直後はこのように真っ暗な画像になります。)
画像処理では、この画像から出発して強調処理とノイズ処理を行なって天体写真を完成させるわけです。
何を使って強調・ノイズ処理をおこなうか？
強調処理：STFのAutoStretchとHDRMT
天体写真の強調処理は、コントラストの強調／構造の強調の２つに分類できます（他に彩度の強調がありますが、今回は輝度ノイズだけに着目しているので、それは省きます）。
コントラスト強調については、STF（Screen Transfer Function ）のAutoStretchを用いることにします。これは、Pixinsightの便利な機能の一つで、画像の情報をもとに一定の強度の強調を自動的に行ってくれる機能です。例えば、上の暗い画像にAutostretchを施すと、次の画像のようになります。
この中身はレベル補正を自動的に行なっているだけです。上の結果の場合、中間スライダと暗部側のスライダを、次の図のようにのように調整しています。
次に構造強調にはPixisightのHDRMT（High Dynamic Range Multiscale Transform）を以下の設定で施します。
HDRMTは、wavelet処理を用いてある長さで変化している星雲などの構造に対して、そのコントラストを部分的に強調する手続き（のようです）。Photoshopでの「明瞭度」やハイパス強調と中身はだいたい同じであると推測します。うえのSTFを施した画像に、さらにHDRMTを行った画像がこちら：
星雲の明暗の変化が激しい部分に注目すると、変化がよくわかると思います。
ノイズ処理：MLT (Multiscale Linear Transform)
MLT (Multiscale Linear Transform)はwavelet処理を元にしたPixinsightの非常に有用な機能で，シャープ化やノイズ低減、その他いろいろなことが出来ます。書いている本人も全体は把握できていないのですが、ここではMLTのノイズリダクション機能を利用します。
1, 2, 3, 4, RとあるLayerはwavelet処理を行う長さスケールを表していて、それぞれ画像上の1, 2, 4, 8, 16pxの長さに対応しています。"Amount”はノイズ処理の強さを表し、"Threshold"はノイズ処理を施すピクセルの変化のしきい値を表しています。
いまは輝度ノイズ（ショットノイズ）の低減を考えているので、1pxの長さの輝度変化のみを選択的に低減させればいいはずなので、Layer１のみを選択し"Threshold"は3.0、"Amount"は適宜調整することにします。
次は"Amount=1"としてMLTのノイズ処理前と、処理後の画像を並べたものです。
どうやって比較するか
まずストレッチの前後にノイズ処理を行う「両方やるよ派」の画像を用意します。つまり
両方やるよ派：元画像 → ノイズ処理（MLT,Amount=0.6）→ ストレッチ（STF+HDRMT） → ノイズ処理（MLT,Amount=0.6）
（"Amount" 0.6というのは顧問がいつもMLTを行うときに選んでいる値で、特に意味はありません）。下の写真はそれによって得られた画像です。
この画像を比較対象に用います。ノイズの強さを定量化するために，中央に四角で示した背景部分の輝度について、pixel毎の輝度値の分散を計算し、その値を「ノイズ強度」とみな巣ことにしました（その方法については後述します）。
この「両方やるよ派」の画像の場合、背景の分散はでした（輝度は0を黒、１を白としてます。もともと輝度の揺らぎが小さいので、分散は極端に小さい値になります）。
つぎに「前にやるよ派」と「後にやるよ派」の画像を用意します。つまり
前にやるよ派：元画像 → ノイズ処理（MLT,Amount = X ） → ストレッチ（STF+HDRMT）
後にやるよ派：元画像 → ストレッチ（STF+HDRMT） → ノイズ処理（MLT,Amount= Y ）
と処理を行なった画像の背景のノイズ強度が、「両方やるよ派」に一致するように、それぞれの"Amount" X,Yの値を調整するわけです。ただしSTFとHDRMTのパラメータは全ての画像で共通です。（ちなみにこの画像の場合 X=1.0, Y=0.95で ノイズ強度が「両方やるよ派」に一致しました。）
以上のようにして用意した３つの画像について、右下の四角部分で示した「星雲のコントラスト高い部分」を比較しよう。というわけです。
このようにすれば、一応は同じまな板の上で、それそれの派閥の処理画像結果を比較できると考えます。ただし「星雲のコントラスト高い部分」については、その点数を数値化する良い方法が浮かびませんでしたので、ここでは見た目で判断する事にします。
結果
それぞれの結果を示します。あらかじめお断りしますと、違いはかなり微妙です。
前にやるよ派↓
後にやるよ派↓
両方やるよ派↓
ちょっとわかりにくいので、さらに拡大して部分ごと比べてみます
ここからはどうしても私感になりますが，
前にやるよ派 ＞ 両方やるよ派 ＞ 後にやるよ派
の順で星雲の構造がはっきりしていると見えます。「後にやるよ派」や「両方やるよ派」では，構造強調の後にノイズ処理が入るので，どうしても星雲の輪郭がボケてしまうようです。ただしその差は，ほんよにごくわずかです。
まとめると，
強調処理の強度が同一で、かつ処理後の背景ノイズの強度を揃えるという条件の下では、ノイズ処理は強調処理の前に行ったほうわずかに結果がよい。
となりますね。
実際の画像処理ではマスクを利用して選択的にノイズ処理を行ったりシャープ化をしたりするので、上のわずかな差はほとんど見えなくなってしまうだろうと思われます。しかし最近Zoom会議などを通していろんな方とお話しした範囲では、ノイズ処理をストレッチの前に行うことを避けているという方が多かったです。waveletのlayerを適切に選びさえすれば，ノイズ処理によって構造が失われる心配は小さく，ストレッチ前にノイズ処理を行う方法も捨てたものではないと，言えるかと思います。
おわりに
というわけで、新年早々かなり時間をかけて検証した割には、結論がしょぼくなってしまいました。
最後に、画像の「ノイズ強度」を計算する方法については、Rubyからrmagickというgemをつつかいました。その導入にあたって twitter上で@rnaさん（このブログでよく引用している
Deep Sky Memories の中の方）や@watsonさんにアドバイスをいただきました。お礼申し上げます。
自分用メモとして、その導入のための概要を書いておきます（環境はMacOS Catalina）
- homebrew からrubyをver.3.0.0p0 指定でインストール（普通に"brew install ruby"とやるとver.2.x.xが入るみたい。するとあとのrmagickのインストールでハマった）
- 同じくhomebrewからimagemagickをインストール。(webの情報ではImageMagick v7だとrmagickが対応していないので、v6以前をインストールする、とあるがこれは既に解決済みでrmagickはv7にも対応しているみたい 2021.1.3）
- ”gem install rmagick”としてrmagickをインストール（gemというのはrubyをインストールすると使えるようになる）
- ソースコードは
の下の方にあるコードを参考にした。
また解析する画像がfitsなら、pythonなどを使うのも簡単で良いみたいです。
astropyも使えるようになっておいたほうが、いろいろ便利そう。
