function [ IDX,cell,C ] = makemapd( scoresmatrix,a,b,dt,st)
%Performs k-means as per MAKEMAP. Applies IMFUSE to input scores together in a 
%high-contrast image, and passes through a Sobel edge filter. Performs
%clustering and then excludes Sobel edge pixels by indexing into cluster
%zero. 
%
%st: sobel threshold. Approx 0.1
%
%NB: MAKEMAPX = MAKE hyperspectral MAP with sobel edge eXclusion
%
%The rationale behind MAKEMAPX is that pixels, which straddle two regions of
%discrete spectral behaviour, will span the inter-cluster regions of a
%cluster plot. Although the clustering algorithm will reproducibly
%negotiate the boundaries between clusters, these pixels cannot be reliably
%assigned to any one region. A raw Sobel filter finds maxima in the contrast
%differential, whilst a contrast-stretched MAT2GRAY(IMFUSE(...,'falsecolor'))
%of scores A and B provides the highest contrast for the selected A and B.
%The practical upshot being that the Sobel filter will pick only the edges
%which present in the current scores.
%
%Edge pixels span the inter-cluster region evenly and relatively diffusely,
%such that their inclusion in the clustering algorithm is of minimal
%effect. Exclusion of edges via zeroing or mean-setting them either forces
%a boundary condition at the mean, skews the clustering, or places them into
%a 'valid' nonzero cluster. Exclusion via setting to NaN reduces the size
%of the output matrices C and D, upon which further functions are
%exquisitely reliant. 
%
%Type HELP MAKEMAP for more information on how to use this.

global nrows rg
mapvars=[a,b,dt,st];
fused=imfuse(reshape(scoresmatrix(:,a),nrows,[]),reshape(scoresmatrix(:,b),nrows,[]));
[sobel_filt]=edge(imadjust(rgb2gray(fused),stretchlim(rgb2gray(fused)),[]),'sobel',st,'nothinning');
sobel_linear=reshape(sobel_filt,size(scoresmatrix,1),[]);
sobel2=(sobel_filt-1)*-1;
mp=size(scoresmatrix,1)/100;
tic;
[cell,IDX,C]=dbscan(scoresmatrix(:,[a,b])',dt,mp);
toc
clustmap=IDX(1,:);
clustmap=reshape(clustmap,nrows,[]);
rIDX=mat2gray(clustmap);

cmap=colourgend(IDX);
cmapr=cmap(:,1);
cmapg=cmap(:,2);
cmapr(sobel_linear==1)=0;
cmapg(sobel_linear==1)=0;
cmap(:,1:2)=[cmapr cmapg];

rg=zeros(256,3);
rg(:,1)=linspace(0,1,256);
rg(:,2)=linspace(1,0,256);

figure
imshow(rIDX);
colormap(rg);
str = sprintf('Cluster Map, Scores (%d,%d); dt = %d',a,b,dt);
title(str);

clustmap(sobel_filt==1)=0;
rIDX=mat2gray(clustmap);
rIDX=imfuse(rIDX,sobel2,'blend');
figure
imshow(rIDX);
colormap(hot);
str = sprintf('Sobel Exclusions, Scores (%d,%d); dt = %d, Sobel Threshold = %f',a,b,dt,st);
title(str);

figure
scatter(scoresmatrix(1:100:end,a),scoresmatrix(1:100:end,b),10,cmap(1:100:end,:),'filled');
str = sprintf('Score %d',a);
xlabel(str)
str = sprintf('Score %d',b);
ylabel(str)
str = sprintf('Cluster Plot, Scores (%d,%d), dt=%d',a,b,dt);
title(str)

putvar(sobel_filt,mapvars)
end