Basic Image Analysis
Chengliang Tang, Yuting Ma
Image Analysis Tool: EBImage
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("EBImage")Import/Export Image
Read image file
Read an image file as an 3D array of doubles ranging from 0 to 1. The third dimension is the slot for the three channels: Red, Green and Blue, or RGB.
library("EBImage")
img <- readImage("chicken.jpg")
print(img)## Image
## colorMode : Color
## storage.mode : double
## dim : 3499 2332 3
## frames.total : 3
## frames.render: 1
##
## imageData(object)[1:5,1:6,1]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.1568627 0.1607843 0.1647059 0.1686275 0.1725490 0.1764706
## [2,] 0.1647059 0.1647059 0.1647059 0.1686275 0.1725490 0.1764706
## [3,] 0.1686275 0.1686275 0.1686275 0.1686275 0.1725490 0.1725490
## [4,] 0.1725490 0.1725490 0.1725490 0.1725490 0.1686275 0.1686275
## [5,] 0.1725490 0.1725490 0.1725490 0.1725490 0.1686275 0.1686275Load images from arrays of pixel values
- Example: handwritten zipcode digits in greyscale of size \(16*16\).
zipcode.RDatacontains 1289 images, each of which is flattened into a 256-dim vector.
load("zipcode.RData")
dat[1:5, 1:5]## V1 V2 V3 V4 V5
## 1 -1 -1 -1 -1.000 -1.000
## 2 -1 -1 -1 -1.000 -1.000
## 3 -1 -1 -1 -1.000 -1.000
## 4 -1 -1 -1 -0.929 0.351
## 5 -1 -1 -1 -1.000 -0.869Display image
- With EBImage: display image using an interactive JavaScript viewer (shown in browser) or R’s built-in graphics capabilities
display(img)
- To display in R graphics:
options("EBImage.display"= "raster") - With default graphics package:
n_r <- n_c <- 16
vec <- dat[70,]
M <- matrix(vec, n_r, n_c)
MM <- M[,rev(1:ncol(M))]
par(mfrow=c(1,2))
image(x=1:n_r, y=1:n_c, z=M, axes = FALSE, xlab="", ylab="", col = grey(seq(0, 1, length = 256)))
image(x=1:n_r, y=1:n_c, z=MM, axes = FALSE, xlab="", ylab="", col = grey(seq(0, 1, length = 256)))
par(mfrow=c(1,1))- Conversion to Image object
img_zip <- Image(vec, dim=c(16, 16))
display(img_zip)
Image data representation
EBImage uses a package-specific class Image to store and process images. It extends the R base class array, and all EBImage functions can also be called directly on matrices and arrays.
str(img)## Formal class 'Image' [package "EBImage"] with 2 slots
## ..@ .Data : num [1:3499, 1:2332, 1:3] 0.157 0.165 0.169 0.173 0.173 ...
## ..@ colormode: int 2Image data can be accessed as a plain R array using the imageData accessor,
dim(img)## [1] 3499 2332 3imageData(img)[1:3, 1:6,]## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.1568627 0.1607843 0.1647059 0.1686275 0.172549 0.1764706
## [2,] 0.1647059 0.1647059 0.1647059 0.1686275 0.172549 0.1764706
## [3,] 0.1686275 0.1686275 0.1686275 0.1686275 0.172549 0.1725490
##
## , , 2
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.1176471 0.1215686 0.1254902 0.1294118 0.1333333 0.1372549
## [2,] 0.1254902 0.1254902 0.1254902 0.1294118 0.1333333 0.1372549
## [3,] 0.1294118 0.1294118 0.1294118 0.1294118 0.1333333 0.1333333
##
## , , 3
##
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.08235294 0.08627451 0.09019608 0.09411765 0.09803922 0.10196078
## [2,] 0.09019608 0.09019608 0.09019608 0.09411765 0.09803922 0.10196078
## [3,] 0.09411765 0.09411765 0.09411765 0.09411765 0.09803922 0.09803922The distribution of pixel intensities can be plotted in a histogram, and their range inspected using the range function.
hist(img)
Basic Image Operations
Resize image
img_small <- resize(img, 128, 128)
display(img_small)
Display muptiple images
img_dog <- readImage("dog.jpg")
img_dog <- resize(img_dog, 128, 128)
img_all <- EBImage::combine(img_small, img_dog)
display(img_all, all=TRUE)
img_all2 <- EBImage::combine(img_small, flip(img_small), flop(img_small))
display(img_all2, all=TRUE)
Adjust brightness
img_bright <- img + 0.2
img_dark <- img - 0.2
display(combine(img_bright, img_dark), all=TRUE)
Adjust contrast
img_low <- img * 0.5
img_high <- img * 2
display(combine(img_low, img_high), all=TRUE)
Crop image
display(img[300:450, 50:200,])
Spatial Transformation
img_rotate <- translate(rotate(img, 45), c(50, 0))
display(img_rotate)
Color Management
Transformation to grayscale
- Mode =“gray”/“grey”: converts a Color image into a Grayscale image, using uniform 1/3 RGB weights.
display(channel(img, mode='gray'))
- Mode = “luminance”: luminance-preserving Color to Grayscale conversion using CIE 1931 luminance weights: 0.2126 * R + 0.7152 * G + 0.0722 * B.
display(channel(img, mode='luminance'))
Transformation to other color channel
- Mode= “asred”, “asgreen”, “asblue”: converts a Grayscale image or an array into a Color image of the specified hue.
display(channel(img_zip, mode='asred'))
Filtering
Linear Filtering
A common preprocessing step involves cleaning up the images by removing local artifacts or noise through smoothing. An intuitive approach is to define a window of a selected size around each pixel and average the values within that neighborhood. After applying this procedure to all pixels, the new, smoothed image is obtained. Mathematically, this can be expressed as \[f'(x,y)=\frac{1}{N} \sum_{s=-a}^{a} \sum_{t=-a}^{a} f(x+s, y+t) \] where \(f'(x,y)\) is the value of the pixel at position \((x,y)\), and aa determines the window size, which is \((2a+1)\) in each direction. \(N=(2a+1)^2\) is the number of pixels averaged over, and \(f′\) is the new, smoothed image.
More generally, we can replace the moving average by a weighted average, using a weight function ww, which typically has the highest value at the window midpoint (s=t=0) and then decreases towards the edges. \[ (w∗f)(x,y)=\sum_{s=-\infty}^{+\infty}\sum_{t=-\infty}^{+\infty} w(s,t) f(x+s, y+s)\]
For notational convenience, we let the summations range from \(-\infty\)to \(\infty\), even if in practice the sums are finite and \(w\) has only a finite number of non-zero values. In fact, we can think of the weight function \(w\) as another image, and this operation is also called the convolution of the images \(f\) and \(w\), indicated by the the symbol \(∗\). Convolution is a linear operation in the sense that \(w∗(c_1f_1+c_2f_2)=c_1w∗f_1+c_2w∗f_2\) for any two images \(f_1\), \(f_2\) and numbers \(c_1\), \(c_2\).
In EBImage, the 2-dimensional convolution is implemented by the function filter2, and the auxiliary function makeBrush can be used to generate the weight function. In fact, filter2 does not directly perform the summation indicated in the equation above. Instead, it uses the Fast Fourier Transformation in a way that is mathematically equivalent but computationally more efficient.
w <- makeBrush(size = 31, shape = 'gaussian', sigma = 5)
plot(w[(nrow(w)+1)/2, ], ylab = "w", xlab = "", cex = 0.7)
img_flo <- filter2(img, w)
display(img_flo)
Here we have used a Gaussian filter of width 5 given by sigma. Other available filter shapes include “box” (default), “disc”, “diamond” and “line”, for some of which the kernel can be binary; see ?makeBrush for details.
If the filtered image contains multiple frames, the filter is applied to each frame separately. For convenience, images can be also smoothed using the wrapper function gblur which performs Gaussian smoothing with the filter size automatically adjusted to sigma.
In signal processing the operation of smoothing an image is referred to as low-pass filtering. High-pass filtering is the opposite operation which allows to detect edges and sharpen images. This can be done, for instance, using a Laplacian filter.
- Low-pass linear filtering: to blur images, remove noise…
f_low <- makeBrush(21, shape='disc', step=FALSE)
display(f_low, title='Disc filter')
f_low <- f_low/sum(f_low)
img_lowPass <- filter2(img, filter=f_low)
display(img_lowPass, title='Filtered image')
- High-pass Laplacian filtering: to detect edges, sharpen images
f_high <- matrix(1, nc=3, nr=3)
f_high[2,2] <- -8
img_highPass <- filter2(img, f_high)
display(img_highPass, title='Filtered image')
Median filter
Another approach to perform noise reduction is to apply a median filter, which is a non-linear technique as opposed to the low pass convolution filter described in the previous section. Median filtering is particularly effective in the case of speckle noise, and has the advantage of removing noise while preserving edges.
The local median filter works by scanning the image pixel by pixel, replacing each pixel by the median on of its neighbors inside a window of specified size. This filtering technique is provided in EBImage by the function medianFilter. We demonstrate its use by first corrupting the image with uniform noise, and reconstructing the original image by median filtering.
l = length(img)
n = l/10
pixels = sample(l, n)
img_noisy = img
img_noisy[pixels] = runif(n, min=0, max=0.5)
display(img_noisy)
img_median = medianFilter(img_noisy, 1)
display(img_median)
Thresholding
Adaptive thresholding
The idea of adaptive thresholding is that, compared to straightforward thresholding, the threshold is allowed to be different in different regions of the image. In this way, one can anticipate spatial dependencies of the underlying background signal caused, for instance, by uneven illumination or by stray signal from nearby bright objects.
Adaptive thresholding works by comparing each pixel’s intensity to the background determined from a local neighbourhood. This can be achieved by comparing the image to its smoothed version, where the filtering window is bigger than the typical size of objects we want to capture.
The function thresh returns the binary image resulting from the comparison between an image and its filtered version with a rectangular window.
img_bengal <- readImage("Bengal_1.jpg")
display(img_bengal)
img_bengal_bw <- channel(img_bengal, mode="gray")
img_seg1 <- thresh(img_bengal_bw, w=30, h=30, offset=0.06) # {f = matrix(1, nc=2*w+1, nr=2*h+1) ; f=f/sum(f) ; x>(filter2(x, f)+offset)}
img_seg2 <- thresh(img_bengal_bw, w=60, h=60, offset=0.06)
img_seg3 <- thresh(img_bengal_bw, w=60, h=60, offset=0.1)
img_seg <- combine(img_seg1, img_seg2, img_seg3)
display(img_seg, all=TRUE)