Python+OpenCVで流星群画像を仕分ける

この記事は約13分で読めます。

梅雨っぽくなってきました。曇りが続くと画像処理が捗りますね。

私(つのだ)は天文業界に足を踏み入れたばかりの若輩者のため、先人の皆様が天文ブログやTwitter、YouTubeなどで天文画像処理について取り上げられたものを頼りに勉強させていただいております。

今回はプライベートで撮影した流星群撮影画像を捌く内容をコードを交えてご紹介します。
12月の撮影後何ヶ月も放置した画像をブログを理由にようやく処理したのはよいものの、記事をこれまた何ヶ月も下書きのまま寝かし、締切を踏み倒し続けてとうとう半年近く経ってしまいました。

つのだ
つのだ

あ、ブログ管理担当のくぼたくん

くぼた
くぼた

……(圧)

つのだ
つのだ

シ…シメキリ、マモル、タイセツ……

くぼた
くぼた

(にっこり)

2020年ふたご座流星群

コロナ禍にあっても―

毎年友人とふたご座流星群を見に行くことが恒例行事になっています。

コロナ禍でしたが、屋外で一晩中空を見上げるだけなのでリスクは低そうと判断し、お互いにマスク等々安全に気をつけて決行しました。

自転車仲間が所属されている天文同好会の観望会にお邪魔させていただいていたのですが、今回は観望会自体は無いが個人的には見に行くということだったので、私たちも恒例となった場所で観望することにしました。

毎年、寒空の下で空を見上げながら一年間の振り返りのようなことをするのですが「今年はコロナで遠征が難しかった」とか「自宅だと曇りが長くほとんど稼働しなかった」等々コロナ禍での星撮りの難しさを改めて実感しました。

はじめての FUJIFILM カメラ

今回は11月に発売されたばかりのFUJIFILM X-S10を携えて撮影に臨みました。

X-S10は手ぶれ補正に惹かれて購入した、私にとって初のFUJIFILMカメラでしたが、インターフェイスが一般的なデジタルカメラに寄せられているためかほとんど違和感なく操作することができ、撮りたいものが素直に撮れ、さらにはX-Trans技術によってエッジがシャープな絵が撮れるという事で、私の中で最近お気に入りのカメラになっています。

流星群の撮影に向いている点としてはUSB給電ができるため、PowerDelivery対応のモバイルバッテリーから給電すると一晩余裕で持つという点です。寒空の下だとカメラバッテリーの消耗が激しく「しばらくして様子を見に行くといつの間にか電池が切れていた」ということはよくありますが、USB給電できるカメラだとその心配がありません。

また、カメラにUSB給電している間はスカイメモに電源を供給し続けることができるという思わぬ副産物もありました。通常、スカイメモ単体にモバイルバッテリーからUSB給電をすると電力消費が少なすぎるのか一定時間が経過するとバッテリーが勝手にスリープモードに入ってしまう問題があり、モバイルバッテリーの選定に悩んでいたのです。最近はX-S10、スカイメモ、モバイルバッテリー、三脚の4点セットだけ持って星撮りに行くスタイルが私の中で定番となっています。

ステライメージはFUJIFILM X-TransのRAWには未対応なのですが、流星群の撮影はJPEGしか撮らないので支障はありません。が、オリオン座付近を撮った絵を上司に見せたら赤いのがこんなに写るの?と驚いていたので、やっぱりRAW対応しないとなぁという思いを個人的に新たにしました。

オリオン座付近
FUJIFILM X-S10 / XF16mm F1.4 R WR
ISO 6400 / F2 / 6秒 / 2枚 / ステライメージ9で比較明合成

今回のふたご群の撮影は、諸事情があって撮影開始が0時前と遅かったのもあり全部で約1300枚。大した量ではないので全部見て片してしまうのが良いのかもしれませんが、徹夜明けにその体力がなく(しかも今回は翌日が平日でした)、後で後でと後回しにし続け、2ヶ月以上放置してしまいました。

今回のふたご群はどうにも短い流星が多く、細かい流星を目を凝らして探すのが億劫だった、というのもありました。

