#' An example of an input shoeprint
#' @format a data frame
"input_example"
#' An example of a reference shoeprint
#' @format a data frame
"reference_example"
#' @title Implementation of a smart sampling algorithm
#'
#' @description Define a smart sampling function which will ensure uniform samples so that we can get better results even with a lower ncross value. This process will divide the input circle area into hexbins and assign probabilities to the bins which will be correlated with the number of points in the bin. This will result in sampling the points from high density areas while still ensuring a more uniform distribution.
#' @name smart_sample
#' @param pts_xy The set of points (for example, the input circle of a shoeprint)
#' @param sample_bins The number of hexbins to sample
#' @param sample_size The number of points within each bin to sample
#' @param xbins Defines the grid of hexbins to divide the overall image (by default, 20x20)
#' @param seed The random seed for reproducing results
#' @return A data frame consisting of the sampled points from the original input circle
#' @export
#' @import hexbin
#' @examples
#' \dontrun{
#' data(input_example)
#'
#' smart_sample(input_example, sample_bins = NULL, sample_size = NULL) # All points from all bins
#' smart_sample(input_example, sample_bins = 30 & sample_size = NULL) # All points from 30 bins
#' smart_sample(input_example, sample_bins = NULL, sample_size = 30) # 30 points from all bins
#' smart_sample(input_example, sample_bins = 30, sample_size = 1) # one point from 30 selected bins
#' }
##############################################################################################################
smart_sample <- function(pts_xy, sample_bins = NULL, sample_size = NULL, xbins = 20, seed = 1) {
## Generate a hexbin object
pts_hxb <- hexbin(pts_xy,xbins=xbins,IDs=TRUE)
## Assign cell ID to data points
pts_xy$cell_id <- pts_hxb@cID
## Get number of points in every bin and sample bins to choose points from
hxb_stats <- data.frame(cell_id=pts_hxb@cell,count=pts_hxb@count)
set.seed(seed)
candidate_hxbs <- sample(hxb_stats$cell_id,ifelse(is.null(sample_bins),length(pts_hxb@cell),min(nrow(hxb_stats),sample_bins)),prob=hxb_stats$count/sum(hxb_stats$count))
## Candidate points and final sample of points
pts_xy <- pts_xy[pts_xy$cell_id %in% candidate_hxbs,]
if(is.null(sample_size)) return(pts_xy[,-3])
set.seed(seed)
return(pts_xy[unlist(tapply(1:nrow(pts_xy),pts_xy$cell_id,function(x) sample(x,min(length(x),sample_size)))),-3])
}
##############################################################################################################
## Defining a vectorized subfunction to get adjacency list for a single vertice in product graph
##
## To compare ot to earlier loops, In this function i1 and j1 are fixed and the adjacency listed for all
## i2's and j2's inside this single function.
## This function use 'expand.grid' function to generate all possible i2,j2 pairs and a vectorize substraction
## to get the edge matches from a single vertice in product graph.
## Parameter b_a is combination of b and a where
## b is single row from the distl_in matrix. Means i1 is fixed
## a is single row from the distl_ref matrix. Means j1 is fixed
## The a length is also cut accordingly to reduce computations by so as to aoid calculations of upper triangulation
## of adjacency matrix
##############################################################################################################
get_edge_vertice <- function(b_a,eps,la,lb)
{
grid<- expand.grid(unlist(b_a[1]),unlist(b_a[2]))
edge_match <- which(abs(grid[,1]-grid[,2])<eps)
if(length(edge_match)==0) return(NULL)
return(edge_match)
}
##############################################################################################################
## Function to return stats on largest found clique
##############################################################################################################
#' @import vec2dtransf
#' @import sp
#' @import graphics
#' @import grDevices
#' @importFrom stats dist
#' @importFrom stats median
get_clique_stats <- function(clique,c_in,circle_in,c_ref,circle_ref,la,lb,plot=TRUE)
{
if (length(clique) < 3) {
cat("Cannot calculate clique, too few points\n")
return(NA)
}
## Subset clique points from input and refrence data
cq_idx <- as.data.frame(t(sapply(clique,function(z) c(i=ifelse(z %% lb==0,(z %/% lb),(z %/% lb)+1),j=ifelse(z %% lb==0,lb,z %% lb)))))
c_in_cq <- c_in[cq_idx$i,]
c_ref_cq <- c_ref[cq_idx$j,]
## Calculate affine transformation parameters i.e Rotation and Translation
affine_tfr_mat <- AffineTransformation(as.data.frame(cbind(c_in_cq,c_ref_cq)))
calculateParameters(affine_tfr_mat)
## Apply affine transformation to input circle and convert classes to spatial points
circle_in_sps <- applyTransformation(affine_tfr_mat,SpatialPoints(circle_in))
c_in_cq_sps <- applyTransformation(affine_tfr_mat,SpatialPoints(c_in_cq))
circle_ref_sps <- SpatialPoints(circle_ref)
## Prepare the match return statistics
circle_in_ref_dist <- apply(spDists(circle_in_sps,circle_ref_sps),1,min)
circle_in_dist_idx <- circle_in_sps@coords[circle_in_ref_dist<2,]
circle_in_ref_dist_idx <- (circle_ref_sps[apply(spDists(circle_in_sps,circle_ref_sps),1,which.min)]@coords)[circle_in_ref_dist<2,]
center_new <- apply(circle_in_ref_dist_idx,2,function(x) ((min(x) + max(x))/2))
radius_new <- max(sqrt((circle_in_ref_dist_idx[,1]-center_new[1])^2+(circle_in_ref_dist_idx[,2]-center_new[2])^2))
## Rotation Angle
R= matrix(affine_tfr_mat@parameters[-c(3,6)],2,byrow=FALSE)
out = svd(R)
s1=sign(out$u[1,1]*out$u[2,2])
s2=sign(out$v[1,1]*out$v[2,2])
out$u = out$u%*%diag(c(s1,1))
out$v = out$v%*%diag(c(s2,1))
out$d = out$d%*%diag(c(s1*s2,1))
theta = atan(out$u[2,1]/out$u[1,1])*180/pi
phi = atan(out$v[1,2]/out$v[1,1])*180/pi
rot_angle = min(abs(theta+phi),180-abs(theta+phi))
if(plot)
{
## 2*2 Plot layout
par(mfrow=c(2,2))
## Plot distance
plot(as.matrix(dist(c_in_cq)),as.matrix(dist(c_ref_cq)),xlab="",ylab="",main="Corresponding Distances Plot")
## Plot matching points
plot(c_ref_cq,col="red",pch=19,main="Matching Points")
points(c_in_cq_sps,col="blue")
## Plot all points
plot( circle_ref,col="red",
xlim=c(min(circle_ref$X,circle_in_sps@coords[,1]),max(circle_ref$X,circle_in_sps@coords[,1])),
ylim=c(min(circle_ref$Y,circle_in_sps@coords[,2]),max(circle_ref$Y,circle_in_sps@coords[,2])),
pch=19,main="All Points"
)
points(circle_in_sps,col="blue",pch=19)
## Plot Close points
plot( circle_in_ref_dist_idx,col="red",pch=19,
xlim=c(min(circle_in_dist_idx[,1],circle_in_ref_dist_idx[,1]),max(circle_in_dist_idx[,1],circle_in_ref_dist_idx[,1])),
ylim=c(min(circle_in_dist_idx[,2],circle_in_ref_dist_idx[,2]),max(circle_in_dist_idx[,2],circle_in_ref_dist_idx[,2])),
main="Close Points"
)
points(circle_in_dist_idx,col="blue",pch=19)
}
gc()
## Return Statistics
mylist <- list(clique_stats = data.frame( clique_size=length(clique),
rotation_angle = rot_angle,
reference_overlap=nrow(circle_in_ref_dist_idx)/nrow(circle_ref),
input_overlap=nrow(circle_in_dist_idx)/nrow(circle_in),
med_dist_euc = round(median(apply((circle_in_ref_dist_idx-circle_in_dist_idx)^2,1,sum)),5),
new_center_x=center_new[1],new_center_y=center_new[2],
new_radius=radius_new),
affine = affine_tfr_mat)
return(mylist)
}
get_os <- function() {
if (grepl("darwin", version$os)) {
return("mac64")
} else if (grepl("linux", version$os)) {
return("lin64")
} else {
return("win64")
}
}
#' Function to perform boosted clique on two input circles
#'
#' @name boosted_clique
#' @param circle_in The input circle
#' @param circle_ref The reference circle
#' @param ncross_in_bins Number of bins in the input circle (See smart_sample)
#' @param xbins_in Number of bins along each axis in the hexbin grid for the input circle
#' @param ncross_in_bin_size Number of points to sample from each bin in the input circle
#' @param ncross_ref_bins Number of bins in the reference circle
#' @param xbins_ref Number of bins along each axis in the hexbin grid for the reference circle
#' @param ncross_ref_bin_size Number of points to sample from each bin in the reference circle
#' @param eps Distance tolerance for declaring an edge match
#' @param seed The random seed for reproducing results
#' @param num_cores The number of processor cores for parallel processing
#' @param plot If TRUE, produce a plot of the clique results
#' @param verbose If TRUE, print out the timing results for each portion of the algorithm
#' @param cl Optionally, a parallel cluster that has already been created
#' @return The statistics for the matching between the two circles
#' @export
#' @import parallel
#' @importFrom stats dist
boosted_clique <- function(circle_in, circle_ref, ncross_in_bins = 30, xbins_in = 20, ncross_in_bin_size = 1, ncross_ref_bins = NULL, xbins_ref = 30, ncross_ref_bin_size = NULL, eps = .75, seed = 1, num_cores = 8, plot = TRUE, verbose = FALSE, cl = NULL) {
start_time_i <- Sys.time()
if(verbose) cat("Preparing circles data for edge matching.\n")
##############################################################################################################
## Sample data using function smart_sample (Not in use as to keep the seed and results same to soyoung's code)
##############################################################################################################
## Smart sample input circle points
c_in <- smart_sample(circle_in,sample_bins=min(ncross_in_bins,nrow(circle_in)),xbins=xbins_in,sample_size=ncross_in_bin_size,seed=seed)
## Smart sample reference circle points or full points
c_ref <- smart_sample(circle_ref,sample_bins=min(ncross_ref_bins,nrow(circle_ref)),xbins=xbins_ref,sample_size=ncross_ref_bin_size,seed=seed)
# set.seed(seed)
# c_in <- as.matrix(circle_in[sample(1:nrow(circle_in),ncross_in_bins),])
# c_ref <- as.matrix(circle_ref[sample(1:nrow(circle_ref),ifelse(is.null(ncross_ref_bins),nrow(circle_ref),ncross_ref_bins)),])
## Data dimesion for later reference within function
la <- nrow(c_in) ## Number of minutiae in image 1
lb <- nrow(c_ref) ## Number of minutiae in image 2
lg <- la * lb ## Number of vertices in the "product graph"
##############################################################################################################
##############################################################################################################
## Prepare distance matrix for adjancey lists calculations
##############################################################################################################
## Pairwise distance matrix between vertices in input and reference data
dist_in <- as.matrix(dist(c_in))
dist_ref <- as.matrix(dist(c_ref))
## Fill NA's to avoid loops
dist_in[row(dist_in)==col(dist_in)] <- NA
dist_ref[row(dist_ref)==col(dist_ref)] <- NA
## Split matrix to list of columns so as to vectorize calculation of adjacency list
distl_in <- as.list(as.data.frame(dist_in))
distl_ref <- as.list(as.data.frame(dist_ref))
dist_l_grid <- expand.grid(distl_ref,distl_in)
## Adding k1 verice 1 name to list so as to get edge match for k2<k1 only and thus avoiding values
## in upper triangulation matrix (Visualising edge match matrix lg*lg)
dist_l_grid$node <- 1:nrow(dist_l_grid)
dist_l_grid$Var2_n <- mapply(function(x,y) unlist(x)[1:unlist(y)],dist_l_grid[,2],as.list(((1:nrow(dist_l_grid)) %/% lb)+1))
dist_l_grid <- dist_l_grid[,c("Var1","Var2_n","node")]
rownames(dist_l_grid) <- 1:nrow(dist_l_grid)
gc()
##############################################################################################################
end_time_i <- Sys.time()
if(verbose) cat(paste("Prepared circles data for edge matching. Took",round(difftime(end_time_i,start_time_i,units="secs")),"Seconds.\n\n\n"))
start_time_al <- Sys.time()
if(verbose) cat("Calculating adjacency list.\n")
##############################################################################################################
## Adjacency list prepasration (Vectorized Code with parallel processing [*nix style])
##############################################################################################################
## Split data to distribute work to cores
dist_l_grid_core_assign <- split(dist_l_grid,floor(seq(0,((num_cores*4)-.01),length.out=nrow(dist_l_grid))))
## Windows
if (!is.null(cl)) {
final_ks <- do.call(c,parLapply(cl, dist_l_grid_core_assign,function(x,eps,la,lb) {gc();apply(x,1,get_edge_vertice,eps=eps,la=la,lb=lb)},eps=eps,la=la,lb=lb))
## Mac/Linux
} else {
final_ks <- do.call(c,mclapply(dist_l_grid_core_assign,function(x,eps,la,lb) {gc();apply(x,1,get_edge_vertice,eps=eps,la=la,lb=lb)},eps=eps,la=la,lb=lb,mc.cores=num_cores))
}
gc()
##############################################################################################################
end_time_al <- Sys.time()
if(verbose) cat(paste("Calculated adjacency list. Took",round(difftime(end_time_al,start_time_al,units="secs")),"Seconds.\n\n\n"))
start_time_mc <- Sys.time()
if(verbose) cat("Finding largest clique.\n")
##############################################################################################################
## Generate graph and get largest cliques
##############################################################################################################
## Prune edge list and write it to disk for parallel clique calculation
fro_node_name <- rep(as.numeric(sapply(names(final_ks),function(x) strsplit(x,"\\.")[[1]][2],USE.NAMES=FALSE)),sapply(final_ks,length))
to_node_name <- unlist(final_ks,use.names=FALSE)
edge_list <- paste(fro_node_name,to_node_name," ")[to_node_name<fro_node_name]
edge_list[length(edge_list)] <- gsub(" "," ",edge_list[length(edge_list)])
edge_list_pmc_format <- c( "%%MatrixMarket matrix coordinate pattern symmetric ",
paste(nrow(dist_l_grid),nrow(dist_l_grid),length(edge_list)," "),
edge_list
)
mytempdir <- tempdir()
writeLines(edge_list_pmc_format, file.path(mytempdir, "edge.mtx"))
ext <- ifelse(get_os() == "win64", ".exe", "")
clique_max <- system(paste0(system.file(package = "shoeprintr", "bin", get_os(), paste0("pmc", ext)), " -f ", file.path(mytempdir, "edge.mtx"), " -a 0"), intern = TRUE)
clique_max <- as.numeric(unlist(strsplit(gsub("Maximum clique: ","",clique_max[grepl("^Maximum clique: ",clique_max)])," ")))
## Cleanup
gc()
##############################################################################################################
end_time_mc <- Sys.time()
if(verbose) cat(paste("Calculated largest clique. Took",round(difftime(end_time_mc,start_time_mc,units="secs")),"Seconds.\n\n\n"))
##############################################################################################################
## Get and return clique stats
##############################################################################################################
return(get_clique_stats(clique_max,c_in,circle_in,c_ref,circle_ref,la=la,lb=lb,plot=plot))
##############################################################################################################
}
#' Adjust it to (0,0) in the left lower corner
#'
#' @name leftcorner_cent
#' @param data The input data
#' @export
leftcorner_cent <- function(data) {
return(apply((-1*data),2,function(x) x-min(x)))
}
##############################################################################################################
## Generate circle on given x ratio,y ratio and radius percentage
##############################################################################################################
get_circle <- function(data,x_rt,y_rt,rad_pct)
{
cir_center <- apply(data,2,function(x) max(x)-min(x))*c(x_rt,y_rt)
cir_rad <- min(apply(data,2,function(x) ((max(x)-min(x))/2)*rad_pct))
return(as.data.frame(data[apply(data,1,function(x,cir_center,cir_rad) sqrt(sum((x-cir_center)^2))<cir_rad,cir_center,cir_rad),]))
}
##############################################################################################################
## Generate circle on given x,y and radius percentage
##############################################################################################################
get_circle_fix <- function(x,y,rad_pct,data)
{
cir_center <- c(x,y)
if (rad_pct <= 1) {
cir_rad <- min(apply(data,2,function(x) ((max(x)-min(x))/2)*rad_pct))
} else {
cir_rad <- rad_pct
}
return(as.data.frame(data[apply(data,1,function(x,cir_center,cir_rad) sqrt(sum((x-cir_center)^2))<cir_rad,cir_center,cir_rad),]))
}
##############################################################################################################
## Generate circle on given x,y and radius
##############################################################################################################
get_circle_rad_fix <- function(x,y,cir_rad,data)
{
cir_center <- c(x,y)
return(as.data.frame(data[apply(data,1,function(x,cir_center,cir_rad) sqrt(sum((x-cir_center)^2))<cir_rad,cir_center,cir_rad),]))
}
#' @title Matches an input and a reference shoeprint
#'
#' @description Function to perform matching of footprints by matching all the candidate circles with input circles. This function performs initial matching between 3 circles on the input print and 27 candidate circles on the reference shoeprint. After the initial cliques, it select 2 best circles for each input circle and performs reinforcement matching on them. In reinforcement matching, full points on reference circle are considered instead of user defined. The function returns the final circle parameters and the statistics of the triangle formed.
#'
#' @name match_print
#' @param print_in The input print
#' @param print_ref The reference print
#' @param circles_input The input circles, specified as a matrix with columns (x, y, rad), or NULL to automatically generate
#' @param circles_reference The reference circles, specified as a matrix with columns (x, y, rad), or NULL to automatically generate
#' @param ncross_in_bins Number of bins in the input circle (See smart_sample)
#' @param xbins_in Number of bins along each axis in the hexbin grid for the input circle
#' @param ncross_in_bin_size Number of points to sample from each bin in the input circle
#' @param ncross_ref_bins Number of bins in the reference circle
#' @param xbins_ref Number of bins along each axis in the hexbin grid for the reference circle
#' @param ncross_ref_bin_size Number of points to sample from each bin in the reference circle
#' @param max_rotation_angle The maximum rotation angle, in degrees, for inclusion in the best circle matches
#' @param eps Distance tolerance for declaring an edge match
#' @param seed The random seed for reproducing results
#' @param num_cores The number of processor cores for parallel processing
#' @param plot If TRUE, produce a plot of the clique results
#' @param verbose If TRUE, print out the timing results for each portion of the algorithm
#' @return The statistics for the matching between the two circles
#' @import dplyr
#' @import ggplot2
#' @importFrom gridExtra grid.arrange
#' @importFrom stats dist
#' @export
#' @examples
#' \dontrun{
#' data(input_example)
#' data(reference_example)
#'
#' ## Transform all the points to have (0,0) at lower left corner.
#' print_in <- leftcorner_cent(reference_example)
#' print_ref <- leftcorner_cent(reference_example)
#'
#' ## Perform Print Match
#' print_stats <- match_print(print_in, print_ref,
#' ncross_in_bins = 30, xbins_in = 20, ncross_in_bin_size = 1,
#' ncross_ref_bins = NULL, xbins_ref = 30, ncross_ref_bin_size = NULL,
#' eps = .75, seed = 1, num_cores = parallel::detectCores(),
#' plot = TRUE, verbose = FALSE)
#' print_stats
#' }
match_print <- function(print_in, print_ref, circles_input = NULL, circles_reference = NULL, ncross_in_bins = 30, xbins_in = 20, ncross_in_bin_size = 1, ncross_ref_bins = NULL, xbins_ref = 30, ncross_ref_bin_size = NULL, max_rotation_angle = 365, eps = .75, seed = 1, num_cores = 8, plot = FALSE, verbose = FALSE) {
##############################################################################################################
## Cut circles on input print and reference print to perform matching
## 3 on Input circle and 27 on Reference circle
##############################################################################################################
## Generate 1 circles on input shoeprint
if (is.null(circles_input)) {
circles_dims <- apply(print_in,2,function(x) max(x)-min(x))
circle_centers <- matrix(c( .25,.8, ## 1/4 on x, 4/5 on y
.25,.3, ## 1/4 on x, 1/5 on y
.7,.6 ## 3/4 on x, 1/2 on y
),3,byrow=TRUE,dimnames=list(NULL,c("x","y"))) %*% diag(circles_dims)
circle_centers <- cbind(circle_centers,c(.4,.4,.4))
} else {
circle_centers <- circles_input
}
circles_in <- apply(circle_centers,1,function(x,data) get_circle_fix(x[1],x[2],x[3],data),data=print_in)
## Generate candidate circles on reference shoeprint
if (is.null(circles_reference)) {
rad_pct <- .5
x_radius <- min(apply(print_ref,2,function(x) ((max(x)-min(x))/2)*rad_pct))
x_ratios <- c(.25,.5,.75)
y_ratios <- x_radius*(1:round((max(print_ref)-min(print_ref))/x_radius))/(max(print_ref)-min(print_ref))
xy_ratios <- expand.grid(x_ratios,y_ratios)
circles_ref <- apply(xy_ratios,1,function(xy_rt,rad_pct,data) get_circle(data,xy_rt[1],xy_rt[2],rad_pct),rad_pct=rad_pct,data=print_ref)
} else {
circles_ref <- apply(circles_reference,1,function(x,data) get_circle_fix(x[1],x[2],x[3],data),data=print_ref)
}
##############################################################################################################
##############################################################################################################
##############################################################################################################
## Perform initial matching on circles
## 3 on Input circle and ~27 on Reference circle
##############################################################################################################
match_grid <- expand.grid(circles_ref,circles_in)
## Use a cluster for parallel computation if running on windows
cl <- NULL
if (get_os() == "win64") cl <- makeCluster(num_cores, renice = 0)
## Loop over all circle to circle comparisons
match_result <- lapply(1:nrow(match_grid), function(match_idx) {
cat(paste("Matching circle pair",match_idx,"out of", nrow(match_grid), "\n"))
## Get the circles
circle_in <- match_grid[match_idx,2][[1]]
circle_ref <- match_grid[match_idx,1][[1]]
## Affine transformation requires 3 control points
if (nrow(circle_ref) > 2 && nrow(circle_in) > 2) {
myboost <- boosted_clique(circle_in, circle_ref,
ncross_in_bins=ncross_in_bins,xbins_in=xbins_in,ncross_in_bin_size=ncross_in_bin_size,
ncross_ref_bins=ncross_ref_bins,xbins_ref=xbins_ref,ncross_ref_bin_size=ncross_ref_bin_size,
eps=eps,seed=seed,num_cores=num_cores,plot=plot,verbose=verbose,cl=cl
)
return(myboost$clique_stats)
} else {
cat("Skipping circle pair", match_idx, "due to no points contained\n")
return(NA)
}
})
##############################################################################################################
##############################################################################################################
##############################################################################################################
## Select best 2 circles against each of input circle for reinforcement learning
## Best criterion is whichever have largest input circle overlap
##############################################################################################################
match_result <- do.call(rbind,match_result)
match_result$circle1 <- rep(1:length(circles_in), each = length(circles_ref))
match_result$circle2 <- rep(1:length(circles_ref), times = length(circles_in))
## Take only ones with a low enough rotation angle
match_result_best <- match_result %>%
filter(rotation_angle <= max_rotation_angle) %>%
group_by(circle1) %>%
arrange(desc(input_overlap)) %>%
slice(1:2)
#best_2_match <- tapply(match_result$input_overlap,rep(1:length(circles_in), each=length(circles_ref)),function(x) order(x,decreasing=TRUE)[1:2])
#best_2_match <- mapply(function(x,y,z) x+(y*z),best_2_match,(1:length(circles_in))-1,MoreArgs=list(z=length(circles_ref)),SIMPLIFY=FALSE)
#best_2_match_params <- match_result[unlist(best_2_match),c("new_center_x","new_center_y")]
best_2_match_params <- match_result_best %>%
ungroup() %>%
select(new_center_x, new_center_y) %>%
as.data.frame()
circles_ref_reinf <- apply(best_2_match_params,1,function(x,rad_pct,data) get_circle_fix(x[1],x[2],rad_pct,data),rad_pct=.6,data=print_ref )
circles_in_reinf <- rep(circles_in, length.out = nrow(match_result_best))
##############################################################################################################
##############################################################################################################
##############################################################################################################
## Perform reinforcement learning
##############################################################################################################
match_result_reinf <- lapply(1:length(circles_in_reinf), function(match_idx_reinf) {
cat(paste("Reinforcement matching circle pair",match_idx_reinf,"out of", length(circles_in_reinf), "\n"))
myboost <- boosted_clique(circle_in=circles_in_reinf[[match_idx_reinf]],circle_ref=circles_ref_reinf[[match_idx_reinf]],
ncross_in_bins=ncross_in_bins,xbins_in=xbins_in,ncross_in_bin_size=ncross_in_bin_size,ncross_ref_bins=NULL,xbins_ref=30,ncross_ref_bin_size=NULL,
eps=eps,seed=seed,num_cores=num_cores,plot=plot,verbose=verbose,cl=cl
)
return(myboost$clique_stats)
})
if (!is.null(cl)) stopCluster(cl)
##############################################################################################################
## Prepare plot and return output
##############################################################################################################
match_result_reinf <- do.call(rbind,match_result_reinf)
match_result_reinf$circle1 <- match_result_best$circle1
match_result_reinf$circle2 <- match_result_best$circle2
#best_match <- tapply(match_result_reinf$input_overlap,rep(1:length(circles_in), each=2),function(x) order(x,decreasing=TRUE)[1],simplify=FALSE)
#best_match <- mapply(function(x,y,z) x+(y*z),best_match,(1:length(circles_in))-1,MoreArgs=list(z=2),SIMPLIFY=FALSE)
#circles_match <- match_result_reinf[unlist(best_match),]
#best_match_params <- match_result_reinf[unlist(best_match),c("new_center_x","new_center_y","new_radius")]
reinf_result_best <- match_result_reinf %>%
#filter(rotation_angle <= max_rotation_angle) %>%
group_by(circle1) %>%
arrange(desc(input_overlap)) %>%
slice(1)
circles_match <- reinf_result_best %>%
ungroup() %>%
as.data.frame()
best_match_params <- circles_match %>%
select(new_center_x, new_center_y, new_radius)
circles_ref_out <- apply(best_match_params,1,function(x,data) get_circle_rad_fix(x[1],x[2],x[3],data),data=print_ref)
circles_match <- cbind(circle_centers[,1:2, drop = FALSE],circles_match[,c("new_center_x","new_center_y","new_radius","rotation_angle","input_overlap")])
names(circles_match) <- c("Fixed_circle_X","Fixed_circle_Y","Match_circle_X","Match_circle_Y","Match_circle_Radius","Rotation_angle","Input_circle_overlap_pct")
## Congurent Triangle output
Fixed_triangle <- as.matrix(dist(circles_match[,c("Fixed_circle_X","Fixed_circle_Y")]))
Match_Triangle <- as.matrix(dist(circles_match[,c("Match_circle_X","Match_circle_Y")]))
p1 <- ggplot(data = as.data.frame(print_in), aes(x = x, y = y)) +
geom_point() +
geom_point(data = circles_in[[1]], color = "red") +
theme_bw()
if (length(circles_in) > 1) p1 <- p1 + geom_point(data = circles_in[[2]], color = "yellow")
if (length(circles_in) > 2) p1 <- p1 + geom_point(data = circles_in[[3]], color = "green")
p2 <- ggplot(data = as.data.frame(print_ref), aes(x = x, y = y)) +
geom_point() +
geom_point(data = circles_ref_out[[1]], color = "red") +
theme_bw()
if (length(circles_ref_out) > 1) p2 <- p2 + geom_point(data = circles_ref_out[[2]], color = "yellow")
if (length(circles_ref_out) > 2) p2 <- p2 + geom_point(data = circles_ref_out[[3]], color = "green")
## Plot
#dev.off()
#par(mfrow=c(1,2))
#plot(print_in,pch=19)
#points(circles_in[[1]],col="red")
#points(circles_in[[2]],col="yellow")
#points(circles_in[[3]],col="green")
#plot(print_ref,pch=19)
#points(circles_ref_out[[1]],col="red")
#points(circles_ref_out[[2]],col="yellow")
#points(circles_ref_out[[3]],col="green")
try(grid.arrange(p1, p2, ncol = 2))
fix1 <- Fixed_triangle[row(Fixed_triangle)<col(Fixed_triangle)]
fix2 <- Match_Triangle[row(Match_Triangle)<col(Match_Triangle)]
if (length(fix1) == 0) fix1 <- NA
if (length(fix2) == 0) fix2 <- NA
old_result <- cbind( circles_match,
data.frame( Triangle_sides = c("A-B","A-C","B-C")[1:nrow(circles_match)],
Fixed_circle_side_length = fix1,
Match_circle_side_length = fix2
)
)
new_result<- sum_result(old_result)
## Return Stats
return(list(old_result, new_result))
}
sum_result<-function(data){
R1<-mean(data[,7])
R2<-mean(abs(data[,9]-data[,10]))
R3<-sd(data[,6])
Result<-c(R1,R2,R3)
return(Result)
}
#' @title Matches a center circle input and a reference shoeprint
#'
#' @description Function to perform matching of center footprints by matching all the candidate circles with input circles. This function performs initial matching between 3 circles on the input print and 27 candidate circles on the reference shoeprint. After the initial cliques, it select 2 best circles for each input circle and performs reinforcement matching on them. In reinforcement matching, full points on reference circle are considered instead of user defined. The function returns the final circle parameters and the statistics of the triangle formed.
#'
#' @name centercircle_match
#' @param input The input print
#' @param reference The reference print
#' @param output The output obtained from the match_print function
#'
#' @export
centercircle_match<-function(input, reference, output){
if (nrow(output) < 3) stop("Must have three circles in the output")
#center 1 - X
cx.in.1<-mean(c(output[1,1],output[2,1]))
#center 1 - Y
cy.in.1<-mean(c(output[1,2],output[2,2]))
#center 2 - X
cx.in.2<-mean(c(output[2,1],output[3,1]))
#center 2 - Y
cy.in.2<-mean(c(output[2,2],output[3,2]))
#center 3 - X
cx.in.3<-mean(c(output[3,1],output[1,1]))
#center 3 - Y
cy.in.3<-mean(c(output[3,2],output[1,2]))
#center 4 - X
cx.in.4<-mean(c(output[1,1],output[2,1],output[3,1]))
#center 4 - Y
cy.in.4<-mean(c(output[1,2],output[2,2],output[3,2]))
#estimated circle location
#center 1 - X
cx.re.1<-mean(c(output[1,3],output[2,3]))
#center 1 - Y
cy.re.1<-mean(c(output[1,4],output[2,4]))
#center 2 - X
cx.re.2<-mean(c(output[2,3],output[3,3]))
#center 2 - Y
cy.re.2<-mean(c(output[2,4],output[3,4]))
#center 3 - X
cx.re.3<-mean(c(output[3,3],output[1,3]))
#center 3 - Y
cy.re.3<-mean(c(output[3,4],output[1,4]))
#center 4 - X
cx.re.4<-mean(c(output[1,3],output[2,3],output[3,3]))
#center 4 - Y
cy.re.4<-mean(c(output[1,4],output[2,4],output[3,4]))
# step 2 : 4 of center circles matching
input_circles <- matrix(c(cx.in.1, cy.in.1, 40), nrow = 1, ncol = 3)
reference_circles <- matrix(c(cx.re.1, cy.re.1, 55), nrow = 1, ncol = 3)
cc1 <- match_print(input, reference, circles_input = input_circles, circles_reference = reference_circles,
ncross_in_bins = 30, xbins_in = 20, ncross_in_bin_size = 1,
ncross_ref_bins = NULL, xbins_ref = 30, ncross_ref_bin_size = NULL,
eps = .75, seed = 1, num_cores = parallel::detectCores(),
plot = FALSE, verbose = TRUE)
input_circles <- matrix(c(cx.in.2, cy.in.2, 40), nrow = 1, ncol = 3)
reference_circles <- matrix(c(cx.re.2, cy.re.2, 55), nrow = 1, ncol = 3)
cc2 <- match_print(input, reference, circles_input = input_circles, circles_reference = reference_circles,
ncross_in_bins = 30, xbins_in = 20, ncross_in_bin_size = 1,
ncross_ref_bins = NULL, xbins_ref = 30, ncross_ref_bin_size = NULL,
eps = .75, seed = 1, num_cores = parallel::detectCores(),
plot = FALSE, verbose = TRUE)
input_circles <- matrix(c(cx.in.3, cy.in.3, 40), nrow = 1, ncol = 3)
reference_circles <- matrix(c(cx.re.3, cy.re.3, 55), nrow = 1, ncol = 3)
cc3 <- match_print(input, reference, circles_input = input_circles, circles_reference = reference_circles,
ncross_in_bins = 30, xbins_in = 20, ncross_in_bin_size = 1,
ncross_ref_bins = NULL, xbins_ref = 30, ncross_ref_bin_size = NULL,
eps = .75, seed = 1, num_cores = parallel::detectCores(),
plot = FALSE, verbose = TRUE)
input_circles <- matrix(c(cx.in.4, cy.in.4, 40), nrow = 1, ncol = 3)
reference_circles <- matrix(c(cx.re.4, cy.re.4, 55), nrow = 1, ncol = 3)
cc4 <- match_print(input, reference, circles_input = input_circles, circles_reference = reference_circles,
ncross_in_bins = 30, xbins_in = 20, ncross_in_bin_size = 1,
ncross_ref_bins = NULL, xbins_ref = 30, ncross_ref_bin_size = NULL,
eps = .75, seed = 1, num_cores = parallel::detectCores(),
plot = FALSE, verbose = TRUE)
comp<-c('1-2','2-3','3-1','1-2-3')
re<-rbind(cc1[[1]][1,1:7],cc2[[1]][1,1:7],cc3[[1]][1,1:7],cc4[[1]][1,1:7])
Result<-data.frame(comp,re)
return(Result)
}
#' @title Starting plot to position three local areas of shoe Q
#'
#' @description Function to plot shoe Q and shoe K together. In shoe Q, three circles are colored.
#'
#' @name start_plot
#' @param input The input print (shoe Q)
#' @param reference The reference print (shoe K)
#' @param input_circle Centers and radius for three circles that we fix in the input print (shoe Q)
#'
#' @export
start_plot<-function(input, reference, input_circle){
cx1<-input_circle[1,1]
cx2<-input_circle[2,1]
cx3<-input_circle[3,1]
cy1<-input_circle[1,2]
cy2<-input_circle[2,2]
cy3<-input_circle[3,2]
r1<-input_circle[1,3]
r2<-input_circle[2,3]
r3<-input_circle[3,3]
P1<-ggplot(data.frame(input), aes(x=x, y=y))+ geom_point(data=data.frame(input), aes(x=x, y=y), color='black',size=0.1) +
geom_point(data=data.frame(int_inside_center(data.frame(input), r1, nseg, cx1,cy1)),color="red",size=0.1)+
gg_circle(r1, xc=cx1, yc=cy1, color="red") +
geom_point(data=data.frame(int_inside_center(data.frame(input), r2, nseg, cx2,cy2)),color="orange",size=0.1)+
gg_circle(r2, xc=cx2, yc=cy2, color="orange") +
geom_point(data=data.frame(int_inside_center(data.frame(input), r3, nseg, cx3, cy3)),color="green",size=0.1)+
gg_circle(r3, xc=cx3, yc=cy3, color="green")
P2<-ggplot(data.frame(reference), aes(x=x, y=y))+ geom_point(data=data.frame(reference), aes(x=x, y=y), color='black',size=0.1)
return(P1+P2)
}
#' @title Detecting points within the circle
#'
#' @description Function to find points within the circle
#'
#' @name int_inside_center
#' @param Data The input print
#' @param r radius of the circle
#' @param c1 x value of the circle
#' @param c2 y value of the circle
#'
#' @export
int_inside_center<-function(Data, r, nseg, c1,c2){
nseg=360
x <- c1+r*cos(seq(0,2*pi, length.out=nseg))
y <- c2+r*sin(seq(0,2*pi, length.out=nseg))
circle<-data.frame(x,y)
shoe_cc<-subset(Data, (c1-r)<Data[,1] & Data[,1]<(c1+r) & (c2-r)<Data[,2] & Data[,2]<(c2+r))
intersect_points<-NULL
CIR<-(shoe_cc[,1]-c1)^2+(shoe_cc[,2]-c2)^2
shoe_cc<-cbind(shoe_cc, CIR)
inside_cir<-subset(shoe_cc, shoe_cc[,3]<r^2)
x<-inside_cir[,1]
y<-inside_cir[,2]
int_pts_cir<-cbind(x,y)
return(unique(int_pts_cir))
}
#' @title Adding the circlular area
#'
#' @description Function to find the circular area
#'
#' @name gg_circle
#' @param r radius of the circle
#' @param xc x value of the circle
#' @param yc y value of the circle
#'
#' @export
#'
gg_circle<-function(r, xc, yc, color="black", fill=NA, ...) {
x <- xc + r*cos(seq(0, pi, length.out=100))
ymax <- yc + r*sin(seq(0, pi, length.out=100))
ymin <- yc + r*sin(seq(0, -pi, length.out=100))
annotate("ribbon", x=x, ymin=ymin, ymax=ymax, color=color, fill=fill, ...)
}
#' @title Find the initial circle position by quantile of x and y ranges
#'
#' @description Function to find the initial circle position by quantile of x and y ranges. Centers are found as (30%, 80%), (20%, 40%), (70%, 70%) with radius of 50 for three circles.
#'
#' @name initial_circle
#' @param input input impression
#'
#' @export
#'
#'
initial_circle<-function(input){
circles_dims <- apply(input, 2, function(x) max(x) - min(x))
circle_centers1 <- matrix(c(0.3, 0.75, 0.2, 0.25, 0.75, 0.65), 3,
byrow = TRUE, dimnames = list(NULL, c("x", "y"))) %*% diag(circles_dims)
circle_centers2 <- circle_centers1 + matrix(c(min(input[,1]), min(input[,1]), min(input[,1]), min(input[,2]),min(input[,2]),min(input[,2])), 3, byrow=FALSE)
circle_centers3 <- cbind(circle_centers2, c(50, 50, 50))
input_circles<-circle_centers3
return(input_circles)
}
#' @title Cut x and y coordinate points less than 1% quantile and larger than 99% quantile
#'
#' @description Function to cut x and y coordinate points less than 1% quantile and larger than 99% quantile
#'
#' @name focus_data2
#' @param data data with x and y coordinate values
#'
#' @export
#'
#'
focus_data2<-function(data){
data.new<-subset(data, quantile(data[,1], 0.01)<data[,1] &
data[,1]<quantile(data[,1], 0.99) &
quantile(data[,2], 0.01)<data[,2] &
data[,2]<quantile(data[,2], 0.99))
data.new<-data.frame(data.new)
names(data.new)<-c("x","y")
return(data.new)
}
#' @title Reverse right side shoe to look like left side shoe
#'
#' @description Function to transform coordinate values to look like from the left side shoe
#'
#' @name rev_R
#' @param Data data with x and y coordinate values from right side of shoe
#'
#' @export
#'
rev_R<-function(Data){
newx<-max(Data[,1])-Data[,1]
rev_Data<-data.frame(newx, Data[,2])
names(rev_Data)<-c("x", "y")
return(rev_Data)
}
#' @title Draw multiple plots
#'
#' @description Draw multiple plots
#'
#' @name multiplot
#'
#' @export
#'
multiplot<-function(..., plotlist=NULL, file, cols=1, layout=NULL) {
library(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
#' @title First plot before starting matching to confirm the circle position.
#'
#' @description Draw edge coordiates with designated circles. Circle will be colored in both input and reference
#'
#' @name step0_plot
#' @param input Input image for fixing three circles
#' @param reference Input image for finding correspondence
#' @param input_circles Three circles which will be fixed in the input
#'
#' @export
#'
step0_plot<-function(input, reference, input_circles){
cx1<-input_circles[1,1]
cx2<-input_circles[2,1]
cx3<-input_circles[3,1]
cy1<-input_circles[1,2]
cy2<-input_circles[2,2]
cy3<-input_circles[3,2]
r1<-input_circles[1,3]
r2<-input_circles[2,3]
r3<-input_circles[3,3]
P1<-ggplot(data.frame(input), aes(x=x, y=y))+ geom_point(data=data.frame(input), aes(x=x, y=y), color='black',size=0.1) +
geom_point(data=data.frame(int_inside_center(data.frame(input), r1, nseg, cx1,cy1)),color="red",size=0.1)+
gg_circle(r1, xc=cx1, yc=cy1, color="red") +
geom_point(data=data.frame(int_inside_center(data.frame(input), r2, nseg, cx2,cy2)),color="orange",size=0.1)+
gg_circle(r2, xc=cx2, yc=cy2, color="orange") +
geom_point(data=data.frame(int_inside_center(data.frame(input), r3, nseg, cx3, cy3)),color="green",size=0.1)+
gg_circle(r3, xc=cx3, yc=cy3, color="green")
P2<-ggplot(data.frame(reference), aes(x=x, y=y))+ geom_point(data=data.frame(reference), aes(x=x, y=y), color='black',size=0.1) +
geom_point(data=data.frame(int_inside_center(data.frame(reference), r1, nseg, cx1,cy1)),color="red",size=0.1)+
gg_circle(r1, xc=cx1, yc=cy1, color="red") +
geom_point(data=data.frame(int_inside_center(data.frame(reference), r2, nseg, cx2,cy2)),color="orange",size=0.1)+
gg_circle(r2, xc=cx2, yc=cy2, color="orange") +
geom_point(data=data.frame(int_inside_center(data.frame(reference), r3, nseg, cx3,cy3)),color="green",size=0.1)+
gg_circle(r3, xc=cx3, yc=cy3, color="green")
return(multiplot(P1,P2,cols=2))
}
#' @title Match input and reference with given circle information in input. The searching area is confined.
#'
#' @description Find corresponding areas in reference for fixed circles in input. To reduce the time, it confines the area for candidate circles in reference image.
#' For the first circle, it confines the area into left toe, for the second circle, it searches only left bottom, for the third circle, it searches right toe area of reference shoe.
#'
#' @name match_print_subarea
#' @param input Input image for fixing three circles
#' @param reference Input image for finding correspondence
#' @param input_circles Three circles which will be fixed in the input
#' @param max_rotation_angle maximum rotation angle that we allow for comparisons.
#'
#' @export
#'
match_print_subarea<-function(input, reference, input_circles, max_rotation_angle){
#original : c(0.3, 0.85, 0.25, 0.3, 0.75, 0.75)
if (is.null(input_circles)) {
circles_dims <- apply(input, 2, function(x) max(x) - min(x))
circle_centers1 <- matrix( c(0.3, 0.85, 0.25, 0.3, 0.75, 0.75), 3, byrow = TRUE, dimnames = list(NULL, c("x", "y"))) %*% diag(circles_dims)
circle_centers2 <- circle_centers1 + matrix(c(min(input[,1]), min(input[,1]), min(input[,1]), min(input[,2]),min(input[,2]),min(input[,2])), 3, byrow=FALSE)
circle_centers3 <- cbind(circle_centers2, c(50, 50, 50))
input_circles<-circle_centers3
} else {
input_circles <- input_circles
}
ref<-data.frame(reference)
ref_len_y<-(max(reference[,2])-min(reference[,2]))
ref_len_x<-(max(reference[,1])-min(reference[,1]))
ref_top<-subset(ref, y>0.5*ref_len_y+min(reference[,2]))
ref_bottom<-subset(ref, y<=0.5*ref_len_y+min(reference[,2]))
location_ref<-list()
location_ref[[1]]<-subset(ref_top, x<(ref_len_x/2+min(reference[,1]))) # ref_top_left
location_ref[[2]]<-subset(ref_bottom, x<(ref_len_x/2+min(reference[,1]))) #ref_bottom_left
location_ref[[3]]<-subset(ref_top, x>=(ref_len_x/2+min(reference[,1]))) #ref_top_right
location_ref[[4]]<-subset(ref_bottom, x>=(ref_len_x/2+min(reference[,1]))) #ref_bottom_right
nseg=360
in.cx<-NULL
in.cy<-NULL
in.r<-NULL
final.cx<-NULL
final.cy<-NULL
final.r<-NULL
MM<-NULL
FM<-NULL
ref_loc<-NULL
rd_score<-NULL
WM<-NULL
ref_loc<-NULL
rd_score<-NULL
wrong_mat<-NULL
WW<-NULL
wrong_set<-NULL
K<-matrix(c(2,3,4,1,3,4,1,2,4,1,2,3), nrow=4, byrow=T)
for ( k in 1:3){
print(paste("circle",k,"matching"))
in.cx[k]<-input_circles[k,1]
in.cy[k]<-input_circles[k,2]
in.r[k]<-input_circles[k,3]
circle_in<-data.frame(int_inside_center(data.frame(input), in.r[k], nseg, in.cx[k], in.cy[k]))
R<-NULL
r_ref<-(input_circles[1,3]+15)
ref_loc<-location_ref[[k]]
ref_loc_len_x<-(max(ref_loc$x)-min(ref_loc$x))
ref_loc_len_y<-(max(ref_loc$y)-min(ref_loc$y))
if(min(ref_loc$y)+r_ref<max(ref_loc$y)-r_ref){
ref_cdd_y<-seq(min(ref_loc$y)+r_ref, max(ref_loc$y), by = r_ref)} else {
ref_cdd_y<-seq(min(ref_loc$y)+r_ref, max(ref_loc$y)+r_ref, by = r_ref)
}
ref_cdd_x<-(min(ref_loc$x)+max(ref_loc$x))/2
for ( j in 1:length(ref_cdd_y)){
circle_cdd_ref<-data.frame(int_inside_center(data.frame(ref), r_ref, nseg, ref_cdd_x, ref_cdd_y[j]))
if(nrow(circle_cdd_ref)>(nrow(circle_in)*0.2)){
M<-try(boosted_clique(circle_in, circle_cdd_ref, ncross_in_bins = 30, xbins_in = 20,
ncross_in_bin_size = 1, ncross_ref_bins = NULL, xbins_ref = 30,
ncross_ref_bin_size = NULL, eps = 0.75, seed = 1, num_cores = parallel::detectCores()-1,
plot = FALSE, verbose = FALSE, cl = NULL))
if (sum(is.na(M[[1]]))<1) {
M2<-M$clique_stats
Mat<-paste0('mat',j)
R<-rbind(R,cbind(Mat,M2))}
}
}
R1<-subset(R,rotation_angle<max_rotation_angle)
if (nrow(R1)!=0){
new.cx<-R1[which.max(R1$input_overlap),7]
new.cy<-R1[which.max(R1$input_overlap),8]
new.r<-R1[which.max(R1$input_overlap),9]+15} else{
new.cx<-R[which.max(R$input_overlap),7]
new.cy<-R[which.max(R$input_overlap),8]
new.r<-R[which.max(R$input_overlap),9]+15
}
circle_ref<-data.frame(int_inside_center(data.frame(reference), new.r, nseg, new.cx, new.cy))
step2_mat<-boosted_clique(circle_in, circle_ref, ncross_in_bins = 30, xbins_in = 20,
ncross_in_bin_size = 1, ncross_ref_bins = NULL, xbins_ref = 30,
ncross_ref_bin_size = NULL, eps = 0.75, seed = 1, num_cores = parallel::detectCores()-1,
plot = FALSE, verbose = FALSE, cl = NULL)$clique_stats
final.cx[k]<-step2_mat[,6]
final.cy[k]<-step2_mat[,7]
final.r[k]<-step2_mat[,8]
MM<-rbind(MM,step2_mat)
#wrong_set<-rbind(location_ref[[K[k,1]]],location_ref[[K[k,2]]],location_ref[[K[k,3]]] )
#count<-0
#repeat{
# w_idx<-sample(1:nrow(wrong_set),1)
# rd.cx<-as.numeric(wrong_set[w_idx,][1])
# rd.cy<-as.numeric(wrong_set[w_idx,][2])
# new.r<-55
# rd_circle_ref<-data.frame(int_inside_center(data.frame(reference), new.r, nseg, rd.cx, rd.cy))
# count<-count+1
# if(nrow(rd_circle_ref)>0.3*nrow(circle_in) & nrow(rd_circle_ref)<1.3*nrow(circle_in) | count==20){
# break}
#}
#wrong_mat<-boosted_clique(circle_in, rd_circle_ref, ncross_in_bins = 30, xbins_in = 20,
#ncross_in_bin_size = 1, ncross_ref_bins = NULL, xbins_ref = 30,
#ncross_ref_bin_size = NULL, eps = 0.75, seed = 1, num_cores = parallel::detectCores(),
#plot = FALSE, verbose = FALSE, cl = NULL)$clique_stats
#WW<-rbind(WW,wrong_mat)
}
Input_X<-input_circles[,1]
Input_Y<-input_circles[,2]
Comp<-c('1-2','1-3','2-3')
d_in_1<-sqrt((Input_X[1]-Input_X[2])^2+(Input_Y[1]-Input_Y[2])^2)
d_in_2<-sqrt((Input_X[1]-Input_X[3])^2+(Input_Y[1]-Input_Y[3])^2)
d_in_3<-sqrt((Input_X[2]-Input_X[3])^2+(Input_Y[2]-Input_Y[3])^2)
Euc_input_dist<-c(d_in_1,d_in_2,d_in_3)
Reference_X<-MM[,6]
Reference_Y<-MM[,7]
d_ref_1<-sqrt((Reference_X[1]-Reference_X[2])^2+(Reference_Y[1]-Reference_Y[2])^2)
d_ref_2<-sqrt((Reference_X[1]-Reference_X[3])^2+(Reference_Y[1]-Reference_Y[3])^2)
d_ref_3<-sqrt((Reference_X[2]-Reference_X[3])^2+(Reference_Y[2]-Reference_Y[3])^2)
Euc_ref_dist<-c(d_ref_1,d_ref_2,d_ref_3)
Reference_radius<-MM[,8]
#FM<-data.frame(Input_X,Input_Y,Reference_X,Reference_Y,Reference_radius,MM[,c(1:5)],Comp,Euc_input_dist,Euc_ref_dist, WW)
FM<-data.frame(Input_X,Input_Y,Reference_X,Reference_Y,Reference_radius,MM[,c(1:5)],Comp,Euc_input_dist,Euc_ref_dist)
P1<-ggplot(data.frame(input), aes(x=x, y=y))+ geom_point(data=data.frame(input), aes(x=x, y=y), color='black',size=0.1) +
geom_point(data=data.frame(int_inside_center(data.frame(input), in.r[1], nseg, in.cx[1], in.cy[1])),color="red",size=0.1)+
gg_circle(in.r[1], xc=in.cx[1], yc=in.cy[1], color="red") +
geom_point(data=data.frame(int_inside_center(data.frame(input), in.r[2], nseg, in.cx[2], in.cy[2])),color="yellow",size=0.1)+
gg_circle(in.r[2], xc=in.cx[2], yc=in.cy[2], color="yellow") +
geom_point(data=data.frame(int_inside_center(data.frame(input), in.r[3], nseg, in.cx[3], in.cy[3])),color="green",size=0.1)+
gg_circle(in.r[3], xc=in.cx[3], yc=in.cy[3], color="green")
P2<-ggplot(data.frame(reference), aes(x=x, y=y))+ geom_point(data=data.frame(reference), aes(x=x, y=y), color='black',size=0.1) +
geom_point(data=data.frame(int_inside_center(data.frame(reference), final.r[1], nseg, final.cx[1],final.cy[1])),color="red",size=0.1)+
gg_circle(final.r[1], xc=final.cx[1], yc=final.cy[1], color="red") +
geom_point(data=data.frame(int_inside_center(data.frame(reference), final.r[2], nseg, final.cx[2],final.cy[2])),color="yellow",size=0.1)+
gg_circle(final.r[2], xc=final.cx[2], yc=final.cy[2], color="yellow") +
geom_point(data=data.frame(int_inside_center(data.frame(reference), final.r[3], nseg, final.cx[3],final.cy[3])),color="green",size=0.1)+
gg_circle(final.r[3], xc=final.cx[3], yc=final.cy[3], color="green")
try(multiplot(P1, P2, cols=2))
return(FM)
}
#' @title Calculate sensitivity and specificity at different threshold for drawing ROC curves
#'
#' @description Using prediction by some methods, calculate sensitivity and specificity at different threshold for drawing ROC curves.
#'
#' @name ROC.data_col
#' @param nsame number of mates
#' @param ndiff number of non-mates
#' @param thres Threshold for ROC curve
#' @param Method method name
#'
#' @export
#'
ROC.data_col<-function(pred_prob, nsame, ndiff, thres, Method){
Sen<-NULL
Spec<-NULL
one_minus_Spec<-NULL
one_minus_Sen<-NULL
ntotal=nsame+ndiff
for (i in 1:length(thres)){
A<-sum(pred_prob[1:nsame]>=thres[i]) # true same - pred same : A
C<-sum(pred_prob[1:nsame]<thres[i]) # true same - pred diff : C
B<-sum(pred_prob[(nsame+1):ntotal]>=thres[i]) # true diff - pred same : B
D<-sum(pred_prob[(nsame+1):ntotal]<thres[i]) # true diff - pred diff : D
Sen[i]<-A/(A+C)
Spec[i]<-D/(B+D)
one_minus_Spec[i]<-(1-Spec[i])
one_minus_Sen[i]<-(1-Sen[i])
}
D=data.frame(thres, Sen=Sen, one_minus_Spec=one_minus_Spec, Spec=Spec, one_minus_Sen, Method=Method)
return(D)
}
#' @title Calculate AUC, EER, optimal threshold for ROC curves
#'
#' @description With the list type of data containing sensitivity and sepecificity values from many thresholds, it calcuates AUC, EER, optiimal threshold
#' to give the minumum summation of FPR and FNR.
#'
#' @name auc_eer_opt.t_generation
#' @param ROC_list List of sensitivity and sepecificity from many different methods
#'
#' @export
#'
auc_eer_opt.t_generation<-function(ROC_list){
auc<-list()
eer<-list()
opt.t<-list()
one_minus_Spec<-list()
one_minus_Sen<-list()
for (i in 1:length(ROC_list)){
auc[[i]]<-with(ROC_list[[i]], simple_auc(Sen, one_minus_Spec))
eer[[i]]<-c(eer.ROC(unique(ROC_list[[i]][, -1])))
opt.t[[i]]<-ROC_list[[i]]$thres[which.max(ROC_list[[i]]$Sen + ROC_list[[i]]$Spec)]
one_minus_Spec[[i]]<-ROC_list[[i]][which(ROC_list[[i]]$thres == opt.t[[i]]),]$one_minus_Spec
one_minus_Sen[[i]]<-ROC_list[[i]][which(ROC_list[[i]]$thres == opt.t[[i]]),]$one_minus_Sen
}
auc_combined<-c(do.call(rbind, auc))
eer_combined<-c(do.call(rbind, eer))
opt.t_combined<-c(do.call(rbind, opt.t))
one_minus_Spec_c<-c(do.call(rbind, one_minus_Spec))
one_minus_Sen_c<-c(do.call(rbind, one_minus_Sen))
name<-c("RF-plus-POC-R","RF", "% Overlap on Q","POC-R","POC-P","FMTC")
auc_eer_data_summary<-data.frame(name, auc_combined,eer_combined,opt.t_combined,one_minus_Spec_c,one_minus_Sen_c)
auc_eer_data_summary2<-auc_eer_data_summary[order(auc_eer_data_summary$auc_combined, decreasing = TRUE), ]
return(auc_eer_data_summary2)
}
#' @title Calculate AUC
#'
#' @description Calculate AUC with information of TPR and FPR
#'
#' @name simple_auc
#' @param TPR True positive rate
#' @param FPR False positive rate
#'
#' @export
#'
simple_auc<-function(TPR, FPR){
# inputs already sorted, best scores first
dFPR <- c(diff(FPR), 0)
dTPR <- c(diff(TPR), 0)
sum(TPR * dFPR) + sum(dTPR * dFPR)/2
}
#' @title Calculate EER
#'
#' @description Calculate EER with ROC data containing FPR and FNR
#'
#' @name eer.ROC
#' @param ROCdata True positive rate
#'
#' @export
#'
eer.ROC<-function(ROCdata){
ROC<-ROCdata
missrates <-ROC$one_minus_Sen
farates <- ROC$one_minus_Spec
# Find the point on the ROC with miss slightly >= fa, and the point
# next to it with miss slightly < fa.
dists <- missrates - farates;
idx1 <- which( dists == min( dists[ dists >= 0 ] ) );
idx2 <- which( dists == max( dists[ dists < 0 ] ) );
stopifnot( length( idx1 ) == 1 );
stopifnot( length( idx2 ) == 1 );
stopifnot( abs( idx1 - idx2 ) == 1 );
# Extract the two points as (x) and (y), and find the point on the
# line between x and y where the first and second elements of the
# vector are equal. Specifically, the line through x and y is:
# x + a*(y-x) for all a, and we want a such that
# x[1] + a*(y[1]-x[1]) = x[2] + a*(y[2]-x[2]) so
# a = (x[1] - x[2]) / (y[2]-x[2]-y[1]+x[1])
x <- c( missrates[idx1], farates[idx1] );
y <- c( missrates[idx2], farates[idx2] );
a <- ( x[1] - x[2] ) / ( y[2] - x[2] - y[1] + x[1] );
eer <- x[1] + a * ( y[1] - x[1] );
return(eer)}
#' @title Find three circles
#'
#' @description Find three circles with different positions for degraded impressions
#'
#' @name initial_circle2
#' @param input Questioned impression
#'
#' @export
#'
initial_circle2<-function(input){
circles_dims <- apply(input, 2, function(x) max(x) - min(x))
circle_centers1 <- matrix(c(0.3, 0.8, 0.2, 0.3, 0.7, 0.7), 3,
byrow = TRUE, dimnames = list(NULL, c("x", "y"))) %*% diag(circles_dims)
circle_centers2 <- circle_centers1 + matrix(c(min(input[,1]), min(input[,1]), min(input[,1]), min(input[,2]),min(input[,2]),min(input[,2])), 3, byrow=FALSE)
circle_centers3 <- cbind(circle_centers2, c(50, 50, 50))
input_circles<-circle_centers3
return(input_circles)
}
#' @title Draw plot of two shoe impressions with each designated circles
#'
#' @description Draw plot of questioned and reference shoe with each designated circles.
#'
#' @name step0_Q_K_plot
#' @param input Questioned impression (Q)
#' @param reference Known impression (K)
#' @param circle_in Circle information for Q
#' @param circle_ref Circle information for K
#'
#' @export
#'
step0_Q_K_plot<-function(input, reference, circle_in, circle_ref){
cx1<-circle_in[1,1]
cx2<-circle_in[2,1]
cx3<-circle_in[3,1]
cy1<-circle_in[1,2]
cy2<-circle_in[2,2]
cy3<-circle_in[3,2]
r1<-circle_in[1,3]
r2<-circle_in[2,3]
r3<-circle_in[3,3]
P1<-ggplot(data.frame(input), aes(x=x, y=y))+ geom_point(data=data.frame(input), aes(x=x, y=y), color='black',size=0.1) +
geom_point(data=data.frame(int_inside_center(data.frame(input), r1, nseg, cx1,cy1)),color="red",size=0.1)+
gg_circle(r1, xc=cx1, yc=cy1, color="red") +
geom_point(data=data.frame(int_inside_center(data.frame(input), r2, nseg, cx2,cy2)),color="orange",size=0.1)+
gg_circle(r2, xc=cx2, yc=cy2, color="orange") +
geom_point(data=data.frame(int_inside_center(data.frame(input), r3, nseg, cx3, cy3)),color="green",size=0.1)+
gg_circle(r3, xc=cx3, yc=cy3, color="green")
cx1<-circle_ref[1,1]
cx2<-circle_ref[2,1]
cx3<-circle_ref[3,1]
cy1<-circle_ref[1,2]
cy2<-circle_ref[2,2]
cy3<-circle_ref[3,2]
r1<-circle_ref[1,3]
r2<-circle_ref[2,3]
r3<-circle_ref[3,3]
P2<-ggplot(data.frame(reference), aes(x=x, y=y))+ geom_point(data=data.frame(reference), aes(x=x, y=y), color='black',size=0.1) +
geom_point(data=data.frame(int_inside_center(data.frame(reference), r1, nseg, cx1,cy1)),color="red",size=0.1)+
gg_circle(r1, xc=cx1, yc=cy1, color="red") +
geom_point(data=data.frame(int_inside_center(data.frame(reference), r2, nseg, cx2,cy2)),color="orange",size=0.1)+
gg_circle(r2, xc=cx2, yc=cy2, color="orange") +
geom_point(data=data.frame(int_inside_center(data.frame(reference), r3, nseg, cx3,cy3)),color="green",size=0.1)+
gg_circle(r3, xc=cx3, yc=cy3, color="green")
return(multiplot(P1,P2,cols=2))
}
#' @title Match input and reference with given circle information in input. Circle in the reference will be also fixed.
#'
#' @description Match between circle 1 in the input impression vs. circle 1 in the referece impression. Repeat this for circle 2 and circle 3.
#' This matching function is for comparing degraded impressions to see how worse the method perform when comparing similar regions.
#'
#'
#' @name edge_fix_circle_mat
#' @param input Input image for fixing three circles
#' @param reference Input image for finding correspondence
#' @param input_circles Three circles which will be fixed in the input
#' @param ref_circles Three circles which will be fixed in the reference too.
#'
#' @export
#'
edge_fix_circle_mat<-function(input, referece, input_circles, ref_circles){
MM<-NULL
in.cx<-NULL
in.cy<-NULL
in.r<-NULL
ref.cx<-NULL
ref.cy<-NULL
ref.r<-NULL
step2_mat<-list()
final.cx<-NULL
final.cy<-NULL
final.cr<-NULL
FM<-NULL
max_rotation_angle<-30
nseg=360
for ( k in 1:3){
print(paste("circle",k,"matching"))
in.cx[k]<-input_circles[k,1]
in.cy[k]<-input_circles[k,2]
in.r[k]<-input_circles[k,3]
circle_in<-data.frame(int_inside_center(data.frame(input), in.r[k], nseg, in.cx[k], in.cy[k]))
ref.cx[k]<-ref_circles[k,1]
ref.cy[k]<-ref_circles[k,2]
ref.r[k]<-ref_circles[k,3]
circle_ref<-data.frame(int_inside_center(data.frame(reference), ref.r[k], nseg, ref.cx[k],ref.cy[k]))
step2_mat<-boosted_clique(circle_in, circle_ref, ncross_in_bins = 30, xbins_in = 20,
ncross_in_bin_size = 1, ncross_ref_bins = NULL, xbins_ref = 30,
ncross_ref_bin_size = NULL, eps = 0.75, seed = 1, num_cores = parallel::detectCores()/2,
plot = TRUE, verbose = FALSE, cl = NULL)$clique_stats
final.cx[k]<-step2_mat[,6]
final.cy[k]<-step2_mat[,7]
final.r[k]<-step2_mat[,8]
MM<-rbind(MM,step2_mat)
}
Input_X<-input_circles[,1]
Input_Y<-input_circles[,2]
Comp<-c('1-2','1-3','2-3')
d_in_1<-sqrt((Input_X[1]-Input_X[2])^2+(Input_Y[1]-Input_Y[2])^2)
d_in_2<-sqrt((Input_X[1]-Input_X[3])^2+(Input_Y[1]-Input_Y[3])^2)
d_in_3<-sqrt((Input_X[2]-Input_X[3])^2+(Input_Y[2]-Input_Y[3])^2)
Euc_input_dist<-c(d_in_1,d_in_2,d_in_3)
Reference_X<-MM[,6]
Reference_Y<-MM[,7]
d_ref_1<-sqrt((Reference_X[1]-Reference_X[2])^2+(Reference_Y[1]-Reference_Y[2])^2)
d_ref_2<-sqrt((Reference_X[1]-Reference_X[3])^2+(Reference_Y[1]-Reference_Y[3])^2)
d_ref_3<-sqrt((Reference_X[2]-Reference_X[3])^2+(Reference_Y[2]-Reference_Y[3])^2)
Euc_ref_dist<-c(d_ref_1,d_ref_2,d_ref_3)
Reference_radius<-MM[,8]
#FM<-data.frame(Input_X,Input_Y,Reference_X,Reference_Y,Reference_radius,MM[,c(1:5)],Comp,Euc_input_dist,Euc_ref_dist, WW)
FM<-data.frame(Input_X,Input_Y,Reference_X,Reference_Y,Reference_radius,MM[,c(1:5)],Comp,Euc_input_dist,Euc_ref_dist)
P1<-ggplot(data.frame(input), aes(x=x, y=y))+ geom_point(data=data.frame(input), aes(x=x, y=y), color='black',size=0.1) +
geom_point(data=data.frame(int_inside_center(data.frame(input), in.r[1], nseg, in.cx[1], in.cy[1])),color="red",size=0.1)+
gg_circle(in.r[1], xc=in.cx[1], yc=in.cy[1], color="red") +
geom_point(data=data.frame(int_inside_center(data.frame(input), in.r[2], nseg, in.cx[2], in.cy[2])),color="yellow",size=0.1)+
gg_circle(in.r[2], xc=in.cx[2], yc=in.cy[2], color="yellow") +
geom_point(data=data.frame(int_inside_center(data.frame(input), in.r[3], nseg, in.cx[3], in.cy[3])),color="green",size=0.1)+
gg_circle(in.r[3], xc=in.cx[3], yc=in.cy[3], color="green")
P2<-ggplot(data.frame(reference), aes(x=x, y=y))+ geom_point(data=data.frame(reference), aes(x=x, y=y), color='black',size=0.1) +
geom_point(data=data.frame(int_inside_center(data.frame(reference), final.r[1], nseg, final.cx[1],final.cy[1])),color="red",size=0.1)+
gg_circle(final.r[1], xc=final.cx[1], yc=final.cy[1], color="red") +
geom_point(data=data.frame(int_inside_center(data.frame(reference), final.r[2], nseg, final.cx[2],final.cy[2])),color="yellow",size=0.1)+
gg_circle(final.r[2], xc=final.cx[2], yc=final.cy[2], color="yellow") +
geom_point(data=data.frame(int_inside_center(data.frame(reference), final.r[3], nseg, final.cx[3],final.cy[3])),color="green",size=0.1)+
gg_circle(final.r[3], xc=final.cx[3], yc=final.cy[3], color="green")
try(multiplot(P1, P2, cols=2))
return(FM)
}
#' @title Summarize MC-COMP macthing across three circles
#'
#' @description Get average of five similarity features and standard deviaitions of rotation angle estimations across three circle matching.
#'
#'
#' @name summarize_list_match
#' @param list_mat Matching result as list type
#' @param mat_class Class either KM or KNM
#' @param size the size of the shoe
#'
#' @export
#'
summarize_list_match<-function(list_mat, mat_class, size){
mat<-list_mat
class<-mat_class
rbind_mat<-do.call(base::rbind, mat)
number<-nrow(rbind_mat)/3 #10
rbind_mat$abs_Euc_diff<-abs(rbind_mat$Euc_input_dist-rbind_mat$Euc_ref_dist)
rbind_mat$Mat<- class
rbind_mat$comp<-factor(rep(1:number, each=3))
sum_mat<-rbind_mat[,c(6,7,8,9,10,14,15,16)]
#average score
ave_mat <- ddply(sum_mat, c("comp"), summarise,
class=mat_class,
size=size,
mean_clique_size=mean(clique_size),
sd_rotation_angle = sd(rotation_angle, na.rm = TRUE),
mean_reference_overlap=mean(reference_overlap),
mean_input_overlap=mean(input_overlap),
mean_med_dist_euc=mean(med_dist_euc),
mean_Euc_diff=mean(abs_Euc_diff)
)
return(ave_mat)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.