diff --git a/NCD/CudaMex.mexw64 b/NCD/CudaMex.mexw64
new file mode 100644
index 0000000000000000000000000000000000000000..591dc985ca02398dea6b266fa0a1498d5195e0a7
Binary files /dev/null and b/NCD/CudaMex.mexw64 differ
diff --git a/NCD/NCD.m b/NCD/NCD.m
new file mode 100644
index 0000000000000000000000000000000000000000..da5fbbcc84597eedb8c2ad7d7ab7ce9bef0b1a9e
--- /dev/null
+++ b/NCD/NCD.m
@@ -0,0 +1,22 @@
+function d=NCD(rg,index)
+
+% whole thing
+
+imX = getCombinedImage(rg);
+
+% imagesc(imX)
+% drawnow
+% axis image
+GX = getCount(imX,index);
+
+gx = min([rg.bCount]);
+
+gExclude =[];
+for i=1:length(rg)
+ rgExclude=rg;
+ rgExclude(i)=[];
+ imExclude = getCombinedImage(rgExclude);
+ gExclude(i)=getCount(imExclude,index);
+end
+
+d = (GX - gx)/(max(gExclude));
diff --git a/NCD/PreProcess.m b/NCD/PreProcess.m
new file mode 100644
index 0000000000000000000000000000000000000000..105db46055ca25fd139491fca5aa8f01124f769d
--- /dev/null
+++ b/NCD/PreProcess.m
@@ -0,0 +1,14 @@
+function imf=PreProcess(im)
+
+% h = fspecial('gaussian',7,0.5);
+% imfilt = imfilter(im, h, 'symmetric');
+% for j=1:50
+% imfilt = imfilter(imfilt, h, 'symmetric');
+% end
+
+imfilt = imgaussfilt(im);
+imhfreq = max((im - imfilt), zeros(size(im)));
+imf = medfilt2(imhfreq,[3 3]);
+imf=stdfilt(imf,getnhood(strel('disk',5)));
+imf = mat2gray(imf);
+
diff --git a/NCD/compress_jp2.m b/NCD/compress_jp2.m
new file mode 100644
index 0000000000000000000000000000000000000000..91d8fd7023ad5605bf9c6cff6045f223ebfc744a
--- /dev/null
+++ b/NCD/compress_jp2.m
@@ -0,0 +1,130 @@
+% Compress using jpeg2000
+% pairwise NCD
+
+clear;
+close all;
+tic;
+% system('del im*.dat*');
+
+input = importdata_phase();
+%% Try all the images
+% % Category 1
+% input{1} = rgb2gray(imread('2014-01-24\Cobblestone category 1 image c_32X.tif'));
+% input{2} = rgb2gray(imread('2014-01-24\Cobblestone category 1 image b_32X.tif'));
+% input{3} = rgb2gray(imread('2014-01-24\Cobblestone category 1 image a _32X.tif'));
+% input{4} = rgb2gray(imread('2014-01-24\Category 1\Catgeroy 1_32X_A.tif'));
+%
+% % Category 2
+% input{5} = rgb2gray(imread('2014-01-24\Cobblestone category 3 image a_32X.tif'));
+% input{6} = rgb2gray(imread('2014-01-24\Cobblestone category 2 image b_32X.tif'));
+% input{7} = rgb2gray(imread('2014-01-24\Category 2\Catgeroy 2_32X_A.tif'));
+% input{8} = rgb2gray(imread('2014-01-24\Category 2\Catgeroy 2_32X_B.tif'));
+% input{9} = rgb2gray(imread('2014-01-24\Category 2\Catgeroy 2_32X_C.tif'));
+%
+% % Category 3
+% input{10} = rgb2gray(imread('2014-01-24\Cobblestone category 2 image a_32X.tif'));
+% input{11} = rgb2gray(imread('2014-01-24\Cobblestone category 3 image b_32X.tif'));
+% input{12} = rgb2gray(imread('2014-01-24\Category 3\Category 3_32X_A.tif'));
+% input{13} = rgb2gray(imread('2014-01-24\Category 3\Category 3_32X_B.tif'));
+% input{14} = rgb2gray(imread('2014-01-24\Category 3\Category 3_32X_C.tif'));
+%
+% % Category 5
+% input{15} = rgb2gray(imread('2014-01-24\Cobblestone category 5 image c_32X.tif'));
+% input{16} = rgb2gray(imread('2014-01-24\Cobblestone category 5 image b_32X.tif'));
+% input{17} = rgb2gray(imread('2014-01-24\Cobblestone category 5 image a_32X.tif'));
+% % letters =['A' 'B' 'C' 'D' 'E' 'F' 'G' 'H' 'I' 'J' 'K' 'L' 'M'];
+% % letters =['A' 'B' 'C' 'D' 'F' 'G' 'H' 'I' 'J' 'L' 'M'];
+% letters = ['C' 'J' 'M'];
+% % letters = ['C' 'M'];
+% for i=1:numel(letters)
+% filename = ['2014-01-24\Category 5\Catgory 5_32X_' letters(i) '.tif'];
+% input{17+i} = rgb2gray(imread(filename));
+% end
+
+
+%% Computation
+
+for i=1:size(input,2)
+% input{i} = mat2gray(input{i});
+% input{i} = imresize(input{i}, 0.5);
+% input{i} = (stdfilt(input{i}, getnhood(strel('disk',3))));
+ [m n] = size(input{i});
+ input{i} = input{i}(1:m/2, 1:n/2);
+ input{i} = CudaMex('ContrastEnhancement', input{i}, [5 5 1],[3 3 1]);
+% if size(input{i},3) == 3
+% input{i} = rgb2gray(input{i});
+% end
+% % input{i} = imadjust(input{i});
+ input{i} = mat2gray(stdfilt(input{i}, getnhood(strel('disk',5))));
+
+ fname = ['im' num2str(i) '.jp2'];
+ imwrite(input{i}, fname, 'jp2', 'Mode', 'lossless','Tilesize', size(input{i})/2);
+ fileinfo = dir(fname);
+ C(i) = fileinfo.bytes;
+end
+
+% comb1 = [input{1}; input{2}];
+% imwrite(comb1, 'comb1.jp2');%, 'Mode', 'lossless');
+% fileinfo = dir('comb1.jp2');
+% s1 = fileinfo.bytes;
+% mn = min(C(1),C(2));
+% mx = max(C(1),C(2));
+% ncd(1,2) = (s1 - mn)/mx;
+%
+% comb2 = [input{1}; input{3}];
+% imwrite(comb2, 'comb2.jp2');%, 'Mode', 'lossless');
+% fileinfo = dir('comb2.jp2');
+% s2 = fileinfo.bytes;
+% mn = min(C(1),C(3));
+% mx = max(C(1),C(3));
+% ncd(1,3) = (s2 - mn)/mx;
+%
+% comb3 = [input{1}; input{4}];
+% imwrite(comb3, 'comb3.jp2');%, 'Mode', 'lossless');
+% fileinfo = dir('comb3.jp2');
+% s3 = fileinfo.bytes;
+% mn = min(C(1),C(4));
+% mx = max(C(1),C(4));
+% ncd(1,4) = (s3 - mn)/mx;
+format long;
+
+for i=1:size(input, 2)
+ for j=1:size(input,2)
+ comb = [input{i} input{j}];
+ fname = ['comb' num2str(i) '_' num2str(j) '.jp2'];
+ imwrite(comb, fname,'jp2' , 'Mode', 'lossless','Tilesize', size(input{3}));
+ fileinfo = dir(fname);
+ s(i,j) = fileinfo.bytes;
+
+% ncd(i,j) = (size(CE{k},1) - min(C(i),C(j)))/max(C(i),C(j));
+ end
+end
+
+for i=1:size(input, 2)
+ for j=1:size(input,2)
+ mn = min(C(i),C(j));
+ mx = max(C(i),C(j));
+ ncd(i,j) = (s(i,j) - mn)/mx;
+ end
+end
+[kGap Gap S idx] = GapSpectral(ncd, 7);
+[idx2 Y] = SpectralCluster(ncd, 2);
+[idx4 Y] = SpectralCluster(ncd, 4);
+[idx5 Y] = SpectralCluster(ncd, 5);
+% for i=1:size(input, 2)
+% for j=1:size(input,2)
+% comb3(:,:,1) = input{i};
+% comb3(:,:,2) = input{j};
+% comb3(:,:,3) = zeros(size(input{i}));
+% fname = ['comb3' num2str(i) '_' num2str(j) '.jp2'];
+% imwrite(comb3, fname,'jp2', 'Mode', 'lossless');
+% fileinfo = dir(fname);
+% s3(i,j) = fileinfo.bytes;
+% mn = min(C(i),C(j));
+% mx = max(C(i),C(j));
+% ncd3(i,j) = (s3(i,j) - mn)/mx;
+% % ncd(i,j) = (size(CE{k},1) - min(C(i),C(j)))/max(C(i),C(j));
+% end
+% end
+%
+toc;
diff --git a/NCD/getCombinedImage.m b/NCD/getCombinedImage.m
new file mode 100644
index 0000000000000000000000000000000000000000..27008cea70e751757f45d6febc743462ac1b9bbb
--- /dev/null
+++ b/NCD/getCombinedImage.m
@@ -0,0 +1,12 @@
+function im = getCombinedImage(rg);
+im=[];
+for i=1:1 %length(rg)
+ rowIm=[];
+ for j=1:length(rg)
+ rowIdx = mod(i+j-2,length(rg))+1;
+% lay(i,j)=rowIdx;
+ rowIm=[rowIm rg(rowIdx).im];
+ end
+ im=[im;rowIm];
+end
+% im=im2bw(im,graythresh(im));
\ No newline at end of file
diff --git a/NCD/getCount.m b/NCD/getCount.m
new file mode 100644
index 0000000000000000000000000000000000000000..28e36e6e0482d128e49941e21c823832e0fe8c1a
--- /dev/null
+++ b/NCD/getCount.m
@@ -0,0 +1,14 @@
+function count = getCount(im,parkey)
+
+if nargin<2
+ tag='';
+else
+ tag=num2str(parkey);
+end
+im=mat2gray(im);
+im=gray2ind(im);
+fname = ['e:\scratch\scratch_' tag '.jp2'];
+imwrite(im,fname,'mode','lossless');
+info=imfinfo(fname);
+count=info.FileSize;
+
diff --git a/NCD/goNCD.m b/NCD/goNCD.m
new file mode 100644
index 0000000000000000000000000000000000000000..47d4cd898b937cec7dde9781795d44785ac62a72
--- /dev/null
+++ b/NCD/goNCD.m
@@ -0,0 +1,53 @@
+
+if isempty(gcp('nocreate'))
+ parpool(36)
+end
+ROOT = 'F:\Raw\Images\Temple\RPE\2014-01-24\'
+tic
+Trellis=[];
+for category=1:5
+ droot = [ROOT 'Category ' num2str(category) '\'];
+ flist = dir([droot '*.tif']);
+ for ff=1:length(flist)
+
+ im = imread([droot flist(ff).name]);
+ im = rgb2gray(im);
+% im=im2double(im);
+
+% im = CudaMex('ContrastEnhancement', im, [5 5 1],[3 3 1]);
+% im = mat2gray(stdfilt(im, getnhood(strel('disk',5))));
+%
+
+ im=gpuArray(im);
+ im=PreProcess(im);
+ im=gather(im);
+ [m n] = size(im);
+ im = im(1:m/2, 1:n/2);
+ nt=[];
+ nt.im=im;
+ nt.category = category;
+ nt.bCount=-1;
+ Trellis=[Trellis nt];
+ end
+end
+
+parfor i=1:length(Trellis)
+ Trellis(i).bCount = getCount(Trellis(i).im,i);
+end
+
+
+d=[];
+parfor i=1:length(Trellis)
+ tTrain = Trellis;
+ tTrain(i)=[];
+ tTest = Trellis(i);
+ for nClass=1:5
+ idx = find([tTrain.category]==nClass);
+ d(i,nClass)=NCD([tTrain(idx) tTest],i)-NCD([tTrain(idx)],i);
+ end
+end
+
+[mm idx]=min(d,[],2);
+length(find(idx==[Trellis.category]'))
+
+toc
\ No newline at end of file