一种基于最小显色误差激光打印机模型绿噪声半色调算法

摘要

本发明公开了一种基于最小显色误差激光打印机模型绿噪声半色调算法，属于图像硬拷贝复制领域的半色调网点生成方法。由于打印机具有点增益和点丢失的非线性扭曲存在，需要使用打印机模型对非线性扭曲进行模拟。而目前的打印机模型通常使用密度或者反射率来换算测量结果，由于密度和反射率并不符合人眼视觉灰度的线性变化，所以通常需要多次打印测量才能确定打印机模型的参数。本发明将打印输出结果转换为孟塞尔明度值，计算打印机模型的最优参数，得到最小显色误差的激光打印机模型，并应用到绿噪声半色调算法中。只需要一次打印测量，所生成的最优参数对应的打印机模型就较好地抑制了打印机的非线性扭曲。

说明书

技术领域

本发明属于数字图像处理领域，涉及一种半色调处理方法，具体涉及一种基于最小显色误差激光打印机模型绿噪声半色调算法。

背景技术

由于墨粉物理扩散及静电荷吸附能力不足，打印机会出现点增益和点丢失现象，这两种现象又被称为打印机非线性扭曲。打印机非线性扭曲将会导致输出图像的灰度级并级、相对明暗关系变化等情况。为了抑制打印机的非线性扭曲，国内外研究者通过模拟打印机打印点的特征及打印机打印机理，建立打印机模型。

目前的打印机模型，通常是将测量值转换为密度值或者反射率，通过比较转换后的密度值、反射率与原图像的灰度值的关系来确定打印机模型。但由于密度值和反射率并不符合人眼视觉灰度变化情况，所以并不能直接通过最小方差的方法求出打印机模型的最优参数，而需要反复打印测量，通过人眼来判定最优参数。

发明内容

为了解决上述技术问题，本发明提出一种基于最小显色误差激光打印机模型绿噪声半色调算法，将打印输出结果转换为孟塞尔明度值，计算打印机模型的最优参数，得到最小显色误差的激光打印机模型，并应用到绿噪声半色调算法中，只需要一次打印测量，所生成的最优参数对应的打印机模型就较好地抑制了打印机的非线性扭曲。

本发明所采用的技术方案是：一种基于最小显色误差激光打印机模型绿噪声半色调算法，其特征在于，包括以下步骤：

步骤1：通过计算测试图像的最小显色误差，来获得对应的最优参数T₁、T₂组合；其具体实现包括以下子步骤：

步骤1.1：绘制包含由黑到白等灰度值变化的20个色块的灰度图像作为测试图像I_s；

步骤1.2：选用绿噪声半色调算法对测试图像I_s进行半色调处理，获得半色调处理结果图像I_h；

步骤1.3：将半色调处理结果图像I_h使用目标激光打印机进行打印输出，获得输出结果图像I_p；

步骤1.4：使用色度计或密度计测量输出结果图像I_p，将测量数据转换为灰度值，获得20个色块测量灰度值数组arrPG；其中所述的将测量数据转换为灰度值，其具体实现包括以下子步骤：

步骤1.4.1：将测量数据转换为CIE1931XYZ色度系统中的亮度因数Y；其中：

若使用色度计对打印输出结果进行测量，则测量得到打印输出结果的亮度因数Y；

若使用密度计对打印输出结果进行测量，则根据光学反射密度定义公式推导出的密度D与亮度因数Y的函数关系，计算出亮度因数Y：

$>D=\lg [Y_0/\left(Y{\rho }_0\right)]\Rightarrow Y=\frac{Y_0}{{10}^D{\rho }_0}$>    (式壹)

其中，Y₀为孟塞尔系统第10级明度的亮度因数，Y₀=102.75，ρ₀为孟塞尔系统第10级明度的反射率，ρ₀=1，D为某区域的光学反射密度，Y为该区域的亮度因数；

步骤1.4.2：将亮度因数Y转换为孟塞尔明度V，根据亮度因数Y与孟塞尔明度V的函数关系，计算孟塞尔明度V：

V=2.217Y^0.352-1.324    (式贰)

步骤1.4.3：将孟塞尔明度V转换为灰度值G，

