연금술사

요약

Motivation

1. 염색정규화(Stain normalization)은 타깃(target image, 또는 reference image)라고 불리는 이미지를 정하고, 이 타깃이미지에 맞춰서 소스이미지의 염색톤을 변화하는 것을 의미합니다.

2. 그러나, 염색정규화에서 이런 타깃이미지를 선정하는 과정이 매우 작위적이고(=관찰자마다 다르며), 타깃이미지를 어떤 것을 고르냐에 따라 인공지능 모델의 성능이 크게 차이가 납니다(아래 Fig1은 어떤 이미지를 고르냐에따라서 성능이 크게 좌우된다는 예시를 표하기위한 그림입니다).

Method: 염색정규화과정에 필요한 trainable parameters을 인공지능 모델에 포함하여 학습

저자들은 데이터로부터 염색정규화에 필요한 파라미터를 모델에 추가합니다. 이 레이어명은 "LStrainNorm"이라고 합니다.

방법은 크게 "입력이미지 전달" -> "각 색공간에서 정규화" -> "다시 RGB로 변경" -> "집계" 입니다.

1. 입력이미지는 RGB 색상에 있는 이미지입니다. $\mathbf{x} \in \mathbb{R}^{H \times H \times 3}$

2. 색 공간에서 정규화

   2.1. 색공간에서 정규화하기위해서, 임의의 색공간 $s$로 변환합니다. 논문내에서는 $\mathbf{x}'^{s} = Transform(\mathbf{x};s)$로 표기합니다.

   2.2. 채널별로 정규화를 진행합니다. $\hat{x}'_{ijk} = \gamma _{k} \frac{x'_{ijk} - \mu_{k}}{\sigma_{k} + \epsilon} + \beta_{k}$

     - 여기서 $\mu_{k}$은  채널별로 이미지의 평균 픽셀값을 의미합니다. 또한, $\sigma_{k}$도 픽셀의 표준편차값을 의미합니다.

     - $\gamma_{k}, \beta_{k}$은 trainable parameter로 스케일(scale)과 값증감(shift)에 활용됩니다. Batch normalization에도 활용되는 값입니다.

     - 이전의 염색정규화과정에서는 이런 값들은 타깃이미지를 정함으로써, 결정되는데 이 논문에서는 이 수치들이 학습가능하게됩니다.

     - 단 완전히 학습값에만 의존하지 않도록 $\gamma_{k} = \gamma_{k}^{int} + \Delta \gamma_{k}^{opt}$로 분리합니다.

    - 그리고,  $ \Delta \gamma_{k}^{opt} $학습하여 조절할수 있도록 합니다. beta에 대해서도 동일히 작업합니다.

  2.3. 스케일링을 합니다. 정규화한 값이 너무 크거나 작거나하지않도록 [0, 1] 사이의 값으로 clip합니다. (식 9)

3. 어텐션을 이용한 컬러 공강 앙상블

  - 위의 과정을 컬러스페이스 s에 대해서 다양하게 진행합니다. 즉, 컬러공간(HSV, CEILAB 등)에 대해서 수행해서 Attention값을 구합니다.

  - Attention 값은 $Att(Q^{s}, K^{s}, V^{s})= softmax(\frac{Q^{s} \cdot   (K^{s})^{T}}{d})\cdot   V^{s}$ 로 계산합니다.

 - Self-attention을 계산할 때 사용되는 Q, K, V은 1x1 kernel을 2D conv하여 구합니다. 즉, 1x1으로 하니 채널축으로만 선형결합한 값입니다.

 - $Q^{s} = Flatten \circ  MaxPool \circ  Conv_{1}^{s}(\hat{ \mathbb{x}^{s}})$ 으로 컬러공간별로 convolution을 따로둡니다. (식11~12)

 - $K^{s} = Flatten \circ  MaxPool \circ  Conv_{2}^{s}(\hat{ \mathbb{x}^{s}})$: 마찬가지로 K도 컬러공간별로 convolution branch을 따로둡니다. 이 QK을 계산하면 patch내의 픽셀들사이의 global context을 얻는것입니다.

 - $V^{s} = Flatten \circ  MaxPool(\hat{ \mathbb{x}^{s}})$: Value은 풀링만하여 넣습니다.  위에서 계산한 global context와 곱합니다. 결과적으로 $(B, C, hw) \cdot (B, hw, hw)$을 하는모양입니다.

 - 이 결과를 원본 (B, C, W, H)와 곱해주기위해서 unflatten(B, C, w, h) 후에 업샘플하여 사이즈를 맞춰주고 (B, C, W, H) 원소별 곱하여 정규화된 이미지를 구합니다.

 - 이 이미지를 LAB, HSV 등 여러 이미지에 대해서 가중합을 하여 최종 정규화된 이미지를 획득합니다. 식(14)

