Local Contrast-Based Pixel Ordering for Exact Histogram Specification

Histogram equalization is one of the basic image processing tasks for contrast enhancement, and its generalized version is histogram specification, which accepts arbitrary shapes of target histograms including uniform distributions for histogram equalization. It is well known that strictly ordered pixels in an image can be voted to any target histogram to achieve exact histogram specification. This paper proposes a method for ordering pixels in an image on the basis of the local contrast of each pixel, where a Gaussian filter without approximation is used to avoid the duplication of pixel values that disturbs the strict pixel ordering. The main idea of the proposed method is that the problem of pixel ordering is divided into small subproblems which can be solved separately, and then the results are merged into one sequence of all ordered pixels. Moreover, the proposed method is extended from grayscale images to color ones in a consistent manner. Experimental results show that the state-of-the-art histogram specification method occasionally produces false patterns, which are alleviated by the proposed method. Those results demonstrate the effectiveness of the proposed method for exact histogram specification.


Introduction
Histograms represent the distribution of pixel values in images, and can offer various image statistics [1] to us, i.e., histograms have useful information for image processing applications such as image enhancement, compression, and segmentation [2][3][4]. Therefore, histograms are the basis for numerous spatial domain processing techniques [5], which include intensity transformation as a special case where the smallest neighborhood of size 1 × 1 is used. Histogram equalization, which theoretically transforms an input image into the corresponding output image having a uniform histogram, is one of the intensity transformation techniques for contrast enhancement [5]. More generally, the method for generating an output image that has a specified target histogram is called histogram specification or histogram matching [5], in which histogram equalization is included as a special case [6]. It is pointed out that histogram equalization may cause over enhancement [7]. In such situations, histogram specification will be a good candidate for substitute methods because it accepts any shape of histograms. In this paper, although we would like to use a Gaussian distribution constantly for the target histogram in histogram specification to simplify our explanation, other choices such as bimodal or multimodal distribution may work well.
For digital images that are composed of pixels with quantized values, it is usual that the number of pixels in an image is greater than the number of possible quantized pixel values. Therefore, there are a large number of pixels that share the same quantized pixel value in a digital image, which means that exact histogram specification for digital images is an ill-posed problem [8]. For the effective utilization of all possible quantized pixel values and the suppression of the occurrence of false contours, a number of research works have been conducted for more than four decades since Hall published his work on discrete distribution transformation [9]. In addition to these technical needs, it is our genuine curiosity to make and see exactly histogram-equalized or specified images because the solution of the ill-posed problem is not unique. A promising approach to exact histogram specification is to obtain a meaningful strict ordering of all pixels in an image. Once we obtain the ordering of pixels, we can assign a limited number of pixel values to the ordered pixels so that the resultant histogram coincides with a specified target histogram.
The simplest way of pixel ordering is random ordering [10]. However, the random ordering is not a meaningful ordering. Coltuc et al. [11] described the process of exact histogram specification when a strict ordering is given, and proposed an ordering method, which is called the local mean method (LM), based on a filter bank of multiple moving average filters. Wan and Shi [12] proposed a wavelet-based exact pixel ordering algorithm (WA) for both exact histogram specification and image enhancement, which takes into account not only local mean intensity values, but also local edge information. Nikolova and Steidl [13] proposed a fast ordering algorithm for exact histogram specification based on a variational approach, which is equivalent to an iterative nonlinear filtering, and demonstrated that their variational ordering method outperforms both LM [11] and WA [12]. They also applied their algorithm to hue and range preserving enhancement of color images [14], where the histogram of an intensity image is exactly equalized or specified, and then an affine transform is applied to each color in a hue and range preserving manner.
In this paper, we propose a local contrast-based pixel ordering method for exact histogram specification, which computes the local contrast of each pixel by using a Gaussian filter [15] without approximation, and compare the proposed method with the state-of-theart method by Nikolova and Steidl [13]. Furthermore, the proposed method for grayscale images is extended to that for color images in a consistent manner. Experimental results demonstrate that Nikolova and Steidl's method can produce false patterns in flat areas of images, which are alleviated by the proposed method. For color images, it is experimentally shown that the proposed method outputs exactly histogram-equalized or specified images with closer hue to original images than that obtained by a separable method that utilizes Nikolova and Steidl's method for exact histogram equalization and specification of each color channel. These results demonstrate the effectiveness of the proposed method compared with the state-of-the-art method.
The rest of this paper is organized as follows: Section 2 provides a brief review of some related works. Section 3 briefly summarizes conventional histogram equalization. Section 4 summarizes conventional histogram specification and Nikolova and Steidl's exact histogram specification [13]. Section 5 proposes a local contrast-based exact histogram specification method, which is extended to that for color images. Section 6 shows experimental results, where the proposed method is compared with the state-of-the-art methods by Nikolova and Steidl [13] and Ramos et al. [16]. Finally, Section 7 concludes this paper.

