Analysis of X-ray microtomography image by Python

Introduction

In this article, I will explain how to analyze a cross-sectional image obtained by X-ray microtomography using Python. (Not a reconstruction from a transparent image)

Specifically, targeting the μ-CT cross-sectional image of carbonate rock as shown below,

Perform segmentation considering three-dimensional connectivity as shown below.

Finally, we get such a 3D image.

What is X-ray microtomography?

Probably 99% of people who are watching Qiita do not know this method, so I will explain it briefly. You all know the following devices.

(Image quoted from https://jp.medical.canon/general/CT_Comparison)

This device is commonly referred to as CT. CT is an abbreviation for Computed Tomography, which is a device that irradiates the human body with X-rays to obtain a cross-sectional image. It is very convenient because you can get an image that looks like a slice of the human body and you can see the internal structure that cannot be seen from the outside.

(Image taken from https://en.wikipedia.org/wiki/CT_scan)

Naturally, there is a demand for this "want to see the inside" other than the human body. For example, it is not known from the outside how many cavities (nests) are inside the casting, but the cavities greatly affect the characteristics of the cast part, so this method of inspection is useful.

In industrial applications, unlike medical applications, radiation damage does not have to be considered, and powerful X-rays can be converged and used, so high-resolution cross-sectional images can be obtained. Such a technique is called X-ray microtomography.

X-ray tomography is a very convenient technique. However, what you basically get is a series of cross-sectional images as shown above. From here, for example, if you want to analyze the 3D structure of a cavity, you need to use special software or write your own program. This time I'll do this in Python.

Download analysis data

This time, we will use the dataset of X-ray tomography images of carbonate rocks published in the following papers.

Pak, T., Archilha, N., Mantovani, I. et al. An X-ray computed micro-tomography dataset for oil removal from carbonate porous media. Sci Data 6, 190004 (2019). https://doi.org/10.1038/sdata.2019.4

Since the amount of free-to-use X-ray tomography data is surprisingly small, I am very grateful that it is published in an open access journal. Image data can be downloaded from the following site.

Pak, Tannaz; Archilha, Nathaly Lopes; Mantovani, Iara Frangiotti; Moreira, Anderson Camargo; Butler, Ian B. (2019): X-ray computed micro-tomography dataset for nanoparticle-based oil removal from porous media. figshare. Collection. https://doi.org/10.6084/m9.figshare.c.4190738.v1

This time, I used "Rock-Reconstructed.raw" published on the above site for analysis. It is about 845MB (lighter as CT image data).

I think the extension ".raw" is an unfamiliar extension, but if you use ImageJ, you can import it as an Image Sequence. Select the file with File => Import => Raw and set as follows. (Details are explained in the paper, so please check if you want to use another image)

Hopefully the image can be read as below.

After that, save it in an appropriate folder with File => Save As => Image Sequence.

Image analysis with OpenCV

The data obtained by X-ray microtomography can now be saved as a series of images. Next, we will perform noise removal and binarization with OpenCV.

This time, we will analyze only the center of the image.

import numpy as np
import matplotlib.pyplot as plt
import cv2
import glob

#Load image path
path = glob.glob('./Rock-reconstructed/*.tif')

#Loading images
images = []
for n,i in enumerate(path):
    img = cv2.imread(i,0)
    images.append(img[100:400,100:400])

An example of the loaded image is shown below.

To be honest, it's a pretty beautiful image, so I don't need it, but I'll try to remove noise. Use OpenCV Non-local means denoising.

# Non-Noise removal using local means
test = images[200]
dn_test = cv2.fastNlMeansDenoising(test,None,10,7,21)

A comparison with and without denoising is shown below.

Binarization using OpenCV is done as follows. Use Otsu's automatic binarization.

#Binarize the original image by Otsu's method
th,th_test = cv2.threshold(test,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)

The comparison with and without noise removal looks like this.

Since the small holes have disappeared, noise removal is unnecessary ... Let's go without this time.

Apply this process to a series of images at once. Then extract the area with cv2.findContours and then exclude the cavities below 10px.

th_denoised = [] #Put images in this list

for i in images:

    #Noise removal(None this time)
    #dn = cv2.fastNlMeansDenoising(i,None,10,7,21)
    #Binarization
    __,th_dn = cv2.threshold(i,0,255,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)
    #Cavity area detection
    i__,cnts,__ = cv2.findContours(th_dn,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)

    #Extract only the cavity area of 25px or more
    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)

    #Add to list
    th_denoised.append(img)

Labeling of concatenated areas with cc3d

Now you have a hollow area in the cross section. Next, we will label the continuous cavities in consideration of the three-dimensional connection structure. A library called connected-components-3d is used for this.

connected-components-3d

See the page above for the detailed algorithm. You can install it using pip as follows.

pip install connected-components-3d

Labeling is very easy. And it's fast.

import cc3d

connectivity = 6
th_denoised = np.array(th_denoised)
labels_out = cc3d.connected_components(th_denoised,connectivity=connectivity)

connectivity indicates the number of adjacent voxels needed to be considered the same area. For example, 6 refers to ± 1 adjacent voxels in each of x, y, and z directions. Besides 6, it seems that you can choose 18, 26 (see documentation).

For the pixels of each image, the number corresponding to the area label is assigned from 1. (0 is the background) Therefore, for example, the number of areas can be obtained as follows.

#Output the number of areas
numbers = np.max(labels_out) - 1
print(numbers)
>> 10001

Now let's look at the distribution of cavity sizes. It can be analyzed as follows.

#Calculation of cavity volume
areas = []
for num in range(np.max(labels_out)):
    count = np.count_nonzero(labels_out == num)
    areas.append(count)

#Display of 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()

The result looks like this. The X axis is logarithmic.

Apparently, the distribution of cavity volume (number of pixels) can be represented by a lognormal distribution. To represent it in a concrete size, you need the image scale and the slice width values in the z direction.

Color the cavity area with RGB

In this way, cc3d makes it easy to separate and analyze cavity regions. However, when viewing it as an image, it is easier to see it by color-coding each area with RGB.

Here is a program that does this:

import random

##Randomly select the RGB color corresponding to each label
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]

##Color each cavity area
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)

##Make sure the cavities are colored correctly
fig = plt.figure(figsize=(5,5),dpi=100,facecolor='white')
plt.imshow(void_colored[0])
plt.show()

#Save image
for num,img in enumerate(void_colored):
    name = str(num).zfill(4)
    cv2.imwrite('./labbeled_rock/'+str(name)+'.png',img)

The areas colored in this way are as follows.

3D display by ImageJ

I'd like to use Python for 3D display, but I can't find a good method, so I'll use ImageJ.

Import the image of the coloring area saved earlier with ImageJ => File => Import => Image Sequence. You can display in 3D with Plugins => Volume Viewer in the loaded state.

Since I do not know the slice width in the Z direction, the result of 3D display after adjusting the scale appropriately is as follows.

At first glance, you can see that the well-connected cavity areas are labeled with the same color. It's working!

Finally

There was a person close to me who was doing X-ray tomography analysis, and I helped a little, so I wrote a memorandum. This article is seriously like a material researcher.

It's so niche that it's unlikely to pierce Qiita's readership ... We look forward to helping you!