% Calculate the mean and Standard Deviation of the images for each R, G, B components
% We first need to find out the noise effects of the digital scanner to our analysis
% We scan the same image 5 times and estimate the average R, G, B, values


clear;
% load the image
Fly41A = imread('fly41A.jpg', 'jpg');
Fly41B = imread('fly41B.jpg', 'jpg');
Fly41C = imread('fly41C.jpg', 'jpg');
Fly41D = imread('fly41D.jpg', 'jpg');
Fly41E = imread('fly41E.jpg', 'jpg');
Fly41F = imread('fly41F.jpg', 'jpg');
Fly41G = imread('fly41G.jpg', 'jpg');
Fly41H = imread('fly41H.jpg', 'jpg');
Fly41I = imread('fly41I.jpg', 'jpg');
Fly41J = imread('fly41J.jpg', 'jpg');

% Noise analysis


[ M, N ] = size( Fly41A(:, :, 1) );

% Statistical data of the images

% RGB mean values
AvgFly41A_R = mean2( Fly41A(:, :, 1) );
AvgFly41A_G = mean2( Fly41A(:, :, 2) );
AvgFly41A_B = mean2( Fly41A(:, :, 3) );

AvgFly41B_R = mean2( Fly41B(:, :, 1) );
AvgFly41B_G = mean2( Fly41B(:, :, 2) );
AvgFly41B_B = mean2( Fly41B(:, :, 3) );

AvgFly41C_R = mean2( Fly41C(:, :, 1) );
AvgFly41C_G = mean2( Fly41C(:, :, 2) );
AvgFly41C_B = mean2( Fly41C(:, :, 3) );

AvgFly41D_R = mean2( Fly41D(:, :, 1) );
AvgFly41D_G = mean2( Fly41D(:, :, 2) );
AvgFly41D_B = mean2( Fly41D(:, :, 3) );

AvgFly41E_R = mean2( Fly41E(:, :, 1) );
AvgFly41E_G = mean2( Fly41E(:, :, 2) );
AvgFly41E_B = mean2( Fly41E(:, :, 3) );

AvgFly41F_R = mean2( Fly41F(:, :, 1) );
AvgFly41F_G = mean2( Fly41F(:, :, 2) );
AvgFly41F_B = mean2( Fly41F(:, :, 3) );

AvgFly41G_R = mean2( Fly41G(:, :, 1) );
AvgFly41G_G = mean2( Fly41G(:, :, 2) );
AvgFly41G_B = mean2( Fly41G(:, :, 3) );

AvgFly41H_R = mean2( Fly41H(:, :, 1) );
AvgFly41H_G = mean2( Fly41H(:, :, 2) );
AvgFly41H_B = mean2( Fly41H(:, :, 3) );

AvgFly41I_R = mean2( Fly41I(:, :, 1) );
AvgFly41I_G = mean2( Fly41I(:, :, 2) );
AvgFly41I_B = mean2( Fly41I(:, :, 3) );

AvgFly41J_R = mean2( Fly41J(:, :, 1) );
AvgFly41J_G = mean2( Fly41J(:, :, 2) );
AvgFly41J_B = mean2( Fly41J(:, :, 3) );



AVG_R = [ AvgFly41A_R AvgFly41B_R AvgFly41C_R AvgFly41D_R AvgFly41E_R AvgFly41F_R AvgFly41G_R AvgFly41H_R AvgFly41I_R AvgFly41J_R]
AVG_G = [ AvgFly41A_G AvgFly41B_G AvgFly41C_G AvgFly41D_G AvgFly41E_G AvgFly41F_G AvgFly41G_G AvgFly41H_G AvgFly41I_G AvgFly41J_G]
AVG_B = [ AvgFly41A_B AvgFly41B_B AvgFly41C_B AvgFly41D_B AvgFly41E_B AvgFly41F_B AvgFly41G_B AvgFly41H_B AvgFly41I_B AvgFly41J_B]

figure(1);
subplot(3,1,1);
plot( AVG_R );
title(' Average values of Red ');
grid;
subplot(3,1, 2);
plot( AVG_G );
title(' Average values of Green ');
grid;
subplot(3, 1,3);
plot( AVG_B);
title(' Average values of Blue ');
grid;

% Percentage of deviation
temp1 = mean( AVG_R );
temp2 = mean( AVG_G );
temp3 = mean( AVG_B );

for counter = 1: 10
   
AvgDeviation_R(1, counter) = AVG_R(1, counter) / temp1 ;
AvgDeviation_G(1, counter) = AVG_G(1, counter) / temp2 ;
AvgDeviation_B(1, counter) = AVG_B(1, counter) / temp3 ;

