[PYTHON] Implementation of ML-EM method, cross-section reconstruction algorithm for CT scan
Introduction
I had the opportunity to have a CT scan for a medical examination. I was wondering how to reconstruct the image of the fault, so I will investigate and implement it.
The image reconstruction algorithm for X-ray CT has traditionally used filter-corrected back projection (FBP), but it seems that there is a limit to the image quality that can be obtained when the X-ray dose is reduced. Therefore, although the calculation time is longer than that of FBP, I would like to try the ML-EM method, which is one of the statistical methods for obtaining clear images.
The ML-EM method is an algorithm that creates projection data again from an image reconstructed based on the actually measured projection data, and sequentially updates it so that it approaches the actually measured value as much as possible.
reference
I referred to the following.
Hiroyuki Shinohara, Basics of Tomography Method 32nd ML-EM Method and OS-EM Method
Toshiyuki Tanaka, Commentary: Special Feature: Basics and New Developments of Measurement Technology Using Electromagnetic Field Modeling Principles / Current Status of X-ray CT and Further Improvement of Image Quality / sicejl / 56/11 / 56_874 / _pdf)
algorithm
First of all, I will write the flow of the algorithm. After that, I will explain the meaning of the words.
- Create the kth projection from the kth image. The first image uses an image with all 1 values
- Find the ratio between the kth projection and the measured projection data. Measured projection data / (1) projection data
- Back-project the ratio
- Multiply the kth image by the back-projected image and update to the k + 1th image.
We will explain projection and back projection, and then implement the ML-EM method.
Measurement data projection
The measurement data is a shadow in which a part of the light is absorbed by the human body by irradiating X-rays from one direction.
CT scans make this projection from different angles. Here, I created the data taken by dividing 360 degrees into 200 divisions (in 1.8 degree increments).
The program that prepared the measurement data is as follows.
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import cv2
def rotate_center(image, angle):
"""
Rotate around the center of the image
image :image
angle :Angle to rotate[deg]
"""
h, w = image.shape[:2]
affine = cv2.getRotationMatrix2D((w/2.0, h/2.0), angle, 1.0) #Coordinates of the center of the image,The angle you want to rotate,Expansion ratio
return cv2.warpAffine(image, affine, (w, h))
def circle_mask(img):
"""
Mask the others with 0, leaving the area of the inscribed circle
"""
x = np.linspace(-1.0,1.0,img.shape[0])
X,Y = np.meshgrid(x,x)
distance_map = np.sqrt(X**2+Y**2) #Distance from the center
img[distance_map>=1.0] = 0.0
return img
def forward_projection(img,theta):
rot_img = rotate_center(img, theta)
return np.sum(rot_img,axis=0)
theta_array = np.linspace(0.0,360.0,200,endpoint=False)
true_img = cv2.imread("CT_img.png ",cv2.IMREAD_GRAYSCALE) #CT image The goal is to restore this
true_img = circle_mask(true_img) #Circular mask
forward_projection_array_meas = []
for theta in theta_array:
forward_projection_array_meas.append(forward_projection(true_img,theta))
forward_projection_array_meas = np.asarray(forward_projection_array_meas)
fig,axes = plt.subplots()
for i in range(200):
axes.plot(forward_projection_array_meas[i])
#In the reference PDF, it was an image
fig,axes = plt.subplots()
axes.imshow(forward_projection_array_meas,aspect=3)
In the paper I referred to, the projection data was arranged side by side and imaged, so I will post it here for comparison.
Back projection
Back projection is the operation of laying out projected data in tiles so that it has the same dimensions as the original image, and adding them together.
Implementation of ML-EM method
It is implemented according to the algorithm of the ML-EM method.
def rotate_center(img, angle):
"""
Rotate around the center of the image
image :image
angle :Angle to rotate[deg]
"""
h, w = img.shape[:2]
affine = cv2.getRotationMatrix2D((w/2.0, h/2.0), angle, 1.0) #Coordinates of the center of the image,The angle you want to rotate,Expansion ratio
rot_img = cv2.warpAffine(img, affine, (w, h))
return rot_img
def circle_mask(img):
"""
Mask the others with 0, leaving the area of the inscribed circle
"""
x = np.linspace(-1.0,1.0,img.shape[0])
X,Y = np.meshgrid(x,x)
distance_map = np.sqrt(X**2+Y**2) #Distance from the center
img[distance_map>=1.0] = 0.0
return img
def forward_projection(img, theta_array):
"""
Repeat the operation of rotating and projecting the kth image
"""
img = circle_mask(img)
forward_projection_array = []
for theta in theta_array:
rot_img_sum = np.sum(rotate_center(img, theta),axis=0)
forward_projection_array.append(rot_img_sum)
forward_projection_array = np.asarray(forward_projection_array)
return forward_projection_array
def back_projection(img_shape, theta_array, forward_projection_array):
"""
Back projection
img_shape :Image size
theta_array :Angle data
forward_projection_array :Measurement data
"""
back_pro_img = np.zeros(img_shape,dtype=np.float64)
for i,theta in enumerate(theta_array):
tile_img = np.tile(forward_projection_array[i],[img_shape[0],1]) #axis=Spread in the 0 direction
tile_img = rotate_center(tile_img, -theta) #rotation
back_pro_img += tile_img
back_pro_img /= theta_array.shape[0] #Divide by the number of superpositions
return back_pro_img
Originally, when the difference between the k-1st and kth images became less than or equal to, the iteration was stopped, but here, for the sake of simplicity, it was set to about 60 times.
init_img = np.ones_like(true_img)
forward_projection_array_k = forward_projection(init_img, theta_array)+1.0e-12 # 1.1 to avoid 0div creating the kth projection from the kth image.0e-12 added
forward_pro_ratio = forward_projection_array_meas/forward_projection_array_k # 2.Find the ratio between the kth projection and the measured projection data. Measured projection data/(1)Projection data
back_pro_img = back_projection(true_img.shape, theta_array, forward_pro_ratio) #3.Back-project the ratio
k_img = init_img * back_pro_img #4.Multiply the k-th image by the back-projected image, and k+Update to the first image
for i in range(60):
forward_projection_array_k = forward_projection(k_img, theta_array)+1.0e-12
forward_pro_ratio = forward_projection_array_meas / forward_projection_array_k
back_pro_img = back_projection(true_img.shape, theta_array, forward_pro_ratio)
k_img *= back_pro_img
result
Original image
1st iteration
60th repetition
It's a little blurry than the original image, but it's well reproduced.
Summary
We implemented the ML-EM method, which is one of the CT image reconstruction methods, created measurement data from sample images, and reconstructed it.
Normally, when the coordinates after rotation are non-integer values, it is necessary to divide the brightness into multiple pixels. In the bibliography, it corresponds to the variable C in the formulation of the EM-ML method. Since the rotation of the image in openCV is interpolated by the bilinear method, the sum of the brightness of the image is not maintained before and after the conversion. However, I think it's not a big problem for me to try it as a hobby.
I would like to try it with actual data, but I can't buy an X-ray CT device. If you mix agar with a little flour and take a picture while illuminating it with a backlight, I think you can get projection data like an X-ray CT image. I will do it if I have a chance.
Recommended Posts