Welcome to foamquant documentation

Important

This project is under construction. New functionalities are constantly added to this package. The project is directly related to the study [Schott2023].

foamquant is a toolbox specifically created for processing time series of 3D images of evolving liquid foams by using open source libraries Scikit-Image [vanderWalt2014] and SPAM [Stamani2020].

The objective is a greater accessibility of time-resolved liquid foam 3D analysis tools for the foam-physicis and food scientists communitites. We propose the following documentation: API and Jupyter Notebooks.

Installation: neither conda, nor pip installations are available yet.

Main dependencies: SPAM and scikit-image

About

Important

This project is under construction. New functionalities are constantly added to this package. The project is directly related to the study [Schott2023].

Overview

Current package structure. The functions in red are not yet included in the package.

Process

Functions for processing batch liquid foam images, all the steps from the raw-reconstructed images to the labelled images :

• Remove background (homogeneization)

• Phase segmentation (binarization)

• Masking (cylindrical or region of interest)

• Remove small objects and holes (volume threshold)

• Bubble segmentation (watershed)

• Remove edge bubbles (edge of a mask if provided)

From the left to the right: raw, phase segmented and bubble segmented images.

FromBinary

Functions to quantify the liquid fraction from a batch of phase segmented images.

The liquid fraction along a cartesian mesh can be returned structured or unstructured.

FromLabelled

Functions to quantify the bubbles regions properties from a batch of labelled images.

The function save the regions properties in a .csv

Tracking

Functions to track the bubbles and their properties from a batch of labelled images.

The color (from green to black) indicates the time index. The red points are the lost tracking positions.

Passage and Average

Functions to convert scalar, vectorial or tensorial properties from cartesian to cylindrical and spherical, and perform time/space averages.

In this example the displacement field is first expressed in a cylindrical basic and then averaged.

Two ways of measuring the internal strain field

• Shape field, defined in [Graner2008] and first used in [Raufaste2015]

• Texture field, defined in [Graner2008]

Label traking

The tracking method was inspired by ID-track presented in [Ando2013].

Tracking of five bubbles, showing various tracked properties: elastic internal strain, number of neighbours, velocity, and volume.

References

[vanderWalt2014]

19. van der Walt et al., scikit-image: Image processing in Python. PeerJ 2:e453 (2014) https://doi.org/10.7717/peerj.453

[Stamani2020]

Stamati et al., (2020). spam: Software for Practical Analysis of Materials. Journal of Open Source Software, 5(51), 2286, https://doi.org/10.21105/joss.02286

[Ando2013]

Andò,E. et al., Experimental micromechanics: grain-scale observation of sand deformation, Géotechnique Letters 2, 107–112, (2012) https://doi.org/10.1680/geolett.12.00027

[Hall2010]

19. 1. Hall et al., Discrete and continuum analysis of localised deformation in sand using X-ray μCT and volumetric digital image correlation. Géotechnique, 60(5), 315-322, (2010) https://doi.org/10.1680/geot.2010.60.5.315

[Graner2008] (1,2)

6. Graner et al., Discrete rearranging disordered patterns, part I: Robust statistical tools in two or three dimensions, Eur. Phys. J. E 25, 349–369 (2008) https://doi.org/10.1140/epje/i2007-10298-8

[Raufaste2015]

Raufaste, C. et al., Three-dimensional foam flow resolved by fast X-ray tomographic microscopy, EPL, 111, 38004, (2015) https://doi.org/10.1209/0295-5075/111/38004

[Schott2023]

6. Schott et al., Three-dimensional liquid foam flow through a hopper resolved by fast X-ray microtomography, Soft Matter, (2023) https://doi.org/10.1039/d2sm01299e

API

Process

FoamQuant.Process.BubbleSegmentation(image, SigSeeds=1, SigWatershed=1, watershed_line=False, radius_opening=None, verbose=False, esti_min_dist=None, compactness=None)

Perform the bubble watershed segmentation

Parameters:

• image (int numpy array) – 3D image

• SigSeeds (int) – Optional, Gaussian filter Sigma for the seeds

• SigWatershed (int) – Optional, Gaussian filter Sigma for the watershed

• watershed_line (Bool) – Optional, If True keep the 0-label surfaces between the segmented bubble regions

• radius_opening (None or int) – Optional, if not None, perform a radius opening operation on the labelled image with the given radius

• verbose (Bool) – Optional, if True, print progression steps of the segmentation

• esti_min_dist (None or float) – Optional, min distance between the seeds

• compactness (None or float) – Optional, compactness of the diffusion

Returns:

int numpy array

FoamQuant.Process.BubbleSegmentation_Batch(nameread, namesave, dirread, dirsave, imrange, verbose=False, n0=3, endread='.tiff', endsave='.tiff', writeparameters=False, Binning=None, SigSeeds=1, SigWatershed=1, watershed_line=False, esti_min_dist=None, compactness=None, radius_opening=None, ITK=False, ITKLevel=1)

Uses BubbleSegmentation function batchwise

Parameters:

• nameread (str) – read image name

• namesave (str) – saved image name

• dirread (str) – read image directory

• dirread – saved image directory

• imrange (list or numpy array) – image indexes array

• verbose (Bool) – if True, print the progression

• n0 (int) – number of digit for the saving index, default is 3

• endread (str) – read image file extension, default is ‘.tiff’

• endsave (str) – saved image file extension, default is ‘.tiff’

• writeparameters (str) – saved in a text file the segmentation parameters

• Binning (str) – saved image file extension, default is ‘.tiff’

• SigSeeds (int) – Optional, Gaussian filter Sigma for the seeds

• SigWatershed (int) – Optional, Gaussian filter Sigma for the watershed

• watershed_line (Bool) – Optional, If True keep the 0-label surfaces between the segmented bubble regions

• esti_min_dist (None or float) – Optional, min distance between the seeds

• compactness (None or float) – Optional, compactness of the diffusion

• radius_opening (None or int) – Optional, if not None, perform a radius opening operation on the labelled image with the given radius

• ITK (Bool) – If True, the ITKwatershed from SPAM is used

• ITKLevel (float) – Optional, default is 1

FoamQuant.Process.FastLocalThickness_Batch(nameread, namesave, dirread, dirsave, imrange, verbose=False, endread='.tiff', endsave='.tiff', n0=3, WalThickness=True, Separation=False, scale=1)

Uses localthickness function batchwise. IMPORTANT, please refer to Dahl, V. A. and Dahl A. B. work: Git-link February 2023: https://github.com/vedranaa/local-thickness.git

Parameters:

• nameread (str) – read image name

• namesave (str) – saved image name

• dirread (str) – read image directory

• dirread – saved image directory

• imrange (list or numpy array) – image indexes array

• verbose (Bool) – if True, print the progression

• endread (str) – read image file extension, default is ‘.tiff’

• endsave (str) – saved image file extension, default is ‘.tiff’

• n0 (int) – number of digit for the saving index, default is 3

• WalThickness (Bool) – If True, save the wall thickness (zeros in the binary image)

• Separation (Bool) – If True, save the separation thickess (ones in the binary image)

• scale (float) – Optional downscaling factor, default is 1

FoamQuant.Process.GetContacts(image, image_noedge, maximumCoordinationNumber=20, returnCoordImage=False)

Return labels [0], centroids [1], coordinations [2] of the no-edge image, (and coordination image [3])

Parameters:

• image (int numpy array) – full 3D image

• image_noedge (int numpy array) – 3D image with removed label at the edges

• maximumCoordinationNumber (int) – the maximum coordination number, default 20

• returnCoordinationImage (Bool) – if True, additionally returns image_noedge coordination image

Returns:

labels [0], centroids [1], coordinations [2] of the no-edge image, (and coordination image [3])

FoamQuant.Process.GetContacts_Batch(nameread, nameread_noedge, namesave, dirread, dirread_noedge, dirsave, imrange, verbose=False, endread='.tiff', endread_noedge='.tiff', endsave='.tiff', n0=3, save='all', maximumCoordinationNumber=20)

Uses GetContacts function batchwise

Parameters:

• nameread (str) – read image name

• namesave (str) – saved image name

• dirread (str) – read image directory

• dirread – saved image directory

• imrange (list or numpy array) – image indexes array

• verbose (Bool) – if True, print the progression

• endread (str) – read labeled image file extension, default is ‘.tiff’

• endread_noedge (str) – read labeled no-edge image file extension, default is ‘.tiff’

• endsave (str) – saved image file extension, default is ‘.tiff’

• n0 (int) – number of digit for the saving index, default is 3

• save (str or list of str) – ‘all’: save all, ‘coord’:save the coordination image, ‘image’: save contact image, table’: save the contact table or ‘pair’: save the contact pairs

• maximumCoordinationNumber (int) – the maximum coordination number, default 20

FoamQuant.Process.MaskCyl(image, rpercent=None)

Create a cylindrical mask of the size of the image along the Z axis

Parameters:

image (float numpy array) – 3D image

Returns:

int numpy array

FoamQuant.Process.Masking_Batch(nameread, namesave, dirread, dirsave, imrange, verbose=False, n0=3, endread='.tiff', endsave='.tiff')

Uses MaskCyl function batchwise

Parameters:

• nameread (str) – read image name

• namesave (str) – saved image name

• dirread (str) – read image directory

• dirread – saved image directory

• imrange (list or numpy array) – image indexes array

• verbose (Bool) – if True, print progression steps of the cleaning

• n0 (int) – number of digit for the saving index, default is 3

• endread (str) – read image file extension, default is ‘.tiff’

• endsave (str) – saved image file extension, default is ‘.tiff’

FoamQuant.Process.PhaseSegmentation(image, method='ostu_global', th=0.5, th0=0.3, th1=0.7, returnOtsu=False)

Perform the phase segmentation

Parameters:

• image (float numpy array) – 3D image

• method (str) – Optional, method for segmenting the phase, either ‘simple’ for simple threshold, ‘ostu_global’ for a global Otsu threshold, ‘niblack’, ‘sauvola’, or ‘sobel’, Default is ‘ostu_global’

• th (float) – Optional, given threshold for ‘simple’ method

• th0 (float) – Optional, given low threshold for ‘sobel’ method

• th1 (float) – Optional, given high threshold for ‘sobel’ method

• returnotsu (Bool) – Optional, if True, returns Otsu threshold for ‘ostu_global’ method

Returns:

int numpy array and float

FoamQuant.Process.PhaseSegmentation_Batch(nameread, namesave, dirread, dirsave, imrange, method='ostu_global', th=None, th0=None, th1=None, returnOtsu=False, verbose=False, ROIotsu=False, n0=3, endread='.tiff', endsave='.tiff')

Uses PhaseSegmentation function batchwise

Parameters:

• nameread (str) – read image name

• namesave (str) – saved image name

• dirread (str) – read image directory

• dirread – saved image directory

• imrange (list or numpy array) – image indexes array

• method (str) – Optional, method for segmenting the phase, either ‘simple’ for simple threshold, ‘ostu_global’ for a global Otsu threshold, ‘niblack’, ‘sauvola’, or ‘sobel’, Default is ‘ostu_global’

• th (float) – Optional, given threshold for ‘simple’ method

• th0 (float) – Optional, given low threshold for ‘sobel’ method

• th1 (float) – Optional, given high threshold for ‘sobel’ method

• returnOtsu (Bool) – Optional, if True, returns Otsu threshold for ‘ostu_global’ method

• verbose (Bool) – if True, print progression steps of the cleaning

• ROIotsu (list or numpy array) – list of length 6 defining the region of interest for determining the single Otsu threshold such as [zmin,zmax,ymin,ymax,xmin,xmax]

• n0 (int) – number of digit for the saving index, default is 3

• endread (str) – read image file extension, default is ‘.tiff’

• endsave (str) – saved image file extension, default is ‘.tiff’

FoamQuant.Process.RemoveBackground(image, method='white_tophat', radius=5)

Remove grey-scale image low frequency background

Parameters:

• image (float numpy array) – 3D image

• method (str) – method for removing the background, either ‘white_tophat’:white tophat filter or ‘remove_gaussian’: remove the Gaussian filtered image

• radius (int) – white_tophat kernel radius or sigma gaussian filter radius

Returns:

float numpy array

FoamQuant.Process.RemoveBackground_Batch(nameread, namesave, dirread, dirsave, imrange, method='white_tophat', radius=5, verbose=False, Binning=None, n0=3, endread='.tiff', endsave='.tiff')

Uses RemoveBackground function batchwise

Parameters:

• nameread (str) – read image name

• namesave (str) – saved image name

• dirread (str) – read image directory

• dirread – saved image directory

• imrange (list or numpy array) – image indexes array

• method (str) – method for removing the background, either ‘white_tophat’:white tophat filter or ‘remove_gaussian’: remove the Gaussian filtered image

• radius (int) – white_tophat kernel radius or sigma gaussian filter radius

• verbose (Bool) – if True, print progression steps of the cleaning

• Binning (None or int) – the binning number

• n0 (int) – number of digit for the saving index, default is 3

• endread (str) – read image file extension, default is ‘.tiff’

• endsave (str) – saved image file extension, default is ‘.tiff’

FoamQuant.Process.RemoveBackground_BatchHome(series, rawdir, prossdir, imrange, method='white_tophat', radius=5, weight=0.1, crop=None, bottom=None, verbose=False, Binning=None, n0rf=3, n0rs=3, n0w=3)

FoamQuant.Process.RemoveEdgeBubble(image, mask=None, rpercent=None, masktopbottom=None, returnmask=False, verbose=False)

Remove the bubbles on the image edges and in intersection with the masks (if given)

Parameters:

• image (int numpy array) – 3D image

• mask (None, Bool or int numpy array) – 3D image, if given, removes also the labels at the intersection with the mask