요약

SimCLR에서는 Instance discrimiation이 학습목적으로 원본이미지를 서로 다른 증강방법을 이용해도, 서로 같은 임베딩이 되게끔 유도한다. CSI은 SimCLR방법과 매우 유사한데, 원본이미지를 증강한 경우만 OOD(Out Of Distribution)으로 학습하는 방법이다.

• 이미지 증강방법들의 집합 $S:=\{S_{0}, S_{i}, ..., S_{K-1}\}$
• 동일 이미지 반환: $I=S_{0}$

여러 이미지 증강방법S로부터 하나를 뽑아, 이미지 모든 이미지를 증강하고(원본반환 포함)이를 SimCLR을 돌림. 이 과정을 여러 증강방법에 대해서 반복함. 

추가적인 학습 테스크로, 증강된 이미지가, 어떤 증강방법으로 이용되었는지를 분류하는 방법도 진행

최종학습 Objective은 증강에 따른 이미지를 OOD Image로 SImCLR하는 목적과 증강방법론 분류모델을 혼합하여 예측하게 됨

• Diagnostic slide: FFPE(Formalin-fixed Paraffin-embedded)한 슬라이드. 
• Tissue slide: Frozen slide. 수술장에서 동TCGA에서 TS#, BS# 라고 표기

주의: Diagnostic slide여도 Fronzen 처리 후에 FFPE하는 경우 Frozen permanent라고 부르려 그냥 FFPE하는 경우보다 dissection 등의 artifact가 더 있을 수 있음

Motivation

1. Vanillar Attention MIL은 IID(independent and identical distribution)에 기반하고 있어서, 패치별 연관성을 고려하지 못함.

2. 이에 대한 대안으로 Transformer을 사용하는 경우도 계산량 때문에, 시퀀스 길이가 짧아야해서 WSI분석에는 적용이 어려움.

예를 들어, 아래의 그림과 같이 Vanillar Attention인 경우(d), $\sum \alpha_{i} * h_{i}$와 같이 attention weight * instance embedding으로 가중합으로 하게되는데, 이때 attention weight은 i의 인스턴스만이고려된 가중치임(=attention weight은 attention score로 유래되고, i 인스턴스만의 정보로만 계산됨). 달리말하면, 예측에 어떤 패치가 중요한지만을 가중치로 둠. 반면 Self-attention pooling을 하면, 국지적인 연관성도 하나의 정보로  확인해서, 이 국지적 연관성이 높다는 것을 하나의 가중치로 확인함.

Method:  TransMIL은 TPT Module로 이뤄짐

1. TPT Module의 중요 특징이 2가지

첫 째로, WSI 패치로 나눌 때, 패치가 N x N으로 연산되지 않고, $\sqrt{N} \times \sqrt{N}$으로 연산되도록 하는 Transformer구조를 가짐. 여기서 국지적인 연관성이 고려됨. * TPT Module에 대한 Full term이 논문내에서는 주어지지 아니함.

둘 째로, Convolution을 가변적으로 진행하기 위한 PPEG(Pyramid Position Encoding Generator). PPEG은 이미지를 N by N kernel 사이즈를 달리하여 Convolution을 적용함. 이 PPEG의 이점은 서로 다른 위치정보를 압축하고, 혼합할 수 있음.

