Grundlegende ITK-Verwendung mit Python gelernt

Notification

Ich habe diesen Artikel als Memo meiner Studie geschrieben. Die Referenzquelle ist dieses Notizbuch. https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization

ITK Eine Bibliothek, die für die Bildverarbeitung im biomedizinischen Bereich unverzichtbar ist und auf Bildvisualisierungstechnologie spezialisiert ist. Es wird in Slicer und ImageJ verwendet.

Umgebung

Google Colab (Starten Sie ein beliebiges Notizbuch.)

Beispieldaten

Laden Sie zuerst die Beispieldaten mit Subversion von GitHub herunter.

!sudo apt install subversion
#Holen Sie sich die Daten
#In der Github-URL enthalten/tree/Meister/Durch Kofferraum ersetzen
!svn checkout https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization/trunk/data

ITK-Installation

#Installieren Sie ITK
!pip install itk

Image Filtering

Was hier zu lernen

• Informieren Sie sich über das von ITK verwendete Bildverarbeitungs-Pipeline-Modell. --Bild und Netz werden zur Darstellung von Datenobjekten verwendet.
• Während der Bearbeitung des Datenobjekts wird ein neues Datenobjekt erstellt.
• Prozesse werden als Prozessobjekte kategorisiert und verfügen über Quell-, Filter- und Mapperobjekte. --source (Leser usw.) ruft die Daten ab und filtert die Objekte, um neue Daten zu generieren. --mapper akzeptiert Daten zur Ausgabe in Dateien oder andere Systeme.

Der Begriff Filter wird häufig in drei Typen verwendet:

--Datenpipeline Daten- und Prozessobjekte werden normalerweise mit den Methoden SetInput () und GetOutput () ineinander konvertiert. Es werden keine neuen Ausgaben oder neuen Pixeldaten für eine Verarbeitung erstellt, bis die Update () -Methode am Ende der Pipeline (Prozess oder Datenobjekt) aufgerufen wird.

• Pipeline-Updates Mehrdimensionale Bilddaten sind groß und können noch größer gemacht werden. ITK behebt dieses Problem durch seine Daten-Streaming-Funktionen.
• Streamen Das Streaming erfolgt durch Aufteilen des Bildes in nicht überlappende Bereiche am Ende der Pipeline. Verbreiten Sie dann die Pipeline als angeforderte Region.

Es gibt drei Arten von Bildregionen in ITK.

--BufferedRegion: Bereich der im Speicher gespeicherten Pixel --LargestPossibleRegion: Maximal möglicher Bereich des Bildes --RequestedRegion: Region, die während des Streamings von einem einzelnen Verarbeitungspfad angefordert wird (BufferedRegion und LargestPossibleRegion sind gleich oder größer als RequestedRegion).

Versuchen.

import numpy as np
import itk
import matplotlib.pyplot as plt

#Üben Sie die räumliche Filterung mit Pacman
fileName = 'data/PacMan.png'
reader = itk.ImageFileReader.New(FileName=fileName)
#Versuchen Sie es mit einem Gaußschen Filter
smoother = itk.RecursiveGaussianImageFilter.New(Input=reader.GetOutput())

print("reader's Output: %s" % reader.GetOutput())
print("smoother's Output: %s" % smoother.GetOutput())

#Aktualisieren, um Verarbeitungsergebnisse zu generieren()Vergiss nicht.
smoother.Update()

print("reader's Output: %s" % reader.GetOutput())
print("smoother's Output: %s" % smoother.GetOutput())

image = reader.GetOutput()#Bestätigung des Originals
plt.subplot(1,2,1),plt.title("original"),plt.imshow(itk.GetArrayFromImage(image))
smoothed = smoother.GetOutput()#Nach dem Glätten
plt.subplot(1,2,2),plt.title("smoothing"),plt.imshow(itk.GetArrayFromImage(smoothed))

Wenn Sie mit einem Pipeline-Objekt aktualisiert haben (obwohl Sie dies bereits getan haben), ist das Verarbeitungsergebnis das neueste

