Add single photon imaging demo

ch08.r

@@ -47,6 +47,27 @@ lines(theta, L_df$"12", lwd=6)
 title(expression(paste("L(", theta, " | S = 12)")))
 grid()
 
+## ML estimation for single-photon imaging
 
+# R code
+library(imager)
+lambda = as.data.frame(load.image("cameraman.tif"))
+lambda = xtabs(value ~ x+y, data=lambda)
+T = 100
+x = c()
+for (i in 1:T) x = append(x, rpois(length(lambda), lambda))
+x = array(x, c(256, 256, 100))
+y = (x>=1)
+mu = apply(y, c(1,2), mean)
+lambdahat = -log(1-mu)
+fig1 = x[,,1]
+
+# Flip matrix since image reads the matrix bottom up
+flip_matrix = function(m) m[,nrow(m):1]
+
+image(flip_matrix(fig1), col=gray.colors(255))
+
+image(flip_matrix(lambdahat), col=gray.colors(255))
 
 #############
+