MR Image preprocessing (1) : histogram equalization

여태까지 나는 단순히 medical image도 natural image의 연장선으로 생각하고 문제를 해결하려고 하였다. 따라서 pixel간의 intensity normalization을 제외하고는 ( memory error을 해결하기 위해 필요시 resizing도 시켜주었다.) 별다른 전처리 작업을 해주지 않았다. 그렇게 계속해서 어떻게 model의 아키텍쳐 부분을 변경할까를 고민하던 도중, 왜 medical image만의 특성을 활용하고 있지 않는지에 대한 의문이 들었다. natural image와는 다르게 MR image는 이미지를 얻어내는 vender별로, 혹은 receive coil의 배치에 따라 signal의 magnitude가 약간씩 달라질 수 있고, 거기에 따른 intensity값이 달라질 수 있다. 물론 이미지를 load한 이후 normalization을 해주긴 하지만, 그전에 (물론 어떤 task를 수행할것인지에 따라) 보다 근본적인 전처리 작업이 필요할 수 있으며 따라서 오늘 포스트는 이러한 전처리 작업 중 하나인 adaptive histogram equalization에 대해 소개하고 이를 어떻게 implementation 해주는지에 대해 이야기하고자 한다.

 

Histogram Equalization

https://docs.opencv.org/3.4/d4/d1b/tutorial_histogram_equalization.html

Histogram equalization은 image contrast를 향상시키기 위해 사용되는 이미지 프로세싱 테크닉이다. 먼저 histogram equalization 등장배경에 대해 설명하자면, MR 이미지에서 비슷한 intensity를 갖는 픽셀들이 모여있는 경우 자세한 내용을 보기 어렵다. 따라서 histogram equalization은 비슷한 intensity가 모여 있는 것을 intensity들이 골구로 퍼지게 해준다. 결국 별로 차이가 안나던 pixel (혹은 voxel)의 intensity 값을 바꿔 차이가 나게 해주는 것이다. (애초에 histogram이 intensity 별로 픽셀들의 개수를 bar로 나타내주는 것이기 때문에 이러한 bar들을 평준화한다는 것은 결국 픽셀들의 intensity 값을 조정해준다는 것) 

Histogram equalization를 수행하기 위한 핵심 공식은 아래와 같다. 여기서 cdf(v)는 해당 pixel intensity의 CDF를 의미하며, L은 총 pixel intensity의 길이로 보통 266이다. (일반적으로 0~255의 pixel intensity를 가지므로)

Adaptive Histogram Equalization

이 histogram equalization에는 하나의 마이너한 문제점이 있다. 기존의 HE는 전체 region에 대해서 intensity값을 재조정해주는 것이기 때문에, 이미지의 특정 구역이 다른 곳과 비교했을 때 다른 분포를 갖고 있는 경우, HE를 적용하게 되면 이미지가 왜곡된다. 따라서 전체 이미지가 아닌 여러개의 작은 region(patch)에 대해서 각각 histogram equalization을 해준다.

 

1. Pytorch implementation 방법

http://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/Python_html/70_Data_Augmentation.html

70_Data_Augmentation

def histogram_equalization(image, 
                           min_target_range = None, 
                           max_target_range = None,
                           use_target_range = True):
    '''
    Histogram equalization of scalar images whose single channel has an integer
    type. The goal is to map the original intensities so that resulting 
    histogram is more uniform (increasing the image's entropy).
    Args:
        image (SimpleITK.Image): A SimpleITK scalar image whose pixel type
                                 is an integer (sitkUInt8,sitkInt8...
                                 sitkUInt64, sitkInt64).
        min_target_range (scalar): Minimal value for the target range. If None
                                   then use the minimal value for the scalar pixel
                                   type (e.g. 0 for sitkUInt8).
        max_target_range (scalar): Maximal value for the target range. If None
                                   then use the maximal value for the scalar pixel
                                   type (e.g. 255 for sitkUInt8).
        use_target_range (bool): If true, the resulting image has values in the
                                 target range, otherwise the resulting values
                                 are in [0,1].
    Returns:
        SimpleITK.Image: A scalar image with the same pixel type as the input image
                         or a sitkFloat64 (depending on the use_target_range value).
    '''
    arr = sitk.GetArrayViewFromImage(image)
    
    i_info = np.iinfo(arr.dtype)
    if min_target_range is None:
        min_target_range = i_info.min
    else:
        min_target_range = np.max([i_info.min, min_target_range])
    if max_target_range is None:
        max_target_range = i_info.max
    else:
        max_target_range = np.min([i_info.max, max_target_range])

    min_val = arr.min()
    number_of_bins = arr.max() - min_val + 1
    # using ravel, not flatten, as it does not involve memory copy
    hist = np.bincount((arr-min_val).ravel(), minlength=number_of_bins)
    cdf = np.cumsum(hist)
    cdf = (cdf - cdf[0]) / (cdf[-1] - cdf[0])
    res = cdf[arr-min_val]
    if use_target_range:
        res = (min_target_range + res*(max_target_range-min_target_range)).astype(arr.dtype)
    return sitk.GetImageFromArray(res)

#cast the images to int16 because data[0] is float32 and the histogram equalization only works 
#on integer types.
output = histogram_equalization(sitk.Cast(inputimg,sitk.sitkInt32))

AHE를 나의 brain data에 대해서 적용해보았을 때의 이미지 결과는 아래와 같다. 오른쪽이 processing을 하기 전, 왼쪽이 processing을 한 이후이다. 전반적으로 object이 존재하지 않는 pixel에 대해서 노이즈가 많이 낀 것을 볼 수 있다. 찾아보니 AHE 단점 중 의미없는 지역의 pixel intensity값이 enhance되어 이미지 왜곡이 생길 수 있다고 한다.

이렇게 전처리한 데이터를 3dResNet34에 feed했을 때 MAE는 3.5->3.4로 0.1 줄었다.