See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/224149953

Algorithms for color look-up-table (LUT) design via joint optimization of node
locations and output values

Conference Paper  in  Acoustics, Speech, and Signal Processing, 1988. ICASSP-88., 1988 International Conference on · April 2010
DOI: 10.1109/ICASSP.2010.5495312 · Source: IEEE Xplore

CITATIONS READS

6 1,307

2 authors, including:

Raja Bala
Samsung
101 PUBLICATIONS   618 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Computer Vision and Imaging in Intelligent Transportation Systems, Pub. Wiley-IEEE Press View project

All content following this page was uploaded by Raja Bala on 18 August 2014.

The user has requested enhancement of the downloaded file.

ALGORITHMS FOR COLOR LOOK-UP-TABLE (LUT) DESIGN VIA JOINT OPTIMIZATION
OF NODE LOCATIONS AND OUTPUT VALUES

Vishal Monga Raja Bala

Department of Electrical Engineering, Xerox Research Center Webster,

Pennsylvania State University, 800 Phillips Rd., Webster, NY,
University Park, PA 16802, USA 14580, USA

ABSTRACT storage. For most applications, processor RAM and cache memory
Real-time processing constraints entail that non-linear color trans- constraints render storage of such full-resolution LUTs impractical.
forms be implemented using multi-dimensional look-up-tables (LUT). Further, the dimensionality of color devices is likely to increase in
Further, relatively sparse LUTs (with efficient interpolation) are em- the future with devices employing more than 3 or 4 colors, e.g. 6
ployed in practice because of storage and memory constraints. Much and higher colorant printing [2].
research has been devoted towards optimizing “nodes” (or equiva- Hence real-world color look-up-tables are built with a sparser
lently partitioning the input color space) of this color LUT based on sampling of the input color space where the sampled values are re-
the curvature of the color transform to be processed through the LUT. ferred to as the “nodes”. Efficient interpolation models [3] are used
Likewise, for a given LUT structure, the optimization of transform to determine output values corresponding to input color variables
output values has been suggested so as to minimize interpolation er- that do not coincide with the nodes of the color LUT. The problem of
ror in an expected sense even if the values stored in the LUT do “node selection” for color LUTs has been extensively researched. To
not agree with true transform output values. This paper presents a some extent, the problem may be likened to quantization [4] where
principled algorithmic approach to combine the merits of these two once again an optimal partitioning of the input color space is desired
complementary techniques. The error (cost) function does not ex- to approximate a range of colors via a reduced palette. The crucial
hibit joint convexity over the multidimensional variable sets of node difference is in the optimization cost function. For color LUT de-
locations and corresponding output values which makes this opti- sign, it is desired to optimize the partitioning or equivalently node
mization particularly challenging. The paper makes two significant selection so as to effectively approximate a multi-dimensional color
contributions: 1.) for the case of simplex interpolation, a cost func- transform processed through the LUT. The curvature of the trans-
tion is formulated that exhibits separable convexity in its arguments form, hence plays a dominant role. As prior work has shown, truly
and enables an efficient alternating convex optimization algorithm, optimal node selection in multiple-dimensions reduces to problems
and 2.) in the aforementioned framework, for fixed node outputs, the in the N P -class and is computationally infeasible. Much attention
optimization of node locations is split into a primary and an auxil- has hence been devoted to sub-optimal and practical algorithms [5],
iary optimization, which greatly improves the quality of the solution [6], [7] that trade-off complexity vs. quality of approximation.
over traditional alternatives where node locations are directly opti- A slightly less explicitly investigated aspect of color LUT design
mized. Preliminary experiments show remarkable improvements in is the problem of output value optimization [8]. That is, while storing
color transform accuracy over what is obtained by individually opti- the exact value of the multi-dimensional color transform will lead to
mizing just the node locations or output values. zero error at the nodes of the LUT, in an expected sense the error
may still be a lot higher as contributed by input color values whose
1. INTRODUCTION output is determined via an interpolation. It is therefore desirable to
optimize the output transform values at these node locations so as to
Device image processing pipelines, such as those in a digital cam- minimize error over a statistically representative training set of input
era or printer often entail the application of color transformations to colors under a given model of color interpolation.
images and objects. These transformations are complex non-linear While the two problems of node location vs. node value opti-
functions making it impractical to process large images in real-time. mization have been investigated individually, we argue and subse-
It is therefore common to implement such mathematical transforms quently demonstrate in this work that significant benefits are to be
via multidimensional look-up-tables (LUTs). made by their joint exploitation. We establish first the fundamental
The most general definition of a color LUT is a sparse data difficulty of this task by formulating the cost function which is to
structure that maps a multi-dimensional input color variable to an- be minimized as a function of node locations as well as output val-
other multi-dimensional output color variable by storing a cleverly ues. To make the problem tractable we particularly analyze the case
sampled subset of all possible color inputs (and corresponding out- of simplex interpolation in arbitrary dimensions. This allows us to
puts) from a multi-dimensional space. The qualifier “sparse” is of map the set of multi-dimensional node locations to a (constrained)
paramount significance here. Well known device color spaces such weight matrix which manifests itself in an attractive analytical form
as RGB, CM Y K as well as perceptually significant device inde- in the cost function. It is then observed that the cost function of in-
pendent color spaces such as CIELAB [1] are of ≥ 3 dimensions terest is not jointly convex over the weight matrix induced by the
with each dimension requiring a bit-depth of 8−bits or higher for node locations and the multi-dimensional variable of node values.
adequate representation. A 16−bit per separation, full-resolution This inhibits the direct application of standard numerical techniques
CM Y K → CIELAB LUT will therefore need 216×4 = 16 GB of in convex optimization [9]. The cost function however exhibits sep-

