OpenCVとNumPyでヒストグラムを作成する

# はじめに
私は撮影術に疎く、[「トライXで万全」「逆光は勝利」「世はなべて三分の一」「ピーカン不許可」「頭上の余白は敵だ」](https://shogakukan-comic.jp/book?isbn=9784098601905)といった偏った知識しか持っていない。
そんな私だから何年か前にちょいとお高いデジカメを買ったとき、液晶画面に表示されるグラフが何なのかをわかっていなかった。
だが、OpenCVの画像処理をいろいろ楽しんできた今ならばわかる。それは輝度のヒストグラムだ。

# Pythonでヒストグラムを作成する
ヒストグラムは別に画像処理に特化した技術ではない。統計の手法の一つとして広く使われており、QC七つ道具の一つとして名を馳せている。
PythonではOpenCVやnumpyがヒストグラムを求める関数を持っている。

統計の世界では棒グラフの数や幅をいくつとすべきかの議論がなされることがあるが、画像の輝度のヒストグラムならば「0から255まで（255を含む）の256階級、幅1で」とするのが普通だろう。

今回の画像は適度にRGBが偏っているレナ嬢とする。

|lenna.png|
|:-:|
|![lenna.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520216/bfbb717f-e3eb-7a09-476c-8ab8d00f1f9e.png)|

## OpenCV
OpenCVでは **cv2.calcHist()** を使う。必須な引数のみ書くと

```Python:
hist = cv2.calcHist(images, channels, mask, histSize, ranges)
```
となる。`images`、`channels`といった引数名からわかるように、これらはリストを指定する。タプルでも可。

### 引数

 - images　　入力画像のリスト。`[image]`と要素1のリストとして指定する。`np.uint8`のほか`np.float32`でも可。
 - channels　カラーチャンネル。グレースケール画像に対しては`[0]`を、BGRに対しては`[0][1][2]`を指定する。
 - mask　　マスク画像。特に指定しない場合は`None`とする。こんなのオプショナルにすればいいのに。
 - histSize　区間の数。ビンと呼ばれることもある。256階調を個別に算出するなら[256]とする。
 - ranges　　カウントする値の範囲。画像の全輝度に対しておこなうときは[0, 256]とする。[0, 255] ではない。

### 出力
 - hist　　ヒストグラム。画像数が1のとき、`[[0.000e+00], [0.000e+00], …]`といった2次元numpy配列を返す。各要素はfloat32。

この戻り値を見る限りではいくら`images`がリストだといっても複数の画像を指定して複数のヒストグラムを一挙に取得することはできないように思える。


### サンプル
RGB画像の各カラーチャンネルに対しヒストグラムを作る例。
OpenCVだから色の並びはBGRであることに注意。

```Python:
import cv2 
import matplotlib.pyplot as plt

filename = "lenna.png"
image = cv2.imread(filename)
COLORS = ["blue","green","red"]

for i, color in enumerate(COLORS):
    hist = cv2.calcHist([image], [i], None, [256], [0, 256])
    plt.plot(hist, color=color)

plt.show()
```

これにより、レナ画像は赤みが強いことが確認できた。
![lenna.hist.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520216/b38e9d3f-651e-7244-ddf4-5f7058fe8ce2.png)


## NumPy
numpyでは **np.histogram()** を使う。

```Python:
hist, bin_edges = np.histogram(a, bins=10, range=None, normed=None, weights=None, density=None)
```
という使い方をする。

### 引数（一部のみ）

 - a　　　入力行列。事前に一次元化する事例が多いが実際は**計算の際に1次元化される。**
 - bins　ビン数。省略可能で初期値は10。数値でもリストやタプルでも可（何なら文字列でも）。
 - range　ビンの範囲。省略可能で初期値はaの最小値から最大値まで。

### 出力
 - hist　　　　ヒストグラム。一次元のnumpy配列で各要素は整数。
 - bin_edges　ビンの各範囲。


### サンプル1　rangeを省略

`bins`を指定し`range`を指定しない状態で描画してみる。

```Python:
for i, color in enumerate(COLORS):
    hist, bin_edges = np.histogram(image[:, :, i], bins=256)
    plt.plot(hist, color=color)
```

|np.histogram() `bins`のみ指定|cv2.calcHist()|
|:-:|:-:|
|![np_hist_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520216/e3234d12-0380-c5cc-0436-5ef9575451e4.png)|![lenna.hist.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520216/b38e9d3f-651e-7244-ddf4-5f7058fe8ce2.png)|

赤に注目してほしい。レナ画像は全体的に赤みが強く輝度50以下の点は存在しないはずだが、今回のグラフでは正規のグラフの50～255の範囲が0～255に拡大してグラフ化されている。これが「rangeを省略すると最小値から最大値までになる」の結果だ。
グラフが千切り状態になっているのは、「最小値から最大値まで」を「0から255まで256階級」に拡大した際に、たとえば「100.1から100.9まで」といった中途半端な範囲が発生しているためだと考えられる。輝度は整数値だから「100.1から100.9まで」の範囲内にある点の数は0個、というわけ。

### サンプル2　rangeを指定

次に`bins`と`range`を指定してみる。

```Python:
for i, color in enumerate(COLORS):
    hist, bin_edges = np.histogram(image[:, :, i], bins=256, range=(0, 256))
    plt.plot(hist, color=color)
```

|np.histogram() `bins`と`range`を指定|cv2.calcHist()|
|:-:|:-:|
|![np_hist_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520216/34439dd0-a722-d56e-0d95-1f3faacb4e7c.png)|![lenna.hist.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/520216/b38e9d3f-651e-7244-ddf4-5f7058fe8ce2.png)|

まだまだ道は長いと思われていたが、さっそくほぼ等しいグラフが描けた。
これは本当に正しいのだろうか。

## 正しさを確認する
レナ画像のサイズは512×512で262144個のピクセルを持つ。ヒストグラムを取れば、ビン数にかかわらず度数の合計は262144になるはずだ。
また、この画像には赤の輝度=255のピクセルは

```python:
np.sum(image[:, :, 2]==255)    # 112
```

から112個あることがわかる。

### 例1　range=(0, 254)
明らかに正しくない例として`range=(0, 254)`と指定してみよう。ビン数は関係ないのでデフォ値のままとする。

```python:
hist, bin_edges = np.histogram(image[:, :, 2], range=(0, 254))
len(hist)        # 10　ビン数（デフォ値）
sum(hist)        # 262032　度数の合計
```

ここで出てきた262032は全ピクセル数262144から輝度255のピクセルの数112を引いた値。つまり輝度255のピクセルがカウントされていないことを意味している。

### 例2　bins=256, range=(0, 255)
次に`bins=256, range=(0, 255)`と指定してみよう。

```python:
hist, bin_edges = np.histogram(image[:, :, 2], bins=256, range=(0, 255))
len(hist)        # 256　ビン数
sum(hist)        # 262144　度数の合計
len(bin_edges)   # 257 bin_edgesの数
bin_edges
[  0.           0.99609375   1.9921875    2.98828125   3.984375
   4.98046875   5.9765625    6.97265625   7.96875      8.96484375
    中略
 249.0234375  250.01953125 251.015625   252.01171875 253.0078125
 254.00390625 255.        ]
```
度数の合計はピクセル数と等しく、すべてのピクセルがカウントされていることがわかる。

ところで、bin_edgesは`(length(hist)+1)`個の要素を持つ。
最初に0があって最後に255がある。最後とは256番目ではなく257番目のことだ。ビン数256で区切っているのに輝度255はそれから外れた257番目のクラスターに属している。これはよろしくない。
また、各階級が正確な1間隔ではなく割り切れない数字になっているのもよろしくない。

### 例3　bins=256, range=(0, 256)
次に`bins=256, range=(0, 256)`と指定してみよう。

```python:
hist, bin_edges = np.histogram(image[:, :, 2], bins=256, range=(0, 256))
len(hist)        # 256　ビン数
sum(hist)        # 262144　度数の合計
len(bin_edges)   # 257 bin_edgesの数
bin_edges
[  0.   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.
    中略
 238. 239. 240. 241. 242. 243. 244. 245. 246. 247. 248. 249. 250. 251.
 252. 253. 254. 255. 256.]
```

これでようやく255という値が256番目のクラスターに属するようになった。また、各階級の幅が1であることも確認できた。

このあたりは[こちらのサイト](https://campkougaku.com/2019/11/23/numpy-histogram/)に目を通しておいたほうが良いだろう。

https://campkougaku.com/2019/11/23/numpy-histogram/



## OpenCVのヒストグラムとNumpyのヒストグラムの比較
cv2.calcHist()にはnp.histogram()にあった`bin_edges`の出力がない。二つのヒストグラムが完全に一致するかどうか確認しよう。
出力のdtypeとshapeをあらためて表記するとこうだ。

||cv2.calcHist()|np.histogram()|→|統一フォーマット|
|:-|:-|:-|:-:|:-|
|**dtype**|dtype('float32')|dtype('int64')|→|dtype('float32')にする|
|**shape**|(256, 1)|(256,)|→|(256,)にする|

このあたりを統一して、内容について完全一致かどうか調べよう。ついでに処理時間も調べよう。
処理時間測定の部分はもっと洗練された書き方があるようだが本題がら外れるのでここではベタ書きしている。

```Python:
import cv2
import numpy as np
import matplotlib.pyplot as plt
import time

COLORS = ["blue","green","red"]
filename = "lenna.png"
image = cv2.imread(filename)

for i, color in enumerate(COLORS):
    print("-" * 20)
    print(color)

    # cv2.calcHist() リストでなくタプルで指示してみる　良い子は真似するな
    start_time = time.time()
    hist_cv = cv2.calcHist((image,), (i,), None, histSize=(256,), ranges=(0, 256))
    end_time = time.time()
    print("cv2.calcHist()", end_time - start_time)

    # np.histogram() rangeをリストで指示してみる
    start_time = time.time()
    hist_np, bin_edges = np.histogram(image[:, :, i], bins=256, range=[0, 256])
    end_time = time.time()
    print("np.histogram()", end_time - start_time)

    hist_cv_sameformat = hist_cv.flatten()                      # (256, 1) -> (256,)
    hist_np_sameformat = hist_np.astype(np.float32)             # int64 -> float32
    if np.array_equal(hist_cv_sameformat, hist_np_sameformat):  # 本当はarray_equalはdtypeを揃える必要はない
        print("完全一致")    
    else:
        print("不一致あり")   
```

結果はこう。

```
--------------------
blue
cv2.calcHist() 0.0020308494567871094
np.histogram() 0.005982637405395508
完全一致
--------------------
green
cv2.calcHist() 0.000997781753540039
np.histogram() 0.003998756408691406
完全一致
--------------------
red
cv2.calcHist() 0.0010008811950683594
np.histogram() 0.003998994827270508
完全一致
```

OpenCVとnumpyで同じ結果を得ることができた。また、OpenCVのほうがだいぶ速いことがわかった。
引数も、名称はともかく、使われ方はほぼ同じであることもわかった。片方が(0, 255)でもう片方が[0, 256]とかだったら目も当てられないところだ。


# 終わりに

np.histogram() の`bin_edges`の解説は省略するつもりだったが、これに触れたことによりrangeが不適切だったときの挙動が理解できるようになった。こういうのは自分で手を動かさないと身につかないものだ。