function [mat2] = lsmb (mat2,mat)

viewing_distance = 300; % millimeters
sample_spacing = 0.4; % millimeters (0.0635 corresponds to 400dpi)
N=256; %(size of the image)
f=pi/180/asin(1/sqrt(1+viewing_distance^2))*sqrt(((0:2*N-1)'*ones(1,2*N)).^2+(ones(2*N,1)*(0:2*N-1)).^2)/2*N/sample_spacing;
MTF = 2.2*(0.192+0.114*f).*exp(-(0.114*f).^1.1); % Daly model of the eye
pointspread = ifft2(MTF);
mat_ols = [[mat zeros(N,N)];zeros(N,2*N)];
mat2_ols = [[mat2 zeros(N,N)];zeros(N,2*N)];
original = fft2(mat_ols);
halftone = fft2(mat2_ols);
original_HVS = ifft2(original.*MTF);
original_HVS = original_HVS(1:N,1:N);
halftone_HVS = ifft2(halftone.*MTF);
halftone_HVS = halftone_HVS(1:N,1:N);
erreur_HVS = sum(sum((original_HVS-halftone_HVS).^2));
%cross correlations
cpp = ifft2(abs(MTF.^2));
cpe = ifft2(MTF.*conj((original-halftone).*MTF));
%cpp = xcorr2(pointspread);
%cpe = xcorr2(pointspread,erreur_HVS);
err=real(erreur_HVS);
err_last = err/2;

while ((err-err_last)/err>0.1)
    err_last=err;
    i=1;
    j=1;
    % Top left corner    
    a0=-mat2(i,j)*2+1;
    de(1)=est_err_toggle(a0,cpp,cpe,i,j,N);
    de(2)=est_err_swap(a0,2,cpp,cpe,i,j,N);
    de(3)=est_err_swap(a0,3,cpp,cpe,i,j,N);
    de(4)=est_err_swap(a0,4,cpp,cpe,i,j,N);
    de(5:9)=0;
    [e,k]=min(de);
    if min(de)<0
        if k==1
            mat2(i,j)=1-mat2(i,j);
            err=err+de;
            cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i:2*N-i,N-j+1:2*N-j);
        else
            [mat2,m,n]=swap(i,j,k,mat2);
            err=err+de;
            cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j)+a0*cpp(N-m+1:2*N-m,N-n+1:2*N-n);
        end
    end
    % First row
    for j=2:N-1
        a0=-mat2(i,j)*2+1;
        de(1)=est_err_toggle(a0,cpp,cpe,i,j,N);
        de(2)=est_err_swap(a0,2,cpp,cpe,i,j,N);
        de(3)=est_err_swap(a0,3,cpp,cpe,i,j,N);
        de(4)=est_err_swap(a0,4,cpp,cpe,i,j,N);
        de(5)=est_err_swap(a0,5,cpp,cpe,i,j,N);
        de(6)=est_err_swap(a0,6,cpp,cpe,i,j,N);
        de(7:9)=0;
        [e,k]=min(de);
        if min(de)<0
            if k==1
                mat2(i,j)=1-mat2(i,j);
                err=err+de;
                cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j);
            else
                [mat2,m,n]=swap(i,j,k,mat2);
                err=err+de;
                cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j)+a0*cpp(N-m+1:2*N-m,N-n+1:2*N-n);
            end
        end
    end
    j=N;
    % Top right corner
    a0=-mat2(i,j)*2+1;
    de(1)=est_err_toggle(a0,cpp,cpe,i,j,N);
    de(4)=est_err_swap(a0,4,cpp,cpe,i,j,N);
    de(5)=est_err_swap(a0,5,cpp,cpe,i,j,N);
    de(6)=est_err_swap(a0,6,cpp,cpe,i,j,N);
    de(7:9)=0;
    de(2:3)=0;
    [e,k]=min(de);
    if min(de)<0
        if k==1
            mat2(i,j)=1-mat2(i,j);
            err=err+de;
            cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j);
        else
            [mat2,m,n]=swap(i,j,k,mat2);
            err=err+de;
            cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j)+a0*cpp(N-m+1:2*N-m,N-n+1:2*N-n);
        end
    end
    
    
    for i=2:N-1
        % First column
        j=1;
        a0=-mat2(i,j)*2+1;
        de(1)=est_err_toggle(a0,cpp,cpe,i,j,N);
        de(2)=est_err_swap(a0,2,cpp,cpe,i,j,N);
        de(3)=est_err_swap(a0,3,cpp,cpe,i,j,N);
        de(4)=est_err_swap(a0,4,cpp,cpe,i,j,N);
        de(8)=est_err_swap(a0,8,cpp,cpe,i,j,N);
        de(9)=est_err_swap(a0,9,cpp,cpe,i,j,N);
        de(5:7)=0;
        [e,k]=min(de);
        if min(de)<0
            if k==1
                mat2(i,j)=1-mat2(i,j);
                err=err+de;
                cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j);
            else
                [mat2,m,n]=swap(i,j,k,mat2);
                err=err+de;
                cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j)+a0*cpp(N-m+1:2*N-m,N-n+1:2*N-n);
            end
        end
        % entire image
        for j=2:N-1
            a0=-mat2(i,j)*2+1;
            de(1)=est_err_toggle(a0,cpp,cpe,i,j,N);
            de(2)=est_err_swap(a0,2,cpp,cpe,i,j,N);
            de(3)=est_err_swap(a0,3,cpp,cpe,i,j,N);
            de(4)=est_err_swap(a0,4,cpp,cpe,i,j,N);
            de(5)=est_err_swap(a0,5,cpp,cpe,i,j,N);
            de(6)=est_err_swap(a0,6,cpp,cpe,i,j,N);
            de(7)=est_err_swap(a0,7,cpp,cpe,i,j,N);
            de(8)=est_err_swap(a0,8,cpp,cpe,i,j,N);
            de(9)=est_err_swap(a0,9,cpp,cpe,i,j,N);
            [e,k]=min(de);
            if min(de)<0
                if k==1
                    mat2(i,j)=1-mat2(i,j);
                    err=err+de;
                    cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j);
                else
                    [mat2,m,n]=swap(i,j,k,mat2);
                    err=err+de;
                    cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j)+a0*cpp(N-m+1:2*N-m,N-n+1:2*N-n);
                end
            end
        end
        j=N-1;
        % last column
        a0=-mat2(i,j)*2+1;
        de(1)=est_err_toggle(a0,cpp,cpe,i,j,N);
        de(4)=est_err_swap(a0,4,cpp,cpe,i,j,N);
        de(5)=est_err_swap(a0,5,cpp,cpe,i,j,N);
        de(6)=est_err_swap(a0,6,cpp,cpe,i,j,N);
        de(7)=est_err_swap(a0,7,cpp,cpe,i,j,N);
        de(8)=est_err_swap(a0,8,cpp,cpe,i,j,N);
        de(9)=0;
        de(2:3)=0;
        [e,k]=min(de);
        if min(de)<0
            if k==1
                mat2(i,j)=1-mat2(i,j);
                err=err+de;
                cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j);
            else
                [mat2,m,n]=swap(i,j,k,mat2);
                err=err+de;
                cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j)+a0*cpp(N-m-1:2*N-m,N-n-1:2*N-n);
            end
        end
        i
    end
    
    i=N;
    j=1;
    % Bottom left corner
    a0=-mat2(i,j)*2+1;
    de(1)=est_err_toggle(a0,cpp,cpe,i,j,N);
    de(2)=est_err_swap(a0,2,cpp,cpe,i,j,N);
    de(8)=est_err_swap(a0,8,cpp,cpe,i,j,N);
    de(9)=est_err_swap(a0,9,cpp,cpe,i,j,N);
    de(3:7)=0;
    [e,k]=min(de);
    if min(de)<0
        if k==1
            mat2(i,j)=1-mat2(i,j);
            err=err+de;
            cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j);
        else
            [mat2,m,n]=swap(i,j,k,mat2);
            err=err+de;
            cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j)+a0*cpp(N-m+1:2*N-m,N-n+1:2*N-n);
        end
    end
    % Last row
    for j=2:N-1
        a0=-mat2(i,j)*2+1;
        de(1)=est_err_toggle(a0,cpp,cpe,i,j,N);
        de(2)=est_err_swap(a0,2,cpp,cpe,i,j,N);
        de(6)=est_err_swap(a0,6,cpp,cpe,i,j,N);
        de(7)=est_err_swap(a0,7,cpp,cpe,i,j,N);
        de(8)=est_err_swap(a0,8,cpp,cpe,i,j,N);
        de(9)=est_err_swap(a0,9,cpp,cpe,i,j,N);
        de(3:5)=0;
        [e,k]=min(de);
        if min(de)<0
            if k==1
                mat2(i,j)=1-mat2(i,j);
                err=err+de;
                cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j);
            else
                [mat2,m,n]=swap(i,j,k,mat2);
                err=err+de;
                cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j)+a0*cpp(N-m+1:2*N-m,N-n+1:2*N-n);
            end
        end
    end
    j=N;
    % Bottom right corner
    a0=-mat2(i,j)*2+1;
    de(1)=est_err_toggle(a0,cpp,cpe,i,j,N);
    de(6)=est_err_swap(a0,6,cpp,cpe,i,j,N);
    de(7)=est_err_swap(a0,7,cpp,cpe,i,j,N);
    de(8)=est_err_swap(a0,8,cpp,cpe,i,j,N);
    de(2:5)=0;
    de(9)=0;
    [e,k]=min(de);
    if min(de)<0
        if k==1
            mat2(i,j)=1-mat2(i,j);
            err=err+de;
            cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j);
        else
           [mat2,m,n]=swap(i,j,k,mat2);
            err=err+de;
            cpe(N:2*N-1,N:2*N-1) = cpe(N:2*N-1,N:2*N-1) - a0*cpp(N-i+1:2*N-i,N-j+1:2*N-j)+a0*cpp(N-m+1:2*N-m,N-n+1:2*N-n);
        end
    end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%  Auxilary functions %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [de] = est_err_toggle(a0,cpp,cpe,i,j,N)
