# em.saltykov: Expectation Maximization (EM) algorithm In unfoldr: Stereological Unfolding for Spheroidal Particles

## Description

Estimation of the empirical sphere diameter distribution

## Usage

 `1` ```em.saltykov(y, bin, maxIt = 32) ```

## Arguments

 `y` vector of observed absolute frequencies of circle diameters `bin` non-decreasing vector of class limits `maxIt` maximum number of iterations to be used

## Details

The function performs the EM algorithm, see reference below.

## Value

vector of count data of absolute frequenties of sphere diameters

## References

Ohser, J. and Muecklich, F. Statistical analysis of microstructures in materials science J. Wiley & Sons, 2000

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216``` ```## Not run: ## Comment: Simulate a Poisson sphere system, ## intersect, discretize and display results library(unfoldr) library(rgl) library(plotrix) spheres <- function(spheres, box=NULL, draw.box=FALSE, draw.axis=FALSE, ...) { xyz <- apply(sapply(spheres, "[[", "center"),1,function(x) x) sizes <- unlist(lapply(spheres,function(x) x\$r)) spheres3d(xyz,radius=sizes,...) if(draw.axis) { axes3d(c('x','y','z'), pos=c(1,1,1), tick=FALSE) #title3d('','','x','y','z') } axes3d(edges = "bbox",labels=TRUE,tick=FALSE,box=TRUE,nticks=0, expand=1.0,xlen=0,xunit=0,ylen=0,yunit=0,zlen=0,zunit=0) if(draw.box) { x <- box\$xrange[2] y <- box\$yrange[2] z <- box\$zrange[2] c3d.origin <- rgl::translate3d(rgl::scale3d(rgl::cube3d(col="gray", alpha=0.1), x/2,y/2,z/2),x/2,y/2,z/2) rgl::shade3d(c3d.origin) } } col <- c("#0000FF","#00FF00","#FF0000","#FF00FF","#FFFF00","#00FFFF") ################################################################# ## `beta` distribution for radii ################################################################# lam <- 50 ## parameter beta distribution (radii) theta <- list("size"=list("shape1"=1,"shape2"=10)) # simulation bounding box box <- list("xrange"=c(-1,4),"yrange"=c(-1.5,3.5),"zrange"=c(0,2)) ## simulate and return full spheres system ## intersect with XZ plane and return full list of intersection profiles ## section plane xy S <- simPoissonSystem(theta,lam,size="rbeta",box=box,type="spheres", intersect="full",n=c(0,0,1),dz=0,pl=1) # check resulting distribution length(S\$S) summary(sapply(S\$S,"[[","r")) theta\$size[[1]]/(theta\$size[[1]]+theta\$size[[2]]) # mean ## interior spheres: ## the ones which intersect one of the lateral planes (without top/bottom planes) ## showing spheres with color intersect notIn <- sapply(S\$S,function(x) !attr(x,"interior")) spheres(S\$S[notIn],box,FALSE,TRUE,color=col) # not intersecting In <- sapply(S\$S,function(x) attr(x,"interior")) spheres(S\$S[In],box,color="gray") ## full sphere system # open3d() # spheres(S\$S,box,FALSE,TRUE,color=col) # planes3d(0,0,1,0,col="darkgray",alpha=1) ## draw intersections sp <- S\$sp id <- sapply(sp,"[[","id") open3d() spheres(S\$S[id],box,FALSE,TRUE,color=col) planes3d(0,0,1,0,col="darkgray",alpha=1) XYr <- t(sapply(sp,function(s) cbind(s\$center[1],s\$center[2],s\$r))) # centers x <- XYr[,1] y <- XYr[,2] r <- XYr[,3] win <- attr(sp,"win") xlim <- win[[1]] ylim <- win[[2]] dev.new() plot(x,y,type="n",xaxs="i", yaxs="i", xlab="x",ylab="y",xlim=xlim,ylim=ylim) for(i in 1:nrow(XYr)) draw.circle(x[i],y[i],r[i],nv=100,border=NULL,col="black") ## digitize inersections ## `win` can also be omitted if not different ## from original itersection window (derived from the box) W <- digitizeProfiles(sp, delta=0.01, win=win) dim(W) dev.new() image(1:nrow(W),1:ncol(W),W,col=gray(1:0)) ################################################################# ## Exact simulation of spheres with log normal radii ################################################################# lam <- 100 ## parameter rlnorm distribution (radii) theta <- list("size"=list("meanlog"=-2.5,"sdlog"=0.5)) # simulation bounding box box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5)) ## simulate only 3D system S <- simPoissonSystem(theta,lam,size="rlnorm",box=box,type="spheres", intersect="original", perfect=TRUE, pl=101) ## show open3d() spheres(S[1:2000],box,TRUE,TRUE,color=col) ## check! mean(log(sapply(S,"[[","r"))) sd(log(sapply(S,"[[","r"))) ################################################################# ## Planar section ################################################################# # planar section of exact simulated `rlnorm` sphere system: # returns diameters for those section profiles # which have their centers inside the intersection window sp <- planarSection(S,d=2.5,intern=TRUE,pl=1) hist(sp) summary(sp) mean(log(sp/2)) sd(log(sp/2)) ################################################################# ## General intersection, all objects (inter=FALSE) ################################################################# SP <- intersectSystem(S, 2.5, n=c(0,0,1), intern=FALSE, pl=1) ## show in 3D id <- sapply(SP,"[[","id") open3d() spheres(S[id],box,TRUE,color=col) planes3d(0,0,-1,2.5,col="black",alpha=1) ## 2D sections XYr <- t(sapply(SP,function(s) cbind(s\$center[1],s\$center[2],s\$r))) # centers x <- XYr[,1] y <- XYr[,2] r <- XYr[,3] xlim <- box\$xrange ylim <- box\$yrange dev.new() plot(x,y,type="n",xaxs="i", yaxs="i", xlab="x",ylab="y",xlim=xlim,ylim=ylim) for(i in 1:nrow(XYr)) draw.circle(x[i],y[i],r[i],nv=100,border=NULL,col="black") # digitize W <- digitizeProfiles(SP, delta=0.01) dev.new() image(1:nrow(W),1:ncol(W),W,col=gray(1:0)) ################################################################# ## Unfolding ################################################################# ret <- unfold(sp,nclass=25) ## Point process intensity cat("Intensities: ", sum(ret\$N_V)/25, "vs.",lam,"\n") ## original diameters r3d <- unlist(lapply(S,function(x) 2*x\$r)) rest3d <- unlist(lapply(2:(length(ret\$breaks)), function(i) rep(ret\$breaks[i],sum(ret\$N_V[i-1])))) op <- par(mfrow = c(1, 2)) hist(r3d[r3d<=max(ret\$breaks)], breaks=ret\$breaks, main="Radius 3D", freq=FALSE, col="gray",xlab="r") hist(rest3d, breaks=ret\$breaks,main="Radius estimated", freq=FALSE, col="gray", xlab="r") par(op) ################################################################# ## Update intersection: find objects which intersect bounding box ################################################################# idx <- updateIntersections(S) sum(!idx) # objects intersecting id <- which( idx != 1) open3d() spheres(S[id],box,TRUE,TRUE,color=col) ################################################################# ## user-defined simulation function ################################################################# lam <- 50 theta <- list("p1"=-2.5,"p2"=0.5) myfun <- function(p1,p2) { c("r"=rlnorm(1,p1,p2)) } box <- list("xrange"=c(0,5),"yrange"=c(0,5),"zrange"=c(0,5)) S <- simPoissonSystem(theta,lam,rjoint=myfun,box=box,type="spheres",pl=1) r <- sapply(S,"[[","r") mean(log(r)) sd(log(r)) open3d() spheres(S,box,TRUE,TRUE,color=col) ## End(Not run) ```

