---
title: "Mathematical Morphology"
author: "D. Renard"
date: "17 January 2017"
output: 
  rmarkdown::html_vignette:
    fig_width: 5
    fig_height: 4
vignette: >
  %\VignetteIndexEntry{Mathematical Morphology}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

This paper is dedicated to a demonstration of the procedures of Mathematical Morphology implemented in RGeostats.

```{r Loading_library,include=FALSE}
knitr::opts_chunk$set(echo = TRUE,fig.align='center')
library(RGeostats)
```

These features will all be demonstrated on a regular 2-D grid (even if they are also available in 3-D). All the graphics will be displayed with a constant aspect ratio equal to 1.
```{r Grid creation}
constant.define("asp",1)
grid = db.create(nx=c(500,300))
```

We first generate a non-conditional simulation which will serve as the starting information stored on the grid. This simulation is performed using the Turning Band simulation Method, using a Model composed of an isotropic cubic structure. Its range is equal to 30 (to be compared with the grid extension 500 by 300)
```{r Initial Simulation}
model = model.create(vartype="Cubic",range=40)
grid = simtub(,grid,model,nbtuba=1000,seed=3213)
grid = db.rename(grid,grid$natt,"Simulation")
plot(grid)
```

The first task is to truncate the previous Gaussian variable (which has a symmetric spread of values around 0). The result is a binary image with 1 when the initial variable was positive, and o otherwise.
```{r Binary Image}
grid = morpho(grid,vmin=0,oper="thresh")
grid = db.rename(grid,grid$natt,"Binary")
plot(grid)
```

Staring from this binary image, we perform an erosion of this image, using a structuring element which corresponds to a cross with an extension of 1 pixel.
```{r Erosion}
radius = 1
option = 0
grid = db.locate(grid,"Simulation","z")
grid = morpho(grid,vmin=0,oper="erode",option=option,radius=radius)
grid = db.rename(grid,grid$natt,"Erosion")
plot(grid)
```

For comparison sake, we can perform a dilation instead, using the same structuring element.
```{r Dilation}
grid = db.locate(grid,"Simulation","z")
grid = morpho(grid,vmin=0,oper="dilate",option=option,radius=radius)
grid = db.rename(grid,grid$natt,"Dilation")
plot(grid)
```

An interesting feature is obtained by comparing the results of the Dilation to the one of the Erosion: it highlight the edge.
```{r Edge}
grid = db.add(grid,Dilation-Erosion)
grid = db.rename(grid,grid$natt,"Edge")
plot(grid)
```

The opening operation consists in an erosion followed by a dilation.
```{r Opening}
grid = db.locate(grid,"Simulation","z")
grid = morpho(grid,vmin=0,oper="open",option=option,radius=radius)
grid = db.rename(grid,grid$natt,"Opening")
plot(grid)
```

while the closing is the inverse sequence, i.e. a dilation followed by an erosion.
```{r Closing}
grid = db.locate(grid,"Simulation","z")
grid = morpho(grid,vmin=0,oper="close",option=option,radius=radius)
grid = db.rename(grid,grid$natt,"Closing")
plot(grid)
```

All previous operations can be run again changing the structuring element from cross to block. The next figure demonstrates the erosion with this new option.
```{r Erosion Block}
radius = 1
option = 1
grid = db.locate(grid,"Simulation","z")
grid = morpho(grid,vmin=0,oper="erode",option=option,radius=radius)
grid = db.rename(grid,grid$natt,"Erosion.Block")
plot(grid)
```


Each one of the previous applications can be performed several times in sequence. The following figure shows the erosion (with the block structuring element) iterated 5 times.
```{r Erosion iterated}
niter=5
grid = db.locate(grid,"Simulation","z")
grid = morpho(grid,vmin=0,oper="erode",option=option,radius=radius,
              niter=niter)
grid = db.rename(grid,grid$natt,"Erosion_iterated")
plot(grid,title=paste0("Erosion (x",niter,")"))
```

Starting from the last binary image (resulting from the Erosion iterated), we now wish to distinguish the different connected components.
```{r Connect components}
grid = db.locate(grid,"Erosion_iterated","z")
grid = morpho(grid,vmin=0.5,oper="ccsize")
grid = db.rename(grid,grid$natt,"Connected components")
plot(grid)
```

The last feature is to calculate the distance to the edge (of the set where initial values were 1). The black color is used to represent pixels where initial values were 0.
```{r Distance to Edge}
grid = db.locate(grid,"Erosion_iterated","z")
grid = morpho(grid,vmin=0.5,oper="dist")
grid = db.rename(grid,grid$natt,"Distance to Edge")
plot(grid)
```