[PYTHON] Implementierung der ML-EM-Methode, Querschnittsrekonstruktionsalgorithmus für CT-Scan

Einführung

Ich hatte die Möglichkeit, mich einem CT-Scan für medizinische Untersuchungen zu unterziehen. Ich habe mich gefragt, wie ich das Bild des Fehlers rekonstruieren kann, also werde ich es untersuchen und implementieren.

Der Bildrekonstruktionsalgorithmus für die Röntgen-CT hat traditionell eine filterkorrigierte Rückprojektion (FBP) verwendet, aber es scheint, dass die Bildqualität, die erhalten werden kann, wenn die Röntgendosis reduziert wird, begrenzt ist. Obwohl die Berechnungszeit länger als die von FBP ist, möchte ich daher die ML-EM-Methode ausprobieren, eine der statistischen Methoden, um klare Bilder zu erhalten.

Die ML-EM-Methode ist ein Algorithmus, der aus einem Bild, das auf der Grundlage der tatsächlichen Messprojektionsdaten rekonstruiert wurde, erneut Projektionsdaten erstellt und diese nacheinander aktualisiert, so dass sie sich dem tatsächlichen Messwert so weit wie möglich nähern.

Referenz

Ich bezog mich auf Folgendes.

Hiroyuki Shinohara, Grundlagen der Fehlerbildgebungsmethode 32. ML-EM-Methode und OS-EM-Methode

Toshiyuki Tanaka, Kommentar: Besonderheit: Grundlagen und neue Entwicklungen der Messtechnik unter Verwendung von Prinzipien zur Modellierung elektromagnetischer Felder / Aktueller Status der Röntgen-CT und weitere Verbesserung der Bildqualität / sicejl / 56/11 / 56_874 / _pdf)

Algorithmus

Zunächst werde ich den Ablauf des Algorithmus schreiben. Danach werde ich die Bedeutung der Wörter erklären.

1. Erstellen Sie die k-te Projektion aus dem k-ten Bild. Das erste Bild verwendet ein Bild mit allen 1 Werten
2. Ermitteln Sie das Verhältnis zwischen der k-ten Projektion und den gemessenen Projektionsdaten. Gemessene Projektionsdaten / (1) Projektionsdaten
3. Projizieren Sie das Verhältnis zurück
4. Multiplizieren Sie das k-te Bild mit dem zurückprojektierten Bild und aktualisieren Sie es auf das k + 1-te Bild.

Wir werden die Projektion und Rückprojektion erklären und dann die ML-EM-Methode implementieren.

Projektion von Messdaten

Die Messdaten sind ein Schatten, in dem ein Teil des Lichts vom menschlichen Körper durch Bestrahlung von Röntgenstrahlen aus einer Richtung absorbiert wird.

CT-Scans machen diese Projektion aus verschiedenen Winkeln. Hier habe ich die Daten erstellt, indem ich 360 Grad in 200 (in Schritten von 1,8 Grad) geteilt habe.

Das Programm, das die Messdaten vorbereitet hat, ist wie folgt.

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

def rotate_center(image, angle):
    """
Drehen Sie um die Bildmitte
    image :Bild
    angle :Drehwinkel[deg]
    """
    h, w = image.shape[:2]
    affine = cv2.getRotationMatrix2D((w/2.0, h/2.0), angle, 1.0) #Koordinaten der Bildmitte,Der Winkel, den Sie drehen möchten,Expansionsverhältnis
    return cv2.warpAffine(image, affine, (w, h))

def circle_mask(img):
    """
Maskieren Sie die anderen mit 0 und lassen Sie den Bereich des Beschriftungskreises frei
    """
    x = np.linspace(-1.0,1.0,img.shape[0])
    X,Y = np.meshgrid(x,x)
    distance_map = np.sqrt(X**2+Y**2) #Entfernung vom Zentrum
    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-Bild Das Ziel ist es, dies wiederherzustellen
true_img = circle_mask(true_img) #Kreismaske

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

#Es war ein Bild im Referenz-PDF
fig,axes = plt.subplots()
axes.imshow(forward_projection_array_meas,aspect=3)

In dem Artikel, auf den ich mich bezog, wurden die Projektionsdaten nebeneinander angeordnet und abgebildet, daher werde ich sie hier zum Vergleich veröffentlichen.

