goNCD.m

%--------------------------------------------------------------------------
% goNCD.m
% main script for evaluating RPE cobblestone morphology using NCD
% (c) 2015-2016 Andrew R. Cohen 
%
%
%       This is the source code for the paper:
%
%  R. Joshi, W. Mankowski, M. Winter, J. S. Saini, T. A. Blenkinsop, J. H. Stern, S. Temple and 
%       A. R. Cohen, "Automated measurement of cobblestone morphology for 
%       characterizing Stem cell-Derived Retinal Pigment Epithelial cell cultures".
%
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     The software is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     in file "gnu gpl v3.txt". If not, see <http://www.gnu.org/licenses/>.
%
% 
% Andrew Cohen
% acohen 'at' coe.drexel.edu
% http://bioimage.coe.drexel.edu
%
%--------------------------------------------------------------------------

if isempty(gcp('nocreate'))
    parpool(36)
end
%  TODO - set path to images
ROOT = '...\RPE\2014-01-24\'
tic

if ~exist('.\scratch','dir')
    mkdir('.\scratch');
end
    
    
% pre-process image settings for NCD w/ new matlab filtering (no cudamex)
FS=1,SF=10,MF=5

% read in and pre-process images
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=gpuArray(im); % can remove this - just for speed
        im=PreProcess(im,FS,MF,SF);
        im=gather(im);
        
        nt=[];
        nt.im=im;
        nt.category = category;
        nt.bCount=-Inf; % init'ed in parallel below
        Trellis=[Trellis nt];
    end
end
% initialize image counts
parfor i=1:length(Trellis)
    Trellis(i).bCount = getCount(Trellis(i).im,i);
end

% cache distances for training sets
dCategory=[];
parfor i=1:5
    ci=find([Trellis.category]==i);
    dCategory(i) = NCD(Trellis(ci),i);
end

% compute distances to each category for each image
d=[];
parfor i=1:length(Trellis)
    tTrain = Trellis;
    % remove test element from training set
    tTrain(i)=[];
    % this is our test element
    tTest = Trellis(i);
    for nClass=1:5
        idx = find([tTrain.category]==nClass);
        % compute NCD(Ax)-NCD(A) for each A, assign to minimum
        %         d(i,nClass)=NCD([tTest tTrain(idx)],i)-NCD([tTrain(idx)],i);
        % avoid recomputing NCD(A)
        if nClass==tTest.category
            d(i,nClass)=dCategory(nClass)-NCD([tTrain(idx)],i);
        else
            d(i,nClass)=NCD([tTest tTrain(idx)],i)-dCategory(nClass);
        end
    end
end
% done - score results
[mm idx]=min(d,[],2);
ncorrect=length(find(idx==[Trellis.category]'))

toc