【画像処理】暗視カメラ出来たかも♬

前回は明度自動調整をやってみて、一応の成果はあったが、提供されているアプリに比べるとどうも得られる画像が負けていた。
ということで、今回は少し中身に踏み込んで改善したので記事にしておく。
まず、最終的に自動的に得られるようになった画像を貼っておく。

前回の画像

元画像

やったこと

・改善の目安をつける--暗い画像とは
・gamma調整画像とは
・いろいろな空間で明るい画像と暗い画像のヒストグラムを書いてみる
・equalizeHistを使う
・ヒストグラムを外挿・補間する

・改善の目安をつける--暗い画像とは

画像の暗さはそもそも照明と反射率と感度の関係で決まっている。
暗い画像は、そのどれかが通常より減少している。
そう考えると、これは適当な処理を施せば完全復元が可能なように思える。
もともと画像は、三原色で記述できる。しかし、暗さはその三原色だけでは表現できない。つまり、暗さの原因は色以外の輝度に原因の端がある。
とはいえ、その輝度が色表現空間の何なのかを見出せば、ほぼ完全復元の可能性が高まると言える。

・gamma調整画像とは

前回も試しているが、輝度調整と言えばgammaがよく知られているので、以下のコードでgamma調整画像を書いてみる。
※なお、上記リンクはガンマで画像を調整している説明で本記事とは関係ありません

gamma_img.py

import cv2
import numpy as np
import matplotlib.pyplot as plt

img_src = cv2.imread('s1.jpg', 1)
for gm in range(10,50,2):
    gamma=gm/10
    lookUpTable = np.zeros((256, 1), dtype = 'uint8')
    for i in range(256):
        lookUpTable[i][0] = 255 * pow(float(i) / 255, 1.0 / gamma)
    # ルックアップテーブルによるガンマ補正
    img_gamma = cv2.LUT(img_src, lookUpTable)
    plt.imshow(img_gamma)
    plt.title('gamma_'+str(gamma)+'.jpg')
    cv2.imwrite('gamma_'+str(gamma)+'.jpg', img_gamma)
    plt.show()

実際、40枚出してみると代表的なものは以下のものが得られます。
$pow(x,y)=x^y$なので、gamma=1のときはlookUpTableの各要素は１になり、同じ数値を返します。
ということで、使えそうな画像はgamma=1.4～2.6あたりがよさそうです。




しかし、やはりよい再生画像に比べるとなんとなくボケているようです。
この変換は何をしているのでしょうかということで、img = cv2.cvtColor(img, cv2.COLOR_BGR2YCR_CB)してヒストグラムを計算すると以下のようになります。
gamma=1.8

gamma=1.0(元画像)

ということで、ほぼ輝度に対応するY成分をこのヒストグラムのより大きな値（255側）にLookUpTableに基づいて変換しているのです。
以下のグラフがこの上式でgammaを変化させて計算されたlookUpTableです。
直線から徐々にカーブがきつくなって、小さな数字から大きな数字への写像になっているのが分かります。

【参考】
・Python OpenCVの基礎　LUTで画像の暗部を持ち上げてみる
この変換を適当な色空間に対して最適に実施すればより鮮明で綺麗な色合いの画像が得られるものと考えられます。
ということで、今回は、まずいろいろな色空間で明るい画像と暗い画像について、ヒストグラムを書いて比較してみました。

・いろいろな空間で明るい画像と暗い画像のヒストグラムを書いてみる

ここではOpenCVで使える色空間でやってみた。

#-*- coding:utf-8 -*-
import cv2
import numpy as np
from matplotlib import pyplot as plt

# 入力画像を読み込み
img = cv2.imread("s1.jpg")
b, g, r = cv2.split(img)

def calc_hist(b,g,r):
    hist_r = cv2.calcHist([r],[0],None,[256],[0,255])
    hist_g = cv2.calcHist([g],[0],None,[256],[0,255])
    hist_b = cv2.calcHist([b],[0],None,[256],[0,255])
    return hist_r, hist_g, hist_b

def plot_hist(hist_r,hist_g,hist_b,file='0'):
    # グラフの作成
    a=plt.figure()
    plt.xlim(0, 255)
    plt.plot(hist_r, "-r", label="Red"+str(255))
    plt.plot(hist_g, "-g", label="Green"+str(255))
    plt.plot(hist_b, "-b", label="Blue"+str(255))
    plt.xlabel("Pixel value", fontsize=20)
    plt.ylabel("Number of pixels", fontsize=20)
    plt.legend()
    plt.grid()
    plt.savefig('original_'+str(file)+'.jpg')
    plt.show()
    plt.close()

def processing(img,file):
    b, g, r = cv2.split(img)
    hist_r, hist_g, hist_b=calc_hist(b,g,r)
    plot_hist(hist_r, hist_g, hist_b,file)

# 入力画像を読み込み
file1="s1"
img = cv2.imread(file1+'.jpg')
processing(img,file1)
file2="s7"
img7 = cv2.imread(file2+'.jpg')
processing(img7,file2)

img = cv2.cvtColor(img, cv2.COLOR_BGR2YCR_CB)
processing(img,file1+'_YCrCb')
img7 = cv2.cvtColor(img7, cv2.COLOR_BGR2YCR_CB)
processing(img7,file2+'_YCrCb')

img = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
processing(img,file1+'_HSV')
img7 = cv2.cvtColor(img7, cv2.COLOR_BGR2HSV)
processing(img7,file2+'_HSV')

img = cv2.cvtColor(img, cv2.COLOR_BGR2YUV)
processing(img,file1+'_YUV')
img7 = cv2.cvtColor(img7, cv2.COLOR_BGR2YUV)
processing(img7,file2+'_YUV')

以下省略
この図を見ると、rgb空間ではそれぞれの値が全体的に255側にシフトしているのが分かります。また、YCｒCｂ空間ではYはgamma変換と同じく255側に広がり、さらにCrやCbもより分布が広がって、100～150の間が100側に幅が広がっているのが分かります。
ということで、この各要素毎の変換のLookUpTableをリーズナブルに求められれば、暗い画像から明るく鮮明な画像への変換の問題は解決すると言えます。
※残念ながら、今回の記事ではこの３次元のLookUpTableを求める問題の解はありませんm_(・・)_m
　DL利用含め、次回以降の課題とします

・equalizeHistを使う

これは、分布を広げるために平坦化関数を利用するというアイディアです。
※参考はありませんが、多く使われているようです
結果は以下のような画像が得られます。はっきり言って、上記二つは妥協してしまうレベルですが、やはり輝度が出すぎで色相がいまいちな気がします。
equalizeHist_treat.jpg

equalizeHist_rgb.jpg

equalizeHist_yuv.jpg

equalizeHist_hsv.jpg

createCLAHE.jpg

おまけに掲載しているプログラムでその心臓部は以下の部分です。
rgbでは以下のとおり、それぞれをcv2.equalizeHist(r)しているだけです。

def histogram_equalize(img,size):
    img = cv2.resize(img, dsize=size, interpolation=cv2.INTER_CUBIC)
    b, g, r = cv2.split(img)
    red = cv2.equalizeHist(r)
    green = cv2.equalizeHist(g)
    blue = cv2.equalizeHist(b)
    return cv2.merge((blue, green, red))

この場合のヒストグラムは以下のとおりです。
equalizeHist_rgb_rgb

equalizeHist_rgb_YCrCb

そして、もう一方の処理は以下のとおり、一度bgrに対して平坦化をした後に、createCLAHE処理をLAB空間のL成分に実施しています。なお、平坦化は（4，4）で実施しています。

def histogram_equalize_treat(img,size):
    img = cv2.resize(img, dsize=size, interpolation=cv2.INTER_CUBIC)
    b, g, r = cv2.split(img)
    red = cv2.equalizeHist(r)
    green = cv2.equalizeHist(g)
    blue = cv2.equalizeHist(b)
    bgr = cv2.merge((blue, green, red))
    lab = cv2.cvtColor(bgr, cv2.COLOR_BGR2LAB)
    lab_planes = cv2.split(lab)
    clahe = cv2.createCLAHE(clipLimit=2.0,tileGridSize=(4,4))
    lab_planes[0] = clahe.apply(lab_planes[0])
    lab = cv2.merge(lab_planes)
    bgr = cv2.cvtColor(lab, cv2.COLOR_LAB2BGR)
    return bgr