$>G=\frac{V-V_p}{V_s-V_p}$>    (式叁)

其中，V为待转换为孟塞尔明度，V_p为纸白对应的孟塞尔明度，V_s为100%网点打印效果对应的孟塞尔明度；

步骤1.5：针对半色调处理结果图像I_h，在[0,1]范围内以0.1为步长，不断地改变式肆中参数将单个像素分割成10×10个小点，计算单个像素曝光时在各个小点的光能量分布值P(x_i,y_i)：

$>P(x_i,y_i)={\exp }^{-\partial [{(x_i-x_c)}^2+{(y_i-y_c)}^2]/l^2}$>    (式肆)

其中，(x_c,y_c)为像素点(m,n)的曝光中心，(x_i,y_i)为该像素曝光中心周围的各个小点，l为激光打印机两个相邻设备打印点的中心距，l=10；

步骤1.6：统计曝光像素及其八邻域的各个小点的光能量，计算曝光像素本身区域、相邻区域、对角区域的平均光能量值A_center、A_nerb、A_diag，其通用计算公式为：

$>A=\frac{{\int }_{y_i-l/2}^{y_i+l/2}{\int }_{x_i-l/2}^{x_i+l/2}P(x,y)\mathrm{dxdy}}{l\times l}$>    (式伍)

其中，(x_i-l/2,x_i+l/2)、(y_i-l/2,y_i+l/2)分别为单位面积区域的横纵坐标取值范围；

步骤1.7：逐点判断半色调处理结果图像I_h中的各个像素的八邻域是否为打印点，计算各个像素区域上的叠加光能量L：

L=uA_center+vA_nerb+wA_diag    (式陆)

其中，u由当前像素是否被曝光决定，A_center为当前像素在本身像素点区域的平均光能量值，v为当前像素的邻近区域中将被曝光的像素个数，A_nerb为邻近方向的打印像素点在当前像素区域的平均光能量值，w为当前像素的对角区域中将被曝光的像素个数，A_diag为对角的打印点在当前像素区域的平均光能量值；

步骤1.8：在[0,3]范围内以0.1为步长，依次改变式柒中参数T₁、T₂，代入式柒中计算各个像素的模拟输出灰度值p

$>p=\left(\begin{array}{cc}0& L<T_1\\ \frac{L-T_1}{T_2-T_1}& T_1\le L<T_2\\ 1& T_2\le L\end{array}\right)$>    (式柒)

其中，如果在变换阈值过程中，T₁>T₂，则此次变换跳过，直接进入下一次变换；

步骤1.9：计算同一灰度色块内的所有像素对应的模拟灰度值p的平均值，获得20个灰度色块的模拟灰度值数组arrIG；

步骤1.10：计算测量灰度值数组arrPG和模拟灰度值数组arrIG之间的均方根误差RMSE：

$>\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{n}{\Sigma }}{({\mathrm{arrPG}}_i-{\mathrm{arrIG}}_i)}^2}{n}}$>    (式捌)

其中，n为灰度色块的个数，n=20；

步骤1.11：判断，RMSE是否是最小值？

如果否，则回转执行所述的步骤1.5；

如果是：则顺序执行下述步骤1.12；

步骤1.12：选择最小显色误差RMSE_min对应的参数T₁、T₂作为最优参数；

步骤2：将最优参数T₁、T₂载入到绿噪声半色调算法中，对图像进行半色调处理；其具体实现包括以下子步骤：

步骤2.1：代入最优参数并根据式肆及式伍分别计算单个打印点在本身像素区域、邻近区域和对角区域的光能量值A_center、A_nerb和A_diag，并保存；

步骤2.2：设待处理图像为一幅大小为M×N的灰度图像，x_m,n表示待处理图像的像素点(m,n)处的灰度值，其中x_m,n∈[0,1]，0为黑，1为白；将周围八邻域已处理像素点的误差值与误差扩散滤波器卷积得到误差值，当前点的输入像素的灰度值与误差值之差作为输入阈值量化器灰度值u_m,n：

$>u_{m,n}=x_{m,n}-\underset{k,j\in S}{\Sigma }{F}_{k,j}{e}_{m-k,n-j}$>    （式玖）

