あいさつ

久々のブログ更新になります．
自作ソフトウェアに PixInsight (以下PI) における STF, HT, DBE 相当の機能を追加しました．

STF

ScreenTransferFunction（STF）は，PI をお使いの方であればお馴染みかと思います．
画像のデータではなく，伝達関数（Transfer Function，データを色や透明度などの(主に)視覚情報へ伝達するための関数）を変更することで，画像の色味などを変えるプロセスです．
本来伝達関数は大変奥の深い概念ではありますが，普通の RGB 画像やモノクロ画像でやることは基本的にシャドウ，中間調，ハイライトの3つのパラメタ値を変更し，それら値を元にできる Midtone Transfer Function (MTF) を変更することになります．
MTF は PI のドキュメントや蒼月城さんの解説動画にあるように，データの値が中間調と等しい時にちょうど0.5となるようにいい感じにストレッチする関数で，次のような式で計算されます．

$$ \rm{MTF}(\it{x}\rm{;} \it{m}\rm{)} = \begin{cases}
0 \hspace{40px} (x = 0) \
\frac{1}{2} \hspace{40px} (x = m)\
1 \hspace{40px} (x = 1)\
\frac{(m-1)x}{(2m - 1)x - m} \hspace{40px} (\rm{otherwise})
\end{cases} $$
PI のドキュメントより

というわけで，現在のソフトに RGB それぞれのカラーバーを追加，中間調の変更機能も追加して，MTF で変換するようにしました．
また，PI の STF プロセスにはいい感じに画像を見やすくオートストレッチする機能もあります．中で何をやっているのかは蒼月城さんが Twitter に書いてくれていました．

これと，PCL のソースコードを参考にすれば，ぶっちゃけ実装することは可能になります．

できました．M42 が一気に見やすくなったところは，オートストレッチをかけているところになります．

HT

Histogram Transformation（HT） は STF と密接に関係しています．STF では伝達関数を変更しましたが，HT では STF での結果になるようにデータの方を変更します．基本的な考え方は STF での伝達関数の適用と同じと言っていいです．
PI ではもう少しパラメタがありますが，今回はそれらは無視して，RGB チャネルそれぞれのシャドウ，中間調，ハイライトの合計9つのパラメタ値をもとにデータを変更します．

完成形はこのようになります．左下のウィンドウでは先ほどと同様に STF で伝達関数を変更するのみですが，右下のウィンドウでは実際のデータのほうが STF の変更に沿うように変換されています．右下のウィンドウのヒストグラムがちゃんと変更されているところに注目していただければ，データが変わっているということがおわかりかと思います．ユーザは9つのパラメタ値と連動する変数を上のビジュアルプログラム上に追加し，それを HT のノードに渡すようにすることで，このような処理を実現します．
STF と HT については変更箇所が多いのでソースコードを載せることはしません．

DBE

DynamicBackgroundEstimation (DBE) は STF や HT とは毛色が違うものの，よく使われるプロセスです．以前 takashi さんが Python で DBE するプログラムを作られていたので，似た感じの処理ができるものを自分も作ることにしました．
本質的にやることは比較的単純で，

  1. ユーザ指定の点群から付近の20*20くらいのボックスを切り出してボックス内の中央値を計算する
  2. 中央値に合うように画像を曲面フィッティングしてバックグラウンドモデルとする
  3. 背景を元画像から引く（または割る）

