%% MATLAB Recipes for Earth Sciences - Chapter 8

%% Section 8.3
clear

%%
unconform1 = imread('unconform.jpg');

%%
whos

%%
unconform1(50:55,50:55,1)

%%
imshow(unconform1)

%%
unconform2 = rgb2gray (unconform1);

%%
imshow(unconform2)

%%
imhist(unconform2)

%%
unconform3 = histeq(unconform2);

%%
imshow(unconform3)

%%
imwrite(unconform3,'unconform3.jpg')

%%
imfinfo('unconform3.jpg')

%%
[x,map] = rgb2ind(unconform1,16);
imshow(unconform1), figure, imshow(x,map)

%% Section 8.4
clear

%%
I1 = hdfread('naivasha.hdf','VNIR_Band3N','Fields','ImageData');
I2 = hdfread('naivasha.hdf','VNIR_Band2','Fields','ImageData');
I3 = hdfread('naivasha.hdf','VNIR_Band1','Fields','ImageData');

%%
naivasha_rgb = cat(3,I1,I2,I3);

%%
imshow(naivasha_rgb);

%%
imwrite(naivasha_rgb,'naivasha.tif','tif')

%% Section 8.5
clear

%%
I1 = hdfread('naivasha.hdf','VNIR_Band3N','Fields','ImageData');
I2 = hdfread('naivasha.hdf','VNIR_Band2','Fields','ImageData');
I3 = hdfread('naivasha.hdf','VNIR_Band1','Fields','ImageData');

%%
naivasha_rgb = cat(3,I1,I2,I3);

%%
upperleft = [-0.319922, 36.214332];
upperright = [-0.400443, 36.770406];
lowerleft = [-0.878267, 36.096003];
lowerright = [-0.958743, 36.652213];

%%
inputpoints(1,:) = upperleft;
inputpoints(2,:) = lowerleft;
inputpoints(3,:) = upperright;
inputpoints(4,:) = lowerright;
inputpoints = fliplr(inputpoints);

%%
basepoints(1,:) = [1,4200];
basepoints(2,:) = [1,1];
basepoints(3,:) = [4100,4200];
basepoints(4,:) = [4100,1];

%%
tform = cp2tform(inputpoints,basepoints,'affine');

%%
[newnaivasha_rgb,x,y] = imtransform(naivasha_rgb,tform);

%%
X = 36.096003 : (36.770406-36.096003)/8569 : 36.770406;
Y = 0.319922 : ( 0.958743- 0.319922)/8400 : 0.958743;

%%
imshow(newnaivasha_rgb,'XData',X,'YData',Y), axis on, grid on
xlabel('Longitude'), ylabel('Latitude')
title('Georeferenced ASTER Image')

%%
print -djpeg70 -r600 naivasha_georef.jpg

%% Section 8.6
clear

%%
data = minput('naivasha_georef.jpg')

%% Section 8.7
clear

%%
I = imread('varves.tif');

%%
imshow(I), axis on

%%
[x,y] = ginput;

%%
ix = 5 * size(I,2) / sqrt((y(2,1)-y(1,1))^2+(x(2,1)-x(1,1))^2);
iy = 5 * size(I,1) / sqrt((y(2,1)-y(1,1))^2+(x(2,1)-x(1,1))^2);