Related Work
In this section, we review recent related works on histogram equalization and specification. Recently, Pallavi et al. [17] reviewed state-of-the-art techniques for image enhancement, in which the techniques of histogram equalization are also included such as contrast limited adaptive histogram equalization (CLAHE) [18] and brightness preserving bi-histogram equalization (BBHE) [19]. Trongtirakul and Agaian [20] proposed a weighted histogram equalization using entropy of probability density function, which outperformed related methods including a low-rank regularized retinex model (LR3M) [21]. More recently, Zhang et al. [22] proposed an unsupervised low-light image enhancement method via a histogram equalization prior (HEP) based on the observation that the feature maps of histogram equalization enhanced image and the ground truth are similar.
Hussain et al. [23] proposed a locally transformed histogram-based technique for dark image enhancement, which does not get affected from the over-enhancement problem. Balado [24] examined an optimum exact histogram specification problem and the inverse problem, and presented their closed-form performance analyses. Ramos et al. [16] proposed two algorithms for histogram specification and quantile transformation of data without local information. As one of the state-of-the-art algorithms for histogram specification, we implemented the first algorithm of Ramos et al. in this paper as described below. However, their algorithm does not produce an output image having the histogram being exactly identical to the target one.
The above review of recent related works draws a conclusion that Nikolova and Steidl's algorithm [13] is still one of the state-of-the-art algorithms for exact histogram specification.

Histogram Equalization
Let F = [ f ij ] be a grayscale image for i = 1, 2, . . . , m and j = 1, 2, . . . , n, where f ij ∈ {0, 1, . . . , L} denotes the pixel value at the position (i, j) on the image plane with m rows and n columns. The maximum pixel value L is typically given by L = 2 8 − 1 = 255 in an 8-bit grayscale. Then, the histogram of the pixel values in F is given by h = [h 0 , h 1 , . . . , h L ] whose kth element is computed by for k = 0, 1, . . . , L, where δ k, f ij denotes the Kronecker delta function defined by δ k, f ij = 0 for k = f ij , and 1 for k = f ij . Note that the histogram h is not normalized, i.e., ∑ L k=0 h k = mn = 1. The corresponding cumulative histogram is given by H = [H 0 , H 1 , . . . , H L ] whose lth element is defined by for l = 0, 1, . . . , L, and recursively computed by H l = H l−1 + h l for l = 1, 2, . . . , L with H 0 = h 0 . Histogram equalization converts a pixel value f ij into where f min = min i,j { f ij }, f max = max i,j { f ij }, and the 'round' operator rounds a given argument toward the nearest integer, and φ HE h denotes the function of converting pixel values in histogram equalization. Figure 1 shows an example of histogram equalization, where an original image in Figure 1a is converted into the image in Figure 1d, and their histograms and cumulative ones are shown in Figure 1b,c and Figure 1e,f, respectively, where the horizontal axis of each graphic denotes the pixel value ranging from 0 to 255, and the vertical axis denotes the (cumulative) number of pixels for each (cumulative) histogram. As shown in Figure 1e, the resultant histogram is not equalized actually, but the cumulative one is linearized as shown in Figure 1f. That is, to be more precise, the conventional histogram equalization is a cumulative histogram linearization.

