Maximum-Likelihood Expectation-Maximization

https://www.openaccessjournals.com/articles/image-reconstruction-for-petct-scanners-past-achievements-and-future-challenges-11017.html

 

PET에서의 MLEM은 Poisson 통계 모델에서의 MLE를 EM 알고리즘으로 푸는 것임.

• FBP에 비해, 느리지만 Poisson noise 억제에 강함
• OSEM, MAP-EM 등으로 속도 및 성능 개선 가능

2023.10.25 - [.../Math] - [Math] Poisson Distribution (포아송분포)

[Math] Poisson Distribution (포아송분포)

---

Poisson Likelihood 모델 설정

각 검출기 $i$에서의 측정값 $y_i$ 는 포아송 분포를 따른다고 가정:

$$y_i \sim \text{Poisson} \left( \hat{y}_i = \sum_j p_{ij} \lambda_j \right)$$

이때 전체 log-likelihood 함수는 다음과 같음:

$$\log \mathcal{L}(\boldsymbol{\lambda}) = \sum_i \left[ y_i \log \left( \sum_j p_{ij} \lambda_j \right) - \left( \sum_j p_{ij} \lambda_j \right) - \log y_i! \right]$$

참고로 $\log y_i!$는 $\lambda_j$ 에 독립이므로 생략 가능함.

이를 통해 다음의 miximization problem으로 정의할 수 있음:

$$\underset{\boldsymbol{\lambda} \geq 0}{\text{maximize}} \quad Q(\boldsymbol{\lambda}) = \sum_i \left[ y_i \log \left( \sum_j p_{ij} \lambda_j \right) - \sum_j p_{ij} \lambda_j \right]$$

• 이는 utility function이므로 최대화를 수행.

EM 방식으로는 다음 surrogate function을 최대화!: 이 과정은 Jensen's inequality에 기반

$$\lambda_j^{k+1} = \lambda_j^k \cdot \frac{1}{\sum_i p_{ij}} \sum_i p_{ij} \cdot \left( \frac{y_i}{\sum_{j'} p_{ij'} \lambda_{j'}^k} \right)$$

요약하면 다음과 같음:

• objective function : $\log \mathcal{L}(\boldsymbol{\lambda})$
• 최적화 방법: Expectation-Maximization (EM)
• solution update equation: multiplicative update (양수 유지 보장)
• 정규화 항: $\sum_i p_{ij}$ - sensitivity normalization

2022.06.02 - [.../Math] - [ML] Likelihood (우도, 기대값)

[ML] Likelihood (우도, 기대값)

---

Mimum Likehood Expection Miximization 수식 **

$$
\lambda_j^{k+1} = \frac{\lambda_j^k}{\sum_i p_{ij}} \sum_i p_{ij} \cdot \left( \frac{y_i}{\sum_{j'} p_{ij'} \lambda_{j'}^k} \right)
$$

where,

• $\lambda_j^k$ : $k$번째 iteration(반복)에서 픽셀 $j$에 대한 이미지 추정값
• $y_i$ : $i$번째 검출기에서의 측정값 (projection bin data, sinogram)
• $p_{ij}$ : system matrix(시스템 행렬)의 요소. 픽셀 $j$가 검출기 $i$에 기여하는 정도
• $\sum_{j'} p_{ij'} \lambda_{j'}^k $ : $i$ 검출기에 대한 예측값 - $\hat{y}_i$

---

MLEM - Iterative Reconstruction

위의 수식은 MLEM 알고리즘의 업데이트 식임:

1.예측 측정값 계산 (forward projection):

$$\hat{y}_i = \sum_{j'} p_{ij'} \lambda_{j'}^k$$

2.실제와 예측의 비율 계산:

$$\frac{y_i}{\hat{y}_i}$$

3.백프로젝션 수행:

$$\sum_i p_{ij} \cdot \left( \frac{y_i}{\hat{y}_i} \right)$$

4.정규화 후 업데이트:
$$\lambda_j^{k+1} = \lambda_j^k \cdot \frac{1}{\sum_i p_{ij}} \sum_i p_{ij} \cdot \left( \frac{y_i}{\hat{y}_i} \right)$$

 

https://gist.github.com/dsaint31x/5a46553d69055b93dcacc249c043bef5

pet_mlem_simul.ipynb

---

Note:

• Poisson likelihood 기반임! (측정치는 counter의 결과이므로 항상 양의 정수로)
• 분자: 측정값과 예측값의 비율
• 분모: 픽셀별 시스템 응답의 총합으로 정규화
• 음수 없음, 수렴 보장 (단, 느림)
• forward / backward projection 반복 구조

---

같이보면 좋은 자료

2025.09.01 - [정리필요./PET, MRI and so on.] - Radon Transform and Inverse Radon Transform-FBP

Radon Transform and Inverse Radon Transform-FBP