function [vWeights,vMetric] = GetWeights_AnisotropicGaussian( cPoints, cPointsIdxs,cCenter, cExtraParam, FixedMetric )

%
% function [vWeights,vMetric] = GetWeights_AnisotropicGaussian( cPoints, cPointsIdxs,cCenter, cExtraParam )
%
% IN:
%   cPoints     : M points in N dimensions in a M by N matrix (typically the whole dataset)
%   cPointsIdxs : K indices (into cPoints) of the points involved in the weights computation
%   cCenter     : Row vector representing the point from which to compute the weights to the cPoints(cPointsIdxs,:)
%   cExtraParam : contains the fields:
%                   Delta : width of the gaussian weight to superimpose to the computed distances
%                   Epsilon : radius of the ellipsoid
%  FixedMetric  : (default) empty array; if not empty, the metric is fixed and given by N-by-N matrix 
%
% OUT:
%   vWeights    : K row vector, i-th entry being the weight from cCenter to cPoints(cPointsIdxs(lk),:)
%   vMetric     : Compute local pca metric. 
%                       vMetric.Cov is covariance matrix of the points.
%                       vMetric.InvCov is pseudo-inverse of vMetric.Cov
%
% USES:
%   svd
%

vMetric = []; vWeights = [];
lNumberOfPoints = length(cPointsIdxs);

% Check if metric is preset (Added, AL 4/24/04)  % NEW!
if isempty(FixedMetric)  % estimate the covariance matrix
        
    % Center the points around cCenter
    lCenteredPoints = cPoints(cPointsIdxs,:);
    %lMean = mean(cPoints,1);
    for lk = 1:lNumberOfPoints
        lCenteredPoints(lk,:) = lCenteredPoints(lk,:) - cCenter;%lMean;
    end;
    
    % Compute the covariance matrix
    try
        lCov = lCenteredPoints'*lCenteredPoints;
        % Compute the distances
        lInvCov = pinv(lCov);
        %lCenteredCenter = cCenter-lMean;
        lDistances=zeros(1,lNumberOfPoints);
        for lk = 1:lNumberOfPoints
            lDistances(lk) = sqrt((lCenteredPoints(lk,:))*(lInvCov)*(lCenteredPoints(lk,:))');
        end;
                
        % Compute the weights
        vWeights = exp(-(lDistances/cExtraParam.Delta).^2);
        % No diffusion outside the ellipsoid of radius cExtraParam.Epsilon
        vWeights(find(lDistances>cExtraParam.Epsilon))=0;
        vMetric.Cov = lCov;
        vMetric.InvCov = lInvCov;
    catch
        vMetric.Cov = lCov;
        vMetric.InvCov = [];
        Weights = zeros(lNumberOfPoints,1);
    end;
    
    %figure(50);cla;plot(cPoints(:,1),cPoints(:,2),'ko');hold on;plot(cCenter(1),cCenter(2),'ro');hold on;scatter(cPoints(cPointsIdxs,1),cPoints(cPointsIdxs,2),15,lDistances,'filled');colorbar;drawnow;pause;
    
else % Use preset metric   NEW!
    
    for lk = 1:lNumberOfPoints
        lCenteredPoint = cPoints(lk,:) - cCenter;
        lDistances(lk) = sqrt(lCenteredPoint*FixedMetric*lCenteredPoint');
    end;
    
    % Compute the weights
    vWeights = exp(-(lDistances/cExtraParam.Delta).^2);
    % No diffusion outside the ellipsoid of radius cExtraParam.Epsilon
    vWeights(find(lDistances>cExtraParam.Epsilon))=0;
    
end

return;