Opening and Preparing Images

The examples in the generators section have numerous examples on how to generate artificial images, but ultimately most people will be working with real images obtained from X-ray tomography or serial sectioning. This tutorial will cover a few different aspects about how to open these images for using a Jupyter notebook, or a Python workflow in general.

Let’s import the packages and functions that we’ll be needing:

import imageio.v2 as imageio
import porespy as ps
import numpy as np
from skimage.io import imread_collection
from matplotlib.pyplot import subplots

[12:36:53] ERROR    PARDISO solver not installed, run `pip install pypardiso`. Otherwise,          _workspace.py:56
                    simulations will be slow. Apple M chips not supported.                                         

Opening a Folder of 2D tif Images

Using ImageJ

ImageJ is the ultimate image processing tool kit. It has many functions and filters for doing quantitative analysis, usually added via plug-ins (see note below), but it’s also very good at converting between file types and formats.

If you have a folder full of 2D images which are either microscope images obtained by serial sectioning or reconstructions of backprojections from X-ray tomography, then you’ll want to convert these into a single tif image. tif supports 3D images, so all the individual 2D images can be stacked together into a single file.

This can be done in ImageJ using “File/Import/Image Sequence” as shown in the screenshots below. Note that the images must be named with a numerical value indicating its order in the sequence, like “slice_001.tif”.

A dialog box will appear asking you about some of the details, which are usually correct by default:

Finally, you’ll want to save this new 3d tif image, which can be done with “save as”:

ImageJ and Fiji

ImageJ has a plugin architecture to add functionality. At some point it was realized that most users install numerous plugins, so Fiji was created, which is ImageJ with a large selection of the most popular plugins preinstalled. The name Fiji is actually a recursive acronym for “Fiji Is Just ImageJ”

Directly in Python

If you’d rather not get manually involved, then it’s also possible to do this in Python. The scikit-image package comes with a handy funciton for just this purpose. It only takes a few minutes to get used to how it works:

seq = imread_collection("../../docs/_static/images/Sample 1/*.tif")

Now let’s inspect the result. We can see a list of files loaded in the files attribute:

seq.files

[]

And we can access the loaded images from the seq object as if it were a list:

fig, ax = subplots(figsize=[5, 5])
ax.imshow(seq[0]);

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[4], line 2
      1 fig, ax = subplots(figsize=[5, 5])
----> 2 ax.imshow(seq[0]);

File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/skimage/io/collection.py:301, in ImageCollection.__getitem__(self, n)
    298     raise TypeError('slicing must be with an int or slice object')
    300 if type(n) is int:
--> 301     n = self._check_imgnum(n)
    302     idx = n % len(self.data)
    304     if ((self.conserve_memory and n != self._cached) or
    305             (self.data[idx] is None)):

File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/skimage/io/collection.py:357, in ImageCollection._check_imgnum(self, n)
    355     n = n % num
    356 else:
--> 357     raise IndexError(f"There are only {num} images in the collection")
    358 return n

IndexError: There are only 0 images in the collection

We can convert this list into a 3D image as follows:

im3d = np.zeros([*seq[0].shape, len(seq)])
for i, im in enumerate(seq):
    im3d[..., i] = im

And now we can use porespy to visualize it:

fig, ax = subplots(figsize=[5, 5])
ax.imshow(ps.visualization.show_planes(im3d));

Convert to Binary

In this tutorial we will not do anything to clean or denoise these images, as that deserves it’s own tutorial. Here we’ll just apply a threshold to convert a boolean with True indicating solid (the fibers above) and False indicating void.

First we need to look at the range of values in the image:

fig, ax = subplots(figsize=[5, 5])
ax.hist(im3d.flatten(), bins=25, edgecolor='k');

From the above histogram we can see that the values span from 22000 to 34000. We can just search for a value which does a decent job of segmenting the solid and void by trial and error:

fig, ax = subplots(figsize=[5, 5])
ax.imshow(im3d[..., 0] < 30000);

Because the image is noisy the result is not perfect. But let’s apply this threshold to the entire image now, then visualize the 3D stack:

im = im3d < 30000
fig, ax = subplots(figsize=[5, 5])
ax.imshow(ps.visualization.sem(im, axis=2));

Opening a 3D Stack

Assuming that the 2D images have been turned into a 3D image, as described above, it is straightforward to load them:

im2 = imageio.volread('../../docs/_static/images/Sample 1.tif')
ps.imshow(im2, axis=0);

Saving a 3D Image

Once you have created a nice image that you want to save, perhaps for visualization in Paraview, you can also use the imageio package:

imageio.volsave('my image.tif', np.array(im, dtype=np.int8))