de = cpp(N,N)-2*a0*cpe(N+i-1,N+j-1);


function [de] = est_err_swap(a0,k,cpp,cpe,i,j,N)
switch k
case 2
    de = cpp(N,N)-2*a0*cpe(N+i-1,N+j-1)+2*a0*cpe(N+i-1,N+j)-2*cpp(N,N+1);
case 3
    de = cpp(N,N)-2*a0*cpe(N+i-1,N+j-1)+2*a0*cpe(N+i,N+j)-2*cpp(N+1,N+1);
case 4
    de = cpp(N,N)-2*a0*cpe(N+i-1,N+j-1)+2*a0*cpe(N+i,N+j-1)-2*cpp(N+1,N);
case 5
    de = cpp(N,N)-2*a0*cpe(N+i-1,N+j-1)+2*a0*cpe(N+i,N+j-2)-2*cpp(N+1,N-1);
case 6
    de = cpp(N,N)-2*a0*cpe(N+i-1,N+j-1)+2*a0*cpe(N+i-1,N+j-2)-2*cpp(N,N-1);
case 7
    de = cpp(N,N)-2*a0*cpe(N+i-1,N+j-1)+2*a0*cpe(N+i-2,N+j-2)-2*cpp(N-1,N-1);
case 8
    de = cpp(N,N)-2*a0*cpe(N+i-1,N+j-1)+2*a0*cpe(N+i-2,N+j-1)-2*cpp(N-1,N);
