Analysis of X-ray microtomography image by Python


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. stk_0001_labbeled_rock.gif

Finally, we get such a 3D image. Volume_Viewer_1.png

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.

636863294260796061IE.png (Image quoted from

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.

1280px-Computed_tomography_of_human_brain_-_large-1.png (Image taken from

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

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.

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)

コメント 2020-02-27 225010.png

Hopefully the image can be read as below.

コメント 2020-02-27 225158.png

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)

An example of the loaded image is shown below. ダウンロード (9).png

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. ダウンロード (11)-1.png

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. ダウンロード (11)-2.png

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)
    __,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

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.


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

#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.xlabel('size of vacancy (px)',fontsize=12)

The result looks like this. The X axis is logarithmic. ダウンロード (13).png

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),
        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 = 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')

#Save image
for num,img in enumerate(void_colored):
    name = str(num).zfill(4)

The areas colored in this way are as follows. ダウンロード (14).png

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. Volume_Viewer_1.png

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


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!

Recommended Posts

Analysis of X-ray microtomography image by Python
Grayscale by matrix-Reinventor of Python image processing-
Example of 3D skeleton analysis by Python
Image processing by matrix Basics & Table of Contents-Reinventor of Python image processing-
[OpenCV / Python] I tried image analysis of cells with OpenCV
Static analysis of Python programs
Image processing by python (Pillow)
Image capture of firefox using python
Expansion by argument of python dictionary
Behavior of python3 by Sakura's server
Clash of Clans and image analysis (3)
Introduction to image analysis opencv python
Story of power approximation by Python
Straight line drawing by matrix-Inventor's original research of Python image processing-
Image processing with Python 100 knocks # 4 Binarization of Otsu (discriminant analysis method)
Explanation of production optimization model by Python
Basics of binarized image processing with Python
Python: Basics of image recognition using CNN
[Learning memo] Basics of class by python
Conditional branching of Python learned by chemoinformatics
Python: Application of image recognition using CNN
Pandas of the beginner, by the beginner, for the beginner [Python]
Drawing with Matrix-Reinventor of Python Image Processing-
[Image classification] Facial expression analysis of dogs
Matrix Convolution Filtering-Reinventor of Python Image Processing-
Time variation analysis of black holes using python
Introduction of Python
python image processing
Extract dominant color of image by k-means clustering
Color page judgment of scanned image with python
Data analysis python
[Various image analysis with plotly] Dynamic visualization with plotly [python, image]
[Python] Comparison of Principal Component Analysis Theory and Implementation by Python (PCA, Kernel PCA, 2DPCA)
A memorandum of extraction by python bs4 request
Affine transformation by matrix (scaling, rotation, shearing, movement) -Reinventor of Python image processing-
Basics of Python ①
Image processing? The story of starting Python for
Judge Yosakoi Naruko by image classification of Tensorflow.
Super (concise) summary of image classification by ArcFace
Thorough comparison of three Python morphological analysis libraries
Copy of python
Analysis of shared space usage by machine learning
Sentiment analysis of large-scale tweet data by NLTK
Image processing by Python 100 knock # 11 smoothing filter (average filter)
Static analysis of Python code with GitLab CI
A well-prepared record of data analysis in Python
Image of closure
[Python + OpenCV] Whiten the transparent part of the image
Introduction of Python
[Language processing 100 knocks 2020] Summary of answer examples by Python
Extract the table of image files with OneDrive & Python
Analysis of financial data by pandas and its visualization (2)
Explanation of the concept of regression analysis using python Part 2
[Python] [Word] [python-docx] Simple analysis of diff data using python
List of posts related to optimization by Python to docker
Summary of Python articles by pharmaceutical company researcher Yukiya
[Python] Specify the range from the image by dragging the mouse
Get the image of "Suzu Hirose" by Google image search.
Challenge principal component analysis of text data with Python
Story of image analysis of PDF file and data extraction
List of Python code used in big data analysis