diff options
Diffstat (limited to 'r/scaleSaliency.R')
| -rw-r--r-- | r/scaleSaliency.R | 246 |
1 files changed, 246 insertions, 0 deletions
diff --git a/r/scaleSaliency.R b/r/scaleSaliency.R new file mode 100644 index 0000000..80172bc --- /dev/null +++ b/r/scaleSaliency.R @@ -0,0 +1,246 @@ +# -*- coding: utf-8 -*- +# +# Scale Saliency algorithm +# +# Reference: +# [1] T. Kadir & M. Brady, Saliency, Scale and Image Description. +# 2001, International Journal of Computer Vision, 45(2), 83-105 +# +# Aaron LI +# 2015/07/29 +# + + +# Calculate Shannon entropy from histogram of probability densities. +# Arguments: +# phist - histogram of probability density +# Return: +# entropy value +calcShannonEntropy <- function(phist) { + # Helper function to calculate the entropy + # Arguments: + # p - probability density of each bin + # Return: + # p * log(p) if p > 0, otherwise 0 + plogp <- function(p) { + if (p < 1e-10) { + return(0.0) + } else { + return(p * log(p)) + } + } + # calculate the entropy + entropy <- -sum(sapply(phist, plogp)) + return(entropy) +} + + +# Generate the template circular regions for each scale. +# Arguments: +# scale_min - minimum scale (radius of circle) +# scale_max - maximum scale +# Return: +# list of matrix which represent the regions of interest (with value of TRUE) +circleROI <- function(scale_min, scale_max) { + rows <- 2 * scale_max + 1 + cols <- rows + rc <- (rows + 1) / 2 # central row + cc <- (cols + 1) / 2 # central col + roi <- list() + for (s in scale_min:scale_max) { + radius2 <- s^2 + m <- matrix(0, nrow=rows, ncol=cols) + roi[[paste("scale", s, sep="")]] <- + ifelse(((row(m)-rc)^2 + (col(m)-cc)^2) <= radius2, + TRUE, FALSE) + } + return(roi) +} + + +# Calculate the scale saliencies for the 1D case: scalar image +# Arguments: +# img - input *scalar* image +# scale_min - minimum scale (pixels of radius of circle) +# scale_max - maximum scale (NOTE: scale_max >= scale_min+2) +# nbins - how many bins used for histogram +# progressbar - whether to display the progress bar +# Return: +# 6-column data.frame contains the scale saliencies results +calcScaleSaliency1D <- function(img, scale_min, scale_max, nbins, + progressbar=TRUE) { + # check scale range first: must have at least 3 scales + stopifnot(scale_max >= scale_min+2) + # get number of rows and columns + rows <- nrow(img) + cols <- ncol(img) + # determine the saliency calculation region of the image + # FIXME: how to deal with the boundaries??? + row_begin <- scale_max + 1 + col_begin <- scale_max + 1 + row_end <- rows - scale_max + col_end <- cols - scale_max + # templates of regions for each scale + roi <- circleROI(scale_min, scale_max) + # R data frame to store the saliency results + scaleSaliency <- data.frame(row=numeric(0), col=numeric(0), + scale=numeric(0), entropy=numeric(0), + disimilarity=numeric(0), saliency=numeric(0)) + # determine the breakpoints for histogram + hist_breaks <- (0:nbins) * (max(img) - min(img))/nbins + min(img) + if (progressbar) { + # progress bar + pb <- txtProgressBar(min=row_begin, max=row_end, style=3) + } + for (ri in row_begin:row_end) { + if (progressbar) { + # update progress bar + setTxtProgressBar(pb, ri) + } + for (ci in col_begin:col_end) { + # filter out the required size of image block, which is + # used to calculate its histogram, entropy, etc. + imgROI <- img[(ri-scale_max):(ri+scale_max), + (ci-scale_max):(ci+scale_max)] + # vectors to store entropies and distances + entropy <- numeric(scale_max-scale_min+1) + distance <- numeric(scale_max-scale_min+1) + # initial probability density for scale of 's-1' + scaleHistPr0 <- rep(0, nbins) + for (s in scale_min:scale_max) { + scaleROI <- roi[[paste("scale", s, sep="")]] + # NOTE: do not use 'breaks=nbins', since the number is a + # suggestion only and breakpoints will be set to 'prtty' + # values in this case. + scaleHist <- hist(imgROI[scaleROI], + breaks=hist_breaks, plot=FALSE) + scaleHistPr <- scaleHist$counts / sum(scaleHist$counts) + # calculate Shannon entropy + entropy[s-scale_min+1] <- calcShannonEntropy(scaleHistPr) + # FIXME: calculate distance of scales??? + distance[s-scale_min+1] <- sum(abs(scaleHistPr-scaleHistPr0)) + # save the probability density of current scale 's' + scaleHistPr0 <- scaleHistPr + } + # smooth the 'distance' vector to reduce the impacts of noise + distance1 <- c(distance[1], distance[1:(length(distance)-1)]) + distance2 <- c(distance[2:length(distance)], + distance[length(distance)]) + distance <- (distance1 + distance + distance2) / 3 + # find the peaks of entropy, and the corresponding scales + peakScale <- c(FALSE, + ((entropy[2:(length(entropy)-1)] > + entropy[1:(length(entropy)-2)]) & + (entropy[2:(length(entropy)-1)] > + entropy[3:length(entropy)])), + FALSE) + #cat("peakScale:", peakScale, "\n") + # calculate the inter-scale saliencies for each entropy peaks + for (s in (scale_min:scale_max)[peakScale]) { + scaleNorm <- s*s / (2*s - 1) + scaleEntropy <- entropy[s-scale_min+1] + disimilarity <- scaleNorm * distance[s-scale_min+1] + saliency <- scaleEntropy * disimilarity + scaleSaliency[nrow(scaleSaliency)+1, ] <- list(ri, ci, s, + scaleEntropy, + disimilarity, + saliency) + } + } + } + if (progressbar) { + # close progress bar + close(pb) + } + return(scaleSaliency) +} + + +# Simple greedy clustering algorithm to filter out salient regions. +# Arguments: +# ssaliency - saliency results from 'calcScaleSaliency*' +# ssaliency_th - inter-scale saliency threshold +# disimilarity_th - disimilarity threshold +# Return: +# clustered & filtered saliency regions +greedyCluster <- function(ssaliency, ssaliency_th, disimilarity_th) { + # filter by global saliency & inter-scale saliency threshold + ssaliency <- ssaliency[((ssaliency$saliency > ssaliency_th) & + (ssaliency$disimilarity > disimilarity_th)), ] + # sort in descending inter-scale saliency + ssaliency <- ssaliency[order(-ssaliency$saliency), ] + # cluster salienct points + clusteredSaliency <- ssaliency[NULL, ] + while (nrow(ssaliency) > 0) { + ss <- ssaliency[1, ] + clusteredSaliency[nrow(clusteredSaliency)+1, ] <- ss + distance2 <- (ssaliency$row - ss$row)^2 + (ssaliency$col - ss$col)^2 + # filter out the points inside the current salient circle + ssaliency <- ssaliency[(distance2 > ss$scale^2), ] + } + return(clusteredSaliency) +} + + +# Plot the image and salient regions with ggplot2 +# Arguments: +# img - input image +# saliency - saliency restults by clusteredSaliency() +plotSalientReg <- function(img, saliency) { + require(reshape2) + require(ggplot2) + plotCircle <- function(xc, yc, radius) { + theta <- seq(0, 2*pi, length.out=100) + gcircle <- annotate("path", + x=xc+radius*cos(theta), + y=yc+radius*sin(theta), + colour="green") + return(gcircle) + } + # plot the image + gp <- ggplot(melt(img), aes(Var2, -Var1, fill=value)) + geom_raster() + # add circles + for (i in 1:nrow(saliency)) { + ss <- saliency[i, ] + gcircle <- plotCircle(ss$col, -ss$row, ss$scale) + gp <- gp + gcircle + } + return(gp) +} + + +# Convert the scale saliency information to DS9 regions. +# +# NOTE: +# However, the rows and columns of the FITS matrix in R correspond +# to the X and Y axes in DS9, which is *swapped*. +# Thus the region width and height correspond to the row range and +# column range, respectively. +# +# Arguments: +# saliency - saliency restults by clusteredSaliency() +# Return: +# vector of DS9 region strings +saliency2region <- function(saliency) { + regions <- with(saliency, + paste("circle(", row, ",", col, ",", scale, ")", + sep="")) + return(regions) +} + + +# Write DS9 region to file with appropriate header information. +# +# Arguments: +# filename - output region file +# region - vector/list of region strings +save.region <- function(filename, region) { + rf <- file(filename, "w") + region.hdr <- c("# Region file format: DS9 version 4.1", "image") + writeLines(region.hdr, rf) + writeLines(region, rf) + close(rf) +} + + +# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=r: # |