% Given a set of grid points that form a patch in 2D (boundary and interior), how to find boundarie points with the correct ordering so that the patch function can be used to draw this region.
function outer_points = get_outer_points(grid_position, indicators, delta_x, delta_y)

    patch_squareform = reshape(indicators,sqrt(length(indicators)), sqrt(length(indicators)));
    contour_squareform = zeros(size(patch_squareform));

    for i=1:size(patch_squareform,1)
        ind_tmp = (patch_squareform(i,:)~=[-1,patch_squareform(i,1:end-1)] | patch_squareform(i,:)~=[patch_squareform(i,2:end),-1]) & patch_squareform(i,:)==1;
        contour_squareform(i,ind_tmp)=1;
    end
    for j=1:size(patch_squareform,2)
        ind_tmp = (patch_squareform(:,j)~=[-1;patch_squareform(1:end-1,j)] | patch_squareform(:,j)~=[patch_squareform(2:end,j);-1]) & patch_squareform(:,j)==1;
        contour_squareform(ind_tmp,j)=1;
    end

    contour_squareform = contour_squareform(:);
    outer_points = grid_position(contour_squareform==1,:);

    % order them
    ordered_points = outer_points(1,:); 
    flag_outer_points = zeros(size(outer_points,1),1);flag_outer_points(1)=1;
    while sum(flag_outer_points)~=length(flag_outer_points)
        dist = sum((repmat(ordered_points(end,:),size(outer_points,1),1) - outer_points).^2,2);
        dist(flag_outer_points==1) = max(dist) + 1;
        [dummy,I] = min(dist);
        if dummy > sqrt(delta_x^2+delta_y^2)*2;
            flag_outer_points(I)=1; continue;
        end
        ordered_points = [ordered_points; outer_points(I,:)]; flag_outer_points(I)=1;
    end
    outer_points = ordered_points;

end