smoother.Update()
#Wenn Sie den Prozess jedoch ändern, müssen Sie erneut aktualisieren.
smoother.SetSigma(10.0)
smoother.Update()
plt.imshow(itk.GetArrayFromImage(smoothed))

Ich glaube auch nicht, dass es viele gibt, aber wenn Sie in der Mitte Änderungen am Reader vornehmen, vergessen Sie nicht, den Reader zu aktualisieren.

# reader.Modified()

#Sie können verschiedene Filter ausprobieren.
# https://itk.org/Doxygen/html/group__IntensityImageFilters.html
fileName = 'data/PacMan.png'
reader = itk.ImageFileReader.New(FileName=fileName)
smoother = itk.MedianImageFilter.New(Input=reader.GetOutput())#Zum Beispiel Median
smoother.SetRadius(5)
smoother.Update()
plt.imshow(itk.GetArrayFromImage(smoother.GetOutput()))

Ergänzung

#Platzieren Sie StreamingImageFilter am Ende der Pipeline, um die Pipeline zu streamen.
#Smoother erzeugt mehrmals eine Ausgabe für jede Streaming-Abteilung des Bildbereichs.
#Der Reader kann nicht gestreamt werden, daher wird die Ausgabe nur einmal generiert.
streamer = itk.StreamingImageFilter.New(Input=smoother.GetOutput())
streamer.SetNumberOfStreamDivisions(4)
reader.Modified()
streamer.Update()
#Sie müssen nur die Schrift einschließlich der Erweiterung angeben
writer = itk.ImageFileWriter.New(Input=smoother.GetOutput())
writer.SetFileName('my_output.mha')
writer.SetNumberOfStreamDivisions(4)
writer.Update()

Image Segmentation

Was hier zu lernen

--Kennen Sie die in ITK verfügbaren Segmentierungsmethoden (First, ConfidenceConnectedImageFilter)

• Verstehen Sie, wie das Leuchten der Region und ihre Parameter das Verhalten ändern
• Verstehen Sie, wie Level-Sets und ihre Parameter ihr Verhalten ändern

In ITK verfügbare Segmentierungsarten

• Level Set
• Wasserscheide (Gradiententropfen von Bildmerkmalen und Analyse entlang der Bereichsgrenzen)
• Verbundene Komponente (eine abstrakte Technik zum Arbeiten mit verbundenen Objekten, die bei der Binärisierung und Beschriftung verwendet werden)
• Label-Voting (zum Zählen von Pixelwerten als Labels und deren Verwendung im Prozess)
• Region Growing --Klassifizierer (Algorithmus für maschinelles Lernen)

Hier stellen wir ein Verarbeitungsbeispiel für Region Growing und Level Set vor.

Region growing

Der grundlegende Ansatz des Regionserweiterungsalgorithmus besteht darin, die Expansion von der Startregion (normalerweise einem oder mehreren Pixeln) aus zu starten, die als innerhalb des segmentierten Objekts betrachtet wird. Pixel neben diesem segmentierten Objektbereich werden ausgewertet, um festzustellen, ob sie als Teil des Objekts betrachtet werden sollen. Wenn sie segmentiert sind, werden sie dem Bereich hinzugefügt und der Vorgang wird fortgesetzt, solange dem Bereich neue Pixel hinzugefügt werden. Der Regionserweiterungsalgorithmus ist das Kriterium, das verwendet wird, um zu bestimmen, ob eine Region Pixel enthält, die Art der Verbindung, die zum Bestimmen von Nachbarn verwendet wird, und der Prozess (oder Prozess), der zum Suchen nach benachbarten Pixeln verwendet wird. Strategie) hängt ab von.

Hier wird die folgende Methode als Beispiel verwendet.

• ConfidenceConnected
• ConnectedThreshold
• IsolatedConnected