となります．
フィッティングはシンプルに最小二乗法を使っています．こちらのブログが詳しいです．
以下のようなコードで実装しました．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
dynamic_background_extraction<IOValue, IOErr>(
    image: Image, points: Roi,
    function_degree: Integer = 1
) -> Image, Image {
    let tag = image.tag();
    let dim = image.scalar().dim();
    let dim = dim.as_array_view();
    let func_plus_one = (*function_degree + 1) as usize;
    let params_right_num = func_plus_one * func_plus_one;
    let mut params_right = Vec::<f32>::with_capacity(params_right_num * dim[0]);
    params_right.resize(params_right_num * dim[0], 0.0);
    let mut params_right = Array::from_vec(params_right).into_shape(
      (dim[0], func_plus_one, func_plus_one)).unwrap();
    let params_left_num = params_right_num * params_right_num;
    let mut params_left = Vec::<f32>::with_capacity(params_left_num * dim[0]);
    params_left.resize(params_left_num * dim[0], 0.0);
    let mut params_left = Array::from_vec(params_left).into_shape(
      (dim[0], func_plus_one, func_plus_one, func_plus_one, func_plus_one)).unwrap();
    let mut model = image.scalar().clone();
    let mut out = image.scalar().clone();
    for i in 0..dim[0] {
        for takes_y in 0..func_plus_one {
            for takes_x in 0..func_plus_one {
                for inner_takes_y in 0..func_plus_one {
                    for inner_takes_x in 0..func_plus_one {
                        let data = points.filterx(image.scalar().slice(s![i, .., ..]));
                        for ((x, y), _) in data {
                            let x = x as f32 / dim[2] as f32;
                            let y = y as f32 / dim[1] as f32;
                            params_left[[i, takes_y, takes_x, 
                              inner_takes_y, inner_takes_x]] += 
                              1.0 
                              * (x as f32).powf((takes_x + inner_takes_x) as f32)
                              * (y as f32).powf((takes_y + inner_takes_y) as f32);
                        }
                    }
                }
                let data = points.filterx(image.scalar().slice(s![i, .., ..]));
                for ((x, y), _) in data {
                    let mut boxdata = Vec::new();
                    for xx in (x-10).max(0)..(x+10).min(dim[2]) {
                        for yy in (y-10).max(0)..(y+10).min(dim[1]) {
                            boxdata.push(model[[i, yy, xx]]);
                        }
                    }
                    boxdata.sort_by(|a, b| a.partial_cmp(b).unwrap());
                    let mid = boxdata.len() / 2;
                    let med = if boxdata.len() % 2 == 0 {
                        (boxdata[mid] + boxdata[mid - 1]) / 2.0
                    } else {
                        boxdata[mid]
                    };
                    boxdata.clear();
                    let x = x as f32 / dim[2] as f32;
                    let y = y as f32 / dim[1] as f32;
                    params_right[[i, takes_y, takes_x]] += 
                      med 
                      * (x as f32).powf(takes_x as f32)
                      * (y as f32).powf(takes_y as f32);
                }
            }
        }
        let params_left_vec = params_left.index_axis(Axis(0), i)
                                          .to_owned().into_raw_vec();
        let params_right_vec = params_right.index_axis(Axis(0), i)
                                          .to_owned().into_raw_vec();
        let a = DMatrix::from_vec(params_right_num, params_right_num, params_left_vec);
        let b = DVector::from(params_right_vec.clone());
        let decomp = a.lu();
        let x = decomp.solve(&b);
        let mut sol = DVector::from(params_right_vec);
        match x {
            Some(vector) => sol = vector,
            None => {
                println!("failure");
            },
        };
        for y in 0..dim[1] {
            for x in 0..dim[2] {
                let xx = x as f32 / dim[2] as f32;
                let yy = y as f32 / dim[1] as f32;
                model[[i, y, x]] = 0.0;
                for takes_x in 0..func_plus_one {
                    for takes_y in 0..func_plus_one {
                        model[[i, y, x]] += 
                          sol[takes_y * func_plus_one + takes_x]
                          * (xx.powf(takes_x as f32))
                          * yy.powf(takes_y as f32);
                    }
                }
                out[[i, y, x]] -= model[[i, y, x]];
            }
        }
    }
                    
    vec![Ok(IOValue::Image(WcsArray::from_array_and_tag(Dimensioned::new(
        out,
        Unit::None,
    ), tag.clone()))), 
        Ok(IOValue::Image(WcsArray::from_array_and_tag(Dimensioned::new(
        model,
        Unit::None,
    ), tag.clone())))]
}

実際に動かしてみると以下のようになります．

DBE

画像は takashi さんのデータを無許可でお借りしています．ごめんなさいm(_ _)m（Github で MIT ライセンスで公開されているデータなのできっと許してくれるはず）まだビジュアルプログラムのノードにオーサーシップを明示することはしていないのですが，その暁には彼の名前は必ず含めるつもりです．他にもテスト用の、被りがあるデータを提供してくださる方を募集しております。^^
上の例では点群は takashi さんのものと同様で，二次曲面でフィッティングしています．左下が元画像，右上が推定されたバックグラウンドモデル，右下がモデルを引き算した結果です．元画像よりも強調していますが，明るさの勾配は除去できているのではないでしょうか．

まとめ

さて，今回は後処理における三種類のプロセスを作ってみました．まだまだデータを仕上げるには必要なプロセスがたくさんありますが，その中身についてちまちま勉強して追実装するのはなかなか楽しいので今後も続けていきたいと思います．