Are there standard methods for computing the direction field of an image?

Question

I'd like to compute a sort of direction field on a 2D image, as (poorly) illustrated from this photoshop mockup. NOTE: This is NOT a vector field as you learn about in differential equations. Instead, this is something that draws along the lines that one would see if they computed level sets of the image.

Are there known methods of obtaining this type of direction field (red lines) of an image? It seems like it almost behaves like the normal to the gradient, but this isn't exactly it, either, since there are places where the gradient is zero and I'd like direction fields at these locations as well.

Answer 1

I was able to find a paper on how to do this for fingerprint processing that went into enough detail that their results were repeatable. It's unfortunately behind a paywall, but here it is for anyone interested and able to access the full text:

Systematic methods for the computation of the directional fields and singular points of fingerprints

---

EDIT: As requested, here is a quick and dirty summary (in Python) of how this is achieved in the above paper.

A naive approach would be to average the gradient in a small square neighborhood around the target pixel, much like the superimposed grid on the image in the question, and then compute the normal. However, if you simply average the gradient, it's possible that opposite gradients in the region will cancel each other (e.g. when computing the orientation along a ridge). Thus, it is common to compute with squared gradients, since gradients pointing in opposite directions would then be aligned. There is a clever formula for the squared gradient based on the original gradient. I won't give the derivation, but here is the formula:

Now, take the sum of squared gradients over the region (modulo some piece-wise defined compensations for the way angles work). Finally, through some arctangent magic, you'll get the orientation field.

If you run the following code on a smooth grayscale bitmap image with the grid-size chosen appropriately and then plot the orientation field O alongside your original image, you'll see how the orientation field more or less gives the angles I asked about in my original question.

from scipy import misc
import numpy as np
import math

# Import the grayscale image
bmp = misc.imread('path/filename.bmp')

# Compute the gradient - VERY important to convert to floats!
grad = np.gradient(bmp.astype(float))

# Set the block size (superimposed grid on the sample image in the question)
blockRadius=5

# Compute the orientation field. Result will be a matrix of angles in [0, \pi), one for each pixel in the original (grayscale) image.

O = np.zeros(bmp.shape)

for x in range(0,bmp.shape[0]):
    for y in range(0,bmp.shape[1]):
        numerator = 0.
        denominator = 0.
        for i in range(max(0,x-blockRadius),min(bmp.shape[0],x+blockRadius)):
            for j in range(max(0,y-blockRadius),min(bmp.shape[0],y+blockRadius)):
                numerator = numerator + 2.*grad[0][i,j]*grad[1][i,j]
                denominator = denominator + (math.pow(grad[0][i,j],2.) - math.pow(grad[1][i,j],2.))
        if denominator==0:
            O[x,y] = 0
        elif denominator > 0:
            O[x,y] = (1./2.)*math.atan(numerator/denominator)
        elif numerator >= 0:
            O[x,y] = (1./2.)*(math.atan(numerator/denominator)+math.pi)
        elif numerator < 0:
            O[x,y] = (1./2.)*(math.atan(numerator/denominator)-math.pi)
            
for x in range(0,bmp.shape[0]):
    for y in range(0,bmp.shape[1]):
        if O[x,y] <= 0:
            O[x,y] = O[x,y] + math.pi
        else:
            O[x,y] = O[x,y]

Cheers!