Überprüfen Sie zunächst, wie "Confidence Connected" verwendet wird. ConfidenceConnected ist eine Methode zum Bestimmen der oberen und unteren Grenzen der Bereichsgrenze, zum Ermitteln der Standardabweichung von den Werten zwischen ihnen und zum Festlegen des zulässigen Bereichs dieser Standardabweichung (wie Hebelwirkung. SD * Multiplikator).

# BrainProtonDensitySlice.Verwenden Sie png
# 'unsigned char'Der Grund für das Laden mit besteht darin, die Datenkonvertierungsmethode später anzuzeigen.
input_image = itk.imread('data/BrainProtonDensitySlice.png', itk.ctype('unsigned char'))
# View the input image
plt.imshow(itk.GetArrayFromImage(input_image))

Überprüfen Sie, wie die Klasse "ConfidenceConnectedImageFilter" verwendet wird. Wenn Sie nicht wissen, wie man es benutzt, schauen Sie nach.

confidence_connected = itk.ConfidenceConnectedImageFilter.New(input_image)
#Überprüfen Sie die Verwendung mit der Hilfefunktion von Python
help(confidence_connected)

#Stellen Sie den vertrauenswürdigen Filterparameter beliebig ein
confidence_connected.SetMultiplier(2.3)#Toleranz ± SD in diesem Beispiel*2.3
confidence_connected.SetNumberOfIterations(5)
confidence_connected.SetInitialNeighborhoodRadius(3)
confidence_connected.SetReplaceValue(255)

Wenn Sie nicht wissen, wie eine Funktion verwendet wird, schlagen Sie sie nach.

# *Multiplier*?
confidence_connected.SetMultiplier?

Stellen Sie den Startpunkt "Startpunkt" der Flächenerweiterung ein. Beachten Sie, dass sich der Ursprung des ITK-Bildobjekts unten links in den Bildkoordinaten befindet. (NumPy ist oben links)

# Set the seed points
confidence_connected.SetSeed([100, 110])
#Ausführen und anzeigen
confidence_connected.Update()
plt.imshow(itk.GetArrayFromImage(confidence_connected))

#Versuchen Sie, das Verarbeitungsergebnis zu verbessern. Wenden Sie einen Medianfilter auf das Originalbild an.
median_filtered_image = itk.median_image_filter(input_image, radius=1)
#Aktualisieren Sie die Parameter erneut. Geben Sie das Bild nach dem Medianfilter ein.
confidence_connected_image = itk.confidence_connected_image_filter(median_filtered_image,
                                                                   seed=[100,110],
                                                                   replace_value=255,
                                                                   multiplier=3.0,
                                                                   number_of_iterations=5,
                                                                   initial_neighborhood_radius=3)
plt.imshow(itk.GetArrayFromImage(confidence_connected_image))

Als nächstes vergleichen wir die drei Glühprozesse der Region. ConnectedThresholdImageFilter vs ConfidenceConnectedmageFilter vs IsolatedConnectedImageFilter

# ConnectedThresholdImageFilter vs ConfidenceConnectedmageFilter vs IsolatedConnectedImageFilter
input_image = itk.imread('data/BrainProtonDensitySlice.png', itk.ctype('unsigned char'))
connected_threshold = itk.ConnectedThresholdImageFilter.New(input_image)
connected_threshold.SetLower(190)
connected_threshold.SetUpper(255)
connected_threshold.SetReplaceValue( 255 )
connected_threshold.SetSeed( [100, 110] );
connected_threshold.Update()

confidence_connected = itk.ConfidenceConnectedImageFilter.New(input_image)
confidence_connected.SetMultiplier(2.3)
confidence_connected.SetSeed( [100, 110] );
confidence_connected.AddSeed( [80, 126] );
confidence_connected.SetNumberOfIterations(5)
confidence_connected.SetInitialNeighborhoodRadius(3)
confidence_connected.SetReplaceValue(255)
confidence_connected.Update()

isolated_connected = itk.IsolatedConnectedImageFilter.New(input_image)
isolated_connected.SetSeed1([98, 112])
isolated_connected.SetSeed2([98, 136])
isolated_connected.SetUpperValueLimit( 245 )
isolated_connected.FindUpperThresholdOff();
isolated_connected.Update()

