半色调技术

这个主要是为了对付作业，对各处找到的资料缝缝补补的整合，因为比较杂乱，已经忘记是在哪儿看到大佬们的文章了，整理一下思路。

背景

半色调图像如常见的印刷品图像，其由浅到深或由淡到浓的变化，是靠网点面积大小或网点覆盖率来表现的。一般用于复制诸如照片之类的连续调原稿时，会采用这种半色调技术，它将图像分成许多点，通过点的不同大小来表现颜色的深浅。
在印刷品的印刷时，印刷机用有限数量的一套油墨(只有黑墨或青、品红、黄、黑墨)来印刷数量不同、大小不同的细小点，印刷品画面上色彩和浓淡就是靠这些细小的点来表示的，由此可以给人眼产生许多灰度级或许多颜色的错觉。当观察印品画面时，网点面积大，颜色就深，称为暗调；网点面积小，颜色就浅，则称为亮调。由于网点在空间上是有一定的距离的，呈离散型分布，并且由于加网的线数总有一定的限制，在图像的层次变化上不能像连续调图像一样实现无级变化，故称加网图像为半色调图像。
目前半色调技术最普遍的处理方式主要有：抖动法、误差扩散法、迭代法。本次主要通过学习抖动法和误差扩散法了解半色调的实现和应用。

算法原理

1.1．抖动法：
抖动法是点处理类方法的一种典型算法，主要分为随机抖动和有序抖动两大类。
这两种算法都需要一个模板，也称为抖动矩阵或阈值矩阵，抖动矩阵不仅决定了当亮度或灰度值减小时网点变成黑点的顺序.而且也决定了半色调图像的质量，所以抖动算法的关键是抖动矩阵的构造。该算法与抖动矩阵进行比较，矩阵中的每个阈值的取值范围是图像的最大灰度值和最小灰度值之间，其数学公式化如下图所示：
在这里插入图片描述

式中f(i,J)代表连续色调图像中的像素点灰度值，r(i,J)代表抖动矩阵的阈值，而h(i,J)代表半色调后的图像灰度值。
1.1.1.随机抖动矩阵：
随机抖动矩阵是通过完全随机产生的，所以半色调后的图像质量常常很不理想，在实际中已经基本不再使用。
1.1.2有序抖动的抖动矩阵：
有序抖动的抖动矩阵是有规律的，具有良好的图像效果和高效的处理速度而被各大打印机厂商采用，有序抖动矩阵主要有两种类型：分散型和聚集型。典型的分散型抖动矩阵是Bayer有序抖动矩阵，而点局部聚簇整体分散是典型的聚集型矩阵，如下图所示：
在这里插入图片描述

Bayer有序抖动矩阵设计：
设计标准图案的算法，是由 Limb 在 1969 年提出的。
先设计出一个2*2矩阵：
在这里插入图片描述

递归关系：
在这里插入图片描述
通过迭代可以获得不同规模大小的有序抖动矩阵，M3被称为Bayer抖动表。
聚簇型矩阵设计：
该类型的模板设计主要的特征是阈值的大小从中心点位置向四周的方向逐渐增大，本次设计中主要比较了4*4大小的cluster和dispersedCL两种聚簇型抖动矩阵。
在这里插入图片描述
规则抖动的优点是算法简单；缺点是图案化有时很明显，这是因为取模运算虽然引入了随机成分，但还是有规律的。另外，点之间进行比较时，只要比标准图案上点的值大就打白点，这种做法并不理想，因为，如果当标准图案点的灰度值本身就很小，而图象中点的灰度只比它大一点儿时，图象中的点更接近黑色，而不是白色。一种更好的方法是将这个误差传播到邻近的象素。

1.2误差扩散法：
在有序抖动处理中，利用了像素点与抖动矩阵比较来判断是否在一个位置放置点，实质是一种点处理过程。通过模板移动的形式会对成像的图像效果带来周期性的网格条纹。在1976年Floyd和Steinberg提出了误差扩散算法，它将半色调加网从“点处理”过渡到“邻域处理”。
误差扩散算法的提出为半色调加网带来了革命性的技术变革，也是半色调技术上的里程碑，并促进了半色调技术的飞速发展。通过误差扩散处理后的半色调图像像素分布各异且无规律性，色调丰富，视觉效果较好。直到目前为止，它依然被视为易于实现且视觉效果较好的半色调技术之一，被广泛应用.
给定阈值t，设原图像为x(i,j)，输出为y(i,j)。对整幅图像按从左到右、从上到下的顺序逐点的依次执行两步操作：
阈值化y(i,j):
h(i,j)=x(i,j)-y(i,j)
b.把量化误差扩散到邻近的未被处理过的点。量化误差是指输入与输出之间的差：量化误差扩散即改变空间上相邻像素x(i,j+1)，x(i+1,j-1)，x(i+1,j))，x(i+1,j+1)的值，将当前像素的量化误差按比例转移并叠加到邻近的像素
换句话说，误差扩散相当于把中值阈值法在每一个点上产生的误差再加到周围的点上，从而保持局部区域的总体灰度基本不变。量化误差扩散的参数可以用一个矩阵描述，称之为误差扩散过滤器，设为P，P={Nr，s}，(r，s)D为P中所有非零元素对应的下标集合。这里的下标r，s是整数，可以为正、为负或为0。这样上述公式可写为：
在这里插入图片描述

