R/step1_selectBF.R

In xhtai/cartridges: Algorithm to Compare Cartridge Case Images

#' Read cartridge case image and find primer region
#'
#' Images must be of the standard format as provided in the NIST Ballistics and
#' Research Toolmark Database. The dimensions are 1944 x 2592 pixels (1944 is
#' the height), and images are on a 255-grayscale. Here we use ring light
#' images. This function finds the primer region, which is the circular region
#' containing the breechface marks and firing pin impression. By default images
#' are cropped to 1919 x 1919 pixels. For more information on the steps involved
#' in finding this region, refer to the GitHub page.
#'
#' This function uses the \code{EBImage} package, which stores the image as a
#' 2592 x 1944 matrix. The primer region that is returned is similarly an
#' \code{EBImage} image, where the primer region is stored in a matrix with
#' binary values. If using this matrix directly on the image matrix in the
#' original orientation (1944 x 2592), note that the should first be transposed.
#'
#'
#' @param fileName location of image file. This function has been tested on
#'   \code{png} files.
#'
#' @return A 1919 x 1919 binary image in \code{EBImage} format. Values 1
#'   indicate the primer region.
#' @examples
#' \dontrun{
#' primerExample <- findPrimer(system.file("extdata", "NBIDE R BF 118.png", package = "cartridges"))
#' }
#'
#' @export

findPrimer <- function(fileName) {
    ebimage <- EBImage::readImage(fileName)
    if (length(dim(ebimage)) == 3) {
        ebimage <- ebimage[, , 1]
    }
    I <- ebimage[337:2255, 13:1931]

    sigma <- 25
    blurred <- EBImage::gblur(I, sigma, radius = 2 * ceiling(2 * sigma) + 1)

    equalized <- makeBinary(blurred)

    # imfill
    tmp <- EBImage::fillHull(equalized)

    # bwlabel -- choose the center region
    L = EBImage::bwlabel(tmp)
    ij <- matrix(c(959, 600, 1200, 959, 959, 780, 780, 1080, 1080,
                   960, 960, 960, 600, 1200, 780, 1080, 780, 1080), ncol = 2)

    Mode <- function(x) {
        ux <- unique(x)
        ux[which.max(tabulate(match(x, ux)))]
    }
    tmp[L != Mode(L[cbind(ij)])] <- 0

    ##### dilate, fill and erode
    se <- EBImage::makeBrush(size = 139, shape = "disc") # should be odd
    tmp <- EBImage::dilate(tmp, se)
    filled <- EBImage::fillHull(tmp)
    se <- EBImage::makeBrush(size = 269, shape = "disc") # should be odd
    outer <- EBImage::erode(filled,se)
    return(outer)
}


#' Find and remove firing pin region
#'
#' After finding the primer region using \code{findPrimer()}, we now remove the
#' firing pin impression, returning a 1919 x 1919 matrix, where areas that are
#' not part of the breechface impression are set to 0. Pixel values range from 0
#' to 255. For more information on the steps involved in finding the firing pin
#' region, refer to the GitHub page.
#'
#' @param fileName location of image file. This function has been tested on
#'   \code{png} files.
#' @param primer binary image in \code{EBImage} format, as returned by
#'   \code{findPrimer}
#'
#' @return A 1919 x 1919 matrix, in the original orientation, i.e. a cropped
#'   version of the original 1944 x 2592 image.
#' @examples
#' \dontrun{
#' FPexample <- findFP(system.file("extdata", "NBIDE R BF 118.png", package = "cartridges"), primer = primerExample)
#' }
#'
#' @export
#' @import methods

