############################################################# # STAT 718A Computer Practical 4 ############################################################# # First of all we shall use R to read in an image and perform # some simple operations x <- read.csv("test.csv") x <- as.matrix (x) # Let's look at some of the data print( x[1:10,1:10] ) # remove the first column y <- x[,-1] ## Let's add some noise to the image, which will demonstrate ## the use of some simple image operations. ## The noise distribution is a mixture of two normals ## with high standard deviations nc<-dim(y)[2] nr<-dim(y)[1] for (i in 1:nr){ for (j in 1:nc){ uu<-runif(1) sigma<-50 if (uu < 0.4){sigma<- 300} z[i,j] <- y[i,j] + rnorm(1)*sigma } } # let's plot the noisy image image(z) # define some working images... z1<- z*0 z2<- z*0 z3<- z*0 z4<- z*0 # mean filter: first order neigborhood # This computes a mean over the 1st order neighbourhood for (i in 2:(nr-1)){ for (j in 2:(nc-1)){ z1[i,j] <- (z[i,j]+z[i-1,j]+z[i+1,j]+z[i,j-1]+z[i,j+1])/5 } } image( z1 ) # mean filter: second order neigborhood # This computes a mean over the 2nd order neighbourhood for (i in 2:(nr-1)){ for (j in 2:(nc-1)){ z2[i,j] <- (z[i,j]+z[i-1,j]+z[i+1,j]+z[i,j-1]+z[i,j+1] + z[i-1,j-1]+z[i+1,j+1]+z[i-1,j+1]+z[i+1,j-1]) / 9 } } image( z2 ) # median filter: second order neigborhood # This computes a median over the 2nd order neighbourhood for (i in 2:(nr-1)){ for (j in 2:(nc-1)){ z3[i,j] <- median(c(z[i,j],z[i-1,j],z[i+1,j],z[i,j-1],z[i,j+1], z[i-1,j-1],z[i+1,j+1],z[i-1,j+1],z[i+1,j-1]) ) } } image( z3 ) # Let's compare the images par(mfrow=c(2,2)) image(z,axes=FALSE) title("noisy image") image(z1,axes=FALSE) title("mean filter 1st order") image(z2,axes=FALSE) title("mean filter 2nd order") image(z3,axes=FALSE) title("median filter") #how close are the filtered images to the truth (image y) library(shapes) insider <- c(2:(nr-1) ) insidec <- c(2:(nc-1) ) norm( z[insider,insidec] - y[insider,insidec] ) norm( z1[insider,insidec] - y[insider,insidec] ) norm( z2[insider,insidec] - y[insider,insidec] ) norm( z3[insider,insidec] - y[insider,insidec] ) # Which is the best image in terms of distance to the truth? # plot the true image! image(y,axes=FALSE) # Laplacian filter: second order neigborhood # This filter emphasizes edges for (i in 2:(nr-1)){ for (j in 2:(nc-1)){ z4[i,j] <- (-8*y[i,j]+y[i-1,j]+y[i+1,j]+y[i,j-1]+y[i,j+1] + y[i-1,j-1]+y[i+1,j+1]+y[i-1,j+1]+y[i+1,j-1]) / 3 } } image( z4 ,axes=FALSE) # Threshold function thresh <- function ( z , t){ zt <- z zt [ z> t] <- 0 zt [ z<= t] <- 1 par(mfrow=c(1,2)) image(z,axes=FALSE) title("image") image(zt,axes=FALSE) title("thesholded") } # Histogram of all gray levels hist(z4) # Summary statistics mean(c(z4)) sd(c(z4)) summary(c(z4)) # good choice of threshold? Try t=100 thresh( z4, 100) ############################################ If you have time try looking at ImageJ. Open the image cerm.tif, apply a mean filter, a median filter and threshold the image. It's often a good idea to `Duplicate' the image to prevent writing on top of your work. Duplicate [Image > Duplicate] mean filter [ Process > Filters > Mean ] Try different choices of radius median filter [Process > Filters > Median] threshold [ Image > Adjust > Threshold ] #############################################