Halftoning算法

3.1 介绍

黑白打印机通过控制黑色墨水的有无来打印灰度图。由于打印机没有灰色的墨水，所以只能通过控制黑白的疏密来得到灰度图，而相应的算法就叫 Haltoning 算法（半调算法）。

3.2 误差扩散 Halftoning 算法

误差扩散算法是 1976 年由 Floyd 和 Steinberg 提出的，也叫 Floyd-Steinberg 算法。这个算法的基础是“误差扩散”（error-diffusion）。当误差累积到一定程度时，就把那个像素设为 1，然后重置误差。如果误差不够大，就保持像素是 0。这样一来就能近似显示出灰色。

下面用数学语言来描述 Halftoning 算法。假设输入图片 $I(m,n)$ 有 $R$ 行 $C$ 列，输出为 $H(m,n)$；用 $E_p(m,n)$ 表示 $(m,n)$ 位置周围的扩散误差；$E_g(m,n)$ 表示在 $(m,n)$ 位置的误差；

扩散误差如何计算呢？我们定义一个误差分布矩阵 $\boldsymbol{C}$，它的元素和为 1，比如：

\[\boldsymbol{C}= \begin{bmatrix} 0.0 & 0.2 & 0.0\\ 0.6 & 0.1 & 0.1 \end{bmatrix}\]

我们用如下公式该矩阵与 $(m,n)$ 周围像素的误差两两相乘再求和，即可得到扩散误差。

\[E_p(m,n) = \sum_{i=1}^I \sum_{j=1}^J C_{ij} E_g[(m-i+1), (n-j+1)]\\ 其中，\sum_i \sum_j C_{ij} = 1，且 C_{11}=0\]

用上面的例子就是：

\[\begin{bmatrix} 0.1 & 0.1 & 0.6\\ 0.0 & 0.2 & 0.0 \end{bmatrix} \times \begin{bmatrix} E_g(m-1,n-2) & E_g(m-1,n-1) & E_g(m-1,n)\\ E_g(m,n-2) & E_g(m,n-1) & \color\red{ E_g(m,n)} \end{bmatrix}\]

注意我们只考虑 $(m,n)$ 位置左上的元素的扩散误差。如果左上没有元素，我们可以设为 0，或着用 $(m,n)$ 处的误差去替代。另外，根据定义式，$\boldsymbol{C}$ 先旋转 180 再乘 $E_p(m,n)$.

那么总误差 $t$ 可以表示为：

\[t = I(m,n)+E_p(m,n)\]

我们取阈值为灰度的最大值，记为 ${\rm maxval}$，则该像素的输出值 $H(m,n)$ 为：

\[H(m,n)= \begin{cases} 0 & t< {\rm maxval}\\ 1 & t\geq {\rm maxval} \end{cases}\]

$(m,n)$ 处的误差 $E_g(m,n)$ 为：

\[E_g(m,n)= \begin{cases} t & H(m,n)=0\\ t-{\rm maxval} & H(m,n)=1 \end{cases}\]

当 $H(m,n)=1$ 时，$t-{\rm maxval}$ 就是重置误差的操作。

用伪代码描述就是：

Define:
    I(R,C) - 有 R 行 C 列的输入图像
    Ep(m,n) - 其他位置传播到 (m,n) 处的误差和 
    Eg(m,n) - (m,n) 处的总误差
    C(i,j) - i 行 j 列的误差分布函数

1. 初始化 Ep(m,n) = Eg(m,n) = 0，两者都有 R 行 C 列
2. 循环 m=1,R
    1. 循环 n=1,C
        1. 计算其他位置传播到 (m,n) 处的误差和得到扩散误差 Ep(m,n)
        2. 将 (m,n) 处的像素值与扩散误差相加: T = I(m,n) + Ep(m,n)
        3. 若 T > threshold
            则执行 7. 和 8.
            否则执行 9. 和 10.
            1. 将像素 (m,n) 设为“开”
            2. 重置 (m,n) 处的误差：
                Eg(m,n) = T - threshold
            3. 将像素 (m,n) 设为“关”
            4. 保留 (m,n) 处的误差：
                Eg(m,n) = T
    2. 结束 n 的循环
3. 结束 m 的循环

代码如下：