otherwise
    de = cpp(N,N)-2*a0*cpe(N+i-1,N+j-1)+2*a0*cpe(N+i-2,N+j)-2*cpp(N-1,N+1);
end



function [mat2,m,n] = swap (i,j,k,mat2)
x=mat2(i,j);
switch k
case 2
    mat2(i,j)=mat2(i,j+1);
    mat2(i,j+1)=x;
    m=i;n=j+1;
case 3
    mat2(i,j)=mat2(i+1,j+1);
    mat2(i+1,j+1)=x;
    m=i+1;n=j+1;
case 4
    mat2(i,j)=mat2(i+1,j);
    mat2(i+1,j)=x;
    m=i+1;n=j;
case 5
    mat2(i,j)=mat2(i+1,j-1);
    mat2(i+1,j-1)=x;
    m=i+1;n=j-1;
case 6
    mat2(i,j)=mat2(i,j-1);
    mat2(i,j-1)=x;
    m=i;n=j-1;
case 7
    mat2(i,j)=mat2(i-1,j-1);
    mat2(i-1,j-1)=x;
    m=i-1;n=j-1;
case 8
    mat2(i,j)=mat2(i-1,j);
    mat2(i-1,j)=x;
    m=i-1;n=j;
otherwise
    mat2(i,j)=mat2(i-1,j+1);
    mat2(i-1,j+1)=x;
    m=i-1;
    n=j+1;
end