[MATLAB 영상처리 활용] Microglial cell extraction/segmentation review - part1

 

 

Microglial cell은 뇌신경계에서 일종의 immune cell로 알려져 있으며 최근에 다양한 시각으로

활발한 연구가 진행중입니다. 일반적인 실험/파이프라인 설계는 two-photon 또는 confocal microscope

등으로 microglia immune response의 dynamics 를 이미징 하고나서, 이미지 분석을 통해 의미있는

데이터를 추출하는 순서로 진행합니다.

 

실험 수행후 수많은 이미지 모음(image-set)을 만들고, 하나의 이미지에서 수많은 유의미한

데이터를 얻을 수 있습니다. 따라서 자동화된 cell 검출 (automatic cell segmentation) 알고리즘/프로그램이

많이 공개되어 왔고 현재에도 계속 발전되고 있습니다. 

 

 

본 포스트에서는 이미 공개된 microglial cell의 검출(segmentation) 알고리즘 묘사를 목표로 합니다.

cell segmentation의 알고리즘으로써 Kozlowsk et al., (출처1) 이 발표한 페이퍼의 methods를 차용합니다.

이 페이퍼가 제시하는 methods는 기초적인 수준의 단순한 이미지 프로세싱을 수행하므로 알고리즘이

직관적으로써 이해가 쉽고, 앞서 본 블로그에 게시한  [MATLAB 영상처리 기초] 내용들을 조합하여 묘사할 수

있습니다. 다만, 저수준의 methods로 인해 기대하는 결과를 얻지 못할수도 있습니다.

본 글에서는 biology 배경없이 철저히 image processing의 관점에서만 이야기 하겠습니다.

 

 

Kozlowsk et al., (출처1) 페이퍼의 Methods를 간략히 소개하자면,

Confocal microscopy based on green fluorescence intensities generated from imaging EGFP

expressing microglia of CX3CR1-EGFP transgenic mice 이미지 데이터를 대상으로 이미지

프로세싱을 수행합니다.

그 절차는,

 

1. 이미지상의 잠재적 microglial cell 후보들을(CP) 자동 검출합니다. 이를 위해 regional maxima intensity

     value를 이용합니다. 대상 이미지의 배경은 검은색(= 낮은 픽셀값)이고 microglial cell로 추정되는 물체는

     높은 intensity value(= 픽셀값)를 갖으므로 시도할 수 있는 방법입니다.

     이것은 MATLAB의 imregionalmax( ) 명령어로 구현할 수 있습니다. 그리고 microglial cell 후보들중

     크기가 50 µ㎡ (50 마이크로 제곱미터) 인 것들을 최종 후보로(LR) 선택 합니다.

   

    페이퍼에 이미지의 resolution이 언급되어 있지만, 이와같은 절대적인 수치(50 µ㎡)는 연구 목적,

    주관적 판단에 따라 그 기준이 다르므로 이 조건은 무시하겠습니다. 그리고 regional maxima intensity

    value의 특성을 이용하여 높은 픽셀값 분포를 갖는 cell을 검출한다는 아이디어는 간단하고 직관적으로

    보이지만 실제로 적용가능한지 검토해 보아야 합니다.

    microglial cell 후보들의 크기 즉, 픽셀합은 MATLAB의 bwarea( ) 등의 함수를 이용할 수 있습니다.

 

 

2. 개별 microglial cell 후보를(LR) 상대로 반복적인 Otsu thresholding을 수행하여 cell mask (CCM)를

    만듭니다. 이때 Otsu thresholding 연산을 위해 여러 level을 설정하여 이것들을 입력으로써 사용 합니다

    - MATLAB의 graythresh( ) 함수의 결과값이 level 입니다.

    ( [MATLAB 영상처리 기초 3] 이미지의 픽셀값 다루기 - Binarization 참조 )

    검출을 목표로 하는 cell의 크기를 TCS라 하며, 앞서 thresholding 연산의 결과로써 cell의 크기가

     TCS ± 100 µ㎡ (100 마이크로 제곱미터) 범위를 만족한다면 이 cell 후보는 '정확한(또는 진짜)

     microglial cell mask' (CM) 가 됩니다.

 

    실제 microglial cell은 3차원 운동을 하므로 2차원 이미지에서는 cell의 z축 움직임(위아래) 에 따라

     - cell이 존재함에도 - thresholding 결과가 달라집니다. 또한 검은색 배경이나 높은 intensity를 갖는

     cell안에 noise 등이 섞여 있음에도 하나의 level (값) 로써 cell mask 를 결정하는 thresholding 방식은

     cell segmentation에 적합하지 않습니다. 

    그럼에도 이미지 프로세싱의 기초 과정중 하나의 소규모 프로젝트로써 적합하다고 생각합니다.

    

 

