SOCR ≫ | DSPA ≫ | Topics ≫ |
In this DSPA Appendix, we present examples of image processing, spectral manipulation, and filtering.
## Warning: package 'knitr' was built under R version 4.1.2
r Biocpkg(“EBImage”) is an R package distributed as part of the Bioconductor project. To install the package, start R and enter:
install.packages("BiocManager")
::install("EBImage") BiocManager
Once the R package(“EBImage”) is installed, it can be loaded by the following command.
library("EBImage")
# install.packages("jpeg")
library(jpeg)
Basic R package(“EBImage”) functionality includes reading, writing, and displaying of images. Images are read using the function readImage
, which takes as input a filename or an URL. To start off, let us load a sample picture distributed with the package.
= system.file("images", "sample.png", package="EBImage")
f = readImage(f) img
The R package (“EBImage”) currently supports three image file formats: jpeg
, png
and tiff
. Additional image formats can be imported using the R GitHubPkg (“aoles/RBioFormats”), which adds support for a wide range of file formats including proprietary microscopy image data and metadata.
Imported images can be visualized by the function display()
.
display(img, method="browser")
In interactive R sessions, display
opens the image in a JavaScript viewer in a web browser tab. Mouse or keyboard shortcuts allow zooming in and out of the image, panning, and cycling through multiple image frames. Images can also be displayed using core R plotting methods, which allows for combining images with other plotting functionality, e.g., adding text labels on top of the image.
display(img, method="raster")
text(x = 20, y = 20, label = "Colorful Parrots", adj = c(0,1), col = "orange", cex = 2)
The graphics displayed in an R device can be saved using R package(“base”) R functions dev.print
or dev.copy
. For example, let’s save our annotated image as a JPEG file and verify its size on disk.
= "parrots.jpg"
filename dev.print(jpeg, filename = filename , width = dim(img)[1], height = dim(img)[2])
png
2
file.info(filename)$size
[1] 34311
The default behavior of display
can be globally changed by setting the "options("EBImage.display") to either
“browser”or
“raster”`. This is useful, for example, to preview images inside RStudio.
It is also possible to read and view color images,
= readImage(system.file("images", "sample-color.png", package="EBImage"))
imgcol display(imgcol)
Images may contain multiple frames, in which case they can be displayed all at once in a grid arrangement by specifying the function argument all = TRUE
.
= readImage(system.file("images", "nuclei.tif", package="EBImage"))
nuc display(nuc, method = "raster", all = TRUE)
Alternatively, single frames can be displayed.
In addition to importing images, images can be exported (saved) to files using the EBImage::writeImage()
. The image that we loaded was a PNG file. To save this image as a JPEG file, the JPEG format allows setting a quality value between 1 and 100 to reflect the desired level of image compression. The default value of the quality
argument of writeImage
is 100. Smaller values yield images of smaller size, but worse resolution (less detail).
writeImage(imgcol, "sample.jpeg", quality = 85)
Similarly, we could have saved the image as a TIFF file and set which compression algorithm we want to use. For a complete list of available parameters see ?writeImage
.
EBImage uses a package-specific class Image
to store and process images. It extends the R base class array
, and all R package(“EBImage”) functions can also be called directly on matrices and arrays. You can find out more about this class by typing ?Image
. Let us peek into the internal structure of an Image
object.
str(img)
Formal class 'Image' [package "EBImage"] with 2 slots
..@ .Data : num [1:768, 1:512] 0.447 0.451 0.463 0.455 0.463 ...
..@ colormode: int 0
..$ dim: int [1:2] 768 512
The .Data
slot contains a numeric array of pixel intensities. We see that in this case the array is two-dimensional, with 768 times 512 elements, and corresponds to the pixel width and height of the image. These dimensions can be accessed using the dim
function, just like for regular arrays.
dim(img)
[1] 768 512
Image data can be accessed as a plain R array
using the imageData
accessor.
imageData(img)[1:3, 1:6]
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.4470588 0.4627451 0.4784314 0.4980392 0.5137255 0.5294118
[2,] 0.4509804 0.4627451 0.4784314 0.4823529 0.5058824 0.5215686
[3,] 0.4627451 0.4666667 0.4823529 0.4980392 0.5137255 0.5137255
The as.array()
method can be used to coerce an Image
to an array
.
is.Image( as.array(img) )
[1] FALSE
The distribution of pixel intensities can be plotted in a histogram, and their range inspected using the range
function.
hist(img)
range(img)
[1] 0 1
A useful summary of Image
objects is also provided by the show
method, which is invoked if we simply type the object’s name.
img
Image
colorMode : Grayscale
storage.mode : double
dim : 768 512
frames.total : 1
frames.render: 1
imageData(object)[1:5,1:6]
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.4470588 0.4627451 0.4784314 0.4980392 0.5137255 0.5294118
[2,] 0.4509804 0.4627451 0.4784314 0.4823529 0.5058824 0.5215686
[3,] 0.4627451 0.4666667 0.4823529 0.4980392 0.5137255 0.5137255
[4,] 0.4549020 0.4666667 0.4862745 0.4980392 0.5176471 0.5411765
[5,] 0.4627451 0.4627451 0.4823529 0.4980392 0.5137255 0.5411765
For a more compact representation without the preview of the intensities array use the print
method with the argument short
set to TRUE
.
print(img, short=TRUE)
Image
colorMode : Grayscale
storage.mode : double
dim : 768 512
frames.total : 1
frames.render: 1
Color images are based on 3 channels (RBG).
print(imgcol, short=TRUE)
Image
colorMode : Color
storage.mode : double
dim : 768 512 3
frames.total : 3
frames.render: 1
They differ from gray-scale images by the property colorMode
and the number of dimensions, 3 (B&W) vs. 4 (color). The colorMode
slot turns out to be convenient when dealing with stacks of images. If it is set to gray-scale
, then the third and all higher dimensions of the array are considered as separate image frames corresponding, for instance, to different z-positions, time points, replicates, etc. On the other hand, if colorMode
is Color
, then the third dimension is assumed to hold different color channels, and only the fourth and higher dimensions are used for multiple image frames. imgcol
contains three color channels, which correspond to the red, green and blue intensities of the photograph. However, this does not necessarily need to be the case, and the number of color channels is arbitrary.
The “frames.total” and “frames.render” fields shown by the object summary correspond to the total number of frames contained in the image, and to the number of rendered frames. These numbers can be accessed using the function numberOfFrames
by specifying the type
argument.
numberOfFrames(imgcol, type = "render")
[1] 1
numberOfFrames(imgcol, type = "total")
[1] 3
Image frames can be extracted using getFrame
and getFrames
. getFrame
returns the i-th frame contained in the image y. If type
is "total"
, the function is unaware of the color mode and returns an xy-plane. For type="render"
the function returns the i-th image as shown by the display function. While getFrame
returns just a single frame, getFrames
retrieves a list of frames which can serve as input to lapply
-family functions. See the “Global thresholding” section for an illustration of this approach.
Let’s look at nuclear/cellular imaging data, which contains 4 total frames that correspond to the 4 separate gray-scale images, as indicated by “frames.render”.
nuc
Image
colorMode : Grayscale
storage.mode : double
dim : 510 510 4
frames.total : 4
frames.render: 4
imageData(object)[1:5,1:6,1]
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.06274510 0.07450980 0.07058824 0.08235294 0.10588235 0.09803922
[2,] 0.06274510 0.05882353 0.07843137 0.09019608 0.09019608 0.10588235
[3,] 0.06666667 0.06666667 0.08235294 0.07843137 0.09411765 0.09411765
[4,] 0.06666667 0.06666667 0.07058824 0.08627451 0.08627451 0.09803922
[5,] 0.05882353 0.06666667 0.07058824 0.08235294 0.09411765 0.10588235
As described in the previous section, the class Image
extends the base class array
and uses colorMode
to store how the color information of the multi-dimensional data should be handled.
The function colorMode
can be used to access and change this property, modifying the rendering mode of an image. For example, if we take a Color
image and change its mode to gray-scale
, then the image won’t display as a single color image anymore but rather as three separate gray-scale frames corresponding to the red, green and blue channels. The function colorMode
does not change the actual content of the image but only changes the way the image is rendered by R package(“EBImage”).
colorMode(imgcol) = Grayscale
display(imgcol, all=TRUE)
Color space conversions between gray-scale
and Color
images are performed using the function channel
. It has a flexible interface which allows to convert either way between the modes, and can be used to extract color channels. Unlike colorMode
, channel
changes the pixel intensity values of the image.
Color
to gray-scale
conversion modes include taking a uniform average across the RGB channels, and a weighted luminescence preserving conversion mode better suited for display purposes.
The asred
, asgreen
and asblue
modes convert a gray-scale image or array into a color image of the specified hue.
The convenience function toRGB
promotes a gray-scale image to RGB color space by replicating it across the red, green and blue channels, which is equivalent to calling channel
with mode set to rgb
. When displayed, this image doesn’t look different from its gray-scale origin, which is expected because the information between the color channels is the same. To combine three gray-scale images into a single rgb image use the function rgbImage
.
The function Image
can be used to construct a color image from a character vector or array of named R colors (as listed by colors()
) and/or hexadecimal strings of the form “#rrggbb” or “#rrggbbaa”.
= matrix(rep(c("red","green", "#0000ff"), 25), 5, 5)
colorMat = Image(colorMat)
colorImg colorImg
Image
colorMode : Color
storage.mode : double
dim : 5 5 3
frames.total : 3
frames.render: 1
imageData(object)[1:5,1:5,1]
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 1 0
[2,] 0 1 0 0 1
[3,] 0 0 1 0 0
[4,] 1 0 0 1 0
[5,] 0 1 0 0 1
display(colorImg, interpolate=FALSE)
Being numeric arrays, images can be conveniently manipulated by any of R’s arithmetic operators. For example, we can produce a negative image by simply subtracting the image from its maximum value.
= max(img) - img
img_neg display( img_neg )
We can also increase the brightness of an image through addition, adjust the contrast through multiplication, and apply gamma correction through exponentiation.
= combine(
img_comb
img,+ 0.3,
img * 2,
img ^ 0.5
img
)
display(img_comb, all=TRUE)
In the example above we have used combine
to merge individual images into a single multi-frame image object.
Furthermore, we can crop and threshold images with standard matrix operations.
= img[366:749, 58:441]
img_crop = img_crop > .5
img_thresh display(img_thresh)
The thresholding operation returns an Image
object with binarized pixel values. The R data type used to store such an image is logical
.
img_thresh
Image
colorMode : Grayscale
storage.mode : logical
dim : 384 384
frames.total : 1
frames.render: 1
imageData(object)[1:5,1:6]
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] FALSE FALSE FALSE FALSE FALSE FALSE
[2,] FALSE FALSE FALSE FALSE FALSE FALSE
[3,] FALSE FALSE FALSE FALSE FALSE FALSE
[4,] FALSE FALSE FALSE FALSE FALSE FALSE
[5,] FALSE FALSE FALSE FALSE FALSE FALSE
For image transposition, use transpose
rather than R’s R package(“base”) function t
. This is because the former one works also on color and multiframe images by swapping its spatial dimensions.
= transpose(img)
img_t display( img_t )
We just saw one type of spatial transformation, transposition, but there are many more, for example translation, rotation, reflection and scaling. The translate()
method moves the image plane by the specified two-dimensional vector in such a way that pixels that end up outside the image region are cropped, and pixels that enter into the image region are set to background.
= translate(img, c(100,-50))
img_translate display(img_translate)
The background color can be set using the argument bg.col
common to all relevant spatial transformation functions. The default sets the value of background pixels to zero which corresponds to black. Let us demonstrate the use of this argument with rotate
which rotates the image clockwise by the given angle.
= rotate(img, 30, bg.col = "white")
img_rotate display(img_rotate)
To scale an image to desired dimensions use resize
. If you provide only one of either width or height, the other dimension is automatically computed keeping the original aspect ratio.
= resize(img, w=256, h=256)
img_resize display(img_resize )
The functions flip
and flop
reflect the image around the image horizontal and vertical axis, respectively.
= flip(img)
img_flip = flop(img)
img_flop
display(combine(img_flip, img_flop), all=TRUE)
Spatial linear transformations are implemented using the general affine
transformation. It maps image pixel coordinates px
using a 3x2 transformation matrix m
in the following way: cbind(px, 1) %*% m
. For example, horizontal sheer mapping can be applied by
= matrix(c(1, -.5, 128, 0, 1, 0), nrow=3, ncol=2)
m = affine(img, m)
img_affine display( img_affine )
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 \(a\) 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 \(w\), 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 r Biocpkg(“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.
= makeBrush(size = 31, shape = 'gaussian', sigma = 5)
w plot(w[(nrow(w)+1)/2, ], ylab = "w", xlab = "", cex = 0.7)
= filter2(img, w)
img_flo 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
.
= gblur(nuc, sigma = 5)
nuc_gblur display(nuc_gblur, all=TRUE )
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.
= matrix(1, nrow = 3, ncol = 3)
fhi 2, 2] = -8
fhi[= filter2(img, fhi)
img_fhi display(img_fhi)
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 of its neighbors inside a window of specified size. This filtering technique is provided in R package(“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.
= length(img)
l = l/10
n = sample(l, n)
pixels = img
img_noisy = runif(n, min=0, max=1)
img_noisy[pixels] display(img_noisy)
= medianFilter(img_noisy, 1)
img_median display(img_median)
Binary images are images which contain only two sets of pixels, with values, say 0 and 1, representing the background and foreground pixels. Such images are subject to several non-linear morphological operations: erosion, dilation, opening, and closing. These operations work by overlaying a mask, called the structuring element, over the binary image in the following way:
erosion: For every foreground pixel, put the mask around it, and if any pixel covered by the mask is from the background, set the pixel to background.
dilation: For every background pixel, put the mask around it, and if any pixel covered by the mask is from the foreground, set the pixel to foreground.
= readImage(system.file('images', 'shapes.png', package='EBImage'))
shapes = shapes[110:512,1:130]
logo display(logo)
= makeBrush(5, shape='diamond')
kern display(kern, interpolate=FALSE)
= erode(logo, kern)
logo_erode= dilate(logo, kern)
logo_dilate
display(combine(logo_erode, logo_dilate), all=TRUE)
Opening and closing are combinations of the two operations above: opening performs erosion followed by dilation, while closing does the opposite, i.e, performs dilation followed by erosion. Opening is useful for morphological noise removal, as it removes small objects from the background, and closing can be used to fill small holes in the foreground. These operations are implemented by opening
and closing
.
In the “Manipulating images” section we have already demonstrated how to set a global threshold on an image. There we used an arbitrary cutoff value. For images whose distribution of pixel intensities follows a bi-modal histogram a more systematic approach involves using the Otsu’s method. Otsu’s method is a technique to automatically perform clustering-based image thresholding. Assuming a bi-modal intensity distribution, the algorithm separates image pixels into foreground and background. The optimal threshold value is determined by minimizing the combined intra-class variance.
Otsu’s threshold can be calculated using the function otsu
. When called on a multi-frame image, the threshold is calculated for each frame separately resulting in an output vector of length equal to the total number of frames in the image.
= otsu(nuc)
threshold threshold
[1] 0.3535156 0.4082031 0.3808594 0.4121094
= combine( mapply(function(frame, th) frame > th, getFrames(nuc), threshold, SIMPLIFY=FALSE) )
nuc_th display(nuc_th, all=TRUE)
Note the use of getFrames
to split the image into a list of individual frames, and combine
to merge the results back together.
The idea of adaptive thresholding is that, compared to straightforward thresholding from the previous section, 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 neighborhood. 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.
= makeBrush(31, "disc")
disc = disc / sum(disc)
disc = 0.05
offset = filter2( nuc, disc )
nuc_bg = nuc > nuc_bg + offset
nuc_th display(nuc_th, all=TRUE)
This technique assumes that the objects are relatively sparsely distributed in the image, so that the signal distribution in the neighborhood is dominated by background. While for the nuclei in our images this assumption makes sense, for other situations you may need to make different assumptions. The adaptive thresholding using a linear filter with a rectangular box is provided by thresh
, which uses a faster implementation compared to directly using filter2
.
display( thresh(nuc, w=15, h=15, offset=0.05), all=TRUE )
# FFT SHIFT
#' This function is useful for visualizing the Fourier transform with the zero-frequency
#' component in the middle of the spectrum.
#'
#' @param img_ff A Fourier transform of a 1D signal, 2D image, or 3D volume.
#' @param dim Number of dimensions (-1, 1, 2, 3).
#' @return A properly shifted FT of the array.
#'
<- function(img_ff, dim = -1) {
fftshift
<- dim(img_ff)[1]
rows <- dim(img_ff)[2]
cols # planes <- dim(img_ff)[3]
<- function(img_ff) {
swap_up_down <- ceiling(rows/2)
rows_half return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)]))
}
<- function(img_ff) {
swap_left_right <- ceiling(cols/2)
cols_half return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half]))
}
#swap_side2side <- function(img_ff) {
# planes_half <- ceiling(planes/2)
# return(cbind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)], img_ff[1:rows, 1:cols, 1:planes_half]))
#}
if (dim == -1) {
<- swap_up_down(img_ff)
img_ff return(swap_left_right(img_ff))
}else if (dim == 1) {
return(swap_up_down(img_ff))
}else if (dim == 2) {
return(swap_left_right(img_ff))
}else if (dim == 3) {
# Use the `abind` package to bind along any dimension a pair of multi-dimensional arrays
# install.packages("abind")
library(abind)
<- dim(img_ff)[3]
planes <- ceiling(rows/2)
rows_half <- ceiling(cols/2)
cols_half <- ceiling(planes/2)
planes_half
<- abind(img_ff[((rows_half+1):rows), (1:cols), (1:planes)],
img_ff 1:rows_half), (1:cols), (1:planes)], along=1)
img_ff[(<- abind(img_ff[1:rows, ((cols_half+1):cols), (1:planes)],
img_ff 1:rows, 1:cols_half, (1:planes)], along=2)
img_ff[<- abind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)],
img_ff 1:rows, 1:cols, 1:planes_half], along=3)
img_ff[return(img_ff)
}else {
stop("Invalid dimension parameter")
} }
str(img)
Formal class 'Image' [package "EBImage"] with 2 slots
..@ .Data : num [1:768, 1:512] 0.447 0.451 0.463 0.455 0.463 ...
..@ colormode: int 0
..$ dim: int [1:2] 768 512
<- matrix(img, nrow = dim(img)[1], ncol = dim(img)[2])
imgMat dim(imgMat)
[1] 768 512
<- fft(imgMat) # fftw2d # Display Re(FT):
ft_imgMat display(fftshift(ft_imgMat))
Warning in storage.mode(x) <- "double": imaginary parts discarded in coercion
<- sqrt(Re(ft_imgMat)^2+Im(ft_imgMat)^2)
mag_ft_imgMat display(mag_ft_imgMat)
# Phase <- atan(Im(img_ff)/Re(img_ff))
<- atan2(Im(ft_imgMat), Re(ft_imgMat))
phase_ft_imgMat # display(phase_ft_imgMat)
# Display FT
::display(log(fftshift(ft_imgMat)), title="FT Magnitude (Cyrillic Alphabet)") EBImage
Warning in storage.mode(x) <- "double": imaginary parts discarded in coercion
# IFT Reconstruction: Magnitude=mag_ft_imgMat AND Phase=phase_ft_imgMat
= mag_ft_imgMat * cos(phase_ft_imgMat)
Real = mag_ft_imgMat * sin(phase_ft_imgMat)
Imaginary = Re(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_imgMat))
ift_ft_imgMat ::display(ift_ft_imgMat, method = "raster", title="(IFT o FT) Magnitude=Img | Phase=Img") EBImage
Remove Fourier coefficients indexed above \(k_{x,max/2}\) and \(k_{y,max/2}\).
<- mag_ft_imgMat
mag_ft_imgMat_1 for (i in 1:dim(mag_ft_imgMat)[1]) {
for (j in 1:dim(mag_ft_imgMat)[2]) {
if (abs(i-(dim(mag_ft_imgMat)[1])/2) > (dim(mag_ft_imgMat)[1])/4 |
abs(j-(dim(mag_ft_imgMat)[2])/2) > (dim(mag_ft_imgMat)[2])/4) {
<- 0
mag_ft_imgMat_1[i,j]
}
}
}
::display(mag_ft_imgMat_1, method = "raster", title="(IFT o LowPassFilter o FT) Magnitude=Img | Phase=Img") EBImage
= mag_ft_imgMat_1 * cos(phase_ft_imgMat)
Real = mag_ft_imgMat_1 * sin(phase_ft_imgMat)
Imaginary = Re(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_imgMat_1))
ift_ft_imgMat_1 ::display(100*ift_ft_imgMat_1, method = "raster", title="(IFT o LowPassFilter o FT) Magnitude=Img | Phase=Img") EBImage
Remove the Fourier coefficients indexed below \(k_{x,max/4}\) and \(k_{y,max/4}\).
# High-pass frequency filtering, remove Fourier coeff's indexed below k_{x,max/4} and k_{y,max/4}
<- mag_ft_imgMat
mag_ft_imgMat_2 for (i in 1:dim(mag_ft_imgMat)[1]) {
for (j in 1:dim(mag_ft_imgMat)[2]) {
if (abs(i-(dim(mag_ft_imgMat)[1])/2) < (dim(mag_ft_imgMat)[1])/4 &
abs(j-(dim(mag_ft_imgMat)[2])/2) < (dim(mag_ft_imgMat)[2])/4) {
<- 0
mag_ft_imgMat_2[i,j]
}
}
}
::display(mag_ft_imgMat_2, method = "raster", title="(IFT o HighPassFilter o FT) Magnitude=Img | Phase=Img") EBImage
= mag_ft_imgMat_2 * cos(phase_ft_imgMat)
Real = mag_ft_imgMat_2 * sin(phase_ft_imgMat)
Imaginary = Re(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_imgMat_2))
ift_ft_imgMat_2 ::display(ift_ft_imgMat_2, method = "raster", title="(IFT o HighPassFilter o FT) Magnitude=Img | Phase=Img") EBImage
Image segmentation performs partitioning of an image, and is typically used to identify objects in an image. Non-touching connected objects can be segmented using the function bwlabel
, while watershed
and propagate
use more sophisticated algorithms able to separate objects which touch each other.
bwlabel
finds every connected set of pixels other than the background, and relabels these sets with a unique increasing integer. It can be called on a thresholded binary image in order to extract objects.
= bwlabel(logo)
logo_label table(logo_label)
logo_label
0 1 2 3 4 5 6 7
42217 1375 2012 934 1957 1135 1697 1063
The pixel values of the logo_label
image range from 0 corresponding to background to the number of objects it contains, which is given by
max(logo_label)
[1] 7
To display the image we normalize it to the (0,1) range expected by the display function. This results in different objects being rendered with a different shade of gray.
display( normalize(logo_label) )
The horizontal gray-scale gradient which can be observed reflects the way bwlabel()
scans the image and labels the connected sets: from left to right and from top to bottom. Another way of visualizing the segmentation is to use the colorLabels
function, which color codes the objects by a random permutation of unique colors.
display( colorLabels(logo_label) )
Some of the nuclei in nuc
are quite close to each other and get merged into one big object when thresholded, as seen in nuc_th
. bwlabel
would incorrectly identify them as a single object. The watershed transformation allows us to overcome this issue. The watershed
algorithm treats a gray-scale image as a topographic relief, or height-map. Objects that stand out of the background are identified and separated by flooding an inverted source image. In case of a binary image its distance map can serve as the input height-map. The distance map, which contains for each pixel the distance to the nearest background pixel, can be obtained by distmap
.
= watershed( distmap(nuc_th), 2 )
nmask display(colorLabels(nmask), all=TRUE)
Let’s try the watershed algorithm on structural MRI data.
library(brainR)
#load the sMRI data (same as in Problem 1)
# 3D sMRI data: http://socr.umich.edu/HTML5/BrainViewer/data/ABIDE_MRI_MPRAGE_peds_defaced.nii.gz
<- "http://socr.umich.edu/HTML5/BrainViewer/data/ABIDE_MRI_MPRAGE_peds_defaced.nii.gz"
brainURL <- file.path(tempdir(), "ABIDE_MRI_MPRAGE_peds_defaced.nii.gz")
brainFile download.file(brainURL, dest=brainFile, quiet=TRUE)
<- readNIfTI(brainFile, reorient=FALSE)
brainVolume
<- dim(brainVolume); brainVolDims brainVolDims
[1] 256 182 256
# get the mid-axial slice index
= brainVolDims[3] / 2; midZ midZ
[1] 128
# plot 2D sMRI slice
image(brainVolume[,, midZ], asp=1, col = hcl.colors(12, "YlOrRd", rev = TRUE)) # Hot-metal
# Apply watershed image segmentation.
= watershed( brainVolume[,, midZ], ext=2 )
sMRI_Watershed_seg display(normalize(sMRI_Watershed_seg), method="raster")
Voronoi tessellation is useful when we have a set of seed points (or regions) and want to partition the space that lies between these seeds in such a way that each point in the space is assigned to its closest seed. This function is implemented in the R package(“EBImage”) by the function propagate
. Let us illustrate the concept of Voronoi tessellation with a basic example. We use the nuclei mask nmask
as seeds and partition the space between them.
The 2D Interactive Voronoi Tessellation webapp provides a hands-on example of dynamic 2D Voronoi tessellation.
= propagate(seeds = nmask, x = nmask, lambda = 100)
voronoiExamp = colorLabels (voronoiExamp)
voronoiPaint display(voronoiPaint)
Only the first frame of the image stack is displayed.
To display all frames use 'all = TRUE'.
The basic definition of Voronoi tessellation, which we have given above, allows for two generalizations:
By default, the space that we partition is the full, rectangular image area, but indeed we could restrict ourselves to any arbitrary subspace. This is akin to finding the shortest distance from each point to the next seed not in a simple flat landscape, but in a landscape that is interspersed by lakes and rivers (which you cannot cross), so that all paths need to remain on the land. propagate
allows for this generalization through its mask
argument.
By default, we think of the space as flat – but in fact it could have hills and canyons, so that the distance between two points in the landscape not only depends on their x- and y-positions but also on the ascents and descents, up and down in z-direction, that lie in between. You can specify such a landscape to propagate
through its x
argument.
Mathematically, we can say that instead of the simple default case (a flat rectangle image with an Euclidean metric), we perform the Voronoi segmentation on a Riemann manifold, which can have an arbitrary shape and an arbitrary metric. Let us use the notation \(x\) and \(y\) for the column and row coordinates of the image, and \(z\) for the elevation of the landscape. For two neighboring points, defined by coordinates \((x, y, z)\) and \((x+dx, y+dy, z+dz)\), the distance between them is given by
\[ ds = \sqrt{ \frac{2}{\lambda+1} \left[ \lambda \left( dx^2 + dy^2 \right) + dz^2 \right] }. \] For \(\lambda=1\), this reduces to \(ds = ( dx^2 + dy^2 + dz^2)^{1/2}\). Distances between points further apart are obtained by summing \(ds\) along the shortest path between them. The parameter \(\lambda\ge0\) has been introduced as a convenient control of the relative weighting between sideways movement (along the \(x\) and \(y\) axes) and vertical movement. Intuitively, if you imagine yourself as a hiker in such a landscape, by choosing \(\lambda\) you can specify how much you are prepared to climb up and down to overcome a mountain, versus sideways walking around it.
When \(\lambda\) is large, the expression becomes equivalent to \(ds = \sqrt{dx^2 + dy^2}\), i. e., the importance of \(dz\) becomes negligible. This is what we did when we used lambda = 100
in our propagate
example.
A more advanced application of propagate
to the segmentation of cell bodies is presented in the “Cell segmentation example” section.
Let’s demonstrate Voronoi tessellation using the 2D sMRI image.
# Try Voronoi tessellation
= thresh(brainVolume[,, midZ], w=100, h=50, offset=0.05)
nmask = opening(nmask, makeBrush(5, shape='disc'))
nmask = fillHull(nmask)
nmask = bwlabel(nmask)
nmask display(nmask, all=TRUE)
= watershed( brainVolume[,, midZ], ext=2 )
sMRI_Watershed_seg <- sMRI_Watershed_seg * nmask
brain display(normalize(brain), method="raster")
= propagate(seeds = brain, x = brain, lambda = 100)
voronoi_sMRI = colorLabels (voronoi_sMRI)
voronoiPaint display(voronoiPaint, method="raster")
The R package(“EBImage”) defines an object mask as a set of pixels with the same unique integer value. Typically, images containing object masks are the result of segmentation functions such as bwlabel
, watershed
, or propagate
. Objects can be removed from such images by rmObject
, which deletes objects from the mask simply by setting their pixel values to 0. By default, after object removal all the remaining objects are relabeled so that the highest object ID corresponds to the number of objects in the mask. The reenumerate
argument can be used to change this behavior and to preserve original object IDs.
= list(
objects seq.int(from = 2, to = max(logo_label), by = 2),
seq.int(from = 1, to = max(logo_label), by = 2)
)= combine(logo_label, logo_label)
logos = rmObjects(logos, objects, reenumerate=FALSE)
z display(z, all=TRUE)
In the example above we demonstrate how the object removal function can be applied to a multi-frame image by providing a list of object indices to be removed from each frame. Additionally we have set reenumerate
to FALSE
keeping the original object IDs.
= function(image) lapply(getFrames(image), function(frame) unique(as.vector(frame)))
showIds
showIds(z)
[[1]]
[1] 0 1 3 5 7
[[2]]
[1] 0 2 4 6
Recall that 0 stands for the background. If at some stage we decide to relabel the objects, we can use for this the standalone function reenumarate
.
showIds( reenumerate(z) )
[[1]]
[1] 0 1 2 3 4
[[2]]
[1] 0 1 2 3
Holes in object masks can be filled using the function fillHull
.
= fillHull(logo)
filled_logo display(filled_logo)
floodFill
fills a region of an image with a specified color. The filling starts at the given point, and the filling region is expanded to a connected area in which the absolute difference in pixel intensities remains below tolerance
. The color specification uses R color names for Color
images, and numeric values for gray-scale
images.
= toRGB(logo)
rgblogo = rbind(c(50, 50), c(100, 50), c(150, 50))
points = c("red", "green", "blue")
colors = floodFill(rgblogo, points, colors)
rgblogo display( rgblogo )
display(floodFill(img, rbind(c(200, 300), c(444, 222)), col=0.2, tolerance=0.2))
Given an image containing object masks, the function paintObjects()
can be used to highlight the objects from the mask in the target image provided in the tgt
argument. Objects can be outlined and filled with colors of given opacities specified in the col
and opac
arguments, respectively. If the color specification is missing or equals NA
it is not painted.
= dim(img)[1:2]
d1 = Image(dim=d1)
overlay = dim(logo_label)-1
d2
= (d1-d2) %/% 2
offset
1]:(offset[1]+d2[1]), offset[2]:(offset[2]+d2[2])] = logo_label
overlay[offset[
= paintObjects(overlay, toRGB(img), col=c("red", "yellow"), opac=c(1, 0.3), thick=TRUE)
img_logo
display( img_logo )
In the example above we have created a new mask overlay
matching the size of our target image img
, and copied the mask containing the “EBImage” logo into that overlay mask. The output of paintObjects
retains the color mode of its target image, therefore in order to have the logo highlighted in color it was necessary to convert img
to an RGB image first, otherwise the result would be a gray-scale image. The thick
argument controls the object contour drawing: if set to FALSE
, only the inner one-pixel wide object boundary is marked; if set to TRUE
, also the outer boundary gets highlighted resulting in an increased two-pixel contour width.
We conclude our vignette by applying the functions described before to the task of segmenting cells. Our goal is to computationally identify and qualitatively characterize the cells in the sample fluorescent microscopy images. Even though this by itself may seem a modest goal, this approach can be applied to collections containing thousands of images, and that need is no longer a modest aim!
We start by loading the images of nuclei and cell bodies. To visualize the cells we overlay these images as the green and the blue channel of a false-color image.
= readImage(system.file('images', 'nuclei.tif', package='EBImage'))
nuc = readImage(system.file('images', 'cells.tif', package='EBImage'))
cel
= rgbImage(green=1.5*cel, blue=nuc)
cells display(cells, all = TRUE)
First, we segment the nuclei using thresh
, fillHull
, bwlabel
and opening
.
= thresh(nuc, w=10, h=10, offset=0.05)
nmask = opening(nmask, makeBrush(5, shape='disc'))
nmask = fillHull(nmask)
nmask = bwlabel(nmask)
nmask
display(nmask, all=TRUE)
Next, we use the segmented nuclei as seeds in the Voronoi segmentation of the cytoplasm.
= opening(cel>0.1, makeBrush(5, shape='disc'))
ctmask = propagate(cel, seeds=nmask, mask=ctmask)
cmask
display(ctmask, all=TRUE)
To visualize our segmentation we use paintObject
.
= paintObjects(cmask, cells, col='#ff00ff')
segmented = paintObjects(nmask, segmented, col='#ffff00')
segmented
display(segmented, all=TRUE)
Let’s use the protocol below to gain some hands-on experience processing a 2D brain image using basic functions, without relying on the EBImage package. Of course, you can also apply all of the methods presented above on this 2D (JPEG) brain image.
# Install and load the JPEG package
# install.packages("jpeg")
library(jpeg)
# Read Image in and display it
# Either download the JPG image locally from https://umich.instructure.com/courses/38100/files/1627149/download?download_frd=1 and then read it in
# mri <- readJPEG("C:/Users/Dinov/Desktop/MRI_ImageHematoma.jpg")
# Or directly load it from the web
<- tempfile(, fileext = ".jpg")
pathToMRI download.file("https://umich.instructure.com/courses/38100/files/1627149/download?download_frd=1", mode = "wb", pathToMRI)
<- readJPEG(pathToMRI)
mri # display(mri)
plot.new(); rasterImage(mri,0,0,1,1)
# Threshold the image
# display( thresh(mri, offset=0.05), all=TRUE )
# plot.new(); rasterImage(EBImage::thresh(mri, offset=0.05),0,0,1,1)
= mri[,,1] > 0.35
mri_thresh = mri[,,1] * mri_thresh
mri_thresh plot.new(); rasterImage(mri_thresh,0,0,1,1)
# Histogram
hist(mri)
# Affine transformation
= matrix(c(1, -.5, 128, 0, 1, 0), nrow=3, ncol=2)
m = EBImage::affine(mri, m)
img_affine # display( img_affine )
plot.new(); rasterImage(img_affine,0,0,1,1)
# Sheer Affine Transformation and translation
= matrix(c(1, -.5, 128, 0, 1, 0), nrow=3, ncol=2)
m = EBImage::affine(EBImage::translate(mri, c(-50, 0)), m)
img_affine # display( img_affine )
plot.new(); rasterImage(img_affine,0,0,1,1)
# Gaussian Kernel smoothing (Convolution)
= EBImage::makeBrush(size = 31, shape = 'gaussian', sigma = 5)
kernel = EBImage::filter2(mri, kernel)
mri_smooth # display(mri_smooth)
plot.new(); rasterImage(mri_smooth,0,0,1,1)
# FFT
<- fft(mri[,,1]) # fftw2d # Display Re(FT):
ft_mriMat ::display(fftshift(ft_mriMat)) EBImage
# FFT Magnitude
<- sqrt(Re(ft_mriMat)^2+Im(ft_mriMat)^2)
mag_ft_mriMat display(mag_ft_mriMat)
#plot.new(); rasterImage(normalize(mag_ft_mriMat),0,0,1,1)
# IFT
<- atan2(Im(ft_mriMat), Re(ft_mriMat)); Real = mag_ft_mriMat * cos(phase_ft_mriMat)
phase_ft_mriMat = mag_ft_mriMat * sin(phase_ft_mriMat)
Imaginary = Re(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_mriMat));
ift_ft_mag_ft_mriMat display(t(ift_ft_mag_ft_mriMat))
# Low- and High-pass Fourier domain filtering
## Low-pass frequency filtering
## Remove Fourier coefficients indexed above $k_{x,max/2}$ and $k_{y,max/2}$.
<- mag_ft_mriMat
mag_ft_mriMat_1 for (i in 1:dim(mag_ft_mriMat)[1]) {
for (j in 1:dim(mag_ft_mriMat)[2]) {
if (abs(i-(dim(mag_ft_mriMat)[1])/2) > (dim(mag_ft_mriMat)[1])/4 |
abs(j-(dim(mag_ft_mriMat)[2])/2) > (dim(mag_ft_mriMat)[2])/4) {
<- 0
mag_ft_mriMat_1[i,j]
}
}
}
# EBImage::display(mag_ft_mriMat_1, method = "raster", title="(IFT o LowPassFilter o FT) Magnitude=Img | Phase=Img")
plot.new(); rasterImage(normalize(mag_ft_mriMat_1),0,0,1,1)
= mag_ft_mriMat_1 * cos(phase_ft_mriMat)
Real = mag_ft_mriMat_1 * sin(phase_ft_mriMat)
Imaginary = Re(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_mriMat_1))
ift_ft_mriMat_1 ::display(t(400*ift_ft_mriMat_1), method = "raster", title="(IFT o LowPassFilter o FT) Magnitude=MRI | Phase=MRI") EBImage
# plot.new(); rasterImage(normalize(ift_ft_mriMat_1),0,0,1,1)
## High-pass frequency filtering
## Remove the Fourier coefficients indexed below $k_{x,max/4}$ and $k_{y,max/4}$.
# High-pass frequency filtering, remove Fourier coeff's indexed below k_{x,max/4} and k_{y,max/4}
<- mag_ft_mriMat
mag_ft_mriMat_2 for (i in 1:dim(mag_ft_mriMat)[1]) {
for (j in 1:dim(mag_ft_mriMat)[2]) {
if (abs(i-(dim(mag_ft_mriMat)[1])/2) < (dim(mag_ft_mriMat)[1])/4 &
abs(j-(dim(mag_ft_mriMat)[2])/2) < (dim(mag_ft_mriMat)[2])/4) {
<- 0
mag_ft_mriMat_2[i,j]
}
}
}
::display(mag_ft_mriMat_2, method = "raster", title="(IFT o HighPassFilter o FT) Magnitude=Img | Phase=Img") EBImage
# plot.new(); rasterImage(normalize(mag_ft_mriMat_2),0,0,1,1)
= mag_ft_mriMat_2 * cos(phase_ft_mriMat)
Real = mag_ft_mriMat_2 * sin(phase_ft_mriMat)
Imaginary = Re(fft(Real+1i*Imaginary, inverse = T)/length(mag_ft_mriMat_2))
ift_ft_mriMat_2 # EBImage::display(t(ift_ft_mriMat_2), method = "raster", title="(IFT o HighPassFilter o FT) Magnitude=Img | Phase=Img")
plot.new(); rasterImage(normalize(ift_ft_mriMat_2),0,0,1,1)
Let’s try to interrogate some artworks from UMMA’s 2020 Curriculum-Collection. This museum collection exhibits a variety of University of Michigan courses that use artwork to enhance the curricula of graduate and undergraduate courses. Each course examines a set of artworks and objects to address the nature of reality, imagination, and vision relative to the current environments and scientific advances. The art exhibit demonstrates diverse and creative examples of learning across the art, science, technology, and biosocial disciplines.
For instance, the painting by Khaled al-Saa’i, “Winter in Ann Arbor”, constructed with natural ink, tempera and gouache on paper, may reveal interesting patterns, shapes, topological structure, and calligraphy. We can interrogate the artwork using the following image-processing protocol:
Further studies may include text recognition, identification of obfuscated structures, or graphical representation of 2D scenes. The readers are encouraged to expand this protocol and apply it to some of the other artworks in these UMMA collection 1 and collection 2.
# UMMA HS650 CC: https://umma.umich.edu/curriculum-collection/health-sciences-650
library(EBImage)
# Load the image from CC
# CC_WhitePoemFigure_Dill <- readImage("https://exchange.umma.umich.edu/media/W1siZiIsIjIwMjEvMDIvMTkvZGJ3dWk2cTVjX2RlZmF1bHQuanBnIl0sWyJwIiwidGh1bWIiLCIxMDAweDEwMDAiXV0?sha=bf87973a898145f8", type="jpeg")
<- readImage("https://exchange.umma.umich.edu/media/W1siZiIsIjIwMjEvMDkvMDMvNXZjaXdkOWdpal9kZWZhdWx0LmpwZyJdLFsicCIsInRodW1iIiwiMTAwMHgxMDAwIl1d?sha=9b8e04fa0a73cbb8", type="jpeg")
CC_WhitePoemFigure_Dill
# Display raw image
display(CC_WhitePoemFigure_Dill, method = "raster", all = TRUE)
# Intensity Histogram
hist(CC_WhitePoemFigure_Dill)
# Generate negative image
<- max(CC_WhitePoemFigure_Dill) - CC_WhitePoemFigure_Dill
img_neg display( img_neg )
# dim: 633 1000 3
# crop and threshold
<- CC_WhitePoemFigure_Dill # [200:578, 250:750, ]
img_crop <- img_crop > 0.55
img_thresh display(img_thresh)
# Low-pass Filter/smooth
<- makeBrush(size = 31, shape = 'gaussian', sigma = 5)
w <- filter2(CC_WhitePoemFigure_Dill, w)
img_flo display(img_flo)
# High-pass filtering Laplace Filter/Edge detection
<- matrix(1, nrow = 3, ncol = 3)
fhi 2, 2] <- -7
fhi[<- filter2(CC_WhitePoemFigure_Dill, fhi)
img_fhi display(img_fhi)
# Denoising
<- medianFilter(CC_WhitePoemFigure_Dill, size=5, cacheSize=20000)
img_median display(img_median)
# Segmentation
<- channel(CC_WhitePoemFigure_Dill,"gray")
grayimage <- otsu(grayimage)
threshold <- combine( mapply(function(frame, th) frame > th,
CC_th getFrames(CC_WhitePoemFigure_Dill),
SIMPLIFY=FALSE) )
threshold, display(CC_th, all=TRUE)
<- makeBrush(51, "disc")
disc <- disc/sum(disc)
disc <- 0.0001
offset <- filter2(CC_WhitePoemFigure_Dill, disc)
CC_bg <- CC_WhitePoemFigure_Dill < CC_bg + offset
CC_th display(CC_th, all=TRUE)
# Watershed segmentation
<- watershed(distmap(CC_th), tolerance=5, ext=20)
nmask display(colorLabels(nmask), all=TRUE)
# Voronoi image tessellation
<- propagate(x = nmask, seeds = nmask, lambda = 0.001)
voronoiExamp <- colorLabels (voronoiExamp)
voronoiPaint display(voronoiPaint, all=TRUE)
<- rgbImage(green=2.0*CC_bg, blue=img_median, red=voronoiPaint[ , , , 1])
CC display(CC, all = TRUE)
= opening(img_fhi<0.6, makeBrush(5, shape='disc'))
ctmask = propagate(img_fhi, seeds=nmask, mask=ctmask)
cmask display(ctmask, all=TRUE)
= opening(img_fhi>0.1, makeBrush(5, shape='disc'))
ctmask = propagate(img_fhi, seeds=nmask, mask=ctmask)
cmask display(ctmask, all=TRUE)
# segmented <- paintObjects(channel(CC_th,"gray"), CC_WhitePoemFigure_Dill, col='#ff00ff')
# segmented <- paintObjects(channel(nmask,"gray"), segmented, col='#ffff00')
# display(segmented, all=TRUE)
= paintObjects(channel(cmask,"gray"), CC_WhitePoemFigure_Dill, col='#ff00ff')
segmented = paintObjects(channel(ctmask,"gray"), segmented, col='#ffff00')
segmented display(segmented, all=TRUE)