Rückprojektion

Die Rückprojektion ist eine Operation, bei der projizierte Daten in Kacheln so angeordnet werden, dass sie die gleichen Abmessungen wie das Originalbild haben, und dann addiert werden.

Implementierung der ML-EM-Methode

Es wird nach dem Algorithmus der ML-EM-Methode implementiert.

def rotate_center(img, angle):
    """
Drehen Sie um die Bildmitte
    image :Bild
    angle :Drehwinkel[deg]
    """
    h, w = img.shape[:2]
    affine = cv2.getRotationMatrix2D((w/2.0, h/2.0), angle, 1.0) #Koordinaten der Bildmitte,Der Winkel, den Sie drehen möchten,Expansionsverhältnis
    rot_img = cv2.warpAffine(img, affine, (w, h))
    return rot_img

def circle_mask(img):
    """
Maskieren Sie die anderen mit 0 und lassen Sie den Bereich des Beschriftungskreises frei
    """
    x = np.linspace(-1.0,1.0,img.shape[0])
    X,Y = np.meshgrid(x,x)
    distance_map = np.sqrt(X**2+Y**2) #Entfernung vom Zentrum
    img[distance_map>=1.0] = 0.0
    return img

def forward_projection(img, theta_array):
    """
Wiederholen Sie den Vorgang des Drehens und Projizierens des k-ten Bildes
    """
    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ückprojektion
    img_shape :Bildgröße
    theta_array :Winkeldaten
    forward_projection_array :Messdaten
    """
    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=In 0-Richtung verteilen
        tile_img = rotate_center(tile_img, -theta) #Drehung
        back_pro_img += tile_img
    back_pro_img /= theta_array.shape[0] #Teilen Sie durch die Anzahl der Überlagerungen
    return back_pro_img

Ursprünglich wurde die Wiederholung gestoppt, wenn der Unterschied zwischen dem k-1. Und dem k-ten Bild kleiner oder gleich wurde, aber der Einfachheit halber wurde er hier auf etwa das 60-fache eingestellt.

init_img = np.ones_like(true_img)
forward_projection_array_k = forward_projection(init_img, theta_array)+1.0e-12 # 1.1, um zu vermeiden, dass 0div die k-te Projektion aus dem k-ten Bild erstellt.0e-12 hinzugefügt
forward_pro_ratio = forward_projection_array_meas/forward_projection_array_k # 2.Finden Sie das Verhältnis zwischen der k-ten Projektion und den gemessenen Projektionsdaten. Gemessene Projektionsdaten/(1)Projektionsdaten
back_pro_img = back_projection(true_img.shape, theta_array, forward_pro_ratio) #3.Projizieren Sie das Verhältnis zurück
k_img = init_img * back_pro_img #4.Multiplizieren Sie das k-te Bild mit dem zurückprojektierten Bild und k+Aktualisiere auf das erste Bild

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

Ergebnis

Original Bild

1. Iteration

2. Wiederholung

Es ist etwas verschwommen als das Originalbild, aber es wird gut reproduziert.

Zusammenfassung

Wir haben die ML-EM-Methode implementiert, eine der CT-Bildrekonstruktionsmethoden, Messdaten aus Probenbildern erstellt und rekonstruiert.

Wenn die Koordinaten nach der Drehung nicht ganzzahlige Werte sind, muss die Helligkeit normalerweise in mehrere Pixel aufgeteilt werden. In der Referenz entspricht es der Variablen C in der Formulierung der EM-ML-Methode. Da die Drehung des openCV-Bildes durch das bilineare Verfahren interpoliert wird, wird die Gesamthelligkeit des Bildes vor und nach der Konvertierung nicht beibehalten. Ich denke jedoch, dass es für mich kein großes Problem ist, es als Hobby auszuprobieren.

Ich würde es gerne mit tatsächlichen Daten versuchen, kann aber kein Röntgen-CT-Gerät kaufen. Wenn Sie ein wenig Mehl mit Agar mischen und ein Bild aufnehmen, während Sie es mit Hintergrundbeleuchtung beleuchten, können Sie Projektionsdaten wie ein Röntgen-CT-Bild erhalten. Ich werde es tun, wenn ich eine Chance habe.