image segmentation using normalized statistics and high-order statistics

 I tried to implement this paper with matlad, but on 2D-grayscal image.

https://www.google.com.eg/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0CBwQFjAA&url=http%253A%252F%252Fwww.ee.nthu.edu.tw%252F~cwlin%252Fpub%252Fiscas02seg.pdf&ei=oVgMVeOZM8zkUceehPgP&usg=AFQjCNEfVXKmgHE4BcnnbPfYQl34cf4H_w&sig2=WnDQ5xkR4Fk4lHV42UU0EQ

 I put a threshold manually to show final segmentation result, where I don't understand how is the probability of object and background calculated to calculate the threshold automatically? I'll be Appreciate if someone explain it to me.

 my code is:
function [] = statistics_seg()
I = imread('t1.png');
I = im2double(rgb2gray(I));
[n m] = size(I);
imm = zeros(size(I));%image mean values
imd = zeros(size(I));%deference image (I-imm) values
imv = zeros(size(I)); %image variance
imv4 = zeros(size(I));%image with 4th variance
norm_statistic = zeros(size(I));%image after normalized statistic
high_statistic = zeros(size(I));%image after 4th variance statistic
f1 = zeros(size(I));
f2 = zeros(size(I));
M = zeros(8,1);%regest the neighbours of pixel
V = zeros(8,1);%regest the neighbours-mean
V4 = zeros(8,1);

 % calculate mean, variance and 4th-order variance
r = 1;
for i = 1:n
    for j = 1:m
        for ii=-1:1
            iwin=ii+i;
            if iwin<1;
                iwin=1;
            end
            if iwin>n;
                iwin=n;
            end
            for jj=-1:1
                jwin=jj+j;
                if jwin<1
                    jwin=1;
                end
                if jwin>m
                    jwin=m;
                end
                if r<= 8
                M(r,1) = I(iwin,jwin);
                else
                    r = 1;
                end
            end            
        end
        %calculate the mean value
        m1 = ((sum(M))+I(i,j))/9;
        imm(i,j) = m1;
        %calculate the variance
        for r = 1:8
        V(r,1) = (M(r,1) - m1)^2;
        V4(r,1) = (M(r,1) - m1)^4;
%         MO3(r,1) = (M(r,1) - m1)^3;
%         MO4(r,1) = (M(r,1) - m1)^4;
        end
        x = (I(i,j)-m1)^2;
        imv(i,j) = (sum(V)+x)/9;
        imdev(i,j) = (imv(i,j))^0.5;%image standard deviation
        imd(i,j) = I(i,j) - imm(i,j);
        norm_statistic(i,j) = imd(i,j)/imv(i,j);

        %find 4th-order variance
        x4 = (I(i,j)-m1)^4;
        imv4(i,j) = (sum(V4)+x4)/9;
        high_statistic(i,j) = imd(i,j)/imv4(i,j);
    end
end
%thresholded the output images
med1 = median(median(norm_statistic))^0.7;
med2 = median(median(high_statistic))^0.7;
for i = 1:n
    for j = 1:m
        if norm_statistic(i,j)>= med1
            f1(i,j) = 1;
        end
        if high_statistic(i,j)>= med2
            f2(i,j) = 1;
        end
    end
end

 imshow(I), title('input image');
figure, imagesc(norm_statistic), title('normalized statistic image');
figure, imagesc(imv4), title('4th-order variance image');
figure, imagesc(high_statistic), title('high statistic image');
figure, imshow(f1), title('final1 thresholded image');
figure, imshow(f2), title('final2 thresholded image');
end

my image

thanks