본문 바로가기

Developement/C/C++

의료영상용 CLAHE 적용에 암부 보정하기.

 방사선을 이용하는 의료영상에 기본적인 후보정 같은걸 안하면 대부분 알아보기 힘든 수준의 결과만 얻게 됩니다. 이런걸 인간이 시각적으로 차이를 구별하고, 더 나아가 의료진이 임상적인 판단을 돕기 위해서 LC(Local Contrast) 를 증가 시키는 방법을 사용하는데, 아마 이중 가장 오래되고 효과적이면서 많이 쓰이는 것이 Contrast Limited Adaptive Histogram Equalization (CLAHE) 알고리즘이 아닐까 합니다.


 그런데 이 CLAHE 알고리즘을 그대로 사용하면 다음과 같은 문제가 생길 수 있습니다.


CLAHE 적용된 16비트 영상을 8비트로 조정한 영상


이 영상이 뭐가 문제인지는 바로 아래의 원본 영상을 대조 해 보면 차이를 알 수 있습니다. (빨간색 점선 안이 검게 타들어 가듯이 레벨이 전체적으로 변해 있는 것이 문제 입니다)


원본 16bit 영상을 8bit 로 조절한 영상


 원본대비, CLAHE 를 적용 하게 되면, 각 분할되는 영역의 복잡도와 그 부분의 레벨 평균에 따른 차이 때문에 전체 영상으로 합쳐 졌을때 부분적으로 음역이 생기고, 이 형태는 나눠진 블럭의 크기와 모양에 따라 차이가 발생하게 됩니다. 그래서 변화가 크게 생기는 오브젝트의 외곽부를 따라 음영처럼 그늘이 지게 되는 문제가 생깁니다.


 이런 문제를 해결 하기 위해서 고려 해 볼 것은 바로 Shading correction 알고리즘이 되겠습니다. 이 알고리즘을 만들게 된 원래 논문은 'Correction of uneven illumination (vignetting) in digital microscopy images' 에서 착안 하였습니다. 작은 구경의 카메라로 인체 내부를 촬영 하게 되면 빛의 균일성(MTF)이 렌즈 외곽으로 갈수록 떨어지는데, 이를 보통 비네팅이라 불립니다. 이를 보정하기 위해서 빛의 양(illumination) 을 구하고, 이를 Gaussian blur 를 적용한 보정 map 을 구하여 원본 영상에 보상하는 원리 입니다.

 이때 원래 대로라면 Gaussian blur 를 통한 shade map 을 얻고, 이를 통해 원본 영상의 음영진 부분을 보상하는 방법을 사용해도 되겠습니다만, 사실 이렇게 하면 엄청나게 느립니다. Gaussian blur 의 특성상 재귀적호출 이 적용되는 크기에 따라 반복적으로 발생하고, 이는 과도한 CPU time 을 소모 하게 됩니다. 그래서 저의 경우는 Gaussian Blur 를 사용하지 않고, 이전에 구현해서 쓰고 있는 resize engine 을 사용하여 매우 빠르면서 효과적인 결과를 도출 하였는데, 이 방법은 다음과 같습니다.