978-1-4244-4296-6/10/$25.00 ©2010 IEEE 998 ICASSP 2010

arable convexity as a function of the vector of the node output val-
ues for a fixed weight matrix (node locations), and vice-versa. The
paper then proposes an iterative algorithm called RESOLVE (RE-
peated Sequential Optimization of Node Location and ValuE) which
involves repeated solving of constrained and unconstrained convex
problems in an alternating manner. Experimental results for the
problem of designing a CM Y K → CIELAB LUT for a color
xerographic printer show dramatic improvements in approximation (a) Uniform node locations and (b) Node locations optimized
outputs coinciding with transform
accuracy by the proposed algorithmic framework over individual op- outputs
timization of node locations or output values.

2. ILLUSTRATION OF THE UNDERLYING CONCEPT

The benefit of jointly optimizing node locations and node values can
be understood by examining a 1-D example. With attention to Fig.
1 (a), the blue curve showns a 1-D function to be approximated. It
is sampled to form nodes at equal intervals (yellow circles), and ap- (c) Node values optimized (d) Node locations and values op-
proximated by a piecewise linearly interpolated function (shown in timized
red). The shaded areas provide a sense for the approximation error
introduced by the interpolation. For a fixed number of nodes, one Fig. 1. Illustration of the benefits of joint optimization of node loca-
approach to reduce approximation error is to optimize the node lo- tions and output values for a 1 − D function.
cations. Since the interpolation is piecewise linear, intuition tells us
that nodes should be placed more densely in regions where the func-
and carefully chosen to adequately represent transform characteris-
tion exhibits high curvature, and sparsely where the function exhibits
tics. Also, note that the above notation implicitly assumes that we
low curvature (or strong linearity). Such an approach would result
are considering color transforms that are Rn → R maps - this is
in an approximation as shown in Fig. 1 (b). Qualitatively, we readily
not a conceptual restriction but rather merely for ease of exposition
see that the approximation error has been reduced. Mathematically
and notational convenience as will be apparent in the following de-
the formulation and solution of a global optimization is non-trivial,
scription. An Rn → Rm map can be seen as a concatenation of m
especially in higher dimensions. However approximate greedy tech-
Rn → R mappings.
niques based on local curvature, or more generally an ”importance
In the ensuing discussion, a multi-dimensional input color that
function”: can be devised [5].
forms a node of the look-up-table will be denoted as xj nd , j =
Another consideration is how the node values are derived. Thus
1, 2, .., M . Let X nd denote the set of the M input nodes, and the
far, the node values are obtained by evaluating the original function
corresponding output node values as yj nd which may be collected
at the node locations. While this produces zero approximation error
(in an ordered fashion) into a vector of output node-values ynd .
at the nodes, it is not necessarily the optimal strategy for arbitrary in-
As illustrated in the previous section, both the node locations (in
puts [6]. This is illustrated in Fig. 1 (c) where node values are chosen
the input color-space) as well as the corresponding output node val-
to minimize an overall error, rather than to precisely match the func-
ues affect the ability of the LUT to approximate the true transform.
tion at the node locations. If the error criterion is mean squared error,
In particular, the output node values do not have to be identical to
then it possible to obtain an optimal set of node values by solving a
the true transform value but may be optimized so as to yield lower
standard least-squares problem, as explained in Section 4.
error overall for a given interpolation model.
While prior techniques have focused on solving for either node Analytically, node locations as well as node values are desired
location or node value, we propose an approach that combines their to be optimized so as to minimize the following cost function:
complementary merits and attempts to optimize both to minimize a
chosen error criterion. This is shown qualitatively in Fig. 1 (d).
T
2
C(X nd , ynd ) = (yi − fˆ(xi , X nd , y nd )) (2)
i=1
3. NOTATION AND PROBLEM SET-UP
where fˆ(xi , X nd , ynd ) = ˆ
xj nd ∈N (xi ) wj yj . That is, f () is
nd

