R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Interface to the Growing Region algorithm

calcGR

R Documentation

Interface to the Growing Region algorithm

Description

Call the GRalgo function to perform the Growing Region algorithm.

Usage

calcGR(contrast, W, seed, sigma_max, range = c(-Inf, +Inf), range.seed = c(-Inf, +Inf),
         breaks = 100, rescale = FALSE, iter_max = 100, sd.robust = FALSE,
		 keep.lower = FALSE, keep.upper = FALSE, verbose = TRUE,
		 history.sigma = FALSE, history.step = FALSE, history.front = FALSE)

Arguments

`contrast`	the contrast value of each observation. numeric vector. REQUIRED.
`W`	the neighbourhood matrix. dgCMatrix. REQUIRED.
`seed`	the index of the initial seeds or a binary indicator of the initial seeds. positive integer vector or logical vector. REQUIRED.
`sigma_max`	the maximum admissible value for the variability of the group contrast. positive numeric. REQUIRED.
`range`	the range of acceptable contrast values for the growing region group. numeric vector of size 2.
`range.seed`	the range of acceptable contrast values for the seeds. numeric vector of size 2.
`breaks`	the break points or the number of break points to use to categorize the contrast distribution. numeric vector or `postive integer`.
`rescale`	should the contrast be scaled ? logical.
`iter_max`	the maximum number of iterations for the expansion of the growing region. `postive integer`.
`sd.robust`	should the median absolute deviation be used to estimte the variability of the group contrast, or the standard deviation ? logical.
`keep.lower`	should removing observations with high intensity of the region be forbidden ? logical.
`keep.upper`	should removing observations with low intensity of the region be forbidden ? logical.
`verbose`	should the execution of the function be traced ? logical.
`history.sigma`	should the values of sigma be recorded ? logical.
`history.step`	should the number of observations included in the growing region set be recorded ? logical.
`history.front`	should the propagation front of the GR set be recorded ? logical.

Details

FUNCTION:
This implementation of the Growing Region algorithm was been proposed by (Revol et al. 1997).

Value

An list containing :

[[GR]] : the index of the observations in the growing region. integer vector.
[[test.break]] : whether the GR algorithm was interrupted an during execution. logical.
[[iter]] : the number of the last iteration of the algorithm. integer.
[[test.id]] : whether the GR set has stabilised during the last iteration. logical.
[[sigma]] : if history.sigma was set to TRUE, the value of the homogeneity criterion at the begining and the end of each step (in columns) for all steps (in row). numeric matrix.
[[history_GR]] : if history.step was set to TRUE, the step when each GR observation was included in the GR set. integer vector.
[[breaks]] : if history.front was set to TRUE, the values used to categorize the contrast. numeric vector.

References

Chantal Revol and Michel Jourlin. A new minimum varance region growing algorithm for image segmentation. Pattern Recognition Letters, 18(3):249-258,1997.

Examples

## load a MRIaggr object
data(MRIaggr.Pat1_red, package = "MRIaggr")

calcThresholdMRIaggr(MRIaggr.Pat1_red, param = c("TTP_t0","MTT_t0"), threshold = 1:10,
                     name_newparam = c("TTP.th_t0","MTT.th_t0"),
                     update.object = TRUE, overwrite = TRUE)
					 
## display raw parameter
multiplot(MRIaggr.Pat1_red, param = "TTP.th_t0", num = 3, numeric2logical = TRUE,
          index1 = list(coords = "MASK_DWI_t0", outline = TRUE))

## extract raw parameter, coordinates and compute the neighbourhood matrix
carto <- selectContrast(MRIaggr.Pat1_red, num = 3, hemisphere = "lesion",
                         param = c("TTP.th_t0","TTP_t0","MASK_DWI_t0"))
coords <- selectCoords(MRIaggr.Pat1_red, num = 3, hemisphere = "lesion")
W <- calcW(coords, range = sqrt(2))$W

## the seed is taken to be the point with the largest TTP in the lesion mask
indexN <- which(carto$MASK_DWI_t0 == 1)
seed <- indexN[which.max(carto[indexN,"TTP_t0"])] 

## Display step by step the GR algorithm with sigma = 1
for(iter in c(0,1,2,5,10)){
  resGR1 <- calcGR(contrast = carto$TTP.th_t0, W = W, 
                   seed = seed, sigma_max = 1, iter_max = iter, verbose = FALSE)
  
  multiplot(MRIaggr.Pat1_red, param = "TTP.th_t0", num = 3,hemisphere = "lesion", legend = FALSE,
         breaks = seq(0,10,0.1), numeric2logical = TRUE, cex = 2,
		 main = paste("iteration=",iter," - slice ", sep = ""),          
         index1 = list(coords = coords[resGR1$GR,], pch = 20, cex = 1),
         index2 = list(coords = coords[seed,], pch = 20, cex = 1)
  )
}