3. 다음의 세가지 세부조건을 모두 만족한다면 적합한 microglial cell 로써 인정됩니다. 

    i) microglial cell은 정해진 크기의 cell mask 밖으로 넘어가면 안됩니다. 

    ii) cell soma의 크기(픽셀합)는 16.7 µ㎡ (16.7 마이크로 제곱미터) 보다 커야하며 cell mask의 50%

        이상을 차지한다.

    iii) 하나의 cell에는 하나의 soma만이(CSM) 존재 합니다. 

    

   여기에서도 마찬가지로 cell soma의 크기 16.7 µ㎡ 를 기준으로 삼는것을 본 글에서는 무시하겠습니다.

 

 

페이퍼에서는 위의 일련과정들을 하나의 figure 로써 제시함으로써 영상처리 프로세스에 대한 이해를 돕습니다.

 

 

The image processing strategy for segmenting cell shape of individual microglia (출처1)

 

 

Microglial cell extraction/segmentation - 1 microglial cell 후보군 찾기

 

이제는 실제의 microglial cell 이미지를 가지고 cell segmentation을 수행 하겠습니다.

시작하기에 앞서 주의할 점은, 위의 페이퍼 Kozlowsk et al., (출처1) 에서 사용하는 이미지를 구할 수

없어서 다른 페이퍼의 이미지로 대체합니다. 즉, 위의 페이퍼에서 제시하는 methods를 다른 페이퍼에서

공개하는 이미지에 적용하는 것입니다.

 

  

IBA1_EDF_#2_SN_L_3_Image 3_EDF_ch00 (출처2)

 

 

위 대상 이미지의 상세정보는 Anwer et al., (출처 2) 의 Materials and methods 부분을 읽어보시길

바랍니다. ( 위의 이미지는 본 블로그 맨밑 출처 2에서 제시하는 이미지로써 이곳에 가셔서 다운로드 

받을 수 있습니다. )

 

우선 기본적인 픽셀값 분포를 보기 위해 대상 이미지를 실수형 타입의 grayscale로 변환시킨후 히스토그램을

보겠습니다. 참고로, 배경이미지는 밝은색의 영향으로 픽셀값은 높을것이고 cell의 검정색을 띄고 있어

픽셀값이 낮을것입니다.

 

 

img = imread('IBA1_EDF_#2_SN_L_3_Image 3_EDF_ch00.tif');
img_gray = im2double(rgb2gray(img));

figure;
subplot(1,2,1); imhist(img_gray);
subplot(1,2,2); imshow(img_gray,[]);

 

 

 

 

 

위의 좌측 히스토그램에서 보듯이 픽셀값 분포는 주로 0.7 - 0.9 사이에 집중되어 있으며 이는

밝은 배경 때문임을 짐작할수 있습니다. 그리고 우측 그림 (그림을 클릭하여 확대) 에서 보듯이

배경이 되는 부분의 index는 0.8196, microglial cell은 0.1176 으로써 이미지가 반전되어 있습니다.

 

Cell segmentation을 수행함에 있어 Regional maxima intensity value 특성 즉, 배경의 낮은 픽셀값과

대비되는 microglial cell의 높은 픽셀값 분포를 위해, 이미지의 픽셀값 분포를 반전시켜야 합니다.

그리고 이는 간단하게도 1 - pixel value 연산을 함으로써 - 예를들어 배경의 index는 1 - 0.8196 = 0.1804,

cell은 1 - 0.1176 = 0.8824 - 이미지의 픽셀분포를 반전시킬 수 있습니다.

(픽셀분포의 반전 과정을 생략하고 regional minima intensity 를 적용해도 될것입니다 - MATLAB의

 imregionalmin() 함수)

 

 

Cell segmentation methods의 첫번째 과정인 imregionalmax( ) 함수를 적용한 microglial cell 후보(CP)

자동 검출을 위해 우선 임의의 cell을 하나 지정하여 작은 구역을 crop 하겠습니다

- [MATLAB 영상처리 기초 4] 이미지의 픽셀값 다루기 - Crop, 선형연산 참조.

그리고 임의지정한 cell에 대하여 imregionalmax( ) 명령어를 적용해 보겠습니다. 

 

 

임의 cell(세포) 지정

 

 

좌: 임의로 지정한 cell , 우: 임의지정한 cell에 imregionalmax( ) 를 적용한 결과

 

img_cell1 = img_gray(684:885, 992:1259);
img_cell1_regionMax = imregionalmax(img_cell1);
figure; 
subplot(1,2,1); imshow(img_cell1,[]);
subplot(1,2,2); imshow(img_cell1_regionMax,[]);

 

 

위 왼쪽 이미지는 배경과 cell의 픽셀분포를 반전시킨후 임의선택한 cell의 모습이고,

오른쪽은 imregionalmax( )를 적용한 결과모습 입니다. imregionalmax( ) 명령어는 grayscale의

이미지에 대해 local maxima를 찾아 black & white (= binary 이미지) 결과를 반환합니다.

예를들어 픽셀값의 분포가 아래와(출처 3) 같고 배경의 픽셀값은 10, 배경사이에 높은값 22, 33,

