Neugebauer Models for Color Error Diffusion Halftoning

In this paper, we propose a method for halftoning color images based on an error diffusion technique, a color design criterion and Neugebauer models for expressing colors. For a natural extension of the conventional method for grayscale error diffusion to its color version, we first reformulate grayscale error diffusion with a one-dimensional Neugebauer model. Then we increase the dimension of the model to derive a color error diffusion method based on a three-dimensional Neugebauer model in RGB (red, green and blue) color space. Moreover, we propose a sparse Neugebauer model based on a color design criterion, or the minimal brightness variation criterion (MBVC), from which we derive a sparse Neugebauer model-based error diffusion method. Experimental results show that color halftone images produced by the proposed methods preserve the color contents in original continuous-tone images better than that by conventional color error diffusion methods. We also demonstrate that the proposed sparse method reduce halftone noise better than the state-of-the-art method based on MBVC.


Introduction
Digital halftoning is the process of converting a continuous-tone image into a pattern of a limited number of gray levels or colors. Recently, the techniques for digital halftoning have attracted the interest of researchers in the field of smart paper technology [1]. Smart paper is another name for the e-paper or electronic paper which is the electronic ink display technology designed to mimic the appearance of ordinary ink on paper. Electronic papers reflect light like ordinary papers, and can hold texts and images indefinitely without drawing electricity, while allowing the image to be changed later [2]. Therefore, it is said to be used for applications such as e-books, electronic newspapers, portable sings and foldable display [3]. Hsu et al. [4] mentioned the movie "Minority Report", in which an electronic newspaper appears, and remarked that those technologies in the near future are becoming reality now. Kao et al. [5] proposed an integrated driving system for playing videos/animations on the electrophoretic display (EPD). The EPD system contains a real-time signal processing hardware engine for video halftoning based on Jarvis error diffusion. The proposed algorithm in this paper is a kind of error diffusion algorithms such as Jarvis error diffusion; therefore, the proposed algorithm can be applied to the EPD system by replacing the conventional error diffusion algorithm with the proposed one.
Error diffusion [6] is one of the typical digital halftoning techniques. Floyd and Steinberg [7] originated a grayscale error diffusion algorithm. A simple extension of grayscale error diffusion to color is to apply grayscale error diffusion to each color channel in a color image separately. These kinds of extensions are referred to as separable methods [8]. Shaked et al. formulated a color design criterion named the minimal brightness variation criterion (MBVC), and developed an ink relocation post-process applicable to arbitrary color halftone algorithms [9] and a color diffusion algorithm which decreases halftone noise successfully [10].
Neugebauer model for color representation can be used to infer the gamut boundaries of printed colorants on paper [11] by printers such as conventional CMYK (cyan, magenta, yellow and black) printers. Monga et al. [12] proposed an extension of tone-dependent error diffusion to color, which employs a Neugebauer printer model with a simple correction by Yule and Nielsen [13]. Inoue et al. [14] applied the Neugebauer model to the estimation of spectral reflectance.
In this paper, we propose a grayscale Neugebauer model for representing an arbitrary grayscale value by a linear combination of black and white, and apply it to grayscale error diffusion to derive a grayscale Neugebauer model-based method. Then we show the equivalence between conventional grayscale error diffusion and the Neugebauer model-based method. After that, we extend the grayscale Neugebauer model-based error diffusion method to its color version by extending the dimension of color space from 1D to 3D. Therefore, the proposed color error diffusion method can be expected as a good candidate for the color extensions of a standard grayscale error diffusion technique. Moreover, we propose a sparse Neugebauer model for improving the quality of halftone images. Experimental results show the effectiveness of the proposed color error diffusion method on real image data including standard images for image processing researches.
The rest of this paper is organized as follows: Section 2 summarizes grayscale error diffusion, and shows the equivalence between conventional grayscale error diffusion and Neugebauer model-based method. Section 3 extends grayscale Neugebauer model to its color version. Section 4 proposes a color Neugebauer model-based error diffusion method. Section 5 proposes a sparse color Neugebauer model-based error diffusion method. Section 6 shows experimental results for demonstrating the performance of the proposed methods. Finally, Section 7 concludes this paper.

Grayscale Error Diffusion
Let F = [ f ij ] be a grayscale image, where f ij ∈ [0, 1] for i = 1, 2, . . . , m and j = 1, 2, . . . , n denotes the pixel value at the position (i, j) in F. Then conventional error diffusion method is described as follows.

