knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "man/figures/README-", out.width = "100%" )
Analyze bullet striations using nonparametric methods
You can install the released version of bulletxtrctr from CRAN with:
install.packages("bulletxtrctr")
And the development version from GitHub with:
# install.packages("devtools") devtools::install_github("heike/bulletxtrctr")
Striae comparisons between bullets are based on land-to-land comparisons.
library(dplyr) library(bulletxtrctr) library(x3ptools) library(randomForest) library(ggplot2) library(readr) library(nbtrd) # devtools::install_github("csafe-isu/nbtrd")
bulletxtrctr
only works on x3p files. See package x3ptools
at https://heike.github.io/x3ptools/ for ways to convert different file formats into x3p standard files.
The NIST Research Ballistics Toolmarks data base (NRBTD)[https://tsapps.nist.gov/NRBTD/Studies/Search] provides access to scans of bullets and cartridge cases from various case studies. In this tutorial, we'll work with two bullets from a single barrel of the Hamby 252 data set. Links to the 12 scans of bullet lands in x3p format are provided in the hamby252demo
object.
These commands will read in the bullets directly from the NRBTD repository, without downloading the files into your working directory:
b1 <- read_bullet(urllist = hamby252demo[[1]]) b2 <- read_bullet(urllist = hamby252demo[[2]])
Instead, we could also download the files into a folder named "data" in our working directory. This is shown in the sequence of commands below:
if (!dir.exists("README_files/data")) { dir.create("README_files/data") } if (!file.exists("README_files/data/Bullet1/Hamby252_Barrel1_Bullet1_Land1.x3p")) { NRBTDsample_download("README_files/data") } b1 <- read_bullet("README_files/data/Bullet1") b2 <- read_bullet("README_files/data/Bullet2")
Combine the results into a single data frame:
b1$bullet <- 1 b2$bullet <- 2 b1$land <- 1:6 b2$land <- 1:6 bullets <- rbind(b1, b2)
We expect data to be recorded at the micron level. The scans posted give measurements in meters:
bullets$x3p[[1]] %>% x3p_get_scale() summary(as.vector(bullets$x3p[[1]]$surface.matrix))
Change measurements to microns:
bullets <- bullets %>% mutate( x3p = x3p %>% purrr::map(.f = x3p_m_to_mum) )
bullets$x3p[[1]] %>% x3p_get_scale() summary(as.vector(bullets$x3p[[1]]$surface.matrix))
Save resolution of the scans in a variable (we will be using this piece of information later:
resolution <- bullets$x3p[[1]] %>% x3p_get_scale() resolution
We are working under the assumption that the scans are aligned such that the base of the bullet are at the bottom (y = 0) of the image, and the land engraved area is displayed left to right from groove to groove, i.e. we are assuming that (0,0) is in the bottom left corner of the image. In scans where no adjustment was made for the barrel's twist (not recommended) the twist will be visible in the image.
image_x3p(bullets$x3p[[1]], file = "man/figures/temp-before.png")
knitr::include_graphics("man/figures/temp-before.png")
The raw scan needs to be flipped such that the heel is along the bottom of the image rather than along the left hand side.
# turn the scans such that (0,0) is bottom left bullets <- bullets %>% mutate( x3p = x3p %>% purrr::map(.f = function(x) x %>% rotate_x3p(angle = -90) %>% y_flip_x3p()) )
Scan after the transformation: a clear right twist is visible in the right slant of striae and grooves:
image_x3p(bullets$x3p[[1]], file = "man/figures/temp-after.png")
knitr::include_graphics("man/figures/temp-after.png")
bullets <- bullets %>% mutate( crosscut = x3p %>% purrr::map_dbl(.f = x3p_crosscut_optimize) ) # now extract the crosscuts bullets <- bullets %>% mutate( ccdata = purrr::map2(.x = x3p, .y = crosscut, .f = x3p_crosscut) )
Visualize the cross cuts:
crosscuts <- bullets %>% tidyr::unnest(ccdata) crosscuts %>% ggplot(aes(x = x, y = value)) + geom_line() + facet_grid(bullet~land, labeller="label_both") + theme_bw()
Note the rather strange cross cut for land 6 in bullet 1. We can look at the scan - and find quite pronounced tank rash. However, the extraction of the land is at a height of 375, which is not as much affected by the rash as the base of the bullet or the top of the scanning area.
Scan of land 6 on bullet 1. The land is affected by quite pronounced tank rash:
filter(bullets, land==6, bullet==1)$x3p[[1]] %>% x3p_add_hline(yintercept = 375, size = 10, color = "#e6bf98") %>% image_x3p(file="man/figures/bullet1-land6.png")
knitr::include_graphics("man/figures/bullet1-land6.png")
bullets <- bullets %>% mutate( grooves = ccdata %>% purrr::map(.f = cc_locate_grooves, method = "middle", adjust = 30, return_plot = TRUE) )
Visualize that the grooves are identified correctly (at least enough to not distort the final result):
gridExtra::grid.arrange( bullets$grooves[[1]]$plot, bullets$grooves[[2]]$plot, bullets$grooves[[3]]$plot, bullets$grooves[[4]]$plot, bullets$grooves[[5]]$plot, bullets$grooves[[6]]$plot, bullets$grooves[[7]]$plot, bullets$grooves[[8]]$plot, bullets$grooves[[9]]$plot, bullets$grooves[[10]]$plot, bullets$grooves[[11]]$plot, bullets$grooves[[12]]$plot, ncol = 6 )
bullets <- bullets %>% mutate( sigs = purrr::map2( .x = ccdata, .y = grooves, .f = function(x, y) { cc_get_signature( ccdata = x, grooves = y, span1 = 0.75, span2 = 0.03) }) )
signatures <- bullets %>% select(source, sigs) %>% tidyr::unnest(sigs) signatures %>% filter(!is.na(sig),!is.na(raw_sig)) %>% ggplot(aes(x = x)) + geom_line(aes(y = raw_sig), colour = "grey70") + geom_line(aes(y = sig), colour = "grey30") + facet_wrap(~source, ncol = 6) + ylim(c(-5,5)) + theme_bw()
bullets$bulletland <- paste0(bullets$bullet,"-", bullets$land) lands <- unique(bullets$bulletland) comparisons <- data.frame( expand.grid(land1 = lands, land2 = lands), stringsAsFactors = FALSE) comparisons <- comparisons %>% mutate( aligned = purrr::map2(.x = land1, .y = land2, .f = function(xx, yy) { land1 <- bullets$sigs[bullets$bulletland == xx][[1]] land2 <- bullets$sigs[bullets$bulletland == yy][[1]] land1$bullet <- "first-land" land2$bullet <- "second-land" sig_align(land1$sig, land2$sig) }) )
Aligned signatures of two matching lands:
subset(comparisons, land1=="2-4" & land2 =="1-2")$aligned[[1]]$lands %>% tidyr::gather(sigs, value, sig1, sig2) %>% ggplot(aes(x = x, y = value, colour = sigs)) + geom_line() + theme_bw() + scale_color_brewer(palette = "Dark2")
Some features are based on aligned signatures:
comparisons <- comparisons %>% mutate( ccf0 = aligned %>% purrr::map_dbl(.f = function(x) extract_feature_ccf(x$lands)), lag0 = aligned %>% purrr::map_dbl(.f = function(x) extract_feature_lag(x$lands)), D0 = aligned %>% purrr::map_dbl(.f = function(x) extract_feature_D(x$lands)), length0 = aligned %>% purrr::map_dbl(.f = function(x) extract_feature_length(x$lands)), overlap0 = aligned %>% purrr::map_dbl(.f = function(x) extract_feature_overlap(x$lands)) )
Other features need an evaluation of striation marks between two aligned signatures:
comparisons <- comparisons %>% mutate( striae = aligned %>% purrr::map(.f = sig_cms_max, span = 75) )
comparisons <- comparisons %>% mutate( cms_per_mm = purrr::map2(striae, aligned, .f = function(s, a) { extract_feature_cms_per_mm(s$lines, a$lands, resolution=resolution) }), matches0 = striae %>% purrr::map_dbl(.f = function(s) { bulletxtrctr:::extract_helper_feature_n_striae(s$lines, type = "peak", match = TRUE) }), mismatches0 = striae %>% purrr::map_dbl(.f = function(s) { bulletxtrctr:::extract_helper_feature_n_striae(s$lines, type = "peak", match = FALSE) }) )
comparisons <- comparisons %>% mutate( bulletA = gsub("([1-2])-([1-6])","\\1",land1), bulletB = gsub("([1-2])-([1-6])","\\1",land2), landA = gsub("([1-2])-([1-6])","\\2",land1), landB = gsub("([1-2])-([1-6])","\\2",land2) )
comparisons <- comparisons %>% mutate( features = purrr::map2(.x = aligned, .y = striae, .f = extract_features_all, resolution = resolution) ) comparisons <- comparisons %>% mutate( legacy_features = purrr::map(striae, extract_features_all_legacy, resolution = resolution) ) legacy <- comparisons %>% tidyr::unnest(legacy_features) # scale features before using them in the random forest, legacy features can be used out of the box features <- comparisons %>% select(land1, land2, ccf0, bulletA, bulletB, landA, landB, features) %>% tidyr::unnest(features) # quick visualization: features %>% ggplot(aes(x = landA, y = landB, fill = ccf)) + geom_tile() + scale_fill_gradient2(low = "grey80", high = "darkorange", midpoint = 0.5) + facet_grid(bulletB~bulletA, labeller = "label_both") + xlab("Land A") + ylab("Land B") + theme(aspect.ratio = 1)
require(randomForest) # use the scaled features, but the unscaled names ... features <- features %>% mutate( cms = cms_per_mm, matches = matches_per_mm, mismatches = mismatches_per_mm, non_cms = non_cms_per_mm ) features$rfscore <- predict(rtrees, newdata = features, type = "prob")[,2] features %>% ggplot(aes(x = landA, y = landB, fill = rfscore)) + geom_tile() + scale_fill_gradient2(low = "grey80", high = "darkorange", midpoint = .5) + facet_grid(bulletB~bulletA, labeller = "label_both") + xlab("Land A") + ylab("Land B") + theme(aspect.ratio = 1)
bullet_scores <- features %>% group_by(bulletA, bulletB) %>% tidyr::nest() bullet_scores <- bullet_scores %>% mutate( bullet_score = data %>% purrr::map_dbl( .f = function(d) max(compute_average_scores(land1 = d$landA, land2 = d$landB, d$rfscore))) ) bullet_scores %>% select(-data)
bullet_scores <- bullet_scores %>% mutate( data = data %>% purrr::map( .f = function(d) { d$samesource = bullet_to_land_predict( land1 = d$landA, land2 = d$landB, d$rfscore, difference=0.1) d }) ) features <- bullet_scores %>% tidyr::unnest(data) features %>% ggplot(aes(x = landA, y = landB, fill = rfscore, colour=samesource)) + geom_tile() + scale_fill_gradient2(low = "grey80", high = "darkorange", midpoint = .5) + scale_colour_manual(values = c("grey80", "darkorange")) + geom_tile(size = 1, data = features %>% filter(samesource)) + facet_grid(bulletB~bulletA, labeller = "label_both") + xlab("Land A") + ylab("Land B") + theme(aspect.ratio = 1)
An interactive interface for doing comparisons is available https://oaiti.org/apps/bulletmatcher/
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.