44, 45 가 보입니다.  

 

 

 

 

이 픽셀값 분포에 대해 local maxima를 구한다면, 22, 33, 45의 경우에만 local maxima로써 결과값 '1' 로써

치환됩니다. 여기에서 local maxima를 구하는 방법은 이미지 픽셀값 분포를 상대로 1차 도함수의 변곡점을

찾는것입니다. 따라서 44의 경우 이웃 픽셀의 값 45로 인해 local maxima로써 선택되지 못한것입니다.

Local maxima, minima 알고리즘에 대한 자세한 설명은 이 기술문서(출처 4)를 읽어 주시기 바랍니다.

 

 

 

 

그렇다면 위에서 임의선택한 cell의 픽셀값 분포는 어떨까요? Grayscale 이미지의 픽셀값 분포를

컬러로 mapping 하여 확인하여 보겠습니다. 그리고 cell의 soma를 수직(이미지의 130번째 열)과

수평(70번째 행)으로 가로지르는 임의의 픽셀범위를 설정하여 픽셀값의 분포 또한 plot 하여 보겠습니다. 

 

 

 

좌: 픽셀값의 변화(수직 방향), 우: 픽셀값의 변화(수평 방향)

 

 

Cell의 soma 부분은 높은 intensity를 갖고 주변으로 갈수록 낮아짐을 볼 수 있습니다.

그리고 soma를 중심으로 수직, 수평 범위의 픽셀값 분포를 그래프로 보자면 역시나 soma의

중심부가 가장 높은 intensity를 갖습니다. 그래프의 x축은 n번째 픽셀, y축은 픽셀값을 나타내며,

검정색 ' * ' marker는 local maxima를 나타냅니다 (그래프를 클릭하여 확대).

 

임의지정한 cell에 imregionalmax( ) 를 적용한 결과 이미지(밑의 좌측)의 하얀색 픽셀이 곧 soma의

intensity plot(밑의 우측)의 검정색 ' * ' marker 임을 나타냅니다 (그래프를 클릭하여 확대). 

 

 

좌: 임의지정한 cell에 imregionalmax( ) 를 적용한 결과, 우: 60번째 행의 픽셀값 변화(수평 방향)

 

 

앞에서 임의선택 하여 crop한 이미지 뿐만 아니라 이미지 전체에 대한 픽셀값 분포가 위 그래프와

같을것 임을 알 수 있으며, 따라서 microglial cell 후보군 찾기를 위해 regional maxima intensity

value의 특성을 이용한다는 아이디어는 실용성이 없을것으로 생각됩니다.

다만 추정하기로는 페이퍼에서 제시하는 입력 이미지가 이미 binarization된 이미지에 regional maxima

intensity value 를 연산하여 cell 후보군 찾기를 시도한 것이 아닌가 생각됩니다.

 

 

Cell 컬러 mapping 이미지의 수직, 수평방향 하얀색 선들은 manual로 그린것이며, colorbar( )과 plot( )

명령어의 세부설정은 MATLAB의 help 에서 검색하여 주의깊게 읽어보시길 권해드립니다.

 

 

figure; imagesc(img_cell1); 
c = colorbar(); 
c.Limits = [0, 1.0];
c.FontSize = 15;
set(gca,'xtick',[])
set(gca,'ytick',[])

ver_line = img_cell1(:,130);
hor_line = img_cell1(70,:);
ver_bw = img_cell1_regionMax(:,130);
ver_marker_pos = find(ver_bw);
hor_bw = img_cell1_regionMax(70,:);
hor_marker_pos = find(hor_bw);

figure; 
subplot(1,2,1); plot(ver_line,'r','LineWidth',4,'Marker','*',...
    'MarkerIndices',ver_marker_pos,'MarkerEdgeColor','k','MarkerFaceColor','k',...
    'MarkerSize',12); ylim([0 1.0]);
legend('vertical line','FontSize',18);
set(gca,'FontSize',20);
subplot(1,2,2); plot(hor_line,'b','LineWidth',4,'Marker','*',...
    'MarkerIndices',hor_marker_pos,'MarkerEdgeColor','k','MarkerFaceColor','k',...
    'MarkerSize',12); ylim([0 1.0]);
legend('horizontal line','FontSize',18);
set(gca,'FontSize',20);

 

 

 

 

 

출처1: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0031814

An Automated Method to Quantify Microglia Morphology and Application to Monitor Activation State Longitudinally In Vivo

 

 

 

출처2: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0284480#

A comparison of machine learning approaches for the quantification of microglial cells in the brain of mice, rats and non-human

본 글에서 사용된 microglial cell 이미지의 저작권 이용 설명

 

 

 

출처3: https://www.mathworks.com/help/images/ref/imregionalmax.html

Regional maxima - MATLAB imregionalmax

 

 

 

출처4: https://www.mathworks.com/help/symbolic/maxima-minima-and-inflection-points.html

Maxima, Minima, and Inflection Points - MATLAB & Simulink Example