clear
close all

%I1 = imread('pcr11sep1323h2000ml40s.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 B/W Image');
level = graythresh(maskedI);
bwI = im2bw(maskedI, 0.6*level);
figure
imshow(bwI)
title('B/W Image');
bwOpenedI = bwareaopen(bwI, 32);
figure
imshow(bwOpenedI);
title('B/W Opened Image')

cc = bwconncomp(bwOpenedI);
imgStats = regionprops(cc, bwOpenedI, 'Area', 'Orientation');
a = [imgStats(:).Area];
o = [imgStats(:).Orientation];
regAreaIndex = find(a == max(a));
regAngle = imgStats(regAreaIndex).Orientation;
rotBWOpenedI = imrotate(bwOpenedI, -1.0*(regAngle+88.5));
figure
imshow(rotBWOpenedI)
title('Rotated B/W Opened Image');
rotOrigI = imrotate(origI, -1.0*(regAngle+88.5));
figure
imshow(rotOrigI)
title('Rotated Original Image');

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];

centRC = [centR; centC]';
sort1CentRC = sortrows(centRC,1);
for i=1:10
    tempSort = sortrows(sort1CentRC((i*8)-7:(i*8)+1,:),2);
    for j=1:8
        sortCentArray{i,j} = [tempSort(j,:)];
    end    
end

%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