Saddle and non-max suppression

saddle_nonmax.py

import numpy as np
from time import time

def getSaddlePoints(img_gray):
  img = np.array(img_gray, dtype=np.float32)

  gx, gy = np.gradient(img)  # row / col -> y / x
  gxx, gxy = np.gradient(gx)
  gyx, gyy = np.gradient(gy)

  gmag = np.sqrt(gx*gx + gy*gy)
  gphase = np.arctan2(gy,gx)
  gphase[gmag < 20] = 0

  S = gxx*gyy - gxy**2
  S = -S
  S[S<0] = 0
  S_suppress = nonmax_supress(S)

  thresh = 0
  count = np.count_nonzero(S_suppress)
  while count > 1000:
    S_suppress[S_suppress<thresh] = 0
    count = np.count_nonzero(S_suppress)
    thresh += 100

  return np.argwhere(S_suppress)

def nonmax_supress(arr, winsize=5):
  """Given 2D image, return suppressed nonmax pixels using window"""
  a = time()
  N = arr.shape[0] * arr.shape[1]
  arr_s = arr.copy()

  for i,j in np.argwhere(arr):
  # for i in range(arr.shape[0]):
  #   for j in range(arr.shape[1]):
      if arr[i,j] <= 0:
        continue
      win = arr[max(0, i-winsize) : min(i+winsize, arr.shape[0]-1),
                max(0, j-winsize) : min(j+winsize, arr.shape[1]-1)]
      
      if arr[i,j] < win.max():
        arr_s[i,j] = 0


  print("nonmax suppress of %d pixels in %g seconds" % (N, time()-a))
  return arr_s