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)
M=8;
f=pi/180/asin(1/sqrt(1+viewing_distance^2))*sqrt(((0:M-1)'*ones(1,M)).^2+(ones(M,1)*(0:M-1)).^2)/M/sample_spacing;
MTF = 2.2*(0.192+0.114*f).*exp(-(0.114*f).^1.1); % Daly model of the eye
pointspread = ifft2(MTF);

xMin = -5:4%minutes of visual angle - corresponding to the viewing distance and the sample spacing
X = xMin(ones(1,length(xMin)),:); Y = X';
D = X.^2 + Y.^2; D = D.^0.5;

ps = 0.952*exp(-2.59*abs(D).^1.36) + 0.048*exp(-2.43*abs(D).^1.74);

a = conv2(mat,ps);
original_HVS = a(6:6+N;6:6+N);
a = conv2(mat2,ps);
halftone_HVS = a(6:6+N;6:6+N);

erreur_HVS = sum(sum((original_HVS-halftone_HVS).^2));
%cross correlations
cpp = xcorr2(ps);
cpe = xcorr2(ps,erreur_HVS);
%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+e;
            cpe(i:2*M+i-1,j:2*M+j-1) = cpe(i:2*M+i-1,j:2*M+j-1) - a0*cpp;
        else
            [mat2,m,n]=swap(i,j,k,mat2);
            err=err+e;
            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+e;
                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+e;
                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+e;
            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+e;
            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+e;
                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+e;
                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+e;
                    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+e;
                    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+e;
                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+e;
                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+e;
            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+e;
            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+e;
                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+e;
                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+e;
            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+e;
            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