其中，S是误差扩散滤波器F(·)的系数模板范围，k、j为F(·)的系数模板的行列值；

步骤2.3：采用迟滞滤波器H(·)对像素点(m,n)处相邻的已处理像素点的模拟灰度值进行权值运算，再乘以迟滞系数h作为反馈，加入输入阈值量化器灰度值u_m,n中，则像素点(m,n)的输出值b_m,n：

$>b_{m,n}=Q\left(u_{m,n}\right)=\left(\begin{array}{cc}1& u_{m,n}-h\underset{p,q\in R}{\Sigma }{H}_{p,q}{b}_{m-p,n-q}^{\prime }>0.5\\ 0& \mathrm{else}\end{array}\right)$>    （式拾）

其中，R是迟滞滤波器H(·)的系数模板范围，输出值1为打印点，0为不打印点；p、q为H(·)的系数模板的行列值；

步骤2.4：将光能量值A_center、A_nerb和A_diag代入式陆，计算像素点(m,n)区域上的光能量叠加值，其中式陆中的u、v、w的取值由像素点(m,n)八邻域中已处理点的输出值决定。

$>u=\left(\begin{array}{cc}1& b_{m,n}=0\\ 0& b_{m,n}=1\end{array}\right)$>    (式十一)

    （式十二）

    （式十三）

步骤2.5：将最优参数T₁、T₂代入式柒中计算当前点的模拟灰度值b′_m,n；

步骤2.6：计算当前点的误差值e_m,n为模拟像素值与输入阈值量化器灰度值的差值e_m,n：

e_m,n=b′_m,n-u_m,n    (式十四)

步骤2.7：按步骤2.1至步骤2.6所述，从待处理图像左上角点开始处理，按左到右、从上到下的顺序，直至到待处理图像右下角点结束。

本发明的创新点在于：

（1）将测量数据转换为孟塞尔明度值，将孟塞尔明度转换为模拟灰度；

（2）通过循环修改参数T₁、T₂，求最小误差，以对应的参数作为最优参数；

（3）将最优参数模拟得到的灰度值代替原始灰度值计算误差；

本发明将打印输出结果转换为孟塞尔明度值，计算打印机模型的最优参数，得到最小显色误差的激光打印机模型，并应用到绿噪声半色调算法中，只需要一次打印测量，所生成的最优参数对应的打印机模型就较好地抑制了打印机的非线性扭曲。

附图说明

图1：是本发明的基于最小显色误差最优参数计算的流程图。

图2：是本发明的基于最小显色误差激光打印机模型绿噪声半色调算法的流程图

图3：是本发明实施例的亮度因数Y与孟塞尔明度V的转换关系示意图。

图4：是本发明实施例的待处理的20个灰度色块。

图5-1：是本发明实施例中所采用的误差扩散滤波器F(·)的系数模板。

图5-2：是本发明实施例中所采用的迟滞滤波器H(·)的系数模板。

图6：是本发明实施例的lena测试图像。

图7：是本发明实施例的由普通绿噪声半色调算法处理并打印扫描的lena测试图。

图8：是本发明实施例的由基于最小显色误差激光打印机模型绿噪声半色调算法处理并打印扫描的lena测试图。

具体实施方式

下面以Pantum2000激光打印机为目标打印机、至冠牌高白复印纸进行打印输出、测量设备为Eyeone色度测量仪，以lena测试图为待处理图像为例，以普通绿噪声半色调算法为对比，使用1200dpi的高分分辨率扫描仪对半色调打印输出图像进行扫描，以此对本发明作进一步地描述。

请见图1、图2，本发明所采用的技术方案是：一种基于最小显色误差激光打印机模型绿噪声半色调算法，其特征在于，包括以下步骤：

步骤1：通过计算测试图像的最小显色误差，来获得对应的最优参数T₁、T₂组合；其具体实现包括以下子步骤：

步骤1.1：请见图4，绘制600×600像素由0到0.95每0.05灰度变化的20个灰度色块作为测试图像I_s，其中0表示黑，1表示白；

步骤1.2：选用绿噪声半色调算法对测试图像I_s进行半色调处理，获得半色调处理结果图像I_h；

