Finding a pin in a haystack - PCA image filtering

I found the following post regarding the anomalous metal object observed in a Curiosity Rover photo to be fascinating - specifically, the clever ways that some programmers used for filtering the image for the object. The following answer on mathematica.stackexchange.com was especially illuminating for its use of a multivariate distribution to describe the color channels for a test region of "sand". This distribution was subsequently used to assess if the rest of the image colors belonged to the same distribution.

I tried a different approach using a Principal Component Analysis (PCA) filter (above), also based on a region of sand. I believe the PCs can be understood in the following way: the PCs represent dominant rgb colors (below), while the loadings indicate the color intensity.

 

The first PC is obviously the main color of the sand and explains 99.949 % of the variance in the colors. Both shadowed areas and sun-lit areas are fairly equally masked out by subtracting a reconstructed image based on PC1, since they are of similar color with differing intensities. What remains are the non-sand regions.

The plot at the top, shows the root mean square error (RMS) of each pixel after PC1 has been filtered out. Furthermore, by using the distribution of the Sand RMS, one can gauge whether other pixels are likely to belong to this distribution. By using a probability > 0.0001 cutoff, the following plot shows the areas that fall outside in terms of RMS (in black).

Or, alternately, a plot of the probabilities themselves.

The last figure requires the image.scale() function found here.

Code to reproduce:

#image source - "http://www.nasa.gov/images/content/694811main_pia16225-43_full.jpg"
require(biOps)
source("image.scale.R") # function found here: http://menugget.blogspot.de/2011/08/adding-scale-to-image-plot.html
 
full.img <- readJpeg("curiosity.jpg")
full.dim <- dim(full.img)
sand.rows <- c(1, 400)
sand.cols <- c(1000, 1345)
obj.rows <- c(700, 1000)
obj.cols <- c(670, 1000)
 
sand.img <- imagedata(full.img[sand.rows[1]:sand.rows[2], sand.cols[1]:sand.cols[2], ])
sand.dim <- dim(sand.img)
obj.img <- imagedata(full.img[obj.rows[1]:obj.rows[2], obj.cols[1]:obj.cols[2], ])
obj.dim <- dim(obj.img)
 
sand.xs <- c(sand.cols[1], sand.cols[2], sand.cols[2], sand.cols[1])
sand.ys <- c(full.dim[1]-sand.rows[2], full.dim[1]-sand.rows[2], full.dim[1]-sand.rows[1], full.dim[1]-sand.rows[1])
obj.xs <- c(obj.cols[1], obj.cols[2], obj.cols[2], obj.cols[1])
obj.ys <- c(full.dim[1]-obj.rows[2], full.dim[1]-obj.rows[2], full.dim[1]-obj.rows[1], full.dim[1]-obj.rows[1])
 
 
#RGB vals
sand.rgb <- matrix(c(sand.img), nrow=prod(sand.dim[c(1:2)]), ncol=3)
colnames(sand.rgb) <- c("r", "g", "b")
sand.rgb <- as.data.frame(sand.rgb)
 
full.rgb <- matrix(c(full.img), nrow=prod(dim(full.img)[c(1:2)]), ncol=3)
colnames(full.rgb) <- c("r", "g", "b")
full.rgb <- as.data.frame(full.rgb)
 
obj.rgb <- matrix(c(obj.img), nrow=prod(dim(obj.img)[c(1:2)]), ncol=3)
colnames(obj.rgb) <- c("r", "g", "b")
obj.rgb <- as.data.frame(obj.rgb)
 
 
 
#PCA model
S <- prcomp(~ r + g + b, data=sand.rgb, center=FALSE, SCALE=FALSE)
S$sdev[1]^2 / sum(S$sdev^2) * 100 # 1st PC explains 99.9% of data
 
#prediction
usePC <- 1
 
PRED <- predict(S, sand.rgb)
RECON <- as.matrix(PRED[,usePC]) %*% t(as.matrix(S$rotation[,usePC]))
sand.rms <- sqrt(rowMeans((sand.rgb - RECON)^2))
dim(sand.rms) <- sand.dim[1:2]
sand.rms <- t(sand.rms)[,rev(seq(sand.dim[1]))]
image(sand.rms)
 