plt.subplot(131)
plt.title("ConnectedThreshold")
plt.imshow(itk.GetArrayFromImage(connected_threshold))
plt.subplot(132)
plt.title("ConfidenceConnected")
plt.imshow(itk.GetArrayFromImage(confidence_connected))
plt.subplot(133)
plt.title("IsolatedConnected")
plt.imshow(itk.GetArrayFromImage(isolated_connected))
plt.tight_layout()
plt.show()

Fast Marching Segmentation(LevelSet) Versuchen Sie als Nächstes die Fast Marching-Segmentierung der LevelSet-Methode. Wenn der Zielbereich der Segmentierung eine einfache Struktur ist, kann ein schneller Expansionsalgorithmus verwendet werden, der als schnelles Marschieren bezeichnet wird. Vereinfachen Sie das Bild mit einem räumlichen Filter, bevor Sie es versuchen.

#Finden Sie heraus, welche Bilddatentypen der Glättungsfilter verarbeiten kann
itk.CurvatureAnisotropicDiffusionImageFilter.GetTypes()

#So wandeln Sie Bilddaten in einen unterstützten Datentyp um
Dimension = input_image.GetImageDimension()
FloatImageType = itk.Image[itk.ctype('float'), Dimension]
caster = itk.CastImageFilter[input_image, FloatImageType].New()
caster.SetInput(input_image)

#Wenden Sie einen Glättungsfilter an.
smoothed_image = itk.curvature_anisotropic_diffusion_image_filter(caster, 
                                                                  time_step=0.125, 
                                                                  number_of_iterations=5,
                                                                  conductance_parameter=3.0)

plt.imshow(itk.GetArrayFromImage(smoothed_image))

#Führen Sie den Kantenerkennungsfilter so aus, wie er ist
gradient = itk.gradient_magnitude_recursive_gaussian_image_filter(smoothed_image, sigma=1.0)
plt.imshow(itk.GetArrayFromImage(gradient))

#Führen Sie den Sigmoidfilter weiter aus

basin_minimum = 2.25
border_minimum = 3.75

alpha = - (border_minimum - basin_minimum) / 3.0
beta = (border_minimum + basin_minimum) / 2.0

sigmoid = itk.sigmoid_image_filter(gradient, output_minimum=0.0, output_maximum=1.0, alpha=alpha, beta=beta)
plt.imshow(itk.GetArrayFromImage(sigmoid))

#Hochgeschwindigkeitsmarsch
fast_marching = itk.FastMarchingImageFilter.New(sigmoid)

seed_value = 0.0
NodeType = itk.LevelSetNode[itk.ctype('float'), Dimension]
node = NodeType()
node.SetIndex([100,110])
node.SetValue(seed_value)

NodeContainerType = itk.VectorContainer[itk.ctype('unsigned int'), NodeType]
nodes = NodeContainerType.New()
nodes.Initialize()
nodes.InsertElement(0, node)

fast_marching_image = itk.fast_marching_image_filter(sigmoid, trial_points=nodes, stopping_value=80)
plt.imshow(itk.GetArrayFromImage(fast_marching_image))

Binarisieren Sie das Bild.

time_threshold = 100 #Schwellenwert einstellen
thresholded_image = itk.binary_threshold_image_filter(fast_marching_image,
                                                      lower_threshold = 0,
                                                      upper_threshold=time_threshold,
                                                      outside_value=0,
                                                      inside_value=255)
plt.imshow(itk.GetArrayFromImage(thresholded_image))

Versuchen Sie von hier aus Shape Detection Level Set. Die Verarbeitungsgeschwindigkeit wird verbessert, aber die Krümmung ist nicht gut. Die Pipeline sieht folgendermaßen aus:

Erstellen Sie zunächst ein ausgeglichenes Bild im unteren Teil der Abbildung.