TPT Module은 아래와 같은 알고리즘입. 크게는 1) 정사각형으로 시퀀싱을 만드는 단계, 2)멀티해드어텐션 단계, 3) PPEG적용단계), 4) 멑리해드어텐션 단계, 5) Slide level예측

그림으로 보면 단순함. 원래는 패치수가 30개면, 30 x 30의 Attention score (alignment score)을 구한후에 Attention weight을 계산하게 되는데, 여기서는 30x30이 아니라, 루트30 x 루트 30으로 계산함.  패치 수가 30개면, 이의 제곱근에 가장가까운 윗 정수는 6으로 계산됨. 즉, 이 행렬은 6x6의 사이즈를 만드는것이 목적임.

이렇게 하는 이유는 Attention score을 N x N으로 만들면 어짜피, 상삼각부분은 i,j의 패치와 i,j의 위치에 해당하는 연관성은 거의 동일할 것이기, 더 단순하게 $\sqrt{N} \times \sqrt{N}$으로 하려는 것으로 보임.

하지만, 이렇게 가공하면, 정방이 안나오기에, 정방이 나올수 있도록, 부족한 패치는 앞의 패치에서 따와 붙여줌. 그렇게 만든 패치의 모음이 $H_{s}^{f}$

이 방법을 이용해서 Sefl-attention을 구하면 됨. 하지만, 조금 더 나아가서 Nystrom Method라는 방법으로 O(n)의 연산을 가지도록 하고, 이렇게하면 토큰수가 엄청 커져도 효율적인 연산이 가능.

한편, PPEG은 위에서 만든 $H_{f}$을 여러 convolution kernel을 동시에 적용 후에,  padding하여 집계하는 방식으로 이뤄짐. 각 kernel을 통과한 featuremap은 단순히 합산으로 연산되고, 이를 Flattenting하여 다시 패치토큰으로 얻어냄.

Motivation

• 한 조직슬라이드를 보면, 같은 패턴(예, Stroma, Cancer, Epithelium)의 조직슬라이드 여러 번 등장한다는 것
• 여러 패턴의 특징을 집합으로해서, 분포로 특징화하면 슬라이드를 표현하기에 더 좋을 것 이라는 가설

Method

방법론은 크게 2가지로 구분됩니다. 1) 패치를 어떻게 프로토타입이랑 연계할 것인가, 2) 슬라이드 임베딩 과정

1. 패치를 임베딩: $z_{N_{j}}^{j} $

2. 슬라이드 임베딩: $z^{j}_{WSI} = \left[ \sum_{n=1}^{N_j} \varphi_j \left( z_{j \, n}, h_1 \right), \, \cdots , \, \sum_{n=1}^{N_j} \varphi_j \left( z_{j \, n}, h_C \right) \right]$

 - 패치의 임베딩인 집합을 패치의 개수가 더 작은 C개의 프로토타입의 집합으로 표현하고자합니다.프로토타입이란 개념을 이용합니다. 여기서 프로토타입은 각 패치의 특징을 대표할만한 벡터를 의미합니다. 예를 들어, 염증을 대표할만한 패치의 벡터, 암 패치를 대표할만한 벡터를 의미합니다. 

  -  이 패치 임베딩의 집합인 WSI은 수식으로 $Z_{WSI}^{j} \in \mathbb{R}^{C \cdot M}$ 으로 표현합니다. 이는 패치1부터 N까지를 프로토타입1에 대해서 유사도를 만들어 모두 합 합니다. 그리고 이 과정을 모든 프로토타입C까지를 진행하여, 벡터로 표현합니다. 결과적으로 M개의 벡터가 프로토타입개수인 C개까지하여 CM차원의 벡터를 얻을 수 있습니다. 슬라이드내 패치 수($N_{j}$)가 가변적이더라도, 이 축에 대해서 합계를 하기 때문에, $CM$의 차원으로 항상 고정인 벡터가 얻어집니다.

  - $\phi^{j}(\cdot, \cdot)$: 각 패치와 프로토타입의 유사도를 새로운 M차원으로 표현합니다. 이 매핑 함수를 정의하기위해서 GMM(Gausian mixture)을 이용합니다. 각 패치 임베딩($z_{n}^{j}$)가 GMM에서 얻어졌다라고 가정합니다. 즉, 패치임베딩은 여러 가우시안분포중 하나의 분포에서 생성될 확률이 있다고 봅니다. 

