[PYTHON] Implémentation de la méthode ML-EM, algorithme de reconstruction en coupe pour scanner

introduction

J'ai eu l'occasion de subir un scanner pour des examens médicaux. Je me demandais comment reconstruire l'image de la faute, je vais donc l'examiner et la mettre en œuvre.

L'algorithme de reconstruction d'image pour la tomodensitométrie aux rayons X a traditionnellement utilisé la rétroprojection corrigée par filtre (FBP), mais il semble qu'il y ait une limite à la qualité d'image qui peut être obtenue lorsque la dose de rayons X est réduite. Par conséquent, bien que le temps de calcul soit plus long que celui de FBP, j'aimerais essayer la méthode ML-EM, qui est l'une des méthodes statistiques pour obtenir des images claires.

La méthode ML-EM est un algorithme qui crée à nouveau des données de projection à partir d'une image reconstruite sur la base des données de projection de mesure réelles, et les met à jour séquentiellement afin qu'elle se rapproche le plus possible de la valeur de mesure réelle.

référence

J'ai évoqué ce qui suit.

Hiroyuki Shinohara, Basics of Fault Imaging Method 32nd ML-EM Method and OS-EM Method

Toshiyuki Tanaka, Commentaire: Dossier spécial: Bases et nouveaux développements de la technologie de mesure utilisant les principes de modélisation du champ électromagnétique / État actuel de la tomodensitométrie à rayons X et amélioration de la qualité de l'image / sicejl / 56/11 / 56_874 / _pdf)

algorithme

Tout d'abord, j'écrirai le flux de l'algorithme. Après cela, j'expliquerai le sens des mots.

1. Créez la kième projection à partir de la kième image. La première image utilise une image avec les 1 valeurs
2. Trouvez le rapport entre la k-ième projection et les données de projection mesurées. Données de projection mesurées / (1) données de projection
3. Rétroprojection du ratio
4. Multipliez la kème image par l'image rétroprojectée et mettez à jour la k + 1ème image.

Nous allons expliquer la projection et la rétroprojection, puis implémenter la méthode ML-EM.

Projection des données de mesure

Les données de mesure sont une ombre dans laquelle une partie de la lumière est absorbée par le corps humain en irradiant des rayons X dans une direction.

Les tomodensitogrammes font cette projection sous différents angles. Ici, j'ai créé les données prises en divisant 360 degrés en 200 (par incréments de 1,8 degrés).

Le programme qui a préparé les données de mesure est le suivant.

import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import cv2

def rotate_center(image, angle):
    """
Rotation autour du centre de l'image
    image :image
    angle :Angle de rotation[deg]
    """
    h, w = image.shape[:2]
    affine = cv2.getRotationMatrix2D((w/2.0, h/2.0), angle, 1.0) #Coordonnées du centre de l'image,L'angle que vous souhaitez faire pivoter,Taux d'expansion
    return cv2.warpAffine(image, affine, (w, h))

def circle_mask(img):
    """
Masquez les autres avec 0 en quittant la zone du cercle inscrit
    """
    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 du centre
    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) #Image CT Le but est de restaurer cette
true_img = circle_mask(true_img) #Masque circulaire

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])

#C'était une image dans le PDF de référence
fig,axes = plt.subplots()
axes.imshow(forward_projection_array_meas,aspect=3)

Dans l'article auquel j'ai fait référence, les données de projection ont été disposées côte à côte et imagées, je vais donc les publier ici pour comparaison.

Rétroprojection

La rétroprojection est une opération dans laquelle les données projetées sont disposées en mosaïques de manière à avoir les mêmes dimensions que l'image d'origine, puis ajoutées ensemble.

Mise en œuvre de la méthode ML-EM

Il est implémenté selon l'algorithme de la méthode ML-EM.

def rotate_center(img, angle):
    """
Rotation autour du centre de l'image
    image :image
    angle :Angle de rotation[deg]
    """
    h, w = img.shape[:2]
    affine = cv2.getRotationMatrix2D((w/2.0, h/2.0), angle, 1.0) #Coordonnées du centre de l'image,L'angle que vous souhaitez faire pivoter,Taux d'expansion
    rot_img = cv2.warpAffine(img, affine, (w, h))
    return rot_img

def circle_mask(img):
    """
Masquez les autres avec 0 en quittant la zone du cercle inscrit
    """
    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 du centre
    img[distance_map>=1.0] = 0.0
    return img

def forward_projection(img, theta_array):
    """
Répétez l'opération de rotation et de projection de la kème 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):
    """
Rétroprojection
    img_shape :Taille de l'image
    theta_array :Données d'angle
    forward_projection_array :Données de mesure
    """
    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=Propagation dans la direction 0
        tile_img = rotate_center(tile_img, -theta) #rotation
        back_pro_img += tile_img
    back_pro_img /= theta_array.shape[0] #Divisez par le nombre de superpositions
    return back_pro_img

À l'origine, lorsque la différence entre la k-1ère et la kième image devenait inférieure ou égale à, la répétition était arrêtée, mais ici, par souci de simplicité, elle était réglée à environ 60 fois.

init_img = np.ones_like(true_img)
forward_projection_array_k = forward_projection(init_img, theta_array)+1.0e-12 # 1.1 pour éviter que 0div crée la kème projection à partir de la kème image.0e-12 ajouté
forward_pro_ratio = forward_projection_array_meas/forward_projection_array_k # 2.Trouvez le rapport entre la kième projection et les données de projection mesurées. Données de projection mesurées/(1)Données de projection
back_pro_img = back_projection(true_img.shape, theta_array, forward_pro_ratio) #3.Rétro-projeter le ratio
k_img = init_img * back_pro_img #4.Multipliez la k-ème image par l'image rétroprojectée, et k+Mettre à jour vers la première 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

résultat

Image originale

1ère itération

60e répétition

Elle est un peu floue que l'image d'origine, mais elle est bien reproduite.

Résumé

Nous avons mis en œuvre la méthode ML-EM, qui est l'une des méthodes de reconstruction d'images CT, créé des données de mesure à partir d'images d'échantillons et les avons reconstruites.

Normalement, lorsque les coordonnées après rotation sont des valeurs non entières, il est nécessaire de diviser la luminosité en plusieurs pixels. Dans la référence, il correspond à la variable C dans la formulation de la méthode EM-ML. Puisque la rotation de l'image openCV est interpolée par la méthode bilinéaire, la luminosité totale de l'image n'est pas maintenue avant et après la conversion. Cependant, je pense que ce n'est pas un gros problème pour moi de l'essayer comme passe-temps.

J'aimerais l'essayer avec des données réelles, mais je ne peux pas acheter un appareil de tomodensitométrie à rayons X. Si vous mélangez un peu de farine avec de la gélose et que vous photographiez en éclairant avec un rétroéclairage, je pense que vous pouvez obtenir des données de projection comme une image CT aux rayons X. Je le ferai si j'ai une chance.