RunClustering: Clustering of a dynamical image sequence

Description

Clusters dynamics of an image sequence. Assumes prior sequence denoising.

The clustering procedure is an iterative procedure. At each step, the (available) children of the voxel associated to a largest neighborhoods (result of the denoising step or made of its children as a result of the getChildren function) are used to build a robust cluster. The center of the latter is compared to previously build ones using the MultiTestH0 procedure. The clusters with equivalent centers are merged together until no further merging are possible. The resulting cluster is added to the cluster list if the number of outliers generated by the robust cluster procedure does not exceed the original number of children. If not the voxel is added to the closest existing cluster.

Further details about the clustering procedure can be found in the references.

Usage

RunClustering(data.array, denois, nproc = 1, min.size = 1, alpha = 0.05,
  do.children.first = FALSE)

Arguments

data.array

a (2D or 3D)+T array containing the dynamic sequence of images (the dataset). The last dimension is the time.

denois

the result of the denoising procedure RunDenoising.

nproc

a numeric value indicating the number of processors used for parallel computation. Default set to 1 (no parallelization).

min.size

a numeric value indicating the smallest size of the neighborhood for a voxel to be clusterized. Default set to 1 (no limitation).

alpha

a numeric value indicating the global level of the multitest. Default set to 0.05.

do.children.first

an boolean. If TRUE compute children list for all voxels before starting the clusterization which will use these lists as neighborhoods.If FALSE (default) neighborhood are those resulting from the denoising step.

Value

a list containing:

cluster, a list of vectors---each vector contains the voxel indexes of one cluster.
centers, a matrix whose columns contain the average dynamics of each cluster.
bad.voxels, vector of non clusterized voxel indexes.

References

Rozenholc, Y. and Reiss, M. (2012) Preserving time structures while denoising a dynamical image, Mathematical Methods for Signal and Image Analysis and Representation (Chapter 12), Florack, L. and Duits, R. and Jongbloed, G. and van~Lieshout, M.-C. and Davies, L. Ed., Springer-Verlag, Berlin

Lieury, T. and Pouzat, C. and Rozenholc, Y. (submitted) Spatial denoising and clustering of dynamical image sequence: application to DCE imaging in medicine and calcium imaging in neurons

Examples

Run this code

# NOT RUN {
# }
# NOT RUN {
    library(DynClust)
    
    ## use fluorescence calcium imaging of neurons performed with Fura 2 excited at 340 nm
    data('adu340_4small',package='DynClust')
    
    ## Gain of the CCD camera:
    G <- 0.146
    ## readout variance of the CCD camera:
    sro2 <- (16.4)^2
    ## Stabilization of the variance to get a normalized dataset (variance=1)
    FT <- 2*sqrt(adu340_4small/G + sro2)
    FT.range = range(FT)
    
    ## launches the denoising step on the dataset with a statistical level of 5%
    FT.den.tmp <- RunDenoising(FT,1,mask.size=NA,nproc=2)
    
    ## launches the clustering step on the dataset with a statistical level of 5%
    FT.clust.tmp  <- RunClustering(FT,FT.den.tmp,nproc=2)
    n.cluster <- length(FT.clust.tmp$clusters)
    print(paste(n.cluster,'clusters using variance set to',sqrt(FT.den.tmp$var),'^2'))
    
    ## get the classified version of the data array and the map of the clusters
    FT.clust.res <- GetClusteringResults(FT,FT.den.tmp,FT.clust.tmp)
    
    ## plotting results of the clusterization
    par(mfrow=c(2,2))
    image(FT.clust.res$clust.map,col=rainbow(n.cluster))
    title('Cluster map')
    matplot(FT.clust.res$clust.center,col=rainbow(n.cluster),type="l",lwd=0.1,lty=1)
    title('Cluster centers')

    ## and more: original and clustered slices at time 50
    image(FT[,,50],zlim=FT.range,col=grey(seq(0,1,length=n.cluster)))
    title('Original sequence at time 50')
    image(FT.clust.res$clust.array[,,50],zlim=FT.range,col=grey(seq(0,1,length=n.cluster)))
    title('Clustered sequence at time 50')

    ####################################################################################
    ## reapply clustering with twice the nominal variance: forces stronger clustering ##
    ####################################################################################

    ## launches the denoising step on the dataset with a statistical level of 5%
    FT.den.tmp <- RunDenoising(FT,2,mask.size=NA,nproc=2)

    ## launches the clustering step on the dataset with a statistical level of 5%
    FT.clust.tmp  <- RunClustering(FT,FT.den.tmp,nproc=2)
    n.cluster <- length(FT.clust.tmp$clusters)
    print(paste(n.cluster,'clusters using twice the nominal variance'))
    
    ## get the classified version of the data array and the map of the clusters
    FT.clust.res <- GetClusteringResults(FT,FT.den.tmp,FT.clust.tmp)
    
    ## plotting results of the clusterization
    par(mfrow=c(2,2))
    image(FT.clust.res$clust.map,col=rainbow(n.cluster))
    title('Cluster map')
    matplot(FT.clust.res$clust.center,col=rainbow(n.cluster),type="l",lwd=0.1,lty=1)
    title('Cluster centers')

    ## and more: original and clustered slices at time 50
    image(FT[,,50],zlim=FT.range,col=grey(seq(0,1,length=n.cluster)))
    title('Original sequence at time 50')
    image(FT.clust.res$clust.array[,,50],zlim=FT.range,col=grey(seq(0,1,length=n.cluster)))
    title('Clustered sequence at time 50')
# }