こうして見てくると、やはり平坦化では限界があることが分かります。
ということで、いよいよ外挿・補間をやってみようと思います。

・ヒストグラムを外挿・補間する

結局、rgb空間でそれぞれの色のヒストグラム関数を外側に適切に伸ばして間を補間すれば、解決するということになります。
これは、通常以下の関数で実施するようです。
※以下の解釈は以下を参考にしています

img1 = cv2.split(img) #rgb要素に分解
for s in range(0,3): #rgbの各要素に対して以下を実行
    hist,bins = np.histogram(img1[s].flatten(),256,[0,256]) #[0,256]区間のヒストグラムの要素を256階調でヒストグラム化
    cdf = hist.cumsum() #その要素までの累積値を配列で返す
    cdf_normalized = cdf * hist.max()/ cdf.max() #累積値とヒストグラムの最大値で規格化
    cdf_m = np.ma.masked_equal(cdf,0) #累積値0の領域をマスク
    cdf_m = (cdf_m - cdf_m.min())*255/(cdf_m.max()-cdf_m.min()) #累積値0の領域を255に広げる
    cdf = np.ma.filled(cdf_m,0).astype('uint8') #cdf_mに0を入れてcdfとする
    img1[s] = cdf[img1[s]] #cdfでimg1[s]を拡張する
img2 = cv2.merge(img1)

・numpy.cumsum
・numpy.ma.masked_equal
・numpy.ma.filled
・Image Enhancement
ということで、以下のおまけのhist_equalize.pyで上記のたぶん一番綺麗な画像が得られました。
ｒとｇ、ｂの要素化するヒストグラムを変更しています。
※ここが恣意的でさらにそれぞれの要素の分布を反映する数値を自動で計算すべきですがやっていません

for s in range(0,3):
    if s==0:
        hist,bins = np.histogram(img3[s].flatten(),256,[0,125])
    else:    
        hist,bins = np.histogram(img3[s].flatten(),256,[0,155])

    cdf = hist.cumsum()
    cdf_normalized = cdf * hist.max()/ cdf.max()
    lst = str(list[s])
    ax1[s].plot(cdf_normalized, color = 'b',label='cdf_'+str(lst))
    ax1[s].hist(img3[s].flatten(),256,[0,256], color = 'r',label='histogram_'+str(lst))
    ax1[s].set_xlim([0,256])
    ax1[s].legend(loc = 'upper left')

    cdf_m = np.ma.masked_equal(cdf,0)
    cdf_m = (cdf_m - cdf_m.min())*255/(cdf_m.max()-cdf_m.min())
    cdf = np.ma.filled(cdf_m,0).astype('uint8')
    img3[s] = cdf[img3[s]]
img4 = cv2.merge(img3)

この場合のヒストグラムは以下のようになります。
若干、ノイズがのっているものの下のimg5（参考サイトのアプリ）とほぼ似た画像が得られています。
※ノイズは上記の平滑化で消えそうですが、鮮度は落ちるのでやりません

また、この場合に作成された３つの変換テーブルcdfは以下のようになっています。

まとめ

・暗視カメラ向け。。暗い画像の輝度・鮮明化を実施した
・ヒストグラムの平坦化でもある程度の画像は得られる。
・外挿・補間をすることによりさらに鮮明な画像が得られる

・問題領域は明確になったので、さらなるシステムティックな画像復元（通常温度の太陽光下の）をやろうと思う
・画像復元は既に解決済みなようだが、コンボリューションの手法で触れや揺らぎからの復元ができるようなのでやってみようと思う
・DLを利用したさらに簡単で改善された手法をやって試してみようと思う

おまけ

hist_equalize.py

import cv2
import numpy as np
from matplotlib import pyplot as plt

def plot_hist(frame, frame1,k=0,l=0,frame_name='frame', frame1_name='frame1'):
    #ヒストグラムの表示
    ax[k,l].plot(frame,label=str(frame_name))
    ax[k,l].plot(frame1,label=str(frame1_name))
    ax[k,l].set_xlim([0, 256])
    ax[k,l].set_ylim([0,max(max(frame),max(frame1))])  #15000])
    ax[k,l].legend()
    #plt.show()

