function [f_flatfield,nearest_neighbors]=flatfieldfill(m_flatfield,mask,dist)
% Function to "fill" a region of an image with data from points outside the
% region
%    [image_out,nearest_neighbors]=flatfieldfill(image_in,mask,dist)
% where
% image_in and image_out are the input and exit images
% mask is an image the same size as image_in, 
%        where the pixels to be filled have value 0
% dist is the number of pixels to use for the averaging outside the mask
% nearest_neighbors is a matrix 
if size(mask)~=size(m_flatfield)
    errordlg('Flatfield and mask are not the same dimensions. Please check')
end
f_flatfield=m_flatfield;
nearest_neighbors=0*mask;
[xmax,ymax]=size(m_flatfield);
[m,n]=size(mask);
msgbox('Will commence filling of the image. This may take some time')
progressbar
for ii=1:xmax
   
    for jj=1:ymax
        if mask(ii,jj)==1
            non_masked_points0=[ii jj];
        else
            non_masked_points0=[];
        end
        %find "dist" higher non-masked pixel in the x-direction
        non_masked_points1=find_non_masked_pixels(mask,ii,jj,'+x',dist);
        %find "dist" lower non-masked pixel in the x-direction
        non_masked_points2=find_non_masked_pixels(mask,ii,jj,'-x',dist);
        %find "dist" higher non-masked pixel in the y-direction
        non_masked_points3=find_non_masked_pixels(mask,ii,jj,'+y',dist);
        %find "dist" lower non-masked pixel in the y-direction
        non_masked_points4=find_non_masked_pixels(mask,ii,jj,'-y',dist);
        %putting them all in a matrix
        try
            non_masked_points=[non_masked_points0;non_masked_points1;...
            non_masked_points2;non_masked_points3;non_masked_points4];
        catch
            ii,jj
        end
        if ~isempty(non_masked_points) ...
                && ~isempty(non_masked_points1) && ~isempty(non_masked_points2)...
                && ~isempty(non_masked_points3) && ~isempty(non_masked_points4)
            %tranforming pixel array to 1D array
            non_masked_points1D=non_masked_points(:,1)+(non_masked_points(:,2)-1)*n;
            %non_masked_points1D=non_masked_points1D(:,1);
            % pixel value is the average of all those found pixels
            pix_val=mean(m_flatfield(non_masked_points1D));
            f_flatfield(ii,jj)=pix_val;
            [nearest_neighbors(ii,jj),kk]=size(non_masked_points1D);
        else
            f_flatfield(ii,jj)=NaN;
        end
    end
    progressbar(ii/xmax)
    %[ii jj nearest_neighbors(ii,jj)]
end

function pts=find_non_masked_pixels(mask,ii,jj,direction,find_num)

coun=0;
[m,n]=size(mask);
pts=[];
switch direction
    case '+x'
        incr=1;
        iincr=ii;
        while coun~=-1
            iincr=iincr+incr;
            if iincr>m % have we moved beyond the high-end border of the mask
                coun=-1;
            else
                if mask(iincr,jj)==1 %is this a non-masked pixel..if so we use it
                    coun=coun+1;
                    pts=[pts; iincr,jj];
                    if coun==find_num, coun=-1;end
                else %no it was a masked pixel let's just go on
                    %just some space we may need
                end
            end
        end
    case '-x'
        incr=-1;
        iincr=ii;
        while coun~=-1
            iincr=iincr+incr;
            if iincr<1 % have we moved beyond the low-end border of the mask
                coun=-1;
            else
                if mask(iincr,jj)==1 %is this a non-masked pixel..if so we use it
                    coun=coun+1;
                    pts=[pts; iincr,jj];
                    if coun==find_num, coun=-1;end
                else %no it was a masked pixel let's just go on
                    %just some space we may need
                end
            end
        end
case '+y'
        incr=1;
        iincr=jj;
        while coun~=-1
            iincr=iincr+incr;
            if iincr>n % have we moved beyond the high-end border of the mask
                coun=-1;
            else
                if mask(ii,iincr)==1 %is this a non-masked pixel..if so we use it
                    coun=coun+1;
                    pts=[pts; ii,iincr];
                    if coun==find_num, coun=-1;end
                else %no it was a masked pixel let's just go on
                    %just some space we may need
                end
            end
        end
    case '-y'
        incr=-1;
        iincr=jj;
        while coun~=-1
            iincr=iincr+incr;
            if iincr<1 % have we moved beyond the low-end border of the mask
                coun=-1;
            else
                if mask(ii,iincr)==1 %is this a non-masked pixel..if so we use it
                    coun=coun+1;
                    pts=[pts; ii,iincr];
                    if coun==find_num, coun=-1;end
                else %no it was a masked pixel let's just go on
                    %just some space we may need
                end
            end
        end       
end