PRED <- predict(S, obj.rgb)
RECON <- as.matrix(PRED[,usePC]) %*% t(as.matrix(S$rotation[,usePC]))
obj.rms <- sqrt(rowMeans((obj.rgb - RECON)^2))
dim(obj.rms) <- obj.dim[1:2]
obj.rms <- t(obj.rms)[,rev(seq(obj.dim[1]))]
image(obj.rms)
 
PRED <- predict(S, full.rgb)
RECON <- as.matrix(PRED[,usePC]) %*% t(as.matrix(S$rotation[,usePC]))
full.rms <- sqrt(rowMeans((full.rgb - RECON)^2))
dim(full.rms) <- full.dim[1:2]
full.rms <- t(full.rms)[,rev(seq(full.dim[1]))]
image(full.rms)
 
 
 
#plot 1
ZLIM <- range(full.rms)
PAL <- colorRampPalette((c("grey10", "red", "yellow", "white")), bias=1)
NCOL=100
 
WIDTHS <- c(2,2,2)
HEIGHTS <- c(2,2,2,2)
OMI <- c(0.1,0.1,0.3,0.1)
RES <- 200
WIDTH <- sum(WIDTHS) + sum(OMI[c(2,4)])
HEIGHT <- sum(HEIGHTS) + sum(OMI[c(1,3)])
LAYOUT <- matrix(c(1,1,2,1,1,3,4,4,5,4,4,6), nrow=4, ncol=3, byrow=TRUE)
 
png("curiosity_w_filter.png", width=WIDTH, height=HEIGHT, res=RES, units="in")
#x11(width=WIDTH, height=HEIGHT)
layout(LAYOUT, widths=WIDTHS, heights=HEIGHTS, respect=TRUE)
#layout.show(6)
par(mar=c(1,1,1,1), omi=OMI)
 
plot(full.img)
mtext("Original Image", side=3, line=0)
polygon(x=sand.xs, y=sand.ys, lwd=2, lty=2, border="white")
polygon(x=obj.xs, y=obj.ys, lwd=2, lty=2, border="white")
text(x=mean(sand.xs), y=min(sand.ys) + c(max(sand.ys)-min(sand.ys))*0.9, labels="SAND", col="white")
text(x=mean(obj.xs), y=min(obj.ys) + c(max(obj.ys)-min(obj.ys))*0.9, labels="OBJECT", col="white")
 
plot(sand.img)
mtext("SAND", side=3, line=0)
plot(obj.img)
mtext("OBJECT", side=3, line=0)
 
image(x=seq(full.dim[2]), y=seq(full.dim[1]), z=full.rms, zlim=ZLIM, col=PAL(NCOL), yaxt="n", xaxt="n", xlab="", ylab="", bty="n", asp=1)
mtext("Sand PCA1 - Filtered Image", side=3, line=0)
polygon(x=sand.xs, y=sand.ys, lwd=2, lty=2, border="white")
polygon(x=obj.xs, y=obj.ys, lwd=2, lty=2, border="white")
text(x=mean(sand.xs), y=min(sand.ys) + c(max(sand.ys)-min(sand.ys))*0.9, labels="SAND", col="white")
text(x=mean(obj.xs), y=min(obj.ys) + c(max(obj.ys)-min(obj.ys))*0.9, labels="OBJECT", col="white")
 
image(x=seq(sand.dim[2]), y=seq(sand.dim[1]), z=sand.rms, zlim=ZLIM, col=PAL(NCOL), yaxt="n", xaxt="n", xlab="", ylab="", bty="n", asp=1)
mtext("SAND", side=3, line=0)
image(x=seq(obj.dim[2]), y=seq(obj.dim[1]), z=obj.rms, zlim=ZLIM, col=PAL(NCOL), yaxt="n", xaxt="n", xlab="", ylab="", bty="n", asp=1)
mtext("OBJECT", side=3, line=0)
 
