#############################################################
#              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 ] 


#############################################