% Copyright (c) 2015 The Johns Hopkins University/
% Applied Physics Laboratory. (JHU/APL)
% 
% This material was developed by JHU/APL under U.S. Government contract.
% The U.S. Government has Unlimited Rights in this material per
% DFARS 252.227-7013/7014, as appropriate.
% 
% For any other permission, please contact the Legal Office at JHU/APL.
%
function [ imgFilled, raster, N] = rasterData(xyIn, cData, outSize, validMask, method, fillValue)

        % METHOD:
        %   1:MAX (default), 2:MIN, 3:MEAN

        if ~exist('method', 'var')  || isempty(method)
            method = 1;
        end


		if ~exist('fillValue', 'var') || isempty(fillValue)
			fillValue = nan;
		end
		
		
        binSpacing = 1;
        z = cData;
        
        Xedges = 0.5:binSpacing:(outSize(2) + 0.5);
        Yedges = 0.5:binSpacing:(outSize(1) + 0.5);

        [N,Xedges,Yedges,binX,binY] = histcounts2(xyIn(:,1), xyIn(:,2), Xedges, Yedges);        
        unbinned = binX == 0;
              
        if isnumeric(method)
            switch method
                case 1  % MAX value
                    raster  = accumarray([binX(~unbinned) binY(~unbinned)], z(~unbinned), [numel(Xedges), numel(Yedges)]-1, @max);
                case 2 % MIN value
                    raster  = accumarray([binX(~unbinned) binY(~unbinned)], z(~unbinned), [numel(Xedges), numel(Yedges)]-1, @min);
                case 3 % AVG value
                   raster  = accumarray([binX(~unbinned) binY(~unbinned)], z(~unbinned), [numel(Xedges), numel(Yedges)]-1);
                   raster = raster./N;
                otherwise
                    error('rasterData.m: Unknown method specified (%g)', method)
            end
            raster(raster == 0) = fillValue;
            
        elseif isa(method, 'functio_nhandle')
            raster  = accumarray([binX(~unbinned) binY(~unbinned)], z(~unbinned), [numel(Xedges), numel(Yedges)]-1, method);
             raster(raster == 0) = fillValue;
        else
            error('rasterData.m: Unknown method specified')
        end
       
        raster = raster';
        
        if exist('validMask', 'var') && all(size(validMask) == size(raster))
            imgFilled = inpaint(raster, validMask);
        else
            imgFilled = raster;
        end
end


function imgFilled = inpaint(inImage, BW)

    BW = double(BW);
    
    imgFilled = inImage;

    toFill = isnan(imgFilled) & BW;

    d = bwdist(~BW);
    f = ones(3); f(2,2) = 0;

    notFilled = 1;
    maxNeigbors =  imfilter(BW, f);

    maxItr = 100; itr = 0;
    while (notFilled > 0) && (itr < maxItr)
        itr = itr+1;

        img = imgFilled;
        img(isnan(img)) = 0;

        count  = imfilter(double(~isnan(imgFilled)), f);
        filling = imfilter(img, f) ./ count;

        pctNeighbors = count./maxNeigbors;

        for p = 0.8:(-0.2):0.4
            fill = pctNeighbors >= p & toFill;
            if any(any(fill))
                break
            end
        end

        if ~any(any(fill))
           notFilled = 0;
           continue;
        end
        imgFilled(fill) = filling(fill);

        toFill = isnan(imgFilled) & BW;

        if ~any(toFill(:))
           notFilled = 0; 
        end
        notFilled = notFilled+1; 
    end

end