### Example output

```Warning messages:
1: In rgl.init(initValue, onlyNULL) : RGL: unable to open X11 display
2: 'rgl.init' failed, running with 'rgl.useNULL = TRUE'.

Attaching package: ‘plotrix’

The following object is masked from ‘package:rgl’:

mtext3d

Size parameters: mx=1.000000, sdx=10.000000

Spheres simulation with `rbeta` (perfect=0):
Mean number: 50.000000 (Box volume 50.000000)
Set label 'N'.

Done. Simulated 2534 spheres.
Converting spheres ...
Converting 108 discs.
Getting plane indices: [0 1 ]
[1] 2534
Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
0.0000087 0.0277244 0.0677158 0.0914636 0.1312227 0.5309968
[1] 0.09090909
null
2
[1] 500 500
dev.new(): using pdf(file="Rplots1.pdf")
Size parameters: mx=-2.500000, sdx=0.500000
Cumulative sum of probabilities: 0.896168, 0.996196, 0.999949, 1.000000 ( mu=139.482809 )

Spheres simulation with `rlnorm` (perfect=1):
Mean number: 100.000000 (exact simulation: 139.482809)
Set label 'N'.

perfect
Done. Simulated 13978 spheres.
Converting spheres ...
null
3
[1] -2.473883
[1] 0.5009234
Converting 486 discs.
Getting plane indices: [0 1 ]
Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
0.005738 0.102507 0.160404 0.179879 0.234117 0.658794
[1] -2.580125
[1] 0.6209039
Converting 542 discs.
Getting plane indices: [0 1 ]
null
4
dev.new(): using pdf(file="Rplots2.pdf")
dev.new(): using pdf(file="Rplots3.pdf")
Intensities:  105.895 vs. 100
Warning: stack imbalance in '<-', 2 then 1
[1] 2740
null
5
Spheres simulation by `myfun`
Box volume: 125.000000, lam: 50.000000, number of spheres: 6100
Done. Simulated 6100 spheres.
Converting spheres ...
[1] -2.497105
[1] 0.4907834
null
6
```

unfoldr documentation built on Sept. 25, 2021, 1:07 a.m.