function [J,thresh_plot] = segment_particle_RegionGrow(I1,x1,y1,parameters)

    % REGION GROWING FUNCTION

    % Stop growth if it hits the side of the image
    I3 = I1 < 0; I3(:,1) = 1; I3(:,end) = 1; I3(1,:) = 1; I3(end,:) = 1;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Growing particle criterion
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Region growing code
    err = 0;
    x = x1(1); y = y1(1);
    reg_size = 1; % Number of pixels in region
    J = false(size(I1));
    J(x,y) = 1;
    thresh_min = 1; thresh_max = 100000; thresh_step = 1; min_maxdist = 1;

    % Calculate local density weighting function
    clear w1
    if parameters.f_ceiling
        w1(1:parameters.nhood_size^2) = ceil(exp(1:parameters.nhood_size...
            ^2).^-1/(exp(parameters.nhood_median)^-1));
        w1(w1 > 1) = ((w1(w1 > 1) - 1)./parameters.w1_reduce) + 1;
    else
        w1(1:parameters.nhood_size^2) = (exp(1:parameters.nhood_size...
            ^2).^-1/(exp(parameters.nhood_median)^-1));
        w1(w1 > 1) = ((w1(w1 > 1) - 1)./parameters.w1_reduce) + 1;
        w1(w1 < 1) = 1 - ((1 - w1(w1 < 1))./parameters.w1_reduce);
    end

    i = thresh_min;
    thresh_count = 1;
    while i < thresh_max
        thresh_count = thresh_count + 1; i = i + thresh_step;
    end
    thresh_plot = zeros(thresh_count,5);

    % Start regiogrowing until distance between region and possible new pixels become
    % higher than a certain threshold

    decreasing = 0;
    J_max = J;
    
    % Start of the region growing loop.  Continually increase the
    % iterations until the intensity difference of the region is optimized
    for iter = 1:thresh_count

        % What is the maximum intensity difference between region and
        % boundary pixels
        reg_maxdist = thresh_min + (iter - 1)*thresh_step;      
        if min_maxdist > reg_maxdist
            criteria_check = 0;
            thresh_flag = 0;
        else
            criteria_check = 1;
            thresh_flag = 1;
        end

        % If there are neighboring pixels that can be added, then execute
        % this loop 
        while criteria_check && ~err

            % B is the logical binary image of the boundary pixels for
            % possible addition
            B = logical(imdilate(J,ones(3))-J);

            % Add weighting function that kicks in after a certain number
            % of pixels have been grown
            h = ones(parameters.nhood_size);
            B1 = imfilter(uint8(J),h);
            B1(~B) = 0;
            W = double(zeros(size(B1)));
            for i = 1:parameters.nhood_size^2
                W(B1 == i) = w1(i);
            end

            % Add pixel with intensity nearest to the mean of the region,
            % to the region
            D = abs(double(I1) - mean(double(I1(J))));
            if reg_size > parameters.min_size
                D = D .* W; 
                % multiply all values by the local density weighting
                % function if the region size is greater than the specified
                % minimum region size 
            end
            D(~B) = 0;
           
            % Go ahead and add the pixel if it is below the intensity
            % threshold and neighbors the region.  Also, update the size of
            % the region.
            D1 = D == min(D(B));
            D1(~B) = 0;
            Dn = find(D1);
            if D(Dn(1)) <= reg_maxdist
                J(Dn(1)) = 1;
                reg_size = reg_size + 1;
            else
                min_maxdist = D(Dn(1));
                criteria_check = 0;
            end

            % If the filling holes option is selected, then fill holes
            % within the region at a predetermined interval to insure
            % continuity of the region
            if (parameters.fill_holes_flag && mod(reg_size,parameters.fill_holes_mod) == 0) ...
                    || criteria_check == 0
                J1 = imfill(J,'holes');
                J = J1;
                reg_size = sum(J(:));
            end

            % Show the progress at predetermined intervals
            if parameters.output_flag && mod(reg_size,parameters.output_mod) == 0
                figure(1), imshow(J,[])
            end

            % THRESHOLD CRITERIA
            % If the region size is too big, stop growing
            if reg_size > parameters.max_size, criteria_check = 0; end
            % If the region has grown into the edges of the image, stop
            if sum(J(I3)) > 0, criteria_check = 0; err = 1;...
                if parameters.output_flag, figure(1), imshow(J,[]), end 
            end

        end

        % If the region has grown this iteration and there were no
        % errors...
        if thresh_flag && ~err
            % Calculate the mean intensity areas that will be used for the
            % stopping criterion
            if parameters.whole_region
                % Mean intensity of the whole region
                mean_region = mean(I1(J));  
            else
                % Mean intensity of a boundary layer inside of particle    
                J3 = imdilate((imdilate(J,ones(parameters.boundary_layer))...
                    -J),ones(parameters.boundary_layer)); 
                J3(~J) = 0;
                mean_region = mean(I1(logical(J3)));
            end
            % Mean intensity of the boundary
            J3 = logical(imdilate(J,ones(parameters.boundary_layer))-J);
            mean_boundary = mean(I1(J3));   % boundary outside region
            % Store the values
            thresh_plot(iter,:) = [reg_maxdist, sum(J(:)), ...
                mean_boundary, mean_region, abs(mean_boundary-mean_region)];

            % This little routine produces images at each increasing
            % threshold value to see the progression of the region growing
            % technique
            if parameters.image_growth_flag
                clear J1
                Jtmp = I1; Jtmp(bwperim(J)) = 255;
                J1(:,:,1) = Jtmp; 
                Jtmp = I1; Jtmp(bwperim(J)) = 0;
                J1(:,:,2) = Jtmp; J1(:,:,3) = Jtmp;
                imwrite(J1,strrep(parameters.image_growth_file,'.tif',...
                    strrep('_i.tif','_i',['_' num2str(iter,'%03.0f')])),'tif')
                clear Jtmp J1
            end

            % This routine will shut down the region growing algorithm
            % after a certain number of decreasing intensity difference
            % values following an intensity maximum
            if thresh_plot(iter,5) >= max(thresh_plot(1:iter-1,5))
                J_max = J;
                decreasing = 0;
            else
                decreasing = decreasing + 1;
            end
            if decreasing >= parameters.decreasing_max, err = 1; end
        
        % ... or if nothing has happened to the region
        elseif ~thresh_flag
            thresh_plot(iter,:) = thresh_plot(iter-1,:);
        end
    end

    % The optimum segmentation is the maximum intensity difference, i.e.,
    % that was stored in variable J_max.
    J = J_max;

    % If wanted, the routine can show a figure with the final segmented
    % region superimposed on the original image
    if parameters.trace_boundary
        trace_boundary(I1,J)
    end

    % This option will plot the threshold plots that include the intensity
    % means for the region and outside the region as well as the intensity
    % differences as a function of the threshold value
    if parameters.threshplot
        nn = find(thresh_plot(:,1)~=0);
        figure
        plot(thresh_plot(nn,1),thresh_plot(nn,2),'-ko','MarkerFaceColor','k'), hold on
        h = get(gcf,'CurrentAxes');
        set(h,'FontName','times','FontSize',24);
        axis square
        xlabel({'Threshold Value'}, ...
            'fontsize',24,...
            'fontweight','b', 'FontName','times')
        ylabel({'Region Size (pixels)'}, ...
            'fontsize',24,...
            'fontweight','b', 'FontName','times')

        figure
        plot(thresh_plot(nn,1),thresh_plot(nn,5),'-ro','MarkerFaceColor','r'), hold on
        plot(thresh_plot(nn,1),thresh_plot(nn,3),'-bo','MarkerFaceColor','b')
        plot(thresh_plot(nn,1),thresh_plot(nn,4),'-ko','MarkerFaceColor','k')
        h = get(gcf,'CurrentAxes');
        set(h,'FontName','times','FontSize',24);
        axis square
        xlabel({'Threshold Value'}, ...
            'fontsize',24,...
            'fontweight','b', 'FontName','times')
        ylabel({'Mean Intensity Value'}, ...
            'fontsize',24,...
            'fontweight','b', 'FontName','times')
    end
        
end