步骤1.3：将半色调处理结果图像I_h使用Pantum2000激光打印机进行打印输出，获得输出结果图像I_p；

步骤1.4：使用Eyeone色度测量仪测量输出结果图像I_p，并按以下方法将测量数据转换为灰度值，获得20个色块测量灰度值数组arrPG：其中所述的将测量数据转换为灰度值，其具体实现包括以下子步骤：

步骤1.4.1：将测量数据转换为CIE1931XYZ色度系统中的亮度因数Y；其具体实现过程为：根据光学反射密度定义公式推导出的密度D与亮度因数Y的函数关系，计算出亮度因数Y：

$>D=\lg [Y_0/\left(Y{\rho }_0\right)]\Rightarrow Y=\frac{Y_0}{{10}^D{\rho }_0}$>    (式壹)

其中，Y₀为孟塞尔系统第10级明度的亮度因数，Y₀=102.75，ρ₀为孟塞尔系统第10级明度的反射率，ρ₀=1，D为某区域的光学反射密度，Y为该区域的亮度因数；

步骤1.4.2：将亮度因数Y转换为孟塞尔明度V，请见图3，根据亮度因数Y与孟塞尔明度V的函数关系，计算孟塞尔明度V：

V=2.217Y^0.352-1.324    (式贰)

步骤1.4.3：将孟塞尔明度V转换为灰度值G，

$>G=\frac{V-V_p}{V_s-V_p}$>    (式叁)

其中，V为待转换为孟塞尔明度，V_p为纸白对应的孟塞尔明度，V_s为100%网点打印效果对应的孟塞尔明度；

步骤1.5：针对半色调处理结果图像I_h，在[0,1]范围内以0.1为步长，不断地改变式肆中参数将单个像素分割成10×10个小点，计算单个像素曝光时在各个小点的光能量分布值P(x_i,y_i)：

$>P(x_i,y_i)={\exp }^{-\partial [{(x_i-x_c)}^2+{(y_i-y_c)}^2]/l^2}$>    (式肆)

其中，(x_c,y_c)为像素点(m,n)的曝光中心，(x_i,y_i)为该像素曝光中心周围的各个小点，l为激光打印机两个相邻设备打印点的中心距，l=10；

步骤1.6：统计曝光像素及其八邻域的各个小点的光能量，计算曝光像素本身区域、相邻区域、对角区域的平均光能量值A_center、A_nerb、A_diag，其通用计算公式为：

$>A=\frac{{\int }_{y_i-l/2}^{y_i+l/2}{\int }_{x_i-l/2}^{x_i+l/2}P(x,y)\mathrm{dxdy}}{l\times l}$>    (式伍)

其中，(x_i-l/2,x_i+l/2)、(y_i-l/2,y_i+l/2)分别为单位面积区域的横纵坐标取值范围；

步骤1.7：逐点判断半色调处理结果图像I_h中的各个像素的八邻域是否为打印点，计算各个像素区域上的叠加光能量L：

L=uA_center+vA_nerb+wA_diag    (式陆)

其中，u由当前像素是否被曝光决定，A_center为当前像素在本身像素点区域的平均光能量值，v为当前像素的邻近区域中将被曝光的像素个数，A_nerb为邻近方向的打印像素点在当前像素区域的平均光能量值，w为当前像素的对角区域中将被曝光的像素个数，A_diag为对角的打印点在当前像素区域的平均光能量值；

步骤1.8：在[0,3]范围内以0.1为步长，依次改变式柒中参数T₁、T₂，代入式柒中计算各个像素的模拟输出灰度值p

$>p=\left(\begin{array}{cc}0& L<T_1\\ \frac{L-T_1}{T_2-T_1}& T_1\le L<T_2\\ 1& T_2\le L\end{array}\right)$>    (式柒)

其中，如果在变换阈值过程中，T₁>T₂，则此次变换跳过，直接进入下一次变换；

步骤1.9：计算同一灰度色块内的所有像素对应的模拟灰度值p的平均值，获得20个灰度色块的模拟灰度值数组arrIG；

步骤1.10：计算测量灰度值数组arrPG和模拟灰度值数组arrIG之间的均方根误差RMSE：