dev.off()
 
 
#plot2
fun <- function(x){x <- abs(x); rgb(x[1], x[2], x[3])}
COL <- apply(S$rotation, 2, fun)
png("PC_colors.png", width=2, height=1, units="in", res=200)
par(mar=c(0,0,0,0), cex=0.5)
plot(1:3, rep(1,3), col=COL, pch=".", asp=1, cex=100, xlim=c(0.5,3.5), ylim=c(0.5,1.5), xaxt="n", yaxt="n", xlab="", ylab="")
text(1:3, rep(1,3),labels=paste(paste("PC", 1:3, sep=""), "\n", round(S$sdev^2 / sum(S$sdev^2) * 100, 3), "%"), col="white")
dev.off()
########################
 
#Probability of being sand
require(MASS)
 
fit <- fitdistr(c(sand.rms), "gamma") # fit a gamma distribution to sand rms
 
ps <- c(1 - rev(10^seq(-5,0,,100)))
qs <- quantile(c(sand.rms), prob=ps)
d <- dgamma(qs, shape=fit$estimate[1], rate=fit$estimate[2], log = FALSE)
hist(sand.rms, n=100, freq=FALSE)
lines(d ~ qs)
 
hist(full.rms, n=100, freq=FALSE, ylim=c(0,0.5))
lines(d ~ qs)
 
#probability of rms belonging to sand rms distribution
prob.sand <- dgamma(full.rms, fit$estimate[1], rate=fit$estimate[2], log = FALSE)
 
###position of object
image(x=seq(full.dim[2]), y=seq(full.dim[1]), z=prob.sand)
#click on object
POS <- locator(1)
 
#replace by c(0,1) at 1e-4 cutoff
prob.sand2 <- prob.sand
prob.sand2 <- replace(prob.sand2, which(prob.sand2 > 1e-4), 1)
prob.sand2 <- replace(prob.sand2, which(prob.sand2 != 1), 0)
 
 
#plot prob cutoff
WIDTHS <- c(5)
HEIGHTS <- c(4)
OMI <- c(0.2,0.2,0.2,0.2)
RES <- 200
WIDTH <- sum(WIDTHS) + sum(OMI[c(2,4)])
HEIGHT <- sum(HEIGHTS) + sum(OMI[c(1,3)])
COLS <- c(1, rgb(abs(S$rotation[1,1]), abs(S$rotation[1,2]), abs(S$rotation[1,3])))
 
png("curiosity_sand_prob_cutoff.png", width=WIDTH, height=HEIGHT, res=RES, units="in")
#x11(width=WIDTH, height=HEIGHT)
 
par(mar=c(0,0,0,0), omi=OMI)
image(x=seq(full.dim[2]), y=seq(full.dim[1]), z=prob.sand2, col=COLS, asp=1, xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
points(POS, cex=4, col="white", lwd=2)
 
mtext("'Sand' (Prob. > 0.0001)", side=3, line=0)
 
dev.off()
 
 
#plot prob
PAL <- colorRampPalette(rev(c("grey10", "red", "yellow", "white")), bias=1)
BREAKS <- 10^seq(-7, 0)
 
WIDTHS <- c(5)
HEIGHTS <- c(4,1)
OMI <- c(0.2,0.2,0.2,0.2)
RES <- 200
WIDTH <- sum(WIDTHS) + sum(OMI[c(2,4)])
HEIGHT <- sum(HEIGHTS) + sum(OMI[c(1,3)])
LAYOUT <- matrix(c(1,2), nrow=2, ncol=1, byrow=TRUE)
 
png("curiosity_sand_prob.png", width=WIDTH, height=HEIGHT, res=RES, units="in")
#x11(width=WIDTH, height=HEIGHT)
layout(LAYOUT, widths=WIDTHS, heights=HEIGHTS, respect=TRUE)
#layout.show(2)
 
par(mar=c(1,1,1,1), omi=OMI)
image(x=seq(full.dim[2]), y=seq(full.dim[1]), z=prob.sand, col=PAL(length(BREAKS)-1), breaks=BREAKS, xaxt="n", yaxt="n", xlab="", ylab="")
points(POS, cex=4, col="cyan", lwd=2)
box()
 
par(mar=c(4,1,0,1))
image.scale(1, breaks=BREAKS, log="x", col=PAL(length(BREAKS)-1), xaxt="n", yaxt="n", xlab="", ylab="")
axis(1, at=BREAKS)
box()
mtext("'Sand' Probability", side=1, line=3)
 
dev.off()

Created by Pretty R at inside-R.org