• rpercent (float) – If mask is True, will create a cylindrical mask with rpercent of the half-size of the image

• masktopbottom (array) – list as [zmin, zmax] for top and bottom edges masking

• returnmask (Bool) – if True, aditionally returns the mask

• verbose (Bool) – Optional, if True, print progression steps of the segmentation

Returns:

int numpy array image or [image, mask] arrays if returnmask is True

FoamQuant.Process.RemoveEdgeBubble_Batch(nameread, namesave, dirread, dirsave, imrange, verbose=False, endread='.tiff', endsave='.tiff', n0=3, maskcyl=False, rpercent=None, masktopbottom=None)

Uses RemoveEdgeBubble function batchwise

Parameters:

• nameread (str) – read image name

• namesave (str) – saved image name

• dirread (str) – read image directory

• dirread – saved image directory

• imrange (list or numpy array) – image indexes array

• verbose (Bool) – if True, print the progression

• endread (str) – read image file extension, default is ‘.tiff’

• endsave (str) – saved image file extension, default is ‘.tiff’

• n0 (int) – number of digit for the saving index, default is 3

• maskcyl (Bool) – if True, create a cylindrical mask

• rpercent (float) – If mask is True, will create a cylindrical mask with rpercent of the half-size of the image

• masktopbottom (array) – list as [zmin, zmax] for top and bottom edges masking

FoamQuant.Process.RemoveSpeckle(image, method='median', radius=1, weight=0.1)

Remove speckle from the image

Parameters:

• image (float numpy array) – 3D image

• method (str) – method for removing the speckle, either ‘median’, ‘gaussian’ or ‘tv_chambolle’

• radius (int) – kernel radius or sigma gaussian filter radius

• weight (int) – weight for tv_chambolle

Returns:

float numpy array

FoamQuant.Process.RemoveSpeckleBin(image, RemoveObjects=True, RemoveHoles=True, BinClosing=False, ClosingRadius=None, GiveVolumes=False, verbose=False, Vminobj=None, Vminhole=None)

Remove small objects and holes in binary image

Parameters:

• image (int numpy array) – 3D image

• RemoveObjects (Bool) – if True, removes the small objects

• RemoveHoles (Bool) – if True, removes the small holes

• BinClosing (Bool) – if True, perform a binnary closing of radius ClosingRadius

• ClosingRadius (int) – radius of the binnary closing

• GiveVolumes (Bool) – if True, returns in addition the used min volume thresholds for objects and holes

• verbose (Bool) – if True, print progression steps of the cleaning

• Vminobj (int) – if given the min volume threshold for the objects is not computed, Vminobj is used as the threshold

• Vminhole (int) – if given the min volume threshold for the holes is not computed, Vminobj is used as the threshold

Returns:

int numpy array, int, int

FoamQuant.Process.RemoveSpeckleBin_Batch(nameread, namesave, dirread, dirsave, imrange, verbose=False, RemoveObjects=True, RemoveHoles=True, Cobj=0.5, Chole=0.5, BinClosing=False, ClosingRadius=None, n0=3, endread='.tiff', endsave='.tiff')

Uses RemoveSpeckleBin function batchwise

Parameters:

• nameread (str) – read image name

• namesave (str) – saved image name

• dirread (str) – read image directory

• dirread – saved image directory

• imrange (list or numpy array) – image indexes array

• verbose (Bool) – if True, print the progression

• RemoveObjects (Bool) – If True, remove the small objects with a volume below Cobj*MaxObjVol

• RemoveHoles (Bool) – If True, remove the small holes with a volume below Chole*MaxholeVol

• Cobj (float) – volume thresholding coefficient for removing the small objects

• Chole (float) – volume thresholding coefficient for removing the small holes

• BinClosing (float) – volume thresholding coefficient for removing the small holes

• n0 (int) – number of digit for the saving index, default is 3

• endread (str) – read image file extension, default is ‘.tiff’

• endsave (str) – saved image file extension, default is ‘.tiff’

FromBinary

FoamQuant.FromBinary.LiqFrac_Batch(nameread, namesave, dirread, dirsave, imrange, TypeGrid='Global', Nz=None, Ny=None, Nx=None, Nr=None, Naz=None, crop=None, Mask=False, verbose=False, endread='.tif', endsave='.pkl', n0=3, structured=True)

Read 3D binary images and save liquid fraction informations in series (for loop). Save liquid fraction dictionary as pickle: {“crop”, “1D, 2D or 3D grid”,”lf”}

Parameters:

• series (str) – series name

• series – read image directory

• savedir (str) – save liquid fraction directory

• imrange (int numpy array) – image index range

• TypeGrid (str) – Optional, type of method: ‘Global’ for global liquid fraction, ‘CartesMesh’ for cartesian grid liquid fraction, ‘CylMesh’ for cylindrical grid liquid fraction

• Nz (int) – Optional, number of sub-regions along z

• Ny (int) – Optional, number of sub-regions along y

• Nx (int) – Optional, number of sub-regions along x

• Nr (int) – Optional, number of sub-regions along r

• crop ([zmin, zmax, ymin, ymax, xmin, xmax] int array) – Optional, Study crop region inside the image

• Mask (int numpy array) – Optional, 3D image

• verbose (int numpy array) – if True or 1, print progress image by image. If 10, print additionally the extracted liquid fraction dictionary

Returns:

None

FoamQuant.FromBinary.LiqFrac_CartesMesh(image, Nz, Ny, Nx, crop=None, Mask=None, structured=True)

Return a 3D zyx grid with its corresponding non-overlapping subvolume (cuboids) liquid fraction

Parameters:

• image (int numpy array) – 3D image

• Nz (int) – number of sub-regions along z

• Ny (int) – number of sub-regions along y

• Nx (int) – number of sub-regions along x

• crop ([zmin, zmax, ymin, ymax, xmin, xmax] int array) – Study crop region inside the image

• Mask (int numpy array) – 3D image

Returns:

int numpy array

FoamQuant.FromBinary.LiqFrac_Glob(image, Nz, Nr, crop=None, Mask=None)

Return the global liquid fraction of an image

Parameters:

• image (int numpy array) – 3D image

• crop ([zmin, zmax, ymin, ymax, xmin, xmax] int array) – Study crop region inside the image

• Mask (int numpy array) – 3D image

Returns:

int numpy array

FromLabelled

FoamQuant.FromLabelled.RegionProp(image, field=False)

Return basic region properties from a labeled image: labels [0], centroids [1], volumes [2], (inertia components [3])

Parameters:

• image (int numpy array) – 3D image

• IncludeInertia (Bool) – if True, also return inertial components

Returns:

array of labels [0], centroids [1], volumes [2], (inertia components [3])

FoamQuant.FromLabelled.RegionProp_Batch(nameread, namesave, dirread, dirsave, imrange, verbose=False, endread='.tif', endsave='.tsv', n0=3, field=False)

Run RegionProp function on a batch of images and save the outputs as .tsv

Parameters:

• readdir (str) – Labeled images folder

• readdir – folder to save the .tsv doc

• readend (str) – tiff image saving end, default is ‘.tiff’

• imrange (int array) – list of image indexes

• IncludeInertia (Bool) – if True, also return inertial components

• verbose (Bool) – if True, print verbose including the number of labels

FromContact

FoamQuant.FromContact.ReadContactTable(nameread, dirread, imrange, verbose=False, endread='.tsv', n0=3, maximumCoordinationNumber=20)

Read Contact tables and return a dictionary

Parameters:

• nameread (str) – read image name

• dirread (str) – read image directory

• imrange (list or numpy array) – image indexes array

• verbose (Bool) – if True, print progression

• endread (str) – read RegionProperties file extension, default is ‘.tsv’

• n0 (int) – number of digit for the saving index, default is 3

• maximumCoordinationNumber (int) – number of digit for the saving index, default is 3

Returns:

contact table dictionary {‘lab’,’lab_noedge’,’Z’, ‘z’,’y’,’x’, ‘labs’, ‘conts’}

FoamQuant.FromContact.Texture(Table, verbose=False)

From a contact table dictionary, compute the texture

Parameters:

• Table (dict) – contact table dictionary

• verbose (Bool) – if True, print progression

Returns:

lab withedge, lab withoutedge, centroid, radius, S1,S2,S3, S1z,S1y,S1x, S2z,S2y,S2x, S3z,S3y,S3x, U1,U2,U3, U, type

FoamQuant.FromContact.Texture_Batch(nameread, namesave, dirread, dirsave, imrange, verbose=False, endsave='.tsv', n0=3, field=False)

Run Texture function on a batch of images and save the outputs as .tsv

Parameters:

• readdir (str) – Labeled images folder

• readdir – folder to save the .tsv doc

• readend (str) – tiff image saving end, default is ‘.tiff’

• imrange (int array) – list of image indexes

• IncludeInertia (Bool) – if True, also return inertial components

• verbose (Bool) – if True, print verbose including the number of labels

Tracking

FoamQuant.Tracking.Combine_Tracking(nameread, dirread, imrange, verbose=False, endread='.tsv', n0=3)

Combine the LabelTracking result tables for a full time-series tracking of the labels present in the first image

Parameters:

• nameread (str) – read image name

• dirread (str) – read image directory

• imrange (list or numpy array) – image indexes array

• verbose (Bool) – if True, print progression

• endread (str) – read RegionProperties file extension, default is ‘.tsv’

• n0 (int) – number of digit for the saving index, default is 3

Returns:

combined dictionary {‘lab’,’labtransl’,’laborig’,’z’,’y’,’x’, ‘dz’, ‘dy’, ‘dx’, ‘vol’, ‘rad’, ‘area’, ‘sph’, ‘volfit’, ‘U’, ‘type’, ‘Utype’}

FoamQuant.Tracking.LabelTracking(Prop1, Prop2, searchbox=[-10, 10, -10, 10, -10, 10], Volpercent=0.05)

Return tracked properties from RegionProperties tables

Parameters:

• Prop1 (dic type) – RegionProperties table at time 1

• Prop1 – RegionProperties table at time 2

• searchbox (searchbox size (1,6) array corresponding to [zmin,zmax,ymin,ymax,xmin,xmax], default is [-10,10,-10,10,-10,10]) – method for removing the background, either ‘white_tophat’:white tophat filter or ‘remove_gaussian’: remove the Gaussian filtered image

• Volpercent (float) – allowed volume variation percentage from one image to the next, default is 0.05 (5 %)

Returns:

numpy arrays: Found label, Found coordinate, Number of candidates, Lost label, Lost coordinates, Volume,Radius,Area,Sphericity,Volume fit, U, type

FoamQuant.Tracking.LabelTracking_Batch(nameread, namesave, dirread, dirsave, imrange, verbose=False, endread='.tsv', endsave='.tsv', n0=3, searchbox=[-10, 10, -10, 10, -10, 10], Volpercent=0.05)

Run LabelTracking batchwise

Parameters:

• nameread (str) – read image name

• namesave (str) – saved image name

• dirread (str) – read image directory

• dirsave (str) – saved image directory

• imrange (list or numpy array) – image indexes array

• verbose (Bool) – if True, print progression

• endread (str) – read RegionProperties file extension, default is ‘.tsv’

• endsave (str) – saved Tracking file extension, default is ‘.tsv’

• n0 (int) – number of digit for the saving index, default is 3

• searchbox (searchbox size (1,6) array corresponding to [zmin,zmax,ymin,ymax,xmin,xmax], default is [-10,10,-10,10,-10,10]) – method for removing the background, either ‘white_tophat’:white tophat filter or ‘remove_gaussian’: remove the Gaussian filtered image

• Volpercent (float) – allowed volume variation percentage from one image to the next, default is 0.05 (5 %)

Returns:

numpy arrays: Found label, Found coordinate, Number of candidates, Lost label, Lost coordinates, Volume,Radius,Area,Sphericity,Volume fit, U, type

FoamQuant.Tracking.Read_LabelTracking(nameread, dirread, imrange, verbose=False, endread='.tsv', n0=3)

Read the full batch of LabelTracking result tables

Parameters:

• nameread (str) – read image name

• dirread (str) – saved image name

• imrange (list or numpy array) – image indexes array

• verbose (Bool) – if True, print progression

• endread (str) – read RegionProperties file extension, default is ‘.tsv’

• n0 (int) – number of digit for the saving index, default is 3

Returns:

dictionary {‘lab1’,’z1’,’y1’,’x1’,’lab2’,’z2’,’y2’,’x2’, ‘Count’, ‘dz’, ‘dy’, ‘dx’, ‘Count’,’vol1’,’rad1’,’area1’,’sph1’,’volfit1’,’U1’,’type1’,’Utype1’,’vol2’,’rad2’,’area2’,’sph2’’volfit2’,’U2’’type2’’Utype2’}

FoamQuant.Tracking.Save_Tracking(combined, namesave, dirsave, verbose=False, endsave='.csv')

Save the combined tracking dictionary as a table

Parameters:

• combined (dict) – combined tracking dictionary

• namesave (str) – save image name

• dirsave (str) – save image directory

• verbose (Bool) – if True, print progression

• endread (str) – read RegionProperties file extension, default is ‘.csv’

Passage

FoamQuant.Passage.Azimuth(X, Y)

Return the aziuthal angle of a point from its coordinates x and y

Parameters:

• X (float) – coordinate x

• Y (float) – coordinate y

Returns:

float, cylindrical azimuthal angle phi

FoamQuant.Passage.Cartesian2Cylindrical_Point(Coord_Cartesian, CoordAxis=[0, 1008, 1008], CylAxisZ=[1, 0, 0], CylAxisY=[0, 1, 0], CylAxisX=[0, 0, 1])

