clear
close all

origI = imread('..\HI images\pcr11sep1523h2000ml40s.tif');
%origI = imread('fakeArray.tiff');
imshow(origI);
title('Original Image');

bkgdI = imopen(origI, strel('disk',15));
figure
surf(double(bkgdI(1:8:end,1:8:end)));
set(gca,'ydir','reverse');
title('Background Estimate');

unifBkgdI = origI - bkgdI;
figure
imshow(unifBkgdI);
title('Image w/Uniform Background');

adjustedI = imadjust(unifBkgdI);
figure
imshow(adjustedI);
title('Adjusted Intensity');

bwMask = imread('bwMaskFile.tiff');
bwMask = bwMask / 65535;
maskedI = bwMask .* adjustedI;
figure
imshow(maskedI);
title('Masked Adjusted Image');

level = graythresh(maskedI);
P = 32
figure
hBWI = gcf;
figure
hBWOpenedI = gcf;
while( (level>0) & (P>0) )
    bwI = im2bw(maskedI, level);
    figure(hBWI)
    imshow(bwI)
    title('B/W Image');
    bwOpenedI = bwareaopen(bwI, P);
    figure(hBWOpenedI)
    imshow(bwOpenedI);
    title('B/W Opened Image')
    fprintf('Threshold value P value %f  %f\n', level, P);
    s = input('Enter new threshold and P values (neg value to end): ', 's');
    s2 = textscan(s, '%f %f');
    level = s2{1};
    P = s2{2};
end

figure
test = 'y';
fprintf('Use the cursor to draw a polygon around unwanted spots\n');
maskedBWOpenedI = bwOpenedI;
while test == 'y'
    mask = roipoly(maskedBWOpenedI);
    maskedBWOpenedI = ~mask .* maskedBWOpenedI;
    imshow(maskedBWOpenedI);
    title('Masked BW Opened Image');
    test = input('Again/reset/done (y/r/d)?\n ', 's');
    if test == 'r'
        maskedBWOpenedI = bwOpenedI;
        test = 'y';
    end
end

cc = bwconncomp(bwOpenedI);
imgStats = regionprops(cc, bwOpenedI, 'Area', 'Orientation');
a = [imgStats(:).Area];
o = [imgStats(:).Orientation];
regAreaIndex = find(a == max(a));
regAngle = imgStats(regAreaIndex).Orientation;
rotationAngle = -1.0 * (regAngle + 88.5);
rotBWOpenedI = imrotate(bwOpenedI, rotationAngle);
figure
imshow(rotBWOpenedI)
title('Rotated B/W Opened Image');
rotOrigI = imrotate(origI, -1.0*(regAngle+88.5));
figure
imshow(rotOrigI)
title('Rotated Original Image');
newAngle = 0.00001;
while newAngle ~= 0.0
    newAngle = input('Rotate more? (+/- degrees, 0 = done)\n ', 's');
    if newAngle ~= 0
        rotationAngle = rotationAngle + newAngle;
        rotOrigI = imrotate(rotOrigI, rotationAngle);
    end
end

ccRot = bwconncomp(rotBWOpenedI);
imgStatsRot = regionprops(ccRot, rotOrigI, 'Area', 'WeightedCentroid', 'MeanIntensity');
temp = [imgStatsRot.WeightedCentroid];
centC = temp(1:2:2*length(imgStatsRot));
centR = temp(2:2:2*length(imgStatsRot));
areas = [imgStatsRot.Area];
maxAreaIndex = find(areas == max(areas));
centC = [centC(1:maxAreaIndex-1) centC(maxAreaIndex+1:end)];
centR = [centR(1:maxAreaIndex-1) centR(maxAreaIndex+1:end)];
meanInt = [imgStatsRot.MeanIntensity];
z = sort(centR);
for i = 1:11
    row(i,:) = z((i-1)*8+1:(i*8));
end
for i = 1:11
    for j= 1:8
        col(i,j) = centC( find(centR == row(i,j)) );
    end
end

for i=1:11
    for j = 1:8
        c = floor(col(i,j));
        r = floor(row(i,j));
        mean4X4Vals(i,j) = mean([rotOrigI(r,c) rotOrigI(r,c+1) rotOrigI(r+1,c) rotOrigI(r+1,c+1)])
    end
end
for i=1:11
    for j=1:8
        c = round(col(i,j));
        r = round(row(i,j));
        sum = 0.0;
        for k = -1:1
            for l = -1:1
                sum = sum + double(rotOrigI(r+k,c+l));
            end
        end
        mean9X9Vals(i,j) = sum / 9.0
    end
end

% for i=1:length(centC)
%     c = floor(centC(i));
%     r = floor(centR(i));
%     mean4X4Vals(i) = mean([rotOrigI(r,c) rotOrigI(r,c+1) rotOrigI(r+1,c) rotOrigI(r+1,c+1)]);
% end
% for i=1:length(centC)
%     c = round(centC(i));
%     r = round(centR(i));
%     sum = 0.0;
%     for j = -1:1
%         for k = -1:1
%             sum = sum + double(rotOrigI(r+j,c+k));
%         end
%     end
%     mean9X9Vals(i) = sum / 9.0;
% end