## Not run: 
## GR with sigma = 2.2
resGR2 <- calcGR(contrast = carto$TTP.th_t0, W = W, 
                 seed = seed, sigma_max = 2.2, iter_max = 50,
                 history.step = TRUE, history.front = TRUE)

## display 
# display the GR over the raw contrast
multiplot(MRIaggr.Pat1_red, param = "TTP.th_t0", num = 3, hemisphere = "lesion", legend = FALSE,
             breaks = seq(0,10,0.1), numeric2logical = TRUE, cex = 2,            
             index1 = list(coords = coords[resGR2$GR,], pch = 20, cex = 1)
)

# display the step of inclusion in GR group for each observation
multiplot(coords[resGR2$GR,],
             resGR2$history.step,breaks = 0:10,
             index1=list(coords = coords[seed,]),
             palette = rainbow(10)
)

# display the front propagation 
multiplot(coords[resGR2$GR,],
             resGR2$Mfront[,7],
             index1 = list(coords = coords[seed,])
)

## End(Not run)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(MRIaggr)
Loading required package: Rcpp
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/MRIaggr/GRalgo-calcGR.Rd_%03d_medium.png", width=480, height=480)
> ### Name: calcGR
> ### Title: Interface to the Growing Region algorithm
> ### Aliases: calcGR
> ### Keywords: functions
> 
> ### ** Examples
> 
> ## load a MRIaggr object
> data(MRIaggr.Pat1_red, package = "MRIaggr")
> 
> calcThresholdMRIaggr(MRIaggr.Pat1_red, param = c("TTP_t0","MTT_t0"), threshold = 1:10,
+                      name_newparam = c("TTP.th_t0","MTT.th_t0"),
+                      update.object = TRUE, overwrite = TRUE)
Step 1 : keep both hemipheres, keep CSF 
Step 2 : thresholding **********
allocContrast[MRIaggr] : Cartographies "TTP.th_t0" "MTT.th_t0" 
                         have been allocated 
> 					 
> ## display raw parameter
> multiplot(MRIaggr.Pat1_red, param = "TTP.th_t0", num = 3, numeric2logical = TRUE,
+           index1 = list(coords = "MASK_DWI_t0", outline = TRUE))
> 
> ## extract raw parameter, coordinates and compute the neighbourhood matrix
> carto <- selectContrast(MRIaggr.Pat1_red, num = 3, hemisphere = "lesion",
+                          param = c("TTP.th_t0","TTP_t0","MASK_DWI_t0"))
> coords <- selectCoords(MRIaggr.Pat1_red, num = 3, hemisphere = "lesion")
> W <- calcW(coords, range = sqrt(2))$W
> 
> ## the seed is taken to be the point with the largest TTP in the lesion mask
> indexN <- which(carto$MASK_DWI_t0 == 1)
> seed <- indexN[which.max(carto[indexN,"TTP_t0"])] 
> 
> ## Display step by step the GR algorithm with sigma = 1
> for(iter in c(0,1,2,5,10)){
+   resGR1 <- calcGR(contrast = carto$TTP.th_t0, W = W, 
+                    seed = seed, sigma_max = 1, iter_max = iter, verbose = FALSE)
+   
+   multiplot(MRIaggr.Pat1_red, param = "TTP.th_t0", num = 3,hemisphere = "lesion", legend = FALSE,
+          breaks = seq(0,10,0.1), numeric2logical = TRUE, cex = 2,
+ 		 main = paste("iteration=",iter," - slice ", sep = ""),          
+          index1 = list(coords = coords[resGR1$GR,], pch = 20, cex = 1),
+          index2 = list(coords = coords[seed,], pch = 20, cex = 1)
+   )
+ }
> 
> ## Not run: 
> ##D ## GR with sigma = 2.2
> ##D resGR2 <- calcGR(contrast = carto$TTP.th_t0, W = W, 
> ##D                  seed = seed, sigma_max = 2.2, iter_max = 50,
> ##D                  history.step = TRUE, history.front = TRUE)
> ##D 
> ##D ## display 
> ##D # display the GR over the raw contrast
> ##D multiplot(MRIaggr.Pat1_red, param = "TTP.th_t0", num = 3, hemisphere = "lesion", legend = FALSE,
> ##D              breaks = seq(0,10,0.1), numeric2logical = TRUE, cex = 2,            
> ##D              index1 = list(coords = coords[resGR2$GR,], pch = 20, cex = 1)
> ##D )
> ##D 
> ##D # display the step of inclusion in GR group for each observation
> ##D multiplot(coords[resGR2$GR,],
> ##D              resGR2$history.step,breaks = 0:10,
> ##D              index1=list(coords = coords[seed,]),
> ##D              palette = rainbow(10)
> ##D )
> ##D 
> ##D # display the front propagation 
> ##D multiplot(coords[resGR2$GR,],
> ##D              resGR2$Mfront[,7],
> ##D              index1 = list(coords = coords[seed,])
> ##D )
> ## End(Not run)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>