Mathematical Morphology¶
This file is meant to demonstrate the use of gstlearn for Morphological Operations
import numpy as np
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
gdoc.setNoScroll()
Setting some global variables
# Set the Global Options
verbose = True
flagGraphic = True
# Define the Space Dimension
ndim = 2
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim)
# Set the Seed for the Random Number generator
gl.law_set_random_seed(5584)
nx = ny = 100
nxy = [nx, ny]
ntotal = nx * ny
dx = dy = 0.01
In several usages, we will need a VectorDouble dimensionned to the total number of pixels. This is created next.
localVD = gl.VectorDouble(ntotal)
Generating an initial square grid covering a 1 by 1 surface (100 meshes along each direction).
grid = gl.DbGrid.create(nxy, [dx, dy])
model = gl.Model.createFromParam(gl.ECov.CUBIC, 0.1, 1.0)
gl.simtub(None, grid, model)
if verbose:
grid.display()
if flagGraphic:
gp.plot(grid, "Simu")
gp.decoration(title="Initial Data Set")
Data Base Grid Characteristics ============================== Data Base Summary ----------------- File is organized as a regular grid Space dimension = 2 Number of Columns = 4 Total number of samples = 10000 Grid characteristics: --------------------- Origin : 0.000 0.000 Mesh : 0.010 0.010 Number : 100 100 Variables --------- Column = 0 - Name = rank - Locator = NA Column = 1 - Name = x1 - Locator = x1 Column = 2 - Name = x2 - Locator = x2 Column = 3 - Name = Simu - Locator = z1
Basic operations¶
We retreive the newly simulated variable (called Simu) from the grid Db into a local Vector (called tab). This vector is then transformed by thresholding and loaded into an image object (called image). This object is very efficient as each pixel is stored into a single bit. A secondary Image object (called image2) is created and will be used in subsequent diadic operations.
vmin = -1
vmax = +1
image2 = gl.BImage(nxy)
tab = grid.getColumn("Simu")
image = gl.morpho_double2image(nxy, tab, vmin, vmax)
volume = gl.morpho_count(image)
print("Grain Volume =", volume, " /", ntotal, "(pixels)\n")
Grain Volume = 6920 / 10000 (pixels)
For visualization (and i the current version), we must first convert the image into a vector and load it into a grid.
if flagGraphic:
gl.morpho_image2double(image, 0, 1.0, 0.0, localVD)
iuid = grid.addColumns(localVD, "Initial Image", gl.ELoc.Z)
res = gp.plot(grid)
The next step interchanges grain and pore
gl.morpho_negation(image, image2)
if flagGraphic:
gl.morpho_image2double(image2, 0, 1.0, 0.0, localVD)
iuid = grid.addColumns(localVD, "Negative Image", gl.ELoc.Z)
res = gp.plot(grid)
Basic Mophological Image transformations¶
We start with the initial image and perform an erosion. The second argument defines the type of structuring element: either Cross (0) or Block (1)
gl.morpho_erosion(0, [1, 1], image, image2)
if flagGraphic:
gl.morpho_image2double(image2, 0, 1.0, 0.0, localVD)
iuid = grid.addColumns(localVD, "Erosion - Cross", gl.ELoc.Z)
res = gp.plot(grid)
We check the result of ersosion when choosing the Block structuring element
gl.morpho_erosion(1, [1, 1], image, image2)
if flagGraphic:
gl.morpho_image2double(image2, 0, 1.0, 0.0, localVD)
iuid = grid.addColumns(localVD, "Erosion - Block", gl.ELoc.Z)
res = gp.plot(grid)
We now perform the dilation of the Initial image (only the Cross structuring element will be used in the next paragraphs)
gl.morpho_dilation(0, [1, 1], image, image2)
if flagGraphic:
gl.morpho_image2double(image2, 0, 1.0, 0.0, localVD)
iuid = grid.addColumns(localVD, "Dilation - Cross", gl.ELoc.Z)
res = gp.plot(grid)
Combining the elementary operations (Erosion and Dilation), we can perform directly an opening
gl.morpho_opening(0, [1, 1], image, image2)
if flagGraphic:
gl.morpho_image2double(image2, 0, 1.0, 0.0, localVD)
iuid = grid.addColumns(localVD, "Opening - Cross", gl.ELoc.Z)
res = gp.plot(grid)
And the closing
gl.morpho_closing(0, [1, 1], image, image2)
if flagGraphic:
gl.morpho_image2double(image2, 0, 1.0, 0.0, localVD)
iuid = grid.addColumns(localVD, "Closing - Cross", gl.ELoc.Z)
res = gp.plot(grid)
Connected components¶
Starting from the Initial image (stored in the Data Base as Reference), we now wish to determine the connected components
gl.morpho_negation(image, image)
if flagGraphic:
gl.morpho_image2double(image, 0, 1.0, 0.0, localVD)
iuid = grid.addColumns(localVD, "Reference", gl.ELoc.Z)
gp.plot(grid)
gp.decoration(title="Starting image used for Connected Components")
Compute the connected components
cc = gl.morpho_labelling(0, 0, image, np.nan)
if flagGraphic:
iuid = grid.addColumns(cc, "Connect Components by Rank", gl.ELoc.Z)
res = gp.plot(grid)
cc = gl.morpho_labelling(0, 1, image, np.nan)
if flagGraphic:
iuid = grid.addColumns(cc, "Connect Components by Volume", gl.ELoc.Z)
res = gp.plot(grid)
Some shortcuts¶
The procedure is made easier for basic morphological operations, using directly the method morpho. This method operates from a variable stored in a Db (organized as a Grid) and returns the result in the same Db.
grid.setLocator("Reference", gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.NEGATION, 0.5, 1.5, verbose=False)
if flagGraphic:
res = gp.plot(grid)
grid.setLocator("Reference", gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.CC, 0.5, 1.5)
if flagGraphic:
gp.plot(grid)
gp.decoration(title="Connected Components by Rank")
grid.setLocator("Reference", gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.CCSIZE, 0.5, 1.5)
if flagGraphic:
gp.plot(grid)
gp.decoration(title="Connected Components by Volume")
Calculation of distance to the edge of the grain
grid.setLocator("Simu", gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.DISTANCE, radius=[1, 1], verbose=False)
if flagGraphic:
gp.plot(grid, "*DISTANCE")
gp.decoration(title="Distance to the Edge")
Calculation of the 2-D angle of the gradient
grid.setLocator("Simu", gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.ANGLE, radius=[1, 1], verbose=False)
if flagGraphic:
gp.plot(grid, "*ANGLE")
gp.decoration(title="Angle of Gradient")
Calculation of Gradient components
grid.setLocator("Simu", gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.GRADIENT)
if flagGraphic:
fig, axs = gp.init(1, 2, figsize=(10, 5))
axs[0].raster(grid, "Morpho.Simu.GRADIENT.1")
axs[0].decoration(title="Gradient along X")
axs[1].raster(grid, "Morpho.Simu.GRADIENT.2")
axs[1].decoration(title="Gradient along Y")
Smoothing the input image
grid.setLocator("Simu", gl.ELoc.Z, 0, True)
neighI = gl.NeighImage([3, 3])
dum = grid.smooth(neighI)
if flagGraphic:
gp.plot(grid, "Smooth*")
gp.decoration(title="Smoothed Image")
Testing the Bitmap Image printout¶
nx = 10
ny = 12
nxy = [nx, ny]
localVD = gl.VectorDouble(nx * ny)
grid = gl.DbGrid.create(nxy)
model = gl.Model.createFromParam(gl.ECov.CUBIC, 6, 1.0)
gl.simtub(None, grid, model)
if flagGraphic:
gp.plot(grid)
gp.decoration(title="Initial Data Set")
vmin = -1
vmax = +1
image2 = gl.BImage(nxy)
tab = grid.getColumn("Simu")
image = gl.morpho_double2image(nxy, tab, vmin, vmax)
gl.morpho_image2double(image, 0, 1.0, 0.0, localVD)
iuid = grid.addColumns(localVD, "Initial Image", gl.ELoc.Z)
res = gp.plot(grid)
image
Array dimension = 2 - Dimension #1 : 10 - Dimension #2 : 12 1234567890 1 1110010011 2 1110011011 3 1110011111 4 1110011111 5 1110001111 6 1111001111 7 1111111110 8 1111111100 9 1110111110 0 0110011110 1 0110000111 2 0011011111
bstrfmt = gl.BImageStringFormat("+", ".", [2, 3])
image.display(bstrfmt)
Array dimension = 2 - Dimension #1 : 10 - Dimension #2 : 12 34567890 4 .++..... 5 .+++.... 6 ..++.... 7 .......+ 8 ......++ 9 .+.....+ 0 .++....+ 1 .++++... 2 ..+.....