\chapter{Matlab Code}
\label{ch:mathlab}

\section{Main Code}
\begin{verbatim}
clear
clc

load DataFile1.mat
I=D1;
[n1,n2]=size(D1);
sigma1=10;
sigma2=30;

figure(1)
imagesc(I(:,:));
title('Initial Image');
saveas(figure(1),['Initial Image DataFile5.png']);

for i=1:6;
     for j=1:n1;
        for k=1:n2;
            mid = mean1(j,k,I);
            if rand>=0.5;
                I(j,k) = random('normal',mid,sigma1,1,1);
            else
                I(j,k)=I(j,k);
            end;
        end;
    end;

    W=random('normal',0,sigma2,n1,n2);
    D=I+W;

I2 = Gibbs(I,D,sigma1,sigma2);
I=I2;

figure(2)
subplot(2,3,i);
imagesc(I); 
figname = sprintf('Image of sweep  %d',i+1);
title (figname);
saveas(figure(2),['Pictures of sweep' int2str(i) ' of DataFile1.mat (sigma=10).png']);
end;

\end{verbatim}

\section{MH Code}
\begin{verbatim}
function [mid] = mean1(j,k,x)
          [n1,n2]=size(x);
            if (j==1) && (k==1)
               mid=(x(j,k+1)+x(j+1,k))/2;
            end;

            if (j==1) && (k==n2)
               mid=(x(1,k-1)+x(j+1,k))/2;
            end;

            if (j==n1) && (k==1)
               mid=(x(j-1,k)+x(j,k+1))/2;
            end;

            if (j==n1) && (k==n2)
               mid=(x(j,k-1)+x(j-1,k))/2;
            end;

            if (j==1 && k~=1 && k~=n2)
               mid=(x(j+1,k)+x(j,k-1)+x(j,k+1))/3;
            end;

            if (j==n1 && k~=1 && k~=n2)
               mid=(x(j-1,k)+x(j,k-1)+x(j,k+1))/3;
            end;

            if (j~=1 && j~=n1 && k==1)
               mid=(x(j-1,k)+x(j+1,k)+x(j,k+1))/3;
            end;

            if (j~=1 && j~=n1 && k==n2)
               mid=(x(j-1,k)+x(j+1,k)+x(j,k-1))/3;
            end;

            if (j~=1 && j~=n1) && (k~=1 && k~=n2)
               mid=(x(j-1,k)+x(j+1,k)+x(j,k+1)+x(j,k+1))/4;
            end;
\end{verbatim}

\section{Gibbs Code}
\begin{verbatim}
function [gib] = Gibbs(I,D,sigma1,sigma2)
[n1,n2]=size(I);
    for j=1:n1;
        for k=1:n2;
 		mid=mean1(j,k,I);
        mu = (mid/sigma1+D(j,k)/sigma2)*(1/((1/sigma1)^2 +(1/sigma2)^2));
        sd = sqrt(1/((1/sigma1)^2+(1/sigma2)^2));
        gib(j,k)=random('normal',mu,sd,1,1);
        end
    end
\end{verbatim}