## R script to confirm that the data set has been created okay. create.txt.files <- FALSE if (create.txt.files) { attach("~/mosaics/data/w81s1.Rda") write.table(w81s1.on, file="w81s1on.txt", quote=F, col.names=FALSE, row.names=FALSE) write.table(w81s1.of, file="w81s1of.txt", quote=F, col.names=FALSE, row.names=FALSE) write(w81s1.w, file="w81s1w.txt", ncolumns=4) detach() } ## Read in data files. w81s1.on <- read.table("w81s1on.txt") w81s1.of <- read.table("w81s1of.txt") w81s1.w <- read.table("w81s1w.txt") w81s1.soma.r <- 9 postscript(file="w81s1_field.ps", horiz=F, width=7, height=7) par(mar=c(2,2,1,1)) plot(NA, ylim=c(1000,0), xlim=c(0,750), asp=1) symbols(w81s1.on, circles=rep(w81s1.soma.r, dim(w81s1.on)[1]), add=TRUE, inches=F, pch=19, fg='green' ) symbols(w81s1.of, circles=rep(w81s1.soma.r, dim(w81s1.of)[1]), add=T, inches=F, asp=1,pch=19, fg='orangered' ) rect(w81s1.w[1], w81s1.w[3], w81s1.w[2], w81s1.w[4], lty=2) dev.off() simple.xnnd <- function(on, of=NULL) { ## Do a simple XNND analysis. ## PTS is a (N,3) matrix, each row being (X, Y, C) where C is the class ## of the cell (type 1 or 2). if (is.null(of)) { ## wanting an autocorrelogram of <- on auto <- TRUE } else { auto <- FALSE } n1 <- dim(on)[1]; n2 <- dim(of)[1] res <- matrix(0, nrow=n2, ncol=3) for (i in 1:n2) { min.d <- 9999999; min.j <- -1 for (j in 1:n1) { if (auto && (i==j) ) { next; } d2 <- ((of[i,1] - on[j,1])^2) + ((of[i,2] - on[j,2])^2) if (d2 < min.d) { min.d <- d2; min.j <- j } } res[i,] <- c(sqrt(min.d), min.j, i) } colnames(res) <- c("NND", "type 1", "type 2") overall.min <- which.min(res[,1]) cat("overall min\n") cat(paste(res[overall.min,]),"\n") res } simple.xnnd(w81s1.on, w81s1.of) ##source('~/langs/R/simple_nnd.r') simple.nnd <- function(x, y) { ## Do a simple NND analysis, standalone, no C or anything fancy. npts <- length(x) res <- matrix(0, nrow=npts, ncol=2) for (i in 1:npts) { min.d <- 9999999; min.j <- -1 for (j in 1:npts) { if (!(i==j)) { d2 <- ((x[i] - x[j])^2) + ((y[i] - y[j])^2) if (d2 < min.d) { min.d <- d2; min.j <- j } } } res[i,] <- c(sqrt(min.d), min.j) } res } simple.hist <- function(cells, label) { nnd <- simple.nnd(cells[,1], cells[,2]) mu <- mean(nnd[,1]); s <- sd(nnd[,1]); ri <- mu/s title <- sprintf("%s N %d mean %.2f sd %.1f RI %.1f", label, dim(cells)[1], mu, s, ri) hist(nnd[,1], main=title, breaks=seq(from=0, to=150, by=5), xlab="distance to nearest neighbour (um)") } postscript(file="w81s1_nnd.ps", width=5, height=8) par(mfrow=c(3,1)) simple.hist(rbind(w81s1.of, w81s1.on), "All") simple.hist(w81s1.on, "On") simple.hist(w81s1.of, "Off") dev.off()