% Segment based on area, Region Growing;
clear  all; close all; clc
[fileName,pathName] = uigetfile('*.*','Please select an image');%檔案筐，選擇檔案
if(fileName)
    fileName = strcat(pathName,fileName);
    fileName = lower(fileName);%一致的小寫字母形式
else 
    J = 0;%記錄區域生長所分割得到的區域
    msgbox('Please select an image');
    return; %退出程式
end

I = imread(fileName);
if( ~( size(I,3)-3 ))
    I = rgb2gray(I);%轉化為單通道灰度圖
end
I = im2double(I); %影象灰度值歸一化到[0,1]之間
Ireshape = imresize(I,[600,800]);
I = Ireshape(51:475,200:699);
gausFilter = fspecial('gaussian',[5 5],0.5);
I = imfilter(I,gausFilter,'replicate');

%種子點的互動式選擇
if( exist('x','var') == 0 && exist('y','var') == 0)
    subplot(2,2,1),imshow(I,[]);
    hold on;
    [y,x] = getpts;%滑鼠取點  回車確定
    x = round(x(1));%選擇種子點
    y = round(y(1));
end

if( nargin == 0)
    reg_maxdist = 0.1;
    %nargin是matlab程式碼編寫中常用的一個技巧，主要用於計算當前主函式的輸入引數個
    %數，一般可以根據nargin的返回值來確定主函式輸入引數的預設值。在實現中，如果
    %使用者輸入的引數個數為零，那麼預設為0.2
end
J = zeros(size(I)); % 主函式的返回值，記錄區域生長所得到的區域
Isizes = size(I);
reg_mean = I(x,y);%表示分割好的區域內的平均值，初始化為種子點的灰度值
reg_size = 1;%分割的到的區域，初始化只有種子點一個
neg_free = 10000; %動態分配記憶體的時候每次申請的連續空間大小
neg_list = zeros(neg_free,3);
%定義鄰域列表，並且預先分配用於儲存待分析的畫素點的座標值和灰度值的空間，加速
%如果影象比較大，需要結合neg_free來實現matlab記憶體的動態分配
neg_pos = 0;%用於記錄neg_list中的待分析的畫素點的個數
pixdist = 0;
%記錄最新畫素點增加到分割區域後的距離測度
%下一次待分析的四個鄰域畫素點和當前種子點的距離
%如果當前座標為（x,y）那麼通過neigb我們可以得到其四個鄰域畫素的位置
neigb = [ -1 0;
          1  0;
          0 -1;
          0  1];
 %開始進行區域生長，當所有待分析的鄰域畫素點和已經分割好的區域畫素點的灰度值距離
 %大於reg_maxdis,區域生長結束
 
 while (pixdist < 0.06 && reg_size < numel(I))
     %增加新的鄰域畫素到neg_list中
     for j=1:4
         xn = x + neigb(j,1);
         yn = y + neigb(j,2);
         %檢查鄰域畫素是否超過了影象的邊界
         ins = (xn>=1)&&(yn>=1)&&(xn<=Isizes(1))&&(yn<=Isizes(1));
         %如果鄰域畫素在影象內部，並且尚未分割好；那麼將它新增到鄰域列表中
         if( ins && J(xn,yn)==0)
             neg_pos = neg_pos+1;
             neg_list(neg_pos,:) =[ xn, yn, I(xn,yn)];%儲存對應點的灰度值
             J(xn,yn) = 1;%標註該鄰域畫素點已經被訪問過 並不意味著，他在分割區域內
         end
     end
    %如果分配的記憶體空問不夠，申請新的記憶體空間
    if (neg_pos+10>neg_free)
        neg_free = neg_free + 100000;
        neg_list((neg_pos +1):neg_free,:) = 0;
    end
    %從所有待分析的畫素點中選擇一個畫素點，該點的灰度值和已經分割好區域灰度均值的
    %差的絕對值時所待分析畫素中最小的
    dist = abs(neg_list(1:neg_pos,3)-reg_mean);
    [pixdist,index] = min(dist);
    %計算區域的新的均值
    reg_mean = (reg_mean * reg_size +neg_list(index,3))/(reg_size + 1);
    reg_size = reg_size + 1;
    %將舊的種子點標記為已經分割好的區域畫素點
    J(x,y)=2;%標誌該畫素點已經是分割好的畫素點
    x = neg_list(index,1);
    y = neg_list(index,2);
%     pause(0.0005);%動態繪製
%     if(J(x,y)==2)
%     plot(x,y,'r.');
%     end
    %將新的種子點從待分析的鄰域畫素列表中移除
    neg_list(index,:) = neg_list(neg_pos,:);
    neg_pos = neg_pos -1;
 end
 
 J = (J==2);%我們之前將分割好的畫素點標記為2
 hold off;
 subplot(2,2,2),imshow(J);
 J = bwmorph(J,'dilate');%補充空洞
 subplot(2,2,3),imshow(J);
 subplot(2,2,4),imshow(I+J);