function dd = datadensity(x,y,varargin)
% datadensity Computes the data density (points/area) of scattered points
%
% Usage:
%   dd = datadensity(x,y,'method','squares','radius',0.2)
%   dd = datadensity(x,y,'method','voronoi')
%
% Input:
%   (x,y) -  coordinates of points
%   method - either 'squares','circles', or 'voronoi'
%       default = 'voronoi'
%   radius - Equal to the circle radius or half the square width. 
%          The default range is max(range(x),range(y))*0.05
%
% Output:
%   dd - Data density
%
%Copy-Left, Alejandro Sanchez, 2012
if nargin==0
    dd = rundemo;
    return
end
method = 'vo';
r = [];
for k=1:2:length(varargin)-1
    switch lower(varargin{k})
        case 'method'
            method = varargin{k+1};
        case 'radius'
            r = varargin{k+1};    
        otherwise
            warning(['Unidentified Property: ',varargin{k}])
    end
end
if isempty(r)
    r = max(range(x),range(y))*0.05;
end
%Correct data if necessary
x = x(:); y = y(:);
%Asuming x and y match
idat = isfinite(x);
x = x(idat); y = y(idat);
%Initialize
Ld = length(x);
dd = zeros(Ld,1);
%Calculate Data Density
switch method(1:2)
    case 'sq'  %---- Using squares ----
        for k=1:Ld
            dd(k) = sum( x>(x(k)-r) & x<(x(k)+r) & y>(y(k)-r) & y<(y(k)+r) );
        end %for
        area = (2*r)^2;
        dd = dd/area;
    case 'ci'
        for k=1:Ld
            dd(k) = sum( sqrt((x-x(k)).^2 + (y-y(k)).^2) < r );
        end
        area = pi*r^2;
        dd = dd/area;
    case 'vo'  %----- Using voronoi cells ------
        [v,c] = voronoin([x,y]);     
        for k=1:length(c) 
            %If at least one of the indices is 1, 
            %then it is an open region, its area
            %is infinity and the data density is 0
            if all(c{k}>1)   
                a = polyarea(v(c{k},1),v(c{k},2));
                dd(k) = 1/a;
            end %if
        end %for
    otherwise
        error('Invalid Method')
end %switch
%Relocate variables and place NaN's
dd(idat) = dd;
dd(~idat) = NaN;
return
%-----------------------------------------------------
function dd = rundemo
%The two lines below can be interchanged
% method = 'voronoi'; 
method = 'squares'; r = [];
x = randn(1000,1);
y = randn(1000,1);
clf
hold on
if strcmpi(method,'voronoi')
    dd = datadensity(x,y,'method',method);    
elseif strcmpi(method,'squares')
    dd = datadensity(x,y,'method',method,'radius',r);
end
%Bidimensional running mean filter
N = 100;
xi = repmat(linspace(min(x(:)),max(x(:)),N),N,1);
yi = repmat(linspace(min(y(:)),max(y(:)),N)',1,N);
F = TriScatteredInterp(x,y,dd); %Construct interpolant
zi = F(xi,yi); %Interpolate scattered data
n = 5; %running mean filter width
c = ones(n(1),1)/n(1); %running mean filter
zi(isnan(zi)) = 0;
zi = conv2(c,c,zi,'same'); 
dds = interp2(xi,yi,zi,x,y);
contour(xi,yi,zi)
gsp(x,y,dds,5)
hold off
colorbar
return
%------------------------------------------------------
function varargout = gsp(x,y,c,ms)
%Graphs scattered poits
map = colormap;
ind = fix((c-min(c))/(max(c)-min(c))*(size(map,1)-1))+1;
%much more efficient than matlab's scatter plot
h = zeros(size(map,1),1);
for k=1:size(map,1) 
    if any(ind==k)
        h(k) = line('Xdata',x(ind==k),'Ydata',y(ind==k), ...
            'LineStyle','none','Color',map(k,:), ...
            'Marker','.','MarkerSize',ms);
    end
end
h(h==0) = [];
if nargout==1
    varargout = h; 
end
return

