clear all
clf

% image processing
x = imread("2246-2.tif");
figure(1)

imshow(x)

disk = strel('disk', 10);
opened_image = imopen(x, disk);
% imshow(io)

erroded_image = imerode(x, disk);
reconstructed_eroded_image = imreconstruct(erroded_image, x);
% imshow(iobr)

closed_image = imclose(opened_image, disk);
% imshow(ioc)

dilated_image = imdilate(reconstructed_eroded_image, disk);
reconstructed_dilated_image = imreconstruct(imcomplement(dilated_image), imcomplement(reconstructed_eroded_image));
% imshow(iobrcbr)

final = reconstructed_dilated_image;

figure(2)
title("Relative Height")
imshow(final)

[row, col] = size(final);
minpix = 0;
maxpix = 0;

% min/max/average determination
for a = 1:row
    for b = 1:col
        if (final(a,b) > maxpix)
            maxpix = final(a,b);
        elseif (final(a,b) < minpix)
            minpix = final(a,b);
        end
        
    end
end

tablepixel = zeros((maxpix - minpix),2);

for (h = 1:(maxpix - minpix))
    tablepixel(h+1,1) = minpix + h;
end


for a = 1:row
    for b = 1:col
        rowvalue = final(a,b);
        tablepixel(rowvalue, 2) = tablepixel(rowvalue, 2) + 1;
    end
end

maxrelpix = 0;

for a = 0:(maxpix-minpix)
    if (tablepixel((maxpix-minpix)-a,2) > (0.005*row*col))
        maxrelpix = tablepixel((maxpix-minpix)-a,1);
        break
    end
    
    
end


finalchange = process(final, 3.2); 
[row, col] = size(finalchange);

% 
% % self-correcting percentile based calculations on empirical data
% % percentile values determined manually, and double checked after
% % calculations.
% % enable if dealing with bulbous code
%
condition = true;
coverage = 0.6;

while(condition)
    
    
    for a = 1:row
        for b = 1:col
            temp = (maxrelpix - double(final(a,b)))*100/maxrelpix;
            alpha = [4.3 5.5 10]; % values derived from experimental data.
            
            if (temp < alpha(1))
                finalchange(a,b) = 256;
                
            end
            
            
            if (temp > alpha(2) && temp < alpha(3))
                finalchange(a,b) = 256;
            elseif (temp > alpha(1))
                finalchange(a,b) = 0;
            end
            
            
        end
    end
    
    validation = 0;
    for a = 1:row
        for b = 1:col
        if (finalchange(a,b)== 255 || finalchange(a,b) == 256)
            validation = validation + 1;
            
            end
        
        end
    end
    
    validationarea = validation/(row*col);
    validationerror = (validationarea-coverage)/coverage;
    
    if ((abs(validationerror) < 0.10))     
        condition = false;
    else
        coverage = coverage + validationerror*coverage;
        alpha = alpha + alpha*validationerror;
    end
    
end


figure(3)
title("Binary Gradient")
imshow(finalchange)


% Cell Area Determination

figure(4)
imshow(finalchange)
title("Binary Gradient - Overlayed Area")
[centers,radii] = imfindcircles(finalchange,[30 50],'ObjectPolarity','bright', ...
    'Sensitivity',0.97)
viscircles(centers, radii,'Color','b')


[o, p] = size(radii);

averageradii = sum(radii,1)/o;
cellarea = pi*(averageradii)^2

% Cell Counter

cellpix = 0;
for a = 1:row
    for b = 1:col
        if (finalchange(a,b)== 255 || finalchange(a,b) == 256)
            cellpix = cellpix + 1;
            
        end
        
    end
end


totalcoverage = cellpix/(row*col)
cellcount = round(cellpix/cellarea)










