Corneal Disc Image Analysis

Setup

To test the image analysis capabilities of the ‘imager’ R package, I decided to download an image of a corneal disc and experiment with it.

setwd('~/Desktop/Research_Volunteering/Corneal Image Analysis/')

# download and load appropriate packages
list.of.packages <- c("dplyr", "ggplot2", "imager")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

# Loading relevant packages
lapply(list.of.packages, library, character.only = TRUE)

# load image
im <- load.image("OCT of corneal disc.png")
plot(im)

# inverted image (might be useful for segmentation)
im_neg = max(im) - im
plot(im_neg)

Exploring imager Package Tutorials

Morphological Operations and Correcting for Varying Luminance

# Thesholding (normal)
par(mfrow = c(3,1))
threshold(im,"80%") %>% plot
threshold(im,"85%") %>% plot
threshold(im,"90%") %>% plot

# Thesholding (inverted)
threshold(im_neg,"20%") %>% plot
threshold(im_neg,"15%") %>% plot
threshold(im_neg,"10%") %>% plot

# Removing any luminance trends with R's built-in lm function so that a fixed threshold works better.
# The first step is to convert the image data to a data.frame, then fit a model:
df <- as.data.frame(im)

m <- lm(value ~ x + y,data=df) #linear trend
summary(m)

## 
## Call:
## lm(formula = value ~ x + y, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.15632 -0.07907 -0.03028  0.04578  0.86355 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.585e-01  4.764e-04  332.73   <2e-16 ***
## x           -6.827e-06  5.350e-07  -12.76   <2e-16 ***
## y           -3.845e-04  1.679e-06 -229.03   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1181 on 431841 degrees of freedom
## Multiple R-squared:  0.1086, Adjusted R-squared:  0.1086 
## F-statistic: 2.631e+04 on 2 and 431841 DF,  p-value: < 2.2e-16

layout(t(1:2))
im.f <- im-fitted(m)
plot(im,main="Before")
plot(im.f,main="After trend removal")

# Morphological operations (based on a threshold of 85%)
im.t <- threshold(im.f,"85%")
plot(im.t)

# Finding average pixel intensity of pixels above threshold
# Convert threshold to pixset
px <- as.pixset(1-im.t) 
layout(t(1:1))

plot(im)
fill(px,5) %>% clean(3) %>% highlight

# taking the pixel set from a threshold image and converting it into a data frame to select pixels from the original picture.
px_to_image <- as.cimg(px)
plot(px_to_image)

px.image_to_df <- as.data.frame(px_to_image)
head(px.image_to_df)

##   x y value
## 1 1 1     1
## 2 2 1     1
## 3 3 1     1
## 4 4 1     1
## 5 5 1     1
## 6 6 1     1

# in the pixset, a value of 0 is black and a value of 1 is white
table(px.image_to_df$value)

## 
##      0      1 
##  66311 365533

# If the pixset is 1, then assign the value of the original image to zero, otherwise, leave the value alone
values_filtered <- NA
for(i in 1:width(df)){
  if(px.image_to_df[i,3] == 1){
    values_filtered[i] = 0
  }else{
   values_filtered[i] = df[i,3] 
  }
}

# creating a dataframe with only the values above threshold
df_filtered <- mutate(df, values_filtered) %>% select(-value) %>% rename(value = values_filtered)
head(df_filtered)

##   x y value
## 1 1 1     0
## 2 2 1     0
## 3 3 1     0
## 4 4 1     0
## 5 5 1     0
## 6 6 1     0

im_filtered <- as.cimg(df_filtered)

# plots demontrating that the below threshold pixels were removed
par(mfrow = c(2,1))
plot(im)
plot(im_filtered)

# Average pixel intensity for pixels above threshold
average_above_threshold <- sum(df_filtered$value)/sum(ifelse(df_filtered$value!=0,1,0))
average_above_threshold 

## [1] 0.339059

Canny Edge Detector

# Getting rid of noise using Gaussian filtering
im_edge <- grayscale(im) %>% isoblur(2) #2 pix. blur

# Computing the image gradient:
gr <- imgradient(im_edge,"xy")
plot(gr,layout="row")

#Compute the gradient magnitude via:
mag <- with(gr,sqrt(x^2+y^2))
plot(mag)

# Use non-maxima suppression to clean up  picture. 
threshold(mag) %>% plot

#Going along the (normalised) gradient
#Xc(im) is an image containing the x coordinates of the image
nX <- Xc(im) + gr$x/mag 
nY <- Yc(im) + gr$y/mag

#nX and nY are not integer values, so we can't use them directly as indices.
#We can use interpolation, though:
val.fwd <- interp(mag,data.frame(x=as.vector(nX),y=as.vector(nY)))

# We can naturally also go backwards along the gradient:
nX <- Xc(im) - gr$x/mag 
nY <- Yc(im) - gr$y/mag
val.bwd <- interp(mag,data.frame(x=as.vector(nX),y=as.vector(nY)))

# Given these two sets of values, non-maxima suppression comes down to 
# killing all values that aren’t larger than their two neighbours along the flow:
throw <- (mag < val.bwd) | (mag < val.fwd)
mag[throw] <- 0
plot(mag)

# Hysteresis 
# The ultimate goal is to classify each pixel as edge/non-edge, 
# which means picking a threshold. Instead of picking a single threshold, 
# Canny suggested picking a double threshold t1,t2, with t1<t, where points with 
# magnitudes above t2 are called “strong edges” points with magnitudes 
# above t1 are called “weak edges”

# strong threshold
t2 <- quantile(mag,.96)
# weak threshold 
t1 <- quantile(mag,.90)
layout(t(1:2))

strong <- mag>t2
plot(strong,main="Initial set of strong edges")
weak <- mag %inr% c(t1,t2)
plot(weak,main="Initial set of weak edges")

layout(t(1:1))

# Hysteresis rescues such weak edges progressively. The usual implementation is iterative.
overlap <- grow(strong,3) & weak 
strong.new <- strong | overlap
plot(strong.new,main="New set of strong edges")