In this discussion, we will strive to optimize a LUT that employs an estimate of the true transform f () for the training point xi and
simplex interpolation [10] in n−dimensions. For generality, we will is formed by the weighted combinations (an interpolation) of node-
not enforce the restriction that the LUT nodes lie on a regular lat- output values corresponding to the input node locations that form the
tice. The input to our algorithm includes a training set S of multi- vertices of the simplex N (xi ) that encloses xi . A 2-D example is
dimensional inputs and outputs obtained from the “true” color trans- shown in Fig. 2.
form f (), where S is given by: The above formulation suggests that it is possible to write the
functional approximation as a linear operator on the vector of node
output values ynd . In particular, the optimization problem of interest
S = {xi , yi }, i = 1, 2, .., T, xi ∈ Rn , yi ∈ R (1) may now be written as:

Or succinctly, S = {X tr , ytr } where X tr is the set of all train- (X nd opt , ynd opt ) = arg min ytr − WX nd .ynd 2 (3)
(X nd ,ynd )
ing inputs and ytr is the (suitably ordered) T × 1 vector of all cor-
responding training outputs. The key requirement of the training set That is, the goal is to solve for a set of node locations X nd and
is that it should be a large enough sampling of the input color space corresponding node values ynd such that the interpolation error over

999
Algorithm 1 SOLVE
1: Using the training set S, estimate the curvature of the transform
and employ a known sub-optimal technique [7], [6], [5], that
allocates more nodes in regions of high transform curvature. Let
∗
X nd denote the resulting node locations.
∗
2: For fixed node locations X nd (a corresponding fixed WX nd ∗
is induced), and for T > M , (3) reduces to a well known un-
constrained convex problem, the solution to which is given by:

