——标题:“图像数据和空间模式分析的基础知识在*R*”作者:“Andrzej oleovic和Wolfgang Huber”日期:7月21日| BioC2015输出:ioslides_presentation: fig_retina: null css: style.css宽屏幕:假书目:BioC2015Oles。bib nocite: | @Pau2010， @Haralick1979， @Jones2005 vignette: > %\VignetteIndexEntry{R} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8}——' ' ' {r回声= FALSE,消息= FALSE}库(BioC2015Oles)库(knitr)美元opts_chunk组(错误= FALSE) set.seed (7) .dpi = 100 ' ' ' # #这个工作坊的目标——学习如何读和写图片的* R * - * R *学习图像表示,如何操作,理解如何应用过滤器和转换图像-应用这些技能细胞显微图像分割和特征提取,探索空间分布的位置细胞# # EBImage # # #的形象*R*的处理和分析工具箱。图像文件的读取和写入。交互式图像查看器。图像操作，转换和过滤。对象检测和特征提取
自*Bioconductor* 1.8 (2006)
最初开发者:Oleg Sk2021欧洲杯体育投注开户lyar Wolfgang Huber Mike Smith Gregoire Pau 贡献者:Joseph Barry Bernd Fischer Ilia Kats Philip A. Marais ! [EBImage-logo](图片/ logo.png)
让我们开始吧!{r, message=FALSE, fig.width=768/。dpi, fig.height = 512 /。dpi, dpi =。dpi/2}库(EBImage) f =系统。file("images"， "sample.png"， package="EBImage") img = readImage(f) display(img)##读取和显示图像
图片可以从本地文件或url中读取。' ' ' {r, fig.width = 480 /。dpi, fig.height = 138 /。dpi, dpi =。dpi, eval=FALSE} bioc = readImage("//www.andersvercelli.com/images/logo/jpg/bioconductor_logo_rgb.jpg") display(bioc)' ' '
! [RBioFormats](图片/ RBioFormats-logo.png) *EBImage*支持JPEG, PNG和TIFF文件格式。要读取专有显微镜图像数据和元数据，请使用*[RBioFormats](https://github.com/aoles/RBioFormats)*。
###显示图像-交互式JavaScript查看器- *R*的内置绘图设备默认的' display '方法可以通过' options(EBImage.display) '设置。' ' ' {r, eval=FALSE}选项(EBImage. eval=FALSE)。Display = "光栅")“# #添加文本标签的“{r, fig.width =暗(img) [1 l] /。dpi fig.height =暗(img) [2 l] /。dpi, dpi =。dpi/2, results='hide'} display(img, method = " rster ") text(x = 20, y = 20, label = "Parrots"， adj = c(0,1)， col = "orange"， cex = 2) filename = "parrot .jpg" dev.print(jpeg, filename = filename, width = dim(img)[1]， height = dim(img)[2])' ' ' {r} file.size(文件名)支持的文件格式:JPEG, PNG和TIFF。{r} writeImage(img， "sample.jpeg"， quality = 85) writeImage(img， "sample.tiff") writeImage(img， "sample_compressed.tiff"， compression = "deflate") files = list。data.frame(row.names= Files, size=file.size(Files))##图像表示
多维像素强度数组- (x, y) - (x, y， **z**) z堆栈- (x, y， **t**)延时- (x, y， **c**)通道- (x, y, c, z, t，…) (样本)!(图片/ sample.png)
!(复制)(图片/ replicates.png) ! [sample-rgb](图片/ sample-rgb.png)
##图像表示' ' ' {r} str(img) getClassDef("Image") dim(img)##图像摘要{r} img imageData(img)[1:3, 1:6] ' ' ' ##图像直方图' ' ' {r, fig.width=4, fig.height=4} hist(img) range(img)“# #彩色图像的“{r, fig.width =暗(img) [1 l] /。dpi fig.height =暗(img) [2 l] /。dpi, dpi =。Dpi /2，} f = system。file("images"， "sample-color.png"， package="EBImage") imgcol = readImage(f) display(imgcol) print(imgcol, short = TRUE)' ' ' ##图像堆栈' ' ' {r, echo=FALSE} nuc = readImage(system。file("images"， "nuclei.tif"， package="EBImage"))' ' ' {r, fig.width =暗(nuc) [1 l] /。dpi fig.height =暗(nuc) [2 l] /。dpi, dpi =。dpi/2} nuc = readImage(system. dpi;file("images"， "nuclei.tif"， package="EBImage")) print(nuc, short = TRUE) display(nuc)““# #图像栈”“{r, fig.width =暗(nuc) [1 l] /。dpi fig.height =暗(nuc) [2 l] /。dpi, dpi =。dpi} display(nuc, method = "raster", all = TRUE) ``` ## Manipulating images Being numeric arrays, images can be conveniently manipulated by any of *R*'s arithmetic operators. ### Cropping ```{r, fig.width=384L/.dpi, fig.height=384L/.dpi, dpi=.dpi/2} img = img[366:749, 58:441] ``` ### Inversion ```{r negative, fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2} img_neg = max(img) - img display(img_neg) ``` ## Manipulating images ### Brightness, contrast, and gamma correction ```{r arithmetic, fig.width=(4*dim(img)[1L]+100)/.dpi, fig.height=(dim(img)[2L]+40)/.dpi, dpi=.dpi/2} img_comb = combine( img, img + 0.3, img * 2, img ^ 0.5 ) display( tile(img_comb, numberOfFrames(img_comb), lwd = 20, fg.col = "white") ) ``` ## Tresholding images ### Binary images Images which contain only two sets of pixels, which represent the background and the foreground pixels. ```{r, fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2} img_thresh = img > .5 display(img_thresh) ``` ## Tresholding images ### Binary images Images which contain only two sets of pixels, which represent the background and the foreground pixels. ```{r} img_thresh ``` ## Spatial transformations ### Translation ```{r translate, fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2} img_translate = translate(img, v = c(100, -50)) display(img_translate) ``` ## Spatial transformations ### Rotation ```{r rotate-pre, echo=FALSE} img_rotate = rotate(img, 30) ```{r rotate, fig.width=dim(img_rotate)[1L]/.dpi, fig.height=dim(img_rotate)[2L]/.dpi, dpi=.dpi/2} img_rotate = rotate(img, angle = 30, bg.col = "white") display(img_rotate) ``` ## Spatial transformations ### Scaling ```{r resize, fig.width=512/.dpi, fig.height=256/.dpi, dpi=.dpi/2} img_resize = resize(img, w = 512, h = 256) display(img_resize) ```{r resize2, fig.width=256/.dpi, fig.height=256/.dpi, dpi=.dpi/2} img_resize = resize(img, 256) display(img_resize) ``` ## Spatial transformations ### Vertical and horizontal reflection ```{r flipflop, fig.width=dim(img)[1L]/.dpi, fig.height=dim(img)[2L]/.dpi, dpi=.dpi/2} display( flip(img) ) display( flop(img) ) ``` ## Affine transformation Described by a 3x2 transformation matrix $\textbf{m}$
$ $ \ {bmatrix}开始x ' _1 & y ' _1 \ \ ' _2 & y ' _2 \ \ \ vdots & \ vdots \ \ \ {bmatrix} = {bmatrix} \开始结束x_1 & y_1 & 1 \ \ x_2 & y_2 & 1 \ \ \ vdots & \ vdots & \ ddots \ \ \ {bmatrix}结束\ * \ textbf {m} $ $
$ $ \ textbf {m} _{\文本翻译}{}= {bmatrix} \开始1 0 \ \ & 0 & 1 \ \ t_x & t_y \ \ \ {bmatrix}, \四\ textbf {m} _{\文本{旋转}}= {bmatrix} \ \开始cosθ}{\ & \ sinθ}{\ \ \ - \ sinθ}{\ & \ cosθ}{\ \ \ & 0 \ \ \ {bmatrix}, \四\ textbf {m} _{\文本{缩放}}= {bmatrix}开始\ W & 0 \ \ 0 & H \ \ & 0 \ \ \ {bmatrix} $ $
仿射变换例子:水平垂直映射' ' ' {r, fig.width=dim(img)[1L]/。dpi fig.height =暗(img) [2 l] /。dpi, dpi =。Dpi /2}角度= pi/6 m = matrix(c(1， -sin(角度)，sin(角度)*dim(img)[2]/ 2,0,1,0)， nrow = 3, ncol = 2) m img_affine = affine(img, m) display(img_affine)' ' ' ##图像转置' ' ' {r转置，fig.width=dim(imgcol)[2L]/。dpi fig.height =暗(imgcol) [1 l] /。dpi, dpi =。Dpi /2} imgcol_t =转置(imgcol)显示(imgcol_t)线性滤波器
移除局部工件或噪声—矩形窗口中的平均值
$$ f'(x,y) = \frac{1}{N} \sum_{s=-a}^{a}\sum_{t=-a}^{a} f(x+s, y+t)，$$其中$f(x,y)$是位置$(x, y)$的像素值，$ a$决定窗口大小，每个方向为$2a+1$。$N=(2a+1)^2$是平均的像素数，$f'$是新的平滑图像。
泛化:使用权重函数$w$进行加权平均
$ $ (w * f) (x, y) = \ sum_ {s = - \ infty} ^ {+ \ infty} \ sum_ {t = - \ infty} ^ {+ \ infty} w (s, t) \ f (x + s y + s) $ $ -卷积图像f和w美元,美元由$ * $ -线性表示:对于任意两个图像f,美元₂和数字c₁美元,美元$ c₂$ $ $ w * (c_1f_1 + c_2f_2) = c_1w * f + c_2w *₂$ $
权重函数可以使用“makeBrush”生成。{r, echo=FALSE} opts_knit$set(global.par=TRUE)' ' ' {r, echo=FALSE} .par=par(mai = c(。45， .75, 0.05, 0.05)){r makeBrush, fig.width=3.8, fig.height=3.5, dev="svg"，} w = makeBrush(size = 31, shape = 'gaussian'， sigma = 5) plot(w[(nrow(w)+1)/2，]， ylab =" w"， xlab =" ")其他可用的笔刷形状:“盒子”(默认)，“圆盘”，“菱形”，“线”##低通滤波二维线性卷积由“filter2”实现(使用FFT)例子:使用宽度为5的高斯滤波器平滑图像' ' ' {r, echo=FALSE} opts_knit$set(global.par=FALSE)' ' ' {r lopass, fig.width=dim(img)[1L]/。dpi fig.height =暗(img) [2 l] /。dpi, dpi =。Dpi /2} img_flo = filter2(img, w) display(img_flo)' ' ' ' {r高通，fig.width=dim(img)[1L]/。dpi fig.height =暗(img) [2 l] /。dpi, dpi =。Dpi /2} fhi = matrix(1, nrow = 3, ncol = 3) fhi[2,2] = -8 fhi img_fhi = filter2(img, fhi) display(img_fhi)补偿背景信号的空间依赖性' ' ' {r, fig.width=(3*dim(nuc)[1L]+80)/。dpi, fig.height =(暗(nuc) (2 l) + 40) /。dpi, dpi =。dpi/2} disc = makeBrush(21， "disc") disc = disc / sum(disc) nuc = getFrame(nuc, 1) nuc_bg = filter2(nuc, disc) offset = 0.02 nuc_thresh = (nuc - nuc_bg) > offset img_comb = combine(nuc, nuc_bg, nuc_thresh) display(tile(img_comb, numberOfFrames(img_comb)， lwd = 20, fg.;Col = "白色"))非线性噪声去除技术-在散斑噪声的情况下特别有效-保留边缘' ' ' {r medianFilter, fig.width=2*dim(img)[1L]/。dpi fig.height =暗(img) [2 l] /。dpi, dpi =。dpi/2, eval=TRUE} l= length(img) n = l/10 img_noise = img img_noise [sample(l, n)] = runif(n, min = 0, max = 1) img_median = medianFilter(img_noise, size = 1) display(combine(img_noise, img_median)， all=TRUE)' ' '
使用常数时间算法实现[@Perreault2007]。
形态运算{。二值图像的非线性滤波-侵蚀:对于每一个成品像素，把蒙版放在它周围，如果蒙版下的任何像素来自bg，将这些像素设置为bg．-扩张:每bg像素，把蒙版放在它周围，如果蒙版下的任何像素来自成品，将这些像素设置为成品．' ' ' {r echo=FALSE} leaf = readImage(系统。file('images'， 'leaf.png'， package='BioC2015Oles')) .leafDim = BioC2015Oles:::tilesDim(dim(leaf)， 3)' ' ' {r, fig.width = .leafDim[1] /。dpi fig.height = .leafDim[2] /。dpi, dpi =。dpi out.width = 2 *。leafDim[1]， results='hide'}file('images'， 'leaf.png'， package='BioC2015Oles')) kern = makeBrush(size = 3) morph = combine(leaf, erosion (leaf, kern)， dilate(leaf, kern)) displayTiles(morph)' ' ' ##形态操作{。Noiterpolation} - **开放**侵蚀之后的膨胀从背景中移除小物体- **关闭:**膨胀之后的侵蚀填补了前景中的空洞' ' ' {r, fig.width=. leafdim[1]/。dpi fig.height = .leafDim[2] /。dpi, dpi =。dpi out.width = 2 *。leafDim[1]} morph = combine(叶子，打开(叶子，kern)，关闭(叶子，kern)) displayTiles(morph)细胞生物学中的应用HeLa细胞的两个通道荧光显微镜图像。' ' ' {r} dna = readImage(system。file("images"， "nuclei.tif"， package="EBImage")) print(dna, short=TRUE) tub = readImage(system. image . txt)file("images"， "cells.tif"， package="EBImage")) print(tub, short=TRUE)##在细胞生物学中的应用{r, fig.width=2*dim(nuc)[1L]/。dpi fig.height =暗(nuc) [2 l] /。dpi, dpi =。dpi/3} display(combine(getFrame(dna, 3)， getFrame(tub, 3))， all=TRUE){r，图.width=BioC2015Oles:::tilesDim(dim(dna)， numberOfFrames(dna))[1L]/。dpi, fig.height=BioC2015Oles:::tilesDim(dim(dna)， numberOfFrames(dna))[2L]/。dpi, dpi =。dpi/3} rgb = rgbImage(绿色= 1.5 * tub，蓝色= dna) displayTiles(rgb) ``` ## Nuclei segmentation - Adaptive thresholding of the DNA channel - Morphological opening and filling of holes - Distance map computation and watershed transformation ### Adaptive thresholding and morphological filtering ```{r, fig.width=2*dim(nuc)[1L]/.dpi, fig.height=dim(nuc)[2L]/.dpi, dpi=.dpi/2} nmaskt = thresh(dna, w = 15, h = 15, offset = 0.05) nmaskf = fillHull( opening(nmaskt, makeBrush(5, shape='disc')) ) display( combine(getFrame(nmaskt, 3), getFrame(nmaskf, 3)), all=TRUE ) ``` ## Nuclei segmentation ### Distance map computation ```{r, fig.width=dim(nuc)[1L]/.dpi, fig.height=dim(nuc)[2L]/.dpi, dpi=.dpi/2} dmap = distmap(nmaskf) range(dmap) display(normalize(dmap), frame = 3) ``` ## Nuclei segmentation ### Watershed transformation ```{r, fig.width=2*dim(nuc)[1L]/.dpi, fig.height=dim(nuc)[2L]/.dpi, dpi=.dpi/2} nmask = watershed(dmap, tolerance = 2) display( combine( toRGB( getFrame(dna, 3) ), colorLabels( getFrame(nmask, 3) ) ), all=TRUE ) ``` `bwlabel` - a simpler and faster function for segmenting connected objects ## Cytoplasm segmentation Identification of cell bodies by Voronoi partitioning using the segmented nuclei as seeds. ```{r, fig.width=2*dim(nuc)[1L]/.dpi, fig.height=dim(nuc)[2L]/.dpi, dpi=.dpi/2} cmaskt = closing( gblur(tub, 1) > 0.105, makeBrush(5, shape='disc') ) cmask = propagate(tub, seeds=nmask, mask=cmaskt, lambda = 0.001) display( combine( toRGB( getFrame(cmaskt, 3) ), colorLabels( getFrame(cmask, 3) ) ), all=TRUE ) ``` ## Voronoi tesselation Voronoi diagram is a partitioning of a plane into regions based on distance to some specific points of that plane (seeds). ### Voronoi tesselation on a Riemann manifold - arbitrary shape (provided in `mask`) - arbitrary metric (dependent on image gradient $z$)
点$(x, y, z)$和点$(x+dx, y+dy, z+dz)$之间的距离大概{$ $ ds = \ \压裂{2}{\λ+ 1}\离开[\λ\离开(dx ^ 2 + dy ^ 2 \右)+ dz ^ 2 \]} $ $$\lambda$控制了欧几里得距离项的权重——当$\lambda$较大时，$ds$趋于欧几里得距离——对于较小的$\lambda$， $ds$趋于图像的强度梯度
##最终分割' ' ' {r, fig.width=2*dim(nuc)[1L]/。dpi fig.height = 2 *暗(nuc) [2 l] /。dpi, dpi =。dpi/2} display(paintObjects(nmask, paintObjects(cmask, rgb, col = "magenta"， thick = TRUE)， col = "yellow"， thick = TRUE)， all = TRUE){r, fig.width=2*dim(nuc)[1L]/. {r, fig.width=2*dim(nuc)[1L]/。dpi fig.height = 2 *暗(nuc) [2 l] /。dpi, dpi =。dpi/2, eval=TRUE} st = stackObjects(cmask, rgb) display(st, all =TRUE)定量细胞描述符-像素强度的基本统计。基本形状特征(面积、周长、半径)“computeFeatures。形状-矩(质心，偏心率，…)“computeFeatures。moment ' - Haralick纹理特征' computeFeatures。haralick` ```{r} head(computeFeatures。Shape (cmask[，，1]， tub[，，1]))##自动化细胞表型使用多元分析方法，我们可以:-检测细胞亚群(聚类)-将细胞分类为预定义的细胞类型或表型(分类)-基于亚群的频率比较不同的生物条件
!(管道)(图片/ pipeline.png)
淋巴结-由B细胞、T细胞、树突状细胞(DC)和巨噬细胞组成的免疫过滤器-癌症分期的临床意义(*前哨*淋巴结)
' ' ' {r, echo=FALSE} .淋巴结= readImage(系统。file("images"， "testscan_1_2_RGB99-4525D.jpg"， package = "BioC2015Oles"))' ' ' {r,呼应= FALSE, fig.width =暗(.lymphnode) [1 l] / (4 * .dpi), fig.height =暗(.lymphnode) (2 l) / (4 * .dpi), dpi =。dpi}显示器(.lymphnode)' ' '
乳腺癌淋巴结活检[@Setiadi2010]。细胞分型是通过提供特定特征的蛋白质抗体染色细胞来完成的。
##淋巴结活检' ' {r} data(brcalymphnode, package = "BioC2015Oles") head(brcalymphnode) nrow(brcalymphnode) table(brcalymphnode$class)细胞空间分布T细胞(左)和肿瘤细胞(右)的*x*和*y*位置图{r, fig.width=7, fig.height=3.5} par(mfrow = c(1,2)， mar = c(4.5, 4.5, 0.5, 0.5)) xlim = range(brcalymphnode$x) ylim = range(brcalymphnode$y) cols = c(' T_cells ' = "dodgerblue4"， ' Tumor ' = "darkmagenta") for(i in seq_along(cols)) plot(子集(brcalymphnode, class==names(cols)[i])[， c("x"， "y")]， pch = "."， asp = 1, xlim = xlim, ylim = ylim, col = cols[i])标记的空间点过程:-位于数学空间中的孤立点集(' xlim '， ' ylim ') -标记有特定属性的点(因子'类')’‘{r、消息= FALSE}图书馆(“spatstat”)ln = (brcalymphnode,购买力平价(x = x, y = y,标志着=类,xrange = xlim yrange = ylim)) ln的“# #凸包‘’{r, fig.width = 2, fig.height = 2} chull = convexhull (ln)标准(mar = c(0, 0, 0, 0))情节(窗口(ln),主要= NULL, asp = 1)情节(chull lty = 2,坳=“lightblue”,添加= TRUE)窗口(ln) = chull ln ' ' ' # #一阶效应:强度——点躺在一个圆的面积一美元在点p = (x, y)美元$ $ $ $ N (p)-局部强度$ $ \λ(p) = \ lim_ {\ rightarrow 0} \压裂{E [N (p)]}{一}$ $-对于*平稳*过程，即在整个区域均匀$$\lambda(p) = \text{const.}$$那么区域A的强度与面积成正比$$E(N(\cdot, A)) = \lambda A$$核平滑的点模式强度估计[@Diggle1985]
$$\hat{\lambda}(p) = \sum_i e(p_i) K(p-p_i)$$，其中$K$为高斯平滑核，$e(p_i)$为边缘修正因子。
“‘{r密度。ppp1, fig.width=5.25, fig.height=2.65}密度= solist(' Diggle's edge correction ' =密度(子集(ln, marks=="Tumor")， Diggle = TRUE)， ' No edge correction ' =密度(子集(ln, marks=="Tumor")， edge = FALSE)) plot(密度，等于。ribbon = TRUE, col = topo。颜色，main = ""){r, fig.width=9.5, fig.height=2.8} rr = relrisk(ln, sigma = 250) plot(rr，等于。ribbon = TRUE, col = topo。colors, nrows = 1, main = "") ``` ## Second order effects: clustering Spatial clustering (or anticlustering) of points: whether and to what extent neighboring points are more/less dense than expected by chance. ### Poisson process
-稳定的强度$\lambda$ -在非重叠区域中物体的出现之间没有依赖关系- $N(p, A)$遵循速率$\lambda A$的泊松分布
###最近邻居(NN)距离分布函数$G$从一个典型的随机点到它的最近邻居的距离$W$的累积分布函数。对于泊松过程$$G(r) = P(W \leq r) = 1-e^{-\ λ \pi r^2}$$## NN距离函数的估计*G* ' ' ' {r Gest} gln = Gest(ln) gln ' ' ' ## NN距离函数的估计*G* ' ' ' ' {r, fig.width=4, fig.height=4, message=FALSE} library(RColorBrewer) par(mar = c(4.5, 4.5, 0.5, 0.5)) plot(gln, lty = 1, col = brewer。pal(4， "Set1")， main = "")' ' '
-有限细胞大小的影响-大距离聚类
Ripley的*K*函数对于一个静止点过程$\ λ K(r)$是在一个给定的，随机选择的点$r$的距离内的期望附加点的数量。推广到非齐次点过程
$ $ K_ {inhom}}{\文本(r) = \ sum_ {i, j} {\ mathbf 1} _ {d le r} (p_i p_j) \ \ * \压裂{e (p_i p_j, r)}{\λ(x_i) \λ(x_j)}, $ $
其中$d(p_i, p_j)$是点$p_i$和点$p_j$之间的距离，$e(p_i, p_j, r)$是边修正因子。对于估计和可视化更有用的是$L$-函数
大概{$ $ L (r) = \ \压裂{K (r)}{\π}}。$ $
对于泊松点模式，理论值为$L(r) = r$。##非齐次的估计*L*-function ' ' ' {r, eval=FALSE} Lln = Linhom(子集(ln, marks=="T_cells")) Lln ' ' ' {r, package =" BioC2015Oles") Lln ' ' ' ##非齐次的估计*L*-function ' ' ' {r, fig.width=4.5, fig.height=4.5} par(mar = c(4.5, 4.5, 0.5, 0.5)) plot(Lln, lty = 1, col = brewer。pal(3， "Set1")， main = "")' ' ' ' ' ' ' ' ' '对相关函数描述点密度如何作为距离参考点的函数而变化。
$ $ g (r) = \压裂{1}{2 \πr} \压裂{dK}{}博士(r) $ $
对于平稳泊松过程$g(r) = 1$。
- $g(r) < 1$——点间抑制- $g(r) > 1$——聚类
##对相关函数' ' ' {r, eval=FALSE} pcfln = pcf(Kinhom(子集(ln, marks=="T_cells"))) plot(pcfln, lty = 1, log =" x"){r, fig.width=4.5, fig.height=4.5, echo=FALSE} par(mar = c(4.5, 4.5, 0.5, 0.5)) data(pcfln, package = "BioC2015Oles") plot(pcfln, lty = 1, log = "x"， main = "")##引用{.small}