clc;
clear;



%處理農田樣本
crI=imread('CR.jpg');
crr=crI(:,:,1);crg=crI(:,:,2);crb=crI(:,:,3);%提取圖像三個波段
crrr=im2double(crr);crgg=im2double(crg);crbb=im2double(crb);%將黑白圖像轉換成雙精度格式并灰度值歸一化為[0,1]之間
crro=[];
crgo=[];
crbo=[];
[crlen crlen2]=size(crrr);%計算rr矩陣的寬度和長度
crx=[];

%將矩陣轉換成一維的矩陣
for i=1:crlen
    crro=[crro,crrr(i,:)];   %cat
end
for i=1:crlen
    crgo=[crgo,crgg(i,:)];
end
for i=1:crlen
    crbo=[crbo,crbb(i,:)];
end

crx=[crro;crgo;crbo];    %將三個波段的矩陣合并為一個矩陣


crlength=length(crx);

%處理海水樣本
seaI=imread('sea.jpg');
sear=seaI(:,:,1);seag=seaI(:,:,2);seab=seaI(:,:,3);%提取圖像三個波段
searr=im2double(sear);seagg=im2double(seag);seabb=im2double(seab);%將黑白圖像轉換成雙精度格式并灰度值歸一化為[0,1]之間
searo=[];
seago=[];
seabo=[];
[sealen sealen2]=size(searr);%計算rr矩陣的寬度和長度
seax=[];
seax1=[];
%將300*300的矩陣轉換成一維的矩陣
for i=1:sealen
    searo=[searo,searr(i,:)];   %cat
end
for i=1:sealen
    seago=[seago,seagg(i,:)];
end
for i=1:sealen
    seabo=[seabo,seabb(i,:)];
end

seax=[searo;seago;seabo];    %將三個波段的矩陣合并為一個矩陣


sealength=length(seax);

%處理蘆葦樣本
reedI=imread('reed.jpg');
reedr=reedI(:,:,1);reedg=reedI(:,:,2);reedb=reedI(:,:,3);%提取圖像三個波段
reedrr=im2double(reedr);reedgg=im2double(reedg);reedbb=im2double(reedb);%將黑白圖像轉換成雙精度格式并灰度值歸一化為[0,1]之間
reedro=[];
reedgo=[];
reedbo=[];
[reedlen reedlen2]=size(reedrr);%計算rr矩陣的寬度和長度
reedx=[];
reedx1=[];
%將矩陣轉換成一維的矩陣
for i=1:reedlen
    reedro=[reedro,reedrr(i,:)];   %cat
end
for i=1:reedlen
    reedgo=[reedgo,reedgg(i,:)];
end
for i=1:reedlen
    reedbo=[reedbo,reedbb(i,:)];
end

reedx=[reedro;reedgo;reedbo];    %將三個波段的矩陣合并為一個矩陣


reedlength=length(reedx);

%處理綠地樣本
grI=imread('GR.jpg');
grr=grI(:,:,1);grg=grI(:,:,2);grb=grI(:,:,3);%提取圖像三個波段
grrr=im2double(grr);grgg=im2double(grg);grbb=im2double(grb);%將黑白圖像轉換成雙精度格式并灰度值歸一化為[0,1]之間
grro=[];
grgo=[];
grbo=[];
[grlen grlen2]=size(grrr);%計算rr矩陣的寬度和長度
grx=[];
grx1=[];
%將矩陣轉換成一維的矩陣
for i=1:grlen
    grro=[grro,grrr(i,:)];   %cat
end
for i=1:grlen
    grgo=[grgo,grgg(i,:)];
end
for i=1:grlen
    grbo=[grbo,grbb(i,:)];
end

grx=[grro;grgo;grbo];    %將三個波段的矩陣合并為一個矩陣


grlength=length(grx);

%處理養殖場樣本
culI=imread('cul.jpg');
culr=culI(:,:,1);culg=culI(:,:,2);culb=culI(:,:,3);%提取圖像三個波段
culrr=im2double(culr);culgg=im2double(culg);culbb=im2double(culb);%將黑白圖像轉換成雙精度格式并灰度值歸一化為[0,1]之間
culro=[];
culgo=[];
culbo=[];
[cullen cullen2]=size(culrr);%計算rr矩陣的寬度和長度
culx=[];
culx1=[];
%將矩陣轉換成一維的矩陣
for i=1:cullen
    culro=[culro,culrr(i,:)];   %cat
end
for i=1:cullen
    culgo=[culgo,culgg(i,:)];
end
for i=1:cullen
    culbo=[culbo,culbb(i,:)];
end

culx=[culro;culgo;culbo];    %將三個波段的矩陣合并為一個矩陣


cullength=length(culx);

%處理灘涂樣本
mudI=imread('mudflat.jpg');
mudr=mudI(:,:,1);mudg=mudI(:,:,2);mudb=mudI(:,:,3);%提取圖像三個波段
mudrr=im2double(mudr);mudgg=im2double(mudg);mudbb=im2double(mudb);%將黑白圖像轉換成雙精度格式并灰度值歸一化為[0,1]之間
mudro=[];
mudgo=[];
mudbo=[];
[mudlen mudlen2]=size(mudrr);%計算rr矩陣的寬度和長度
mudx=[];
mudx1=[];
%將矩陣轉換成一維的矩陣
for i=1:mudlen
    mudro=[mudro,mudrr(i,:)];   %cat