def plot_hist1(frame, k=0,frame_name='frame'):
    #ヒストグラムの表示
    ax[k].plot(frame,label=str(frame_name))
    ax[k].set_xlim([0, 256])
    ax[k].set_ylim([0,max(frame)])  #15000])
    ax[k].legend()
    #plt.show()

def something(frame,YCrCb=0):
    imgOrg = frame  #cv2.imread('61.jpg', 1)
    #BGRをYCrCbに変換します
    orgYCrCb = cv2.cvtColor(imgOrg, cv2.COLOR_BGR2YCR_CB)
    #輝度のヒストグラムを作成
    histOrgY = cv2.calcHist([orgYCrCb], [YCrCb], None, [256], [0, 256]) #0:Y 1:Cr 2:Cb
    return histOrgY

def something2(frame, frame1,YCrCb=0):
    imgOrg = frame  #cv2.imread('61.jpg', 1)
    imgLut = frame1  #cv2.imread('LUT.jpg', 1)
    #BGRをYCrCbに変換します
    orgYCrCb = cv2.cvtColor(imgOrg, cv2.COLOR_BGR2YCR_CB)
    lutYCrCb = cv2.cvtColor(imgLut, cv2.COLOR_BGR2YCR_CB)
    #輝度のヒストグラムを作成
    histOrgY = cv2.calcHist([orgYCrCb], [YCrCb], None, [256], [0, 256]) #0:Y 1:Cr 2:Cb
    histLutY = cv2.calcHist([lutYCrCb], [YCrCb], None, [256], [0, 256])
    return histOrgY,histLutY

#fig, ax = plt.subplots(1, 2, figsize=(12, 6))
#fig1, ax1 = plt.subplots(1, 3, figsize=(12, 4))
img = cv2.imread('s1.jpg',1)
img5 = cv2.imread('s7.jpg',1)
plt.imshow(img5)
plt.title("img5")
cv2.imwrite('img5.jpg', img5)
plt.show()
plt.imshow(img)
plt.title("img")
cv2.imwrite('img.jpg', img)
plt.show()

lab = cv2.cvtColor(img, cv2.COLOR_BGR2YCrCb)
fig1, ax1 = plt.subplots(1, 3, figsize=(12, 4))
img1 = cv2.split(lab)
list =('l','a','b')
for s in range(0,1):
    hist,bins = np.histogram(img1[s].flatten(),256,[0,100])
    cdf = hist.cumsum()
    cdf_normalized = cdf * hist.max()/ cdf.max()
    lst = str(list[s])
    ax1[s].plot(cdf_normalized, color = 'b',label='cdf_'+str(lst))
    ax1[s].hist(img1[s].flatten(),256,[0,256], color = 'r',label='histogram_'+str(lst))
    ax1[s].set_xlim([0,256])
    ax1[s].legend(loc = 'upper left')
    cdf_m = np.ma.masked_equal(cdf,0)
    cdf_m = (cdf_m - cdf_m.min())*255/(cdf_m.max()-cdf_m.min())
    cdf = np.ma.filled(cdf_m,0).astype('uint8')
    img1[s] = cdf[img1[s]]
plt.show()

img1[1]=img1[1]
img1[2]=img1[2]
img1 = cv2.merge(img1)
img2 = cv2.cvtColor(img1, cv2.COLOR_YCrCb2BGR)
plt.imshow(img2)
plt.title("img2")
cv2.imwrite('img2.jpg', img2)
plt.show()

fig, ax = plt.subplots(1, 3, figsize=(12, 4))
histOrgY=something(img2,0)
plot_hist1(histOrgY,0, "img2")
histOrgY=something(img2,1)
plot_hist1(histOrgY,1, "img2")
histOrgY=something(img2,2)
plot_hist1(histOrgY,2,"img2")
#cv2.imwrite('img2_hist.jpg', ax)
plt.show()

img3 = cv2.split(img)
list =('r','g','b')
hist_r = np.bincount(img3[0].ravel(),minlength=256)
hist_g = np.bincount(img3[1].ravel(),minlength=256)
hist_b = np.bincount(img3[2].ravel(),minlength=256)
# グラフの作成
plt.xlim(0, 255)
plt.plot(hist_r, "-r", label="Red")
plt.plot(hist_g, "-g", label="Green")
plt.plot(hist_b, "-b", label="Blue")
plt.xlabel("Pixel value", fontsize=20)
plt.ylabel("Number of pixels", fontsize=20)
plt.legend()
plt.grid()
plt.show()

