Analyze bullet striations using nonparametric methods

Comparing lands from two bullets

Striae comparisons between bullets are based on land-to-land comparisons.

  1. Load libraries
library(dplyr)
library(bulletxtrctr)
library(x3ptools)
library(randomForest)
library(ggplot2)
library(readr)
library(nbtrd) # devtools::install_github("csafe-isu/nbtrd")
  1. 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")
## 6 files found. Reading ...
b2 <- read_bullet("README_files/data/Bullet2")
## 6 files found. Reading ...

Combine the results into a single data frame:

We expect data to be recorded at the micron level. The scans posted give measurements in meters:

## [1] 1.5625e-06
## [1] 1.5625e-06
summary(as.vector(bullets$x3p[[1]]$surface.matrix))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##       0       0       0       0       0       0   24829

Change measurements to microns:

## [1] 1.5625
## [1] 1.5625
summary(as.vector(bullets$x3p[[1]]$surface.matrix))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.513 117.626 166.723 155.933 199.429 216.341   24829

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")

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.

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")

  1. Get the ideal cross sections

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")

  1. Get the groove locations

Visualize that the grooves are identified correctly (at least enough to not distort the final result):

  1. Extract signatures
signatures <- bullets %>% select(source, sigs) %>% tidyr::unnest()
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()

  1. Align signatures and extract features

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")
## Warning: Removed 37 rows containing missing values (geom_path).

Some features are based on aligned signatures:

Other features need an evaluation of striation marks between two aligned signatures:

comparisons <- comparisons %>% mutate(
  cms_per_mm = purrr::map2(striae, aligned, .f = function(s, a) {
    extract_feature_cms_per_mm(s$lines, a$lands, resolution=1.5625)
  }),
  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)
)
  1. Extract Features
comparisons <- comparisons %>% mutate(
  features = purrr::map2(.x = aligned, .y = striae, .f = extract_features_all, resolution = 1.5625)
)

comparisons <- comparisons %>% mutate(
  legacy_features = purrr::map(striae, extract_features_all_legacy, resolution = 1.5625)
)

comparisons <- comparisons %>% tidyr::unnest(legacy_features) 
# scale features before using them in the random forest, legacy features can be used out of the box


# quick visualization:
comparisons %>% 
  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)

  1. Get Score predictions for each land to land comparison
comparisons$rfscore <- predict(bulletr::rtrees, newdata = comparisons, type = "prob")[,2]

comparisons %>% 
  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)

  1. Determine bullet-to-bullet scores
## # A tibble: 4 x 3
##   bulletA bulletB bullet_score
##   <chr>   <chr>          <dbl>
## 1 1       1              0.982
## 2 2       1              0.681
## 3 1       2              0.681
## 4 2       2              0.989
  1. Use bullet-to-bullet scores to predict land to land scores

An interactive interface for doing comparisons is available https://oaiti.org/apps/bulletmatcher/