Return unstructured cylindrical coordinates array (N,3) from unstructured cartesian coordinates (N,3)

Parameters:

• Coord_Cartesian (numpy array) – unstructured cartesian coordinates, expressed in the image basis [z,y,x]

• CoordAxis (numpy array) – Coordinates of the axis of rotation in the image basis [z,y,x]

• CylAxisZ (numpy array) – Direction vector of the axis of rotation Z expressed in the image basis [z,y,x]

• CylAxisY (numpy array) – Direction vector of the axis Y expressed in the image basis [z,y,x]

• CylAxisX (numpy array) – Direction vector of the axis X expressed in the image basis [z,y,x]

Returns:

unstructured cylindrical coordinates (N,3), expressed in the cylindrical basis [r,phi,z]

FoamQuant.Passage.Cartesian2Cylindrical_Tensor(Coord_Cartesian, T_Cartesian, CoordAxis=[0, 1008, 1008], CylAxisZ=[1, 0, 0], CylAxisY=[0, 1, 0], CylAxisX=[0, 0, 1])

Return unstructured cylindrical coordinates (N,3) and tensors (N,3,3) from unstructured cartesian coordinates (N,3) and tensors (N,3,3)

Parameters:

• Coord_Cartesian (numpy array) – unstructured cartesian coordinates, expressed in the image basis [z,y,x]

• CoordAxis (numpy array) – Coordinates of the axis of rotation in the image basis [z,y,x]

• CylAxisZ (numpy array) – Direction vector of the axis of rotation Z expressed in the image basis [z,y,x]

• CylAxisY (numpy array) – Direction vector of the axis Y expressed in the image basis [z,y,x]

• CylAxisX (numpy array) – Direction vector of the axis X expressed in the image basis [z,y,x]

Returns:

unstructured cylindrical coordinates [0] and tensors [1], expressed in the cylindrical basis [r,phi,z]

FoamQuant.Passage.Cartesian2Cylindrical_Vector(Coord_Cartesian, V_Cartesian, CoordAxis=[0, 1008, 1008], CylAxisZ=[1, 0, 0], CylAxisY=[0, 1, 0], CylAxisX=[0, 0, 1])

Return unstructured cylindrical coordinates (N,3) and vectors (N,3) from unstructured cartesian coordinates (N,3) and vectors (N,3)

Parameters:

• Coord_Cartesian (numpy array) – unstructured cartesian coordinates, expressed in the image basis [z,y,x]

• CoordAxis (numpy array) – Coordinates of the axis of rotation in the image basis [z,y,x]

• CylAxisZ (numpy array) – Direction vector of the axis of rotation Z expressed in the image basis [z,y,x]

• CylAxisY (numpy array) – Direction vector of the axis Y expressed in the image basis [z,y,x]

• CylAxisX (numpy array) – Direction vector of the axis X expressed in the image basis [z,y,x]

Returns:

unstructured cylindrical coordinates [0] and vectors [1], expressed in the cylindrical basis [r,phi,z]

FoamQuant.Passage.Cartesian2Spherical_Point(Coord_Cartesian, CoordAxis=[0, 1008, 1008], CylAxisZ=[1, 0, 0], CylAxisY=[0, 1, 0], CylAxisX=[0, 0, 1])

Return unstructured spherical coordinates array (N,3) from unstructured cartesian coordinates (N,3)

Parameters:

• Coord_Cartesian (numpy array) – unstructured cartesian coordinates, expressed in the image basis [z,y,x]

• CoordAxis (numpy array) – Coordinates of the axis of rotation in the image basis [z,y,x]

• CylAxisZ (numpy array) – Direction vector of the axis of rotation Z expressed in the image basis [z,y,x]

• CylAxisY (numpy array) – Direction vector of the axis Y expressed in the image basis [z,y,x]

• CylAxisX (numpy array) – Direction vector of the axis X expressed in the image basis [z,y,x]

Returns:

unstructured cylindrical coordinates (N,3), expressed in the cylindrical basis [r,phi,z]

FoamQuant.Passage.Cartesian2Spherical_Tensor(Coord_Cartesian, T_Cartesian, CoordAxis=[0, 1008, 1008], CylAxisZ=[1, 0, 0], CylAxisY=[0, 1, 0], CylAxisX=[0, 0, 1])

Return unstructured spherical coordinates (N,3) and tensors (N,3,3) from unstructured cartesian coordinates (N,3) and tensors (N,3,3)

Parameters:

• Coord_Cartesian (numpy array) – unstructured cartesian coordinates, expressed in the image basis [z,y,x]

• CoordAxis (numpy array) – Coordinates of the axis of rotation in the image basis [z,y,x]

• CylAxisZ (numpy array) – Direction vector of the axis of rotation Z expressed in the image basis [z,y,x]

• CylAxisY (numpy array) – Direction vector of the axis Y expressed in the image basis [z,y,x]

• CylAxisX (numpy array) – Direction vector of the axis X expressed in the image basis [z,y,x]

Returns:

unstructured cylindrical coordinates [0] and tensors [1], expressed in the cylindrical basis [r,phi,z]

FoamQuant.Passage.Cartesian2Spherical_Vector(Coord_Cartesian, V_Cartesian, CoordAxis=[0, 1008, 1008], CylAxisZ=[1, 0, 0], CylAxisY=[0, 1, 0], CylAxisX=[0, 0, 1])

Return unstructured spherical coordinates (N,3) and vectors (N,3) from unstructured cartesian coordinates (N,3) and vectors (N,3)

Parameters:

• Coord_Cartesian (numpy array) – unstructured cartesian coordinates, expressed in the image basis [z,y,x]

• CoordAxis (numpy array) – Coordinates of the axis of rotation in the image basis [z,y,x]

• CylAxisZ (numpy array) – Direction vector of the axis of rotation Z expressed in the image basis [z,y,x]

• CylAxisY (numpy array) – Direction vector of the axis Y expressed in the image basis [z,y,x]

• CylAxisX (numpy array) – Direction vector of the axis X expressed in the image basis [z,y,x]

Returns:

unstructured cylindrical coordinates [0] and vectors [1], expressed in the cylindrical basis [r,phi,z]

FoamQuant.Passage.CylCoord(Coord_Cartesian, CoordAxis=[0, 1008, 1008], CylAxisZ=[1, 0, 0], CylAxisY=[0, 1, 0], CylAxisX=[0, 0, 1])

Return the cylindrical coordinate of a point from its cartesian coordinates

Parameters:

• Coord_Cartesian (numpy array) – cartesian coordinates in the image basis [z,y,x]

• CoordAxis (numpy array) – Coordinates of the axis of rotation in the image basis [z,y,x]

• CylAxisZ (numpy array) – Direction vector of the axis of rotation Z expressed in the image basis [z,y,x]

• CylAxisY (numpy array) – Direction vector of the axis Y expressed in the image basis [z,y,x]

• CylAxisX (numpy array) – Direction vector of the axis X expressed in the image basis [z,y,x]

Returns:

cylindrical coordinates [r,phi,z]

FoamQuant.Passage.CylR(X, Y)

Return the cylindrical radius r of a point from its coordinates x and y

Parameters:

• X (float) – coordinate x

• Y (float) – coordinate y

Returns:

float, cylindrical radius r

FoamQuant.Passage.Polar(Z, spheR)

Return the spherical polar angle of a point from its coordinates x, y and z

Parameters:

• X (float) – coordinate x

• Y (float) – coordinate y

• Z (float) – coordinate z

Returns:

float, spherical radius r

FoamQuant.Passage.Pzyx2phithetar(Azi, Polar)

Return the spherical passage matrix, from the spherical aziuthal and polar angles Azi and Polar

Parameters:

• Azi (float) – azimuthal angle Azi

• Polar (float) – polar angle Polar

Returns:

Passage matrix from (z,y,x) to (phi,theta,r), (3,3) numpy array

FoamQuant.Passage.Pzyx2zphir(Azi)

Return the cylindrical passage matrix from the cylindrical aziuthal angle Azi

Parameters:

phi (float) – azimuthal angle Azi

Returns:

Passage matrix from (z,y,x) to (z,phi,r), (3,3) numpy array

FoamQuant.Passage.SpheCoord(Coord_Cartesian, CoordAxis=[0, 1008, 1008], CylAxisZ=[1, 0, 0], CylAxisY=[0, 1, 0], CylAxisX=[0, 0, 1])

Return the spherical coordinates of a point from its cartesian coordinates

Parameters:

• Coord_Cartesian (numpy array) – cartesian coordinates in the image basis [z,y,x]

• CoordAxis (numpy array) – Coordinates of the axis of rotation in the image basis [z,y,x]

• CylAxisZ (numpy array) – Direction vector of the axis of rotation Z expressed in the image basis [z,y,x]

• CylAxisY (numpy array) – Direction vector of the axis Y expressed in the image basis [z,y,x]

• CylAxisX (numpy array) – Direction vector of the axis X expressed in the image basis [z,y,x]

Returns:

spherical coordinates [r,phi,z]

FoamQuant.Passage.SpheR(X, Y, Z)

Return the spherical radius of a point from its coordinates x, y and z

Parameters:

• X (float) – coordinate x

• Y (float) – coordinate y

• Z (float) – coordinate z

Returns:

float, spherical radius r

Average

FoamQuant.Average.Grid_Pavg(Coord, P, Range, N, NanFill=True, verbose=False, structured=True)

Return averaged scalar field over a 3D grid

Parameters:

• Coord (numpy array) – unstructured point coordinates (N,3)

• P (numpy array) – unstructured scalar array (N,1)

• Range (numpy array) – Averaging 3D ranges, such as [zmin,zmax,ymin,ymax,xmin,xmax]

• N (numpy array) – Number of points along each dimension, such as [Nz,Ny,Nx]

• NanFill (Bool) – If True fill with NaN when the Count is equal to 0 (else fill with zeros)

• verbose (Bool) – If True, print the averaging Cartesian grid

• structured (Bool) – If True, return a structured scalar averaged field (Nz,Ny,Nx,1), else return an unstructured field (Nz*Ny*Nx,1)

Returns:

[0] [grid0,grid1,grid2], [1] Grid coordinates, [2] Mean scalar field, [3] Std scalar field, [4] Count

FoamQuant.Average.Grid_Tavg(Coord, T, Range, N, NanFill=True, verbose=False, structured=True)

Return averaged tensor field over a 3D grid

Parameters:

• Coord (numpy array) – unstructured point coordinates (N,3)

• P (numpy array) – unstructured tensor array (N,3,3)

• Range (numpy array) – Averaging 3D ranges, such as [zmin,zmax,ymin,ymax,xmin,xmax]

• N (numpy array) – Number of points along each dimension, such as [Nz,Ny,Nx]

• NanFill (Bool) – If True fill with NaN when the Count is equal to 0 (else fill with zeros)

• verbose (Bool) – If True, print the averaging Cartesian grid

• structured (Bool) – If True, return a structured tensor averaged field (Nz,Ny,Nx,3,3), else return an unstructured field (Nz*Ny*Nx,3,3)

Returns:

[0] [grid0,grid1,grid2], [1] Grid coordinates, [2] Mean tensor field, [3] Std tensor field, [4] Count

FoamQuant.Average.Grid_Vavg(Coord, V, Range, N, NanFill=True, verbose=False, structured=True)

Return averaged vector field over a 3D grid

Parameters:

• Coord (numpy array) – unstructured point coordinates (N,3)

• P (numpy array) – unstructured vector array (N,3)

• Range (numpy array) – Averaging 3D ranges, such as [zmin,zmax,ymin,ymax,xmin,xmax]

• N (numpy array) – Number of points along each dimension, such as [Nz,Ny,Nx]

• NanFill (Bool) – If True fill with NaN when the Count is equal to 0 (else fill with zeros)

• verbose (Bool) – If True, print the averaging Cartesian grid

• structured (Bool) – If True, return a structured vector averaged field (Nz,Ny,Nx,1), else return an unstructured field (Nz*Ny*Nx,3)

Returns:

[0] [grid0,grid1,grid2], [1] Grid coordinates, [2] Mean vector field, [3] Std vector field, [4] Count

Figure

FoamQuant.Figure.Cut3D(image, zcut=False, ycut=False, xcut=False, showcuts=False, showaxes=False, cmap='gray', interpolation=None, figblocksize=5, returnfig=False, vmin=None, vmax=None, printminmax=False, colorbars=False)

Plot a 3x1 figure showing three orthogonal cross-sections of the 3D image.

Parameters:

• image (int) – 3D numpy array

• zcut (int or False) – Optional z cut value

• ycut (int or False) – Optional y cut value

• xcut (int or False) – Optional x cut value

• showcuts (Bool) – Optional plot the orthogonal cuts

• showaxes (Bool) – Optional plot the axes

• cmap (str or cmap type) – Optional the color map used for the cuts, Default cmap = ‘gray’

• interpolation (str or None) – Optional type of interpolation, Default interpolation = None

• figblocksize (float) – Optional size of the subfigure, Default figblocksize = 5

• returnfig (Bool) – Optional, if should return the figure, if not returns None

• vmin (Bool) – Optional, min value for the color range

• vmax (Bool) – Optional, max value for the color range

• printminmax (Bool) – Optional, print min and max for the whole image and the three projections

Returns:

None or fig type

FoamQuant.Figure.Histogram(image, histtitle=False)

Plot a 1x1 grey value histogram.

Parameters:

image (numpy array) – 3D image.

Returns:

None

FoamQuant.Figure.LinCmap(vmin=0, vmax=10, first_color='b', last_color='r', verbose=True)

Creates a linear colormap for matplotlib.

Parameters:

• vmin (float) – min value for the color-range

• vmax (float) – max value for the color-range

• first_color_black (str or matplotlib color) – first color

• last_color_black (str or matplotlib color) – last color

• verbose (Bool) – If True, prints the number of labels and shows the colormap.

Returns:

matplotlib colormap

FoamQuant.Figure.Proj3D(image, showaxes=False, cmap='gray', interpolation=None, figblocksize=5, returnfig=False, vmin=None, vmax=None, printminmax=False, colorbars=False)

Plot a 3x1 figure showing three orthogonal projections of the 3D image.

Parameters:

• image (numpy array) – 3D image

• showaxes (Bool) – Optional plot the axes

• cmap (str or cmap type) – Optional the color map used for the projections, Default cmap = ‘gray’

• interpolation (str or None) – Optional the type of interpolation, Default interpolation = None

• figblocksize (float) – Optional size of the subfigure, Default figblocksize = 5

• returnfig (Bool) – Optional if should return the figure, if not returns None

• vmin (Bool) – Optional, min value for the color range

• vmax (Bool) – Optional, max value for the color range

• printminmax (Bool) – Optional, print min and max for the whole image and the three projections

Returns:

None or fig

FoamQuant.Figure.RandomCmap(nlabels, type='bright', first_color_black=True, last_color_black=False, verbose=True, GiveRange=None)

Creates a random colormap for matplotlib. Reference: copied from https://github.com/delestro/rand_cmap

Parameters:

• nlabels (int) – Number of labels (size of colormap)

• type (str) – ‘bright’ for strong colors, ‘soft’ for pastel colors

• first_color_black (Bool) – Option to use first color as black, True or False

• last_color_black (Bool) – Option to use last color as black, True or False

• verbose (Bool) – Prints the number of labels and shows the colormap. True or False

Returns:

matplotlib colormap

FoamQuant.Figure.ellipse_plot(Fig, X, Y, Vect, Val, scale_factor=1, mirror=False)

Plot internal 2D structured strain field map.

Parameters:

• Fig (figure type) – matplotlib.pyplot figure

• X (numpy array) – structured (N,M) horyzontal positions

• Y (numpy array) – structured (N,M) vertical positions

• Vect (numpy array) – structured in-plane eigenvectors [U+x,U+y,U-x,U-y], (N,M,4)

• Val (numpy array) – structured in-plane eigenvalues [U+,U-,Uphi], (N,M,3)

• scale_factor (float) – Optional, by default 1

• mirror (Bool) – Optional, if True the plot is horyzontally flipped

Returns:

None

Movie

VTK

FoamQuant.VTK.json_rand_dictionary(Ncolors, namecmap, first_color_black=True)

Save a json random colormap to be used with ParaView or Tomviz.

Parameters:

• Ncolors (int) – Number of labels (size of colormap)

• type (str) – ‘bright’ for strong colors, ‘soft’ for pastel colors

• first_color_black (Bool) – Option to use first color as black, True or False

Returns:

None

FoamQuant.VTK.writeGlyphsVTK(coordinates, pointData, fileName='spam.vtk')

Write a plain text glyphs vtk.

Parameters:

• coordinates ((N,3) array of float) – Coordinates of the centre of all n glyphs

• pointData (dict {'field1name':field1,'field2name':field2, ...}) – (N,1) arrays for scalar values, (N,3) for vector values and (N,3) for tensor values

• fileName (int) – Optional name of the output file. By default=’.vtk’

Returns:

str index

Helper

FoamQuant.Helper.RangeList(i1, i2, verbose=False)

Return an index range

Parameters:

• i1 (int) – first index

• i2 (int) – last index

• verbose (Bool) – If True, print the range

Returns:

int numpy array

FoamQuant.Helper.strindex(i, n0)

Return str index written on n0 digit

Parameters:

• i (int) – index

• n0 (int) – number of 0 digit

Returns:

str index

Examples

Liquid foam flow: basic analysis

Binder tutorial:

Jupyter example:

Part 1: Image processing

# Import Libraries
import numpy as np
import matplotlib.pyplot as plt
from tifffile import imread, imsave
import skimage.measure
import pickle as pkl
import os

# Import FoamQuant library
from FoamQuant import *

# Set matplotlib default font size
plt.rc('font', size=20)

# Create the processing pipeline
ProcessPipeline = ['P1_Raw','P2_PhaseSegmented','P3_Cleaned','P4_BubbleSegmented','P5_BubbleNoEdge']

for Pi in ProcessPipeline:
    if  os.path.exists(Pi):
        print('path already exist:',Pi)
    else:
        print('Created:',Pi)
        os.mkdir(Pi)

path already exist: P1_Raw
path already exist: P2_PhaseSegmented
path already exist: P3_Cleaned
path already exist: P4_BubbleSegmented
path already exist: P5_BubbleNoEdge

A) The raw image

# Read/Save image names and directories
nameread = 'Raw_'
namesave = 'PhaseSegmented_'
dirread = ProcessPipeline[0]+'/'
dirsave = ProcessPipeline[1]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

# Read the first image of the series
Raw = imread(dirread+nameread+strindex(imrange[0], 3)+'.tiff')

# Show a 3D-cut view of the volume
Cut3D(Raw,
      showcuts=True,
      showaxes=True,
      figblocksize=7,
      zcut=15,       # tune this parrameter if you wish
      ycut=False,    # tune this parrameter if you wish
      xcut=False)    # tune this parrameter if you wish

/gpfs/offline1/staff/tomograms/users/flosch/Rheometer_Jupyter/FoamQuant/Figure.py:78: UserWarning: This figure was using constrained_layout==True, but that is incompatible with subplots_adjust and or tight_layout: setting constrained_layout==False.
  plt.tight_layout()

B) Phase segmentation

# Otsu simple threshold phase segmentation of the whole series
th = PhaseSegmentation_Batch(nameread,
                             namesave,
                             dirread,
                             dirsave,
                             imrange,
                             method='ostu_global',
                             returnOtsu=True,
                             verbose=True,
                             n0=3,
                             endread='.tiff',
                             endsave='.tiff')

PhaseSegmented_ 1: done

PhaseSegmented_ 2: done

PhaseSegmented_ 3: done

PhaseSegmented_ 4: done

PhaseSegmented_ 5: done

PhaseSegmented_ 6: done

PhaseSegmented_ 7: done

PhaseSegmented_ 8: done

PhaseSegmented_ 9: done

PhaseSegmented_ 10: done

# Otsu thresholds
print('Otsu thresholds:',th)

Otsu thresholds: [137, 140, 138, 141, 141, 141, 143, 143, 139, 139]

Let’s see the result…

# Read the first image of the series
Seg = imread(dirsave+namesave+strindex(imrange[0], 3)+'.tiff')

zcut=15       # tune this parrameter if you wish
ycut=False    # tune this parrameter if you wish
xcut=False    # tune this parrameter if you wish
cmap='jet'    # tune this parrameter if you wish: e.g. 'bone'

# Show a 3D-cut view of the volume
Cut3D(Seg, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)             # Phase segmented image
Cut3D((Seg>0)*Raw, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)     # Phase segmented image * Raw image
Cut3D((1-Seg)*Raw, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)     # (1-Phase segmented image) * Raw image

C) Remove small holes & regions

# Read/Save image names and directories
nameread = 'PhaseSegmented_'
namesave = 'Cleaned_'
dirread = ProcessPipeline[1]+'/'
dirsave = ProcessPipeline[2]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

Remove all holes and objects with: - Vobj < Cobj * max(Vobj) - Vhole < Chole * max(Vhole)

Since in liquid foam images, the liquid and gas phases both consist of unique large regions, Cobj and Chole can be strict (large thresholds). All the other smaller regions are often due to imaging artefacts.

# remove holes and objects
RemoveSpeckleBin_Batch(nameread,
                       namesave,
                       dirread,
                       dirsave,
                       imrange,
                       verbose=True,
                       endread='.tiff',
                       endsave='.tiff',
                       n0=3,
                       Cobj=0.1,  # tune this parrameter if you wish
                       Chole=0.1) # tune this parrameter if you wish

Before: Nobj 9
After: Nobj 1
Before: Nobj 20
After: Nobj 1
First image (vox): maxObj 13524383 maxHole 2351568
Thresholds (vox): thrObj 1352438 thrHole 235157

Before: Nhol 9
After: Nhol 1
Before: Nhol 20
After: Nhol 1
Cleaned_001: done

Before: Nhol 6
After: Nhol 1
Before: Nhol 21
After: Nhol 1
Cleaned_002: done

Before: Nhol 9
After: Nhol 1
Before: Nhol 26
After: Nhol 1
Cleaned_003: done

Before: Nhol 8
After: Nhol 1
Before: Nhol 28
After: Nhol 1
Cleaned_004: done

Before: Nhol 3
After: Nhol 1
Before: Nhol 22
After: Nhol 1
Cleaned_005: done

Before: Nhol 4
After: Nhol 1
Before: Nhol 41
After: Nhol 1
Cleaned_006: done

Before: Nhol 8
After: Nhol 1
Before: Nhol 19
After: Nhol 1
Cleaned_007: done

Before: Nhol 6
After: Nhol 1
Before: Nhol 23
After: Nhol 1
Cleaned_008: done

Before: Nhol 1

/gpfs/offline1/staff/tomograms/users/flosch/Rheometer_Jupyter/FoamQuant/Process.py:199: UserWarning: Only one label was provided to remove_small_objects. Did you mean to use a boolean array?
  image = remove_small_objects(label(image), min_size=Vminobj)

After: Nhol 1
Before: Nhol 19
After: Nhol 1
Cleaned_009: done

Before: Nhol 8
After: Nhol 1
Before: Nhol 26
After: Nhol 1
Cleaned_010: done

## Let's see the result...
# Read the first image of the series
Cleaned = imread(dirsave+namesave+strindex(imrange[1], 3)+'.tiff')
# Show a 3D-cut view of the volume
Cut3D(Cleaned, showcuts=True, showaxes=True)

… we cannot see much like this

Let’s check again the number of objects and holes in the images

for imi in imrange:
    # Read the "non-cleaned" image
    NoCleaned = imread(dirread+nameread+strindex(imi, 3)+'.tiff')
    # regprops of obj and holes
    regions_obj=skimage.measure.regionprops(skimage.measure.label(NoCleaned))
    regions_holes=skimage.measure.regionprops(skimage.measure.label(NoCleaned<1))
    # number of obj and holes
    print(nameread+strindex(imi, 3),'Number of objects:',len(regions_obj), 'Number of holes:',len(regions_holes))

PhaseSegmented_001 Number of objects: 9 Number of holes: 20
PhaseSegmented_002 Number of objects: 6 Number of holes: 21
PhaseSegmented_003 Number of objects: 9 Number of holes: 26
PhaseSegmented_004 Number of objects: 8 Number of holes: 28
PhaseSegmented_005 Number of objects: 3 Number of holes: 22
PhaseSegmented_006 Number of objects: 4 Number of holes: 41
PhaseSegmented_007 Number of objects: 8 Number of holes: 19
PhaseSegmented_008 Number of objects: 6 Number of holes: 23
PhaseSegmented_009 Number of objects: 1 Number of holes: 19
PhaseSegmented_010 Number of objects: 8 Number of holes: 26

for imi in imrange:
    # Read the "cleaned" image
    Cleaned = imread(dirsave+namesave+strindex(imi, 3)+'.tiff')
    # regprops of obj and holes
    regions_obj=skimage.measure.regionprops(skimage.measure.label(Cleaned))
    regions_holes=skimage.measure.regionprops(skimage.measure.label(Cleaned<1))
    # number of obj and holes
    print(namesave+strindex(imi, 3),'Number of objects:',len(regions_obj), 'Number of holes:',len(regions_holes))

Cleaned_001 Number of objects: 1 Number of holes: 1
Cleaned_002 Number of objects: 1 Number of holes: 1
Cleaned_003 Number of objects: 1 Number of holes: 1
Cleaned_004 Number of objects: 1 Number of holes: 1
Cleaned_005 Number of objects: 1 Number of holes: 1
Cleaned_006 Number of objects: 1 Number of holes: 1
Cleaned_007 Number of objects: 1 Number of holes: 1
Cleaned_008 Number of objects: 1 Number of holes: 1
Cleaned_009 Number of objects: 1 Number of holes: 1
Cleaned_010 Number of objects: 1 Number of holes: 1

D) Labelled images

# Read/Save image names and directories
nameread = 'Cleaned_'
namesave = 'BubbleSeg_'
dirread = ProcessPipeline[2]+'/'
dirsave = ProcessPipeline[3]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

# Segment the bubbles with default parrameters
# for more parrameters, try help(BubbleSegmentation_Batch)
BubbleSegmentation_Batch(nameread,
                         namesave,
                         dirread,
                         dirsave,
                         imrange,
                         verbose=True,
                         endread='.tiff',
                         endsave='.tiff',
                         n0=3)

Path exist: True
Distance map: done
Seeds distance map: done
Seeds: done
Watershed distance map: done
Watershed: done
BubbleSeg_001: done

Distance map: done
Seeds distance map: done
Seeds: done
Watershed distance map: done
Watershed: done
BubbleSeg_002: done

Distance map: done
Seeds distance map: done
Seeds: done
Watershed distance map: done
Watershed: done
BubbleSeg_003: done

Distance map: done
Seeds distance map: done
Seeds: done
Watershed distance map: done
Watershed: done
BubbleSeg_004: done

Distance map: done
Seeds distance map: done
Seeds: done
Watershed distance map: done
Watershed: done
BubbleSeg_005: done

Distance map: done
Seeds distance map: done
Seeds: done
Watershed distance map: done
Watershed: done
BubbleSeg_006: done

Distance map: done
Seeds distance map: done
Seeds: done
Watershed distance map: done
Watershed: done
BubbleSeg_007: done