#Das Binärbild ist Weiße Materie.Verwenden Sie png
binary_mask = itk.imread('data/WhiteMatter.png', itk.ctype('float'))

smoothed_image = itk.smoothing_recursive_gaussian_image_filter(binary_mask)

#Abstufungsverarbeitung
rescaler = itk.IntensityWindowingImageFilter[FloatImageType,FloatImageType].New(smoothed_image)
rescaler.SetWindowMinimum(0)
rescaler.SetWindowMaximum(255)
#Invertieren Sie den Eingangspegel, um intern einen negativen Wert zu haben
rescaler.SetOutputMinimum(5.0)
rescaler.SetOutputMaximum(-5.0)

rescaler.Update()
plt.imshow(itk.GetArrayFromImage(rescaler))

Erstellen Sie als Nächstes das Feature-Image in der Abbildung.

#Verwenden Sie das geladene Bild wie oben beschrieben vor
gradient = itk.gradient_magnitude_recursive_gaussian_image_filter(caster, sigma=1.0)

basin_minimum = 1
border_minimum = 5.0

alpha = - (border_minimum - basin_minimum) / 3.0
beta = (border_minimum + basin_minimum) / 2.0

#Funktionsbild in der Abbildung
sigmoid = itk.sigmoid_image_filter(gradient, output_minimum=0.0, output_maximum=1.0, alpha=alpha, beta=beta)
plt.imshow(itk.GetArrayFromImage(sigmoid))

#shape_detection_level_Set-Verarbeitung ausführen
shape_detection_image = itk.shape_detection_level_set_image_filter(rescaler,
                                                                   feature_image=sigmoid,
                                                                   maximum_r_m_s_error=0.001,
                                                                   number_of_iterations=100,
                                                                   propagation_scaling=1.0,
                                                                   curvature_scaling=2.0)
#Binarisierung
thresholded_image = itk.binary_threshold_image_filter(shape_detection_image,
                                                      lower_threshold=-1e7,
                                                      # Cut at the zero set
                                                      upper_threshold=0.0,
                                                      outside_value=0,
                                                      inside_value=255)
plt.imshow(itk.GetArrayFromImage(thresholded_image))

Lassen Sie uns einige Parameter in Bezug auf die Krümmungs- und Ausbreitungsbedingungen des eingestellten Pegels in den eingestellten Formerkennungspegel einführen. Min Basin Min Boundary Propagation Scaling Curvature Scaling Number of Iterations

def explore_shape_detection_image_parameters(basin_minimum, boundary_minimum, propagation_scaling,
                                             curvature_scaling, number_of_iterations):

    input_image = itk.imread('data/BrainProtonDensitySlice.png', itk.ctype('unsigned char'))

    # Prepare the initial level set image
    binary_mask = itk.imread('data/WhiteMatter.png', itk.ctype('float'))

    smoothed_image = itk.smoothing_recursive_gaussian_image_filter(binary_mask)

    Dimension = input_image.GetImageDimension()
    FloatImageType = itk.Image[itk.ctype('float'), Dimension]
    caster = itk.CastImageFilter[input_image, FloatImageType].New()
    caster.SetInput(input_image)

    # Prepare the speed image
    gradient = itk.gradient_magnitude_recursive_gaussian_image_filter(caster, sigma=1.0)

    alpha = - (boundary_minimum - basin_minimum) / 3.0
    beta = (boundary_minimum + basin_minimum) / 2.0

    sigmoid = itk.sigmoid_image_filter(gradient, output_minimum=0.0, output_maximum=1.0, alpha=alpha, beta=beta)

    rescaler = itk.IntensityWindowingImageFilter[FloatImageType,FloatImageType].New(smoothed_image)
    rescaler.SetWindowMinimum(0)
    rescaler.SetWindowMaximum(255)
    # Invert the input level set to have negative values inside
    rescaler.SetOutputMinimum(5.0)
    rescaler.SetOutputMaximum(-5.0)

    rescaler.Update()
    
    shape_detection_image = itk.shape_detection_level_set_image_filter(rescaler,
                                                                       feature_image=sigmoid,
                                                                       maximum_r_m_s_error=0.001,
                                                                       number_of_iterations=number_of_iterations,
                                                                       propagation_scaling=propagation_scaling,
                                                                       curvature_scaling=curvature_scaling)

    thresholded_image = itk.binary_threshold_image_filter(shape_detection_image,
                                                          lower_threshold=-1e7,
                                                          # Cut at the zero set
                                                          upper_threshold=0.0,
                                                          outside_value=0,
                                                          inside_value=255)
    plt.imshow(itk.GetArrayFromImage(thresholded_image))
    plt.title(str(basin_minimum)+"-"+str(boundary_minimum)+"-"+str(propagation_scaling)+"-"+str(curvature_scaling)+"-"+str(number_of_iterations))
    plt.show()
    return

explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=5., 
                                         propagation_scaling=1., curvature_scaling=2., number_of_iterations=100)
explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=5., 
                                         propagation_scaling=1., curvature_scaling=2., number_of_iterations=500)
explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=5., 
                                         propagation_scaling=1., curvature_scaling=4., number_of_iterations=100)
explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=20., 
                                         propagation_scaling=1., curvature_scaling=1., number_of_iterations=100)

Image Registration

Was hier zu lernen

• Wissen, wie ITK Berechnungen im physischen Raum durchführt
• Verstehen, warum ITK die Registrierung im physischen Raum anstelle des Pixelraums durchführt
• Berühren Sie die Komponenten des ITK-Registrierungsframeworks und kennen Sie die verfügbaren Werte

#Werden die Zieldaten vor der Ausrichtung im gleichen geometrischen Raum kalibriert? Wenn Sie ein Resampling benötigen, tun Sie es.
#Ich frage mich, ob ich es mit einem IJ-Sample machen werde. zukünftige Arbeit
# https://imagej.nih.gov/ij/images/
# https://imagej.nih.gov/ij/images/pet-series/
# https://imagej.nih.gov/ij/images/ct-scan.zip

#Laden Sie die zu verwendenden Daten (Referenzbild)
PixelType = itk.ctype('float')
fixedImage = itk.imread('data/BrainProtonDensitySliceBorder20.png', PixelType)
plt.imshow(itk.array_view_from_image(fixedImage))

#Bild, das nicht in Position ist
movingImage = itk.imread('data/BrainProtonDensitySliceShifted13x17y.png', PixelType)
plt.imshow(itk.array_view_from_image(movingImage))

#Sammeln Sie die für die Konvertierung erforderlichen Informationen
#Anzahl der Dimensionen
Dimension = fixedImage.GetImageDimension()
#Bildtyp
FixedImageType = type(fixedImage)
MovingImageType = type(movingImage)
#Art der zu konvertierenden Daten
TransformType = itk.TranslationTransform[itk.D, Dimension]
#Initialisierung der Untertasse
initialTransform = TransformType.New()
#Definieren Sie eine Funktion, die optimiert wird, um Ausrichtungsfehler zu minimieren
optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
        LearningRate=4,
        MinimumStepLength=0.001,
        RelaxationFactor=0.5,
        NumberOfIterations=200)
#Ausrichtungsalgorithmus einstellen (mittlerer quadratischer Fehler)
metric = itk.MeanSquaresImageToImageMetricv4[
    FixedImageType, MovingImageType].New()
#Initialisieren Sie das Ausrichtungsobjekt
registration = itk.ImageRegistrationMethodv4.New(FixedImage=fixedImage,
        MovingImage=movingImage,
        Metric=metric,
        Optimizer=optimizer,
        InitialTransform=initialTransform)
#Ausrichtungseinstellungen
#Konvertierungstyp des bewegten Bildes
movingInitialTransform = TransformType.New()
initialParameters = movingInitialTransform.GetParameters()
initialParameters[0] = 0 #x
initialParameters[1] = 0 #y
movingInitialTransform.SetParameters(initialParameters)
registration.SetMovingInitialTransform(movingInitialTransform)
#Kriterienumrechnungstyp
identityTransform = TransformType.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)

