pNumberOfTopCoordinates = 100,
pNormalizationNorm = 2,
pCorrelationMatrixThres = 0.0,
pNumberOfEigenVectors = 40,

% Thresholding on coordinates
fprintf('Thresholding on coordinates...\n');
lXTop=zeros(size(cX));
lXTopNorm=zeros(size(cX));
for lk = 1:size(cX,1);
    [lSort,lSortIdxs]=sort(cX(lk,:));
    lXTop(lk,lSortIdxs(size(cX,2)-pNumberOfTopCoordinates:size(cX,2)))=cX(lk,lSortIdxs(size(cX,2)-pNumberOfTopCoordinates:size(cX,2)));
    lXTopNorm(lk,:)=lXTop(lk,:)/norm(lXTop(lk,:),pNormalizationNorm);
end;

clear lXTop;

lXTopNormSparse=sparse(lXTopNorm); clear lXTopNorm;

% Compute correlation
fprintf('Computing correlations...\n');
lCorr=lXTopNormSparse*lXTopNormSparse';
lCorr=lCorr/max(max(lCorr));

% Threshold coorelation
fprintf('Thresholding correlations...\n');
lCorrVector=lCorr(:);figure;plot(sort(full(lCorrVector)));axis tight;title('Distribution function of correlation');
lCorrThres=full(lCorr);
lCorrThres(find(abs(lCorr)<median(lCorrVector)))=0;
lCorrThresSparse=sparse(lCorrThres);

fprintf('Creating diffusion matrix...\n');
lCorrThresDSparse=sum(lCorrThresSparse,2);
lCorrThresDSparseInv=lCorrThresDSparse.^(-0.5);

% Compute the diffusion matrix
lA=diag(lCorrThresDSparseInv)*lCorrThresSparse*diag(lCorrThresDSparseInv);
lA = (lA+lA')/2;            % Force simmetry to precision

% Compute the eigenvectors
fprintf('Computing eigenvectors...\n');
lOpts.issym=1;lOpts.isreal=1;
[lEigenVecs,lEigenVals]=eigs(lA,pNumberOfEigenVectors,'LA',lOpts);

figure;plot(diag(lEigenVals));title('EigenValues');
figure;scatter3(lEigenVecs(:,2),lEigenVecs(:,3),lEigenVecs(:,4),10,cLabels(:,2),'filled');title('EigenImage (2,3,4)');
figure;scatter3(lEigenVecs(:,5),lEigenVecs(:,6),lEigenVecs(:,7),10,cLabels(:,2),'filled');title('EigenImage (5,6,7)');

vEigenVecs=lEigenVecs;

% Display results with classification
while(true)    
    h=figure;scatter(vEigenVecs{1}(:,2),vEigenVecs{1}(:,3),20,lLaplacianClusterIdxs{1},'filled');
    k = waitforbuttonpress;
    point1 = get(gca,'CurrentPoint');    % button down detected
    finalRect = rbbox;                   % return figure units
    point2 = get(gca,'CurrentPoint');    % button up detected
    point1 = point1(1,1:2);              % extract x and y
    point2 = point2(1,1:2);
    p1 = min(point1,point2);             % calculate locations
    offset = abs(point1-point2);         % and dimensions
    x = [p1(1) p1(1)+offset(1) p1(1)+offset(1) p1(1) p1(1)];
    y = [p1(2) p1(2) p1(2)+offset(2) p1(2)+offset(2) p1(2)];
    hold on
    axis manual
    plot(x,y,'r','linewidth',5)    
    names(find( (vEigenVecs(:,2)>p1(1)) & (vEigenVecs(:,3)>p1(2)) & (vEigenVecs(:,2)<p1(1)+offset(1)) & (vEigenVecs(:,3)<p1(2)+offset(2))),:)
    pause
    close(h);
end;