Getting Started

This pages demonstrates how the core-ct library can be used, focusing on tools to visualize the core. In this notebook, we will go over an example workflow for visualizing and analyzing a core CT scan, including:

  • Creating a Core from dicom files

  • Displaying various views of the core

  • Taking a single slice from a core

  • Trimming unwanted space from the slice

  • Creating a brightness trace plot from a core slice

All code on this page is from a Jupyter notebook that can be downloaded here

Setup

1. Import Necessary Libraries

import os
import matplotlib.pyplot as plt
from core_ct import importers
from core_ct import visualize

2. Set Behavior of Matplotlib

The IPython kernel allows us to use magic commands to customize the behavior of our notebook. We will use a magic command to control the behavior of matplotlib.

  • To create non-interactive plots that appear inside the notebook, use %matplotlib inline. This is the default behavior, so you do not need to state it explicitly.

  • To create interactive plots that open in a new window, use %matplotlib. You may need to write this command twice to get the desired behavior. This also provides a basic GUI to tweak subplot spacing.

  • To create interactive plots that appear inside the notebook, use %matplotlib widget. This may cause more lag when interacting with the plots than opening it in a separate window.

You can also add magic commands throughout the notebook if you want different behaviors for different plots. Note that you cannot change the GUI backend mid-notebook, but you can change between inline and interactive. A full list of matplotlib magic options can be found here.

# Choose your desired plotting behavior by commenting/uncommenting lines in this cell

# %matplotlib inline

%matplotlib widget

# %matplotlib 
# %matplotlib

Load the Core Data

For this demo, we will be using sediment core CT scan data from the 2023 paper “The life and death of a subglacial lake in West Antarctica.” by Siegfried, Venturelli et al. [1]. The data used in this example can be downloaded here. For the curious, the full dataset can be downloaded from Zenodo, which you can find here. (Uncompressed size ~32 GB) [2].

Load this data using the function dicom from the importers module.

dir_path = os.path.join("siegfried_venturelli_core")
# the force option ignores files that produce errors (non-DICOMs)
# the ignore_file_extensions allows us to read files without the '.dcm' extension
my_core = importers.dicom(dir = dir_path, force = True, ignore_file_extensions = True)

Displaying the Core

Before we do any processing of the core, we want to get some information about it, including:

  • the dimensions of the 3D data array of the core scan’s pixels

  • how the core is oriented, that is, how the axes of the core scan correspond to the axes of the data array. We can use the display_core method to show us orthogonal views of the core for each axis. If you are using matplotlib’s interactive display, you can experiment with tweaking the subplot spacing.

print("The shape of the core is: ", my_core.shape()) 

# display orthogonal views of the core sliced along each axis
fig, (ax1, ax2, ax3) = visualize.display_core(my_core)
fig.set_figwidth(8)
plt.show()

The shape of the core is:  (1560, 1560, 34)

Notice that the core appears strangely stretched. This is because the pixel dimensions (the size of each pixel in mm along each axis) are unequal. We can choose to plot the core view in mm instead of pixels to remove the distortion by setting mm = True.

fig, (ax1, ax2, ax3) = visualize.display_core(my_core, mm = True)
fig.set_figwidth(8)
plt.show()

From these visualizations, we can see that axis 0 is the long axis of the core (down the page), axis 1 goes across the page, and axis 2 goes into the page. We also see that clearly we want to collapse along axis 2 to get a nice cross-section of the entire core. But first, we want to trim away unwanted data.

Trimming the Core

There’s a lot of extra space surrounding the Core which we want to trim away. Since this core is a cylinder, we will use the trim_radial function. We want to trim about axis 0, since that is the axis that gives us a coronal slice (circular core cross section).

If we don’t want to trim a core radially, numpy style indexing can be used. More information on that can be found here.

trimmed_core = my_core.trim_radial(axis=0, radius=30, y_center=770, z_center=16)
fig, (ax1, ax2, ax3) = visualize.display_core(trimmed_core, mm = True)
fig.set_figwidth(8)
plt.show()

Taking a Single Slice of the Core

Because the core looked reasonably centered, we will simply take a slice at the center of axis 2. We can display the slice using the display_slice method to ensure it looks good.

my_slice = trimmed_core.slice(axis = 2, loc = trimmed_core.shape()[2]//2)
# another way to take the slice is by using the bracket operator like this (up to personal preference):
# my_slice = trimmed_core[:, :, trimmed_core.shape()[2]//2]
visualize.display_slice(my_slice)
plt.show()

Trimming the Slice

The core slice looks like it was taken at a good location, but the sediment core looks like it got disrupted near the top and bottom. We want to trim that off before performing any analysis. We can visualize the possible trim lines before we actually perform the trim to make sure it’s in the right spot.

# horizontal trim
visualize.visualize_trim(slice = my_slice, 
                        axis = 0, 
                        loc_start = 500, 
                        loc_end = 200)
plt.show()

Now that we know where we want to trim our slice, we use the bracket operator ([]) on the Slice object to actually perform the trim. We can use the display_slice method to display the slice along with a colorbar. We can choose to either display the axes in millimeters or pixels. If you are using an interactive window, you can drag the colorbar to change the color mapping.

new_slice = my_slice[500:-200]
visualize.display_slice(new_slice, mm = True)
plt.show()

Creating a Brightness Trace

Now we will use the function display_slice_bt_std to create a graph of the slice next to the plot of the brightness trace across it and the standard deviation of the brightness.

fig, (ax1, ax2, ax3) = visualize.display_slice_bt_std(new_slice)
# you can change different plot elements attached to fig or the axes
fig.suptitle("Sediment Core Brightness Trace")
fig.set_figwidth(8)
plt.show()

Filtering the Slice

If you want to filter the slice to display only a certain density range, for example if you wanted to isolate high-density particles in a core, you can use the filter method on the Slice object. There is also a filter that operates on a Core object as well, that is particularly useful if you want to find the percent volume of the core within a certain density range. Note that filter takes a function that returns a boolean.

# range was approximated using the colorbar of previous plots
filtered_slice = new_slice.filter(lambda x: x > 1400 and x < 3500)
visualize.display_slice(filtered_slice)
plt.show()

References

[1] Siegfried, M. R., Venturelli, R. A., Patterson, M. O., Arnuk, W., Campbell, T. D., Gustafson, C. D., Michaud, A. B., Galton-Fenzi, B. K., Hausner, M. B., Holzschuh, S. N., Huber, B., Mankoff, K. D., Schroeder, D. M., Summers, P., Tyler, S., Carter, S. P., Fricker, H. A., Harwood, D. M., Leventer, A., Rosenheim, B. E., Skidmore, M. L., Priscu, J. C., and the SALSA Science Team. (2023). The life and death of a subglacial lake in West Antarctica. Geology. https://doi.org/10.1130/G50995.1

[2] Siegfried, M. R., Venturelli, R. A., Patterson, M. O., Arnuk, W., Campbell, Gustafson, Chloe D., C. D., Michaud, A. B., Galton-Fenzi, B. K., Hausner, M. B., Holzschuh, S. N., Huber, B., Mankoff, K. D., Schroeder, D. M., Summers, P. T., Tyler, S., Carter, S. P., Fricker, H. A., Harwood, D. M., Leventer, A., Rosenheim, B. E., Skidmore, M. L., Priscu, J. P., and the SALSA Science Team. (2023). Data for Siegfried*, Venturelli*, et al., 2023, Geology (1.0) [Data set]. Zenodo. https://doi.org/10.5281/ZENODO.7597019