end;

AvgDeviation_R
AvgDeviation_G
AvgDeviation_B

figure(2);
subplot( 3, 1, 1);
plot( AvgDeviation_R );
title(' Deviation of Red Average ');
grid;
subplot(3, 1, 2);
plot( AvgDeviation_G );
title(' Deviation of Green Average ');
grid;
subplot(3, 1,3);
plot( AvgDeviation_B );
title(' Deviation of Blue Average ');
grid;

% Standard Deviation of the noise
stdFly41A_R = std2( Fly41A(:, :, 1) );
stdFly41A_G = std2( Fly41A(:, :, 2) );
stdFly41A_B = std2( Fly41A(:, :, 3) );

stdFly41B_R = std2( Fly41B(:, :, 1) );
stdFly41B_G = std2( Fly41B(:, :, 2) );
stdFly41B_B = std2( Fly41B(:, :, 3) );

stdFly41C_R = std2( Fly41C(:, :, 1) );
stdFly41C_G = std2( Fly41C(:, :, 2) );
stdFly41C_B = std2( Fly41C(:, :, 3) );

stdFly41D_R = std2( Fly41D(:, :, 1) );
stdFly41D_G = std2( Fly41D(:, :, 2) );
stdFly41D_B = std2( Fly41D(:, :, 3) );

stdFly41E_R = std2( Fly41E(:, :, 1) );
stdFly41E_G = std2( Fly41E(:, :, 2) );
stdFly41E_B = std2( Fly41E(:, :, 3) );

stdFly41F_R = std2( Fly41F(:, :, 1) );
stdFly41F_G = std2( Fly41F(:, :, 2) );
stdFly41F_B = std2( Fly41F(:, :, 3) );

stdFly41G_R = std2( Fly41G(:, :, 1) );
stdFly41G_G = std2( Fly41G(:, :, 2) );
stdFly41G_B = std2( Fly41G(:, :, 3) );

stdFly41H_R = std2( Fly41H(:, :, 1) );
stdFly41H_G = std2( Fly41H(:, :, 2) );
stdFly41H_B = std2( Fly41H(:, :, 3) );

stdFly41I_R = std2( Fly41I(:, :, 1) );
stdFly41I_G = std2( Fly41I(:, :, 2) );
stdFly41I_B = std2( Fly41I(:, :, 3) );

stdFly41J_R = std2( Fly41J(:, :, 1) );
stdFly41J_G = std2( Fly41J(:, :, 2) );
stdFly41J_B = std2( Fly41J(:, :, 3) );


STD_R = [ stdFly41A_R stdFly41B_R stdFly41C_R stdFly41D_R stdFly41E_R stdFly41F_R stdFly41G_R stdFly41H_R stdFly41I_R stdFly41J_R]
STD_G = [ stdFly41A_G stdFly41B_G stdFly41C_G stdFly41D_G stdFly41E_G stdFly41F_G stdFly41G_G stdFly41H_G stdFly41I_G stdFly41J_G]
STD_B = [ stdFly41A_B stdFly41B_B stdFly41C_B stdFly41D_B stdFly41E_B stdFly41F_B stdFly41G_B stdFly41H_B stdFly41I_B stdFly41J_B]

temp4 = mean( STD_R );
temp5 = mean( STD_G );
temp6 = mean( STD_B );

for counter = 1: 10
StandDeviation_R(1, counter)  = STD_R(1, counter) / temp4 ;
StandDeviation_G(1, counter)  = STD_G(1, counter) / temp5 ;
StandDeviation_B(1, counter)  = STD_B(1, counter) / temp6 ;
end;


StandDeviation_R
StandDeviation_G
StandDeviation_B


figure(3);
subplot( 3,1,1);
plot( STD_R);
title(' Standard Deviation of Red ');
grid;
subplot(3, 1, 2);
plot( STD_G );
title(' Standard Deviation of Green ');
grid;
subplot(3, 1, 3);
plot( STD_B );
title(' Standard Deviation of Blue ');
grid;

figure(4);
subplot( 3,1,1);
plot( StandDeviation_R);
title(' Percentage of Deviation for Standard Deviation of Red ');
grid;
subplot(3, 1, 2);
plot( StandDeviation_G );
title(' Percentage of Deviation for Standard Deviation of Green ');
grid;
subplot(3, 1, 3);
plot( StandDeviation_B );
title(' Percentage of Deviation for Standard Deviation of Blue ');
grid;