∗
Fig. 2. Illustration of simplex interpolation in 2-D carried out over a ynd = (WX nd ∗ T WX nd ∗ )−1 WX nd ∗ T ytr (4)
set of nodes. The training point is denoted by “x”.
Note that the entries of the weight matrix WX nd ∗ are cor-
responding to the training inputs xi . These are determined
the statistically representative training set S = {X tr , ytr } may be uniquely as the barycentric coordinates [10] w.r.t the enclos-
minimized. In (3), ytr is the T × 1 vector of training output values, ing simplex given a simplex topology associated with the set of
∗
WX nd is a T × M matrix of interpolation weights (determined as node locations X nd .
nd ∗ nd ∗
a function of X tr and X nd ) that selectively operates on the entries 3: X and y are determined as the node locations and output
of the vector of node output values ynd . In particular, the i−th row values respectively.
of WX nd , i = 1, 2, .., T contains at most n + 1 entries that are
non-zero (for n-dimensional input nodes) and they are corresponding
to the output node values in ynd of the enclosing simplex, which Algorithm 2 RESOLVE
contribute to the interpolation. 1: Initialization - same as Step 1 of the SOLVE algorithm to yield
∗
node locations X nd .
∗
2: As in Step 2 of SOLVE, optimize node ouput values to get ynd .
4. ALGORITHMS FOR LOOK-UP-TABLE (LUT) DESIGN nd ∗
3: Based on the optimized vector of node output values y in
Step 2, optimize the weight matrix as:
4.1. Cost Function Analysis and Algorithm Development
Examining the cost function in (2), it is not difficult to see that ∗
C(X nd , ynd ) is not jointly convex as a function of X nd and ynd . WX nd opt = arg min ytr − WX nd .ynd 2 (5)
(WX nd )
In fact, for fixed ynd , the function is not even guaranteed to be con-
vex in X nd and as much past research has shown [5], [7] clever subject to:
sub-optimal algorithms are needed to yield good solutions that seek
to minimize interpolation error.
However, for the case of simplex (and more generally any linear) (i) Membership Constraint : wiT ei = 0. i = 1, 2, ..T. (6)
interpolation, the revised formulation of (3) is much more appealing. (ii) Non-negativity : wi,j ≥ 0, i = 1, 2, ..T ; j = 1, 2, ..M. (7)
It is of interest to observe that the cost function now exhibits separa-
ble convexity in WX nd and ynd . That is, for a fixed weight matrix (iii) Interpolation Constraint : WX nd .1 = 1. (8)
WX nd (induced by node locations X nd ) the function is convex in In (6), wi ∈ RM represents the i − th row of WX nd and ei is
ynd and vice-versa. Note that, we are still not claiming convexity a membership vector comprising only of zeros and ones. And in
of the cost function in X nd but only in the corresponding induced (8), 1 ∈ RM is a vector with all its entries equal to 1.
weight matrix WX nd . This separation is crucial to the formulation 4: Map the optimized WX nd opt to a corresponding set of node
as well as the numerical nature of the algorithms presented next. ∗
Two variants are proposed: 1.) a one shot sequential optimiza- locations X nd . In particular, solve:
tion of node locations followed by output value optimization; 2.) an
iterative algorithm that solves constrained and unconstrained convex ∗
X nd = arg min cij (WX nd (i, j)−WX nd opt (i, j))2 (9)
problems in an alternating manner. The first variant is actually con- X nd
ij
ceptually rather simple but sets the tone for the key contribution of
this paper, namely the second variant. where cij = K.eij + . Here eij is the matrix of ones and
zeros formed by collecting the membership vectors in a row-
4.2. Sequential Optimization of node Location and ValuE (SOLVE) wise fashion. Both K, > 0 where in practice, K is chosen
large while is small.
The stepwise description of SOLVE is given in Algorithm 1. 5: Repeat steps 2-4 till no further minimization of the cost function
The algorithm simply uses a state of the art technique for deter- is possible.
mining node locations and identifies that the node output values may
be uniquely optimized by solving a convex least squares problem.
Albeit simplistic, the algorithm may be used when constraints on
LUT structures are severe, (e.g. regular lattices as required by ICC The stepwise description of RESOLVE is given in Algorithm 2. In
profiles [11]) and do not lend themselves to further optimization. step 3 of Algorithm 2, (5) represents a convex optimization problem
with three convex constraints. The first constraint guarantees that the
4.3. REpeated - Sequential Optimization of node Location and training point of interest for which the optimization of weights (and
ValuE (RESOLVE) consequently node locations) is carried out, maintains membership
in the same simplex. ei is a vector of ones and zeros with ones cor-