(참고) https://github.com/rageworx/librawprocessor


  • 원본 영상을 간략화 시킬 Gaussian blur radial 만큼의 양으로 나눈 크기로 down scale 을 수행 하되, 각 손실되는 레벨을 평균적으로 유지할 수 있도록 BiLinear filter 를 적용함.
  • 다시 up scale 을 수행하되, 이때 Gaussina blur 처럼 효과를 얻어 내기 위해 B-Spline filter 를 사용하여, 선형성이 강한 복호값들을 가지도록 함.


 또한 의료영상은 픽셀정보 자체가 컬러 이미지의 RGB 를 통해 얻어 내는 illumination 의 값과 같다고 가정 할 수 있기 때문에, (물론 빛이 아니라 X-Ray 의 양을 말하고, 이를 바로 illumination 으로 생각할 수 있습니다) 다음과 같이 shading correction map 을 만들수 있습니다.



 제가 사용한 의료 영상은 보통 우리가 아는 어두은 값이 0 이 아닌 가장 높은 값이 어두은 Inverse 된 형태를 가지고 있습니다. Float 으로 치면 1.0f 가 가장 어두은 색으로 표현되고, 이는 실제 X-Ray 가 Detector 에 어떠한 저항체를 받지 않고 바로 전달 된 것을 의미 하며, 빛으로 치면 가시광선이 바로 센서에 다다른 것으로 이해 될 수 있습니다. 하지만 X-Ray 는 빛과 달리 반대의 개념으로 이를 이해 해야 하고, 과거 CR(필름을 사용하는 형태)들이 네거티브 필름과 같이 반전되어 있기 때문에 그대로 DR(Digital Radiograph, 디지털 X-Ray 영상) 에서도 동일하게 표현 합니다. 다만, DCM (DICOM) 파일에 포함 될 때엔 실제 데이터를 변환 하여, 0 (또는 0.0f) 가 가장 어두은 색으로 지정 되기도 합니다.


 암부 보정을 위해서는 원래 영상에서 그래도 만들어진 보정값을 써도 되지만, 색상이라 볼수 있는 각 픽셀의 암부 수준을 쉽게 이해 하기 위해서 원본이 반전 되어 있는 것을 다시 원래 어두운 값이 0 인 것으로 변환 해야 합니다. 그래서 위 이미지 에서는 실제 0 에 가까운 값이 밝게 표시 되고, 어두은 값이 1.0f 에 가까운 높은 값이 됩니다.


 이를 test 를 위한 코드를 만들어 보정을 시도 합니다. (GCC 6.4 및, OpenMP 를 사용하는 코드 입니다) , 또한 일정한 데이터 scale 을 유지하기 위해 unsigned short 로 저장되는 레벨을 float 으로 변환 하여 vector를 만들어 사용 하였습니다. 이는 보통 float vect = (float)us_level / 65535.0f; 와 같은 수식으로 std::vecotr< float > data; 형에 채워서 사용하였습니다.



    float* datasrc = (float*)src->data();
    float* datamap = (float*)smap->data();

    #pragma omp parellel for
    for( unsigned cnt=0; cnt< src->size(); cnt++ )
    {
        float corr_f = 0.0f;
        float thres_f = 0.075f;
        float boost_f = 1.2f;

        if ( datamap[ cnt ] < thres_f )
        {
            float prev_f = datasrc[ cnt ];
            float prev_min_f = prev_f * 0.1f;

            corr_f = ( thres_f - datamap[ cnt ] ) * boost_f;

            datasrc[ cnt ] = prev_f - corr_f;

            if( datasrc[ cnt ] > prev_f )
            {
                datasrc[ cnt ] = prev_f;
            }
            else
            if( datasrc[ cnt ] < prev_min_f )
            {
                datasrc[ cnt ] = prev_min_f;
            }
        }
    }


 코드는 매우 간단하게 만들어서 테스트 하였고, 여기엔 자동으로 보정에 대해 보상을 하는 형태는 아직 구현 되지 않았습니다만, CLAHE 를 통해 대부분의 다른 조건의 영상들을 비슷한 형태의 조건으로 만들어 주는 결과를 고려 해 볼때 다른 어떤 영상을 사용해도 비슷한 결과를 얻을 수 있어 보입니다.



 위 방법을 통해서 얻어지는 결과는 위 이미지 처럼 window level 을 조절해서 봤을 때 더욱이 도드라 지는 것을 어느정도 보상되어 진다는 것을 알수 있습니다. 그리고, 이를 다시 background masking  을 통해 원래 영상형태로 만들어 주면 ...



 위와 같은 결과를 얻어 낼 수 있습니다. 다소 CLAHE 를 통해 손실되거나 외곡되는 부분이 있는 편 이지만, 이는 또다른 signal filter 를 통해 완화 할 수 있게 됩니다. 또는 여기에 다시 Dyamic range 보정을 하여 균등한 레벨로 만들어 주기만 해도 별다른 signal processing 없이 anatomy 까지 알아 볼 수 있는 기초 수준을 만들어 내는 듯 합니다.


Drago 알고리즘을 적용한 이미지.



 다만 이 방법으론 전체적인 flatten 한 영상을 만들어 내지는 못 합니다. 특히 L-Spine 부의 임상적 관측에 필요한 요소들이 균일하게 보여 지는 요소부족은 어느정도 있으며, 이는 결국 다른 signal process 가 필요하다는 반증이기도 합니다.


 이 알고리즘은 추후 librawprocessor 에 별도의 기능으로 구현 하여 포함 시키도록 하겠습니다. 또한 librawprocessorMIT License 를 따르는 Open Source 이므로 사용권 획득만 하시면 사용에 제한이 없으므로, 많은 관심을 가져 주시면 감사하겠습니다.