findFP <- function(fileName, primer) {
    if (!requireNamespace("imager", quietly = TRUE) || !requireNamespace("purrr", quietly = TRUE) || !requireNamespace("dplyr", quietly = TRUE)) {
        stop("This function uses a Canny edge detector, and this implementation requires the packages imager, purrr, and dplyr. Please install them. For an example, see help(FPexample).", call. = FALSE)
    }
    ebimage <- EBImage::readImage(fileName)
    if (length(dim(ebimage)) == 3) {
        ebimage <- ebimage[, , 1]
    }
    J <- makeBinary(ebimage)
    J <- J[337:2255, 13:1931]

    # canny
    sigma <- 32
    blurred <- EBImage::gblur(J, sigma)

    tmp <- blurred@.Data
    imagerImage <- imager::as.cimg(tmp)
    out <- cannyEdges(imagerImage*255, q1 = .83*.4, q2 = .83)

    out[primer == 0] <- 0

    if (sum(out) < 2000) {
        thres <- .83
        while (sum(out) < 2000) {
            thres <- thres - .02
            out <- cannyEdges(imagerImage*255, q1 = thres*.4, q2 = thres)
            out[primer == 0] <- 0
        }
    }

    outtmp <- as.matrix(out)
    outebimage <- EBImage::Image(outtmp, colormode = 'Grayscale') # back to ebimage

    se <- EBImage::makeBrush(size = 399, shape = "disc") # should be odd
    tmp <- EBImage::dilate(outebimage, se)
    tmp <- EBImage::fillHull(tmp)
    se <- EBImage::makeBrush(size = 379, shape = "disc") # should be odd
    FP <- EBImage::erode(tmp,se)

    while (sum(FP == 1) < 675000) { #3000000 -- fix this -- percentage of BF to remove
        se <- EBImage::makeBrush(9, 'disc')
        FP <- EBImage::dilate(FP, se)
    }

    # second pass for FP
    newJ <- J
    newJ[primer == 0] <- 1 # 1 instead of 255
    newJ[FP == 1] <- 0

    sigma <- 10 # changed from 2 -- blur more
    blurred <- EBImage::gblur(newJ, sigma)
    tmp <- blurred@.Data
    imagerImage <- imager::as.cimg(tmp)
    tmp <- cannyEdges(imagerImage*255, q1 = .2, q2 = .9)

    outtmp <- as.matrix(tmp) + FP@.Data # combine with earlier FP
    outtmp[outtmp > 1] <- 1
    outebimage <- EBImage::Image(outtmp, colormode = 'Grayscale') # back to ebimage

    se <- EBImage::makeBrush(139, 'disc')
    tmp <- EBImage::dilate(outebimage,se)
    tmp <- EBImage::fillHull(tmp)
    se <- EBImage::makeBrush(139, 'disc')
    FP <- EBImage::erode(tmp, se)

    newI <- ebimage[337:2255, 13:1931]
    newI[FP == 1] <- 0 #FP==1 (white) -- remove
    newI[primer == 0] <- 0 #BW3==0 (black) -- remove (outer regions)

    # convert to matrix, 0-255
    newI <- t(newI@.Data * 255)
    return(newI)
}

# ACKNOWLEDGEMENTS: the rest of the code in this file was adapted from Simon
# Barthelme, "A loop-free canny edge detector"
# (http://dahtah.github.io/imager/canny.html)
#' @importFrom magrittr "%>%"

fillInit <- function(strong) {
    lab <- imager::label(strong, TRUE)*strong
    as.data.frame(lab) %>% dplyr::filter(value > 0) %>% dplyr::group_by(value) %>% dplyr::summarize(x = x[1], y = y[1])
}

#Starts a fill at each successive location, and accumulates the results
rescueFill <- function(strong, weak) {
    v <- strong
    v[weak == 1] <- .9
    loc <- fillInit(strong)
    #Transform the data.frame into a list of locations
    tmp <- dplyr::select(loc, -value)
    loc <- purrr::transpose(tmp)
    #Fold
    out <- purrr::reduce(loc, function(v,l) imager::bucketfill(v, l$x, l$y, color = 1, sigma = .1, high = TRUE), .init = v)
    imager::as.cimg(out == 1)
}

nonmax <- function(gr) {
    mag <- with(gr, sqrt(x^2 + y^2))
    ang <- with(gr, atan2(y, x))
    grs <- list(x = gr$x/mag, y = gr$y/mag)
    X <- imager::Xc(gr$x)
    Y <- imager::Yc(gr$y)
    val.bwd <- imager::interp(mag, data.frame(x = as.vector(X - grs$x),
                                     y = as.vector(Y - grs$y)))
    val.fwd <- imager::interp(mag,data.frame(x = as.vector(X + grs$x),
                                     y = as.vector(Y + grs$y)))

    throw <- (mag < val.bwd) | (mag < val.fwd)
    mag[throw] <- 0
    mag
}


cannyEdges <- function(im, q1, q2){
    mag <- imager::imgradient(im, "xy") %>% nonmax
    t1 <- q1*max(mag)
    t2 <- q2*max(mag)
    strong <- imager::as.cimg(mag > t2)
    weak <- imager::as.cimg(mag >= t1 & mag <= t2)

    rescueFill(strong, weak)

}

makeBinary <- function(image) {
    med <- median(image)
    ret1 <- ifelse(image < med, 0, 1)
    diff1 <- abs(sum(ret1 == 1) - sum(ret1 == 0))
    ret2 <- ifelse(image <= med, 0, 1)
    diff2 <- abs(sum(ret2 == 1) - sum(ret2 == 0))

    if (diff1 <= diff2) return(ret1) else return(ret2)
}