clear;close all;
fcol=['r','g','b']; %初始化色彩
for i =1:3
    T=imread(strcat('IMG',num2str(i),'.jpg'));%讀入圖像，不顯示
    %巴特沃茲濾波，不顯示濾波后圖像
    for k=1:3
        A=T(:,:,k);
        for j=1:300
            B=A(j,:);
            C=double(B); %將整形轉換為雙精度型
            [M,N]=butter(8,0.4); %使用8階，截至頻率是0.4的巴特沃茲濾波器
            D=filter(M,N,C);
            A(k,:)=D;
            T(:,:,k)=A;
       end
    end
   %在濾波后圖像上截取目標圖
   %imshow(T);
%    [x,y]=ginput(4);
%    a=min(x);
%    b=min(y);
%    width=max(x)-min(x);
%    height=max(y)-min(y);
%    I=imcrop(T,[a,b,width,height]);% %在目標圖上求灰度（強度）值
   rect=imcrop(T,[940,500,300,100]);
   %imshow(rect);
   %反應曲線
x=1:300; y=50*ones(1,length(x)); %y=ones(a,b);表示形成a行b列全為1的數列。
                                 %詞句表示y=100的一條長為550的線段。
Cur0=improfile(rect,x,y); %improfile用于沿著圖像中一條直線段路徑或直線路徑計算并繪制其強度（灰度）值。
                             %彩色圖像是將紅綠藍分離開了，個表達獨立的線條圖形。
%三點平均濾波

 Cur1=Cur0(:,:,2)';
 
    for p=2:(length(Cur1)-1)
    Cur2(p)=(Cur1(p-1)+Cur1(p)+Cur1(p+1))/3;
    Cur2(1)=Cur1(1);
    Cur2(length(Cur1))=Cur1(length(Cur1));
end
plot(Cur2,'g');grid on;
hold on;
%將曲線值的轉至存成一個矩陣，方便后面計算
juzhen(i,:)=Cur2;
%尋找峰值，法一：最小值閾值法。
%求最大峰值
%定義存儲最大峰值X坐標的矩陣目前為空矩陣,求綠色圖像的峰值關系。
BRange=[];
k=1;
for t=1:length(Cur2)
   if abs(min(Cur2)-Cur2(t))<10;
       BRange(k)=t;
       k=k+1;
   else
       continue
   end
end
MAX=mean(Cur2(BRange(1):BRange(length(BRange))));
% disp('MAX：') %顯示MAX:這個字符串
% disp(MAX);    %顯示MAX的值
%求最小峰值
Sub=Cur2(1:(BRange(1)-100));
SRange=[];
k=1;
for t=1:length(Sub)
    if abs(min(Sub)-Sub(t))<10;
        SRange(k)=t;
        k=k+1;
    else
        continue
    end 
end
MIN=mean(Cur2(SRange(1):SRange(length(SRange))));
% disp('MIN:');
% disp(MIN);
q=MAX/MIN;
% disp('q:');
% disp(q);
%將q值形成一個矩陣，方便以后調用
Qzhen(i,:)=q;
disp('Qzhen：') %顯示Qzhen:這個字符串
disp(Qzhen);    %顯示Qzhen的值
end
%循環結束
%進行Qzhen的數據分析，方差分析
fangcha=var(Qzhen);
disp('fangcha：') %顯示fangcha:這個字符串
disp(fangcha);    %顯示fangcha的值
bzpiancha=std(Qzhen);
disp('bzpiancha：') %顯示bzpiancha:這個字符串
disp(bzpiancha);    %顯示bzpiancha的值
%將10條曲線取平均作圖
pingjun=zeros(1,299);
for i=1:299
    pingjun(:,i)=mean(juzhen(:,i));
end
figure,plot(pingjun,'b');grid on;
%求平均曲線的峰值
BRange=[];
k=1;
for t=1:length(pingjun)
   if abs(min(pingjun)-pingjun(t))<10;
       BRange(k)=t;
       k=k+1;
   else
       continue
   end
end
MAX=mean(pingjun(BRange(1):BRange(length(BRange))));
% disp('MAX：') %顯示MAX:這個字符串
% disp(MAX);    %顯示MAX的值
%求最小峰值
Sub=pingjun(1:(BRange(1)-100));
SRange=[];
k=1;
for t=1:length(Sub)
    if abs(min(Sub)-Sub(t))<10;
        SRange(k)=t;
        k=k+1;
    else
        continue
    end 
end
MIN=mean(pingjun(SRange(1):SRange(length(SRange))));
% disp('MIN:');
% disp(MIN);
pingjunQ=MAX/MIN;
disp('pingjunQ：') %顯示Qzhen:這個字符串
disp(pingjunQ);    %顯示Qzhen的值