% FIRST STEP: read each image and cut them

% http://rapidfire.sci.gsfc.nasa.gov/subsets/?subset=Japan 
% pixel size 250m
% before tsunami feb 23 2011, band 721
% after tsunami mar 13 2011, band 721

clf; clear; close all
I1=imread('Japan.2011054.aqua.721.250m.tif');
I2=imread('Japan.2011072.aqua.721.250m.tif');
I3=imread('Japan.2011054.aqua.250m.tif');
I4=imread('Japan.2011072.aqua.250m.tif');

%Show original image
figure, imshow(I1)
figure, imshow(I2)
figure, imshow(I3)
figure, imshow(I4)

%cut image. 
%imcrop retrieve coordinate (copy postion M1)
M1 = [4278.5 3100.5 636 1716]
G1 = imcrop(I1,M1) 
G2 = imcrop(I2,M1)
G3 = imcrop(I3,M1) 
G4 = imcrop(I4,M1)

% show cut images
figure, imshow(G1) 
figure, imshow(G2)
figure, imshow(G3) 
figure, imshow(G4)

%Convert it to black & white:
level1 = graythresh(G1); 
level2 = graythresh(G2); 
bw1 = im2bw(G1,level1);
bw2 = im2bw(G2,level2);

%Remove any noise. Pick a good minimum size:
bw1 = bwareaopen(bw1,50000);
bw2 = bwareaopen(bw2,50000);

%SECOND STEP: Inland, Invert image to get a single grain and 
% Invert image again to remove noise
bw1=1-bw1;
bw2=1-bw2;

%Identify all white grains with the connected components 
cc1 = bwconncomp(bw1, 8)
cc2 = bwconncomp(bw2, 8)

% Remove noise inverted image to get a single grain
sbw1 = bwareaopen(bw1, 4000);
sbw2 = bwareaopen(bw2, 4000);


%THIRD STEP: Inland, to get a single grain properties to comparate areas
%Identify all white grains with the connected components function:
cc1 = bwconncomp(sbw1, 8)
cc2 = bwconncomp(sbw2, 8)

% Display grain number 
idx = 1 
grain1 = false(size(sbw1)); 
grain2 = false(size(sbw2)); 
grain1(cc1.PixelIdxList{idx}) = true;
grain2(cc2.PixelIdxList{idx}) = true;

% cut image again to focus un the area tha we want to calculate
M2 = [44.5 127.5 477 1241];
bw1 = imcrop(sbw1,M2);
bw2 = imcrop(sbw2,M2);

%label grains
bw1=bwlabel(bw1,8)
bw2=bwlabel(bw2,8)

figure, imshow(bw1)
%impixelinfo
figure, imshow(bw2)
%impixelinfo

% grain properties
graindata1=regionprops(cc1,'all');
graindata2=regionprops(cc2,'all');
grain_area1=[graindata1.Area];
grain_area2=[graindata2.Area];
Area_inland=grain_area2-grain_area1;



% FOURTH STEP: Single grain dilatation
SE = strel('disk', 150, 8);
dsg1= imdilate(1-sbw1,SE);
dsg2= imdilate(1-sbw2,SE);
figure, imshow(dsg1)
figure, imshow(dsg2)
D=dsg1-(1-sbw1); %dilatation area
figure, imshow(D)


%FIFTH: Ocean Floor: extract information for ocean floor on dilated grain Mask
D2 = logical(D);
D3 = cat(3, D2, D2, D2);
A3 = immultiply(G3,D3);
A4 = immultiply(G4,D3);
figure;imshow(A3);
figure;imshow(A4);

% SIXTH: Ocean Floor: Aply Mahalanobis distance on the mask in both images


filtSize = 25;

if (exist('myPatches.mat'))
    load('myPatches.mat');
else
    % call read patch code
roiPixA = getPixStd(A3, filtSize);
roiPixB = getPixStd(A4, filtSize);

roiPix = cat(1, roiPixA, roiPixB);
myStd = max(std(roiPix));

save('myPatches', 'roiPix', 'myStd');
end

% A3 & A4
% read and check image, get size
    [nRows nCols nLyrs] = size(A3);
    nPix = nRows * nCols;
    
    % smooth image A3
    fprintf('Smoothing image ...\n');
    myFilt = ones(filtSize) ./ filtSize^2;
    smoothImg1 = imfilter(A3, myFilt);
    imgPix1= single(reshape(A3, nPix, 3));
    
    % smooth image A4
    fprintf('Smoothing image ...\n');
    myFilt = ones(filtSize) ./ filtSize^2;
    smoothImg2 = imfilter(A4, myFilt);
    imgPix2= single(reshape(A4, nPix, 3));
    
    % compute distance vector A3
    mahDist1 = mahal(imgPix1, roiPix);
    
    % compute distance vector A4
    mahDist2 = mahal(imgPix2, roiPix);
    
    
    % reshape distance and store in matrix
    myDist1 = reshape(mahDist1, nRows, nCols);
    
    % reshape distance and store in matrix
    myDist2 = reshape(mahDist2, nRows, nCols);
    
    %display the distance
    figure;
    imagesc(myDist1);
    title('Mahalanobis distance from region color');
    %colorbar;
    %impixelinfo
    
    % define how many standard deviations are allowed
    kThresh = 1.2;
    
    % myDist contains the distances of each pixel to mean color
    % if minimum is more than k*std label it 0, else label 1
    fprintf('Thresholding pixels with distance from region color ...\n');
    mySegm1 = false(nRows, nCols);
    mySegm1(myDist1 <= kThresh * myStd) = true;
    
    fprintf('Thresholding pixels with distance from region color ...\n');
    mySegm2 = false(nRows, nCols);
    mySegm2(myDist2 <= kThresh * myStd) = true;

    % get the red, green and blue parts of the image A3
    Ir1 = A3(:,:,1);
    Ig1 = A3(:,:,2);
    Ib1 = A3(:,:,3);

    % get the red, green and blue parts of the image A4
    Ir2 = A4(:,:,1);
    Ig2 = A4(:,:,2);
    Ib2 = A4(:,:,3);

    % label the segmented pixels in red A3
    Ir1(mySegm1==1) = 255;
    Ib1(mySegm1==1) = 0;
    Ig1(mySegm1==1) = 0;
    
% label the segmented pixels in red A4
    Ir2(mySegm2==1) = 255;
    Ib2(mySegm2==1) = 0;
    Ig2(mySegm2==1) = 0;

    % put the red, green, and blue parts back together
    I1 = cat(3,Ir1,Ig1,Ib1);
    I2 = cat(3,Ir2,Ig2,Ib2);
    figure; imshow(I1);
    figure; imshow(I2);
    

% SEVENTH: Calculate area
mycount1=sum(mySegm1(:));
mycount2=sum(mySegm2(:));
Area_offshore=mycount2-mycount1;


%FINAL RESULTS image pixel size 250 m

Area_offshore=mycount2-mycount1

Area_inland=grain_area2-grain_area1


   