%--------------------------------------------------------
%Interpol: a function for interpolation of eye blinks
%in eye pupil data
%Usage:
%interpol(rawdiamet,samplerate)
%
%(c) EJh 2007-2010
%--------------------------------------------------------


function diamet = interpoll(rawdiamet,eyeSamplerate)

%Settings
%windowstart = -1;   %cut out 1 second pre-stim
%windowend = 6;      %until 6 sec post-stim
diffthresh = 20;     %threshold for eyeblink detection
blinkwindowlength = 10;  %number of datapoints discarded after blink
mincontiggood = 10;  %the minimum number of contiguous correct measurement in diameter data
interpolationvals = 2;  %number of datapoints used for interpolation (at start and end of interpolation)


%Clean up the pupil diameter data
%diamet = bpf(rawdiamet,eyeSamplerate,(1/60),-1);
diamet=rawdiamet;

%find all onsets of eyeblinks
errordiamet = vertcat(abs(diff(diamet))>diffthresh,0);

%Now discard all data within blinks, looping over too high diffs
blinkonsets = find(errordiamet);
for i=1:length(blinkonsets)
    errordiamet(blinkonsets(i):blinkonsets(i)+blinkwindowlength)=1;
end
errordiamet=errordiamet(1:length(diamet));

%Now also set all remaining zero-periods in diameter record to error in
%the errordiamet vector
errordiamet=or(errordiamet,rawdiamet==0);

%Now check if there are too small duration pieces of measurement and
%also remove those.
ngood=0;
for i=1:length(diamet)
    if errordiamet(i)==0
        ngood=ngood+1;
    else
        if ngood>0
            if ngood<mincontiggood
                %too little good datapoints; throw them away
                errordiamet(i-ngood:i-1)=1;
            end
            ngood=0;
        end
    end
end

%plot(errordiamet,'k')

%Now interpolate all eye-blinks
%first do the start of the data
if errordiamet(1)==1
    firstgood = find(diff(errordiamet)==-1);
    diamet(1:firstgood(1)) = ...
        mean(diamet(firstgood(1)+1:firstgood(1)+interpolationvals));
    errordiamet(1:firstgood(1)) =0;
end

%Then the end; first make sure there's not one single value at the end
%because this may cause a crash
if errordiamet(end)==0&errordiamet(end-1)==1
    errordiamet(end)=1;
end
if errordiamet(end)==1
    lastgood = find(diff(errordiamet)==1);
    diamet(lastgood(end)+1:end) = ...
        mean(diamet(lastgood(end)-(interpolationvals-1):lastgood(end)));
    errordiamet(lastgood(end)+1:end) =0;
end

%And finally all in between

startinterps=find(diff(errordiamet)==1)+1;
endinterps=find(diff(errordiamet)==-1);

for i=1:length(startinterps)
   
    gaplength=(endinterps(i)-startinterps(i))+1;
    startval = mean(diamet(startinterps(i)-interpolationvals:startinterps(i)-1));
    endval = mean(diamet(endinterps(i)+1:endinterps(i)+interpolationvals));
    datapointgap = (endval - startval)/(gaplength+1);
    %replace values in diameter vector for interpolated ones
    for dpoint=1:gaplength
        diamet(startinterps(i)+dpoint-1) = ...
            startval+dpoint*datapointgap;
    end
end