end
for i=1:mudlen
    mudgo=[mudgo,mudgg(i,:)];
end
for i=1:mudlen
    mudbo=[mudbo,mudbb(i,:)];
end

mudx=[mudro;mudgo;mudbo];    %將三個波段的矩陣合并為一個矩陣


mudlength=length(mudx);

%處理互米花草樣本
huI=imread('humi.jpg');
hur=huI(:,:,1);hug=huI(:,:,2);hub=huI(:,:,3);%提取圖像三個波段
hurr=im2double(hur);hugg=im2double(hug);hubb=im2double(hub);%將黑白圖像轉換成雙精度格式并灰度值歸一化為[0,1]之間
huro=[];
hugo=[];
hubo=[];
[hulen hulen2]=size(hurr);%計算rr矩陣的寬度和長度
hulength=hulen*hulen2
hux=[];
hux1=[];
%將矩陣轉換成一維的矩陣
for i=1:hulen
    huro=[huro,hurr(i,:)];   %cat
end
for i=1:hulen
    hugo=[hugo,hugg(i,:)];
end
for i=1:hulen
    hubo=[hubo,hubb(i,:)];
end

hux=[huro;hugo;hubo];    %將三個波段的矩陣合并為一個矩陣


hulength=length(hux);

%處理房屋樣本
fwI=imread('fw.jpg');
fwr=fwI(:,:,1);fwg=fwI(:,:,2);fwb=fwI(:,:,3);%提取圖像三個波段
fwrr=im2double(fwr);fwgg=im2double(fwg);fwbb=im2double(fwb);%將黑白圖像轉換成雙精度格式并灰度值歸一化為[0,1]之間
fwro=[];
fwgo=[];
fwbo=[];
[fwlen fwlen2]=size(fwrr);%計算rr矩陣的寬度和長度
fwx=[];
fwx1=[];
%將矩陣轉換成一維的矩陣
for i=1:fwlen
    fwro=[fwro,fwrr(i,:)];   %cat
end
for i=1:fwlen
    fwgo=[fwgo,fwgg(i,:)];
end
for i=1:fwlen
    fwbo=[fwbo,fwbb(i,:)];
end
fwx=[fwro;fwgo;fwbo];    %將三個波段的矩陣合并為一個矩陣


fwlength=length(fwx);

%處理海三棱藨草的樣本
sanI=imread('san.jpg');
sanr=sanI(:,:,1);sang=sanI(:,:,2);sanb=sanI(:,:,3);%提取圖像三個波段
sanrr=im2double(sanr);sangg=im2double(sang);sanbb=im2double(sanb);%將黑白圖像轉換成雙精度格式并灰度值歸一化為[0,1]之間
sanro=[];
sango=[];
sanbo=[];
[sanlen sanlen2]=size(sanrr);%計算rr矩陣的寬度和長度
sanx=[];
sanx1=[];
%將矩陣轉換成一維的矩陣
for i=1:sanlen
    sanro=[sanro,sanrr(i,:)];   %cat
end
for i=1:sanlen
    sango=[sango,sangg(i,:)];
end
for i=1:sanlen
    sanbo=[sanbo,sanbb(i,:)];
end

sanx=[sanro;sango;sanbo];    %將三個波段的矩陣合并為一個矩陣


sanlength=length(sanx);

%處理未利用土地的樣本
unI=imread('un.jpg');
unr=unI(:,:,1);ung=unI(:,:,2);unb=unI(:,:,3);%提取圖像三個波段
unrr=im2double(unr);ungg=im2double(ung);unbb=im2double(unb);%將黑白圖像轉換成雙精度格式并灰度值歸一化為[0,1]之間
unro=[];
ungo=[];
unbo=[];
[unlen unlen2]=size(unrr);%計算rr矩陣的寬度和長度
unx=[];
unx1=[];
%將矩陣轉換成一維的矩陣
for i=1:unlen
    unro=[unro,unrr(i,:)];   %cat
end
for i=1:unlen
    ungo=[ungo,ungg(i,:)];
end
for i=1:unlen
    unbo=[unbo,unbb(i,:)];
end

unx=[unro;ungo;unbo];   %將三個波段的矩陣合并為一個矩陣

unlength=length(unx);

%將樣本置于一個一維向量中
%
% P_G=[crgo,seago,reedgo,grgo,culgo,mudgo,hugo,fwgo,sango,ungo];
% P_B=[crbo,seabo,reedbo,grbo,culbo,mudbo,hubo,fwbo,sanbo,unbo];
% P=[P_RP_G,P_B];
 P=[crx,seax,reedx,grx,culx,mudx,hux,fwx,sanx,unx];



%設定目標向量，一維
T=[];
for i=1:100
    T(1,i)=0.1;
    T(2,i)=0.9;
    T(3,i)=0.1;
