% PRML image de-noising

clc;

clear;

close all;

A=imread('a.jpg');%读入名字为a.jpg的图片

imshow(A);

M=;

N=;%把图片处理一下大小

A=imresize(A,[M,N]);

X = rgb2gray(A);

for i=:M%%二值化图像，要根据实验图像来改天阈值，本实验阈值为200

    for j=:N

        if(X(i,j)<)%200为阈值，请您根据实际实验图片改变

            X(i,j)=;

        else

            X(i,j)=;

        end

    end

end

Y=X%把原图X保存在Y里

imshow(X)%显示X

for i=:M%这两个for循环是改变二值化图像10%点的值。

    for j=:N

        if(rand()<0.1)%以百分之10的概率进行改变

            if(Y(i,j)==)

               % Y(i,j)=;

               Y(i,j)=;

            else

                Y(i,j)=;

        end

        end

    end

end

figure;

imshow(Y);%显示带有噪声的图像

YY=zeros(M,N)

for i=:M%把{,}转换为{-，+}

    for j=:N

              if  Y(i,j)==

                YY(i,j)=;

              else

                YY(i,j)=-;

              end

    end

end

%参数设置

beta=1.0;

yita=2.1;

h=;

%step1

    R=YY; %R是要逼近X的图像，YY是噪声图像

%step2

Change=

while Change  %系统扫描法,可以试试随机选点法的效果

    Change=;

    for i=:M-

        for j=:N-

            temp=R(i,j);

            %若这个点状态为+，计算这一点的能量。

            %解释一下为什么是2*beta..因为这个点不仅影响自己的能量函数，也影响了周围四个点的能量函数。

            %所以-*beta=-beat*(..)-beat*(...).第一项相当于它自己的能量函数，第二项相当于周围能量函数的一部分。

            %R(i,j)不会对其他点的能量造成影响，所以只考虑（i,j)点的能量与其组成团的点的部分能量函数即可。

            R(i,j)=;

        Epos=h*R(i,j)-*beta*(R(i+,j)*R(i,j)+R(i-,j)*R(i,j)+R(i,j+)*R(i,j)+R(i,j-)*R(i,j))-yita*R(i,j)*YY(i,j);

             R(i,j)=-;

        Eneg=h*R(i,j)-*beta*(R(i+,j)*R(i,j)+R(i-,j)*R(i,j)+R(i,j+)*R(i,j)+R(i,j-)*R(i,j))-yita*R(i,j)*YY(i,j);

        if Epos~=Eneg

        R(i,j)=1.0*(Epos<Eneg)+(Epos>Eneg)*-1.0;

        else

           R(i,j)=temp;

        end

            if temp~=R(i,j)

                Change=;%若有变化则继续while1的系统扫描，若每个点都没有改变则结束程序

            end

        end

    end

end

for i=:M

    for j=:N

        if(R(i,j)==)

            Y(i,j)=;

        else

            Y(i,j)=;

        end

    end

end

figure;

imshow(Y)

%错误的概率

 disp( ['error rate is %d  ' num2str(sum(sum(Y~=X))/(M*N))])