### R code from vignette source 'PLLPanalysis.Rnw'

###################################################
### code chunk number 1: installPackages (eval = FALSE)
###################################################
## if (!requireNamespace("BiocManager", quietly=TRUE))
##     install.packages("BiocManager")
## BiocManager::install(c("EBImage", "parallel"))


###################################################
### code chunk number 2: loadPackages
###################################################
library("DonaPLLP2013")


###################################################
### code chunk number 3: loadImages
###################################################
xGFP <- readImage(system.file("extdata/cxcr4b_02C1.tif", package="DonaPLLP2013")) 
xRFP <- readImage(system.file("extdata/cxcr4b_02C2.tif", package="DonaPLLP2013")) 


###################################################
### code chunk number 4: displayImages
###################################################
dim(xGFP)
xGFP.proj <- apply(xGFP, c(1, 2), max)
writeImage(normalize(xGFP.proj), "PLLPanalysis-displayImages.jpeg")


###################################################
### code chunk number 5: setUnits
###################################################
Lpx <- c(x=0.1318, y=0.1318, z=1)
Lbox <- Lpx*dim(xGFP)  
Lbox
Leff <- round(5/Lpx["x"])           
Leff


###################################################
### code chunk number 6: backgroundSubtraction
###################################################
xbckGFP <- readImage(system.file("extdata/bgC1.tif", package="DonaPLLP2013"))
xbckRFP <- readImage(system.file("extdata/bgC2.tif", package="DonaPLLP2013"))
xGFP <- xGFP-mean(xbckGFP, na.rm=TRUE)
xRFP <- xRFP-mean(xbckRFP, na.rm=TRUE)


###################################################
### code chunk number 7: gaussianSmoothing
###################################################
Lpsf <- round(0.5/Lpx["x"])     
Lpsfodd <- ifelse(Lpsf%%2 == 0, Lpsf+1, Lpsf)
z <- makeBrush(size=Lpsfodd, shape="gaussian", sigma=Lpsfodd/2)
x2GFP <- filter2(xGFP, filter=z)
x2RFP <- filter2(xRFP, filter=z)
x2GFP.proj <- apply(x2GFP, c(1, 2), max)
writeImage(normalize(x2GFP.proj), "PLLPanalysis-gaussianSmoothing.jpeg")


###################################################
### code chunk number 8: adaptiveThresholding
###################################################
mk <- function(size) makeBrush(size, shape="disc")
mask <- thresh(x2GFP, w=Leff, h=Leff, offset=0.01)
mask <- erode(closing(mask, mk(Lpsfodd)), mk(Lpsfodd-2))
maskSlice <- mask[, , 20]
writeImage(maskSlice, "PLLPanalysis-maskSlice.jpeg")
maskDensity.proj <- apply(mask, c(1, 2), mean)
writeImage(normalize(maskDensity.proj), "PLLPanalysis-maskDensity.jpeg")


###################################################
### code chunk number 9: obtainOrganMask
###################################################
organ.mask <- apply(fillHull(mask), c(1, 2), max) 
organ.mask <- opening(closing(organ.mask, mk(Leff/2)), mk(Leff/2))
organ.labels <- bwlabel(organ.mask)
I <- which(table(organ.labels)[-1] < max(table(organ.labels)[-1])) 
organ.mask <- rmObjects(organ.labels, I)
organ.mask[organ.mask > 0] <- 1
organ.mask <- fillHull(organ.mask)
writeImage(organ.mask, "PLLPanalysis-organMask.jpeg")


###################################################
### code chunk number 10: principalComponentAnalysis
###################################################
organ.coord <- as.data.frame(which(organ.mask == 1, arr.ind=TRUE),
                             stringsAsFactors=FALSE) 
colnames(organ.coord) <- c("x", "y")
organ.coord$x <- (organ.coord$x-1)*Lpx["x"]
organ.coord$y <- -(organ.coord$y-1)*Lpx["y"]
organ.pca <- prcomp(organ.coord)
pc1 <- organ.pca$rotation[, 1]
pc1.theta <- atan(pc1["y"]/pc1["x"])*(180/pi)
pc1.theta


