よしなにスターエンハンスしてみる

目次

  1. 1. あいさつ
  2. 2. モチベーション
  3. 3. 手法
  4. 4. プログラム
  5. 5. 結果
  6. 6. 得られた知見と感想

あいさつ

こんにちは．先日行ったあいあい岬で撮影したアンタレス付近を編集していて，こんなことをやってみましたので，その紹介記事です．

モチベーション

Twitterでフォローしているnagahiroさんのこちらの記事を見て，初めて”Ringing”という名前がついていることを知ったのですが，これは天体写真を編集しているとよくある現象です．例えば，Ps，Lrの「明るさの最小値」やSI8の「スターシャープ」をかけて星雲を強調していると，微光星の周りが暗く落ち込んでしまうということがよく起きます．画面上で見たり小さい写真用紙に印刷する時はあまり気になりませんが，A4写真用紙等に印刷すると主に星雲の中でかなり気になるレベルで影響が出てきます（僕の場合だけかもしれませんが笑）．これをなんとかしたいというのが今回のモチベーションになります．

手法

nagahiroさんはこれをハイパスを使って，かなりうまく補正していますね．私はそういった前提知識がなかったので，別の方法を考えてみました．具体的には以下の通りです．

  • 星雲の中の微光星を見た時，その輝度や彩度は内側から「微光星（めっちゃ明るい，彩度低）」→「Ringing（暗い，彩度低）」→「星雲（微光星ほどではないが明るく，かつ彩度が高い）」のようになっている，という仮定をおきます．
  • よって，テキトーなフィルタを使って周辺の画素を使って埋めてやる（何回か繰り返す）．
  • 具体的にはRingingが起こっている画素を中心とした5×5のフィルタで，彩度が高い画素トップ9くらいを選んで，その平均をとってやる．

    例えばこんな感じの5×5の並びがあった時に，Ringing Pixelは星雲あたりの黄色っぽい画素を使って置き換えてやればいい感じになるんじゃないかな？というアイデアです．

とまあこんな感じになります．それでは実装していきます．

プログラム

実装にはC++（Visual Studio 2017に入ってるコンパイラ）とOpenCVを使います．OpenCVの導入については「OpenCV 4.1.0をVisual Studio 2019から使用する時の手順」あたりを参考にどうぞ．

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
#include "stdafx.h"
#include <opencv2/opencv.hpp>
#include <iostream>
#include <string>
#include <sstream>
#include <algorithm>

using namespace std;
using namespace cv;

#define INF 10000000

int main(int argc, char* argv[])
{
    string filename = "input.tif";
    Mat image = imread(filename, 2 | 4);
    
    //画像読込み
    if (image.data == NULL) { 
        cerr << "Error: Couldn't read file " << filename << endl;
        return -1;
    }
	
    int xsize = image.cols;
    int ysize = image.rows;
    int th = 10000; //輝度のスレッショルド

    Mat output = Mat::zeros(ysize, xsize, CV_16UC3); //16bit tifの場合こう
    double sat[5][5] = { 0.0 }; //彩度
    int rank[5][5] = { 0 }; //彩度のランク(大きいほど明るい)
    ushort maxc, minc;
    maxc = -INF, minc = INF;

    for (int t = 0; t < 30; t++) { //30回操作
        for (int j = 0; j < ysize; j++) {
            for (int i = 0; i < xsize; i++) {
                if (image.at<Vec3w>(j, i)[0] < th  && image.at<Vec3w>(j, i)[1] < th && image.at<Vec3w>(j, i)[2] < th && j > 2 && i > 2 && j < ysize - 2 && i < xsize - 2) { //Ringingが起きているなら
                    int lc, mc;
                    for (lc = 0; lc < 5; lc++) {
                        for (mc = 0; mc < 5; mc++) {
                            rank[lc][mc] = 0;
                        }
                    }
                    lc = mc = 0;
                    // 周辺画素の彩度を計算
                    for (int l = j - 2; l <= j + 2; l++) {
                        for (int m = i - 2; m <= i + 2; m++) {
                            maxc = max(image.at<Vec3w>(l, m)[0], image.at<Vec3w>(l, m)[1]);
                            minc = min(image.at<Vec3w>(l, m)[0], image.at<Vec3w>(l, m)[1]);
                            maxc = max(maxc, image.at<Vec3w>(l, m)[2]);
                            minc = min(minc, image.at<Vec3w>(l, m)[2]);
                            sat[lc][mc] = maxc == 0 ? 0 : ((double)(maxc - minc) / (double)maxc);
                            mc++;
                        }
                        lc++;
                        mc = 0;
                        maxc = -INF, minc = INF;
                    }
                    // 彩度でランキング
                    for (lc = 0; lc < 5; lc++) {
                        for (mc = 0; mc < 5; mc++) {
                            for (int lcc = 0; lcc < 5; lcc++) {
                                for (int mcc = 0; mcc < 5; mcc++) {
                                    if (!(lc == lcc && mc == mcc) && sat[lc][mc] > sat[lcc][mcc]) {
                                        rank[lc][mc]++;
                                    }
                                }
                            }
                        }
                    }

                    lc = mc = 0;
                    // フィルタを使って，周辺画素から置換
                    for (int l = j - 2; l <= j + 2; l++) {
                        for (int m = i - 2; m <= i + 2; m++) {
                            if (rank[lc][mc] >= 16) {
                                for (int k = 0; k < 3; k++) {
                                    output.at<Vec3w>(j, i)[k] += image.at<Vec3w>(l, m)[k] / 9;
                                }
                            }
                            mc++;
                        }
                        mc = 0;
                        lc++;
                    }
                }else { //Ringingが起きていなければ無視
                    output.at<Vec3w>(j, i)[0] = image.at<Vec3w>(j, i)[0];
                    output.at<Vec3w>(j, i)[1] = image.at<Vec3w>(j, i)[1];
                    output.at<Vec3w>(j, i)[2] = image.at<Vec3w>(j, i)[2];
                }
            }
        }
    }
    imwrite("output.tif", output);
    return 0;
}

結果

結果は冒頭のツイートの通りで，落込みの完全解消までには至りませんでしたが，SI8のスターエンハンスに比べ，微光星の主張が激しすぎず，かつ周辺の彩度を保ったまま落込みを補完することができたんじゃないかなあと思います．画面上ではあまり違いが分かりませんが，この後少し微調整して印刷したところ，Ringingに見られるぼつぼつがある程度解消され，かなり有効だということがわかりました．

処理前

処理後

いかがでしょうか，黒く落ちすぎている部分が，周りの星雲の色でもって塗られていることがわかりますね．

得られた知見と感想

プログラムについてですが，16bit tiff形式を読書きするやり方がわからず少し時間を使いました．imreadの引数や行列の定義はプログラムの通りですし，画素へのアクセスはimage.at<Vec3w>(y, x)[k]</code>でいけるっぽいです．
上記のプログラムで何回か(30回)フィルタをかませていますが，繰り返してもあまり意味がありませんでした．それと，今回は輝度だけでRingingが起きている画素かを判断していて，これはよくないですね．例えば暗黒星雲をRingingが起きている画素と判断し，暗黒星雲のある部分だけ明るくなってしまうということが起きました（以下画像参照）．

ハイパスを使ったやり方も楽しそうなので，原理を勉強しつつ今後実装したいと考えています．