void halftoning(unsigned char **in_image,
                unsigned char **out_image,
                int threshold,
                unsigned char one,
                unsigned char zero,
                int rows, int cols)
{
    float **eg, **ep;
    float c[2][3] = {
        {0.0, 0.2, 0.0},
        {0.6, 0.1, 0.1}};
    float sum_p, t;
    int i, j, m, n, xx, yy;

    /***********************************************
      *
      *   Calculate the total propogated error
      *   at location(m,n) due to prior
      *   assignment.
      *
      *   Go through the input image.  If the output
      *   should be one then display that pixel as such.
      *   If the output should be zero then display it
      *   that way.
      *
      *   Also set the pixels in the input image array
      *   to 1's and 0's in case the print option
      *   was chosen.
      *
      ************************************************/

    eg = malloc(rows * sizeof(float *));
    for (i = 0; i < rows; i++)
    {
        eg[i] = malloc(cols * sizeof(float));
        if (eg[i] == NULL)
        {
            printf("\n\tmalloc of eg[%d] failed", i);
        } /* ends if */
    }     /* ends loop over i */

    ep = malloc(rows * sizeof(float *));
    for (i = 0; i < rows; i++)
    {
        ep[i] = malloc(cols * sizeof(float));
        if (ep[i] == NULL)
        {
            printf("\n\tmalloc of ep[%d] failed", i);
        } /* ends if */
    }     /* ends loop over i */

    for (i = 0; i < rows; i++)
    {
        for (j = 0; j < cols; j++)
        {
            eg[i][j] = 0.0;
            ep[i][j] = 0.0;
        }
    }

    /**********************************************
       *
       *   29 February 1988 - Fix to remove a solid 
       *   line at the bottom of each region. Loop 
       *   over ROWS-1 and then draw an extra line.
       *
       **********************************************/

    for (m = 0; m < rows; m++)
    {
        for (n = 0; n < cols; n++)
        {

            sum_p = 0.0;
            for (i = 0; i < 2; i++)
            {
                for (j = 0; j < 3; j++)
                {

                    xx = m - i;
                    yy = n - j;
                    if (xx < 0)
                        xx = 0;
                    if (xx >= rows)
                        xx = rows - 1;
                    if (yy < 0)
                        yy = 0;
                    if (yy >= cols)
                        yy = cols - 1;

                    sum_p = sum_p + c[i][j] * eg[xx][yy];
                } /* ends loop over j */
            }     /* ends loop over i */

            ep[m][n] = sum_p;
            t = in_image[m][n] + ep[m][n];

            /**********************************
               *
               *    Here set the point [m][n]=one
               *
               ***********************************/

            if (t > threshold)
            {
                eg[m][n] = t - threshold;
                out_image[m][n] = one;
            } /* ends if t > threshold  */

            /**********************************
               *
               *    Here set the point [m][n]=zero
               *
               ***********************************/

            else
            { /* t <= threshold */
                eg[m][n] = t;
                out_image[m][n] = zero;
            } /* ends else t <= threshold  */

        } /* ends loop over n columns */
    }     /* ends loop over m rows */

    for (i = 0; i < rows; i++)
    {
        free(eg[i]);
        free(ep[i]);
    }

} /*  ends halftoning  */

运行结果：

Figure 3.2 Input lena Image	Figure 3.3 Output Halftoned lena Image

多级 Halftone

上面我们只有黑白 2 种灰度。如果增加到 $k$ 种灰度，那只需要修改：

\[H(m,n)= \begin{cases} \lfloor \frac{t}{ {\rm maxval} }\cdot k\rfloor & t< {\rm maxval}\\ k-1 & t\geq {\rm maxval} \end{cases}\\ E_g(m,n)= \begin{cases} t-{\rm maxval}\cdot \frac{H(m,n)}{k} & H(m,n)=0..k-2\\ t-{\rm maxval} & H(m,n)=k-1 \end{cases}\]

参考

3.1 “Personal Computer Based Image Processing with Halftoning,” John A Saghri, Hsieh S. Hou, Andrew F. Tescher, Optical Engineering, March 1986, Vol. 25, No. 3, pp. 499-504. 论文链接

扩展资料

上篇图片输入输出

下篇直方图与均衡化