Analyse d'image de microtomographie à rayons X par Python
introduction
Dans cet article, je vais vous expliquer comment analyser l'image en coupe obtenue par microtomographie à rayons X en utilisant Python. (Pas une reconstruction à partir d'une image transmise)
Plus précisément, cibler l'image en coupe μ-CT de la roche carbonatée comme indiqué ci-dessous
Effectuez une segmentation en tenant compte de la connectivité tridimensionnelle comme indiqué ci-dessous.
Enfin, nous obtenons une telle image en trois dimensions.
Qu'est-ce que la microtomographie à rayons X?
Probablement 99% des personnes qui regardent Qiita ne connaissent pas cette méthode, je vais donc l'expliquer brièvement. Vous connaissez tous les appareils suivants.
(L'image est tirée de https://jp.medical.canon/general/CT_Comparison)
Cet appareil est communément appelé CT. CT est une abréviation pour Computed Tomography, qui est un appareil qui irradie le corps humain avec des rayons X pour obtenir une image en coupe. C'est très pratique car vous pouvez obtenir une image qui ressemble à une tranche du corps humain et vous pouvez voir la structure interne qui ne peut pas être vue de l'extérieur.
(Image tirée de https://en.wikipedia.org/wiki/CT_scan)
Bien sûr, il existe une demande autre que le corps humain pour ce «désir de voir à l'intérieur». Par exemple, on ne sait pas de l'extérieur combien de cavités (nids) se trouvent à l'intérieur de la pièce moulée, mais les nids affectent grandement les caractéristiques des pièces coulées, cette méthode d'inspection est donc utile.
Dans les applications industrielles, contrairement aux applications médicales, les dommages causés par les rayonnements ne doivent pas être pris en compte, et des rayons X puissants peuvent être convergés et utilisés, ce qui permet d'obtenir des images en coupe à haute résolution. Une telle technique est appelée micro-tomographie aux rayons X.
La tomographie aux rayons X est une technique très pratique. Cependant, vous obtenez essentiellement une série d'images en coupe, comme indiqué ci-dessus. À partir de là, par exemple, si vous souhaitez analyser la structure tridimensionnelle d'une cavité, vous devez utiliser un logiciel spécial ou écrire votre propre programme. Cette fois, je vais le faire en Python.
Télécharger les données d'analyse
Cette fois, nous utiliserons l'ensemble de données d'images de tomographie aux rayons X de roches carbonatées publiées dans les articles suivants.
Étant donné que la quantité de données de tomographie à rayons X disponibles gratuitement est étonnamment petite, je suis très reconnaissante qu'elles soient publiées dans une revue en libre accès comme celle-ci. Les données d'image peuvent être téléchargées à partir du site suivant.
Cette fois, j'ai utilisé "Rock-Reconstructed.raw" publié sur le site ci-dessus pour l'analyse. Il est d'environ 845 Mo (plus léger que les données d'image CT).
Je pense que l'extension ".raw" est une extension inconnue, mais si vous utilisez ImageJ, vous pouvez l'importer en tant que séquence d'images. Sélectionnez un fichier avec File => Import => Raw et définissez comme suit. (Les détails sont expliqués dans l'article, veuillez donc vérifier si vous souhaitez utiliser une autre image)
Espérons que l'image peut être lue comme ci-dessous.
Après cela, enregistrez-le dans un dossier approprié avec Fichier => Enregistrer sous => Séquence d'images.
Analyse d'image par OpenCV
Les données obtenues par microtomographie à rayons X peuvent désormais être enregistrées sous la forme d'une série d'images. Ensuite, nous effectuerons la suppression du bruit et la binarisation avec OpenCV.
Cette fois, je n'analyserai que le centre de l'image.
import numpy as np
import matplotlib.pyplot as plt
import cv2
import glob
#Charger le chemin de l'image
path = glob.glob('./Rock-reconstructed/*.tif')
#Chargement des images
images = []
for n,i in enumerate(path):
img = cv2.imread(i,0)
images.append(img[100:400,100:400])
Un exemple de l'image chargée est présenté ci-dessous.
Pour être honnête, c'est une très belle image, donc je n'en ai pas besoin, mais je vais essayer de supprimer le bruit. Utiliser OpenCV Non-local signifie le débruitage.
# Non-Élimination du bruit par des moyens locaux
test = images[200]
dn_test = cv2.fastNlMeansDenoising(test,None,10,7,21)
Voici une comparaison avec et sans suppression du bruit.
La binarisation à l'aide d'OpenCV se fait comme suit. Utilisez la binarisation automatique d'Otsu.
#Binariser l'image originale par la méthode d'Otsu
th,th_test = cv2.threshold(test,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
La comparaison avec et sans suppression du bruit ressemble à ceci.
Depuis que les petits trous ont disparu, je me demande si la suppression du bruit est superflue ... Allons-y sans cette fois.
Appliquez ce processus à une série d'images à la fois. Ensuite, extrayez la zone avec cv2.findContours, puis excluez les cavités inférieures à 10px.
th_denoised = [] #Mettre les images dans cette liste
for i in images:
#Suppression du bruit(Aucun cette fois)
#dn = cv2.fastNlMeansDenoising(i,None,10,7,21)
#Binarisation
__,th_dn = cv2.threshold(i,0,255,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)
#Détection de zone de cavité
i__,cnts,__ = cv2.findContours(th_dn,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
#Extraire uniquement la zone de la cavité de 25 pixels ou plus
img = np.zeros_like(i).astype(np.uint8)
for cnt in cnts:
if(cv2.contourArea(cnt) > 10):
img = cv2.drawContours(img, [cnt], 0, 255, -1)
#Ajouter à la liste
th_denoised.append(img)
Étiqueter la zone concaténée avec cc3d
Vous avez maintenant une zone creuse dans la section transversale. Ensuite, nous étiquetons les cavités continues en tenant compte de la structure de connexion tridimensionnelle. Une bibliothèque appelée connected-components-3d est utilisée pour cela.
Voir la page ci-dessus pour l'algorithme détaillé. Vous pouvez l'installer à l'aide de pip comme suit.
pip install connected-components-3d
L'étiquetage est très simple. Et c'est rapide.
import cc3d
connectivity = 6
th_denoised = np.array(th_denoised)
labels_out = cc3d.connected_components(th_denoised,connectivity=connectivity)
La connectivité indique le nombre de boxels de proximité nécessaires pour être considérés comme la même zone. Par exemple, 6 fait référence à ± 1 cases adjacentes dans chacune des directions x, y et z. Outre 6, il semble que 18, 26 peuvent être sélectionnés (voir document).
Pour chaque pixel de l'image, le numéro correspondant à l'étiquette de zone est attribué à partir de 1. (0 est l'arrière-plan) Par conséquent, par exemple, le nombre de zones peut être obtenu comme suit.
#Sortie du nombre de zones
numbers = np.max(labels_out) - 1
print(numbers)
>> 10001
Regardons maintenant la distribution des tailles de cavité. Il peut être analysé comme suit.
#Calcul du volume de la cavité
areas = []
for num in range(np.max(labels_out)):
count = np.count_nonzero(labels_out == num)
areas.append(count)
#Affichage de la distribution
fig = plt.figure(figsize=(7,4),dpi=100,facecolor='white')
plt.hist(areas[1:], bins=np.logspace(0, 6, 20),rwidth=0.8)
plt.xscale('log')
plt.xlabel('size of vacancy (px)',fontsize=12)
plt.ylabel('Amounts',fontsize=12)
plt.show()
Le résultat ressemble à ceci. L'axe X est logarithmique.
Apparemment, la distribution du volume de la cavité (nombre de pixels) peut être représentée par une distribution normale logarithmique. Pour le représenter dans une taille concrète, vous avez besoin de l'échelle de l'image et des valeurs de largeur de tranche dans la direction z.
Colorez la zone de la cavité avec RVB
De cette manière, cc3d facilite la séparation et l'analyse des régions de cavité. Cependant, lorsque vous le visualisez sous forme d'image, il est plus facile de le voir en codant par couleur chaque zone avec RVB.
Voici un programme qui fait cela:
import random
##Sélectionnez au hasard la couleur RVB correspondant à chaque étiquette
param_list = np.linspace(0,np.max(labels_out),np.max(labels_out)+1,dtype=np.uint16)
color_list = {}
for i in param_list:
if(i != 0):
color_list[str(i)] = [random.randint(0,255),
random.randint(0,255),
random.randint(0,255)]
else:
color_list[str(i)] = [0,0,0]
##Colorez chaque zone de cavité
void_colored = []
for img in labels_out:
h,w = img.shape
img_flat = img.flatten()
img_colored = np.zeros((img_flat.shape[0],3),dtype=np.uint8)
non_zero_idx = np.where(img_flat != 0)
for n in non_zero_idx[0]:
img_colored[n] = color_list[str(img_flat[n])]
void_colored.append(img_colored)
void_colored = np.array(void_colored)
void_colored = void_colored.reshape(void_colored.shape[0],h,w,3)
##Assurez-vous que les cavités sont correctement colorées
fig = plt.figure(figsize=(5,5),dpi=100,facecolor='white')
plt.imshow(void_colored[0])
plt.show()
#Enregistrer l'image
for num,img in enumerate(void_colored):
name = str(num).zfill(4)
cv2.imwrite('./labbeled_rock/'+str(name)+'.png',img)
Les zones colorées de cette manière sont les suivantes.
Affichage 3D par ImageJ
J'aimerais utiliser Python pour l'affichage 3D, mais je ne trouve pas de bonne méthode, j'utiliserai donc ImageJ.
Importez l'image de la zone de coloration enregistrée précédemment avec ImageJ => Fichier => Importer => Séquence d'images. Vous pouvez afficher la 3D avec Plugins => Volume Viewer dans l'état chargé.
Comme je ne connais pas la largeur de la tranche dans la direction Z, le résultat de l'affichage 3D après avoir ajusté l'échelle de manière appropriée est le suivant.
À première vue, vous pouvez voir que les zones de cavité bien connectées sont étiquetées avec la même couleur. Ça marche!
finalement
Il y avait une personne près de moi qui faisait des analyses de tomographie par rayons X, et j'ai aidé un peu, alors j'en ai rédigé un mémorandum. Cet article est sérieusement comme un chercheur en matériaux.
C'est tellement une niche qu'elle ne semble pas coller au lectorat de Qiita ... Nous sommes impatients de vous aider!
Recommended Posts