Distance map: done
Seeds distance map: done
Seeds: done
Watershed distance map: done
Watershed: done
BubbleSeg_008: done

Distance map: done
Seeds distance map: done
Seeds: done
Watershed distance map: done
Watershed: done
BubbleSeg_009: done

Distance map: done
Seeds distance map: done
Seeds: done
Watershed distance map: done
Watershed: done
BubbleSeg_010: done

# Create a random colormap
rcmap = RandomCmap(5000)

Number of labels: 5000

Let’s see the result…

# Read the first image of the series
Lab = imread(dirsave+namesave+strindex(imrange[0], 3)+'.tiff')
# Show a 3D-cut view of the volume
Cut3D(Lab,
      showcuts=True,
      showaxes=True,
      cmap=rcmap,
      interpolation='nearest',
      figblocksize=7,           # tune this parrameter if you wish
      zcut=11,                  # tune this parrameter if you wish
      ycut=False,               # tune this parrameter if you wish
      xcut=False)               # tune this parrameter if you wish

E) Visualize the result in Parraview

# Create a .json random colormap that can be used in ParaView
json_rand_dictionary(Ncolors=5000, namecmap='random_cmap.json', dirsave = dirsave, first_color_black=True)

Download your ‘random_cmap.json’ and vizualize your bubble-segmented image in Paraview

Part 2: Quantification

# Create the quantification folders
QuantFolders = ['Q1_LiquidFraction','Q2_RegProps','Q3_Tracking','Q4_MergedTracking']

for Qi in QuantFolders:
    if os.path.exists(Qi):
        print('path already exist:',Qi)
    else:
        print('Created:',Qi)
        os.mkdir(Qi)

path already exist: Q1_LiquidFraction
path already exist: Q2_RegProps
path already exist: Q3_Tracking
path already exist: Q4_MergedTracking

A) Liquid fraction

# Read/Save names and directories
nameread = 'Cleaned_'
namesave = 'LFGlob_'
dirread = ProcessPipeline[2]+'/'
dirsave = QuantFolders[0]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

1) Whole images liquid fraction

# Get the whole images liquid fraction
# (volume percentage of liquid)
LiqFrac_Batch(nameread,
              namesave,
              dirread,
              dirsave,
              imrange,
              TypeGrid='Global',
              verbose=10,
              structured=False)

Path exist: True
Liquid fraction image 1: done
crop:None
LiqFrac:0.14812106324011087
LFGlob_001: done
Liquid fraction image 2: done
crop:None
LiqFrac:0.14536507936507936
LFGlob_002: done
Liquid fraction image 3: done
crop:None
LiqFrac:0.14993783068783068
LFGlob_003: done
Liquid fraction image 4: done
crop:None
LiqFrac:0.14911722096245905
LFGlob_004: done
Liquid fraction image 5: done
crop:None
LiqFrac:0.14338044847568657
LFGlob_005: done
Liquid fraction image 6: done
crop:None
LiqFrac:0.14822921390778535
LFGlob_006: done
Liquid fraction image 7: done
crop:None
LiqFrac:0.14275831443688586
LFGlob_007: done
Liquid fraction image 8: done
crop:None
LiqFrac:0.14482527084908037
LFGlob_008: done
Liquid fraction image 9: done
crop:None
LiqFrac:0.1425321869488536
LFGlob_009: done
Liquid fraction image 10: done
crop:None
LiqFrac:0.14275472411186696
LFGlob_010: done

## Let's see the result...
# Read the liquid fraction of the first image of the series
with open(dirsave + namesave + '001' + '.pkl','rb') as f:
    LF = pkl.load(f)['lf']

print('Whole image liquid fraction:',LF,'%')

Whole image liquid fraction: 0.14812106324011087 %

2) Stuctured liquid fraction in Cartesian subvolumes

# Read/Save image names and directories
nameread = 'Cleaned_'
namesave = 'LFCartesMesh_'
dirread = ProcessPipeline[2]+'/'
dirsave = QuantFolders[0]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

# Get liquid fraction in cartesian subvolumes
# (volume percentage of liquid in each subvolumes)
LiqFrac_Batch(nameread,
              namesave,
              dirread,
              dirsave,
              imrange,
              TypeGrid='CartesMesh',
              Nz=10,        # tune this parrameter if you wish
              Ny=10,        # tune this parrameter if you wish
              Nx=10,        # tune this parrameter if you wish
              verbose=1,
              structured=True)

Path exist: True
LFCartesMesh_001: done
LFCartesMesh_002: done
LFCartesMesh_003: done
LFCartesMesh_004: done
LFCartesMesh_005: done
LFCartesMesh_006: done
LFCartesMesh_007: done
LFCartesMesh_008: done
LFCartesMesh_009: done
LFCartesMesh_010: done

# Read the cartesian grid of liquid fraction
LF=[]
for imi in imrange:
    imifordir = strindex(imi, n0=3)
    with open(dirsave + namesave + imifordir + '.pkl','rb') as f:
        LF.append(pkl.load(f)['lf'])
LF=np.mean(LF,0)

# If "structured=True", the liquid fraction is saved as a 3D mesh by LiqFrac_Batch
# Such as for 3D images, we can reuse Cut3D or Proj3D to vizualise the liquid-fraction meshed volume
fig,ax,neg = Cut3D(LF,
                   vmin=0.11,    # tune this parrameter if you wish
                   vmax=0.22,   # tune this parrameter if you wish
                   cmap='jet',
                   printminmax=True,
                   returnfig=True,
                   figblocksize=7)
# Colorbars
fig.colorbar(neg[1], label=r'$\phi_l$')
fig.colorbar(neg[1], label=r'$\phi_l$')
fig.colorbar(neg[1], label=r'$\phi_l$')

vmin = 0.11 vmax = 0.22
MIN: 0.0943104 MAX: 0.220256
Min0: 0.1130944 Min0: 0.1777792
Min1: 0.101056 Max1 0.20806399999999997
Min2: 0.101056 Max2: 0.21072639999999998

<matplotlib.colorbar.Colorbar at 0x2ba3e379ea90>

# The Proj3D function is similar to Cut3D, but show an average along the 3 directions
fig,ax,neg = Proj3D(LF,
                    vmin=0.11,    # tune this parrameter if you wish
                    vmax=0.22,   # tune this parrameter if you wish
                    cmap='jet',
                    printminmax=True,
                    returnfig=True,
                    figblocksize=7)

# Colorbars
fig.colorbar(neg[0], label=r'$\phi_l$')
fig.colorbar(neg[1], label=r'$\phi_l$')
fig.colorbar(neg[2], label=r'$\phi_l$')

vmin = 0.11 vmax = 0.22
MIN: 0.0943104 MAX: 0.220256
Min0: 0.132816 Min0: 0.16055424
Min1: 0.12508662153846156 Max1 0.1802792426035503
Min2: 0.12507100591715978 Max2: 0.18537196307692308

/gpfs/offline1/staff/tomograms/users/flosch/Rheometer_Jupyter/FoamQuant/Figure.py:140: UserWarning: This figure was using constrained_layout==True, but that is incompatible with subplots_adjust and or tight_layout: setting constrained_layout==False.
  plt.tight_layout()

<matplotlib.colorbar.Colorbar at 0x2ba3db8f08b0>

3) Unstuctured liquid fraction in Cartesian subvolumes

If “structured=False”, the liquid fraction is saved as a 1D array by LiqFrac_Batch

This may be practical if one want to plot the liquid fraction as a function of other parrameters such as the cartesian/cylindrical/spherical coordinates or the bubble deformation, for example.

# structured = False
LiqFrac_Batch(nameread,
              namesave,
              dirread,
              dirsave,
              imrange,
              TypeGrid='CartesMesh',
              Nz=10,      # tune this parrameter if you wish
              Ny=5,       # tune this parrameter if you wish
              Nx=5,       # tune this parrameter if you wish
              verbose=1,
              structured=False)

Path exist: True
LFCartesMesh_001: done
LFCartesMesh_002: done
LFCartesMesh_003: done
LFCartesMesh_004: done
LFCartesMesh_005: done
LFCartesMesh_006: done
LFCartesMesh_007: done
LFCartesMesh_008: done
LFCartesMesh_009: done
LFCartesMesh_010: done

# We can plot the liquid fraction as a function of the z coordinate for the first image
with open(dirsave + namesave + '001' + '.pkl','rb') as f:
    pack = pkl.load(f)
lf = pack['lf']
z = pack['zgrid']

fig, ax = plt.subplots(1,1, figsize = (10, 10))
plt.plot(z, lf,'ko', alpha=0.5)
plt.xlabel(r'$z$ (vox)', fontsize=30)
plt.ylabel(r'$\phi_l$', fontsize=30)
plt.ylim((0,0.36))
plt.grid(True)

B) Individual bubble properties

We are going to extract the individual bubble volume, radius, sphericity, moment of inertial, strain tensor, etc.

# Read/Save names and directories
nameread = 'BubbleSeg_'
namesave = 'Props_'
dirread = ProcessPipeline[3]+'/'
dirsave = QuantFolders[1]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

Get some properties in the given field of view (field=[zmin,zmax,ymin,ymax,xmin,xmax])

• Label and centroid coodinate: ‘lab’,‘z’,‘y’,‘x’

• Volume, equivalent radius, area, sphericity: ‘vol’,‘rad’,‘area’,‘sph’

• Volume from ellipsoid fit: ‘volfit’

• Ellipsoid three semi-axis and eigenvectors: ‘S1’,‘S2’,‘S3’,‘e1z’,‘e1y’,‘e1x’,‘e2z’,‘e2y’,‘e2x’,‘e3z’,‘e3y’,‘e3x’,

• Internal strain components: ‘U1’,‘U2’,‘U3’

• Internal strain von Mises invariant: ‘U’

• Oblate (-1) or prolate (1) ellipsoid:‘type’

# Region properties
RegionProp_Batch(nameread,
                 namesave,
                 dirread,
                 dirsave,
                 imrange,
                 verbose=True,
                 field=[40,220,40,220,40,220], # tune this parrameter if you wish
                 endread='.tiff',
                 endsave='.tsv')

Path exist: True
Props_001: done
Props_002: done
Props_003: done
Props_004: done
Props_005: done
Props_006: done
Props_007: done
Props_008: done
Props_009: done
Props_010: done

# Read the regionprop files
properties = Read_RegionProp(namesave, dirsave, imrange)

# histogram of some extracted properties
prop=['vol','rad','area','sph']
xlab=[r'$Volume$ (vox)',r'$Radius$ (vox)',r'$Area$ (vox)',r'$Sphericity$']

fig, ax = plt.subplots(1,4, figsize = (7*4, 7), constrained_layout=True)
for i in range(4):
    H=ax[i].hist(properties[prop[i]], bins=50)
    ax[i].set_xlabel(xlab[i], fontsize=20)
    ax[i].set_ylabel(r'Frequency', fontsize=20)
    ax[i].grid(True)
    ax[i].set_yscale('log')   # tune this parrameter if you wish

More properties can be extracted from individual images, such as coordination (number of neighbours and contact topology), and individual contact area and orientation.

The SPAM package is great for extracting these properties! Have a look here if you wish to know more: https://ttk.gricad-pages.univ-grenoble-alpes.fr/spam/

C) Velocity field

1) Tracking table

We are going to track the centroid of each bubble from one image to the next

# Read/Save image names and directories
nameread = 'Props_'
namesave = 'Tracking_'
dirread = QuantFolders[1]+'/'
dirsave = QuantFolders[2]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

Track the bubbles from one image to the next: - keep the bubbles candidate if their centroid is inside the box searchbox=[zmin,zmax,ymin,ymax,xmin,xmax] - keep the bubbles candidate if the volume from one image to the next is not changing more than Volpercent (the next volume should be between (1-Volpercent)V and (1+Volpercent)V) - select the bubble with the closest distance

# Tracking
LLostlab, LLostX, LLostY, LLostZ = LabelTracking_Batch(nameread,
                                                       namesave,
                                                       dirread,
                                                       dirsave,
                                                       imrange,
                                                       verbose=False,
                                                       endread='.tsv',
                                                       endsave='.tsv',
                                                       n0=3,
                                                       searchbox=[-10,10,-10,10,-10,10],   # tune this parrameter if you wish
                                                       Volpercent=0.05)              # tune this parrameter if you wish

# Lost tracking: N percentage
# Lost tracking: 2 1.2738853503184715 %

Path exist: True
Lost tracking: 13 2.579365079365079 %
Lost tracking: 9 1.7786561264822136 %
Lost tracking: 7 1.36986301369863 %
Lost tracking: 16 3.11284046692607 %
Lost tracking: 11 2.1825396825396823 %
Lost tracking: 12 2.3483365949119372 %
Lost tracking: 9 1.761252446183953 %
Lost tracking: 11 2.156862745098039 %
Lost tracking: 7 1.3806706114398422 %

# Read the tracking files
tracking = Read_LabelTracking(namesave, dirsave, imrange, verbose=True)

Tracking_001_002 : done
Tracking_002_003 : done
Tracking_003_004 : done
Tracking_004_005 : done
Tracking_005_006 : done
Tracking_006_007 : done
Tracking_007_008 : done
Tracking_008_009 : done
Tracking_009_010 : done

# Convert -1 in np.nan
# and create Coord and v arrays (non-structured coordinate and velocity arrays)
Listx1 = tracking['x1']
Listy1 = tracking['y1']
Listz1 = tracking['z1']

Listx2 = tracking['x2']
Listy2 = tracking['y2']
Listz2 = tracking['z2']