$p(x) = \sum_{k=1}^{K}\pi_{k} \cdot N(x|\mu_{x}, \sigma_{k})$

=> $p(z_{n_j}; \theta_j) = \sum_{c=1}^{C} p(c_{n_j} = c; \theta_j) \cdot p(z_{n_j} \mid c_{n_j} = c; \theta_j)$

=> $p(z_{n_j}; \theta_j) = \sum_{c=1}^{C} \pi_{c_j} \cdot N(z_{n_j}; \mu_{c_j}, \Sigma_{c_j})$

 - $ \pi $: 프로토타입c의 확률을 의미합니다. 여기서는 c번째 가우시안 분포가 선택될 확률을 나타냅니다. 이 확률의 총합은 1입니다.

 - $ \theta $ : 최적화해야할 파라미터를 의미합니다. GMM은 Mixutre probablity ($ \pi $), 평균($ \mu_{c} $), 공분산 행렬($ \Sigma_{c_j} $)의 집합을 의미합니다.

위의 수식을 풀어서 해석해보면, C개의 가우시안분포의확률 x 선택된 가우시안분포에서의 확률밀도 입니다.

1. 첫 번째 항($ p(c_{n_j} = c; \theta_j)  $)은 c번째 가우시안 분포로 선택될 확률을 의미합니다. 

2. 두 번째 항($ p(z_{n_j} \mid c_{n_j} = c; \theta_j) $): c번째 가우시안 분포일때, 확률밀도를 의미합니다. 예를 들어, 염증(c=염증)을 의미하는 가우시안 분포에서 임베딩이 나왔다면, 그 임베딩값이 나왔을 확률 밀도를 계산하는 것입니다.

여기서 중요한건, 혼합확률(Mixture probablity, $\pi _{c} ^{j}$)은 슬라이드 별로 계산합니다. 슬라이드의 특징을 고려해서, 특징을 뽑아서 WSI을 GMM의 파라미터로 설명합니다.

최종적으로, 아래와 같이 표현할 수 있습니다.

$ z_j^{\text{WSI}} = \left[z_j^{\text{WSI},1}, \cdots, z_j^{\text{WSI},C}\right] = \left[\pi b_j^1, \mu b_j^1, \Sigma b_j^1 \,  \cdots, \pi b_j^C, \mu b_j^C, \Sigma b_j^C \, \right]$

Q. GMM에 알고리즘에서 초기화는 어떻게하나?

1. 혼합확률(mixture probability, $\pi_{c} ^{j, (0)}$)은 균등분포로 1/C을 줍니다. 

2. 각 분포의 평균값($ \mu_{c} ^{j, (0)}$): c클레스 프로토타입에 대한 벡터

