Fast corner detection using a look-up table in IDL: how to use a look-up table in IDL for improved speed

Implementing algorithms in IDL usually involve large arrays of data. One technique that can speed up some algorithms in IDL, is to make use of a precomputed look-up table. This essentially allows bypassing a computation step at the expense of increased memory use. This example implements the FAST corner detector algorithm. This is a simple algorithm compared to many other corner detector algorithms. Every pixel is compared with 16 other pixels forming a circle around the pixel in question. Each of the 16 surrounding pixels is evaluated as to whether it is considered similar or different from the reference pixel. Finally, the reference pixel is flagged as a corner if it has at least N consecutive surrounding pixels marked as significantly bigger or significantly smaller.

To test the speed, I used the following code as a reference.

  IDL> im=read_image(filepath('ohare.jpg', subdir=['examples','data']))
  % Loaded DLM: JPEG.
  IDL> im = reform(im[0,*,*])
  IDL> tic & x=fast_corner_detector(im) & toc
  % Compiled module: FAST_CORNER_DETECTOR.
  % Compiled module: ARRAY_INDICES.
  % Time elapsed: 42.885000 seconds.
  IDL> tic & z=fastcornerfinder(im, threshold=50, n_cont=12) & y=array_indices(z,where(z)) & toc
  % Compiled module: FASTCORNERFINDER.
  % Time elapsed: 2.6559999 seconds.
  IDL> array_equal(x,y)
  1
  IDL> help, x, y
  X               LONG      = Array[2, 4078]
  Y               LONG      = Array[2, 4078]

In this case the speed was improved by a factor of 16 compared to the reference code. This is meant to illustrate a technique for making IDL code run faster. I do not guarantee that the implementation is suitable for any specific purpose. The code is listed below:

 ;+
 ; Reference:
 ;  http://www.edwardrosten.com/work/fast.html
 ;-
 function FastCornerFinder, im, threshold=threshold, n_cont=n_cont
  compile_opt idl2, logical_predicate
  common fast_common, lookup, x_shift, y_shift

  ; number of consecutive matches to look for
  n = n_elements(n_cont) eq 0 ? 9b : byte(n_cont)
  if n lt 1 || n gt 16 then message, 'n_cont must be between 1 and 16'
  if n_elements(threshold) eq 0 then threshold = 50
  if n_elements(lookup) eq 0 then begin
    ; one-time common initialization
    ; for r = 3 the shifts are
    x_shift = [-1, 0, 1, 2, 3, 3, 3, 2, 1, 0,-1,-2,-3,-3,-3,-2]
    y_shift = [-3,-3,-3,-2,-1, 0, 1, 2, 3, 3, 3, 2, 1, 0,-1,-2]
    ; lookup table returns the maximum number of consecutive
    ; bits that are set, bitwise shift (ishft) is used.
    lookup = bytarr(65536)
    for i=0, 2^16-1 do begin
      x = i or ishft(i, 16)
      y = x
      for j=0, 15 do begin
        if y eq 0 then break
        y = x and ishft(y, 1)
      endfor
      lookup[i] = j
    endfor
  endif

  ; make an array where the bits represent whether each of the
  ; 16 positions around the circle is significantly different
  ; from the center. Test significantly smaller or bigger.
  bitsmaller = uintarr(size(im, /dimensions))
  bitbigger = uintarr(size(im, /dimensions))
  ; ensure signed pixels, so that subtraction can go negative
  fim = fix(im)
  ; loop over the 16 positions around the circle
  for i=0, 15 do begin
    bitsmaller or= ishft(1us, i) * ((fim - shift(fim, x_shift[i], y_shift[i])) gt fix(threshold))
    bitbigger or= ishft(1us, i) * ((fim - shift(fim, x_shift[i], y_shift[i])) lt -fix(threshold))
  endfor

  ; use the lookup array to convert to number of consecutive bits
  return, (lookup[bitsmaller] ge n) or (lookup[bitbigger] ge n)
end