Coord = []; v=[]
for vali in range(len(Listx1)):
    for i in range(len(Listx1[vali])):
        if Listx1[vali][i]==-1:
            Listx1[vali][i]=np.nan
        if Listy1[vali][i]==-1:
            Listy1[vali][i]=np.nan
        if Listz1[vali][i]==-1:
            Listz1[vali][i]=np.nan
        if Listx2[vali][i]==-1:
            Listx2[vali][i]=np.nan
        if Listy2[vali][i]==-1:
            Listy2[vali][i]=np.nan
        if Listz2[vali][i]==-1:
            Listz2[vali][i]=np.nan
        Coord.append([Listz1[vali][i],
                      Listy1[vali][i],
                      Listx1[vali][i]])
        v.append([Listz2[vali][i]-Listz1[vali][i],
                  Listy2[vali][i]-Listy1[vali][i],
                  Listx2[vali][i]-Listx1[vali][i]])

Coord=np.asarray(Coord)
v=np.asarray(v)

# Create a linear colormap
lincmap = LinCmap(vmin=0, vmax=len(LLostX), first_color="lime", last_color="k")

fig, ax = plt.subplots(1,3, figsize = (45, 15), constrained_layout=True)
# Plot of tracked bubbles (green)
for i in range(len(Listx1)):
    ax[0].plot([Listx1[i],Listx2[i]],[Listy1[i],Listy2[i]],'.-', color=lincmap.to_rgba(i))
    ax[1].plot([Listx1[i],Listx2[i]],[Listz1[i],Listz2[i]],'.-', color=lincmap.to_rgba(i))
    ax[2].plot([Listy1[i],Listy2[i]],[Listz1[i],Listz2[i]],'.-', color=lincmap.to_rgba(i))
# Plot of the lost position (red)
for i in range(len(LLostlab)):
    ax[0].plot(LLostX[i],LLostY[i],'ro')
    ax[1].plot(LLostX[i],LLostZ[i],'ro')
    ax[2].plot(LLostY[i],LLostZ[i],'ro')
# Axes
ax[0].set_xlabel('x', fontsize=20); ax[0].set_ylabel('y', fontsize=20)
ax[1].set_xlabel('x', fontsize=20); ax[1].set_ylabel('z', fontsize=20)
ax[2].set_xlabel('y', fontsize=20); ax[2].set_ylabel('z', fontsize=20)

Text(0, 0.5, 'z')

# in 3D!
# Plot of tracked bubbles (green)
ax = plt.figure(figsize = (10, 10)).add_subplot(projection='3d')
for i in range(len(Listx1)):
    for j in range(len(Listx1[i])):
        ax.plot([Listx1[i][j],Listx2[i][j]],
                [Listy1[i][j],Listy2[i][j]],
                [Listz1[i][j],Listz2[i][j]],
                '.-', color=lincmap.to_rgba(i))
# Plot of the lost position (red)
for i in range(len(LLostlab)):
    ax.plot(LLostX[i],LLostY[i],LLostZ[i],'ro')
# Axes
ax.set(xlim=(40, 220), ylim=(40, 220), zlim=(40, 220),
       xlabel='x', ylabel='y', zlabel='z')

[(40.0, 220.0),
 (40.0, 220.0),
 (40.0, 220.0),
 Text(0.5, 0, 'x'),
 Text(0.5, 0, 'y'),
 Text(0.5, 0, 'z')]

2) Structured averaged flow field

If structured=True, Grid_Vavg function can be used to obtained cartesian 3D structured average flow field (here respectively along z,y,x)

# structured = True
# tune Range and N parrameters if you wish
Lgrid_z, Coordavg_z,Vavg_z,Vstd_z, Count_z = Grid_Vavg(Coord, v, Range=[40,220,40,220,40,220], N=[1,10,10], NanFill=True, verbose=False, structured=True)
Lgrid_y, Coordavg_y,Vavg_y,Vstd_y, Count_y = Grid_Vavg(Coord, v, Range=[40,220,40,220,40,220], N=[10,1,10], NanFill=True, verbose=False, structured=True)
Lgrid_x, Coordavg_x,Vavg_x,Vstd_x, Count_x = Grid_Vavg(Coord, v, Range=[40,220,40,220,40,220], N=[10,10,1], NanFill=True, verbose=False, structured=True)

vmin=0   # tune this parrameter if you wish
vmax=100  # tune this parrameter if you wish

fig, ax = plt.subplots(1,3, figsize = (30, 10), constrained_layout=True)
# Colormesh plot of the number of bubble inside the averaging box
neg1=ax[0].pcolormesh(Lgrid_z[2],Lgrid_z[1], Count_z[0,:,:], cmap = 'jet', shading='nearest', vmin=vmin,vmax=vmax)
neg2=ax[1].pcolormesh(Lgrid_y[2],Lgrid_y[0], Count_y[:,0,:], cmap = 'jet', shading='nearest', vmin=vmin,vmax=vmax)
neg3=ax[2].pcolormesh(Lgrid_x[1],Lgrid_x[0], Count_x[:,:,0], cmap = 'jet', shading='nearest', vmin=vmin,vmax=vmax)
# Averaged flow field
ax[0].quiver(Lgrid_z[2],Lgrid_z[1],Vavg_z[0,:,:,2],Vavg_z[0,:,:,1], pivot='mid')
ax[1].quiver(Lgrid_y[2],Lgrid_y[0],Vavg_y[:,0,:,2],Vavg_y[:,0,:,0], pivot='mid')
ax[2].quiver(Lgrid_x[1],Lgrid_x[0],Vavg_x[:,:,0,1],Vavg_x[:,:,0,0], pivot='mid')
# Axes
ax[0].set_xlabel('x'); ax[0].set_ylabel('y')
ax[1].set_xlabel('x'); ax[1].set_ylabel('z')
ax[2].set_xlabel('y'); ax[2].set_ylabel('z')
# Colorbars
fig.colorbar(neg1, label='Count')
fig.colorbar(neg2, label='Count')
fig.colorbar(neg3, label='Count')

<matplotlib.colorbar.Colorbar at 0x2ba3fba8ee80>

3) Unstructured averaged flow field

If structured=False, Grid_Vavg function can be used to obtained the same average, but unstructured. Here it is convenient for example plotting the velocity along x as a function of the position z

# structured = False
Lgrid, Coordavg,Pavg,Pstd, Count = Grid_Pavg(Coord,
                                             v[:,0],
                                             Range=[40,220,40,220,40,220],  # tune this parrameter if you wish
                                             N=[10,1,1],                    # tune this parrameter if you wish
                                             NanFill=True,
                                             verbose=False,
                                             structured=False)

# Do an average over (time,z,y):
ax = plt.figure(figsize = (10, 10))
for vali in range(len(Listx1)):
    plt.plot(Listz1[vali], Listz2[vali]-Listz1[vali],'k.', alpha=0.1)
# Averaged velocity along x as a function of the position z
plt.errorbar(Lgrid[0], Pavg, yerr=Pstd, capsize=5, color='k')
plt.xlabel('z (vox)', fontsize=20); plt.ylabel('Vy (vox/image)', fontsize=20)
plt.grid(True)

4) From Cartesian to Cylindrical coordinates

# Convert Cartesian to Cylindrical: (z,y,x) -> (r,theta,z)
CoordCyl, VCyl = Cartesian2Cylindrical_Vector(Coord, v, CoordAxis = [0,125,125], CylAxisZ = [1,0,0],CylAxisY = [0,1,0],CylAxisX = [0,0,1])

# Do an average over (time,r,z):
Lgrid, Coordavg,Pavg,Pstd, Count = Grid_Pavg(CoordCyl,
                                             VCyl[:,2],
                                             Range=[0,110,-np.pi,np.pi,40,220],  # tune this parrameter if you wish
                                             N=[1,10,1],                         # tune this parrameter if you wish
                                             NanFill=True,
                                             verbose=False,
                                             structured=False)

# Velocity along z as a function of the azimuthal angular position theta
ax = plt.figure(figsize = (10, 10))
plt.plot(CoordCyl[:,1], VCyl[:,2],'k.', alpha=0.1)
plt.errorbar(Lgrid[1], Pavg, yerr=Pstd, capsize=5, color='k')
plt.xlabel(r'$\theta$ (rad)', fontsize=20); plt.ylabel('Vz (vox/image)', fontsize=20)
plt.grid(True)

D) Combine the tracking tables

Get the individual bubble flow path by combining the subsequent image tracking tables.

# Combine the subsequent trackings over the whole series
combined = Combine_Tracking(namesave, dirsave, imrange, verbose=False, endread='.tsv', n0=3)

# Convert lost bubbles data into np.nan
for axis in ['z','y','x']:
    for i in range(len(combined[axis])):
        for j in range(len(combined[axis][i])):
            if combined[axis][i][j]==-1:
                combined[axis][i][j]=np.nan

# Show the individual paths by a random color
fig, ax = plt.subplots(1,3, figsize = (45, 15), constrained_layout=True)
for i in range(len(combined['x'])):
    ax[0].plot(combined['x'][i],combined['y'][i])
    ax[1].plot(combined['x'][i],combined['z'][i])
    ax[2].plot(combined['y'][i],combined['z'][i])
# Axes
ax[0].set_xlabel('x'); ax[0].set_ylabel('y')
ax[1].set_xlabel('x'); ax[1].set_ylabel('z')
ax[2].set_xlabel('y'); ax[2].set_ylabel('z')

Text(0, 0.5, 'z')

# in 3D!
# Show the individual paths by a random color
ax = plt.figure(figsize = (10, 10)).add_subplot(projection='3d')
for i in range(len(combined['x'])):
    plt.plot(combined['x'][i],combined['y'][i],combined['z'][i], '.-')
# Axes
ax.set(xlim=(40, 220), ylim=(40, 220), zlim=(40, 220),
       xlabel='x', ylabel='y', zlabel='z')

[(40.0, 220.0),
 (40.0, 220.0),
 (40.0, 220.0),
 Text(0.5, 0, 'x'),
 Text(0.5, 0, 'y'),
 Text(0.5, 0, 'z')]

# in 3D!
# Show the non-tracked paths in red
ax = plt.figure(figsize = (10, 10)).add_subplot(projection='3d')
for i in range(len(Listx1)):
    for j in range(len(Listx1[i])):
        ax.plot([Listx1[i][j],Listx2[i][j]],
                [Listy1[i][j],Listy2[i][j]],
                [Listz1[i][j],Listz2[i][j]],
                'r,-')
# Show the tracked paths in green
for i in range(len(combined['x'])):
    ax.plot(combined['x'][i],combined['y'][i],combined['z'][i], 'g.-')

namesave = 'MergedTracking'
dirsave = QuantFolders[3]+'/'

# You can then open it in Paraview! and use the filter "table to points", then "glyph" and "theshold" for example.
Save_Tracking(combined, namesave, dirsave)

For more information see https://foamquant.readthedocs.io/

Liquid foam flow: advance analysis

Not yet available.

In-situ bread baking

Jupyter example:

Welcome!

• In this jupyter notebook, as an example, we are going to analyse 10 croped images of a bread sample, having \(Kondis\) with extra gluten, baked by convection heating and imaged in the middle of the bread (series \(Extra-Convection-Middle\)).

• 5 images at the very beginning of the baking (indexes 1 to 5), and 5 images at the very end of the baking (indexes 6 to 10).

# Import Libraries
import numpy as np
import matplotlib.pyplot as plt
from tifffile import imread, imsave
import skimage.measure
import pickle as pkl
import os
import spam

# Import FoamQuant library
from FoamQuant import *

# Set matplotlib default font size
plt.rc('font', size=20)

Part 1: Image processing

# Create the processing pipeline
ProcessPipeline = ['P1_Raw',
                   'P2_PhaseSegmented',
                   'P3_Cleaned',
                   'P4_PoreSegmented',
                   'P5_PoreNoEdge',
                   'P6_WallThickness',
                   'P7_Contact']

for Pi in ProcessPipeline:
    if  os.path.exists(Pi):
        print('path already exist:',Pi)
    else:
        print('Created:',Pi)
        os.mkdir(Pi)

path already exist: P1_Raw
path already exist: P2_PhaseSegmented
path already exist: P3_Cleaned
path already exist: P4_PoreSegmented
path already exist: P5_PoreNoEdge
path already exist: P6_WallThickness
path already exist: P7_Contact

A) The raw images

# Read/Save image names and directories
nameread = 'Raw_'
namesave = 'PhaseSegmented_'
dirread = ProcessPipeline[0]+'/'
dirsave = ProcessPipeline[1]+'/'
# Images indexes:
# -> 1 to 5 are the 5 first images of the series
# -> 6 to 10 are the 5 last images of the series
imrange = [1,2,3,4,5,6,7,8,9,10]

# Read the first and last image of the series
RawFirst = imread(dirread+nameread+strindex(imrange[0], 3)+'.tif')
RawLast = imread(dirread+nameread+strindex(imrange[-1], 3)+'.tif')

# Show a 3D-cut view of the two volume
Cut3D(RawFirst,
      showcuts=True,
      showaxes=True,
      figblocksize=7,
      zcut=50,       # tune this parrameter if you wish
      ycut=False,    # tune this parrameter if you wish
      xcut=False,    # tune this parrameter if you wish
      cmap='bone')

Cut3D(RawLast,
      showcuts=True,
      showaxes=True,
      figblocksize=7,
      zcut=50,       # tune this parrameter if you wish
      ycut=False,    # tune this parrameter if you wish
      xcut=False,    # tune this parrameter if you wish
      cmap='bone')

Example of raw images in the beginning of baking

Example of raw images in the end of baking

B) Phase segmentation

