Python3 & OpenCV で画像処理を学ぶ[2] 〜 情報可視化のためのヒストグラム基礎実験
1. はじめに
本記事では、画像の持つ特徴を定量的に知るために、その基本的な手段として、基礎的な統計量（平均, 分散など）を踏まえつつ、一次元・二次元ヒストグラムの算出・描画方法について学びます。また、HLVやL*a*b*色空間でヒストグラムを描画することで、画像データが持つ色の情報をうまく可視化することを試みます。
前回↓
はあまりコードを書きませんでしたが、今回は numpy や matplotlib で色々書いていきます。
import numpy as np import scipy.stats as sstats import pandas as pd from matplotlib import pyplot as plt import matplotlib.cm as cm import seaborn as sns import cv2 %matplotlib inline
2. ヒストグラム基礎実験
2.1 明度ヒストグラムと統計量の基礎
ヒストグラムは、ある画像の画素値がどのように分布するかを示すものでした。0 から 255 までの値について(横軸)、その値を取る画素の数(縦軸)を集計してグラフにしたもの、と見ることができます。
(アルマちゃんのドット絵のヒストグラム。ドット絵なので、画素のとる値が飛び飛びなのがよくわかります)
以下、基本的な統計量について確認しておきましょう。
最小値[最大値]は、その名の通り、ある画像の画素値の中で最小[最大]のものです。データの分布の両端の値とも言えます。画像処理的には、最大値と最小値の差を元にコントラストが算出されたりするようです。
平均 mean : $\mu^{2}$, 分散 variance : は、以下の式で表されます。
$k$ 個目の画素の値を$P_k$, 画像全体の画素の数を $N$ として、
平均 : $\mu = \dfrac{1}{N} \sum\limits_{k=1}^{N}P_k$
分散 :
分散の正の平方根を標準偏差 standard deviation : と言います。分散も標準偏差も、平均値から離れた画素が多いほど大きくなる統計量であるため、分布の広さ・狭さを見る指標になります。
中央値 median は、画素値を小さい方から順に並べて、ちょうど真ん中に位置する画素値です。大きい順に並べても同じです。偶数個の場合は、真ん中にある二つの値の平均を取るようです。メディアンフィルタなどで使用されるようです。
最頻値 mode は、画像のヒストグラムの中で一番頻度が高いものです。
データの分布が対称な形（正規分布など）であれば、平均と分散だけでも指標として良いのですが、そうでない場合、平均は外れ値の影響を受けやすいため（例 : 平均年収はしばしば上方の少数の外れ値に引っ張られる）、関心のある目的に応じて、最頻値や中央値が分布の性質を見るために使われたりします。
分布の偏りの方向に応じて、 最頻値≦中央値≦平均値 または 平均値≦中央値≦最頻値 という関係が成り立つようです。
以下、ヒストグラムにこれらの統計量を表示して描画してみました。 赤い太線が平均、その両側の少し離れたところにある赤い点線が平均±標準偏差の範囲。緑の線が中央値、黄色の線が最頻値です（少々見づらいですが……）。
mean : 114.0581962 stddev : 65.6394397318 median : 96.0 mode : 70
mean : 174.468091709 stddev : 50.5730009287 median : 183.0 mode : 244
説明の都合で後出しになりましたが、上記のヒストグラムの表示に使った関数のコードが以下です。だいたい matplotlib です。公式チュートリアルでも、ヒストグラムの表示には matplotlib を使っているようです。OpenCV の図形描画で頑張る方向性もあるようですが、普通にこちらの方が良さそうです。
def grayhist(img, stats=False): """ 入力 : BGR画像, 統計量の表示の有無(Optional) 出力 : グレースケールのimshow, Histgram(+統計量) """ # スタイルの設定。seabornの設定は残るため、常に最初に書いておく sns.set() sns.set_style(style='ticks') # プロット全体と個々のグラフのインスタンス fig = plt.figure(figsize=[15,4]) ax1 = fig.add_subplot(1,2,1) sns.set_style(style='whitegrid') ax2 = fig.add_subplot(1,2,2) # グレースケール画像化→三重配列に戻す img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) img = cv2.cvtColor(img, cv2.COLOR_GRAY2RGB) ax1.imshow(img) img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY) # 一次元配列化 img = np.array(img).flatten() # 本来rangeは[0,255]となるはずだが、255に値が集まりすぎたので弾いた img = img[img!=255] # 軸の刻み目 ax2.set_xticks(np.linspace(0,255,9).astype(int)) # ヒストグラムを計算→色をつける N, bins, patches = ax2.hist(img, range=(0,255), bins=256) for (i,patch) in enumerate(patches): color = cm.gray(bins[i]/256) patch.set_facecolor(color) if stats==True: # 統計量を表示する mean = img.mean() std = np.std(img) median = np.median(img) mode = sstats.mode(img)[0][0] # 統計量のラインをひく ax2.axvline(mean, color='#d95763', linestyle='solid', linewidth=3) ax2.axvline(median, color='#6abe30', linestyle='solid', linewidth=2) ax2.axvline(mode, color='#ba9641', linestyle='solid', linewidth=2) ax2.axvline(mean + std, color='#d95763', linestyle='dashdot', linewidth=1) ax2.axvline(mean - std, color='#d95763', linestyle='dashdot', linewidth=1) # 統計量の説明の文字 ax2.text(mean,N.max()*1.075, "$\mu$",color='#d95763',horizontalalignment='center') ax2.text(median,N.max()*1.18,"median", color='#6abe30',rotation=45) ax2.text(mode,N.max()*1.15,"mode",color='#ba9641',rotation=45) ax2.text(mean+std,N.max()*1.075, "$\mu+\sigma$", color='#d95763',horizontalalignment='center') ax2.text(mean-std,N.max()*1.075, "$\mu-\sigma$", color='#d95763',horizontalalignment='center') fig.tight_layout() plt.show() print("mean : ", mean) print("stddev : ", std) print("median : ", median) print("mode : ", mode) else: fig.tight_layout() plt.show()
実際に実行するときは以下のようになります。
# 2nd arg: 1=BGR wiz = cv2.imread('Wizdomiot.png',1) pdr = cv2.imread('pandra.png',1) alma = cv2.imread('alma.png',1) grayhist(alma) grayhist(pdr, True) grayhist(wiz, True)
2.2 RGB, HSV, L*a*b* 空間でのヒストグラム
先ほどはグレースケールの明暗の値についてヒストグラムを見ました。そこで次に、色情報をヒストグラムで見ることを考えます。
まずは基本ということで、RGBチャンネルごとのヒストグラムを見てみます。
先ほどはヒストグラムの計算にmatplotlibを用いましたが、今回はOpenCVの cv2.calcHist([img], [channel], mask_img, [binsize], ranges) を使ってみます。
wiz = cv2.cvtColor(cv2.imread('Wizdomiot.png',1), cv2.COLOR_BGR2RGB) pdr = cv2.cvtColor(cv2.imread('pandra.png',1), cv2.COLOR_BGR2RGB) lyrith = cv2.cvtColor(cv2.imread('Lyrith.png',1), cv2.COLOR_BGR2RGB) sb = cv2.cvtColor(cv2.imread('SCARY BANQUET1.png',1), cv2.COLOR_BGR2RGB) alma = cv2.cvtColor(cv2.imread('alma.png',1), cv2.COLOR_BGR2RGB)
実際に書いた関数のコードが以下です。
def rgb_hist(rgb_img): sns.set() sns.set_style(style='ticks') fig = plt.figure(figsize=[15,4]) ax1 = fig.add_subplot(1,2,1) sns.set_style(style='whitegrid') ax2 = fig.add_subplot(1,2,2) ax1.imshow(rgb_img) color=['r','g','b'] for (i,col) in enumerate(color): # 各チャンネルのhist # cv2.calcHist([img], [channel], mask_img, [binsize], ranges) hist = cv2.calcHist([rgb_img], [i], None, [256], [0,256]) # グラフの形が偏りすぎるので √ をとってみる hist = np.sqrt(hist) ax2.plot(hist,color=col) ax2.set_xlim([0,256]) plt.show()
実行すると、以下のようになります。
rgb_hist(wiz)
rgb_hist(pdr)
CG制作ソフトウェアのユーザーにとっては、情報として比較的見慣れたものであると思います。
トーンカーブなどで色調補正をするときに、いま明部にどれだけ赤成分があるか、暗部の青成分はどれくらいか、などといったことを見ながら確認したりしますね。
次に、HSV空間とL*a*b*空間でヒストグラムを描くことを考えます。
HSVは色相 Hue, 彩度 Saturation, 明度 Value(Brightness)の三軸からなる色空間で、L*a*b*は L* が輝度, a*が赤-緑成分, b*が黄-青成分を示す補色空間でした。(例えばa*は負値を含む軸で、プラス方向が赤、マイナス方向が緑)
RGB空間に比べて、比較的直感的な軸であったり、人の色の認知に即した形で情報が表現されているため、より画像の色の性質をわかりやすく示してくれるかもしれません。
以下、とりあえずHSVとL*a*b*のそれぞれのチャンネルでヒストグラムを表示する関数のコードを書いてみました。
def axes4(): """ 四分割したプロットを描画するときのための関数。 axesが4つ入ったリストを返す """ sns.set() fig = plt.figure(figsize=[15,15]) sns.set_style(style='ticks') ax1 = fig.add_subplot(2,2,1) sns.set_style(style='whitegrid') ax2 = fig.add_subplot(2,2,2) ax3 = fig.add_subplot(2,2,3) ax4 = fig.add_subplot(2,2,4) axes = [ax1,ax2,ax3,ax4] return axes def tri_hists(rgb_img, color='rgb'): """ 'rgb','hsv','lab'の指定を受け取って、 「入力画像」と「各チャンネルのヒストグラム」の計4枚のプロットを描画する """ axes = axes4() axes[0].set_title('image') axes[0].imshow(rgb_img) if color=='rgb': img_ch = img color=['R','G','B'] elif color=='hsv': img_ch = cv2.cvtColor(rgb_img, cv2.COLOR_RGB2HSV) color=['Hue','Saturation','Brightness'] elif color=='lab': img_ch = cv2.cvtColor(rgb_img, cv2.COLOR_RGB2Lab) color=['L*','a*','b*'] else: return for (i,col) in enumerate(color): # 各チャンネルのhist if col=='Hue': # Hueの値域は[0,180) hist = cv2.calcHist([img_ch], [i], None, [180], [0,180]) else: hist = cv2.calcHist([img_ch], [i], None, [256], [0,256]) hist = np.sqrt(hist) axes[i+1].set_title(col) axes[i+1].plot(hist) plt.show()
実行結果は以下の通りです。
tri_hists(wiz,'hsv')
tri_hists(wiz,'lab')
……いまいちよくわかりませんね。
明るさを表すBrightnessやL*は良いとして、Hueも色相の色と対応させればわかりそうですが、「彩度が 100/255 付近に多く分布している」などと言われても解釈に困ります。
「よくわからない」原因は色々考えられますが、いくつか挙げられるのは、彩度やa*, b*単独では最終的な「色」の想像がつきにくいことであったり、a*,b*軸は単純に馴染みが薄いことなどがあると思います。
そこで、二軸ヒストグラムの描画を試みます。例えば、少なくとも彩度は色相とセットでないと色がわからないため、彩度軸×色相軸を組み合わせたヒストグラムを作り、グラフ上で対応する色とその頻度をセットで見てみることで、（明度を一旦無視した）「色」の頻度分布を明らかにできる可能性があります。
この場合、 (色相, 彩度) のペアに 頻度 を対応させる二変数関数となるので、情報としては三次元となります。それをそのまま立体的な三次元グラフにしても良いのですが、隠れてしまう部分が発生するなどして見にくい可能性があるので、今回は三軸目に明度を割り当てて表現することを考えます。
などと書いていてもイメージがつきにくいので、色相-彩度の二次元ヒストグラムを描画する関数を書いてみました。以下に示します。
def hs_hist(rgb_img): """ [目的] 「横軸に 色相 , 縦軸に 彩度 をとり、 点(H,S)における頻度を 明度 で表現するグラフ」を作りたい。 [実装] 入力のRGB画像をHSVに変換し、(H, S)の二次元ヒストグラムを計算して(H, S, 頻度)の配列を作る。 その配列をHSV2RGBとして変換し、RGB画像としてimshowに表示させる関数。 """ sns.set() sns.set_style(style='ticks') img_hsv= cv2.cvtColor(rgb_img, cv2.COLOR_RGB2HSV) # (H,S) だけにして、H,Sも分離 hs_2d = img_hsv[:,:,0:2].reshape(-1,2) h = hs_2d[:,0] s = hs_2d[:,1] # ヒストグラムのbinの設定。OpenCVにおいてHの値域は0~180である（0-255に納めるためか） hedges = np.arange(0,180) sedges = np.arange(0,255) # 二次元ヒストグラム H, xedges, yedges = np.histogram2d(h,s, bins=(hedges, sedges)) H = H.T # log scaleで偏りを緩和 & 正規化 H_log = np.log(H+1) H_norm = H_log/H_log.max()*255 # (H,S,頻度)の配列にするために、まずH[S]を縦[横]にリピートし、x行y列の配列にする x = H_norm.shape[1] y = H_norm.shape[0] hue_xy = np.repeat(xedges[:-1],y).reshape(x,y).T sat_xy = np.repeat(yedges[:-1],x).reshape(y,x) # depth方向にくっつけて、(H,S,頻度)の配列にする。uint8型でないとcvtColorが受けつけないらしい HS_hist = np.dstack((hue_xy, sat_xy, H_norm)).astype('uint8') HS_hist_im = cv2.cvtColor(HS_hist, cv2.COLOR_HSV2RGB) HS_hist_im = cv2.resize(HS_hist_im,(360,100)) # 以下、元の画像をax1に、ヒストグラムをax2に表示する fig = plt.figure(figsize=[15,4]) ax1 = fig.add_subplot(121) ax2 = fig.add_subplot(122) ax1.imshow(rgb_img) ax1.set_title('image') ax2.set_title('Hue-Saturation Histgram') ax2.set_xlabel('Hue') ax2.set_ylabel('Saturation') ax2.set_xticks(np.linspace(0,360,13)) ax2.set_xlim(0,360) ax2.set_ylim(0,100) ax2.imshow(HS_hist_im,origin='low',interpolation='bicubic') plt.show()
実行すると、以下のようになります。
hs_hist(wiz)
縦軸に彩度、横軸に色相をとり、頻度を明度で表現した図を、左側の元画像と共に右に表示しています。こちらのWizdomiotの例でいうと、魔法陣の青〜紫のグローと、画面下部の緑色の光をヒストグラム上で見て取ることができます。
hs_hist(sb)
このSCARY BANQUETの例では、(H=270~300,S=30~50)付近に月の色の分布が目立ちます。Hue=0付近に分布しているのは、床の低彩度の赤〜橙色や本棚ですね。
hs_hist(lyrith)
最後のLyrithの例は、元画像からは赤〜黄色系統の色が一見目立ちますが、ヒストグラムを見ると緑〜青色成分もそれなりに含まれていることがわかります。実際、画面下部の暗部などをよく見てみると、黄緑っぽい緑成分や紫がかった青成分の存在を認めることができると思います。（他人事のように解説していますが、むろん制作したのは僕なので、意図してそのようにしています）
先ほどと比べて、「どんな色がどのように分布しているか」がわかりやすい情報が得られたのではないかと思います。左に表示した元画像と見比べても、このヒストグラムにはデータとして納得感があるのではないでしょうか。
さて、次も同様にa*b*ヒストグラムを描いてみましょう。どうなるでしょうか。まずは書いたコードを示します。
def ab_hist(rgb_img): """ 基本的にhs_histと同じ """ sns.set() sns.set_style(style='ticks') img_lab = cv2.cvtColor(rgb_img, cv2.COLOR_RGB2Lab) # (a,b) だけにして、a,bも分離 ab_2d = img_lab[:,:,1:3].reshape(-1,2) a = ab_2d[:,0] b = ab_2d[:,1] aedges = np.arange(0,255) bedges = np.arange(0,255) H, xedges, yedges = np.histogram2d(a,b, bins=(aedges, bedges)) H = H.T H_log = np.log(H+1) H_norm = H_log/H_log.max()*255 # (頻度,a,b)の配列とするために、まずa[b]を縦[横]にリピートし、x行y列の配列にする。 x = H_norm.shape[1] y = H_norm.shape[0] a_xy = np.repeat(xedges[:-1],y).reshape(x,y).T b_xy = np.repeat(yedges[:-1],x).reshape(y,x) # depth方向にくっつけて、(頻度,a,b)の配列にする。 ab_hist = np.dstack((H_norm, a_xy, b_xy)).astype('uint8') ab_hist_im = cv2.cvtColor(ab_hist, cv2.COLOR_Lab2RGB) ab_hist_im = cv2.resize(ab_hist_im,(255,255)) fig = plt.figure(figsize=[15,4]) ax1 = fig.add_subplot(121) ax2 = fig.add_subplot(122) ax1.imshow(rgb_img) ax1.set_title('image') ax2.set_title('a*b* Histgram') ax2.set_xlabel('G ← a* → R') ax2.set_ylabel('B ← b* → Y') ax2.axvline(255/2, color='white', linestyle='solid', linewidth=1) ax2.axhline(255/2, color='white', linestyle='solid', linewidth=1) ax2.imshow(ab_hist_im,origin='low',interpolation='bicubic') plt.show()
関数を実行すると以下のようになります。
ab_hist(wiz)
L*a*b*空間は、a*は赤-緑成分、b*は黄-青成分で、正負で補色関係になっている空間でした。すなわち、右のヒストグラムでいうと、十字の真ん中である原点(無彩色)に対して対称に分布する成分があれば、補色の成分が多く含まれていると解釈することができそうです。（軸の目盛りは無視してください……）
上の例を見ると、基本的にa*軸のマイナス(=緑)、b*軸のマイナス(=青)成分に振れており、補色関係にあると言える成分は少ないことがわかります。わかりにくいのですが、画像上部の魔法陣の内側に赤に近い紫成分（a*軸-,b*軸+）が微妙に含まれていて、これが緑成分と補色をなしていると言えなくもないでしょうか……。
ab_hist(sb)
こちらは、基本的にa*軸+,b*軸-の成分、すなわち紫が多いことがわかります。こちらの画像で補色を増やそうと考える場合、例えば建物内部の影になっているところの暗部に黄緑成分をトーンカーブで付加する、などといった案をヒストグラムから挙げることができそうです。
ab_hist(lyrith)
こちらは、基本的にa*軸+,b*軸+の成分、すなわち橙色の成分がかなり多いことがわかります。彩度も幅広く分布しており、先の例に比べて色数はかなり多い方と言えるのではないでしょうか。また、a*軸+,b*軸+の橙色に対応する形で、a*軸-,b*軸-の青〜緑成分が含まれていることが確認できます（こちらは彩度など抑えめではありますが）
3. おわりに
3.1 まとめ
今回はヒストグラムの描き方を学び、また、二次元ヒストグラムの実験を通じて、画像の特徴を可視化する方法について考えることができました。
もちろん、「画素値の頻度を集計する」ヒストグラムは、データの可視化の手段としては比較的シンプルなものであり、「画像の特徴」を考える上で抜け落ちている重要な情報も多いと考えられます。例えば、「明度の高いピクセルが多いこと」がわかっても、それが画像中のどこに多く分布しているか？などの位置に関する情報は、ヒストグラムを通じて知ることはできません。
しかしむしろ、上記で見てきたように、それ自体は単なる集計の可視化であるヒストグラムからでも、様々な観点からデータを見ることでわかることは想像以上に多かったというのが僕の感想です。そして今回取り上げた視点以外にも、ヒストグラムを通じて見つけられる知見はまだまだ多くありうるのではないかと思います。
3.2 今後の課題
今回至らなかった以下の点については、今後の課題としたいと思います。
- 描画されたヒストグラムについて、値の前処理が一貫しなかったこと。
色の分布の形状を見やすくするためにと、出力を見ながら平方根をとったりlogをとったりしましたが、正当性のある補正方法が何なのかは結局わかりませんでした。
あまり綺麗な描画ではなかったと思います。原因として考えられるのは、OpenCVがHの値を[0-180)としていること(色相の分解能が2°どまりになる)や、そもそも変換時に整数しか認められていなかったことによる情報落ちなどでしょうか。また、せっかくなので極座標に変換するなどして円形で見てみればよかったですね。
3.3 次回の方針
次回↓は、ヒストグラムなどを用いつつ、実際に画像を操作・変換する処理の基本を学ぶ予定です。