Histogram Specification
In this section, we summarize conventional histogram specification, which is not an exact method, and then briefly summarize the state-of-the-art method for exact histogram specification by Nikolova and Steidl [13]. After that, we propose a local contrast-based method for exact histogram specification, whose objects to be processed are extended from grayscale images to color images.

Conventional Histogram Specification
Letĥ = [ĥ 0 ,ĥ 1 , . . . ,ĥ L ] be a target histogram into which we want to convert the histogram h of the original image F, where we assume that ∑ L k=0ĥ k = ∑ L k=0 h k for convenience. Then, we can compute the cumulative histogramĤ = [Ĥ 0 ,Ĥ 1 , . . . ,Ĥ L ] fromĥ in the same manner as above. Theoretically, we can describe the relationship between the target histogramĥ and the original one h as where f HS ij denotes the output pixel value of histogram specification corresponding to an input pixel value f ij , and (φ HÊ h ) −1 denotes the inverse function of φ HÊ h which equalizes the target histogramĥ. However, in the digital condition, the functions φ HE h and φ HÊ h become staircase functions, which are not bijective. Therefore, they have no inverse functions.
A prescription for this situation may be the use of linear interpolation as follows: . . , L}, where we define that φ HÊ h (−1) = 0. Then, we can determine the output pixel value of histogram specification by corresponding to an input pixel value f ij . Figure 2 shows an example of histogram specification, where the resultant image in Figure 2a is obtained with a Gaussian target histogram with mean 127.5 and standard deviation 50 from the original image in Figure 1a. Figure 2b shows the obtained histogram, which is not Gaussian, but the cumulative one fits into the target cumulative histogram as shown in Figure 2c, where a blue solid line denotes the target one, and the yellow dashed line denotes the obtained one. That is, to be more precise, the conventional histogram specification is a cumulative histogram specification.

Nikolova and Steidl's Method
In this subsection, we briefly summarize Nikolova and Steidl's fast ordering algorithm for exact histogram specification [13]. They proposed a fixed point algorithm that attains the minimizer of fully smoothed l 1 -TV (total variation) functional as follows: where [25] of F where the superscript T denotes the matrix transpose [26], (Gu) κ denotes the κth component of Gu ∈ R r with r = 2mn − m − n, where R denotes the set of real numbers, and G is a forward difference operator given by where ⊗ denotes the Kronecker product, I n is the n × n identity matrix, D n denotes the forward difference matrix defined by and the function θ is defined by with a constant α > 0, which is a smooth approximation of l 1 norm, and has the derivative These functions are illustrated in Figure 3, where blue lines denote the functions for l 1 norm, and orange lines denote their smooth approximations.
. . , mn, which are combined into a vector equation as u = f − ξ(βG T θ (Gu)), where the functions θ and ξ are applied to each element of their arguments. This is a fixed point equation for u which gives rise to a fixed point algorithm for minimizing J(u, f ) as described in Algorithm 1.

Algorithm 1 Fully smoothed l 1 -TV minimization algorithm [13]
Require: a vectorized image f = [ f 1 , f 2 , . . . , f mn ] T , parameters β > 0 and α > 0 used in functions ξ and θ , the number of iterations T > 0 Ensure: the minimizer u (T) of a fully smoothed l 1 -TV functional J(u, f ) in (6) 1: Initialize the variable u as u (0) ← f 2: for t ← 1, 2, . . . , T do 3: In this algorithm, we set α = 0.05, β = 0.1 and T = 5 according to Nikolova and Steidl [13]. The obtained image from u (T) with the input image f in Figure 1a is shown in Figure 4a. This procedure for computing u (T) from f can be viewed as a nonlinear filter which reduces the quantization noise in f slightly. Therefore, the image in Figure 4a is very close to that in Figure 1a. The pixel values in u (T) can be ordered in a strict way with a high probability [8].

Assume that the pixel values in u
λ mn where λ ι for ι = 1, 2, . . . , mn denote the ordered indices, and the cumulative histogramĤ = [Ĥ 0 ,Ĥ 1 , . . . ,Ĥ L ] is computed from a given target histogramĥ. Then, ifĤ l−1 < ι ≤Ĥ l withĤ −1 = 0 for an index ι and a pixel value l, then the corresponding output pixel value is given by Figure 4b shows the output image with the same Gaussian target histogram as Figure 2, and Figure 4c shows the target and obtained histograms, where we can see that they have the same Gaussian shape, which demonstrates the exactness of Nikolova and Steidl's method.

