1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
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: #
|