###################################################
### code chunk number 11: rotateImages
###################################################
xGFP <- rotate(xGFP, angle=pc1.theta)
xRFP <- rotate(xRFP, angle=pc1.theta)
x2GFP <- rotate(x2GFP, angle=pc1.theta)
x2RFP <- rotate(x2RFP, angle=pc1.theta)
mask <- rotate(mask, angle=pc1.theta)
organ.mask <- rotate(organ.mask, angle=pc1.theta)


###################################################
### code chunk number 12: excludeExteriorPixels
###################################################
organ.mask.rep <- replicate(dim(xGFP)[3], organ.mask)
I1 <- which(organ.mask.rep == 1)
I2 <- which(mask == 1)
I <- intersect(I1, I2)
mask <- array(0, dim=dim(organ.mask.rep))
mask[I] <- 1


###################################################
### code chunk number 13: cropImages
###################################################
I <- which(organ.mask.rep == 1, arr.ind=TRUE)
intRange <- function(x) {rg=range(x); seq(rg[1], rg[2], by=1)}
x.range <- intRange(I[, 1])
y.range <- intRange(I[, 2])
xGFP <- xGFP[x.range, y.range, ]
xRFP <- xRFP[x.range, y.range, ]
x2GFP <- x2GFP[x.range, y.range, ]
x2RFP <- x2RFP[x.range, y.range, ]
mask <- mask[x.range, y.range, ]
organ.mask <- organ.mask[x.range, y.range]
writeImage(organ.mask, "PLLPanalysis-cropImages.jpeg")


###################################################
### code chunk number 14: cropCoordinates
###################################################
Lbox <- Lpx*dim(xGFP)
Lbox


###################################################
### code chunk number 15: runningMedianFunctions
###################################################
getCoordinates <-
function(s, xrange, yrange, zrange) {
    Ix=which(s$x >= min(xrange) & s$x <= max(xrange))
    Iy=which(s$y >= min(yrange) & s$y <= max(yrange))
    Iz=which(s$z >= min(zrange) & s$z <= max(zrange))
    return(list(x=Ix, y=Iy, z=Iz))
}

getCubeIntensity <-
function(x0, x, y, z, spatial, Lcube) {
    crd <- getCoordinates(spatial,
                          xrange=c(x-Lcube/2, x+Lcube/2),
                          yrange=c(y-Lcube/2, y+Lcube/2),
                          zrange=c(z-Lcube/2, z+Lcube/2))
    cube.median <- median(imageData(x0)[crd$x, crd$y, crd$z], na.rm=TRUE)
    return(list(median=cube.median))
}

runningMedian <-
function(x, grid.x, grid.y, grid.z, Lx, Ly, Lz, nCores=2, Lcube) {
    grid <- as.list(data.frame(t(expand.grid(grid.x, grid.y, grid.z))))
    spatial <- list(x=(1:dim(x)[1]-1)*Lx, 
                    y=(1:dim(x)[2]-1)*Ly,
                    z=(1:dim(x)[3]-1)*Lz)
    chooseCores <- function(numCoresWanted) {
        if(.Platform$OS.type == "windows") return(1)
        return(numCoresWanted)
    }
    dataIntensities <-
        mclapply(grid, function(s) getCubeIntensity(x0=x, x=s[1], y=s[2], z=s[3],
                                                    spatial=spatial, Lcube=Lcube), 
                 mc.cores=chooseCores(nCores), 
                 mc.preschedule=FALSE)
    dataIntensities=unlist(dataIntensities)
    result.median=array(dataIntensities,
                        dim=c(length(grid.x), length(grid.y), length(grid.z)))
    return(result.median)
}


###################################################
### code chunk number 16: runningMedianGrid
###################################################
Ljump <- 5
Lcube <- 10
grid.x <- seq(from=0, to=Lbox["x"], by=Ljump)
grid.y <- seq(from=0, to=Lbox["y"], by=Ljump)
grid.z <- seq(from=0, to=Lbox["z"], by=Ljump)


###################################################
### code chunk number 17: runningMedian
###################################################
I <- which(mask == 0)
xGFP.maskOnly <- xGFP
xRFP.maskOnly <- xRFP
xGFP.maskOnly[I] <- NA
xRFP.maskOnly[I] <- NA
resultGFP <- runningMedian(x=xGFP.maskOnly, 
                           grid.x=grid.x, grid.y=grid.y, grid.z=grid.z, 
                           Lx=Lpx["x"], Ly=Lpx["y"], Lz=Lpx["z"], 
                           nCores=4, Lcube=Lcube)
