function labelMap = myGrow2D(dF,labelMap,sigma,growTimes)
    % grow
    [H,W] = size(dF);
%     dF = dF - median(dF(:));
%     dF = dF/max(dF(:));
%     labelMap = dF>thr;
    %% boundary
    candidate_pix = [];
    dw = zeros(9,1);
    dh = zeros(9,1);
    cnt = 1;
    for x = -1:1
        for y = -1:1
            dw(cnt) = x;
            dh(cnt) = y;
            cnt = cnt + 1;
        end
    end

    pix0 = find(labelMap);

    for kk = 1:9
        [ih0,iw0] = ind2sub([H,W],pix0);
        ih1 = max(min(ih0+dh(kk),H),1);
        iw1 = max(min(iw0+dw(kk),W),1);
        pix1 = sub2ind([H,W],ih1,iw1);
        candidate_pix = [pix1;candidate_pix];
    end
    candidate_pix = setdiff(candidate_pix,pix0);
    candidate_pix = unique(candidate_pix);

    % neighborhood
    radius = 2;
    filter = zeros(radius*2+1,radius*2+1);
    for x = -radius:radius
        for y = -radius:radius
            if(x^2 + y^2 <=radius^2)
                filter(x+radius+1,y+radius+1) = 1;
            end
        end
    end

    dF0 = dF;
    dF0(~labelMap) = 0;
    neighbor_record = imfilter(uint8(labelMap),filter);
    summation = imfilter(dF0,filter);
    dF0 = dF;
    dF0(labelMap) = 0;
%     neighbor_record2 = imfilter(uint8(~labelMap),filter);
%     summation2 = imfilter(dF0,filter);
    clear dF0;
    labelMap = false(H,W);
    Max_rounds = growTimes; % 30 can be changed. Just an arbitrary number.
    for nn = 1:Max_rounds   
        summation2 = zeros(H,W);
        neighbor_record2 = ones(H,W);
        [pix0,candidate_pix,neighbor_record,summation,neighbor_record2,summation2] = growRegion2D(pix0,...
            dF,candidate_pix,filter,neighbor_record,summation,neighbor_record2,summation2,sigma);
%         fprintf('Grow %d rounds\n',nn);
    end

    labelMap(pix0) = true;
% 
%     fprintf('Done\n')

    
%     imShow = zeros(2*H,W);
%     imShow(1:H,:) = dF;
%     imShow(H+1:end,:) = labelMap;
%     zzshow(imShow);
end





