天文はかせ三段目

天文はかせ三段目(仮)

仙台高専(旧・宮城高専)天文部の活動記録です

PIのPixelMathで,FAPの対数現像を真似する

前回のエントリ

に引き続き,今回は「対数現像」をPixelMathで真似してみます(「再現」と書かないのは,対数現像の中身がわからないため。本当に一致しているかの確証はないです)。ではまいりましょう

対数現像とは

荒井俊也氏が開発したFlat Aid Pro (FAP)は,通常のフラット補正では取りきれない光害カブリやフラットズレなどを補正することを目的としたソフトです。さらにガイドエラー補正や星のHDR処理(飽和復元合成)など有用な機能が詰め込まれています。800万画素以下なら無償というところも素晴らしいです*1

今回注目する「対数現像」はこのソフトに実装されています。レベル補正に輝度圧縮機能が付属していて,そこで「デジタル現像」と「対数現像」が選択できるようになっています。

しかし,対数現像の詳細は公表されていません。

対数現像は極度に明るい部分は飽和輝度を超えてしまうことがあるものの、星の芯が残りやすいという特徴があります。たとえ最大輝度部分がピンポイントで飽和しても画像として破たんするわけではな いので、星の表現のしやすさではデジタル現像より優れているのではないかというのが私の見解です。 Flat Aide Pro ver1 マニュアルより

この記述と「対数」というキーワードだけが手がかりになります。ここから推測していたわけですが,顧問はいろいろ試してみましたがあまりうまくいきませんでした。。。

Log Stretch

一方で,画像のストレッチに対数関数を用いる手法は,その他のソフトでも見受けられます。PixinsightではAutoHistgramというストレッチ機能内に"Logarithmic Transformation"というのがあります(その中身も不明)。さらに調べていると,天文計算のためのpythonパッケージを集めたサイト"The Astropy Project"”Logstretch”という関数があって,そこでは

The stretch is given by

 y=\displaystyle\frac{\log(ax+1)}{\log(a+1)}~~~~~~(1)

と書かれています。またTwitter上で@dabo_picさんからいただいた情報では,天文研究者の間で使われているSAOImageDS9でもほぼ同じ関数形のストレッチが使われているとのことです。この関数は次のようなトーンカーブと同じです(a=100)の場合

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

本エントリでは,上の関数が「対数現像」の正体であると仮定します。

Pixinsightで「対数現像」を真似する

さて,あらためてFAPの対数現像のウインドウをみてみます

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

FAPでも,対数現像はレベル補正にプラスアルファする形で実装されています。なので(1)式をそのまま適用するだけではうまくいきません。これはデジタル現像のときと同様ですね。そこで,処理前のj番目の輝度値を X_j,処理後のj番目の輝度値を Y_j,暗部側のスライダの輝度値をmとして,

 Y_j=\displaystyle\frac{\log\left\{ a(X_j-m)+1\right\} }{\log\left\{ a(1-m)+1\right\}}~~~~~~(2)

としてみます。また係数 aは強調の度合いを決めるパラメータで,これは明部側スライダの輝度値nからもとまると仮定します。その求め方は,デジタル現像の真似の時と同じようにやります。すなわち次の図の(c)

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

のように,X軸と接する部分での関数の傾きが,輝度圧縮をする前のレベル補正に対応するトーンカーブと一致するように aを決めるわけです。その計算は省略しますが,かんたんな微分の計算より aは方程式

 \displaystyle\frac{a}{\log(a(1-m)+1)}=\frac{1}{n-m}~~~~~~(3)

の解として定まることがわかります。しかし残念ながらこの式は、紙と鉛筆ではaについて解けません(解けないよね?)。

これをいちいち数値的に解くのも面倒なので、aの値は,処理結果をみながらなんとなく決めてもらったほうが早いと思います。以下の適用例でもM13については,そうしております*2

適用例

それでは適用結果を見てみます。PixelMathに入力する数式は次のようにします

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

コピペ用:iif( a*($T-m)+1 > 0, ln(a*($T-m)+1)/ln(a*(1-m)+1),0)
m=0.26, a=300

対数関数は負の値が代入されると不味いので、if文を使って処理しています。

まずはASI294MCで撮影した次のM13の画像に適用してみます。元画像は以下の通り:

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

これに対して、FAPで対数現像した結果は下のようになりました。

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

FAPにて対数現像。値はデフォルトのまま

これに対して、次がPixelMathによる真似です。ただしaの値はFAPの結果に近くなるように調整しています。

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

PixelMathにて,(2)式をa=750として適用

1パラメータ調整で、ほぼ似た結果を出力することができました。

ちなみに、比較のためにSIでのディジタル現像の結果を下に再掲載します。上の結果と比べると対数現像が確かに高輝度部の諧調に優れていることがわかります。

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

SIによるデジタル現像(比較用)

つぎに、デジカメによるのぎょしゃ座周辺の星野写真に、適用してみます。

f:id:snct-astro:20201014211102j:plain
上が元画像。

これに対してFAPでの対数現像を適用した結果がこちら: 

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

FAPによる対数現像:明部スライダ46600,暗部スライダ14600,輝度圧縮強さ160として元画像に適用

ストレッチ具合が弱いのは、明部スライダの値を大きくしすぎたせいだと思います(まあそれは良いでしょう)。これに対して、PixelMathによる真似は以下のようになりました;

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

PixelMathによる再現:m=14600/65535=0.22,a=20((3)式から推定)

違いがほとんどわかりません。

終わりに

というわけで、「対数現像」もリバースエンジニアリングに成功したと考えています。

PIのAutoHistgramの"Logarithmic Transformation"は、まさに対数現像であろうと考えています。ですのでPIユーザーなら、わざわざPixelMathを使う必要もないわけです。

最後に対数現像の「意味」について考えてみたいと思います。

Twitterにて,「対数現像わからんのおー」と呟いておったところ,ぐらすのすち同志に映像メディア情報学会誌の論文を教えてもらいました。その記事のイントロには「人間の明るさ感覚Rと入射光Lとの間にFechnerの法則

R\sim \log L

が成り立つことが知られている云々」との記述があります。一方で(1)式は、線形のグラフにlogを作用させただけです。つまりリニア画像に対して

  • 対数現像は、人間の網膜に入射する光量と、視覚の感度の間の非線形性をモデル化している

といえます*3。一方で、

  • デジタル現像は、フィルムに入射する光量と、現像液の感度の間の非線形性をモデル化している

わけです。どちらにも科学的な根拠があるわけですね。

 

さて、い か が で し た か?

今回のエントリ、役に立ちそうに見えて、実際の画像処理にはほとんど役立たないと思います。ただ、いろいろな処理の理解を深めることはできましたので、顧問は満足しています。思えば初心者の頃、ストレッチといえばGimpトーンカーブをいじることしか知りませんでした。紆余曲折を経て、またトーンカーブ処理の重要性を再認識した次第です。

 

 

*1:顧問はお金を払って,1年ライセンスを使っていました

*2:M13の結果はPIからの出力で、32bitです。これをそのままFAPで読むと16bitに変換されるのですが、その変換のされ方が不明なため、パラメータaを定めることができなかった。

*3:ただ、対数現像の結果がディスプレイに映し出された場合、それを見る人間の網膜の感度はlogになるから、我々が認識するのは\log\{log(ax+1)\}という光で・・・とかいうツッコミはなしで(無しでいいのか?)。