16bit版ヒストグラム平坦化処理

### はじめに

Opencvには画像処理のいろいろな関数がありますが、多くの関数は8bit(0-255)をベースに作られています。研究用途の画像データでは、より解像度の高い16bitが用いられています。
opencvでは実装されていない16bit Grayscaleのヒストグラム平坦化処理について覚えとして記載します。

[OpenCV-Python Tutorials 1 documentation » OpenCV-Pythonチュートリアル » OpenCVを使った画像処理 » ヒストグラム その2: ヒストグラム平坦化](http://labs.eecs.tottori-u.ac.jp/sd/Member/oyamada/OpenCV/html/py_tutorials/py_imgproc/py_histograms/py_histogram_equalization/py_histogram_equalization.html#histogram-equalization)
を参考に16bit版を作成しました。

### 環境

windows 10 Pro
Anaconda
Python 3.9
opencv 4.5

### コードと解説

画像データをopencvで読み込むとnumpyのarray形式になります。
通常の画像ファイルのデータタイプはuint8（符号なし8ビット整数型）になります。
Leenaさんの画像を題材します。
![Lenna_gray.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/165987/5a145077-2d56-83b6-0c6f-cc93d6581f76.png)

```python
# モジュールをインポートして、題材としてLennaさんを読み込みます。
import cv2
import numpy as np
from matplotlib import pyplot as plt

# from PIL import Image

# 画像をグレースケールで読み込む
file_name = 'Lenna_gray.png'
img = cv2.imread(file_name, cv2.IMREAD_ANYDEPTH)
# img = cv2.imread(file_name, cv2.IMREAD_GRAYSCALE | cv2.IMREAD_ANYCOLOR)

print(img.dtype)
# I = Image.open(file_name)
# img = np.array(I)

# >>> uint8
```



```python
# 16bitのデータに変換します。
img_16 = img.copy()
img_16 = img_16.astype('uint16')

# img_16 = img_16*256
print(img_16.dtype)
print(img_16.max())
# >>> uint16
# >>> 245
```


```python
### 8bit版（チュートリアルと同じ）チュートリアルを再現してみます。
def equalize_hist8(img, info=True):
    img2 = img.copy()
    hist,bins = np.histogram(img.flatten(),256,[0,256])

    cdf = hist.cumsum()
    cdf_normalized = cdf * hist.max()/ cdf.max()
    
    if info: 
        plt.imshow(img, cmap = "gray")
        plt.title('Original')
        plt.show()
        
        plt.plot(cdf_normalized, color = 'b')
        plt.hist(img.flatten(),256,[0,256], color = 'r')
        plt.xlim([0,256])
        plt.legend(('cdf','histogram'), loc = 'upper left')
        plt.show()
    

    cdf_m = np.ma.masked_equal(cdf,0)
    cdf_m = (cdf_m - cdf_m.min())*255/(cdf_m.max()-cdf_m.min())
    cdf_m2 = np.ma.filled(cdf_m,0).astype('uint8')
    img2 = cdf_m2[img]
    hist2,bins2 = np.histogram(img2.flatten(),256,[0,256])

    cdf2 = hist2.cumsum()
    cdf2_normalized = cdf2 * hist2.max()/ cdf2.max()

    
    if info: 
        plt.imshow(img2,cmap = "gray")
        plt.title('EqualizeHist')
        plt.show()
            
        plt.plot(cdf2_normalized, color = 'b')
        plt.hist(img2.flatten(),256,[0,256], color = 'r')
        plt.xlim([0,256])
        plt.legend(('cdf','histogram'), loc = 'upper left')
        plt.show()
    
    return img2

_ = equalize_hist8(img, info=True)

```
![Leena8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/165987/ae5fa900-6155-b455-bd8f-f0ef239f0338.png)

![hist08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/165987/bef9a4ad-9cf6-5042-a9a5-9ffa42d0c15a.png)

![Leena_eq8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/165987/09cacbf6-9f58-fe88-c3a0-4795a06ae9c4.png)

![hist_eq8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/165987/ec6ee9a7-bd0a-1355-0a2a-c58d39755264.png)


```python
### 16bit版
# - 配列のデータタイプをuint16に変更
# - 少し時間がかかります。（2分ぐらい） histgramの計算に時間がかかります。

def equalize_hist16(img, info=True):
    img2 = img.copy()
    hist,bins = np.histogram(img.flatten(),65536,[0,65536])

    cdf = hist.cumsum()
    cdf_normalized = cdf * hist.max()/ cdf.max()
    
    if info: 
        plt.imshow(img,cmap = "gray")
        plt.title('Original')
        plt.show()
        
        plt.plot(cdf_normalized, color = 'b')
        plt.hist(img.flatten(),65536,[0,65536], color = 'r')
        plt.xlim([0,65536])
        plt.legend(('cdf','histogram'), loc = 'upper left')
        plt.show()
    
    cdf_m = np.ma.masked_equal(cdf,0)
    cdf_m = (cdf_m - cdf_m.min())*65535/(cdf_m.max()-cdf_m.min())
    cdf_m2 = np.ma.filled(cdf_m,0).astype('uint16')
    img2 = cdf_m2[img]
    
    hist2,bins2 = np.histogram(img2.flatten(),65536,[0,65536])
    cdf2 = hist2.cumsum()
    cdf2_normalized = cdf2 * hist2.max()/ cdf2.max()


    if info:   
        plt.imshow(img2,cmap = "gray")
        plt.title('EqualizeHist')
        plt.show()
         
        plt.plot(cdf2_normalized, color = 'b')
        plt.hist(img2.flatten(),65536,[0,65536], color = 'r')
        plt.xlim([0,65536])
        plt.legend(('cdf','histogram'), loc = 'upper left')
        plt.show()
    
    return img2 

_ = equalize_hist16(img_16, info=True)
```
![leena16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/165987/9903e016-719f-7ec3-9834-ebeceeb5f8d1.png)

![hist16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/165987/85aa33a8-edc6-3b65-ce81-89b5f07f7125.png)

![Leena_eq16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/165987/40b0a72a-fa4a-b5eb-647c-a479d5cb75be.png)

![hist_eq16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/165987/a30d4690-f024-1fe4-8b0f-a549a42e2537.png)


```python
### 自動でDataTypeを確認して計算（8、16の両方に対応できるように変更）

def equalize_hist(img, info=True):
    img2 = img.copy()
    dataType = img.dtype
    print(dataType)
    
    if dataType == 'uint8':
        tval = 256

    elif dataType == 'uint16':
        tval = 65536
    
        
    hist,bins = np.histogram(img.flatten(),tval,[0,tval])

    cdf = hist.cumsum()
    cdf_normalized = cdf * hist.max()/ cdf.max()
    
    if info:
        plt.imshow(img,cmap = "gray")
        plt.title('Original')
        plt.show() 
        
        plt.plot(cdf_normalized, color = 'b')
        plt.hist(img.flatten(),tval,[0,tval], color = 'r')
        plt.xlim([0,tval])
        plt.legend(('cdf','histogram'), loc = 'upper left')
        plt.show()
        
    cdf_m = np.ma.masked_equal(cdf,0)
    cdf_m = (cdf_m - cdf_m.min())*(tval-1)/(cdf_m.max()-cdf_m.min())
    cdf_m2 = np.ma.filled(cdf_m,0).astype(dataType)
    img2 = cdf_m2[img]
    
    
    if info:  
        plt.imshow(img2,cmap = "gray")
        plt.title('EqualizeHist')
        plt.show()   
        
        hist2,bins2 = np.histogram(img2.flatten(),tval,[0,tval])
        cdf2 = hist2.cumsum()
        cdf2_normalized = cdf2 * hist2.max()/ cdf2.max()
        
        plt.plot(cdf2_normalized, color = 'b')
        plt.hist(img2.flatten(),tval,[0,tval], color = 'r')
        plt.xlim([0,tval])
        plt.legend(('cdf','histogram'), loc = 'upper left')
        plt.show()
        
    return img2

_ = equalize_hist(img_16, info=True)
```



### まとめ

ヒスグラムの計算で時間がかかってしまっていますが、16bitに対応するヒストグラム平坦化処理のコードを作成しました。

opencv、SciPyのndimageなどは、基本8bitベースで作られているため、16bitのイメージを入れるとエラーが出たり、8bitにキャストされてしまうことがあるので注意が必要です。