{'fraction_non_alphanumeric': 0.04928048878527213, 'fraction_numerical': 0.011603290725406652, 'mean_word_length': 4.277612784613209, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Statistical mechanics of digital halftoning

We consider the problem of digital halftoning from the view point of statistical mechanics. The digital halftoning is a sort of image processing, namely, representing each grayscale in terms of black and white binary dots. The digital halftoning is achieved by making use of the threshold mask, namely, for each pixel, the halftoned binary pixel is determined as black if the original grayscale pixel is greater than or equal to the mask value and is determined as white vice versa. To determine the optimal value of the mask on each pixel for a given original grayscale image, we first assume that the human-eyes might recognize the black and white binary halftoned image as the corresponding grayscale one by linear filters. The Hamiltonian is constructed as a distance between the original and the recognized images which is written in terms of the threshold mask. We are confirmed that the system described by the Hamiltonian is regarded as a kind of antiferromagnetic Ising model with quenched disorders. By searching the ground state of the Hamiltonian, we obtain the optimal threshold mask and the resulting halftoned binary dots simultaneously. From the power-spectrum analysis, we find that the binary dots image is physiologically plausible from the view point of human-eyes modulation properties. We also propose a theoretical framework to investigate statistical performance of inverse digital halftoning, that is, the inverse process of halftoning. From the Bayesian inference view point, we rigorously show that the Bayes-optimal inverse-halftoning is achieved on a specific condition which is very similar to the so-called Nishimori line in the research field of spin glasses.


I. INTRODUCTION
Recently, a lot of problems of information science and technology have been investigated by several useful tools developed in the research field of statistical mechanics of spin glasses [1][2][3][4]. Statistical mechanics of information is now widely spreading in various subjects such as neural networks [3], error-correcting codes [5][6][7], CDMA multi-user demodulator [8], image processing, etc [3,9]. Especially, image restoration by making use of a graphical model of Markov random fields has been investigated extensively from both analytical and numerical point of views [9]. In the research field of image processing, digital halftoning instead of image restoration, which is defined as a process of generating a pattern of pixels with limited number of colors, especially converting a grayscale image into the binary black and white picture, has been widely used in various practical situations in media, such as the printing of newspapers, fax machines and so forth [10,11].
To achieve the digital halftoning, one needs the strategy to arrange the geometrical-combination of black and white pixels so as to make human-eyes to have a kind of optical illusion. Namely, the halftoning relies on the fact that the human-eyes act as a spatial low-pass filter and can not recognize any specified structure in the part of image dominated by high frequency components.
We plot the shape of these contrast sensitivity functions (1) and (2)  high-frequency regime f ρ ≫ 1, the both contrast sensitivity functions decrease to zero, which means that human eyes cannot recognize high-frequency components in images. In other words, the halftone algorithm should be constructed such that the halftoned binary dots contain relatively high-frequency components to describe the original grayscale revels. This fact is an important guide for us to consider halftoning algorithms from the view point of physiology. Up to now, to achieve fine qualities of the digital halftoning, a lot of techniques, for instance, clustered-dot ordered dither method, threshold mask method [10,11] and blue noise mask method [14,15], error-diffusion method [16], etc. have been proposed and developed by many researchers. However, no attempt has been done to formulate the important problem from the view point of statistical mechanics of information. It is not difficult for us to assume that each pixel for the halftone and the grayscale images is represented by Ising spin and Q-Ising spin (or Potts spin), respectively. Therefore, it seems that statistical mechanical approach is very useful for the problem.
On the other hand, to retrieve the original grayscale image from the halftone binary dots, the so-called inversehalftoning has been also developed especially in scanner technology [17,18]. As such an inverse halftoning strategy, the conventional smoothing filters have been used widely. However, it is worth while for us to look for the alternative or reconsider the problem from the view point of statistical mechanics of disordered spin systems.
In this paper, we propose a statistical-mechanical modeling of both digital halftoning and the inverse process, namely, inverse-halftoning. To attempt to generate the halftoned binary dots, we first introduce a threshold mask matrix with the same size as the original image. For the threshold mask, each pixel of the halftoned image is determined as a black pixel if the value of the original grayscale pixel is greater than the component of the threshold mask matrix and is determined as a white pixel vice versa. Then, we assume that human eyes might recognize the original grayscale image by a linear filter. By taking into account that the distance between the original and the recognized images is written in terms of the threshold mask, we naturally introduce a Hamiltonian to be minimized to determine the optimal threshold mask as the ground state. We explicitly show that the system described by the Hamiltonian is a variant of the antiferromagnetic Ising model with disorders. For a demonstration of our method to generate the optimal threshold mask, computer simulations of simulated annealing [19,20] is carried out. To investigate the statistics of the threshold mask image and the halftoned image, we evaluate the power-spectrum numerically.
On the other hand, in order to investigate statistical performance of the inverse digital halftoning, in our previous study [17], we introduced the standard regularization theory for a Markov random-fields model to represent original grayscale image and construct the inverse process as a kind of image reconstruction of the grayscale image from a given halftone binary dots. Then, we evaluated the statistical performance by making use of Markov chain Monte Carlo (MCMC) simulations and analysis of the infinite-range mean-field model. We also investigated the inverse-halftoning process as a dynamics of disordered spin systems. In this paper, we reconsider the inverse halftoning problem from the Bayesian inference point of view. We show that the Bayes-optimal inverse-halftoning that minimizes the mean square error is achieved on a specific condition which is similar to the so-called Nishimori line [21] in the research field of spin glasses. Finally we discuss the possibility that both halftoning and the inverse-halftoning processes are unified under a single Hamiltonian, namely, we argue whether it is possible for us to obtain the threshold mask, the halftoned and inverse-halftoned images simultaneously by finding the ground state of the spin systems or not.
This paper is organized as follows. In the next Sec. II, we explain our procedure to generate the optimal threshold mask for a given original grayscale image. In this section, the Hamiltonian is introduced as a function of the threshold mask under the assumption that the human-eyes might recognize the halftoned image as the corresponding original grayscale version by a linear filter. In Sec. III, we solve the combinatorial optimization problem described by the Hamiltonian by making use of the simulated annealing. We show the resulting halftoned image for a standard grayscale image. We find that our algorithm generates the binary dots which induce the illusion of a continuous-tone image. We also investigate the power-spectrum statistics of the threshold mask array and the halftoned image. In Sec.IV, we formulate our inverse-halftoning procedure based on the MPM estimation of Bayesian statistics. We show clearly that the Bayes-optimal inverse-halftoning that minimizes the mean-square error is achieved on a specific condition. In Sec. V, we discuss the possibility that both halftoning and the inverse-halftoning processes are unified under a single Hamiltonian. The last section is summary.

II. STATISTICAL-MECHANICAL MODELING OF HALFTONING
We first deal with the 'forward problem' of the digital halftoning. In this section, we first explain the procedure to generate the binary dots to represent the original grayscale image via what we call threshold dither method [10,11]. Then, we explain how one formulates the halftoning processes as a problem of disordered spin systems and why statistical mechanics of information is useful for the problem.

A. Halftoning by threshold dither method
Let us first define the original grayscale image located on the square lattice with size L 1 × L 2 by g = {g x,y ∈ 0, 1, · · · , Q − 1|x = 1, 2, · · · , L 1 , y = 1, 2, · · · , L 2 } . ( To generate the binary dots to represent the grayscale image g, we introduce the threshold mask matrix with the size l 1 × l 2 (l 1 ≤ L 1 , l 2 ≤ L 2 ): t = {t x,y ∈ 0, 1, · · · , Q − 1|x = 1, 2, · · · , l 1 , y = 1, 2, · · · , l 2 } We should keep in mind that each component of the matrix t x,y should be satisfied the following periodic boundary condition: Then, the halftoned image is defined by h = {h x,y ∈ 0, 1|x = 1, 2, · · · , L 1 , y = 1, 2, · · · , L 2 } and each pixel in the h is calculated as where the step function Θ(x) is defined conventionally as For this setup of the halftoning procedure, a pixel of the resulting halftoned binary image is set to one if the pixel of the original grayscale image is greater than or equal to the corresponding pixel of the threshold mask; otherwise the pixel is set to zero. For instance, the digital halftoning is achieved by threshold mask of Bayers' type for Q = 16 with size l 1 × l 2 = 4 × 4: Namely, for a part with size 4 × 4 of the original image g, say, for    g 1,1 g 1,2 g 1,3 g 1,4 g 2,1 g 2,2 g 2,3 g 2,4 g 3,1 g 3,2 g 3,3 g 3,4 g 4,1 g 4,2 g 4,3 g 4,4    =    we have the corresponding block of the halftone image h as In FIG. 2, we show a typical example of the original grayscale image (256 grayscale levels: left panel), the image with reduced grayscales Q = 16: center panel) and the resulting halftone image obtained by the Bayer's threshold mask (9). The size of the image is L 1 × L 2 = 256 × 256 and the number of the grayscale levels is 256. Here we first converted the original 256-grayscales g ′ image to the reduced |N(x, y)| = Q(< 256)-grayscale image g by the following transform: and used ∀ x,y g x,y = C(g Therefore, once we obtain the appropriately threshold matrix t, the halftone binary dots h is determined uniquely, and for this reason, the quality of the halftoning is dependent on how we choose the threshold mask appropriately. In this paper, we consider the case of the threshold mask having the same size as that of the original grayscale image, that is, the case of l 1 = L 1 , l 2 = L 2 , and propose an algorithm to obtain the optimal threshold mask as lowest energy states of the Hamiltonian which is a square distance between the original and the recognized images. We shall introduce it in the next subsection.

B. A recognition model of human vision
As we explained in the previous subsection, our main problem is now to determine the threshold mask t efficiently. For this purpose, we might assume that human eyes might retrieve the grayscale image, let us call, the recognized image: s = {s x,y ∈ 0, 1, · · · , Q − 1|x = 1, 2, · · · , L 1 , y = 1, 2, · · · , L 2 } from the halftoned image h by making use of the following linear filter: where N(x, y) denotes the nearest neighboring pixels of the pixel located at (x, y) and the point (x, y) itself. Apparently, for a two-dimensional square lattice, the ingredients of N(x, y) are (x, y), (x − 1, y + 1), (x, y + 1), (x + 1, y + 1), (x − 1, y), (x + 1, y), (x − 1, y − 1), (x, y − 1), (x + 1, y − 1) and the size is |N(x, y)| = 9. The above choice for the recognition model comes from our assumption that human-eyes might recognize a part of halftone pictures more dark when the density of black pixels located on the part is relatively high. The justification of this simple assumption could be checked by comparing the recognized image s, which is calculated in terms of the resultant halftone binary dots via (14), with the original grayscale image g. We shall discuss the results in Sec. III. Obviously, for the simplest choice, we might set the weight W x,y appearing in (14) as ∀ x,y W x,y = 1. We should notice that for this choice, the spatial structure of the black dots in N(x, y) is not taken into account and the number |N(x, y)| itself determines the grayscale level of the corresponding pixel s x,y in the recognized image s. Therefore, the choice of the effective grayscale level can be controlled by choosing the number of the nearest neighboring pixels |N(x, y)|. Taking into account this limitation, we first reduce 256-grayscale levels to Q(< 256) by (12) and regard the Q-grayscale image as the g.
We also might utilize the other choice such as the following two-dimensional Gaussian-type: For this choice, not only the number |N(x, y)| but also the spatial structure determines the grayscale levels in the recognized image. However, we naturally assume that the spatial structure does not affect the quality of halftoning if the size of the window |N(x, y)| is small enough in comparison with the image size L 2 .
In this paper, we use the following specific choice of the weight W x−i,x−j with size |N(x, y)| = 5 × 5 = 25 as We should keep in mind that each recognized pixel s x,y takes the minimum s x,y = 0 and the maximum s x,y = 255 from the definition (14). Hence, each pixel in the recognized image s calculated by (14) represents 'defective' 256-grayscale levels, that is, ∀ x,y s x,y ∈ {0, 5, 10, · · · , 255}. Then, we might choose a strategy to determine the threshold mask t that minimizes the square distance between the original and the recognized images. Namely, we minimize the energy g − s 2 ≡ L1 x=1 L2 y=1 (g x,y − s x,y ) 2 . We should notice that from the relation (14), the above distance g − s 2 is written in terms of the threshold mask t for a given original image g. By substituting (14) into the distance g − s 2 , we obtain the Hamiltonian of the system as a function of the threshold mask t. Thus, our digital halftoning is now reduced to a combinatorial optimization problem of the following Hamiltonian: The above Hamiltonian H(t|g) is a starting point of statistical-mechanical modeling of the digital halftoning. We should bear in mind that in the above Hamiltonian, t are dynamical variables, whereas, g are quenched disorders.
As we have the Hamiltonian of the system, one can generate statistical ensembles according to standard statistical mechanics. In its thermal equilibrium at temperature β −1 , each possible microscopic state t obeys the following Boltzmann-Gibbs distribution: One of the simplest uses of the Hamiltonian is to regard its minimum energy state as the optimal threshold mask t ≡ {t x,y ∈ 0, 1|x = 1, 2, · · · , L 1 , y = 1, 2, · · · , L 2 }. Namely, we havê which is rewritten in terms of the statistical-mechanical terminology as Calculating the Q L1×L2 -sums: tx,y=0 (· · · ) and taking the zero-temperature limit are performed by Gibbs sampler with temperature annealing procedure during the Monte Carlo steps. The solution might be a good candidate for the suitable threshold mask.
We should mention that our procedure is categorized in the method 'model-based halftoning' which means that the threshold mask is dependent on the original image. In engineering, there exist several such algorithms [22][23][24]. However, we should stress that these studies are completely different from our statistical-mechanical modeling.

C. The corresponding spin system
It is helpful for us to consider the meaning of the Hamiltonian (17) from the view point of disordered spin systems. To see the relationship between halftoning processes and the corresponding disordered spin system, we expand the square of the Hamiltonian (17). Then, we have It should be noted that we canceled a constant term L1 x=1 L2 y=1 (g x,y ) 2 and used the relation Θ(x) = (1 + sgn(x))/2. Thus, our system is nothing but an Ising model defined on the two-dimensional square lattice with antiferromagnetic interactions with random field g on pixel. However, the ground state which minimizes the above Hamiltonian with respect to t is complicated due to the quenched variables as a given original image g appearing in the argument of sign function such as sgn(g ij − t ij ). The lowest energy state in the dynamical variable space of t might be highly degenerated and it means that we can use various techniques developed in the field of statistical physics of disordered spin systems to obtain the lowest energy state and to investigate the lowest energy properties. In the next section, we show some demonstrations to make the halftoned binary dots that imitate the original grayscale image by minimizing the Hamiltonian (17) via simulated annealing.

III. NUMERICAL EXPERIMENTS
In this section, we show some demonstrations to generate the optimal threshold mask and the resulting halftoned image for a given standard grayscale image. We shall demonstrate our algorithm for one of the well-known standard images shown in FIG.4 (upper left). For the original image g with the size L 1 × L 2 = 400 × 400 as input data, we construct the Hamiltonian H(t|g) and carry out the calculation of the threshold maskt given by (20) via simulated annealing with the temperature schedule β ∼ √ t. For the solution of the threshold maskt, we make the halftoned image h according to (6).
In In the lower right panel of FIG. 4, we show the recognized image s. We find that the recognized image is quite similar to the original grayscale image and this result is a justification for our assumption to construct the recognition model (14). From these pictures, we might conclude that our algorithm works very well and the resulting halftoned binary dots look fine to represent the original grayscale levels. However, it is quite important for us to evaluate the performance quantitatively. In following , we evaluate the performance of the halftoning quantitatively.
In FIG. 6 (right), we show the histograms of grayscale levels for the original and recognition images. We find that the both shapes of the histograms have a similarity, however, there is a gap between the histograms in their height. Apparently, if the perfect minimization of the Hamiltonian is achieved, these two histogram should coincide with each other. Some theoretical argument on this issue is given in Appendix A.
As we mentioned, this gap comes from the fact that the grayscale levels in the recognition image are restricted (defected) to s x,y = 0, 5, 10, 15, · · · , 255 due to the definition of the weight W x−i,y−j with size |N(x, y)| ≪ 256. We discuss this issue later. In the right panel of FIG. 6,   From the view point of human visual systems, we check the power-spectrum of the threshold mask and the halftone binary dots. In order to investigate the statistical properties of the model, we next investigate the power spectrum of the threshold mask and the halftone binary dots. Namely, we evaluate where we should replace o i,j by the threshold mask t i,j or the halftoned binary dots h i,j . We should bear in mind that here we consider the case L 1 = L 2 = L.
We show the results in FIG. 7. As well-known, the so-called blue noise mask, there exists principal frequency and below the value, the power spectrum drops to zero, whereas, the high-frequency components remains finite. On the other hand, the threshold mask generated by our algorithm apparently depends on the structure of the original grayscale image. Therefore, it is assumed that some of the spatial structures of the original image, especially, the smoothness leading up to the low frequency components in the power spectrum remains finite. However, from the power spectrum P t (f x , f y ) shown in FIG.8, we find that the low-frequency components almost disappear and highfrequency components remain finite. More convenient spectrum statistics is the so-called radial-averaged power spectrum [10] defined by where R(f ρ ) means the region of annular rings in the Fourier space with the radiuses f ρ and f ρ + ∆ ρ , namely, the width of the annular rings is given by ∆ ρ . N (R(f ρ )) is the number of frequency samples in R(f ρ ). In FIG. 8, we plot the spectrum statistics for both threshold mask and resulting binary dots of the halftoned image. Inset of FIG. 8 shows the same statistics for the threshold mask image. From this figure, we find that the low-frequency components are not dominant and relatively high-frequency components are larger than the low-frequency counter part. This tendency should be plausible from the view point of human-eyes modulation as we explained in FIG. 1.

IV. INVERSE DIGITAL-HALFTONING
In the previous sections, we consider the digital halftoning from the view point of statistical mechanics of spin systems. In engineering perspective, it is also important for us to retrieve the original grayscale image from a given halftoned binary dots when we attempt to capture the halftone image via scanner machines. The inverse problem of the digital halftoning is referred to as inverse-halftoning. In this section, we investigate the statistical performance of the inverse-halftoning by making use of the concept of statistical mechanics of disordered spin systems. Especially, we discuss the condition on which the Bayes-optimal inverse-halftoning is achieved.

A. Bayesian formula
We first provide a Bayesian formula of inverse-halftoning. As we mentioned, the halftoned image h is given pixelwise as h x,y = Θ(g x,y − t x,y ), ∀ x,y . Thus, we might regard the process as a kind of 'channel' which is written in terms of the following conditional probability: for a given threshold mask t. Hence, the halftoning process is regarded as stationary memory-less deterministic erasure channel from the information theoretical point of view. When we choose the ferromagnetic prior, the posterior is obtained as P (σ|h : t) = P t (h|σ)P J (σ) where we defined (x,y),(k,l) (· · · ) as a sum over the nearest neighbouring pixel pairs. From the view point of statistical mechanics, the above posterior is regarded as the Boltzmann distribution of the Q-Ising model in which the mobility of each spin is tightly constrained in the range t x,y ≤ g x,y ≤ Q − 1 for h x,y = 1 and 0 ≤ g x,y < t x,y for h x,y = 0.
For the above posterior, the Maximizer of Posterior Marginal (MPM for short) estimate is given by pixel-wise optimization of the posterior aŝ σ x,y = arg max σx,y P (σ x,y |h : t) where Θ Q (x) stands for a function to convert a real value x to the nearest integer.
In the inverse halftoning, we should solve ∀ x,y h x,y = Θ(g x,y − t x,y ) with respect to g for a given halftoned image h and the mask t.
In the reference [25], one of the authors considered restoration processes of grayscale images by making use of bit-decomposed data. Namely, we generate the Q binary images whose pixel is given by ∀ x,y h (m) x,y = Θ(g i − m), m ∈ {0, · · · , Q − 1} and transmitting them through some noise channels. Then, the paper [25] dealt with the procedure to retrieve the original image g from the degraded Q slices of the binary images x,y |x = 1, · · · , L 1 , y = 1, · · · , L 2 } and η (m) ≡ {η (m) x,y |x = 1, · · · , L 1 , y = 1, · · · , L 2 } stands for the additive noise. If there is no degrading process, it is obvious that the vector is identical to the original grayscale image g. However, in the present inverse-halftoning case, the only information we have is just a single slice h (m) . This fact makes the problem hard to treat. Obviously, this inverse-halftoning is a typical ill-posed problem because there are a lot of candidates to satisfy the equations. In Appendix B, we evaluate two relevant quantities, namely, degree of degeneracy for possible solutions and mutual information to evaluate the difficulties of the problem.
In the next subsection, we discuss the relationship between the so-called Bayes-optimal solution and the Nishimori line established in the research field of spin glasses.

B. Bayes-optimal inverse-halftoning and the condition to achieve it
In the previous studies [18], we investigated the inverse-halftoning on the bases of Markov chain Monte Carlo simulations and analysis of the infinite-range mean-field model. Then, we found several conditions on the hyperparameters appearing in the Hamiltonian that gives a minimum of the mean square error numerically. However, so far, we do not yet provide any mathematically rigorous results on the performance of the inverse-halftoning defined on realistic two dimensional square lattices.
In this section, we attempt to prove that the Bayes-optimal inverse-halftoning is achieved on a specific condition which is similar to the so-called Nishimori line. For the purpose, let us use here an alternative definition of pixel index. Namely, for x, y = 1, 2, · · · , L, we change the index by means of i = x + Ly − L, i = 1, 2, · · · , L 2 ≡ N (conversely, x = (i mod L), y = [i/L]). Then, we consider the true prior and the true inverse-halftoning process as follows.
where we defined the original image g = {g i = 1, · · · , Q|i = 1, · · · , N } and the halftoned image h = {h i = 0, 1|i = 0, · · · , N }. We also used the definition ij to represent all nearest neighbouring pairs on the arbitrary lattice in finite dimension. The sums g (· · · ) and h (· · · ) denote g = i gi=1,··· ,Q (· · · ) and h = i hi=0,1 (· · · ), respectively. It should be noted that the above likelihood is another representation of the following dither method for each pixel For this original image and the halftone process given by the likelihood (29), we naturally use the following posterior: The quantity to be evaluated is the following mean square error D for an arbitrary pixel i: where E g,h [· · · ] ≡ P t (h|g)P J0 (g)(· · · )dgdh and we defined the average of a single pixel over the posterior by To proceed the proof, we should notice that for any stochastic variable A, the following fact is satisfied When we set A = g i − Θ Q ( σ i J ), we have the following inequality: In following, we shall evaluate the lower bound of the D, that is, the right most term of the above inequality. Among the staffs of the right most part in the above inequality, E g,h [g i ] is easily evaluated as where we defined and Therefore, E g,h [g i ] is identical to the local magnetization m On the other hand, the lest of the term appearing in the right hand side of the equation (34) is written by Therefore, the lower bound of the mean square error is evaluated as follows.
where the equality on the last line is satisfied for J = J 0 . We also used the fact that the sign of (m J0,g>t are monotonically increasing function with respect to J. δ = δ 1 , δ 2 is a quantization error due to the Q-generalized step function Θ Q (· · · ) defined as From the argument we presented above, we found that the performance of the inverse-halftoning achieved by the posterior (33) is optimized on the specific condition J = J 0 which is similar to the Nishimori line (point) [21] in the research field of spin glasses.

V. SIMULTANEOUS GENERATION OF MASK, HALFTONE AND RECOGNITION IMAGES
Finally we show that both halftoning and the inverse-halftoning processes are unified under a single Hamiltonian, namely, it is possible for us to obtain the threshold mask, the halftoned and inverse-halftoned images simultaneously by finding the ground state of the spin systems.
In this paper, we proposed the digital halftoning and the inverse-halftoning separately. However, from the form of the Hamiltonian H(t|g), the definition of the halftoned image: ∀ x,y h x,y = Θ(g x,y − t x,y ) and the recognition image: ∀ x,y s x,y = i,j∈N(x,y) W x−i,y−j Θ(g x,y − t x,y ), it is confirmed that these three important images are obtained simultaneously in the single theoretical framework, that is, minimizing the Hamiltonian H(t|g) with respect to the mask t for a given original image g.
Then, we should remember that there is a gap between the histograms of the grayscale levels for the recognition and the original images. This gap might make us hard to accept the recognition image as the inverse-halftoned one as a solution.
As we mentioned before, this gap comes from the fact that the grayscale levels in the recognition image are restricted (defected) to s x,y = 0, 5, 10, 15, · · · , 255 due to the definition of the weight W x−i,y−j with size |N(x, y)| = 5 × 5 = 25 ≪ 256 = Q. To reduce the gap, we might use the following linear filter: where the initial condition s (0) are chosen as ∀ x,y s x,y = i,j∈N(x,y) W x−i,y−j Θ(g x,y − t x,y ) for the minimum energy state of the Hamiltonian H(t|g). Thus, for a given s (0) , we recursively operate the above map (44) n times, and then, one might obtain more plausible image than the ∀ x,y s (0) x,y as the inverse-halftoned image. We show the resulting histogram of the grayscale levels in FIG. 9 and the corrected recognition images in FIG. 10 for n = 1 and 2. We set the number of nearest neighbouring pixels around the pixel at (x, y) as |N ′ (x, y)| = 3 × 3 = 9, namely, N ′ (x, y) = {(x−1, y +1), (x, y +1), (x+1, y +1), (x−1, y), (x, y), (x+1, y), (x−1, y −1), (x, y −1), (x+1, y −1)}. From these figures, we find that the gap between two images actually reduced and the resulting recognition image is improved by 'mixing effect' on the defected grayscale levels s x,y = 0, 5, 10, 15, · · · , 255. However, when we increase the number of iteration n, the mixing effect works too much on the recognition image and it makes the local structure of the image too smooth (see FIG. 9 (right) and FIG. 10 (lower right)).

VI. SUMMARY
In this paper, we proposed a possible statistical-mechanical modeling for digital halftoning. This formalism helps us to understand the problem as a combinatorial optimization which is described by a sort of disordered spin systems. Finding the ground state was achieved by simulated annealing and we found that the resulting binary dots looks well to represent the original grayscale levels. The quality of the binary dots was evaluated from the power-spectrum statistics. We found that the binary dots contain relatively high frequency components which are plausible from the view point of human-eyes modulation properties. We also proposed a theoretical framework to evaluate the statistical performance of the inverse digital-halftoning based on statistical mechanics. From the Bayesian inference view point, we rigorously show that the Bayes-optimal inverse-halftoning is achieved on a specific condition which is very similar to the so-called Nishimori line in the research field of spin glasses. We hope our formulation might be applied to generating of binary halftoned images effectively and evaluating the performance for the inverse-halftoning from halftoned binary images obtained by various algorithms. In Sec. III, we evaluated the performance of halftoning through two different measurements, namely, the histograms of grayscale levels and the power-spectrum via computer simulations. However, it might be helpful for us to evaluate the performance analytically. In this Appendix, we derive the distribution of the recognized pixels, that is, the number of black pixels in a window with |N(x, y)| = Q(≪ 256) analytically. Namely, we consider the distribution of the following quantities: where we consider the case of L 1 = L 2 = L and defined the index i so as to satisfy i = x + Ly − L, namely, i takes i = 1, 2, · · · , L 2 and we set L 2 = N . It is obvious that (x, y) for a given i are given as x = (i mod L), y = [i/L]. Here we consider the case of |N(i)| = Q = O(1). Then, we obtain the general formula for the distribution as follows. where H eff denotes the following effective Hamiltonian: In the above expression, we omitted the term i g 2 i which is independent of the dynamical variable t. We should keep in mind that the normalization ∞ −∞ dµ i P (µ i : g) = 1 is satisfied. It also should be noted that the pixels are fully connected in the window with size |N(i)| = Q. Then, the effective Hamiltonian is reduced to the decoupled form H eff = i H (i) eff and we immediately obtain by using the saddle point method in the limit of N → ∞. f is given by Then, equation (A2) leads to Obviously, the above distribution is dependent on the original image g and such data-averaged distribution is evaluated after slightly complicated algebra as follows.
As a demonstration, we choose uniform grayscale images having a single grayscale g 0 , namely, the distribution of G is given by P (g) = δ(g − g 0 ). Then, the data-average is easily performed as In FIG. 11, we plot the P β,Q (µ) for several values of Q and β and g. From this figure, we find that as temperature β −1 decreases, the distribution changes its shape to the delta function in which the peak is located at g 0 /Q. This means that simulated annealing safely finds the ground state of the Hamiltonian having the zero energy, namely, This immediately reads µ i = g 0 /Q for all index i. Therefore, for a given image having a single grayscale g 0 , the distribution should converge to the delta function with the peak at g 0 /Q if the annealing schedule of β is appropriate. Here we shall consider the inverse process of halftoning as an ill-posed problem. In the inverse halftoning, we should solve ∀ x,y h x,y = Θ(g x,y − t x,y ) with respect to g for a given halftoned image h and the mask t. Obviously, this problem is ill-posed because there are a lot of candidates to satisfy the equations. In this Appendix, we analytically evaluate several relevant quantities to show the difficulty in finding the solution.

Degree of degeneracy for possible solutions
We easily assume that the number of the solutions for a given h and t is exponential order, however, it is helpful for us to evaluate the number more precisely. For the purpose, let us use the definition introduced in the previous section, namely, i = x + Ly − L, i = 1, 2, · · · , L 2 ≡ N (conversely, x = (i mod L), y = [i/L]). Then, the number of the candidates for the solution of the equations ∀ i h i = Θ(g i − t i ) is evaluated for a given a realization of the original image and the threshold mask as follows.
As the number N seems to be an exponential order object, we might rewrite the data average of the logarithm of the number as a self-averaging quantity where we defined the expectation E g,t [· · · ] by Therefore, the above averages could be carried out for a specific choice of the P (g) and P (t). For distributions of the threshold mask, we suppose that in each mask with size Q × Q (Q ≪ N ), each component t i takes a value among 0, · · · , Q − 1 grayscales with equal probability 1/Q. Then, the P (t) is reduced to the product of the effective single site distribution as On the other hand, as a distribution of the original image, we consider snapshots from the infinite-range ferromagnetic Q-Ising model, that is, where Z 0 is a normalization constant for the probability P J0 (g). After simple algebra, we find that the P J0 (g) is rewritten such as P J0 (g) = function as J −1 0 decreases, however, the φ remains a finite value even if J −1 0 → ∞. This means that the number of the candidates for the solutions of the inverse halftoning is always exponential order. Therefore, we need some systematic approach to solve this type of the ill-posed problems. Namely, we need to introduce the ferromagnetic prior P J (σ) to compensate the lack of information as we discussed in Sec. IV.

Mutual information
We next consider the difficulties of retrieving original grayscale images from a slightly different point of view. Here we calculate the mutual information between the original image g and the halftone image h. From the definition of the mutual information, we should evaluate I(g : h) = H(h) − H(h|g), where we defined the entropy H(h) and the conditional entropy H(h|g) as H(h) = − h P (h) log P (h), H(h|g) = − g P J0 (g) h P t (h|g) log P t (h|g), respectively. It should be noted that we used P t (h|g) = i δ(h i , Θ(g i − t i )) (B8) and P (h) = g P t (h|g)P J0 (g).
For the infinite-range model and constant mask t i = t 0 ∀ i , these entropies per are calculated analytically as halftoning process by the dither method with the degraded process by a Gaussian noise with mean a 0 g and the variance a 2 , we also calculate the mutual information per pixel for the Gaussian noise. We choose the same distribution P J0 (g) as the halftone case. We immediately have In FIG. 13, we plot the mutual information per pixel for the Gaussian noise as a function of J 0 for various cases of the signal-noise ratio, a 0 /a.