function [mat2]=halftone_color_visual (type,mat,thresh,nois,serp,BNM_size,sharp)
switch type
case 1
    pixel= [[0;0;0] [255;0;0] [0;255;0] [0;0;255] [255;255;0] [255;0;255] [0;255;255] [255;255;255]];
    pixel_HSV = opponent(pixel')';
    [l1 l2 l3]=size(mat);
    matHSV = opponent(reshape(mat,l1*l2,l3));
    weights = [1,1,1]';
    mat2=matHSV;
    for i=1:l1*l2
        [p,k]=min(abs(matHSV(i,:)'*ones(1,8)-pixel_HSV)'*weights);
        mat2(i,:) = pixel(:,k)'/255;
    end
    mat2=reshape(mat2,l1,l2,l3);
case 2
    white=round(255*rand(size(mat,1),size(mat,2),3));
    mat2=halftone_color_visual(1,mat-white,thresh,noise,serp,BNM_size,sharp);
case 3
    vc=process_color(BNM_size);
    vc2=process(BNM_size);
    h=fspecial('lagrangian');
    m = conv2(mat(:,:,1),h);
    m=m(2:size(m)-1,2:size(m)-1);
    mat(:,:,1)=mat(:,:,1)-sharp*m;
    mat(:,:,1)=max(0,min(255,mat(:,:,1)));
    mat2(:,:,1)=temp(vc(:,:,1),mat(:,:,1));

    m = conv2(mat(:,:,2),h);
    m=m(2:size(m)-1,2:size(m)-1);
    mat(:,:,2)=mat(:,:,2)-sharp*m;
    mat(:,:,2)=max(0,min(255,mat(:,:,2)));
    mat2(:,:,2)=temp(vc(:,:,2),mat(:,:,2));

    m = conv2(mat(:,:,3),h);
    m=m(2:size(m)-1,2:size(m)-1);
    mat(:,:,3)=mat(:,:,3)-sharp*m;
    mat(:,:,3)=max(0,min(255,mat(:,:,3)));
    mat2(:,:,3)=temp(vc(:,:,3),mat(:,:,3));

case 5
    FS = [0 0 0 7 5;3 5 7 5 3;1 3 5 3 1];
    FS = FS/sum(sum(FS));
    mat2=err_diff(FS,mat,thresh);
otherwise
    mat2=mat;
end

function [mat2]=temp (template,mat)

[b1 b2]=size(template);
[l1 l2]=size(mat);
mat0=mat;
mat=zeros(b1*ceil(l1/b1),b2*ceil(l2/b2));
mat(1:l1,1:l2)=mat0;
x=max(max(template));
for i=1:ceil(l1/b1)
    for j=1:ceil(l2/b2)
        mat2(((i-1)*b1+1):(i*b1),((j-1)*b2+1):(j*b2))=(mat(((i-1)*b1+1):(i*b1),((j-1)*b2+1):(j*b2))>((template+1/2/255*(x<1)+1/2*(x>=1))*255/x));
    end
end
mat2=mat2(1:l1,1:l2);


function [mat2]=err_diff(FS,mat,thresh)
[img_r, img_c,d] = size(mat);
xFS_r = size(FS, 1);
xFS_c = fix(size(FS, 2)/2);
pixel= [[0;0;0] [1;0;0] [0;1;0] [0;0;1] [1;1;0] [1;0;1] [0;1;1] [1;1;1]];
pixel_HSV = opponent(pixel')';

temp = zeros(img_r+xFS_r, img_c+2*xFS_c,3);
temp(1:img_r, xFS_c +1:xFS_c+img_c,:) = mat/255;
load monitor;load xyz;
m1 = xyz'*phosphors;
m2=[.3897 .6890 -.0787;-.2298 1.1834 .0464;0 0 1];
m3=[1 1 1;1 -1 1;1 0 -1];
w=[2 0 0;0 .5 0; 0 0 .5];
m4=inv(m3);
m5=inv(m2);
m6=inv(m1);
trans = m1*m2*m3*w*m4*m5*m6/8;
trans1 = 7/16*trans;trans1=[0.6316 -0.1306 0.0323;-0.0430 0.3993 0.0327;-0.0167 -0.1082 0.7379];
trans2 = 1/16*trans;trans2=[-0.1949 0.1289 -0.0242; 0.0817 -0.0730 0.0645;0.0454 0.1585 -0.4017];
trans3 = 5/16*trans;trans3=[0.3598 -0.0549 0.0403; -0.0018 0.2906 0.0173;-0.0080 -0.0895 0.4867];
trans4 = 3/16*trans;trans4=[0.2181 -0.0112 0.0047; 0.0222 0.1515 0.0580;0.0129 0.0213 0.1614];


for ir = 1 : img_r,
    
    for ic = xFS_c+1 : img_c,
        error = reshape(temp(ir, ic,:),1,3);
        temp(ir, ic,:) = (error>thresh/255);
        error = error-reshape(temp(ir, ic,:),1,3);
        temp(ir,ic+1,:) = reshape(temp(ir,ic+1,:),1,3) + error*trans1;
        temp(ir+1,ic+1,:) = reshape(temp(ir+1,ic+1,:),1,3) + error*trans2;
        temp(ir+1,ic,:) = reshape(temp(ir+1,ic,:),1,3) + error*trans3;
        temp(ir+1,ic-1,:) = reshape(temp(ir+1,ic-1,:),1,3) + error*trans4;
    end
    
    temp(ir:ir+xFS_r-1,img_c+1:img_c+xFS_c,:) = ...
        temp(ir:ir+xFS_r-1,img_c+1:img_c+xFS_c,:)+temp(ir+1:ir+xFS_r, 1:xFS_c,:);
    
    for ic = img_c+1 : img_c+xFS_c,
        error = reshape(temp(ir, ic,:),1,3);
        temp(ir, ic,:) = (error>255);
        error = error-reshape(temp(ir, ic,:),1,3);
        temp(ir,ic+1,:) = reshape(temp(ir,ic+1,:),1,3) + error*trans1;
        temp(ir+1,ic+1,:) = reshape(temp(ir+1,ic+1,:),1,3) + error*trans2;
        temp(ir+1,ic,:) = reshape(temp(ir+1,ic,:),1,3) + error*trans3;
        temp(ir+1,ic-1,:) = reshape(temp(ir+1,ic-1,:),1,3) + error*trans4;
    end
    
    temp(ir+1:ir+xFS_r, xFS_c+1:2*xFS_c,:) = ...
        temp(ir+1:ir+xFS_r, xFS_c+1:2*xFS_c,:) + temp(ir:ir+xFS_r-1,  ...
        img_c+xFS_c+1:img_c+2*xFS_c,:);
    
    temp(:, 1:xFS_c,:) = zeros(img_r+xFS_r, xFS_c,3);
    temp(:, img_c+xFS_c+1:img_c+2*xFS_c,:) = zeros(img_r+xFS_r, xFS_c,3);
end
mat2 = temp(1:img_r, xFS_c+1:xFS_c+img_c,:);

function mat1 = opponent(mat)
mat1=mat;
m = [1 1 1;1 -1 1;1 0 -1];
mat1 = mat*m;