clear all;
band1 = imread ( 'ori\band1.gif' );
band2 = imread ( 'ori\band2.gif' );
band3 = imread ( 'ori\band3.gif' );
band4 = imread ( 'ori\band4.gif' );

min1 = min(min(band1));
min2 = min(min(band2));
min3 = min(min(band3));
min4 = min(min(band4));

max1 = max(max(band1));
max2 = max(max(band2));
max3 = max(max(band3));
max4 = max(max(band4));

img1 = uint8 ( ( double(band1) - double(min1) ) * 255 / ( double(max1) - double(min1) ) );
img2 = uint8 ( ( double(band2) - double(min2) ) * 255 / ( double(max2) - double(min2) ) );
img3 = uint8 ( ( double(band3) - double(min3) ) * 255 / ( double(max3) - double(min3) ) );
img4 = uint8 ( ( double(band4) - double(min4) ) * 240 / ( double(max4) - double(min4) ) );

% imwrite ( img1 , 'balanced\balband1.gif', 'bmp' );
% imwrite ( img2 , 'balanced\balband2.gif', 'bmp' );
% imwrite ( img3 , 'balanced\balband3.gif', 'bmp' );
% imwrite ( img4 , 'balanced\balband4.gif', 'bmp' );
% 
% new = zeros ( 512, 512, 3 );
% new(:,:,1) = img1;
% new(:,:,2) = img2;
% new(:,:,3) = img3;
% new = uint8 ( new );
% imwrite ( new , 'balanced\color\colored1.gif', 'bmp' );
% 
% new(:,:,1) = img4;
% imwrite ( new , 'balanced\color\colored2.gif', 'bmp' );
% 
% new(:,:,1) = img1;
% new(:,:,2) = img4;
% imwrite ( new , 'balanced\color\colored3.gif', 'bmp' );
% 
% new(:,:,2) = img2;
% new(:,:,3) = img4;
% imwrite ( new , 'balanced\color\colored4.gif', 'bmp' );

clear band1; clear band2; clear band3; clear band4;
clear min1; clear min2; clear min3; clear min4;
clear max1; clear max2; clear max3; clear max4;
% clear new;

l1 = imread ( 'training\l1.gif' );
land = zeros ( 300, 4 );
w1 = imread ( 'training\w4.gif' );
water = zeros ( 300, 4 );
lc = 1;
wc = 1;
for x = 1 : 512
    for y = 1 : 512
        if ( l1(x,y) == 255 )
            land( lc, 1 ) = img1(x,y);
            land( lc, 2 ) = img2(x,y);
            land( lc, 3 ) = img3(x,y);
            land( lc, 4 ) = img4(x,y);
            lc = lc + 1;
        end
        if ( w1(x,y) == 255 )
            water( wc, 1 ) = img1(x,y);
            water( wc, 2 ) = img2(x,y);
            water( wc, 3 ) = img3(x,y);
            water( wc, 4 ) = img4(x,y);
            wc = wc + 1;
        end
    end
end
clear l1; clear w1;
clear x; clear y;
lc
wc
clear lc; clear wc;

lmean = mean ( land );
wmean = mean ( water );
lcov = cov ( land );
wcov = cov ( water );
lcovInv = inv ( lcov );
wcovInv = inv ( wcov );
dlcov = det ( lcov );
dwcov = det ( wcov );
sqrtdlcov = sqrt ( dlcov );
sqrtdwcov = sqrt ( dwcov );

v  = zeros ( 4 );
p1 = zeros ( 4 );
p2 = zeros ( 4 );
oimg1 = zeros ( 512, 512 );
oimg2 = zeros ( 512, 512 );
oimg3 = zeros ( 512, 512 );

p11 = 0.5;
p12 = 1 - p11;
p21 = 0.3;
p22 = 1 - p21;
p31 = 0.7;
p32 = 1 - p31;

for x = 1 : 512
    for y = 1 : 512
        v = double ( [ img1(x,y) img2(x,y) img3(x,y) img4(x,y) ] );
        p1 = v - wmean;
        t1 = exp ( -0.5 * p1 * wcovInv * p1' ) / sqrtdwcov;
        p2 = v - lmean;
        t2 = exp ( -0.5 * p2 * lcovInv * p2' ) / sqrtdlcov;
        t11 = p11 * t1 ;
        t12 = p12 * t2 ;
        if ( t11 > t12 )
            oimg1(x,y) = 1;
        else
            oimg1(x,y) = 256;
        end
        t21 = p21 * t1 ;
        t22 = p22 * t2;
        if ( t21 > t22 )
            oimg2(x,y) = 1;
        else
            oimg2(x,y) = 255;
        end
        t31 = p31 * t1 ;
        t32 = p32 * t2;
        if ( t31 > t32 )
            oimg3(x,y) = 1;
        else
            oimg3(x,y) = 256;
        end
    end
end

oimg1 = uint8 ( oimg1 );
oimg2 = uint8 ( oimg2 );
oimg3 = uint8 ( oimg3 );

imwrite ( oimg1, 'out\bayes1.gif', 'bmp' );
imwrite ( oimg2, 'out\bayes2.gif', 'bmp' );
imwrite ( oimg3, 'out\bayes3.gif', 'bmp' );

imshow ( oimg2 );