if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("EBImage")
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.1686275
Load images from arrays of pixel values
zipcode.RData
contains 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.869
display(img)
options("EBImage.display"= "raster")
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))
img_zip <- Image(vec, dim=c(16, 16))
display(img_zip)
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 2
Image data can be accessed as a plain R array using the imageData accessor,
dim(img)
## [1] 3499 2332 3
imageData(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.09803922
The distribution of pixel intensities can be plotted in a histogram, and their range inspected using the range function.
hist(img)
img_small <- resize(img, 128, 128)
display(img_small)
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)
img_bright <- img + 0.2
img_dark <- img - 0.2
display(combine(img_bright, img_dark), all=TRUE)
img_low <- img * 0.5
img_high <- img * 2
display(combine(img_low, img_high), all=TRUE)
display(img[300:450, 50:200,])
img_rotate <- translate(rotate(img, 45), c(50, 0))
display(img_rotate)
display(channel(img, mode='gray'))
display(channel(img, mode='luminance'))
display(channel(img_zip, mode='asred'))
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.
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')
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')
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)
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)