本次设计中验证和比较的过滤器模板如下图所示：
在这里插入图片描述

实现代码

clc
clear;
%% 
I=imread('sample.jpeg'); %读取图像
A=rgb2gray(I);
[length,width]=size(A);
%% 抖动法
%随机抖动
RandomJitter=fix(255*rand(8,8));
%% 分散性抖动
% bayer 有序抖动
m1 = [[0 2];[3 1]];
u1=ones(2, 2);
m2=[[4*m1 4*m1+2*u1];[4*m1+3*u1 4*m1+u1]];
u2=ones(4, 4);
m3=[[4*m2 4*m2+2*u2];[4*m2+3*u2 4*m2+u2]];%8*8
bayer3=m3;
u3=ones(8, 8);
m4=[[4*m3 4*m3+2*u3];[4*m3+3*u3 4*m3+u3]];%16*16
bayer4=m4;
% 长这样：Bayer=[0,48,12,60,3,51,15,63;
%        32,16,44,28,35,19,47,31;
%         8,56,4,52,11,59,7,55;
%        40,24,36,20,43,27,39,23;
%        2,50,14,62,1,49,13,61;
%        34,18,46,30,33,17,45,29;
%        10,58,6,54,9,57,5,53;
%        42,26,38,22,41,25,37,21]%聚簇型抖动
%cluster
cluster=[6,7,8,9;5,0,1,10;4,3,2,11;15,14,13,12];%dispersedCL
dispersedCL=[0,4,2,6;12,8,14,10;3,7,1,5;15,11,13,9];clusterImg=dithercount(I,cluster);
dispersedCLimg=dithercount(I,dispersedCL);
bayer3Img=dithercount(I,bayer3);
bayer4Img=dithercount(I,bayer4);
randomditterimg=dithercount(I,RandomJitter);
% subplot(141);
figure(1)
imshow(bayer3Img)
title('bayer8*8') 
hold on
% subplot(142);
figure(2)
imshow(bayer4Img);
title('bayer16*16') 
hold on
% subplot(143);
figure(3)
imshow(clusterImg);
title('cluster')
hold on
% subplot(144);
figure(4)
imshow(dispersedCLimg);
title('dispersedCL')
hold on
figure(5)
imshow(randomditterimg);
title('random')
hold on
%% 误差扩撒
[length,width]=size(A);
zero=zeros(length,1);
B=[zero,A,zero];
zero=zeros(1,width+2);
B=[B;zero];
[length,width]=size(B);
% floyd-steinberg
a=7/16;
b=3/16;
c=5/16;
d=1/16;
for i=1:(length-1)%对每行进行处理for j=2:(width-1)%对每列进行处理if (B(i,j)/2)<127 %进行阈值比较floydimg(i,j)=0;elsefloydimg(i,j) = 255;endwc=(B(i,j)-floydimg(i,j));%定义误差B(i,j+1)=B(i,j+1)+wc*a;B(i+1,j-1)=B(i+1,j-1)+wc*b;B(i+1,j)=B(i+1,j)+wc*c;B(i+1,j+1)=B(i+1,j+1)+wc*d;%误差分配end
end
floydimg=floydimg(1:end-1,2:end-1);
figure(6)
imshow(floydimg);
title('floyd-steinberg')
hold on
%% burkers
[length,width]=size(A);
C=A;
zero=zeros(length,2);
C=[zero,C,zero];
zero=zeros(1,width+2);
burkersimg=zeros(length,width);
for i=3:(length-2)%对每行进行处理for j=2:(width-2)%对每列进行处理if (C(i,j)/1.8)<128 %进行阈值比较burkersimg(i,j)=0;elseburkersimg(i,j) = 255;endwc=(C(i,j)-burkersimg(i,j));%定义误差C(i+1,j)=C(i+1,j)+wc*(8/32);C(i+2,j)=C(i+2,j)+wc*(4/32);C(i-2,j+1)=C(i-2,j+1)+wc*(2/32);C(i-1,j+1)=C(i-1,j+1)+wc*(4/32);C(i,j+1)=C(i,j+1)+wc*(8/32);C(i+1,j+1)=C(i+1,j+1)+wc*(4/32); C(i+2,j+1)=C(i+2,j+1)+wc*(2/32); end
end
burkersimg=burkersimg(:,2:end);
figure(7)
imshow(burkersimg);
title('burkers')
hold on
%% sierra
[length,width]=size(A);
a=5/32;
b=4/32;
c=3/32;
d=2/32;
C=A;
sierraimg=zeros(length,width);
for i=3:(length-2)%对每行进行处理for j=2:(width-2)%对每列进行处理if (C(i,j)/1.5)<127 %进行阈值比较sierraimg(i,j)=0;elsesierraimg(i,j) = 255;endwc=(C(i,j)-sierraimg(i,j));%定义误差C(i+1,j)=C(i+1,j)+wc*a;C(i+2,j)=C(i+2,j)+wc*c;C(i-2,j+1)=C(i-2,j+1)+wc*d;C(i-1,j+1)=C(i-1,j+1)+wc*b;C(i,j+1)=C(i,j+1)+wc*a;C(i+1,j+1)=C(i+1,j+1)+wc*b; C(i+2,j+1)=C(i+2,j+1)+wc*d; C(i-1,j+2)=C(i-1,j+2)+wc*d;C(i,j+2)=C(i,j+2)+wc*c;C(i+1,j+2)=C(i+1,j+2)+wc*d;end
end
figure(8)
imshow(sierraimg);
title('sierra')
hold on
%% stucki
[length,width]=size(A);
a=7/42;
b=5/42;
c=3/42;
d=1/42;
D=A;
stuckiimg=zeros(length,width);
for i=3:(length-2)%对每行进行处理for j=2:(width-2)%对每列进行处理if (D(i,j)/2)<127%进行阈值比较stuckiimg(i,j)=0;elsestuckiimg(i,j) = 255;endwc=(D(i,j)-stuckiimg(i,j));%定义误差D(i+1,j)=D(i+1,j)+wc*a;D(i+2,j)=D(i+2,j)+wc*b;D(i-2,j+1)=D(i-2,j+1)+wc*c;D(i-1,j+1)=D(i-1,j+1)+wc*b;D(i,j+1)=D(i,j+1)+wc*a;D(i+1,j+1)=D(i+1,j+1)+wc*b; D(i+2,j+1)=D(i+2,j+1)+wc*c; D(i-2,j+2)=D(i-2,j+2)+wc*d;D(i-1,j+2)=D(i-1,j+2)+wc*c;D(i,j+2)=D(i,j+2)+wc*b;D(i+1,j+2)=D(i+1,j+2)+wc*c;D(i+2,j+2)=D(i+2,j+2)+wc*d;end
end
figure(9)
imshow(stuckiimg);
title('stucki')
hold on
%% jarvis
a=7/48;
b=5/48;
c=3/48;
d=1/48;
D=A;
jarvisimg=zeros(length,width);
for i=3:(length-2)%对每行进行处理for j=2:(width-2)%对每列进行处理if (D(i,j)/1.8)<127%进行阈值比较jarvisimg(i,j)=0;elsejarvisimg(i,j) = 255;endwc=(D(i,j)-jarvisimg(i,j));%定义误差D(i+1,j)=D(i+1,j)+wc*a;D(i+2,j)=D(i+2,j)+wc*b;D(i-2,j+1)=D(i-2,j+1)+wc*c;D(i-1,j+1)=D(i-1,j+1)+wc*b;D(i,j+1)=D(i,j+1)+wc*a;D(i+1,j+1)=D(i+1,j+1)+wc*b; D(i+2,j+1)=D(i+2,j+1)+wc*c; D(i-2,j+2)=D(i-2,j+2)+wc*d;D(i-1,j+2)=D(i-1,j+2)+wc*c;D(i,j+2)=D(i,j+2)+wc*b;D(i+1,j+2)=D(i+1,j+2)+wc*c;D(i+2,j+2)=D(i+2,j+2)+wc*d;end
end
figure(10)
imshow(jarvisimg);
title('jarvis')
hold on%% 
%抖动法的模板应用
function f=dithercount(image,template)
[length,width,height]=size(image);
if height==1gI=rgb2gray(image);
elsegI=image;
end
bw =zeros(length,width);
[teml,temw]=size(template);
if teml==4m=16;
elseif teml==8m=4;
elseif teml==16m=1;
end
for i=1:lengthfor j=1:widthif (gI(i,j)/m> template(bitand(i, teml-1) + 1, bitand(j,temw-1) + 1))bw(i,j)= 255;elsebw(i,j)= 0;endend
end
f=bw;
end