Python+OpenCVで画像処理

「怠惰であれ」というプログラマのモットーの下「OpenCVの直線検知で流星が含まれる画像を仕分ける」というテーマで1300枚を捌くことにしました。

記事中で紹介しているソースコードはGitHub Gistにまとめてありますので、ご参照ください。

Python 3.9を使用し、記事中で使用するライブラリは以下のようにimportしてあります(型ヒントを使用しているためそのまま動かす場合はPython 3.5以上が必須です)。

import math
import os
import sys
import typing
# 3rd-party libraries
import cv2
import numpy

直線検知で仕分け

まずは以下のload_gray関数で画像をグレイスケールで読み込み、detect_lines関数で画像中の直線を検出します。

def load_gray(filepath: str) -> numpy.array:
    """
    グレイスケール画像読み込み
    :param str filepath: 入力ファイルパス
    :return: グレイスケール画像データ
    """
    img = cv2.imread(filepath, cv2.IMREAD_GRAYSCALE)
    return img

def detect_lines(img: numpy.array) -> typing.List[numpy.array]:
    """
    直線検出
    :param numpy.array img: 入力画像(グレイスケール)
    :return: 検出直線リスト
    """
    ret, thr = cv2.threshold(img, 127, 255, cv2.ADAPTIVE_THRESH_MEAN_C)
    lines = cv2.HoughLinesP(thr, rho=1, theta=math.pi/180, threshold=200, minLineLength=20)
    return lines

これを走らせると、1300枚中400枚近く検知されてしまいました。

result = []
for filepath in tqdm(image_list):
    filepath = str(filepath)
    img = load_gray(filepath)
    lines = detect_lines(img)
    if lines is not None:
        result.append((filepath, lines))

print("detected:{}/{}", len(result), len(image_list))

面積を持った明るい領域の除去

画像を調べてみると、次のような原因によるものでした。

  • 流れてきた雲がスカイよりも明るく閾値を越えてしまう
  • 画像端に写り込んだライトなどが閾値を越えてしまう
雲がかかってしまった画像
左:元画像 中央:閾値処理 右:領域を塗りつぶし
カメラの前を横切ったことでライトが写り込んでしまった画像
左:元画像 中央:閾値処理 右:領域を塗りつぶし

そこで、画像から輪郭検出して面積を求め、一定面積よりも大きい領域をスカイの値で塗りつぶすことにしました。

薄雲の中を流れる流星などは検知できませんが、どのみち見ても残念な気持ちにしかならないと思うので気にしないことにします。

以下のdetect_area関数の引数thresholdは塗りつぶし判定の面積割合で、ここでは画像全体の0.01%よりも大きい面積を持つ輪郭を求めます。fill_area関数では先に求めた輪郭(contours)を含む外接矩形をbuffer_ratio分だけ拡張して画像の中央値(空が大部分を占めている画像なので画像の中央値≒スカイの値とみなす)で塗りつぶします。

ここで領域を若干拡張するのは、雲のように領域の境界付近がグラデーション状になっている画像で、稀にグラデーションが外接矩形をはみ出すことがあったため、念の為で広めに塗りつぶす意図があります。

def detect_area(img: numpy.array, threshold: float = 0.0001) -> typing.List[numpy.array]:
    """
    閾値を超える面積を持つ輪郭の検出
    :param numpy.array img: 入力画像
    :param float threshold: 閾値(画像全体の何%を`(0, 1]`で指定)
    :return: 閾値を超えた面積の領域リスト
    """
    height, width = img.shape
    img_area = width * height
    ret, thr = cv2.threshold(img, 127, 255, cv2.ADAPTIVE_THRESH_MEAN_C)
    contours, hierarchy = cv2.findContours(thr, 1, 2)
    contours = [cnt for cnt in contours if (cv2.contourArea(cnt) / img_area) > threshold]
    return contours

T = typing.TypeVar("T")
def clamp(v: T, min_v: T, max_v: T) -> T:
    """
    clamp関数
    入力値を`[min_v, max_v]`の範囲に収める
    :param T v: 入力値
    :param T min_v: 最小値
    :param T max_v: 最大値
    :return: `[min_v, max_v]`内に収めた値
    """
    return min(max_v, max(v, min_v))