registration.SetNumberOfLevels(1) # number of multi-resolution levels
registration.SetSmoothingSigmasPerLevel([0]) #Glättungsstufe
registration.SetShrinkFactorsPerLevel([1]) #Schrumpfstufe
#Aktualisieren Sie die Ausrichtungseinstellungen und bestätigen Sie
registration.Update()
#Holen Sie sich das Ergebnis des Konvertierungsprozesses
transform = registration.GetTransform()
finalParameters = transform.GetParameters()
#Bewegungsentfernung
translationAlongX = finalParameters.GetElement(0)
translationAlongY = finalParameters.GetElement(1)
#Anzahl der Wiederholungen
numberOfIterations = optimizer.GetCurrentIteration()
#Methodengenauigkeit (kleiner ist besser)
bestValue = optimizer.GetValue()

print("Result = ")
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Iterations    = " + str(numberOfIterations))
print(" Metric value  = " + str(bestValue))

#Bildkomposition
CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(movingInitialTransform)
outputCompositeTransform.AddTransform(registration.GetModifiableTransform())
#Resampling
resampler = itk.ResampleImageFilter.New(Input=movingImage,
        Transform=outputCompositeTransform,
        UseReferenceImage=True,
        ReferenceImage=fixedImage)
#Ergebnisanzeige
resampler.SetDefaultPixelValue(100)

OutputPixelType = itk.ctype('unsigned char')
OutputImageType = itk.Image[OutputPixelType, Dimension]

resampler.Update()

result = resampler.GetOutput()
plt.subplot(1,2,1)
plt.title("registered")
plt.imshow(itk.array_view_from_image(result))
plt.subplot(1,2,2)
plt.title("moving_org")
plt.imshow(itk.array_view_from_image(movingImage))

Lassen Sie uns die Fehlausrichtung durch den Unterschied sehen

difference = itk.SubtractImageFilter.New(Input1=fixedImage,
        Input2=resampler)
resampler.SetDefaultPixelValue(1)

difference.Update()
plt.title("difference")
plt.imshow(itk.array_view_from_image(difference.GetOutput()))

#Fehlausrichtung vor und nach der Registrierung
resampler.SetTransform(identityTransform)
difference.Update()
plt.imshow(itk.array_view_from_image(difference.GetOutput()))

Als nächstes werden verschiedene Arten von Bildern (multimodal) ausgerichtet. Mattes MutualInformation, die häufig für die Ausrichtung in multimodalen Fällen verwendet wird, wird verwendet. Verwenden Sie eine starre Transformation.

#Referenzbild
t1_image = itk.imread('data/DzZ_T1.nrrd', itk.F)
print("Pixel spacing:", t1_image.GetSpacing())
plt.imshow(itk.array_view_from_image(t1_image))

#Drehen des T2-verbesserten Bildes
t2_image = itk.imread('data/DzZ_T2.nrrd', itk.F)
print("Pixel spacing:", t2_image.GetSpacing())
plt.imshow(itk.array_view_from_image(t2_image))

#Ergänzung: Richtung erhalten
pyDirM=itk.GetArrayFromVnlMatrix(t2_image.GetDirection().GetVnlMatrix().as_matrix())
print("Direction matrix:\n", pyDirM)

Die Ausrichtung erfolgt anhand der Korrelation von Mattes Mutual Information, die häufig für die Ausrichtung in multimodalen Fällen verwendet wird.

T1ImageType = type(t1_image)
T2ImageType = type(t2_image)
assert T1ImageType==T2ImageType, "Images must have the same pixel type!"

#Initialisieren Sie den Konvertierungstyp als Starrkörperkonvertierung
TransformTypeR = itk.Rigid2DTransform[itk.D]
initialTransform = TransformTypeR.New()

MetricType = itk.MattesMutualInformationImageToImageMetricv4[
        T1ImageType, T2ImageType]
metric = MetricType.New()