Conventional Method
ij ] be a halftoned binary image of F. Then each pixel value h F ij in H F is given by in raster scan order, where θ denotes a threshold value for binarization, and is given as θ = 1 2 in this paper. This binarization in Equation (1) causes an error as follows: which is diffused to subsequent adjoining pixels as follows: where such a pixel f i+k,j+l at the position (i + k, j + l) is updated to f i+k,j+l by adding the error e ij weighted by predetermined coefficients w kl as shown in Figure 1, where "#" denotes the current pixel being processed, and "-" denotes the past pixel in the raster scan order. For subsequent pixels, the updated values f ij are binarized with Equation (1) instead of the original f ij .

Grayscale Neugebauer Model-Based Method
We present a grayscale Neugebauer model as follows: where K = 0 and W = 1 denote black and white, respectively, and the coefficients defined by can be viewed as the occurrence probabilities of black and white because of p K ij ≥ 0, p W ij ≥ 0 and p K ij + p W ij = 1. In other words, the grayscale Neugebauer model expresses an arbitrary grayscale value f ij as an internally dividing point between black (K) and white (W) as illustrated in Figure 2. Based on the grayscale Neugebauer model in Equation (4), we propose another expression of grayscale error diffusion as follows (the equivalence of these expressions will be discussed in the next subsection): A pixel value h F ij in H F is given by selecting black or white which has larger probability than another as After this binarization, the occurrence probabilities becomẽ using which, we define the errors in occurrence probabilities p K ij and p W ij as which are diffused to the subsequent pixels as follows: where the occurrence probabilities p K i+k,j+l and p W i+k,j+l of the subsequent adjoining pixel (i + k, j + l) are updated to p K i+k,j+l and p W i+k,j+l , respectively. For subsequent pixels, the updated values p K ij and p W ij are used for the binarization in Equation (7) instead of the original p K ij and p W ij .

Equivalence of the Two Methods
Here we show an equivalence of the above two methods for grayscale error diffusion. (1). On the other hand, when h F ij = W = 1 in Equation (7), (1). Therefore, the binarization in Equation (7) is equivalent to that in Equation (1).
Next, we describe the equivalence of error diffusion procedures. Assume that h F ij = K = 0 in Equation (7). Then the occurrence probabilities in Equations (8) and (9) (10) and (11)

and the errors in Equations
respectively. Substituting these errors into Equations (12) and (13), we have from both of which, using the relationships in Equations (5) and (6), we have which coincides with Equation (3) with Equation (2) for h F ij = K = 0. On the other hand, when h F ij = W = 1 in Equation (7), the occurrence probabilities in Equations (8) and (9) (10) and (11) become

and the errors in Equations
respectively. Substituting these errors into Equations (12) and (13), we have from both of which, using the relationships in Equations (5) and (6), we have which also coincides with Equation (3) with Equation (2) for h F ij = W = 1.
Consequently, the grayscale Neugebauer model-based error diffusion in Equations (12) and (13) can be reduced to the conventional error diffusion in Equation (3) with Equation (2), that shows the equivalence of the two methods.

Extending Grayscale Neugebauer Model to Higher Dimensions
The grayscale Neugebauer model illustrated in Figure 2 has a one-dimensional structure, where a midpoint f ij is expressed as a linear interpolant between K and W. This one-dimensional structure of grayscale Neugebauer model can be extended to two-dimensional bilinear interpolation as follows: Let f ij = (r ij , g ij ) be a two-dimensional color vector on an RG (red and green) color plane as illustrated in Figure 3. Then f ij can be expressed as a bilinear interpolant of four basic colors: K (black), R (red), G (green) and Y (yellow) as shown in Figure 3, where a linear interpolant between K and G is expressed as g ij G + (1 − g ij )K for a given green value g ij . Similarly, a linear interpolant between R and Y is expressed as g ij Y + (1 − g ij )R. Using those expressions, f ij on the RG color plane can be expressed as a linear interpolant of g ij G + (1 − g ij )K and g ij Y + (1 − g ij )R for a given red value r ij as follows: where four basic colors are given by K = (0, 0), R = (1, 0), G = (0, 1) and Y = (1, 1), and the coefficients are given by Thus we have a two-dimensional Neugebauer model in Equation (20) with Equations (21)-(24). Next, we further extend the above two-dimensional Neugebauer model to three-dimensional one in order to deal with arbitrary colors in three-dimensional color space. Figure 4 illustrates the derivation of three-dimensional Neugebauer model as a trilinear interpolation of eight basic colors in RGB color space for representing an arbitrary color as a linear combination of the basic colors. Let f ij = (r ij , g ij , b ij ) be a three-dimensional color vector in an RGB color space as shown in Figure 4. Then the orthogonal projection of f ij onto the bottom plane KRYG of the RGB color cube is given by Equation (20) using two-dimensional Neugebauer model. Similarly, another orthogonal projection of f ij onto the top plane BMWC of the RGB color cube is given by ( On the contrary, f ij can be represented as an interpolant of the above two orthogonal projections of f ij as follows: where eight basic colors are given by K = (0, 0, 0), R = (1, 0, 0), G = (0, 1, 0), Y = (1, 1, 0), B = (0, 0, 1), M = (1, 0, 1), C = (0, 1, 1) and W = (1, 1, 1), and the coefficients are given by This three-dimensional Neugebauer model in Equation (25) with Equations (26)-(33) is used for color error diffusion in the next section.

Color Neugebauer Model-Based Method
Analogous to the grayscale Neugebauer model-based method in Section 2.2, we propose a color error diffusion method based on the color Neugebauer model in Section 3.
In color image halftoning, it is common to use a limited number of basic colors instead of true color. Let S = {K, R, G, Y, B, M, C, W} be the set of eight basic colors. Then the color quantization of a pixel color f ij of an RGB color image F to one of the eight basic colors can be described as which is the corresponding pixel color in the color halftone image H F = [h F ij ] of F, and p x ij for x ∈ S are given by Equations (26)-(33). This expression in Equation (34) is a natural extension of grayscale binarization in Equation (7) to color quantization.
After this quantization, the occurrence probabilities becomẽ using which, we define the errors in occurrence probabilities p x ij for x ∈ S as which are diffused to the subsequent pixels as follows: where the occurrence probability p x i+k,j+l of the subsequent adjoining pixel (i + k, j + l) is updated to p x i+k,j+l . For subsequent pixels, the updated values p x ij for x ∈ S are used for the quantization in Equation (34) instead of the original p x ij .
The above method is summarized in Algorithm 1.

Sparse Color Neugebauer Model-Based Method
The above Neugebauer model expresses an arbitrary color as a linear combination of eight basic colors in S. However, this expression is somewhat redundant, because a color in RGB color space exists in a tetrahedron, the vertices of which are selected from S, and therefore, four basic colors are sufficient for expressing an arbitrary color as a linear combination of basic colors.
Shaked et al. [10] proposed a high-quality color error diffusion method based on the MBVC, i.e., to reduce halftone noise, select from within all halftone sets by which the desired color may be rendered, the one whose brightness variation is minimal [9,10].  Figure 5. Although Shaked et al. [9] suggested a different brightness order, K, B, R, M, G, C, Y and W, for a CRT screen, we would like to use the order in Figure 5 for consistency with Shaked's algorithm. That is, selecting a suitable brightness order is a problem-dependent on the display environment of colors.  = (r, g, b). Then this situation can be illustrated as Figure 6, and we would like to express f as a linear combination of S f as: where p 1 , p 2 , p 3 and p 4 are nonnegative coefficients satisfying p 1 + p 2 + p 3 + p 4 = 1. As shown in Figure 6, the vector f − s 1 can be expressed as a linear combination of the vectors s 2 − s 1 , s 3 − s 1 and s 4 − s 1 as follows: Solving this simultaneous equations for p 2 , p 3 and p 4 , we have the following expression of f : Putting p 1 = 1 − p 2 − p 3 − p 4 , we have a sparse Neugebauer model as follows: where q η = 0, and t η ∈ S \ S f for η = 1, 2, 3, 4, i.e., every t η is the element of the complement of S f in S. By replacing the standard Neugebauer model in Algorithm 1 with the above sparse Neugebauer model, we propose another color error diffusion algorithm as Algorithm 2, where the function MBVQ( f ) computes an MBVQ for a given color f by Shaked's algorithm.

Experimental Results
In this section, we demonstrate the performance of the proposed error diffusion methods using real images. Figure 7 shows an example of color image halftoning, where an original color image in Figure 7a is halftoned by conventional separable error diffusion method as shown in Figure 7b. The separable method [8] applies grayscale error diffusion to each color channel separately, and combines the resultant three grayscale halftone images into a single color halftone image. Figure 7c shows the result of the proposed color Neugebauer model-based error diffusion method summarized in Algorithm 1. Although Figure 7b,c look similar to each other from a distance, a close look at these images reveals that Algorithm 1 loses the clarity of the reproduction of contours. Algorithm 2 will improve the quality of contour reproduction as shown in the following results. In order to show the difference between Figure 7b,c clearly, we show the occurrence probabilities of eight basic colors (K = (0, 0, 0), R = (1, 0, 0), G = (0, 1, 0), Y = (1, 1, 0), B = (0, 0, 1), M = (1, 0, 1), C = (0, 1, 1) and W = (1, 1, 1)) in each image in Figure 8, where the vertical and horizontal axes denote the occurrence probability and the basic colors, respectively, where the occurrence probability of a basic color x ∈ S is defined by Blue, orange and green bars denote the original image and the halftoned images by the separable and the proposed methods, respectively. This figure shows that the occurrence probabilities for the separable method deviate from the original ones, on the other hand, that for the proposed method follow the original ones well. To confirm this difference quantitatively, we compute the absolute error of the occurrence probabilities defined by for x ∈ S, where p x Orig and p x Half denote the occurrence probabilities of basic color x in an original and its halftoned images, respectively. Figure 9 shows E(p x Orig , p x Half ) for eight basic colors, where orange and green bars denote that the occurrence probabilities p x Half are computed from the halftoned images by the separable and proposed methods, respectively. This figure shows that the proposed method can decrease the absolute error compared with the separable method.   Figure 10 shows additional results on the standard image database SIDBA [18], which contains twelve color images as shown in Figure 10a. Figures 10b,c show the results of the separable and the proposed methods, respectively. Figure 11 shows the mean absolute error of the occurrence probabilities of the eight basic colors, where the vertical and horizontal axes denote the mean absolute error and the twelve images shown in Figure 10, respectively. The proposed method denoted by green bars achieved smaller errors compared with the separable method denoted by orange bars for all twelve SIDBA images.  Next, we compare the proposed sparse Neugebauer model-based method (Algorithm 2) with Shaked's method [10]. Figure 12 compares conventional separable, Shaked's and the proposed sparse methods, where Figure 12a shows an original color image, which is halftoned by the separable, Shaked's and proposed methods as Figure 12b-d, respectively. Their zoomed parts are also shown at the bottom of them. As shown in Figure 12f, separable method places dark and light pixels side by side with each other here and there, which conflicts with the MBVC. On the other hand, Shaked's method reduces the brightness variation, however, the subtle gradations on the surface of ceramic ware are disturbed in Figure 12g. Compared with Shaked's method, the proposed method achieves both the brightness variation reduction and the subtle gradation reproduction as shown in Figure 12h.
For evaluating the quality of halftoned images in Figure 12 quantitatively, we show the occurrence probabilities of eight basic colors in each image in Figure 13, where the vertical and horizontal axes denote the occurrence probability and the basic colors, respectively. Blue, orange and green bars denote the original image and the halftoned images by Shaked's and the proposed methods, respectively. This figure shows that the occurrence probabilities for the Shaked method deviate from the original ones, on the other hand, that for the proposed sparse method follows the original ones well. To confirm this difference quantitatively, we show the absolute error of the occurrence probabilities in Figure 14, which shows that the proposed sparse method can decrease the absolute error compared with Shaked's method. Figure 15 shows halftoned SIDBA images by Shaked's and the proposed sparse methods, which produce visually similar halftone images to each other.
However, the produced halftone images are not identical as shown in Figure 16, where Figure 16a-d show the zoomed parts of the same subregions of Milkdrop images in Figure 10a,b, and Figure 15a,b, respectively. The original patch in Figure 16a is halftoned by the separable method as shown in Figure 16b, where all eight basic colors are used to reproduce the original colors. In Figure 16c, Shaked's method succeeds in reducing black and white pixels for minimizing the brightness variation. However, in Figure 16c, the second darkest blue and the second lightest yellow pixels are still intermingled spatially. On the other hand, in Figure 16d, the proposed sparse method reduces the spatial mixture of blue and yellow pixels compared with Figure 16c. This result demonstrates that the proposed sparse method (Algorithm 2) can separate different colors spatially on the image plane better than Shaked's method, which improves the clarity of the reproduction of contours, and therefore, Algorithm 2 can also be applicable to the reproduction of textual information as well as Shaked's method.
To confirm the difference between Shaked's and the proposed sparse methods quantitatively, we show the mean absolute error of the occurrence probabilities of the eight basic colors in Figure 17, where the proposed sparse method denoted by green bars achieved smaller errors compared with the Shaked's method denoted by orange bars for all twelve SIDBA images.

Conclusions
In this paper, we proposed a method for color error diffusion by extending the grayscale Neugebauer model to its color version. The color Neugebauer model represents a color by a linear combination of basic colors, where the coefficients can be interpreted as the occurrence probabilities of the basic colors because they satisfy the probability axioms.
We also proposed a sparse Neugebauer model for color error diffusion halftoning. The proposed sparse method is based on the MBVC proposed by Shaked et al. Therefore, the proposed sparse method can also reduce halftone noise as well as Shaked's method.
Experimental results demonstrated that the color halftone images given by the proposed method preserve the occurrence probabilities of the basic colors of the original images well. Compared with Shaked's method, the proposed sparse method achieved a better performance both visually and quantitatively.