clear all
close all
clc

%% initial settings
% set path
dataFolderPath = [pwd,'\'];
fileId = dir([dataFolderPath '\*.tif']);

for i = 1:numel(fileId)
    thisFile = fileId(i).name;
    tmpEu = strmatch('Eu',thisFile);
    if ~isempty(tmpEu)
        if ~exist('imagePathEu','var')
            imagePathEu = [dataFolderPath thisFile];
        else
            error('found more than one Eu file');
        end
    else      
        tmpCo = strmatch('Co',thisFile);
        if ~isempty(tmpCo)
            if ~exist('imagePathCo','var')
                imagePathCo = [dataFolderPath thisFile];
            else
                error('found more than one Co file');
            end
        else
            tmpSi = strmatch('SG',thisFile);
            if ~isempty(tmpSi)
                if ~exist('imagePathSi','var')
                    imagePathSi = [dataFolderPath thisFile];
                else
                    error('found more than one Si file');
                end
            end
        end
    end
end

outputPath = [dataFolderPath 'Ouputs'];

mkdir(outputPath)
outputTxtPath = [outputPath '\result_' date '.txt'];
% set interested regions
ROI_UpperLeft = [350,550];
ROI_LowerRight = [1550,1750];
% set the range of each class
interestedRangeEu = [  8000, 40000;...
                       8000, 40000;...
                       3300,  8000;...
                       8000, 40000;...
                       3000,  7000;...
                      ];
interestedRangeCo = [  1000,  2000;...
                       2000,  5000;...
                       2000,  6000;...
                       6000, 40000;...
                       6000, 40000;...
                    ];
if size(interestedRangeEu,1) ~= size(interestedRangeCo,1)
    error('two ranges have different number of rows');
else
    numClass = size(interestedRangeEu,1);
end
% set the size of ball
radius = 2;

%% processing
% load images
temp = double(imread(imagePathEu));
Eu = temp(ROI_UpperLeft(1):ROI_LowerRight(1),ROI_UpperLeft(2):ROI_LowerRight(2));
clear temp
temp = double(imread(imagePathCo));
Co = temp(ROI_UpperLeft(1):ROI_LowerRight(1),ROI_UpperLeft(2):ROI_LowerRight(2));
clear temp
temp = double(imread(imagePathSi));
Si = temp(ROI_UpperLeft(1):ROI_LowerRight(1),ROI_UpperLeft(2):ROI_LowerRight(2));
clear temp
% band-pass filters
for r = 1:numClass
    thisEuRange = interestedRangeEu(r,:);
    thisCoRange = interestedRangeCo(r,:);
    
    tmpBpEu = (Eu>=min(thisEuRange) & Eu<=max(thisEuRange));
    tmpBpEu = double(imopen(tmpBpEu,strel('disk',radius)));
    
    thisLabel = bwlabel(tmpBpEu); %
    thisProps = regionprops(thisLabel,'Centroid');
    thisCentroid = round([thisProps.Centroid]);
    centerMap = zeros(size(thisLabel));
    for t = 1:numel(thisCentroid)/2
        thisX = thisCentroid(t*2-1);
        thisY = thisCentroid(t*2);
        centerMap(thisY,thisX) = 1;
    end
    ballMap = imdilate(centerMap,strel('disk',radius));
    
    [thisLabel,N] = bwlabel(ballMap);
%     figure,imagesc(thisLabel),axis image;
    
    for n = 1:N
        clear checkList
        thisIdx = find(thisLabel == n);
        thisCoBall = Co(thisIdx);
        % selection rule
        checkList = (thisCoBall>=min(thisCoRange)) & (thisCoBall<=max(thisCoRange));
        if sum(checkList)<length(checkList)*.95
            %bad ball 
            ballMap(thisIdx) = 0;
        end
    end
    figure,imagesc(ballMap),axis image;
    classImg{r} = ballMap;
end

% lookup values in Signal Image
for k = 1:numClass
    clear regMean touchedIdx1 touchedIdx2 touchedIdx
    regMean = [];
    [nextLabel,numLabel] = bwlabel(classImg{k});

    if numLabel>0
        localSingle =0;
        for t = 1:numLabel
            clear regIdx
            regIdx = find(nextLabel == t);
            if ~isempty(regIdx)                
                regMean(end+1) = mean(Si(regIdx));
                localSingle(t)= regIdx(1);%
            end
        end
        localMatrix{k} = localSingle;
    else
        regMean = [];
    end
    
    numLabelMatrix(k) = numLabel;

    if ~isempty(regMean)
        regStd = std(regMean);
        selectedIdx = find(abs(regMean-mean(regMean))<3*regStd);
        newMean = regMean(selectedIdx);
        classProps(k).rangeEu = [interestedRangeEu(k,:)];
        classProps(k).rangeCo = [interestedRangeCo(k,:)];
        classProps(k).regionNumOriginal = numLabel;
        classProps(k).regionsMean = regMean;
        classProps(k).regionNumRefined = numel(newMean);
        classProps(k).classMean = mean(newMean);
        classProps(k).classStd = std(newMean);  
    else
        
        classProps(k).rangeEu = [interestedRangeEu(k,:)];
        classProps(k).rangeCo = [interestedRangeCo(k,:)];
        classProps(k).regionNumOriginal = numLabel;
        classProps(k).regionsMean = 0;
        classProps(k).regionNumRefined = 0;
        classProps(k).classMean = 0;
        classProps(k).classStd = 0;
    end
end


%% output
fid = fopen(outputTxtPath,'w');
fprintf(fid,'**************************Settings***************************\n');
% fprintf(fid,'folder  directory :    %s\', [dataFolderPath ' ']);
fprintf(fid,'region of interest:    (%d, %d)-(%d, %d)\n',ROI_UpperLeft(1), ROI_UpperLeft(2), ROI_LowerRight(1),ROI_LowerRight(2));
fprintf(fid,'Eu image information : %s\n',imagePathEu);
fprintf(fid,'Co image information : %s\n',imagePathCo);
fprintf(fid,'Sg image information : %s\n',imagePathSi);
fprintf(fid,'ball radius       :     %d\n',radius);
fprintf(fid,'*************************************************************\n');
for k = 1:numClass
    fprintf(fid,'\nClass Type     %d\n', k);
    fprintf(fid,'Class Range Eu    %d - %d\n', classProps(k).rangeEu(1),classProps(k).rangeEu(2));
    fprintf(fid,'Class Range Co    %d - %d\n', classProps(k).rangeCo(1),classProps(k).rangeCo(2));
    fprintf(fid,'Ball  Num      %d\n', classProps(k).regionNumRefined);
    fprintf(fid,'Ball  Mean     %5.1f\n', classProps(k).classMean);
    fprintf(fid,'Ball  Std      %5.1f\n', classProps(k).classStd);
    imwrite(classImg{k},[outputPath '\Class Type = ' num2str(k) '.tiff']);
end    
fclose(fid);

display('!!!!!!!!!!!!!!!!!!!done!!!!!!!!!!!!!!!!!!!');

    [row, col] = size(Si);
    
    for kk=1:length(localMatrix)
        Sii = Si;
        regIdx =  localMatrix{kk};
        regMean = classProps(kk).regionsMean; 
        regMean = round(regMean);
        Sii(regIdx) = 65536;
       [ii,jj] = ind2sub([row,col],regIdx);
        figure,imshow(Sii,[]);
        title(['Type ',num2str(kk),' ball'])
        for mm = 1: numLabelMatrix(kk)
          hold on; text(jj(mm),ii(mm),num2str(regMean(mm)),'Color','r') 
        end
        
    end
    
    outputTxtPath_detail = [outputPath '\result_detail' date '.txt'];
    
    fid = fopen(outputTxtPath_detail,'w');
    fprintf(fid,'**************************Settings***************************\n');
  % fprintf(fid,'folder  directory :    %s\', [dataFolderPath ' ']);
    fprintf(fid,'region of interest:    (%d, %d)-(%d, %d)\n',ROI_UpperLeft(1), ROI_UpperLeft(2), ROI_LowerRight(1),ROI_LowerRight(2));
    fprintf(fid,'Eu image information : %s\n',imagePathEu);
    fprintf(fid,'Co image information : %s\n',imagePathCo);
    fprintf(fid,'Sg image information : %s\n',imagePathSi);
    fprintf(fid,'ball radius       :     %d\n',radius);
    fprintf(fid,'*************************************************************\n');
    for k = 1:numClass
       regMean = classProps(k).regionsMean; 
       regMean = round(regMean);
       fprintf(fid,'\nClass Type     %d\n', k);
       fprintf(fid,'Class Range Eu    %d - %d\n', classProps(k).rangeEu(1),classProps(k).rangeEu(2));
       fprintf(fid,'Class Range Co    %d - %d\n', classProps(k).rangeCo(1),classProps(k).rangeCo(2));
       fprintf(fid,'Ball  Num      %d\n', classProps(k).regionNumRefined);
       fprintf(fid,'Mean intensity for each ball \n'); 
         for ii = 1:numLabelMatrix(k)
             fprintf(fid,'%8d',regMean(ii));
             if mod(ii,10)==0
             fprintf(fid,'\n');
             end
         end
   
         fprintf(fid,'\n');
         fprintf(fid,'\n');
      
    end    
   fclose(fid);
    
    