%%
imshow(I,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')

%%
[CX,CY,C] = improfile;

%%
imshow(I,'XData',[0 ix],'YData',[0 iy]), hold on
plot(CX,CY), hold off

%%
figure
plot(CY,C(:,1),'r',CY,C(:,2),'g',CY,C(:,3),'b')
xlabel('Centimeters'), ylabel('Intensity')

%%
imshow(I,'XData',[0 ix],'YData',[0 iy]), hold on
plot(CX,CY), hold off
laminae = ginput;

%%
imshow(I,'XData',[0 ix],'YData',[0 iy])
hold on
plot(CX,CY)
plot(laminae(:,1),laminae(:,2),'ro')
xlabel('Centimeters'), ylabel('Centimeters')
hold off

%%
newlaminae(:,1) = laminae(:,2);
newlaminae(:,2) = 1 : length(laminae(:,2));
age = interp1(newlaminae(:,1),newlaminae(:,2),CY);

%%
subplot(2,1,1), plot(CY,C(:,1),CY,C(:,2),CY,C(:,3))
xlabel('Pixel ID'), ylabel('Intensity'), title('Color vs. Length')
subplot(2,1,2), plot(age,C(:,1),age,C(:,2),age,C(:,3))
xlabel('Years'), ylabel('Intensity'), title('Color vs. Age')

%% Section 8.8´
clear

%%
I1 = imread('grainsize.tif');
ix = 3; iy = 284 * 3 / 367;
imshow(I1,'XData',[0 ix],'YData',[0 iy])
title('Original Image')

%%
I2 = rgb2gray(I1);
imshow(I2,'XData',[0 ix],'YData',[0 iy])
title('Grayscale Image')

%%
I3 = imadjust(I2);
imshow(I3,'XData',[0 ix],'YData',[0 iy])
title('Adjusted Intensity Values')

%%
I4 = imopen(I3,strel('disk',1));
imshow(I4,'XData',[0 ix],'YData',[0 iy])
title('No Background')

%%
I5 = imsubtract(I3,I4);
imshow(I5,'XData',[0 ix],'YData',[0 iy])
title('Background')

%%
I6 = im2bw(I4,0.2);
imshow(I6,'XData',[0 ix],'YData',[0 iy])
title('Binary Image')

%%
I7 = imclearborder(I6);
himage1 = imshow(I6,'XData',[0 ix],'YData',[0 iy]); hold on
set(himage1, 'AlphaData', 0.7);
himage2 = imshow(imsubtract(I6,I7),'XData',[0 ix],'YData',[0 iy]);
set(himage2, 'AlphaData', 0.4);
title('Image Border'), hold off

%%
[B,L] = bwboundaries(I7,'noholes');
imshow(label2rgb(L,@jet,'w','shuffle'),'XData',[0 ix],'YData',[0 iy])
title('Define Objects')

%%
[labeled,numObjects] = bwlabeln(I7,8);
numObjects

%%
graindata = regionprops(labeled,'basic');
grainareas= [graindata(:).Area];
objectareas = 3^2 * grainareas * 367^(-2);

%%
max_area = max(objectareas)
min_area = min(objectareas)
mean_area = mean(objectareas)

%%
clf
v = 0.0005 : 0.0005 : 0.15;
hist(objectareas,v)
xlabel('Grain Size in Millimeters^2')
ylabel('Number of Grains')
axis([0 0.1 0 30])

%%
D = bwdist(~I7,'cityblock');

%%
D = -D;
D(~I7) = -Inf;
 
%%
L2 = watershed(D); 
imshow(label2rgb(L2,@jet,'w','shuffle'),'XData',[0 ix],'YData',[0 iy])
title('Watershed Segmentation')

%%
objects = sortrows(L2(:),1);
max(objects)
clear objectsizes
for i = 2 : max(objects)
 clear individualobject
 individualobject = objects(objects == i);
 objectsizes(i) = length(individualobject);
end
objectsizes = objectsizes';
objectsizes = sortrows(objectsizes,1);
objectsizes = objectsizes(objectsizes~=0);

%%
objectareas = 3^2 * objectsizes * 367^(-2);

%%
max_area = max(objectareas)
min_area = min(objectareas)
mean_area = mean(objectareas)

%%
clf
v = 0.0005 : 0.0005 : 0.15;
hist(objectareas,v)
xlabel('Grain Size in Millimeters^2'),ylabel('Number of Grains')
axis([0 0.1 0 70])

%%
figure
imshow(I1,'XData',[0 ix],'YData',[0 iy])
data = ginput;

%%
data(end+1,:) = data(1,:)

%%
imshow(I1,'XData',[0 ix],'YData',[0 iy]), hold on
plot(data(:,1),data(:,2)), hold off

%%
polyarea(data(:,1),data(:,2))

%% Section 8.9
clear

%%
I1 = imread('lakesediment.jpg');
ix = 1; iy = 1;
imshow(I1,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeter'), ylabel('Centimeter')
title('Original Image')

%%
I2 = rgb2gray(I1);
imshow(I2,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('Grayscale')

%%
I3 = imadjust(I2);
imshow(I3,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('Better Contrast')

%%
I4 = imopen(I3,strel('disk',5));
imshow(I4,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('W/O Background')

%%
I5 = imsubtract(I3,I4);
imshow(I5,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('Background')

%%
I6 = im2bw(I4,0.03);
imshow(I6,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('Only Charcoal')

%%
I7 = im2bw(I4,0.6);
imshow(I7,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('All Objects')

%%
100*sum(sum(I6==0))/sum(sum(I7==0))
