Mathematical Morphology¶

This file is meant to demonstrate the use of gstlearn for Morphological Operations

In [1]:
import numpy as np
import pandas as pd
import sys
import os
import gstlearn as gl
import gstlearn.plot as gp
import gstlearn.document as gdoc
import matplotlib.pyplot as plt

gdoc.setNoScroll()

Setting some global variables

In [2]:
# Set the Global Options
verbose = True
flagGraphic = True
gp.setDefaultGeographic(dims=[5,5])

# 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)
In [3]:
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.

In [4]:
localVD = gl.VectorDouble(ntotal)

Generating an initial square grid covering a 1 by 1 surface (100 meshes along each direction).

In [5]:
grid = gl.DbGrid.create(nxy, [dx,dy])
model = gl.Model.createFromParam(gl.ECov.CUBIC, 0.1, 1.)
gl.simtub(None, grid, model)
if verbose:
    grid.display()
if flagGraphic:
    ax = grid.plot("Simu")
    ax.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.

In [6]:
vmin = -1
vmax = +1
image2 = gl.BImage(nxy)
tab = grid.getColumn("Simu")
image = gl.morpho_double2image(nxy,tab,vmin,vmax)
In [7]:
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.

In [8]:
if flagGraphic:
    gl.morpho_image2double(image, 0, 1., 0., localVD)
    iuid = grid.addColumns(localVD,"Initial Image",gl.ELoc.Z)
    ax = grid.plot()

The next step interchanges grain and pore

In [9]:
gl.morpho_negation(image, image2)

if flagGraphic:
    gl.morpho_image2double(image2, 0, 1., 0., localVD)
    iuid = grid.addColumns(localVD,"Negative Image",gl.ELoc.Z)
    ax = grid.plot()

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)

In [10]:
gl.morpho_erosion(0, [1,1], image, image2)

if flagGraphic:
    gl.morpho_image2double(image2, 0, 1., 0., localVD)
    iuid = grid.addColumns(localVD,"Erosion - Cross",gl.ELoc.Z)
    ax = grid.plot()

We check the result of ersosion when choosing the Block structuring element

In [11]:
gl.morpho_erosion(1, [1,1], image, image2)

if flagGraphic:
    gl.morpho_image2double(image2, 0, 1., 0., localVD)
    iuid = grid.addColumns(localVD,"Erosion - Block",gl.ELoc.Z)
    ax = grid.plot()

We now perform the dilation of the Initial image (only the Cross structuring element will be used in the next paragraphs)

In [12]:
gl.morpho_dilation(0, [1,1], image, image2)

if flagGraphic:
    gl.morpho_image2double(image2, 0, 1., 0., localVD)
    iuid = grid.addColumns(localVD,"Dilation - Cross",gl.ELoc.Z)
    ax = grid.plot()

Combining the elementary operations (Erosion and Dilation), we can perform directly an opening

In [13]:
gl.morpho_opening(0, [1,1], image, image2)

if flagGraphic:
    gl.morpho_image2double(image2, 0, 1., 0., localVD)
    iuid = grid.addColumns(localVD,"Opening - Cross",gl.ELoc.Z)
    ax = grid.plot()

And the closing

In [14]:
gl.morpho_closing(0, [1,1], image, image2)

if flagGraphic:
    gl.morpho_image2double(image2, 0, 1., 0., localVD)
    iuid = grid.addColumns(localVD,"Closing - Cross",gl.ELoc.Z)
    ax = grid.plot()

Connected components¶

Starting from the Initial image (stored in the Data Base as Reference), we now wish to determine the connected components

In [15]:
gl.morpho_negation(image, image)

if flagGraphic:
    gl.morpho_image2double(image, 0, 1., 0., localVD)
    iuid = grid.addColumns(localVD,"Reference",gl.ELoc.Z)
    ax = grid.plot()
    ax.decoration(title="Starting image used for Connected Components")

Compute the connected components

In [16]:
cc = gl.morpho_labelling(0, 0, image, np.nan)

if flagGraphic:
    iuid = grid.addColumns(cc,"Connect Components by Rank",gl.ELoc.Z)
    ax = grid.plot()
In [17]:
cc = gl.morpho_labelling(0, 1, image, np.nan)

if flagGraphic:
    iuid = grid.addColumns(cc,"Connect Components by Volume",gl.ELoc.Z)
    ax = grid.plot()

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.

In [18]:
grid.setLocator("Reference",gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.NEGATION,0.5,1.5,verbose=False)

if flagGraphic:
    ax = grid.plot()
In [19]:
grid.setLocator("Reference",gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.CC,0.5, 1.5)

if flagGraphic:
    ax = grid.plot()
    ax.decoration(title="Connected Components by Rank")
In [20]:
grid.setLocator("Reference",gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.CCSIZE,0.5,1.5)

if flagGraphic:
    ax = grid.plot()
    ax.decoration(title="Connected Components by Volume")

Calculation of distance to the edge of the grain

In [21]:
grid.setLocator("Simu",gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.DISTANCE,radius=[1,1],verbose=False)

if flagGraphic:
    ax = grid.plot("*DISTANCE")
    ax.decoration(title="Distance to the Edge")

Calculation of the 2-D angle of the gradient

In [22]:
grid.setLocator("Simu",gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.ANGLE,radius=[1,1],verbose=False)

if flagGraphic:
    ax = grid.plot("*ANGLE")
    ax.decoration(title="Angle of Gradient")

Calculation of Gradient components

In [23]:
grid.setLocator("Simu",gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.GRADIENT)

if flagGraphic:
    fig, axs = plt.subplots(1, 2, figsize=(10,5))
    axs[0].raster(grid,"Morpho.Simu.GRADIENT.1")
    ax.decoration(title="Gradient along X")
    
    axs[1].raster(grid,"Morpho.Simu.GRADIENT.2")
    ax.decoration(title="Gradient along Y")

Smoothing the input image

In [24]:
grid.setLocator("Simu",gl.ELoc.Z,0,True)
neighI = gl.NeighImage([3,3])
dum = grid.smooth(neighI)

if flagGraphic:
    ax = grid.plot("Smooth*")
    ax.decoration(title="Smoothed Image")

Testing the Bitmap Image printout¶

In [25]:
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.)
gl.simtub(None, grid, model)
if flagGraphic:
    ax = grid.plot()
    ax.decoration(title="Initial Data Set")
In [26]:
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., localVD)
iuid = grid.addColumns(localVD,"Initial Image",gl.ELoc.Z)
ax = grid.plot()
In [27]:
image
Out[27]:
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
In [28]:
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 ..+.....
In [ ]: