In srvanderplas/ShoeAlignR: Image Alignment

library(imager)
library(dplyr)
library(ImageAlignR)
if (!"ShoeSampleData" %in% installed.packages()) {
  devtools::install_github("srvanderplas/ShoeData")
}

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#",
  # prompt = "",
  fig.width = 5, fig.height = 9
)

The R code underlying the functions in this vignette were borrowed with permission from Vinces Gaitan's blog post: https://www.kaggle.com/vicensgaitan/image-registration-the-r-way/notebook

I have modified the code where necessary in order to make it more pipeline-friendly and function oriented.

imlinks <- system.file(package = "ShoeSampleData", "extdata/") %>%
  list.files(pattern = "036285L", full.names = T) %>%
  sort()

clean_shoe_img <- function(im) {
  suppressMessages({
    im_bbox <- im %>%
      imsplit(axis = "c") %>%
      (function(x) is.finite((x[[1]] + x[[2]]) / x[[3]])) %>%
      as.cimg() %>%
      (function(x) x == 1)

    crop.bbox(im, im_bbox) %>%
      grayscale() %>%
      map_halfimg(fun = autocrop) %>%
      crop.borders(nx = 5, ny = 5) %>%
      autocrop() %>%
      threshold() %>%
      shrink(3) %>%
      grow(3) %>%
      autocrop() %>%
      # img_rotate_refit() %>%
      # magrittr::extract2("img") %>%
      grayscale()
  })
}

img_a <- load.image(imlinks[1]) %>% clean_shoe_img()
img_b <- load.image(imlinks[2]) %>% clean_shoe_img()
plot(imlist(img_a, img_b))

We need to pad image a so that it is the same size as image b:

dim(img_a)
dim(img_b)
pad_size <- dim(img_b) - dim(img_a)
img_a <- pad(img_a, nPix = pad_size[1], axes = "x", pos = 1, val = max(img_a)) %>%
  pad(nPix = pad_size[2], axes = "y", pos = 1, val = max(img_a))

We can then overlay the two images to see how far apart they are:

plot(img_a + img_b)

Step 1: Keypoint Detection

The Harris corner detection algorithm has high values (light pixels) in areas where there are well-defined corners in the image.

hkp <- img_a %>%
  harris_corners(sigma = 6)

plot(hkp)

The detector seems to be working fairly well, outlining major features of the shoe print.

Thresholding the image from the previous step at 99% produces the best 1% of keypoints:

hkp_threshold <- hkp %>% threshold("99%")
hkp_threshold %>% plot()

The regions are labeled, and region centers are computed for each separately labeled region

keypoints <- hkp_threshold %>% label() %>% region_centers(bord = 30)

plot(img_a)
points(keypoints$mx, keypoints$my, col = "red")

Streamlined version

The first two steps can be streamlined using the harris_keypoints() function:

hkp <- harris_keypoints(img_a, sigma = 6)
plot(img_a)
points(hkp$centers$mx, hkp$centers$my, col = "red")

Step 2: Image Orientation

Calculating the dominant orientations for the whole image produces:

angles <- img_a %>% oriented_gradients(sigma = 3)
angles

Step 3: Feature Detection

For each angle, we pull features from a 40x40 area around the keypoint. These features will be used to identify points of similarity across the two images.

hkp_a <- harris_keypoints(img_a, sigma = 6)
hkp_b <- harris_keypoints(img_b, sigma = 6)

angles_a <- img_a %>% oriented_gradients(sigma = 3)
angles_b <- img_b %>% oriented_gradients(sigma = 3)

get_kpf <- function(angles, hkp, im) {
  kpa <- data_frame(angle = angles, v = list(hkp$centers)) %>%
    tidyr::unnest(v) %>%
    dplyr::rename(theta = angle, x = mx, y = my) %>%
    mutate(idx = 1:n()) %>%
    rowwise() %>%
    tidyr::nest(-theta, -idx, .key = "v") %>%
    select(-idx)
  purrr::pmap(list(theta = kpa$theta, v = kpa$v), descriptor_orientation, im = grayscale(im)) %>%
    do.call("rbind", .)
}

kpf_a <- get_kpf(angles_a, hkp_a, img_a)
kpf_b <- get_kpf(angles_b, hkp_b, img_b)

Step 4: Match points

Match points are calculated using the K nearest neighbors algorithm, combined with some thresholding by distance.

match_points <- knn_points(kpf_a, kpf_b, hkp_a$centers, hkp_b$centers, show_plot = T)

Step 5: RANSAC

RANSAC is then used to find points that have similar homography.

ransac_points <- ransac(match_points$points_a, match_points$points_b)

par(mfrow = c(1, 2))

plot(img_a)
hkp_a$centers %$% points(mx, my, col = "orange")
points(match_points$points_a[ransac_points$inliers, ], col = "purple", pch = 16)
plot(img_b)
hkp_b$centers %$% points(mx, my, col = "orange")
points(match_points$points_b[ransac_points$inliers, ], col = "purple", pch = 16)

Step 6: Image Warping

The homography can be used to warp one image onto the other:

map_fcn <- map_affine_gen(ransac_points$homography)

img_a_warp <- imwarp(img_a, map_fcn, direction = "backward", boundary = "neumann")

plot(img_a_warp)

We can then overlay the two images:

blank_channel <- as.cimg(img_b > 0 & img_a_warp > 0)
overlaid_images <- imappend(imlist(img_a_warp, blank_channel, img_b), axis = "c")
plot(overlaid_images)