def fill_area(img: numpy.array, contours: typing.List[numpy.array], buffer_ratio: float = 0.01, color: float = None) -> numpy.array:
    """
    領域の外接矩形で塗りつぶす
    :param numpy.array img: 入力画像
    :param contours: 領域リスト
    :param float buffer_ratio: バッファ率 (e.g. 6000 * 0.01 => 60px)
    :param float color: 塗りつぶしの色(未指定の場合は入力画像の中央値で塗りつぶす)
    :return: 塗りつぶし後の画像
    """
    height, width = img.shape
    x_buffer = int(width * buffer_ratio)
    y_buffer = int(height * buffer_ratio)
    # detect fill color
    if color is None:
        color = numpy.median(img)
    # fill bounding rect
    for cnt in contours:
        x, y, w, h = cv2.boundingRect(cnt)
        left = clamp(x - x_buffer, 0, width)
        top = clamp(y - y_buffer, 0, height)
        right = clamp(x + w + x_buffer, 0, width)
        bottom = clamp(y + h + y_buffer, 0, height)
        pts = numpy.asarray([
            [left, top],
            [left, bottom],
            [right, bottom],
            [right, top],
        ])
        img = cv2.fillPoly(img, pts=[pts], color=(color,))
    return img

流星検知処理の完成

以上の塗りつぶし処理を組み合わせて流星の写っている写真を検知するdetect_meteor関数ができました。

def line_length(line: numpy.array) -> float:
    """
    直線の長さの算出
    `cv2.HoughLinesP()`の返り値は`[ [ [start_x, start_y, end_x, end_y] ], ... ]`形式になっている
    :param line: `[ [start_x, start_y, end_x, end_y] ]`であること
    :return: 直線の長さ
    """
    assert len(line) == 1
    sx, sy, ex, ey = line[0]
    dx = ex - sx
    dy = ey - sy
    return math.sqrt(dx * dx + dy * dy)

def detect_meteor(filepath: str, area_threshold: float = 0.0001, line_threshold: float = 100) -> typing.Optional[typing.Tuple[str, typing.List[numpy.array]]]:
    """
    流星の検出
    :param str filepath: 入力画像ファイルパス
    :param float area_threshold: 面積のある領域検知用閾値(`detect_area()`関数`threshold`参照)
    :param float line_threshold: 検出した直線を流星と判定する最小の長さ
    :return: (画像ファイルパス, 検出した直線) or None
    """
    img = load_gray(filepath)
    area_contours = detect_area(img, area_threshold)
    if area_contours:
        img = fill_area(img, area_contours)
    lines = detect_lines(img)
    if lines is not None:
        length = max([line_length(x) for x in lines])
        if length > line_threshold:
            return lines
    return None

処理を走らせると、1300枚中10枚の画像に流星が写っていると検知されました。画像を一つ一つ確認してみると全ての画像に流星が写っています。良さそうです。

result = []
for filepath in tqdm(image_list):
    filepath = str(filepath) # fqdm -> str
    lines = detect_meteor(filepath)
    if lines is not None:
        result.append((filepath, lines))

print("detected:{}/{}", len(result), len(image_list))
検出画像合成結果(4枚中に5本の流星)
FUJIFILM X-S10 / XF16mm F1.4 R WR
ISO 6400 / F2 / 6秒 / 4枚 / ステライメージ9で比較明合成

画像処理はMacBook Air (M1, 2020) 16GB上で行い、6240×4160サイズのJPEG画像 1300枚を処理するのにかかった時間は約5分でした。画像をAmazon Prime Photosにアップしつつバックグラウンドで処理を走らせればアップロードが完了するまでには選別が終わりそうですね。

行き当たりばったりで作った検出処理なため、別の環境・条件での撮影画像ではパラメータの調整をしないとうまくマッチしないかもしれません。もっといい手法等があれば、是非ブログやつぶやき等でご指摘・ご反応いただけると幸いです。

関連情報

タイトルとURLをコピーしました