dykifier.server <- function(input,output,session){
library("inline")
threshold <- reactiveValues()
threshold$pvalue <- 0.02
threshold$lvalue <- -0.01
X <- c()
for(i in 1:length(segments) ){
X = rbind(X,segments[[i]])
}
type <- rep(3, nrow(X) )
type[ rtype == "quartz diorite" ] <- 1
type[ rtype == "basalt" ] <- 2
tcolors <- brewer.pal("Set1", n=3)
ranges <- reactiveValues(x = NULL, y = NULL)
intersection.probs <- reactive({
input$update
a.ext <- isolate( input$a.ext )
an <- 1 / ( sqrt(2*pi)*a.ext)
av <- 2*a.ext^2
withProgress(message = 'Updating probabilities', value = 0, {
ns <- length(segments)
k <- ncol(angles1)
na <- nrow(angles1)
p <- matrix( NA, ncol=5, nrow = k*na )
index <- 1
for(i in 1:na ){
incProgress( 0.5 / na )
si = i %% length(segments)
if(si==0){
si = length(segments)
}
d.ext <- sdevs[si, 1] * isolate( input$d.ext )
dv <- 2*d.ext^2
dn <- 1 / ( sqrt(2*pi)*d.ext )
pd <- dn * exp( - dists[i,]^2 / dv )
pa1 <- av * exp( - (angles1[i, ] - 1)^2 / av )
pa2 <- av * exp( - (angles2[i, ] - 1)^2 / av )
#p[i, ] <- pd*pa1*pa2
p[index:(index+k-1), ] <- cbind(pd*pa1*pa2, rep(i, k),
KNN$nn.index[i, ],
rep(i, k) %% ns,
KNN$nn.index[i, ] %% ns )
index <- index+k
}
p[ p[, 4] == 0, 4 ] = ns
p[ p[, 5] == 0, 5 ] = ns
})
porder <- order(p[,1], decreasing=TRUE)
segs <- matrix(NA, ncol=7, nrow=na )
cpp_segs_src <- '
Rcpp::IntegerVector porder(po);
Rcpp::NumericMatrix p(probs);
Rcpp::NumericMatrix segs(hsegs);
Rcpp::NumericMatrix coords(x);
int np = porder.size();
int na = segs.nrow();
bool *used = new bool[na];
for(int i=0; i<na; i++){
used[i] = false;
}
int index = 0;
for(int i=0; i<np; i++){
int i1 = p(porder[i]-1, 1)-1;
if( used[i1] ){
continue;
}
int i2 = p(porder[i]-1, 2)-1;
if( used[i2] ){
continue;
}
segs(index, 0) = p( porder[i]-1, 0);
segs(index, 1) = coords(i1, 0);
segs(index, 2) = coords(i1, 1);
segs(index, 3) = coords(i2, 0);
segs(index, 4) = coords(i2, 1);
segs(index, 5) = p( porder[i]-1, 3 );
segs(index, 6) = p( porder[i]-1, 4 );
used[i1] = true;
used[i2] = true;
index++;
}
delete[] used;
return segs;
'
cpp_segs <- cxxfunction(
signature(po="integer", probs="matrix",
hsegs="matrix", x="matrix" ),
cpp_segs_src, plugin="Rcpp")
segs <- cpp_segs(porder, p, segs, coords)
# used <- rep(FALSE, na )
# index <- 1
# for(i in 1:length(porder) ){
# incProgress( 0.5 / (k*na) )
# i1 <- p[porder[i], 2]
# if( used[i1] ){
# next
# }
# i2 <- p[porder[i], 3]
# if( used[i2] ){
# next
# }
# segs[index, ] <- c( p[ porder[i], 1], coords[i1, ], coords[i2, ] )
# used[i1] = TRUE
# used[i2] = TRUE
# index <- index+1
# }
# })
segs <- segs[complete.cases(segs), ]
segs[,1] <- segs[,1] / max(segs[,1])
colnames(segs) <- c("Probability", "x1", "y1", "x2", "y2", "sid1", "sid2")
segs
})
intersection.selected <- reactive({
segs <- intersection.probs()
segs <- segs[ which( segs[,1] >= threshold$pvalue ) , ]
segs
})
intersection.selected.connected <- reactive({
segs <- intersection.selected()
id.segs <- rep( NA, length(segments) )
id.ph <- rep( NA, nrow(segs) )
did <- 0
for( i in 1:nrow(segs) ){
cid <- NA
if( !is.na( id.segs[ segs[i,6] ]) ){
cid = id.segs[ segs[i,6] ]
}
if( !is.na( id.segs[ segs[i,7] ]) ){
if( !is.na(cid) ){
#replace with single id
eid <- id.segs[ segs[i,7] ]
ind <- which(id.segs == eid)
id.segs[ind] = cid
ind <- which(id.ph == eid)
if(length(ind) > 0 ){
id.ph[ind] = cid
}
}
else{
cid = id.segs[ segs[i,7] ]
}
}
if( is.na(cid) ){
did <- did + 1
cid <- did
}
id.ph[i] = cid
id.segs[ segs[i,6] ] = cid
id.segs[ segs[i,7] ] = cid
}
ndykes = length( unique(id.segs) )
for( i in 1:length(id.segs) ){
if( is.na( id.segs[i] ) ){
did = did+1
id.segs[i] = did
}
}
dlengths = rep(0, did+1)
for(i in 1:did){
ind <- which(id.ph == i)
if(length(ind) > 0 ){
xd <- (segs[ind,2] - segs[ind,4])^2
yd <- (segs[ind,3] - segs[ind,5])^2
dlengths[i] <- sqrt( sum(xd+yd) )
}
ind <- which(id.segs == i)
if(length(ind) > 0 ){
dlengths[i] <- dlengths[i] + sum(slength[ind] )
}
}
list( ndykes = ndykes, nsegs = length( unique(id.segs) ) - ndykes, id.segs = id.segs,
id.ph = id.ph, dlengths = dlengths )
})
## Dykes plot ##
output$dykes <- renderPlot({
plot(X, pch=".", asp=1, xlab="x", ylab="y")
csegs <- intersection.selected.connected()
for( i in 1:length(segments) ){
if( csegs$dlengths[ csegs$id.segs[i] ] > threshold$lvalue ){
lines( segments[[i]], col=tcolors[type[i]], lwd=3 )
}
}
hsegs <- intersection.selected()
for( i in 1:nrow(hsegs) ){
if( csegs$dlengths[ csegs$id.ph[i] ] > threshold$lvalue){
segments( hsegs[i,2], hsegs[i, 3], hsegs[i,4], hsegs[i,5], col="gray", lwd=2)
}
}
})
observe({
brush <- input$dykes_brush
if (!is.null(brush)) {
ranges$x <<- c(brush$xmin, brush$xmax)
ranges$y <<- c(brush$ymin, brush$ymax)
}
})
output$dykes.zoom <- renderPlot({
plot(X, pch=".", asp=1, xlim = ranges$x, ylim = ranges$y, xlab="x", ylab="y")
csegs <- intersection.selected.connected()
for( i in 1:length(segments) ){
if( csegs$dlengths[csegs$id.segs[i]] > threshold$lvalue ){
lines( segments[[i]], col=tcolors[type[i]], lwd=3 )
}
}
hsegs <- intersection.selected()
for( i in 1:nrow(hsegs) ){
if(csegs$dlengths[csegs$id.ph[i] ] > threshold$lvalue ){
segments( hsegs[i,2], hsegs[i, 3], hsegs[i,4], hsegs[i,5], col="gray", lwd=2)
}
}
})
output$histogram <- renderPlot({
par(mar=c(4,4,1,1) )
hsegs <- intersection.probs()
hist(hsegs[,1], 200, main="", xlab="Probability")
abline(v=threshold$pvalue)
})
ob.hclick <- observeEvent(input$h_click, {
if( !is.null( input$h_click) ){
ev <- input$h_click
threshold$pvalue <- input$h_click$x
}
})
output$lhistogram <- renderPlot({
par(mar=c(4,4,1,1) )
hsegs <- intersection.selected()
xd <- hsegs[,2] - hsegs[,4]
yd <- hsegs[,3] - hsegs[,5]
ls <- sqrt(xd^2 + yd^2)
hist(ls, input$nlbins, main="", xlab="Phantom Segment Lengths")
})
output$dhistogram <- renderPlot({
par(mar=c(4,4,1,1) )
dsegs <- intersection.selected.connected()
hist(dsegs$dlengths, input$ndbins, main="", xlab="Dike Lengths")
abline(v=threshold$lvalue)
})
ob.dhclick <- observeEvent(input$dh_click, {
if( !is.null( input$dh_click) ){
ev <- input$dh_click
threshold$lvalue <- input$dh_click$x
}
})
output$downloadData <- downloadHandler( "phantom-dikes.csv",
content = function(file) {
write.csv(intersection.selected(), file)
}
)
output$ndykes.info <- renderText({
cdykes <- intersection.selected.connected()
paste("Number of connected dykes: ", cdykes$ndykes )
})
output$nsegs.info <- renderText({
cdykes <- intersection.selected.connected()
paste("Number of disconnected segments left: ", cdykes$nsegs )
})
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.