霧除去（Haze Reduction）とは

霧除去はその名の通り、霧がかった画像から白い靄を取り除く画像処理技術です。
霧除去のアルゴリズムは様々なのものが提案されていますが、中でも最も有名なものはHe.Kの手法(Dark Channel法)※だと思います。
この手法は、画像中の霧がかった領域が白っぽくなる特性を利用し、白くない度合い（DarkChannel）を定義しました。
そして、Dark Channelを用いて局所的に補正することで、霧の濃淡の影響を受けない処理を実現しました。

今回はこのDark Channelを用いた霧除去アルゴリズムをPythonで実装していきます。

※He, Kaiming, Jian Sun, and Xiaoou Tang. "Single image haze removal using dark channel prior." IEEE transactions on pattern analysis and machine intelligence 33.12 (2010): 2341-2353.
He KaimingはCVの研究者なら一度は名前を聞いたことがあると思います。Guided FilterやResNetなどを提案しており、CV業界に影響を与えています。確かCVPR2019でベストリサーチャー賞を受賞されてた覚えがあります。

Dark Channelによる霧除去アルゴリズム

多くの霧除去アルゴリズムでは、霧画像の撮像モデルを用いています。
霧が含まれた撮像画像をI、霧が含まれない目的画像をJとした時の撮像モデルは以下になります。

$\displaystyle I(x, y) = J(x, y)t(x, y) + A(1-t(x, y))$

ここで、Aは大気光（推定値なので後で説明します）、tは透過マップ（こっちも後で推定します）になります。
目的画像はJなので、観測画像Iとそれ以外の値が定まれば求める事ができます。

Dark Channelの推定

Dark Channelは以下で定義されます。

$\displaystyle DC(x, y) = \min_{(u, v) \in \Omega(x, y) } \min_{c \in {R, G, B}} I^c(x, y)$

Ωはパッチ範囲です。
この式を言葉で説明すると、 1. 入力画像の各ピクセル毎にR,G,Bの中で最も値の小さいものを選択する 2. パッチ領域（例えば3x3）毎に最も小さい値を選択するという処理になります。
これをPythonで実装すると以下のようになります。

import cv2
import numpy as np
from tensorflow.keras.layers import MaxPooling2D


def MinFilter(img, ksize=3):
    I = 255 - img
    O = MaxPooling2D(ksize, strides=1, padding="same")(I[None,:,:,None])
    O = O.numpy()[0,:,:,0]
    return (255 - O).astype(np.uint8)


if __name__ == "__main__":
    img = cv2.imread("image.jpg")
    DC = np.min(img, axis=-1)  # min_{R, G, B}に相当する処理
    DC = MinFilter(DC) # min_{Ω}に相当する処理

この実装のポイントはmin_{Ω}の部分です。
普通に実装しようとすると２重（もしくは４重）のforループでSliding Window処理を行う必要があります。
しかしながら、numpy arrayにランダムアクセスすると処理がとても重くなってしまいます。
そこで、tensorflowやpytorchのMaxPooling関数を利用し、MinPooling処理を実装しました。

実行例を以下に示します。

入力画像	Dark Channel画像

Dark Channel画像では、輝度値の高い部分（白い所）が霧が濃い場所となります。

大気光Aの推定

大気光は霧領域の中で最も明るい領域と定義されています。具体的な処理としては、まずDark Channel画像の中から輝度値が上位0.1%のピクセルを選択します。
その後、その座標に相当する入力画像のピクセルを選択し、その中で最も高い輝度値を大気光Aとします。

この処理をPythonで実装すると以下のようになります。

def estimate_A(img, DC):
    index = np.where(DC > np.percentile(DC, 99.9))
    A = np.max(img[index])
    return A

np.percentile(arr, N)は、numpy arrayの中から、上位N%の要素を返すメソッドです。とても便利ですが今回初めてその存在に気が付きました。

透過マップの推定

透過マップはDark Channel画像と大気光から求める事ができます。

$\displaystyle t(x, y) \approx 1 - \frac{DC(x, y)}{A}$

w = 0.5
t = 1 - w*(DC/A)

wは補正の強さを決めるハイパーパラメータであり、0から1の実数値を取ります。

目的画像の算出

透過マップ、大気光が求まったので、いよいよ目的画像を推定します。
目的画像は以下の式で求まります。

$\displaystyle J(x, y) = \frac{I(x, y) - A}{\max(t(x, y), t_0)} + A$

t0はゼロ除算を抑制する調整項であり、論文では0.1が設定されています。
画像処理ではこういった除算処理が登場しますが、ゼロ除算処理が曖昧に記述されているケースがよくあります。
本論文は、実装しやすいように記述されており、とても助かりました。

pythonで実装すると以下のようになります。

t[t < 0.1] = 0.1
J = (img - A)/t[:,:,None] + A

実行結果

霧除去アルゴリズムを実行した結果を示します。

入力画像	霧除去画像

横並びにすると改善効果がパッと分からないですね。
なので、連続画像にしてみました。

画像

２枚目と３枚目は改善効果が確認できますね。
多分ハイパーパラメータの設定（Kernel Size、重み）次第ではもう少し改善すると思うので、いい感じのパラメータを探してみたいです。

まとめ

Kaiming先生のHaze Reductionを実装しました。
非学習ベース（フィルタリング処理）で実装されている為、実装してすぐに試せるのは非常に嬉しかったです。
多分、このアルゴリズムに刺さるデータがあると思うので（ダイナミックレンジが深いとか）、使いどころを考えていきたいです。
最後に全体の実装を載せておきます。

全体の実装

import cv2
import numpy as np
from tensorflow.keras.layers import MaxPooling2D
import sys



def MinFilter(img, ksize=3):
    I = 255 - img
    O = MaxPooling2D(ksize, strides=1, padding="same")(I[None,:,:,None])
    O = O.numpy()[0,:,:,0]
    return (255 - O).astype(np.uint8)


def estimate_A(img, DC):
    index = np.where(DC > np.percentile(DC, 99.9))
    A = np.max(img[index])
    return A
    

def reduction(img, w=0.5):
    img = img.astype(np.int64)
    DC = np.min(img, axis=-1)
    #DC = np.max(img, axis=-1)
    DC = MinFilter(DC, 7)
    cv2.imwrite("dark.jpg", DC)
    A = estimate_A(img, DC)
    t = 1 - w*(DC/A)
    t[t < 0.1] = 0.1
    J = (img - A)/t[:,:,None] + A
    return J



if __name__ == "__main__":
    img = cv2.imread(sys.argv[1])
    out = reduction(img, float(sys.argv[2]))
    cv2.imwrite("out.png", out)