# Otsu simple threshold phase segmentation of the whole series
th = PhaseSegmentation_Batch(nameread,
                             namesave,
                             dirread,
                             dirsave,
                             imrange,
                             method='ostu_global',
                             returnOtsu=True,
                             verbose=True,
                             n0=3,
                             endread='.tif',
                             endsave='.tif')

PhaseSegmented_ 1: done

PhaseSegmented_ 2: done

PhaseSegmented_ 3: done

PhaseSegmented_ 4: done

PhaseSegmented_ 5: done

PhaseSegmented_ 6: done

PhaseSegmented_ 7: done

PhaseSegmented_ 8: done

PhaseSegmented_ 9: done

PhaseSegmented_ 10: done

# Otsu thresholds used for the segmentation
print('Otsu thresholds:',th)

Otsu thresholds: [93, 93, 94, 94, 94, 94, 94, 94, 94, 94]

Let’s see the result…

# Read the first and last image of the series
SegFirst = imread(dirsave+namesave+strindex(imrange[0], 3)+'.tif')
SegLast = imread(dirsave+namesave+strindex(imrange[-1], 3)+'.tif')

# Let's see the result for the first image
zcut=50       # tune this parrameter if you wish
ycut=False    # tune this parrameter if you wish
xcut=False    # tune this parrameter if you wish
cmap='bone'   # tune this parrameter if you wish: e.g. 'bone'

# Show a 3D-cut view of the volume
Cut3D(SegFirst, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)                  # Phase segmented image
Cut3D((SegFirst>0)*RawFirst, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)     # Phase segmented image * Raw image
Cut3D((1-SegFirst)*RawFirst, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)     # (1-Phase segmented image) * Raw image

Check the phase segmentation of the first raw image (in the beginning of baking)

# Let's see the result for the last image
zcut=50       # tune this parrameter if you wish
ycut=False    # tune this parrameter if you wish
xcut=False    # tune this parrameter if you wish
cmap='bone'    # tune this parrameter if you wish: e.g. 'bone'

# Show a 3D-cut view of the volume
Cut3D(SegLast, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)             # Phase segmented image
Cut3D((SegLast>0)*RawLast, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)     # Phase segmented image * Raw image
Cut3D((1-SegLast)*RawLast, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)     # (1-Phase segmented image) * Raw image

Check the phase segmentation of the last raw image (in the end of baking)

C) Remove small holes & regions

Due to artefacts the phase segmented images can have speckles, e.g. small pore regions (1, white) and holes (0, black). In our case (bread images), this is problematic for the pore segmentations. One may oversegment the pores due to small holes flowting in the middle of the pore regions.

• Since in the bread images the solid phase consist of a unique region (0, black), one can filter all the small holes.

• However for the small pore regions (1, white), we cannot be sure. They may be speckle artefacts or actual pores in the bread. Therefore they are not going to be filtered. All the pore having a volumes below 3x3x3=27 voxels should be disgarded because of the this resolution limitation.

# Read/Save image names and directories
nameread = 'PhaseSegmented_'
namesave = 'Cleaned_'
dirread = ProcessPipeline[1]+'/'
dirsave = ProcessPipeline[2]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

Remove all holes with: - Vhole < Chole * max(Vhole)

Since in bread images, the solid phase consist of unique large regions, Chole can be strict (large thresholds). All the other smaller regions are often due to imaging artefacts.

# remove holes and objects
RemoveSpeckleBin_Batch(nameread,
                       namesave,
                       dirread,
                       dirsave,
                       imrange,
                       verbose=True,
                       RemoveObjects=False,
                       RemoveHoles=True,
                       BinClosing=False,
                       endread='.tif',
                       endsave='.tif',
                       n0=3,
                       Chole=0.1) # tune this parrameter if you wish

Before: Nobj 6104
After: Nobj 1
Before: Nhol 983
After: Nhol 1
First image (vox): maxObj 4987372 maxHole 22006312
Thresholds (vox): thrObj 2493686 thrHole 2200631

Before: Nhol 1316
After: Nhol 1
Cleaned_001: done

Before: Nhol 1283
After: Nhol 1
Cleaned_002: done

Before: Nhol 1212
After: Nhol 1
Cleaned_003: done

Before: Nhol 1384
After: Nhol 1
Cleaned_004: done

Before: Nhol 1418
After: Nhol 1
Cleaned_005: done

Before: Nhol 1868
After: Nhol 1
Cleaned_006: done

Before: Nhol 1945
After: Nhol 1
Cleaned_007: done

Before: Nhol 1878
After: Nhol 1
Cleaned_008: done

Before: Nhol 1928
After: Nhol 1
Cleaned_009: done

Before: Nhol 2064
After: Nhol 1
Cleaned_010: done

# Read the first image of the series
Seg = imread(dirread+nameread+strindex(imrange[0], 3)+'.tif')
Cleaned = imread(dirsave+namesave+strindex(imrange[0], 3)+'.tif')

zcut=50       # tune this parrameter if you wish
ycut=False    # tune this parrameter if you wish
xcut=False    # tune this parrameter if you wish
cmap='bone'    # tune this parrameter if you wish: e.g. 'bone'

# Show a 3D-cut view of the volume
Cut3D(Seg, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)
Cut3D(Cleaned, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)
Cut3D(Cleaned-Seg, showcuts=True, showaxes=True, figblocksize=7,zcut=zcut,ycut=ycut,xcut=xcut, cmap=cmap)

Check the small holes removal in the first image of the series (in the beginning of baking)

D) Labelled images

# Read/Save image names and directories
nameread = 'Cleaned_'
namesave = 'PoreSeg_'
dirread = ProcessPipeline[2]+'/'
dirsave = ProcessPipeline[3]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

# Segment the pores with ITK watershed (ITK=True) implemented in SPAM (works well for polydispersed and elongated pores)
# if you wish, you can also try the default watershed segmentation, and play with the segmentation parrameters (ITK=False)
BubbleSegmentation_Batch(nameread,
                         namesave,
                         dirread,
                         dirsave,
                         imrange,
                         ITK=True,
                         ITKLevel=1,
                         verbose=True,
                         endread='.tif',
                         endsave='.tif',
                         n0=3)

Path exist: True
PoreSeg_001: done

PoreSeg_002: done

PoreSeg_003: done

PoreSeg_004: done

PoreSeg_005: done

PoreSeg_006: done

PoreSeg_007: done

PoreSeg_008: done

PoreSeg_009: done

PoreSeg_010: done

# Create a random colormap to distinguish the pores
rcmap = RandomCmap(5000)

Number of labels: 5000

Let’s see the result…

# Read the first image of the series
LabFirst = imread(dirsave+namesave+strindex(imrange[0], 3)+'.tif')
LabLast = imread(dirsave+namesave+strindex(imrange[-1], 3)+'.tif')
# Show a 3D-cut view of the volume
Cut3D(LabFirst,
      showcuts=True,
      showaxes=True,
      cmap=rcmap,
      interpolation='nearest',
      figblocksize=7,           # tune this parrameter if you wish
      zcut=50,                  # tune this parrameter if you wish
      ycut=False,               # tune this parrameter if you wish
      xcut=False)               # tune this parrameter if you wish

Cut3D(LabLast,
      showcuts=True,
      showaxes=True,
      cmap=rcmap,
      interpolation='nearest',
      figblocksize=7,           # tune this parrameter if you wish
      zcut=50,                  # tune this parrameter if you wish
      ycut=False,               # tune this parrameter if you wish
      xcut=False)               # tune this parrameter if you wish

Check the pores segmentation in the first image of the series (in the beginning of baking)

Check the pores segmentation in the last image of the series (in the end of baking)

-> To visualize the segmentation result in Paraview

Download your ‘random_cmap.json’ and vizualize your pore-segmented image in Paraview

# Create a .json random colormap that can be used in ParaView
json_rand_dictionary(Ncolors=5000, namecmap='random_cmap.json', dirsave = dirsave, first_color_black=True)

E) Remove the pores at the edges of the image

The pores on the edge of the images (or cut by the mask) are irrelevant for measuring the individual pore properties. To obtain clean statistics, all the pores touching the edges of the image are removed.

# Read/Save image names and directories
nameread = 'PoreSeg_'
namesave = 'PoreNoEdge_'
dirread = ProcessPipeline[3]+'/'
dirsave = ProcessPipeline[4]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

# Remove the pores at the end of the image (default)
# for more parrameters, try help(BubbleSegmentation_Batch)
RemoveEdgeBubble_Batch(nameread,
                       namesave,
                       dirread,
                       dirsave,
                       imrange,
                       verbose=True,
                       endread='.tif',
                       endsave='.tif',
                       n0=3)

Path exist: True
PoreNoEdge_001: done

PoreNoEdge_002: done

PoreNoEdge_003: done

PoreNoEdge_004: done

PoreNoEdge_005: done

PoreNoEdge_006: done

PoreNoEdge_007: done

PoreNoEdge_008: done

PoreNoEdge_009: done

PoreNoEdge_010: done

Let’s see the result…

# Read the first image of the series
NoedgeFirst = imread(dirsave+namesave+strindex(imrange[0], 3)+'.tif')
NoedgeLast = imread(dirsave+namesave+strindex(imrange[-1], 3)+'.tif')
# Show a 3D-cut view of the volume
Cut3D(NoedgeFirst,
      showcuts=True,
      showaxes=True,
      cmap=rcmap,
      interpolation='nearest',
      figblocksize=7,           # tune this parrameter if you wish
      zcut=50,                  # tune this parrameter if you wish
      ycut=False,               # tune this parrameter if you wish
      xcut=False)               # tune this parrameter if you wish

Cut3D(NoedgeLast,
      showcuts=True,
      showaxes=True,
      cmap=rcmap,
      interpolation='nearest',
      figblocksize=7,           # tune this parrameter if you wish
      zcut=50,                  # tune this parrameter if you wish
      ycut=False,               # tune this parrameter if you wish
      xcut=False)               # tune this parrameter if you wish

Check the pores removal on the edges in the first image of the series (in the beginning of baking)

Check the pores removal on the edges in the last image of the series (in the end of baking)

G) Contact images

We can then extract contacts data from the pore-segmented images (both with edges and no-edges are required). The GetContacts function save the coordination table, coordination images and contact table, batchwise.

# Read/Save image names and directories
nameread = 'PoreSeg_'
nameread_noedge = 'PoreNoEdge_'
namesave = 'Contact_'
dirread = ProcessPipeline[3]+'/'
dirread_noedge = ProcessPipeline[4]+'/'
dirsave = ProcessPipeline[6]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

GetContacts_Batch(nameread, nameread_noedge, namesave, dirread, dirread_noedge, dirsave, imrange,
                  verbose=False,
                  endread='.tif',
                  endread_noedge='.tif',
                  endsave='.tif',
                  n0=3,
                  save='all',
                  maximumCoordinationNumber=20)

Path exist: True

# Read the first image of the series
Lab = imread(dirsave+namesave+strindex(imrange[0], 3)+'.tif')
# Show a 3D-cut view of the volume
Cut3D(Lab,
      showcuts=True,
      showaxes=True,
      cmap=rcmap,
      interpolation='nearest',
      figblocksize=7,           # tune this parrameter if you wish
      zcut=False,                  # tune this parrameter if you wish
      ycut=False,               # tune this parrameter if you wish
      xcut=False)               # tune this parrameter if you wish

# Read the last image of the series
Lab = imread(dirsave+namesave+strindex(imrange[-1], 3)+'.tif')
# Show a 3D-cut view of the volume
Cut3D(Lab,
      showcuts=True,
      showaxes=True,
      cmap=rcmap,
      interpolation='nearest',
      figblocksize=7,           # tune this parrameter if you wish
      zcut=False,                  # tune this parrameter if you wish
      ycut=False,               # tune this parrameter if you wish
      xcut=False)               # tune this parrameter if you wish

Check the contacts labelling in the first image of the series (in the beginning of baking)

Check the contacts labelling in the last image of the series (in the end of baking)

F) Fast Local-wall thickness

