stMincut

rm(list=ls())
library(igraph)
library(Matrix)


# Helper functions for graph creation -------------------------------------

getSubscripts <- function(idx,msize)
{
  idx<-idx - 1;
  fsize <- msize[1] * msize[2]
  z <- (idx %/% fsize)
  idx <- idx - (z*fsize)
  y <- idx %/% msize[1]
  x<-idx - (y*msize[1])
  return(c(x+1,y+1,z+1))
}
getNeighbourIndex<-function(idx,msize,dX,dY)
{
  s<-getSubscripts(idx,msize)
  neighbours<-rep(NA,6)
  neighbours[1]<-getIndex(s[1]+1,s[2],s[3],msize,base=TRUE) #voxel below
  neighbours[2]<-getIndex(s[1]+dX,s[2]+1,s[3],msize,base=TRUE) #next column
  neighbours[3]<-getIndex(s[1]+dX,s[2]-1,s[3],msize,base=TRUE) #prev column
  neighbours[4]<-getIndex(s[1]+dY,s[2],s[3]+1,msize,base=TRUE) #next frame
  neighbours[5]<-getIndex(s[1]+dY,s[2],s[3]-1,msize,base=TRUE) #prev frame
  neighbours[6]<-idx
  return(neighbours)
}
getIndex <- function(x,y,z,msize,base=TRUE)
{
  if(base)
  {
    x[x<1]<-NA
  }else{
    x[x<1]<-NA
  }

  if(base)
  {
    x[x>msize[1]]<-msize[1]
    #      x[x>msize[1]]<-NA
  }else{
    x[x>msize[1]]<-NA
  }
  y[y<1]<-NA
  y[y>msize[2]]<-NA
  z[z<1]<-NA
  if(is.na(msize[3]))
  {
    z[z>1]<-NA
  }else{
    z[z>msize[3]]<-NA
  }

  fsize <- msize[1] * msize[2]
  idx = (z-1) * fsize
  idx <- (y - 1) * msize[1] + idx
  idx <- x + idx;
  return(idx)
}
getBaseNeighbourIndex<-function(idx,msize)
{
  s<-getSubscripts(idx,msize)
  if(s[1]!=msize[1])
    #if(s[1]!=1)
  {
    stop('Not a base layer voxel')
  }
  neighbours<-rep(NA,4)
  neighbours[1]<-getIndex(s[1],s[2]+1,s[3],msize,base=TRUE) #voxel next column
  neighbours[2]<-getIndex(s[1],s[2]-1,s[3],msize,base=TRUE) #voxel prev column
  neighbours[3]<-getIndex(s[1],s[2],s[3]+1,msize,base=TRUE) #next frame
  neighbours[4]<-getIndex(s[1],s[2],s[3]-1,msize,base=TRUE) #prev frame
  return(neighbours)
}
# Create artificial image -------------------------------------------------
img<-matrix(0,nrow=6,ncol=5)

szImg<-dim(img)
numEl<-szImg[1]*szImg[2]


img[c(3,5),]<-1

sources_S1<-vector(length=(numEl*6))
targets_S1<-vector(length=(numEl*6))


mnames<-matrix(nrow=numEl,ncol=3)


# build the internal n-links ----------------------------------------------

dX_s1<-1
dY_s1<-1
for(iVoxel in 1:numEl)
  {
  if(iVoxel %% 1000 == 0)
    {
    cat(sprintf('Processing Voxel:%i\n',iVoxel))
    }
  subscript <- getSubscripts(iVoxel,szImg)
  mnames[iVoxel,]<-subscript
  neighbours<-getNeighbourIndex(iVoxel,szImg,dX_s1,dY_s1)
  sources_S1[((iVoxel-1)*6 + 1):(iVoxel*6)]<-iVoxel
  targets_S1[((iVoxel-1)*6 + 1):(iVoxel*6)]<-neighbours
  }

blanks<-is.na(targets_S1)
sources_S1<-sources_S1[!blanks]
targets_S1<-targets_S1[!blanks]

#make the bottom layer (Vb) well connected
sources_base<-vector(length=(numEl / szImg[1])*4)
targets_base<-vector(length=(numEl / szImg[1])*4)
baseVoxels<-seq(from=szImg[1],to=numEl,by=szImg[1])

for(iloc in 1:length(baseVoxels))
  {
  iVoxel=baseVoxels[iloc]
  neighbours<-getBaseNeighbourIndex(iVoxel,szImg)
  sources_base[((iloc-1)*4 + 1):(iloc*4)]<-iVoxel
  targets_base[((iloc-1)*4 + 1):(iloc*4)]<-neighbours
  }
blanks<-is.na(targets_base)
sources_base<-sources_base[!blanks]
targets_base<-targets_base[!blanks]

sources<-c(sources_S1,sources_base)
targets<-c(targets_S1,targets_base)


# create the adjacency matrix ---------------------------------------------
adjM<-sparseMatrix(i=sources,j=targets,x=10000)


# create the inital MRF ---------------------------------------------------
g<-graph.adjacency(adjM,mode='directed',weighted=TRUE)

g<-simplify(g)  #remove the edge loops

g<-set.vertex.attribute(g,'r',1:numEl,mnames[,1])
g<-set.vertex.attribute(g,'c',1:numEl,mnames[,2])
g<-set.vertex.attribute(g,'f',1:numEl,mnames[,3])
g<-set.vertex.attribute(g,'s',1:numEl,1)
g<-set.vertex.attribute(g,'color',1:numEl,'green')
g<-set.vertex.attribute(g,'name',1:vcount(g),sprintf('s=%i,r=%i,c=%i,f=%i',V(g)$s,V(g)$r,V(g)$c,V(g)$f))

# add the source and sink vertices {s,t} ----------------------------------
g<-g+vertex('s',f=0,r=0,c=0,id='s',name='s')
g<-g+vertex('t',f=0,r=0,c=0,id='t',name='t')


# generate the costs for the voxels ---------------------------------------

costS1<-apply(img,2,cumsum)

# normalise costs to be between -1 and +1
costS1<-(costS1-median(costS1))/((max(costS1)-min(costS1))/2)


#voxels with a low cost are likely to be of interest, they should be connected to the source

ri<-which(costS1<0)
g['s',ri,attr='weight']<-abs(costS1[ri])


#join the base layer (Vb) to the source vertex
rVb<-seq(from=szImg[1],to=numEl,by=szImg[1])

g['s',rVb,attr='weight']<-2

#join the other vertices to the terminal vertex
rOther<-setdiff(1:numEl,c(ri,rVb))
g[rOther,'t',attr='weight']<-costS1[rOther]


cut<-stMincuts(g,'s','t')