# File         clusterSpaceFill
# By           Daniel B. Carr
# Copyright    2005, 2006 class use permitted

#  clusterSize is function is used in the main function below

# clusterSize function____________________________________
clusterSize = function(mat){

# input:  merge matrix transposed 
# output:  combined group size for each merge step
   n = ncol(mat)
   vec = rep(0,n)
   for( i in 1:n){
      j = mat[1,i]
      k = mat[2,i]
   # if k <0  both are singletons
   # if 1 singleton it must be j   
      if(k < 0) vec[i] =   2 else    
      if(j < 0) vec[i] =   1+ vec[k] else vec[i] = vec[j]+vec[k]
   }
   return(vec)
}


clusterSpaceFill = function(clusSummary,boxBoundsRLTB=c(1284,1,1024,1),drawHeight=-1){
# Input
# clusSummary      hclust() output list
# boxBoundsRLTB    Bounds for the region to be split, for example in pixels
#                  The order is right, left, top, bottom  	
# drawHeight	 Which boxes have merge heights less than this 

# Output
# points           Locations for plotting the items
# pointsBoxes      Boxes containing the points: right, left, top, bottom
# boxes            Box sides for recursively splitting the space
#                  right, left, top, bottom
# drawBox          Logical to draw boxes for merges less than drawHeight 			

# compute vector of group sizes at each merge
mat = t(clusSummary$merge)
clusSize = clusterSize(mat)   

# allocate storage 
n = ncol(mat)
node.rect = matrix(0,nrow=4,ncol=n)     # location for box sides
p.rect = matrix(0,nrow=4,ncol=n+1)      # box sides for points
p.point = matrix(0,nrow=2,ncol=n+1)     # points (x,y) average of box sides
dis = clusSummary$height
notdraw = rep(T,n)
contrast = c(1,-1,-1,1)
cont = matrix(c(.5,.5,0,0,0,0,.5,.5),nrow=2,byrow=T)
node.rect[,n] = boxBoundsRLTB

#_______________main loop____________________

y = n
for (i in 1:n){
   node = y[1]
   r = node.rect[,node]
   r1 = r
   r2 = r
   y = y[-1]
   if(dis[node] < drawHeight & notdraw[node])notdraw[node] = F  # height is below a bound
   leaf = mat[,node]
   j = leaf[1]
   k = leaf[2]

   if(j < 0) jsize = 1 else jsize = clusSize[j]
   if(k < 0) ksize = 1 else ksize = clusSize[k]

# maybe use total clusSize[node]
   splrat = jsize/(jsize+ksize)

   if(r%*%contrast > 0) {
      r1[1] = r1[2]+ (r[1]-r[2])*splrat
      r2[2] = r1[1]
   } else {
     r1[3] = r1[4]+ (r[3]-r[4])*splrat
     r2[4] = r1[3]
   }

   if(jsize==1){
       j = abs(j)
       p.point[,j] = cont%*%r1
       if(notdraw[node]) p.rect[,j] = r1  
   } else {
       notdraw[j] = notdraw[node]
       node.rect[,j] = r1
   }

   if(ksize==1){
       k = abs(k)
       p.point[,k] = cont%*% r2
       if(notdraw[node]) p.rect[,k] = r2
   } else {
       notdraw[k] = notdraw[node]
       node.rect[,k] = r2
   }

   y = c(y,leaf[leaf>0])
}

return(list(points=t(p.point),pointBoxes=t(p.rect),
            clusterBoxes=t(node.rect),drawBox=!notdraw)) 	
}