Utilisation élémentaire d'ITK apprise avec Python
Notification
J'ai écrit cet article comme mémo de mon étude. La source de référence est ce cahier. https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization
ITK Une bibliothèque indispensable pour le traitement d'images dans le domaine biomédical, spécialisée dans la technologie de visualisation d'images. Il est utilisé dans Slicer et ImageJ.
environnement
Google Colab (Lancez n'importe quel notebook.)
Exemple de données
Commencez par télécharger les exemples de données depuis GitHub à l'aide de subversion.
!sudo apt install subversion
#Obtenez les données
#Inclus dans l'URL Github/tree/Maître/Remplacer par le coffre
!svn checkout https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization/trunk/data
Installation ITK
#Installez ITK
!pip install itk
Image Filtering
Quoi apprendre ici
- Apprenez à connaître le modèle de pipeline de traitement d'image utilisé par ITK. --Image et Mesh sont utilisés pour représenter des objets de données.
- Un nouvel objet de données est créé au cours du processus de manipulation de l'objet de données.
- Les processus sont classés en tant qu'objets de processus et ont des objets source, filtre et mappeur. --source (lecteur, etc.) récupère les données et filtre les objets pour générer de nouvelles données. --mapper accepte les données pour la sortie vers des fichiers ou d'autres systèmes.
Le terme filtre est largement utilisé dans trois types:
- Pipeline de données Les données et les objets de processus sont généralement convertis les uns aux autres à l'aide des méthodes SetInput () et GetOutput (). Aucune nouvelle sortie ou nouvelle donnée de pixel pour aucun traitement ne se produira jusqu'à ce que la méthode Update () soit appelée à la fin du pipeline (processus ou objet de données).
- Mises à jour des pipelines Les données d'image multidimensionnelles sont volumineuses et peuvent être encore plus volumineuses. ITK résout ce problème grâce à ses capacités de streaming de données. --Diffusion Le streaming se fait en divisant l'image en zones non superposées à la fin du pipeline. Puis propagez le pipeline en tant que région demandée.
Il existe trois types de régions d'image dans ITK.
--BufferedRegion: zone de pixels stockée en mémoire --LargestPossibleRegion: Zone maximale possible de l'image --RequestedRegion: Région demandée par un seul chemin de traitement pendant le streaming (BufferedRegion et LargestPossibleRegion sont identiques ou plus grands que RequestedRegion)
Essayer.
import numpy as np
import itk
import matplotlib.pyplot as plt
#Pratiquez le filtrage spatial avec Pacman
fileName = 'data/PacMan.png'
reader = itk.ImageFileReader.New(FileName=fileName)
#Essayez d'utiliser un filtre gaussien
smoother = itk.RecursiveGaussianImageFilter.New(Input=reader.GetOutput())
print("reader's Output: %s" % reader.GetOutput())
print("smoother's Output: %s" % smoother.GetOutput())
#Mettre à jour pour générer des résultats de traitement()N'oublie pas.
smoother.Update()
print("reader's Output: %s" % reader.GetOutput())
print("smoother's Output: %s" % smoother.GetOutput())
image = reader.GetOutput()#Confirmation de l'original
plt.subplot(1,2,1),plt.title("original"),plt.imshow(itk.GetArrayFromImage(image))
smoothed = smoother.GetOutput()#Après le lissage
plt.subplot(1,2,2),plt.title("smoothing"),plt.imshow(itk.GetArrayFromImage(smoothed))
Si vous avez mis à jour avec un objet pipeline (bien que vous l'ayez déjà fait), le résultat du traitement sera le dernier
smoother.Update()
#Cependant, si vous modifiez le processus, vous devez à nouveau mettre à jour.
smoother.SetSigma(10.0)
smoother.Update()
plt.imshow(itk.GetArrayFromImage(smoothed))
De plus, je ne pense pas qu'il y en ait beaucoup, mais si vous apportez des modifications au lecteur au milieu, n'oubliez pas de mettre à jour le lecteur.
# reader.Modified()
#Vous pouvez essayer différents filtres.
# 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())#Par exemple, médiane
smoother.SetRadius(5)
smoother.Update()
plt.imshow(itk.GetArrayFromImage(smoother.GetOutput()))
Supplément
#Placez StreamingImageFilter à la fin du pipeline pour diffuser le pipeline.
#Smoother génère une sortie plusieurs fois pour chaque division de diffusion en continu de la zone d'image.
#Le lecteur ne peut pas être diffusé, donc la sortie n'est générée qu'une seule fois.
streamer = itk.StreamingImageFilter.New(Input=smoother.GetOutput())
streamer.SetNumberOfStreamDivisions(4)
reader.Modified()
streamer.Update()
#Il vous suffit de spécifier l'écriture avec l'extension
writer = itk.ImageFileWriter.New(Input=smoother.GetOutput())
writer.SetFileName('my_output.mha')
writer.SetNumberOfStreamDivisions(4)
writer.Update()
Image Segmentation
Quoi apprendre ici
- Connaître les méthodes de segmentation disponibles dans ITK (First, ConfidenceConnectedImageFilter)
- Comprendre comment la région rayonnante et ses paramètres changent de comportement
- Comprendre comment les level sets et leurs paramètres modifient leur comportement
Types de segmentation disponibles dans ITK
- Level Set
- Bassin versant (chute de gradient des caractéristiques de l'image et analyse le long des limites de la zone) --Composant connecté (une technique abstraite pour travailler avec des objets connectés utilisée dans la binarisation et l'étiquetage)
- Vote d'étiquette (pour compter les valeurs de pixel comme étiquettes et les utiliser dans le processus)
- Region Growing --Classifier (algorithme d'apprentissage automatique)
Ici, nous allons présenter un exemple de traitement pour Region Growing et Level Set.
Region growing
L'approche de base de l'algorithme d'expansion de région consiste à démarrer l'expansion à partir de la région de départ (généralement un ou plusieurs pixels) qui est considérée comme étant à l'intérieur de l'objet segmenté. Les pixels adjacents à cette zone d'objet segmentée sont évalués pour déterminer s'ils doivent être considérés comme faisant partie de l'objet. S'ils sont segmentés, ils sont ajoutés à la zone et le processus se poursuit tant que de nouveaux pixels sont ajoutés à la zone. L'algorithme d'expansion de région est le critère utilisé pour déterminer si une région contient des pixels, le type de connexion utilisé pour déterminer les voisins et le processus (ou processus) utilisé pour rechercher les pixels voisins. Stratégie) dépend de.
Ici, la méthode suivante est utilisée comme exemple.
- ConfidenceConnected
- ConnectedThreshold
- IsolatedConnected
Tout d'abord, vérifiez comment utiliser "Confidence Connected". ConfidenceConnected est une méthode permettant de déterminer les limites supérieure et inférieure de la limite de la zone, de trouver l'écart type par rapport aux valeurs entre elles et de définir la plage autorisée de cet écart type (comme l'effet de levier. SD * Multiplicateur).
# BrainProtonDensitySlice.Utiliser png
# 'unsigned char'La raison du chargement avec est d'afficher la méthode de conversion de données plus tard.
input_image = itk.imread('data/BrainProtonDensitySlice.png', itk.ctype('unsigned char'))
# View the input image
plt.imshow(itk.GetArrayFromImage(input_image))
Vérifiez comment utiliser la classe "ConfidenceConnectedImageFilter". Si vous ne savez pas comment l'utiliser, recherchez-le.
confidence_connected = itk.ConfidenceConnectedImageFilter.New(input_image)
#Vérifiez comment utiliser en utilisant la fonction d'aide de python
help(confidence_connected)
#Définir arbitrairement le paramètre de filtre connecté à la confiance
confidence_connected.SetMultiplier(2.3)#Tolérance, ± SD dans cet exemple*2.3
confidence_connected.SetNumberOfIterations(5)
confidence_connected.SetInitialNeighborhoodRadius(3)
confidence_connected.SetReplaceValue(255)
Si vous ne savez pas comment utiliser une fonction, recherchez-la.
# *Multiplier*?
confidence_connected.SetMultiplier?
Définissez le point de départ "point de départ" de l'expansion de la zone. Notez que l'origine de l'objet ITK Image se trouve en bas à gauche des coordonnées de l'image. (NumPy est en haut à gauche)
# Set the seed points
confidence_connected.SetSeed([100, 110])
#Exécuter et visualiser
confidence_connected.Update()
plt.imshow(itk.GetArrayFromImage(confidence_connected))
#Essayez d'améliorer le résultat du traitement. Appliquez un filtre médian à l'image d'origine.
median_filtered_image = itk.median_image_filter(input_image, radius=1)
#Mettez à nouveau les paramètres à jour. Entrez l'image après le filtre médian.
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))
Ensuite, comparons les trois processus lumineux de la région. 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) Ensuite, essayez la segmentation Fast Marching de la méthode LevelSet. Si la zone cible de segmentation est une structure simple, un algorithme d'expansion rapide appelé marche rapide peut être utilisé. Simplifiez l'image avec un filtre spatial avant d'essayer.
#Découvrez les types de données d'image que le filtre de lissage peut gérer
itk.CurvatureAnisotropicDiffusionImageFilter.GetTypes()
#Comment convertir des données d'image dans un type de données pris en charge
Dimension = input_image.GetImageDimension()
FloatImageType = itk.Image[itk.ctype('float'), Dimension]
caster = itk.CastImageFilter[input_image, FloatImageType].New()
caster.SetInput(input_image)
#Appliquez un filtre lissant.
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))
#Exécuter le filtre de détection de bord tel quel
gradient = itk.gradient_magnitude_recursive_gaussian_image_filter(smoothed_image, sigma=1.0)
plt.imshow(itk.GetArrayFromImage(gradient))
#Continuer à exécuter le filtre sigmoïde
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))
#Marche à grande vitesse
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))
Binarisez l'image.
time_threshold = 100 #Définir le seuil
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))
À partir de là, essayez le jeu de niveaux de détection de forme. La vitesse de traitement est améliorée, mais la courbure n'est pas bonne.
Le pipeline ressemble à ceci:
Commencez par créer une image équilibrée dans la partie inférieure de la figure.
#L'image binaire est la matière blanche.Utiliser png
binary_mask = itk.imread('data/WhiteMatter.png', itk.ctype('float'))
smoothed_image = itk.smoothing_recursive_gaussian_image_filter(binary_mask)
#Traitement des gradations
rescaler = itk.IntensityWindowingImageFilter[FloatImageType,FloatImageType].New(smoothed_image)
rescaler.SetWindowMinimum(0)
rescaler.SetWindowMaximum(255)
#Inverser le niveau d'entrée pour avoir une valeur négative en interne
rescaler.SetOutputMinimum(5.0)
rescaler.SetOutputMaximum(-5.0)
rescaler.Update()
plt.imshow(itk.GetArrayFromImage(rescaler))
Ensuite, créez l'image de fonction dans la figure.
#Pré-traiter en utilisant l'image chargée dans la procédure ci-dessus
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
#Image caractéristique dans la figure
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_Exécuter le traitement de l'ensemble
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)
#Binarisation
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))
Introduisons quelques paramètres liés aux conditions de courbure et de propagation du jeu de niveaux dans le jeu de niveaux de détection de forme. 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
Quoi apprendre ici
- Savoir comment ITK effectue des calculs dans l'espace physique
- Comprendre pourquoi ITK effectue l'enregistrement dans l'espace physique au lieu de l'espace pixel --Touchez les composants du cadre d'enregistrement ITK et connaissez les valeurs disponibles
#Les données cibles sont-elles calibrées dans le même espace géométrique avant l'alignement? Si vous avez besoin d'un rééchantillonnage, faites-le.
#Je me demande si je vais le faire avec un échantillon IJ. travail futur
# https://imagej.nih.gov/ij/images/
# https://imagej.nih.gov/ij/images/pet-series/
# https://imagej.nih.gov/ij/images/ct-scan.zip
#Charger les données à utiliser (image de référence)
PixelType = itk.ctype('float')
fixedImage = itk.imread('data/BrainProtonDensitySliceBorder20.png', PixelType)
plt.imshow(itk.array_view_from_image(fixedImage))
#Image déplacée
movingImage = itk.imread('data/BrainProtonDensitySliceShifted13x17y.png', PixelType)
plt.imshow(itk.array_view_from_image(movingImage))
#Rassemblez les informations nécessaires à la conversion
#Nombre de dimensions
Dimension = fixedImage.GetImageDimension()
#Type d'image
FixedImageType = type(fixedImage)
MovingImageType = type(movingImage)
#Type de données à convertir
TransformType = itk.TranslationTransform[itk.D, Dimension]
#Initialisation de la soucoupe
initialTransform = TransformType.New()
#Définissez une fonction qui optimise pour minimiser les erreurs d'alignement
optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
LearningRate=4,
MinimumStepLength=0.001,
RelaxationFactor=0.5,
NumberOfIterations=200)
#Définir l'algorithme d'alignement (erreur quadratique moyenne)
metric = itk.MeanSquaresImageToImageMetricv4[
FixedImageType, MovingImageType].New()
#Initialiser l'objet d'alignement
registration = itk.ImageRegistrationMethodv4.New(FixedImage=fixedImage,
MovingImage=movingImage,
Metric=metric,
Optimizer=optimizer,
InitialTransform=initialTransform)
#Paramètres d'alignement
#Type de conversion d'image en mouvement
movingInitialTransform = TransformType.New()
initialParameters = movingInitialTransform.GetParameters()
initialParameters[0] = 0 #x
initialParameters[1] = 0 #y
movingInitialTransform.SetParameters(initialParameters)
registration.SetMovingInitialTransform(movingInitialTransform)
#Type de conversion des critères
identityTransform = TransformType.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)
registration.SetNumberOfLevels(1) # number of multi-resolution levels
registration.SetSmoothingSigmasPerLevel([0]) #Niveau de lissage
registration.SetShrinkFactorsPerLevel([1]) #Niveau de rétrécissement
#Mettez à jour les paramètres d'alignement et confirmez
registration.Update()
#Obtenez le résultat du processus de conversion
transform = registration.GetTransform()
finalParameters = transform.GetParameters()
#Distance de déplacement
translationAlongX = finalParameters.GetElement(0)
translationAlongY = finalParameters.GetElement(1)
#Nombre de répétitions
numberOfIterations = optimizer.GetCurrentIteration()
#Précision de la méthode (plus petite est meilleure)
bestValue = optimizer.GetValue()
print("Result = ")
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Iterations = " + str(numberOfIterations))
print(" Metric value = " + str(bestValue))
#Composition d'image
CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(movingInitialTransform)
outputCompositeTransform.AddTransform(registration.GetModifiableTransform())
#Rééchantillonnage
resampler = itk.ResampleImageFilter.New(Input=movingImage,
Transform=outputCompositeTransform,
UseReferenceImage=True,
ReferenceImage=fixedImage)
#Affichage des résultats
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))
Voyons le désalignement par la différence
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()))
#Mauvais alignement avant et après l'enregistrement
resampler.SetTransform(identityTransform)
difference.Update()
plt.imshow(itk.array_view_from_image(difference.GetOutput()))
Ensuite, différents types d'images (multimodales) sont alignés. Mattes MutualInformation, qui est souvent utilisé pour l'alignement dans les cas multimodaux, est utilisé. Utilisez une transformation rigide.
#Image de référence
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))
#Rotation de l'image améliorée T2
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))
#Supplément: obtenir la direction
pyDirM=itk.GetArrayFromVnlMatrix(t2_image.GetDirection().GetVnlMatrix().as_matrix())
print("Direction matrix:\n", pyDirM)
L'alignement est effectué en utilisant la corrélation par Mattes Mutual Information, qui est souvent utilisée pour l'alignement dans les cas multimodaux.
T1ImageType = type(t1_image)
T2ImageType = type(t2_image)
assert T1ImageType==T2ImageType, "Images must have the same pixel type!"
#Initialiser le type de conversion en tant que conversion de corps rigide
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)
#Il peut estimer automatiquement l'échelle, mais il ne semble pas fonctionner avec Mattes Mutual Information.
#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()
#Réglage du décalage
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))
C'est assez bizarre. Peut-être que je devrais ajuster les paramètres. N'ayez pas peur de le répéter une fois de plus en utilisant l'image alignée.
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()))
Produit des informations d'alignement.
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))
Regardons la différence.
difference = itk.AddImageFilter.New(Input1=t1_image,
Input2=resamplerResult)
difference.Update()
plt.imshow(itk.GetArrayFromImage(difference.GetOutput()))
c'est tout.
Reference
- https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization
- https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch4.html#x37-1770004
Recommended Posts