Proposed Method
In this section, we propose a local contrast-based pixel ordering method for exact histogram specification. We first describe the method for grayscale images, and then it is extended to that for color images. We also give a complexity analysis of the proposed method. The main idea of our approach is that the problem of pixel ordering can be divided into 256 subproblems for an 8-bit grayscale image. We solve the downsized subproblems separately to have 256 groups of ordered pixels, and then concatenate them to have a sequence of all ordered pixels.

Local Contrast-Based Exact Histogram Specification for Grayscale Images
As demonstrated above, we can transform the histogram of a digital image into the specified one, when all pixels in the image are ordered in a strict and faithful way [13]. In most cases, pixels in a digital image take a limited number of discrete values, e.g., 256 values are available for 8-bit grayscale images, which is smaller than the number of pixels in the image, e.g., m × n = 256 × 256 = 65, 536 256 for the above image in Figure 1a. Therefore, many pixels share the same pixel value with each other in the image, and the ordering of the pixels with the same values becomes an ill-posed problem [8]. However, it is fortunate for us that we have a rough ordering of 256 groups of pixels in the image, which means that the entire pixel ordering problem can be divided into 256 subproblems that can be solved separately. In this subsection, we propose a method for ordering pixels in each group on the basis of the local contrast of each pixel, and then the separately ordered subgroups are finally merged into an entire pixel ordering.
The first step of the proposed method is Gaussian filtering without approximation, which means that the weighted mean of all pixels is outputted at each pixel, where the weights are given by a Gaussian function as follows: for i = 1, 2, . . . , m and j = 1, 2, . . . , n, where σ denotes the standard deviation positive constant.
be an m × n matrix which expresses the Gaussian-filtered image of F. Then, we can compute F GF by the following matrix operations. The numerator of (10) is computed by where the left matrix G L is given by where denotes the Hadamard product or elementwise product of matrices, and 1 m denotes the m-dimensional column vector of ones. Similarly, the right matrix G R in (11) is given by On the other hand, the denominator of (10) is computed by where E denotes the m × n matrix of ones. Then, we have the Gaussian-filtered image by F GF = N u D e , where denotes the Hadamard division or elementwise division of matrices. Next, we define the local contrast (LC) of each pixel by d ij := f ij − f GF ij which is added to the corresponding pixel as an attribute, and divide the mn pixels in F into (L + 1) groups as S k := {( f ij , d ij ) | f ij = k}, where we have that h k = |S k | for k = 0, 1, . . . , L, i.e., h k is the cardinality of S k . Then, we sort the elements of S k in the ascending order of d ij for k = 0, 1, . . . , L, and concatenate them to have a sequence. As a result, we expect to have the entire pixel ordering µ ι ∈ {1, 2, . . . , mn} for ι = 1, 2, . . . , mn, from which the corresponding pixel position (i ι , j ι ) is given as the integer part (j ι − 1) and remainder (i ι − 1) of the division (µ ι − 1)/m, i.e., (µ ι − 1) = (j ι − 1)m + (i ι − 1), and {d ij } is ordered as follows: For exact histogram specification, we would like to have a strict ordering as described in (15). That is, for different pixels with the same pixel value as For conventional Gaussian filters with finite kernel sizes such as 3 × 3, 5 × 5 and 7 × 7 pixels, the lower bounds of the probability of f GF However, some symmetric patterns may result in the case of To reduce the probability of such cases, we adopt the Gaussian filter with sufficiently large kernel size as described in (10).
For an input image F with a target histogramĥ and its cumulative version H = [Ĥ 0 ,Ĥ 1 , . . . ,Ĥ L ], the proposed method outputs a histogram-specified image where each pixel value is given by f LC iι,jι = l for the index µ ι satisfyingĤ l−1 < ι ≤Ĥ l withĤ −1 = 0 as well as Nikolova and Steidl's method summarized in Section 4.2. The proposed method is summarized in Algorithm 2.
. . , L 4: ι ← 1 5: for k ← 0, 1, . . . , L do 6: Sort the elements of S k in the ascending order of d ij as d i 1 ,j 1 for u ← 1, 2, . . . , |S k | do 8: µ ι ← m(j u − 1) + i u Search for l satisfyingĤ l−1 < ι ≤Ĥ l withĤ −1 = 0 15: f LC i ι ,j ι ← l, where (j ι − 1) and (i ι − 1) are the integer part and remainder of (µ ι − 1)/m 16: end for Additionally, we can make an image whose pixel values indicate ordered numbers where each pixel value is given by f I iι,jι = ι for ι = 1, 2, . . . , mn. For histogram equalization, the number of pixels to which the same pixel value is assigned is given bȳ h = mn/(L + 1) , where x denotes the floor function that gives the largest integer less than or equal to x [27]. That is, we make a target histogramĥ = [ĥ k ] asĥ k =h for k = 0, 1, . . . , L. However, the truncation by the floor function may cause the shortage of total amount in the histogram: ∑ L k=0ĥk < mn. To cover the shortage, if Q = mn − ∑ L k=0ĥk > 0, then we add 1 tô h 0 ,ĥ 1 , . . . ,ĥ Q−1 to have ∑ L k=0ĥk = mn. This modified histogram can be inputted to Algorithm 2 as a target histogram for exact histogram equalization. Figure 5 shows the results of the proposed method, where Figure 5a,b show the results of the histogram equalization and specification with σ = 50 for the original image in Figure 1a, respectively, and Figure 5c shows their histograms, where we can see that both exact histogram equalization and specification are achieved.