$>\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{n}{\Sigma }}{({\mathrm{arrPG}}_i-{\mathrm{arrIG}}_i)}^2}{n}}$>    (式捌)

其中，n为灰度色块的个数，n=20；

步骤1.11：判断，RMSE是否是最小值？

如果否，则回转执行所述的步骤1.5；

如果是：则顺序执行下述步骤1.12；

步骤1.12：选择最小显色误差RMSE_min对应的参数T₁、T₂作为最优参数；本实施例选择最小显色误差RMSE_min=3.28时，对应的参数T₁=0.4、T₂=2.3作为最优参数；

步骤2：将最优参数T₁、T₂载入到绿噪声半色调算法中，对图像进行半色调处理；其具体实现包括以下子步骤：

步骤2.1：代入最优参数并根据式肆及式伍分别计算单个打印点在本身像素区域、邻近区域和对角区域的光能量值A_center、A_nerb和A_diag，并保存；

步骤2.2：设待处理图像为一幅大小为M×N的灰度图像，x_m,n表示待处理图像的像素点(m,n)处的灰度值，其中x_m,n∈[0,1]，0为黑，1为白；将周围八邻域已处理像素点的误差值与误差扩散滤波器卷积得到误差值，当前点的输入像素的灰度值与误差值之差作为输入阈值量化器灰度值u_m,n：

$>u_{m,n}=x_{m,n}-\underset{k,j\in S}{\Sigma }{F}_{k,j}{e}_{m-k,n-j}$>    （式玖）

其中，S是误差扩散滤波器F(·)的系数模板范围，k、j为F(·)的系数模板的行列值，请见图5-1；

步骤2.3：采用迟滞滤波器H(·)对像素点(m,n)处相邻的已处理像素点的模拟灰度值进行权值运算，再乘以迟滞系数h作为反馈，加入输入阈值量化器灰度值u_m,n中，则像素点(m,n)的输出值b_m,n：

$>b_{m,n}=Q\left(u_{m,n}\right)=\left(\begin{array}{cc}1& u_{m,n}-h\underset{p,q\in R}{\Sigma }{H}_{p,q}{b}_{m-p,n-q}^{\prime }>0.5\\ 0& \mathrm{else}\end{array}\right)$>    （式拾）

其中，R是迟滞滤波器H(·)的系数模板范围，输出值1为打印点，0为不打印点；p、q为H(·)的系数模板的行列值，请见图5-2；

步骤2.4：将光能量值A_center、A_nerb和A_diag代入式陆，计算像素点(m,n)区域上的光能量叠加值，其中式陆中的u、v、w的取值由像素点(m,n)八邻域中已处理点的输出值决定。

$>u=\left(\begin{array}{cc}1& b_{m,n}=0\\ 0& b_{m,n}=1\end{array}\right)$>    (式十一)

    （式十二）

    （式十三）

步骤2.5：将最优参数T₁=0.4、T₂=2.3代入式柒中计算当前点的模拟灰度值b′_m,n；

步骤2.6：计算当前点的误差值e_m,n为模拟像素值与输入阈值量化器灰度值的差值e_m,n：

e_m,n=b′_m,n-u_m,n    (式十四)

步骤2.7：按步骤2.1到步骤2.6所述，从待处理图像左上角点开始处理，按左到右、从上到下的顺序，直至到待处理图像右下角点结束。

本实施例中，将测量数据转换为CIE1931XYZ色度系统中的亮度因数Y；也可以使用色度计对打印输出结果进行测量，并且测量得到打印输出结果的亮度因数Y。

请见图6，是本发明实施例的lena测试图像；请见图7，是本发明实施例的由普通绿噪声半色调算法处理并打印扫描的lena测试图；其整体明暗明显偏暗，其头发及帽子上的花的阶调完全并级，无法分辨细节轮廓；请见图8，是本发明实施例的由基于最小显色误差激光打印机模型绿噪声半色调算法处理并打印扫描的lena测试图，其整体明暗还原较符合原图，在暗调部分细节仍能较好地再现，有效地抑制了打印机非线性扭曲。

以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围，因此，凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。