resultRFP <- runningMedian(x=xRFP.maskOnly, 
                           grid.x=grid.x, grid.y=grid.y, grid.z=grid.z, 
                           Lx=Lpx["x"], Ly=Lpx["y"], Lz=Lpx["z"], 
                           nCores=4, Lcube=Lcube)


###################################################
### code chunk number 18: plotFluorescenceChannels
###################################################
GFP.profile=apply(resultGFP, 1, median, na.rm=TRUE)
RFP.profile=apply(resultRFP, 1, median, na.rm=TRUE)
plot(-grid.x, rev(GFP.profile), 
     xlab="Distance from leading edge (microns)", 
     ylab="Fluorescence intensity (a.u.)", 
     type="b", 
     axes="F", 
     pch=1, 
     ylim=c(0, range(GFP.profile, RFP.profile, na.rm=TRUE)[2]), 
     col="darkgreen")
points(-grid.x, rev(RFP.profile), type="b", pch=2, col="red")
axis.at.x <- seq(-200 ,0 , by=25)
axis(1, at=axis.at.x, labels=-axis.at.x)
axis(2)
legend("topright", legend=c("GFP", "RFP"), pch=1:2, col=c("darkgreen", "red"))


###################################################
### code chunk number 19: plotRatio
###################################################
resultRatio <- resultRFP/resultGFP
ratio.profile <- apply(resultRatio, 1, median, na.rm=TRUE)
plot(-grid.x, rev(ratio.profile), 
     xlab="Distance from leading edge (microns)", 
     ylab="RFP/GFP fluorescence intensity ratio (-)", 
     type="b", 
     axes="F")
axis.at.x <- seq(-200, 0, by=25)
axis(1, at=axis.at.x, labels=-axis.at.x)
axis(2)


###################################################
### code chunk number 20: solutionBackgroundSubtraction
###################################################
solGFP <- readImage(system.file("extdata/1_100_mCherryC1.tif",
                    package="DonaPLLP2013"))
solGFPbck <- readImage(system.file("extdata/PBSC1.tif", package="DonaPLLP2013"))
solRFP <- readImage(system.file("extdata/1_100_mCherryC2.tif",
                    package="DonaPLLP2013"))
solRFPbck <- readImage(system.file("extdata/PBSC2.tif", package="DonaPLLP2013"))
solRFP <- solRFP-mean(solRFPbck)
solGFP <- solGFP-mean(solGFPbck)


###################################################
### code chunk number 21: correctSampleRatioByControlRatio
###################################################
solRatio <- mean(solRFP)/mean(solGFP)
resultRatioCorrected <- resultRatio/solRatio


###################################################
### code chunk number 22: wholeSampleRatio
###################################################
print(median(resultRatioCorrected, na.rm=TRUE))


###################################################
### code chunk number 23: varyMembraneThickness
###################################################
membraneRatio <- array(dim=5)
I <- which(mask == 1)
membraneRatio[3] <- median(xRFP[I])/median(xGFP[I])

calcRatioOnModifiedMembrane <- function(xGFP, xRFP, mask, morph=erode, steps=1) {
 for (i in seq_len(steps)) mask <- morph(mask, mk(2))
 I <- which(mask == 1)
 return(median(xRFP[I])/median(xGFP[I]))
}

membraneRatio[1] <- calcRatioOnModifiedMembrane(xGFP, xRFP, mask, erode, 2)
membraneRatio[2] <- calcRatioOnModifiedMembrane(xGFP, xRFP, mask, erode, 1)
membraneRatio[4] <- calcRatioOnModifiedMembrane(xGFP, xRFP, mask, dilate, 1)
membraneRatio[5] <- calcRatioOnModifiedMembrane(xGFP, xRFP, mask, dilate, 2)
plot(1:5, membraneRatio, xlab="Mask ID", 
     ylab="RFP/GFP fluorescence intensity ratio (-)", 
     ylim=c(0, max(membraneRatio)))