Finally, the fast local wall thickness tool developped by Dahl, V. A. and Dahl A. B. (Git-link February 2023: https://github.com/vedranaa/local-thickness.git) can be used to determine the distribution of thicknesses and mean wall thicknes over the whole image.

# If you wish to import the pakage please do as follow
import localthickness as lt

# Read/Save image names and directories
nameread = 'Cleaned_'
namesave = 'WallThickness_'
dirread = ProcessPipeline[2]+'/'
dirsave = ProcessPipeline[5]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

# The localthickness function is used in FastLocalThickness_Batch for batchwise analysis!
FastLocalThickness_Batch(nameread, namesave, dirread, dirsave, imrange,
                         verbose=True,
                         endread='.tif',
                         endsave='.tif',
                         n0=3,
                         WalThickness=True,
                         Separation=True,
                         scale=1)

Path exist: True
WallThickness_001: done

WallThickness_002: done

WallThickness_003: done

WallThickness_004: done

WallThickness_005: done

WallThickness_006: done

WallThickness_007: done

WallThickness_008: done

WallThickness_009: done

WallThickness_010: done

Let’s see the result…

# Read the first and last Local Wall thickness images of the series
WTFirst = imread(dirsave+namesave+strindex(imrange[0], 3)+'_WT.tif')
WTLast = imread(dirsave+namesave+strindex(imrange[-1], 3)+'_WT.tif')

# Show a 3D-cut view of the volume
Cut3D(WTFirst,
      showcuts=True,
      showaxes=True,
      cmap='jet',
      figblocksize=7,           # tune this parrameter if you wish
      zcut=50,                  # tune this parrameter if you wish
      ycut=False,               # tune this parrameter if you wish
      xcut=False)               # tune this parrameter if you wish

Cut3D(WTLast,
      showcuts=True,
      showaxes=True,
      cmap='jet',
      figblocksize=7,           # tune this parrameter if you wish
      zcut=50,                  # tune this parrameter if you wish
      ycut=False,               # tune this parrameter if you wish
      xcut=False)               # tune this parrameter if you wish

Check the first local-wall thickness image of the series (in the beginning of baking)

# Read the first and last Separation images of the series
SEPFirst = imread(dirsave+namesave+strindex(imrange[0], 3)+'_SEP.tif')
SEPLast = imread(dirsave+namesave+strindex(imrange[-1], 3)+'_SEP.tif')
# Show a 3D-cut view of the volume
Cut3D(SEPFirst,
      showcuts=True,
      showaxes=True,
      cmap='jet',
      figblocksize=7,           # tune this parrameter if you wish
      zcut=50,                  # tune this parrameter if you wish
      ycut=False,               # tune this parrameter if you wish
      xcut=False)               # tune this parrameter if you wish

Cut3D(SEPLast,
      showcuts=True,
      showaxes=True,
      cmap='jet',
      figblocksize=7,           # tune this parrameter if you wish
      zcut=50,                  # tune this parrameter if you wish
      ycut=False,               # tune this parrameter if you wish
      xcut=False)               # tune this parrameter if you wish

Check the first separation image of the series (in the beginning of baking)

Check the last separation image of the series (in the end of baking)

Part 2: Quantification

In this second part, we are going to reuse the processed images for quantifying bread properties:

• Porosity \(\phi\)

• Volume \(V\), Shape eigenvalues (\(S_1\),\(S_2\),\(S_3\)), Elongation \(E\)

• Coordination \(Z\)

• Local Wall Thickness \(h_w\)

# Create the quantification folders
QuantFolders = ['Q1_Porosity','Q2_RegProps','Q3_WallThickness']

for Qi in QuantFolders:
    if os.path.exists(Qi):
        print('path already exist:',Qi)
    else:
        print('Created:',Qi)
        os.mkdir(Qi)

path already exist: Q1_Porosity
path already exist: Q2_RegProps
path already exist: Q3_WallThickness

A) Porosity \(\phi\)

# Read/Save names and directories
nameread = 'Cleaned_'
namesave = 'Porosity_'
dirread = ProcessPipeline[2]+'/'
dirsave = QuantFolders[0]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

# Get the whole images liquid fraction
# (volume percentage of liquid)
LiqFrac_Batch(nameread,
              namesave,
              dirread,
              dirsave,
              imrange,
              TypeGrid='Global',
              verbose=1,
              structured=False)

Path exist: True
Porosity_001: done
Porosity_002: done
Porosity_003: done
Porosity_004: done
Porosity_005: done
Porosity_006: done
Porosity_007: done
Porosity_008: done
Porosity_009: done
Porosity_010: done

## Let's see the result...
# Read the liquid fraction of the first image of the series
LPorosity=[]
for imi in imrange:
    with open(dirsave+namesave+strindex(imi,3)+'.pkl','rb') as f:
        SF = pkl.load(f)['lf']
    print(imi, 'Whole image porosity:',round(1-SF,3),'%')
    LPorosity.append(1-SF)

1 Whole image porosity: 0.322 %
2 Whole image porosity: 0.321 %
3 Whole image porosity: 0.323 %
4 Whole image porosity: 0.321 %
5 Whole image porosity: 0.321 %
6 Whole image porosity: 0.578 %
7 Whole image porosity: 0.581 %
8 Whole image porosity: 0.583 %
9 Whole image porosity: 0.585 %
10 Whole image porosity: 0.587 %

# Read/Save names and directories
nameread = 'Cleaned_'
namesave = 'CartesPorosity_'
dirread = ProcessPipeline[2]+'/'
dirsave = QuantFolders[0]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

# structured = False
LiqFrac_Batch(nameread,
              namesave,
              dirread,
              dirsave,
              imrange,
              TypeGrid='CartesMesh',
              Nz=10,      # tune this parrameter if you wish
              Ny=1,       # tune this parrameter if you wish
              Nx=1,       # tune this parrameter if you wish
              verbose=1,
              structured=False)

Path exist: True
CartesPorosity_001: done
CartesPorosity_002: done
CartesPorosity_003: done
CartesPorosity_004: done
CartesPorosity_005: done
CartesPorosity_006: done
CartesPorosity_007: done
CartesPorosity_008: done
CartesPorosity_009: done
CartesPorosity_010: done

# We can plot the liquid fraction as a function of the z coordinate for the first image
fig, ax = plt.subplots(1,1, figsize = (10, 10))
for imi in imrange:
    with open(dirsave+namesave+strindex(imi,3)+'.pkl','rb') as f:
        pack = pkl.load(f)
    lf = pack['lf']
    z = pack['zgrid']
    plt.plot(z, (1-np.asarray(lf))*100,'o-', alpha=0.5, label='Image {0}'.format(imi))
plt.xlabel(r'$z$ (vox)', fontsize=30)
plt.ylabel(r'$Porosity$ ($\%$)', fontsize=30)
plt.grid(True)
plt.legend(fontsize=17)

<matplotlib.legend.Legend at 0x2b63b8aee2b0>

Porosity (in percent) as a function of the vertical position z (in voxels)

B) Individual pores properties

We are going to extract the individual pore volume, radius, sphericity, moment of inertial, strain tensor, etc.

These are the properties we are mainly interested in: - Volume \(V\) - Shape eigenvalues (\(S_1\),\(S_2\),\(S_3\)) - Elongation \(E\)

# Read/Save names and directories
nameread = 'PoreSeg_'
namesave = 'Props_'
dirread = ProcessPipeline[3]+'/'
dirsave = QuantFolders[1]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

Get some properties in the given field of view (field=[zmin,zmax,ymin,ymax,xmin,xmax])

• Label and centroid coodinate: ‘lab’,‘z’,‘y’,‘x’

• Volume, equivalent radius, area, sphericity: ‘vol’,‘rad’,‘area’,‘sph’

• Volume from ellipsoid fit: ‘volfit’

• Ellipsoid three semi-axis and eigenvectors: ‘S1’,‘S2’,‘S3’,‘e1z’,‘e1y’,‘e1x’,‘e2z’,‘e2y’,‘e2x’,‘e3z’,‘e3y’,‘e3x’,

• Internal strain components: ‘U1’,‘U2’,‘U3’

• Internal strain von Mises invariant: ‘U’

• Oblate (-1) or prolate (1) ellipsoid:‘type’

# Region properties
RegionProp_Batch(nameread,
                 namesave,
                 dirread,
                 dirsave,
                 imrange,
                 verbose=True,
                 field=[40,220,40,220,40,220], # tune this parrameter if you wish
                 endread='.tif',
                 endsave='.tsv')

Path exist: True
Props_001: done
Props_002: done
Props_003: done
Props_004: done
Props_005: done
Props_006: done
Props_007: done
Props_008: done
Props_009: done
Props_010: done

# Read the regionprop files
properties_Beg = Read_RegionProp(namesave, dirsave, imrange[0:5])
# Read the regionprop files
properties_End = Read_RegionProp(namesave, dirsave, imrange[5:9])

# histogram of some extracted properties
prop=['vol','S3','S2','S1','U','type']
xlab=[r'$Volume$ (vox)',r'$Shape$ (vox)',r'$Elongation$',r'<- Prolate / Oblate ->']

fig, ax = plt.subplots(1,4, figsize = (7*4, 7), constrained_layout=True)

# Volume distribution of the pores
bins=np.power(10,np.linspace(np.log10(1),np.log10(1e5),100))
H=ax[0].hist(properties_Beg[prop[0]], bins=bins,alpha=0.5)
ax[0].set_xlabel(xlab[0], fontsize=20)
ax[0].set_ylabel(r'Frequency', fontsize=20)
ax[0].grid(True)
ax[0].set_yscale('log')   # tune this parrameter if you wish
ax[0].set_xscale('log')   # tune this parrameter if you wish

# Shape eigenvalues distribution of the pores
bins=np.power(10,np.linspace(np.log10(1),np.log10(60),100))
for i in range(3):
    H=ax[1].hist(properties_Beg[prop[i+1]],
                 bins=bins,
                 alpha=0.5,
                 label=prop[i+1])
ax[1].set_xlabel(xlab[1], fontsize=20)
ax[1].set_ylabel(r'Frequency', fontsize=20)
ax[1].grid(True)
ax[1].set_yscale('log')   # tune this parrameter if you wish
ax[1].set_xscale('log')   # tune this parrameter if you wish
ax[1].legend(fontsize=15)

# Elongation distribution of the pores: shape eig max / eig min
bins=np.power(10,np.linspace(np.log10(1),np.log10(8),100))
H=ax[2].hist(properties_Beg[prop[1]] / properties_Beg[prop[3]],
             bins=bins,
             alpha=0.5)
ax[2].set_xlabel(xlab[2], fontsize=20)
ax[2].set_ylabel(r'Frequency', fontsize=20)
ax[2].grid(True)
ax[2].set_yscale('log')   # tune this parrameter if you wish
ax[2].set_xscale('log')   # tune this parrameter if you wish

# Deviation from a spherical shape: how oblate or prolate are the pore?
bins=np.linspace(-2,2,100)
H=ax[3].hist(properties_Beg[prop[4]] * properties_Beg[prop[5]],
             bins=bins,
             alpha=0.5)
ax[3].set_xlabel(xlab[3], fontsize=20)
ax[3].set_ylabel(r'Frequency', fontsize=20)
ax[3].grid(True)
ax[3].set_yscale('log')   # tune this parrameter if you wish

Individual pore properties histograms

C) Coordination \(Z\)

# Read/Save image names and directories
namesave = 'Contact_'
dirsave = ProcessPipeline[6]+'/'
# Images indexes
imrange = [1,2,3,4,5,6,7,8,9,10]

# Read the first image of the series
Lab = imread(dirsave+'Coordination_'+strindex(imrange[0], 3)+'.tif')
# Show a 3D-cut view of the volume
fig,ax,neg = Cut3D(Lab,
                   cmap='seismic',
                   interpolation='nearest',
                   figblocksize=7,
                   returnfig=True,
                   vmin=0,
                   vmax=20)

fig.colorbar(neg[2], label=r'$Z$')

# Read the last image of the series
Lab = imread(dirsave+'Coordination_'+strindex(imrange[-1], 3)+'.tif')
# Show a 3D-cut view of the volume
fig,ax,neg = Cut3D(Lab,
                   cmap='seismic',
                   interpolation='nearest',
                   figblocksize=7,
                   returnfig=True,
                   vmin=0,
                   vmax=20)

fig.colorbar(neg[2], label=r'$Z$')

vmin = 0 vmax = 20

/gpfs/offline1/staff/tomograms/users/flosch/Old/PSI_2021_Bread/FoamQuant/Figure.py:78: UserWarning: This figure was using constrained_layout==True, but that is incompatible with subplots_adjust and or tight_layout: setting constrained_layout==False.
  plt.tight_layout()

vmin = 0 vmax = 20

<matplotlib.colorbar.Colorbar at 0x2b643aef7340>

Coordination image of the first image of the series (in the beginning of baking)

Coordination image of the last image of the series (in the end of baking)

# Read contact table
TableFirst = ReadContactTable(namesave+'table_', dirsave, imrange[:5], verbose=False)
TableLast = ReadContactTable(namesave+'table_', dirsave, imrange[5:], verbose=False)

# remove the values at the edge
LZnoedge=[]
for Table in [TableFirst, TableLast]:
    Znoedge = []; coordnoedge = []
    for t in range(len(Table)):
        table = Table[t]
        for i in range(len(table['Z'])):
            if table['lab_noedge'][i]>0:
                Znoedge.append(table['Z'][i])
                coordnoedge.append([table['z'][i],table['y'][i],table['x'][i]])
    coordnoedge = np.asarray(coordnoedge)
    LZnoedge.append(Znoedge)

# Coordination histogram before and after removing the edge pores
fig, ax = plt.subplots(1,1, figsize = (7, 7), constrained_layout=True)
H=ax.hist(LZnoedge, bins=21, label=['Beginning','End'], alpha=0.8)
ax.set_xlabel(r'Coordination $Z$', fontsize=20)
ax.set_ylabel(r'Frequency', fontsize=20)
ax.grid(True)
ax.set_yscale('log')
ax.legend(fontsize=15)

<matplotlib.legend.Legend at 0x2b643afaccd0>

Coordination histogram

C) Local Wall Thickness \(h_w\)

def Remove0Array(image):
    import numpy as np
    ZZ,YY,XX = np.shape(image)
    Array0 = np.reshape(image,(ZZ*YY*XX))
    Array=[]
    for i in range(len(Array0)):
        if Array0[i]>0:
            Array.append(Array0[i])
    return Array

# Load first image as list
WTFirstList = Remove0Array(WTFirst)
WTFirstLast = Remove0Array(WTLast)

fig, ax = plt.subplots(1,1, figsize = (7, 7), constrained_layout=True)
bins=np.power(10,np.linspace(np.log10(1),np.log10(13),10))

H=plt.hist(WTFirstList, alpha=0.5, bins=bins, label='beginning')
H=plt.hist(WTFirstLast, alpha=0.5, bins=bins, label='end')
ax.set_xlabel(r'$h_w$ (vox)', fontsize=20)
ax.set_ylabel(r'Frequency', fontsize=20)
ax.grid(True)
ax.set_yscale('log')   # tune this parrameter if you wish
ax.set_xscale('log')   # tune this parrameter if you wish

Local-wall thickness histogram

You have now completed the jupyter example for analysis on bread! We hope this has been useful to you!

For more information on the tools, the references or contacts, have a look on https://foamquant.readthedocs.io

Note

Contacts: Florian Schott and Rajmund Mokso