image=imread('mri1.jpg');
I=rgb2gray(image);
figure,imshow(I),title('原始圖像')
I=double(I)/255;
[M,N]=size(I);
x1=65;y1=56;x2=72;y2=113;
seed1=I(x1,y1);             %將生長起始點(diǎn)灰度值存入seed中
seed2=I(x2,y2); 
Y=zeros(M,N);             %作一個全零與原圖像等大的圖像矩陣Y，作為輸出圖像矩陣
Z=zeros(M,N); 
Y(x1,y1)=I(x1,y1);          %將Y中與所取點(diǎn)相對應(yīng)位置的點(diǎn)設(shè)置為與原圖像相同的灰度
Z(x2,y2)=I(x2,y2);          %將Z中與所取點(diǎn)相對應(yīng)位置的點(diǎn)設(shè)置為與原圖像相同的灰度

sum1=seed1;               %儲存符合區(qū)域生長條件的點(diǎn)的灰度值的和
suit1=1;                   %儲存符合區(qū)域生長條件的點(diǎn)的個數(shù)
count1=1;                 %記錄每次判斷一點(diǎn)周圍八點(diǎn)符合條件的新點(diǎn)的數(shù)目
threshold1=0.05555;         %判斷域值
while count1>0
 s1=0;                    %記錄判斷一點(diǎn)周圍八點(diǎn)時，符合條件的新點(diǎn)的灰度值之和
 count1=0;
 for i1=1:M
   for j1=1:ceil(N/2)
     if Y(i1,j1)~=0
      if (i1-1)>0 & (i1+1)<(M+1) & (j1-1)>1 & (j1+1)<((ceil(N/2))+1)
 %判斷此點(diǎn)是否為圖像邊界上的點(diǎn)
       for u=-1:1                                %判斷點(diǎn)周圍八點(diǎn)是否合和域值條件
        for v=-1:1                               %u,v為偏移量
          if  Y(i1+u,j1+v)==0 & abs(I(i1+u,j1+v)-seed1)<=threshold1    
%判斷是否未存在于輸出矩陣Y，并且為符合域值條件的點(diǎn)
             Y(i1+u,j1+v)=I(i1+u,j1+v);            %符合以上兩條件即將其在Y中與之位置對應(yīng)的點(diǎn)設(shè)置為白場
             count1=count1+1;                    %此次循環(huán)新點(diǎn)數(shù)增1
             s1=s1+I(i1+u,j1+v);                   %此點(diǎn)的灰度之加入s中
          end
        end  
       end
      end
     end
   end
 end
suit1=suit1+count1;                         %將n加入符合點(diǎn)數(shù)計數(shù)器中
sum1=sum1+s1;                           %將s加入符合點(diǎn)的灰度值總合中
seed1=sum1/suit1;                          %計算新的灰度平均值
end
sum2=seed2;                              %儲存符合區(qū)域生長的的灰度值
suit2=1;                                   %儲存符合區(qū)域生長條件的點(diǎn)的個數(shù)
count2=1;                                 %記錄每次判斷一點(diǎn)周圍八點(diǎn)符合條件的新點(diǎn)的數(shù)目
threshold2=0.06255;                         %域值
while count2>0
 s2=0;                                   %記錄判斷一點(diǎn)周圍八點(diǎn)時，符合條件的新點(diǎn)的灰度值之和
 count2=0;
 for i2=1:M
   for j2=ceil(N/2):N
     if Z(i2,j2)~=0
      if (i2-1)>0 & (i2+1)<(M+1) & (j2-1)>(ceil(N/2))& (j2+1)<(N+1)
 %判斷此點(diǎn)是否為圖像邊界上的點(diǎn)
       for u=-1:1                            %判斷點(diǎn)周圍八點(diǎn)是否合和域值條件
        for v=-1:1                                       %u,v為偏移量
          if  Z(i2+u,j2+v)==0 & abs(I(i2+u,j2+v)-seed2)<=threshold2    
%判斷是否未存在于輸出矩陣Y，并且為符合域值條件的點(diǎn)
             Z(i2+u,j2+v)=I(i2+u,j2+v);         %符合以上兩條件即將其在Y中與之位置對應(yīng)的點(diǎn)設(shè)置為白場
             count2=count2+1;                 %此次循環(huán)新點(diǎn)數(shù)增1
             s2=s2+I(i2+u,j2+v);                %此點(diǎn)的灰度之加入s中
          end
        end  
       end
      end
     end
   end
 end
suit2=suit2+count2;                             %將n加入符合點(diǎn)數(shù)計數(shù)器中
sum2=sum2+s2;                               %將s加入符合點(diǎn)的灰度值總合中
seed2=sum2/suit2;                              %計算新的灰度平均值
end
figure,imshow(Y+Z),title('分割后白質(zhì)圖像')