3. 공분산($Sigma_{c}^{j, (0)}$: Indentity matrix

4. 프로토타입 벡터에 대한 초기화: K-means clsutering을 훈련데이터셋에서 해서 

이미지에서 blur(블러, 흐림 현상)을 확인하는 데 사용되는 다양한 알고리즘과 지표들이 있습니다. 주로 이미지의 선명도, 엣지(경계)의 강도를 측정하여 블러를 평가합니다.

이 평가하는 방법은 2가지가 있습니다.

1. 원본과 대비가 가능한 경우: With reference image
2. 원본과 대비가 필요 없는 경우: without reference image

이 포스팅에서 without reference image(no refereince image)인 대표적인 지표 및 알고리즘들을 크게 분류하면 다음과 같습니다.

1. Spatial domain: 방법은 이미지의 픽셀과 인접 픽셀 간의 관계를 계산하여 흐림과 선명한 이미지를 구분하는 기법입니다.
   1. Grayscale Gradient base method: 이미지를 그레이스케일로 변환하여, 인접픽셀과의 그레디언트값(변화도)를 측정합니다. 측정된 변화도가 클수록 또렷한 이미지로 간주됩니다. 대표적으로 Laplacian variance가 있습니다.
2. Spectral-domain: 이미지의 고주파 성분과 저주파 성분을 분석하여 이미지의 선명도를 평가하는 방법입니다. 고주파 성분은 이미지의 선명한 부분과 관련이 있으며, 세부 정보와 경계 정보를 많이 포함하고 있습니다. 반면 저주파 성분은 흐릿한 부분에 해당합니다
3. Learning: 머신러닝을 이용한 방법
4. Combination: 2개 이상의 조합을 이용하는 방법

1. Laplacian Variance (라플라시안 분산)

라플라시안 분산은 라플라시안 커널을 이용해서, 이미지의 2차미분을 구해 분산을 구하는 방법입니다. 이 분산이 뜻하는 바는 픽셀의 흩어짐 정도로, 또렷한 이미지일수록 높은 분산값을 가지며, 흐릿한 이미지일수록 낮은 분산값을 가지게 됩니다. 특정한 임계점을 두어 또렷하다, 흐리다라고 판단할 수 있습니다. 

다음의 장점을 지닙니다.

• 간단하고 빠릅니다.: 계산이 비교적 간단하여 실시간 블러 감지에도 적합합니다.
• 효과적임: 엣지 정보를 기반으로 하므로 다양한 종류의 블러를 효과적으로 감지할 수 있습니다.
• 수학적 직관성: 분산을 이용한 접근 방식은 통계적으로도 타당성이 있습니다.

2차원 이미지에서는 다음과 같이 표현합니다. 2차미분을 구하기 위해서는 Laplacian kernel을 이용해서 구합니다. 

$\nabla^2 I = \frac{\partial^2 I}{\partial x^2} + \frac{\partial^2 I}{\partial y^2}$

이를 메트릭스 연산으로 표현하면, 라플라시안 커널이 됩니다.

$$\begin{bmatrix}
0 & 1 & 0 \\
1 & -4 & 1 \\
0 & 1 & 0
\end{bmatrix}
$$

import cv2
import numpy as np


def cal_laplacian_variance(image_array: np.ndarray):
    """라플라시안 분산을 구함

    Args:
        image_array (np.ndarray): RGB image array

    Returns:
        float: Laplacian variance
    """

    gray_image = cv2.cvtColor(image_array, cv2.COLOR_RGB2GRAY)
    laplacian = cv2.Laplacian(gray_image, cv2.CV_64F)
    return laplacian.var()

라플라시안 분산은 분산이기에 항상 0또는 양의값을 갖게됩니다. 이런 경계면이 모호한 이미지를 바로 식별할 수 있습니다.

또는 경계면이 또렷한(=focal plane)에 맞게 하는 경우의 이미지도 선택할 수 있습니다.

2.  Wavelet-based transform

이미지도 신호기이 때문에, 이를 신호로 처리하는 방법이 있습니다. Spectral-domain에 속하는 방법으로, 이미지를 저주파 및 고주파로 나눕니다. 각각 저주파 고주파는 다시, 저저, 고고, 저고,고저로 나눌 수 있고, 저저를 제외한 나머지 신호 강도를 이용하여 이미지의 선명도를 추출할 수 있습니다.

Blur Detection for Digital Images Using Wavelet Transform*이 논문에서도 저저주파를 제외하고 신호강도를 측정합니다.

def cal_wavelet(image_array):
    # 이미지를 그레이스케일로 변환 (이미지가 컬러일 경우)
    if len(image_array.shape) == 3:
        image = cv2.cvtColor(image_array, cv2.COLOR_RGB2GRAY)

    # 2D 웨이브렛 변환 수행 (Haar 웨이브렛 사용)
    coeffs = pywt.dwt2(image, "haar")  # 'haar' 외에도 다양한 웨이브렛 사용 가능
    LL, (LH, HL, HH) = coeffs

    # 고주파 성분 (LH, HL, HH)의 에너지를 합산하여 선명도 측정
    high_freq_energy = (
        np.sum(np.abs(LH) ** 2) + np.sum(np.abs(HL) ** 2) + np.sum(np.abs(HH) ** 2)
    )

    return high_freq_energy