Radial interpolation for corneal reflex removal with Starburst

Hi everybody,

I'm trying to break down the Starburst algorithm to do some eye tracking. The code is old and needs updating, so I'm trying to rewrite the outdated parts of code to get it to work again, but it is hard as it is very poorly commented (if ever).

I got stuck onto the function that is supposed to remove the corneal reflection previuosly identified on the eye image. I have the position onto the image and the approximate radius of it.

The concept of radial interpolation for removal is very easy: the central point of the region is given the value of the average gray levels of the points onto the perimeter of the corneal reflex. Then the points from the border to the center are linarly interpolated to assign their gray values.

The problem I have is with the two last lines of code, that give a bunch of syntax errors. I can fix the syntax errors so that the function actually produces output, but I'd like to know if what I obtain is reasonable. The code writing is too compact and not enough commented for me to understand if I'm actually getting a reasonable result. Can anybody help me out with this?

I attach here the original code:

function [I] = remove_corneal_reflection(I, crx, cry, crr, angle_delta)

% Input:
% I = input image
% angle_delta = angle step size
% [crx cry] = corneal reflection center
% crr = corneal reflection radius
% Output:
% I = output image with the corneal reflection removed

if crx==0 | cry==0 | crr<=0
   return;
end

[height width] = size(I);

if crx-crr < 1 || crx+crr > width || cry-crr < 1 || cry+crr > height
    fprintf(1, 'Error! Corneal reflection is too near the image border\n');
    return;
end

theta=[0:pi/360:2*pi];
tmat = repmat(theta',[1 crr]);
rmat = repmat([1:crr],[length(theta) 1]);
[xmat,ymat]=pol2cart(tmat,rmat);
xv=reshape(xmat,[1 prod(size(xmat))]);
yv=reshape(ymat,[1 prod(size(ymat))]);

avgmat=ones(size(rmat))*mean(I(round(ymat(:,end)+cry)+(round(xmat(:,end)+crx)-1)*height));
permat=repmat(I(round(ymat(:,end)+cry)+(round(xmat(:,end)+crx)-1)*height),[1 crr]);
wmat=repmat((1:crr)/crr,[length(theta) 1]);
imat=avgmat.*(1-wmat) + permat.*(wmat); 
I(round(yv+cry)+(round(xv+crx-1).*height))=imat;

and this is how I fixed the last lines:

avgmat=ones(size(rmat))*mean(I(round(ymat(:,end)+cry)+(round(xmat(:,end)+crx)-1)*height));
permat=repmat(I(round(ymat(:,end)+cry)+(round(xmat(:,end)+crx)-1)*height),[1 crr]);
permat = double(permat);
wmat=repmat((1:crr)/crr,[length(theta) 1]);
imat=-avgmat.*(wmat-1) + permat.*(wmat); 
I(round(yv+cry)+(round(xv+crx-1).*height))=imat;

[Walter Roberson]

? Those lines do not have a syntax error, and I cannot see how moving the negative would make a difference?

The double() you inserted could potentially be important, depending on the class of the data.

[Beatrice Pazzucconi]

The double I inserted is essential in that avgmat is double and permat is initially uint8. Not putting the double results in error.

Those lines I posted originally do give a syntax error: you cannot do 1 - matrix, results in error, while doing matrix - 1 takes away 1 from all matrix entries.

The point is: was 1 - matrix intented to be like matrix - 1 with negative sign or like I - matrix? Because I don't see how I could build an identity out of a non-square matrix (because wmat in general is not squared), but since we're using regression, it would make sense to me to do so.

[Image Analyst]

Can you attach a sample image that it's supposed to work with along with values for the other input parameters?

[Walter Roberson]

(1 - matrix) should not be an error for any matrix of any of the numeric classes, or the logical class, or the character class, or class sym or symfun. It would be an error for class struct or cell, and for many object oriented classes, but for those classes, doing (matrix - 1) would not be any better.

[Beatrice Pazzucconi]

This is the original image

This is the one I get with the corneal reflection pointed out (the red cross):

The function input parameters are: I = the first image, after being turned into grayscale and smoothed with gaussian filetr wit sigma = 2:

I = rgb2gray(I); % first convert image to grayscale
f = fspecial('gaussian', [ceil(2.5*sigma) ceil(2.5*sigma)], sigma);
I = imfilter(I, f, 'symmetric');

and the points identifying the corneal reflection (from another function), crx, cry, crr correspondind to corneal reflection x and y coordinates and approximate radius (the x and y coordinates are plotted in the second figure).

The image is 720x1280 pixels, and

crx =

679.4065

cry =

310.4634

crr =

    26

The input angle_delta is not used in the original function, so I did not use it in my version. The part I did not post of the original code is the same as the original code If I run the original code:

Error using  .* 
Integers can only be combined with integers of the same class, or scalar 
Error in remove_corneal_reflection (line 58)
imat=avgmat.*(1-wmat) + permat.*(wmat);

This is what I get calling the function directly. If I run the algorithm line by line, it is clear it oringinates from the fact that permat is uint8 and not double. If I correct it, as in my algorithm, it gave the same error, so I though it could be the 1-wmat.

UPDATE - re-running the code today the error is not given anymore. I attach here the resuting image:

So I guess it was maybe a problem with residual variables.