for s in range(0,3):
    if s==0:
        hist,bins = np.histogram(img3[s].flatten(),256,[0,125])
    else:    
        hist,bins = np.histogram(img3[s].flatten(),256,[0,155])

    cdf = hist.cumsum()
    cdf_normalized = cdf * hist.max()/ cdf.max()
    lst = str(list[s])
    ax1[s].plot(cdf_normalized, color = 'b',label='cdf_'+str(lst))
    ax1[s].hist(img3[s].flatten(),256,[0,256], color = 'r',label='histogram_'+str(lst))
    ax1[s].set_xlim([0,256])
    ax1[s].legend(loc = 'upper left')

    cdf_m = np.ma.masked_equal(cdf,0)
    cdf_m = (cdf_m - cdf_m.min())*255/(cdf_m.max()-cdf_m.min())
    cdf = np.ma.filled(cdf_m,0).astype('uint8')
    img3[s] = cdf[img3[s]]

img4 = cv2.merge(img3)
plt.imshow(img4)
plt.title("img4")
cv2.imwrite('img4.jpg', img4)
plt.show()

fig, ax = plt.subplots(2, 3, figsize=(12, 4))
histOrgY,histLutY=something2(img5,img4,0)
plot_hist(histOrgY,histLutY,0,0,"img5", "img4")
histOrgY,histLutY=something2(img5,img4,1)
plot_hist(histOrgY,histLutY,0,1,"img5", "img4")
histOrgY,histLutY=something2(img5,img4,2)
plot_hist(histOrgY,histLutY,0,2,"img5","img4")
#cv2.imwrite('img4_hist.jpg',ax)
plt.show()

equalize.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import cv2
import numpy as np
import pylab as plt

def plot_hist(frame, frame1,k=0,l=0,frame_name='frame', frame1_name='frame1'):
    #ヒストグラムの表示
    ax[k,l].plot(frame,label=str(frame_name))
    ax[k,l].plot(frame1,label=str(frame1_name))
    ax[k,l].set_xlim([0, 256])
    ax[k,l].set_ylim([0,max(max(frame),max(frame1))])  #15000])
    ax[k,l].legend()
    #plt.show()

def something(frame, frame1,YCrCb=0):
    imgOrg = frame  #cv2.imread('61.jpg', 1)
    imgLut = frame1  #cv2.imread('LUT.jpg', 1)
    #BGRをYCrCbに変換します
    orgYCrCb = cv2.cvtColor(imgOrg, cv2.COLOR_BGR2YCR_CB)
    lutYCrCb = cv2.cvtColor(imgLut, cv2.COLOR_BGR2YCR_CB)
    #輝度のヒストグラムを作成
    histOrgY = cv2.calcHist([orgYCrCb], [YCrCb], None, [256], [0, 256]) #0:Y 1:Cr 2:Cb
    histLutY = cv2.calcHist([lutYCrCb], [YCrCb], None, [256], [0, 256])
    return histOrgY,histLutY

def treatise(frame,size):
    frame1 = cv2.resize(frame, dsize=size, interpolation=cv2.INTER_CUBIC)
    gridsize=8
    #moto_file=frame #'sora_7_after.png'
    bgr = frame1  #cv2.imread(s1.jpg,1) #k1.jpg s1.jpg
    lab = cv2.cvtColor(bgr, cv2.COLOR_BGR2LAB)
    lab_planes = cv2.split(lab)
    clahe = cv2.createCLAHE(clipLimit=3.0,tileGridSize=(gridsize,gridsize))
    lab_planes[0] = clahe.apply(lab_planes[0])
    lab = cv2.merge(lab_planes)
    bgr = cv2.cvtColor(lab, cv2.COLOR_LAB2BGR)
    return bgr

def histogram_equalize(img,size):
    img = cv2.resize(img, dsize=size, interpolation=cv2.INTER_CUBIC)
    b, g, r = cv2.split(img)
    red = cv2.equalizeHist(r)
    green = cv2.equalizeHist(g)
    blue = cv2.equalizeHist(b)
    return cv2.merge((blue, green, red))

def histogram_equalize_treat(img,size):
    img = cv2.resize(img, dsize=size, interpolation=cv2.INTER_CUBIC)
    b, g, r = cv2.split(img)
    red = cv2.equalizeHist(r)
    green = cv2.equalizeHist(g)
    blue = cv2.equalizeHist(b)
    bgr = cv2.merge((blue, green, red))
    lab = cv2.cvtColor(bgr, cv2.COLOR_BGR2LAB)
    lab_planes = cv2.split(lab)
    clahe = cv2.createCLAHE(clipLimit=2.0,tileGridSize=(4,4))
    lab_planes[0] = clahe.apply(lab_planes[0])
    lab = cv2.merge(lab_planes)
    bgr = cv2.cvtColor(lab, cv2.COLOR_LAB2BGR)
    return bgr

def histogram_equalize_hsv(img,size):
    img = cv2.resize(img, dsize=size, interpolation=cv2.INTER_CUBIC)
    img = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
    h, s, v = cv2.split(img)
    h = cv2.equalizeHist(h)
    s = cv2.equalizeHist(s)
    v = cv2.equalizeHist(v)
    bgr = cv2.merge((h, s, v))
    return bgr

def histogram_equalize_yuv(img,size):
    img = cv2.resize(img, dsize=size, interpolation=cv2.INTER_CUBIC)
    img = cv2.cvtColor(img, cv2.COLOR_BGR2YUV)
    # equalize the histogram of the Y channel
    img[:,:,0] = cv2.equalizeHist(img[:,:,0])
    # convert the YUV image back to RGB format
    img_output = cv2.cvtColor(img, cv2.COLOR_YUV2BGR)
    return img_output

video_input = cv2.VideoCapture(0)
size=(600,450)  #(2560,1920) #(1280,960) #(640,480)
#cv2.namedWindow("gammma correction", cv2.WINDOW_NORMAL)
while(1):
    #ret, frame2 = video_input.read()
    frame = cv2.imread('s7.jpg', 1)
    #frame1 = cv2.imread('s1_clahe_2_color.jpg', 1)
    frame2 = cv2.imread('s1.jpg', 1)
    cv2.imshow("frame", frame)
    cv2.imshow("frame2", frame2)
    bgr = treatise(frame2,size)
    cv2.imshow("equalizeHist", bgr)
    bgr = histogram_equalize(frame2,size)
    cv2.imshow("equalizeHist_rgb", bgr)
    bgr = histogram_equalize_yuv(frame2,size)
    cv2.imshow("equalizeHist_yuv", bgr)
    bgr = histogram_equalize_hsv(frame2,size)
    cv2.imshow("equalizeHist_hsv", bgr)
    bgr = histogram_equalize_treat(frame2,size)
    cv2.imshow("equalize_treat", bgr)

    fig, ax = plt.subplots(2, 3, figsize=(12, 4))
    histOrgY,histLutY=something(frame, frame2,0)
    plot_hist(histOrgY,histLutY,0,0,"frame", "frame2")
    histOrgY,histLutY=something(frame, frame2,1)
    plot_hist(histOrgY,histLutY,0,1,"frame", "frame2")
    histOrgY,histLutY=something(frame, frame2,2)
    plot_hist(histOrgY,histLutY,0,2,"frame", "frame2")
    #plt.pause(1)
    #plt.close()

    #fig, ax = plt.subplots(1, 3, figsize=(12, 4))
    histOrgY,histLutY=something(frame, bgr,0)
    plot_hist(histOrgY,histLutY,1,0,"frame", "bgr")
    histOrgY,histLutY=something(frame, bgr,1)
    plot_hist(histOrgY,histLutY,1,1,"frame", "bgr")
    histOrgY,histLutY=something(frame, bgr,2)
    plot_hist(histOrgY,histLutY,1,2,"frame", "bgr")
    plt.show()
    plt.close()

    k = cv2.waitKey(30) & 0xff
    if k == 27:
        break

暗視カメラ

はっきり言って、上記の技術使うと真っ暗な中でものが見える。。。ここまで来ると見えればどれでもいいのかもという気がしてくる。むしろ目的や使い方のようだ。
※YUVやHSVが案外いい結果を出している、ある程度明るければCLAHEもいい感じです