Extension to Color Images
In this subsection, we extend the above proposed method for grayscale images to color images. Let C = [c ij ] be an RGB color image, where c ij = [c ij1 , c ij2 , c ij3 ] denotes the RGB color vector of the (i, j)th pixel in C, and c ijk ∈ {0, 1, . . . , L} for k = 1, 2 and 3. Then, we make a histogram into which all RGB values in C are united as follows: For example, Figure 6 shows an RGB color image and its histogram defined in (16).
Letĥ C be a target histogram with the cumulative oneĤ C = [Ĥ C This procedure for color images is summarized in Algorithm 3. Figure 7 shows the results of the proposed method for color images, where Figure 7a,b show the results of the histogram equalization and specification with σ = 50 for the original color image in Figure 6a, respectively, and Figure 7c shows their histograms, where we can also see that both exact histogram equalization and specification are achieved.

Complexity Analysis
In this subsection, we analyze the complexity of the proposed method in Algorithm 2. Let N = mn be the input size. The computation of Gaussian filtering in line 1 is described as matrix calculations. The time complexity of the numerator (11) is estimated as mn 2 + m 2 n = mn(m + n) = N(m + n) ≈ N 3/2 , which is the same as that of the denominator (14). The Hadamard division requires N elementwise divisions. Then, the order of time complexity of Gaussian filtering is O(N 3/2 ). The next main procedure is sort in line 6, which is repeated L + 1 times. Assume that the cardinality |S l | is N L+1 on average. Then, the time complexity of the sorting procedure is estimated as (L + 1) × N L+1 log N L+1 = N[log N − log(L + 1)] from which we have O (N log N). The time complexity of the remaining procedures is bounded by O (N). Consequently, the time complexity of the proposed method is given by O (N 3/2 ).
The space complexities of images F, F GF , D, and F LC are O(N), and that of (L + 1) groups S k for k = 0, 1, . . . , L in line 3 is also O(N). Histogram arraysĥ andĤ require O(L + 1). Selecting the largest one, we can see that the space complexity of the proposed method is O(N).

Experimental Results
In this section, we show experimental results in comparison with Nikolova and Steidl's method [13]. First, we show the results with a synthetic input image shown in Figure 8a. Nikolova and Steidl's method summarized in Section 4.2 outputs the image in Figure 8b, where we can see that vertical stripes appear in the flat areas in the original image. As a result, Nikolova and Steidl's method achieves the equalized histogram shown in Figure 8c. Figure 8d shows the horizontal profile of the output image in Figure 8b, where the vertical and horizontal axes denote the pixel value and the column index j, respectively. Although the pixel values fluctuate with small steps, the pixel values on the left side are larger than that on the right side as well as the original image.  Next, we compare the ability in strict ordering by Nikolova and Steidl's and the proposed methods. Figure 10 compares the filter outputs used in the two methods with the input image in Figure 1a. Nikolova and Steidl's filter in Algorithm 1 outputs u (T) , the values in which are sorted in ascending order in Figure 10a where the vertical and horizontal axes denote the pixel value and index, respectively. The differences between the neighboring pixel values in Figure 10a are very small, and therefore, we take the logarithm of them as shown in Figure 10b. Similarly, we sort the values of the output of the Gaussian filter in (10) as shown in Figure 10c, and take the logarithm of the difference values in Figure 10c as shown in Figure 10d, where we can see that the minimum value is greater than that in Figure 10c, which is preferable for strict ordering because the larger the difference is, the clearer the ordering is. The minimum values of the differences without taking logarithm are 4.23 × 10 −12 and 1.66 × 10 −8 for Figure 10b,d, respectively. The proposed method has a parameter σ in (10), the effect of which is investigated in Figure 11a, where the vertical and horizontal axes denote the minimum difference in the sorted output values of the Gaussian filter and the parameter σ, respectively, for the input image in Figure 1a. In this figure, although we cannot find clear tendency of the minimum difference as a function of σ, the obtained difference values are larger than that of u (T) given by Algorithm 1 that exemplifies the insensitivity of the proposed method to σ. Figure 11b,c show the results of histogram equalization with σ = 1 and 70, respectively. These images demonstrate that the proposed method is insensitive to the values of parameter σ.
We investigate the limitation of the proposed method on the value of σ. Figure 12 shows the results of histogram equalization for the input image in Figure 8a, where Figure 12a,b correspond to σ = 10 8 and 5 × 10 8 , respectively. Figure 12a is similar to Figure 9a. However, in Figure 12b, we observe an artifact. As a result, we estimate that the limitation of the value of σ exists between 10 8 and 5 × 10 8 , and recommend using the value of σ smaller than 10 8 . Figure 13 shows twelve grayscale images in the Standard Image Data-BAse (SIDBA) [28], which are enhanced by Nikolova and Steidl's exact histogram equalization as shown in Figure 14, where the contrast of all images is enhanced well because the histogram of each image is exactly equalized.   Similarly, Figure 15 shows the results by the proposed method, which also equalizes the histogram of each image exactly. Therefore, the output images of both methods are similar to each other. Figure 16a shows the root mean squared error (RMSE) between Nikolova and Steidl's output image and the corresponding output image by the proposed method, where "Text" image shown in (b) get the largest RMSE value, but the difference is small as visualized in (c), where error-free pixels have the neutral gray value 127.
We also compare the results of exact histogram specification with the Gaussian target histogram, which are shown in Figures 17 and 18 for Nikolova and Steidl's and the proposed methods, respectively. In Figure 17, we observe that the obtained contrast is not so high as that of Figure 14 because the Gaussian target histogram suppresses the numbers of extremely dark and bright pixels, and assigns many pixel values to middle-range gray pixels. As a result, the dynamic range of gray pixels is boosted relatively. The output images in Figure 18 have a similar tendency to that in Figure 17 because the proposed method also obtains an exact Gaussian target histogram as well as Nikolova and Steidl's method.   Figure 19a shows the RMSE between the corresponding output images in Figures 17 and 18, where three images "Bui.", "Lig." and "Tex." have relatively large RMSE values. We can find the differences at the bottom of those images as shown in Figure 19b, where the bottom areas of "Bui.", "Lig." and "Tex." are arranged from top to bottom. The corresponding results by Nikolova and Steidl's and the proposed methods are shown in Figure 19c and Figure 19d, respectively, where we can see artifacts in (c) similar to that in Figure 8b. On the other hand, Figure 19d has no such artifacts.  Figure 20 shows twelve color images in the SIDBA image dataset [28], which are enhanced by a separable histogram equalization that applies Nikolova and Steidl's exact histogram equalization algorithm to each color channel separately as shown in Figure 21, where we can see that the contrast of every image is enhanced, but the hue has changed from the original images in Figure 20.   Figure 22 shows the results from the proposed histogram equalization method described in Section 5.2, where we can see that the hue of each image is closer to the original one than that in Figure 21.  Figure 23 shows the RMSE of hue of each output image from the corresponding original image, where the vertical and horizontal axes denote the RMSE value and the images, respectively, and the blue and orange bars denote the separable and proposed methods, from which we observe that the proposed method preserves the original hue better than the separable method. Figure 24 shows the results of a separable histogram specification that applies Nikolova and Steidl's exact histogram specification algorithm to each color channel separately. Similarly to the above results of histogram equalization in Figure 21, here we can see the hue change from the original images.
On the other hand, Figure 25 shows the results by the proposed histogram specification method, where we can also see that the hue of each image is closer to the original one than that in Figure 24. Figure 26 shows the RMSE of hue of each output image from the corresponding original image, from which we also observe that the proposed method preserves the original hue better than the separable method.   The above experimental results demonstrate the effectiveness of the proposed local contrast-based pixel ordering method for exact histogram equalization and specification compared with the state-of-the-art method by Nikolova and Steidl [13]. A drawback of Nikolova and Steidl's algorithm observed in Figure 8b, where a stripe pattern occurred on flat areas, may be a newfound phenomenon. The proposed method can avoid the occurrence of such false patterns. In the above experiments, although we used a Gaussian target histogram for histogram specification constantly, we can choose other candidates such as multimodal distributions without any problems.
For the implementation of Nikolova and Steidl's algorithm described in Algorithm 1, we need to use sparse matrices even if the input image is not so large, e.g., the original image in Figure 1a with 256 × 256 pixels because the difference operator G in (7) becomes a very large matrix with r × mn ≈ 8.56 × 10 9 elements. On the other hand, the proposed method can be implemented without sparse matrices.
Additionally, we implemented Ramos's recent algorithm for histogram specification [16]. Figure 27 shows the result of Ramos's histogram specification, where Figure 27a shows the output image computed from the input image in Figure 1a with the Gaussian target histogram with mean 127.5 and standard deviation 50. Figure 27b shows the obtained histogram of the output image in Figure 27a in a blue line with the target histogram in the orange line. This figure shows that Ramos's algorithm is not an exact histogram specification method, and the obtained histogram is similar to that of the conventional histogram specification shown in Figure 2b. On the other hand, the proposed method gives exactly the same histogram as a prescribed target histogram as shown in Figure 2c. In Ramos's formulation of l p norm minimization, we set p = ∞. For details, please refer to Ramos's original paper [16]. For the implementation of Ramos's algorithm, we need to use sparse matrices with about mn × L elements as well as Nikolova and Steidl's one.
A limitation of the proposed method is that users need to prepare a target histogram for histogram specification. However, it may be difficult to know the optimal shape of the histogram for a given image in advance. The estimation of suitable target histograms for given images will be a subject to be considered in the future.

Conclusions
In this paper, we proposed a method for ordering pixels in an image for exact histogram specification, and compared it with the state-of-the-art algorithm by Nikolova and Steidl [13]. We divided the problem of pixel ordering into small subproblems which can be solved memory-efficiently. This idea makes the problem of exact histogram specification more tractable than ever. The proposed method uses the Gaussian filter without approximation instead of the nonlinear filter used in Nikolova and Steidl's algorithm. We use the local contrast defined with the Gaussian-filtered image as a key to ordering pixels. As a result, it was experimentally demonstrated that we could avoid the occurrence of false patterns that were observed in the results obtained by Nikolova and Steidl's algorithm.
We also extended the proposed method for grayscale images to that for color images. Experimental results showed that the proposed method keeps the hue of the original images better than conventional separable methods of histogram equalization and specification for color images.
Our future work will include the development of a method for making a preferable target histogram for histogram specification depending on a given input image.