end
for i=101:200
    T(1,i)=0.2;
    T(2,i)=0.2;
    T(3,i)=0.9;
end
for i=201:300
    T(1,i)=0.9;
    T(2,i)=0.5;
    T(3,i)=0.2;
end
for i=301:400
    T(1,i)=0.6;
    T(2,i)=0.5;
    T(3,i)=0.7;
end
for i=401:464
    T(1,i)=0.3;
    T(2,i)=0.1;
    T(3,i)=0.1;
end
for i=465:564
    T(1,i)=0.8;
    T(2,i)=0.5;
    T(3,i)=0.1;
end
for i=565:664
    T(1,i)=0.4;
    T(2,i)=0.1;
    T(3,i)=0.7;
end
for i=665:764
    T(1,i)=0.5;
    T(2,i)=0.9;
    T(3,i)=0.9;
end
for i=765:864
    T(1,i)=0.25;
    T(2,i)=0.6;
    T(3,i)=0.9;
end
for i=865:964
    T(1,i)=0.7;
    T(2,i)=0.5;
    T(3,i)=0.1;
end


clc;
net=newcf(minmax(P),[15,3],{'tansig','purelin'},'traingdx');
%當前輸入層權值和閾值
imputweights=net.iw{1,1}
inputbias=net.b{1}  %產生隱含層各神經元的閾值
%當前網路層權值和閾值 
layerweights=net.lw{2,1}
layerbias=net.b{2}

%設置訓練參數
net.trainparam.show=50 %迭代50次顯示一次
net.trainparam.lr=0.05 %學習率為0.05（一般大小在0-1之間）
net.trainparam.mc=0.9;  %動量因子0.9（0-1之間）
net.trainparam.epochs=1000;%迭代次數1000次
net.trainparam.goal=1e-3;%誤差
clc
%調用traingdm算法訓練bp網絡
newnet=train(net,P,T);
clc
%仿真
I=imread('tm2000mask.jpg');

ro=I(:,:,1);
go=I(:,:,2); 
bo=I(:,:,3); 
rr=im2double(ro);%將r,g,b圖像變為雙精度格式，其灰度值歸一到0，1之間
gg=im2double(go);
bb=im2double(bo);
[m,n]=size(rr)
%將圖像變為一維向量形式
x=[];
for i=1:m
    x=[x,rr(i,:)];
end

y=[];
for i=1:m
    y=[y,gg(i,:)];
end

z=[];
for i=1:m
    z=[z,bb(i,:)];
end
tm=[x;y;z];
a=sim(newnet,tm);


rR=[];
R=[a(1,:)];
for i=0:779
   rR=[rR;R((i*1090+1):(i*1090)+1090)];
end



rG=[];
G=[a(2,:)];
for i=0:779
   rG=[rG;G((i*1090+1):(i*1090)+1090)];
end



rB=[];
B=[a(3,:)];
for i=0:779
   rB=[rB;B((i*1090+1):(i*1090)+1090)];
end
cr=0;
sea=0;
reed=0;
gr=0;
cul=0;
mud=0;
hu=0;
fw=0;
san=0;
un=0;
for i=1:780
 for j=1:1090
    if ((0.098<rR(i,j))&(rR(i,j)<0.12))
        cr=cr+1;
    else if ((0.19<rR(i,j))&(rR(i,j)<0.22))
            sea=sea+1;
        else if ((0.88<rR(i,j))&(rR(i,j)<0.92))
            reed=reed+1;
            else if ((0.58<rR(i,j))&(rR(i,j)<0.62))
                    gr=gr+1;
                else if ((0.28<rR(i,j))&(rR(i,j)<0.32))
                        cul=cul+1;
                    else if ((0.78<rR(i,j))&(rR(i,j)<0.82))
                            mud=mud+1;
                        else if ((0.38<rR(i,j))&(rR(i,j)<0.42))
                                hu=hu+1;
                            else if ((0.48<rR(i,j))&(rR(i,j)<0.52))
                                    fw=fw+1;
                                else if ((0.24<rR(i,j))&(rR(i,j)<0.26))
                                        san=san+1;
                                    else if ((0.68<rR(i,j))&(rR(i,j)<0.78))
                                            un=un+1;
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
 end
end
figure,pie([cr,sea,reed,gr,cul,mud,hu,fw,san,un]);
title('各類型土地利用比例圖');
legend('農田','海水','蘆葦','綠地','養殖場','灘涂','互米花草','房屋','海三棱藨草','未利用地');

%將數據轉成0-255格式
R1=abs(rR)*255;
G1=abs(rG)*255;
B1=abs(rB)*255;
%轉成彩圖RGB格式顯示
R2=uint8(R1);
G2=uint8(G1);
B2=uint8(B1);
RGB=cat(3,R2,G2,B2);
figure,imagesc(RGB);


disp('1--分類后紅色波段圖像');
disp('2--分類后綠色波段圖像');
disp('3--分類后藍色波段圖像');
choice=input('請輸入數字，對應要顯示的圖像：');
if choice==1
    figure,imshow(abs(rR))
else if choice==2
        figure,imshow(abs(rG))
    else if choice==3
            figure,imshow(abs(rB))
        end
    end
end