Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Calculates and visualises proportion of fibres that can be extracted/sampled, from a skeletal muscle section leaving all immediate neighbours of sampled fibres in place.
library(deldir)
library(png)
dat = read.delim("results.csv",sep=",",stringsAsFactors=FALSE)
im = readPNG("AVE_mitocyto.png")
h = dim(im)[1]
w = dim(im)[2]
td = data.frame(
xCoord = dat$Value[dat$Channel=="xCoord"],
yCoord = h-dat$Value[dat$Channel=="yCoord"],
Area = dat$Value[dat$Channel=="Area"],
cell_id = 1:length(dat$Value[dat$Channel=="xCoord"])
)
cdat = data.frame(x=td$xCoord,y=td$yCoord)
meanrad = mean(sqrt(td$Area/(2*pi)))
del = deldir(cdat)
del$delsgs$dist = sqrt((del$delsgs$x2-del$delsgs$x1)^2+(del$delsgs$y2-del$delsgs$y1)^2)
adj = matrix(0.0,nrow=nrow(cdat),ncol=nrow(cdat))
for(i in 1:length(del$delsgs$dist)){
distance = del$delsgs$dist[i]
if(distance>5*meanrad) distance = 0.0
adj[del$delsgs$ind1[i],del$delsgs$ind2[i]] = distance
adj[del$delsgs$ind2[i],del$delsgs$ind1[i]] = distance
}
adj = round(adj)
colnames(adj) = td$cell_id
rownames(adj) = td$cell_id
candidates = sample(td$cell_id)
unavailable = c()
sampled = c()
pdf("Samples.pdf",width=16,height=8)
#png("frames/frame%05d.png",width=1600, height=800,pointsize=18)
while(length(candidates)>0){
cellid = candidates[1]
neighb_1 = td$cell_id[adj[,cellid]>0.0]
neighb_2 = c()
for(n in neighb_1){
neighb_2 = unique(c(neighb_2,td$cell_id[adj[,n]>0.0]))
}
neighb_2 = neighb_2[neighb_2!=cellid]
sampled = c(sampled, cellid)
#unavailable = unique(c(unavailable,neighb_1))
unavailable = unique(c(unavailable,neighb_1,neighb_2))
candidates = candidates[!((candidates%in%unavailable)|(candidates%in%sampled))]
}
sdat=cdat[td$cell_id%in%sampled,]
udat=cdat[td$cell_id%in%unavailable,]
op = par(mfrow=c(1,2))
prop = length(sampled)/length(cdat$x)
mlab = paste("Proportion sampled:",formatC(100*prop,3),"(%)")
rang = range(c(cdat$x,cdat$y))
plot(del)
points(sdat$x,sdat$y,pch=16,col="red",cex=0.75)
points(udat$x,udat$y,pch=16,col="blue",cex=0.75)
plot(cdat$x,cdat$y,type="n",xlab="x (px)",ylab="y (px)",cex.axis=1.45,cex.lab=1.45,xlim=c(0,w),ylim=c(0,h),main=mlab)
lim <- par()
rasterImage(im, 0,0,w,h)
points(cdat$x,cdat$y,pch=16,col="yellow",cex=0.75)
points(sdat$x,sdat$y,pch=16,col="red",cex=0.55)
points(udat$x,udat$y,pch=16,col="blue",cex=0.55)
par(op)
#}
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.