scales = itk.OptimizerParameters[itk.D](initialTransform.GetNumberOfParameters())
translation_scale = 1.0 / 1000
scales.SetElement(0, 1.0)
scales.SetElement(1, translation_scale)
scales.SetElement(2, translation_scale)

#Es kann die Skala automatisch schätzen, scheint jedoch nicht mit Mattes Mutual Information zu funktionieren.
#ScalesEstimatorType = itk.RegistrationParameterScalesFromPhysicalShift[MetricType]
#scalesEstimator = ScalesEstimatorType.New(Metric=metric)

optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
        Scales=scales,
        #ScalesEstimator=scalesEstimator,
        #LearningRate=1.0,
        MinimumStepLength=0.001,
        # RelaxationFactor=0.8,
        NumberOfIterations=300)

registration = itk.ImageRegistrationMethodv4.New(FixedImage=t1_image,
        MovingImage=t2_image,
        Metric=metric,
        MetricSamplingPercentage=0.1, # 10%
        Optimizer=optimizer,
        InitialTransform=initialTransform)

initialParameters = initialTransform.GetParameters()
#Offset-Einstellung
initialParameters[0] = 0.0 # rotation in radians
initialParameters[1] = -40.0 # x translation in millimeters
initialParameters[2] = +30.0 # y translation in millimeters

initialTransform.SetParameters(initialParameters)

resampler = itk.ResampleImageFilter.New(Input=t2_image,
        Transform=initialTransform,
        UseReferenceImage=True,
        ReferenceImage=t1_image)
resampler.SetDefaultPixelValue(0)
resampler.Update()

plt.subplot(131)
plt.imshow(itk.GetArrayFromImage(resampler.GetOutput()))
plt.subplot(132)
plt.imshow(itk.GetArrayFromImage(t1_image))

Es ist ziemlich aus. Vielleicht sollte ich die Parameter anpassen. Haben Sie keine Angst, es mit dem ausgerichteten Bild noch einmal zu wiederholen.

registration.SetMovingInitialTransform(initialTransform)

Dimension = t1_image.GetImageDimension()
identityTransform = TransformTypeR.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)

registration.SetNumberOfLevels(1)
registration.SetSmoothingSigmasPerLevel([0])
registration.SetShrinkFactorsPerLevel([1])

try:
    registration.Update()
except RuntimeError as exc:
    print("Exception ocurred:\n", exc, "\n\n")
finally:
    print("Registration finished")
    
CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(registration.GetTransform())
outputCompositeTransform.AddTransform(initialTransform)
    
resamplerResult = itk.ResampleImageFilter.New(Input=t2_image,
        Transform=outputCompositeTransform,
        UseReferenceImage=True,
        ReferenceImage=t1_image)
resamplerResult.SetDefaultPixelValue(0)
resamplerResult.Update()
# final position of T2 image when resampled into T1 image's pixel grid
plt.imshow(itk.GetArrayFromImage(resamplerResult.GetOutput()))

Gibt Ausrichtungsinformationen aus.

transform = registration.GetTransform()
finalParameters = transform.GetParameters()
rotationInRadians = finalParameters.GetElement(0)
translationAlongX = finalParameters.GetElement(1)
translationAlongY = finalParameters.GetElement(2)
rotationInDegrees = rotationInRadians * 360 / 3.141592 # radian to degree
numberOfIterations = optimizer.GetCurrentIteration()

bestValue = optimizer.GetValue()

print("Result = ")
print(" Rotation degr = " + str(rotationInDegrees))
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Rotation rad. = " + str(rotationInRadians))
print(" Iterations    = " + str(numberOfIterations))
print(" Metric value  = " + str(bestValue))

Schauen wir uns den Unterschied an.

difference = itk.AddImageFilter.New(Input1=t1_image,
        Input2=resamplerResult)
difference.Update()
plt.imshow(itk.GetArrayFromImage(difference.GetOutput()))

das ist alles.

Reference

• https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization
• https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch4.html#x37-1770004