1000
responding to node locations (simplex vertices) that do not enclose Node Location Node value # of Nodes Avg. ΔE 95% ΔE
the training point. The second and third constraints account for the Uniform Hi-res LUT 4096 6.2 8.2
fact that the entries in each row of WX nd are interpolation weights SSD as in [5] Hi-res LUT 4096 4.1 7.4
- hence they must be non-negative and sum to one. Uniform Least-squares 4096 5.6 7.89
Step 4 is based on the (fairly non-trivial) observation that the SOLVE-L SOLVE-L 4096 3.11 6.1
mapping from node locations1 X nd to the corresponding weight ma- SOLVE-G SOLVE-G 1800 2.95 6.1
trix WX nd is potentially a many-one mapping with a non-unique RESOLVE RESOLVE 1800 1.89 4.83
inverse. In some cases, it may even be impossible to determine a
set of node locations that exactly yield the optimized WX nd opt . We
choose to therefore conduct a search on node locations X nd via (9) Table 1. Benefits of optimizing both node location and value in
which minimizes a weighted Frobenius norm between the desired approximating a forward printer color transform.
optimal weight matrix obtained via (5) and the one induced by node
locations. The weighting is cleverly contrived so the solution favors
induced weight matrices that are close in structure and content to the leads to the lowest error demonstrating the strength of the alternate
optimized one. convex optimization framework.
It is useful to re-emphasize here that the formulation of the in-
terpolation weight matrix induced by node locations, and the sub- 6. CONCLUSION
sequent spilt into a primary optimization as in (5) and an auxiliary
optimization in (9) to determine node locations is a very significant A principled mathematical approach has been presented for optimiz-
idea in making the joint optimization framework tractable. As is ing the accuracy of a multidimensional color LUT transform. Key
well-known [7], a direct optimization of multi-dimensional node lo- contributions are: 1.) the observation that joint optimization of LUT
cations to minimize (2) is painfully slow and yields rather poor re- nodes and values is superior to separate optimization of either quan-
sults in terms of the quality of minima obtained. tity alone; 2.) the formulation of a cost function that exhibits sepa-
rable convexity with respect to the weighting matrix WX nd , deter-
mined only by node locations X nd , and the vector ynd determined
5. EXPERIMENTAL RESULTS only by node values, thus enabling an alternating convex optimiza-
tion solution; 3) for fixed ynd , the derivation of WX nd via a con-
In order to demonstrate the virtues of optimizing both node locations strained convex optimization in conjunction with an auxiliary La-
and values, we carried out the following experiment. Because the grangian optimization to ensure a WX nd that is induced by actual
availability of the true transform as a compact analytical closed-form node locations. Results demonstrate markedly improved accuracy.
is impractical for problems like this, a very high resolution (334 ) While the application focuses on color transforms, this approach can
CM Y K → CIELAB LUT was first designed by printing and be applied for any LUT-based functional approximation problem.
measuring about 1500 CM Y K patches using state of the art tech-
niques in building printer-models [1]. This extremely hi-resolution 7. REFERENCES
LUT was treated as the true transform.
An independent test set of 216 CM Y K vectors spanning the [1] R. Bala, “Device characterization,” Digital Color Imaging Handbook:
Chapter 5, pp. 687–725, 2003.
printer gamut were processed through the high-resolution LUT to
obtain “ground truth” CIELAB color vectors. The same CM Y K [2] D. Pankow, Tempting the palette: a survey of color printing processes,
color vectors were then processed through a given low-resolution RIT Cary Graphic Arts Press, 2005.
LUT approximation built using a variety of different techniques (in [3] R. Bala and V. Klassen, “Efficient color transformation implementa-
Table 1), and the output was compared against the ground-truth in tion,” Digital Color Imaging Handbook: Chapter 11, pp. 687–725,
terms of the CIE76 ΔE metric [1]. 2003.
Results of the aforementioned experiment are presented in Ta- [4] A. Gersho and R. Gray, Vector Quantization and Signal Compression,
ble 1. The first four rows of Table 1 show approximations via color Springer-Verlag, 1990.
LUTs that are restricted to be regular lattices of size 84 = 4096 [5] V. Monga and R. Bala, “Sort-select-damp: An efficient strategy for
nodes. Uniform lattice with no optimization, individual optimiza- color look-up-table lattice design,” SID/IS&T Color Imaging Confer-
tion of node locations and node values, and the SOLVE algorithm ence, Nov 2008.
instantiated on the lattice (denoted by SOLVE-L) respectively, com- [6] M. J. Baines, “Algorithms for optimal discontinuous piecewise linear
prise the first four rows of Table 1. The SOLVE-G and RESOLVE and constant l2 fits to continuous functions,” Mathematics of Compu-
tation, pp. 645–669, Apr 1994.
entries in rows 5 and 6 do not have the lattice restriction.
The benefits of optimizing both node location and node value vs. [7] J. Z. Chang, J. P. Allebach, and C. A. Bouman, “Sequential linear
interpolation of multidimensional functions,” IEEE Trans. on Image
standard uniform LUT design and individual optimization of either
Processing, pp. 1231–1245, Sep 1997.
LUT node locations or output values are readily apparent. Note also
that if the constraint of a regular lattice is removed, the number of [8] T. Balasubramanian, “Iterative technique for refining color correction
lookup tables,” US Patent 5649072, 1995.
LUT nodes M requires for a given level of accuracy can be signifi-
cantly reduced. Finally, the iterative structure of RESOLVE clearly [9] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge Uni-
versity Press, 2004.
1 In general, conducting efficient search over node locations to obtain de- [10] P. Hemmingway, “n-simplex interpolation,” HP Labs Technical Report,
sirable weight matrices is an open problem. The solution presented here is Nov 2002.
reasonably easily implemented but not necessarily the fastest or numerically [11] INTERNATIONAL COLOR CONSORTIUM, “What is an icc pro-
most stable. We are currently investigating this and will report a detailed file?,” http://www.color.org/iccprofile.xalter.
analysis in future work.

1001

View publication stats