Skip to content

Instantly share code, notes, and snippets.

@dimitrs
Created June 15, 2014 13:40
Show Gist options
  • Star 3 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save dimitrs/d667ccb433d5f2baa8df to your computer and use it in GitHub Desktop.
Save dimitrs/d667ccb433d5f2baa8df to your computer and use it in GitHub Desktop.
A Python implementation of the paper "Accurate Eye Center Location through Invariant Isocentric Patterns".
import cv2
import numpy as np
import matplotlib.pyplot as plt
from skimage import io
img = io.imread('eye.png')
# Left eye
img = np.around(img[0:img.shape[0], 0:img.shape[1]/2])
gray_image = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
gray_image = np.float32(gray_image)
# Note: at higher scales the isophotes curvature might degenerate with the inherent
# effect of losing important structures in the image.
gray_image = cv2.GaussianBlur(gray_image,(9,9),0)
# Calculate Isophotes curvature k where Lx and Ly are the first order derivatives
# of the luminance function L(x,y) in the x and y dimension, respectively.
# The sign of the isophote depends on the intensity of the outer side of the curve,
# so its vector direction will always point towards the highest change in the luminance.
Ly,Lx = np.gradient(gray_image)
Lxy, Lxx = np.gradient(Lx)
Lyy, Lyx = np.gradient(Ly)
Lvv = Ly**2 * Lxx - 2*Lx * Lxy * Ly + Lx**2 * Lyy
Lw = Lx**2 + Ly**2
k = - Lvv / (Lw**1.5)
# The displacement vectors Dx and Dy, point to the estimated position of the centers.
Dx = -Lx * (Lw / Lvv)
Dy = -Ly * (Lw / Lvv)
magnitude_displacement = np.sqrt(Dx**2 + Dy**2)
# The curvedness indicates how curved a shape is.
curvedness = np.absolute(np.sqrt(Lxx**2 + 2 * Lxy**2 + Lyy**2))
# Sometimes the isophote curvature could assume extremely small or big values.
# This indicates that we are dealing with a “straight line” or a “single dot”
# isophote. Since the estimated radius to the isophote center would be too high
# to fall into the centermap or too little to move away from the originating pixel,
# the calculation of the displacement vectors in these extreme cases can simply
# be avoided.
minrad = 2
maxrad = 20
# A negative sign indicates a change in the direction of the gradient
# (i.e. from brighter to darker areas). Therefore, it is possible to
# discriminate between dark and bright centers by analyzing the sign of the
# curvature
center_map = np.zeros(gray_image.shape, gray_image.dtype)
(height, width)=center_map.shape
for y in range(height):
for x in range(width):
if Dx[y][x] == 0 and Dy[y][x] == 0:
continue
if (y + Dx[y][x])>0 and (x + Dy[y][x])>0:
if (x + Dx[y][x]) < center_map.shape[1] and (y + Dy[y][x]) < center_map.shape[0] and k[y][x]<0:
if magnitude_displacement[y][x] >= minrad and magnitude_displacement[y][x] <= maxrad:
center_map[int(y+Dy[y][x])][int(x+Dx[y][x])] += curvedness[y][x]
# Since every vector gives a rough estimate of the center, the accumulator
# is convolved with a Gaussian kernel so that each cluster of
# votes will form a single center estimate.
center_map = cv2.GaussianBlur(center_map,(5,5),0)
# Maximum isocenter (MIC) in the centermap will be used as the most probable
# estimate for the soughtafter location of the center of the eye.
mic = np.unravel_index(np.argmax(center_map), center_map.shape)
# Draw it
fig = plt.figure()
plt.subplot(121),plt.imshow(gray_image, cmap = 'gray'),plt.title('Input')
fig.gca().add_artist(plt.Circle((mic[1],mic[0]),2,color='r'))
plt.subplot(122),plt.imshow(center_map, cmap = 'gray'),plt.title('Output')
plt.show()
from pylab import *
savefig('./eye.jpg')
@dimitrs
Copy link
Author

dimitrs commented Jun 15, 2014

@willard-yuan
Copy link

Is it wrong in line 59?

if (y + Dx[y][x])>0 and (x + Dy[y][x])>0:

I think it should be:

if (x + Dx[y][x])>0 and (y + Dy[y][x])>0:

@asanchez19
Copy link

Hello, I was wondering what's the python version ? I'm getting an error on line 9. And my img is not null.
Thank you in advance.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment