Mathematical Morphology¶

In [1]:
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
    return false;
}

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

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

Setting some global variables

In [3]:
# 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 [4]:
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 [5]:
localVD = gl.VectorDouble(ntotal)

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

In [6]:
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 [7]:
vmin = -1
vmax = +1
image2 = gl.BImage(nxy)
tab = grid.getColumn("Simu")
image = gl.morpho_double2image(nxy,tab,vmin,vmax)
In [8]:
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 [9]:
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 [10]:
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 [11]:
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 [12]:
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 [13]:
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 [14]:
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 [15]:
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 [16]:
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 [17]:
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 [18]:
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 [19]:
grid.setLocator("Reference",gl.ELoc.Z)
dum = grid.morpho(gl.EMorpho.NEGATION,0.5,1.5,verbose=False)

if flagGraphic:
    ax = grid.plot()
In [20]:
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 [21]:
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 [22]:
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 [23]:
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 [24]:
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 [25]:
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 [26]:
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 [27]:
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 [28]:
image
Out[28]:
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 [29]:
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 [ ]: