function [nF, nM, nMF] = centrosome(imgfile, flt, varargin)
% [nF, nM, nMF] = centrosome(imgfile, flt, varargin) 
% 
% [nF, nM, nMF] = centrosome('myfile.TIF', flt, 'countCells', 1, 'plotResults', 0); 
% Leon Peshkin 2005 Harvard Medical School 
countCells = 0; 
plotResults = 0; 
for ndx = 1:2:nargin-2
    switch varargin{ndx}
        case 'countCells'
            countCells = varargin{ndx+1};
        case 'plotResults'
            plotResults = varargin{ndx+1};
        otherwise 
    end
end
TH_F = .90 ;     %  Correlation thresholds for determining Faulty mitosis 
TH_M = .78 ;       %  Normal Mitosis  and 
TH_C = .7 ;         %  cells 
idntyDstTH = 15; % if separated by less then that many pixels, two cells are the same 
img = double(imread(imgfile)); 
[rowsI, colsI] = size(img);
[rowsF, colsF] = size(flt.F1); 
shftR = floor(rowsF/2); shftC = floor(colsF/2); 
shftR = rowsF-1; shftC = colsF-1; 

cmF1 = normxcorr2(flt.F1, img); 
cmF2 = normxcorr2(flt.F2, img); 
cmF3 = normxcorr2(rot90(flt.F2), img); 
cmF4 = normxcorr2(rot90(flt.F1), img); 
cmM1 = normxcorr2(flt.M1, img); 
cmM2 = normxcorr2(flt.M2, img);               % \,
cmM3 = normxcorr2(rot90(flt.M1), img); 
cmM4 = normxcorr2(rot90(flt.M2), img);      %  /
cmM5 = normxcorr2(rot90(flt.M2,2), img);    %  '\
ResXRange = rowsF-shftR:rowsI+rowsF-1-shftR; ResYRange = colsF-shftC:colsI+colsF-1-shftC; 
maskF = (cmF1(ResXRange, ResYRange) > TH_F) | ...  
        (cmF2(ResXRange, ResYRange) > TH_F) | ...
        (cmF4(ResXRange, ResYRange) > TH_F) | ...
        (cmF3(ResXRange, ResYRange) > TH_F) ; 
lociF = bwmorph(maskF,'shrink',Inf);
[lociFx lociFy] = find(lociF);
nF = size(lociFx, 1); 
maskM = (cmM1(ResXRange, ResYRange) > TH_M) | ...
    (cmM2(ResXRange, ResYRange) > TH_M) | (cmM4(ResXRange, ResYRange) > TH_M) | ... 
    (cmM5(ResXRange, ResYRange) > TH_M) | (cmM3(ResXRange, ResYRange) > TH_M);
lociM = bwmorph(maskM,'shrink',Inf);
[lociMx lociMy] = find(lociM);
nM = size(lociMx, 1);   nMF = 0; 
   %  see how many cells were identified both as abnormal and normal mitotic 
if ( (nM > 2) & (nF > 2) )
    DelaunayTriM = delaunay(lociMx, lociMy);    % Auxiliary structure to find a closest point
    if ~isempty(DelaunayTriM)
        clstNdx   = dsearch(lociMx, lociMy, DelaunayTriM, lociFx, lociFy); 
        clstDst   = sqrt(sum(([lociFx'; lociFy'] - [lociMx(clstNdx)'; lociMy(clstNdx)']).^2))';
        nMF = size(find(clstDst < idntyDstTH),1);
        sprintf('%s has %d cells identified as both Abnormal (%d) and Normal (%d) mitotic \n', imgfile, nMF, nF, nM)
        nM = nM - nMF; 
    end
end
if countCells
    cmC1 = normxcorr2(flt.C1, img); 
    cmC2 = normxcorr2(flt.C2, img); 
    cmC3 = normxcorr2(flt.C3, img); 
    cmC4 = normxcorr2(flt.C4, img); 
    cmC5 = normxcorr2(rot90(flt.C1,2), img); 
    maskC = (cmC1(ResXRange, ResYRange) > TH_C) | (cmC2(ResXRange, ResYRange) > TH_C) | ...
            (cmC4(ResXRange, ResYRange) > TH_C) | (cmC3(ResXRange, ResYRange) > TH_C);
    nC = size(find(bwmorph(maskC,'shrink',Inf)), 1);
end 
if plotResults 
    figure; imagesc(img .* ~maskF); title(['Found ' num2str(nF) ' faulty mithosis sites'], 'FontSize', 18)
    set(gcf,'Position', [1 1 colsI rowsI]); 
    figure; imagesc(img .* ~maskM); title(['Found ' num2str(nM) ' normal mithosis sites'], 'FontSize', 18)
    if countCells  
       figure; imagesc(img .* ~maskC); title(['Found ' num2str(nC) ' cells'], 'FontSize', 18) 
    end 
    nMvec=[]; nCvec=[]; nFvec=[]; rng = .6:.05:.95;
    for TH = rng
        maskM = (cmM1(ResXRange, ResYRange) > TH) | (cmM2(ResXRange, ResYRange) > TH) | ... 
             (cmM4(ResXRange, ResYRange) > TH) | (cmM5(ResXRange, ResYRange) > TH) | ... 
             (cmM3(ResXRange, ResYRange) > TH);
        nMvec(end+1) = size(find(bwmorph(maskM,'shrink',Inf)),1); 
        maskF = (cmF1(ResXRange, ResYRange) > TH) |  (cmF2(ResXRange, ResYRange) > TH) | ...
             (cmF4(ResXRange, ResYRange) > TH) | (cmF3(ResXRange, ResYRange) > TH) ; 
        nFvec(end+1) = size(find(bwmorph(maskF,'shrink',Inf)),1); 
        if countCells 
            maskC = (cmC1(ResXRange, ResYRange) > TH) | (cmC2(ResXRange, ResYRange) > TH) | ...
                 (cmC4(ResXRange, ResYRange) > TH) | (cmC3(ResXRange, ResYRange) > TH);
            nCvec(end+1) = size(find(bwmorph(maskC,'shrink',Inf)),1); 
        end
   end 
   figure; hold on, plot(nMvec,nFvec); title(['ROC curve for ' imgfile]), xlabel('# of Mitotic cells'), ylabel('# of Abnormal Mitoses')
   if countCells, plot(nCvec, nMvec + nFvec, 'r-'); end 
end
