Method for matching color in images

Abstract

A method and apparatus for transforming a first color distribution based on a second color distribution. The method includes the steps of determining, for each of the first and second color distributions, a one-dimensional histogram along a direction in a color space; matching the one-dimensional histogram determined for the first color distribution and the one-dimensional histogram determined for the second color distribution so as to generate a transform mapping; transforming the first color distribution based on the generated transform mapping; and repeating the determining, matching, and transforming steps for other directions in the color space until the generated transform mapping converges.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention is described in conjunction with the appended figures, wherein:

FIG. 1 is an example (right image) of color transfer using Reinhard's Color Transfer [15] illustrating how the transfer of a source image (center image) fails to re-synthesize the color scheme of the target image (left image), and that human interaction is required to make the process successful;

FIG. 2 is an example of a one-dimensional (1D) pdf transfer on grayscale pictures with FIG. 2A showing a source and a target image, and FIG. 2B showing the source and target histograms and a matching function that can transform pixels from the source to match the histogram of the target;

FIGS. 3A-3D illustrate a outline flow diagram (FIG. 3A) and various stages of a more detailed flow chart (FIGS. 3B-3D) with the first part of the process being optional, but capable of handling very effectively linear transformations, whilst the second step of the process can handle any non-linear transformation and at the end of the iterations, when no more improvement is visible, the two set of samples share the exact same statistics; while using this process for matching colors in images, a grain artifact can appear and the noise smoothing module is used to address this problem;

FIG. 4 is an example of the new color transfer process as applied to a target image (a) and a source image (b) picture using three iterations (c-e), with the final result (f) having the same color distribution as the target image (a);

FIG. 5 illustrates the result of noise smoothing on grayscale images where two consecutive archive frames (a) and (b) suffer from extreme brightness variation and the corresponding mapping transformation (e) is overstretched resulting in an increased level of noise on the mapped original frame (c) whilst the proposed grain artifact reducer is able to reproduce the noise level of the original picture though leaving the top of the original picture saturated and not retrieved but succeeding in preserving the soft gradient in image (d);

FIG. 6 illustrates the result of noise smoothing for a color picture with the details of the picture being preserved, while the spurious graininess in sky has been washed out; and

FIG. 7 is an example of employing the novel color grading process for matching lighting conditions by using the color properties of the sunset to synthesize the ‘evening’ scene depicted at sunset in accordance with an embodiment of the present invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

Specific details are given in the following description to provide a thorough understanding of the embodiments. However, it will be understood by one of ordinary skill in the art that the embodiments may be practiced without these specific details. For example, circuits may be shown in block diagrams in order not to obscure the embodiments in unnecessary detail. In other instances, well-known circuits, processes, algorithms, structures, and techniques may be shown without unnecessary detail in order to avoid obscuring the embodiments.

Also, it is noted that the embodiments may be described as a process which is depicted as a flowchart, a flow diagram, a data flow diagram, a structure diagram, or a block diagram. Although a flowchart may describe the operations as a sequential process, many of the operations can be performed in parallel or concurrently. In addition, the order of the operations may be re-arranged. A process is terminated when its operations are completed, but could have additional steps not included in the figure. A process may correspond to a method, a function, a procedure, a subroutine, a subprogram, etc. When a process corresponds to a function, its termination corresponds to a return of the function to the calling function or the main function.

Furthermore, embodiments may be implemented by hardware, software, firmware, middleware, microcode, hardware description languages, or any combination thereof. When implemented in software, firmware, middleware or microcode, the program code or code segments to perform the necessary tasks may be stored in a machine readable medium such as storage medium. A processor(s) may perform the necessary tasks.

An embodiment of the present invention is a process for automating the color grading or color balancing task even when the lighting conditions have dramatically changed, as shown in FIG. 1. It resolves many of the problems with previous work and is computationally extremely efficient. Its implementation based on the below description in dedicated or purposely designed circuitry is possible, but an implementation exploiting the Graphics Subsystem of general purpose processors is described as a more economic manner of implementing the present invention.

Considering now the task of transferring the color from target image (left) to source image (center) in FIG. 1. One simple approach is to separately match the red, green and blue channels using the standard 1D histogram matching technique. After one iteration of this matching process, there is no further change in the mapping since the 1D histograms no longer change. Typically this methods yields a poor result, as shown in image on the right of FIG. 1.

One aspect of this invention is to notice that if a colorspace transformation (a rotation of the colorspace) is applied to the color data between each iteration of 1D matching, this causes the entire color surface distribution of the center image of FIG. 1 to eventually warp into the shape of the color distribution of the target image on the left of FIG. 1. The image on the right of FIG. 1 shows the final color transfer appearance in which source image has been colored with the colors in the target image.

It is an important observation for the purpose of present invention to note that even if a random sequence of rotations is used, the process converges eventually to a valid picture with some successful color transfer properties.

In addition, it was found that the transformation need not be a rotation, although other transformations may be less efficient.

The invention proposes in this example that if 1D slices of the multidimensional color space are matched independently, eventually the entire multidimensional surface will be mapped effectively. To understand this aspect of the invention, it may be beneficial to discuss some background on image histograms.

An image histogram can be understood as the frequency distribution of color in the image. The histogram of a gray scale image is 1D because pixels assume only scalar values of gray—typically from 0 to 255 for instance. FIG. 2B shows this. The histogram of a color image is 3D because each pixel has associated with it 3 values for red, green and blue (in RGB space), or luminance, and two color components (in YUV space).

Matching 1D histograms is a well understood process and used in histogram equalization for instance. For example, in case of the two gray scale images shown in FIG. 2A it may be required to transfer the brightness distribution of the target image on the right to the source image one on the left. To do so requires mapping each of the 256 gray scale values in the image on the right to some other grayscale value such that the resulting brightness of source image is the same as target image. A technique for creating this mapping is well known [4]. A process may be represented as the following steps.

1. Set the histogram of the target image to be stored in some memory array H₁ such that the n^th element of that array is H₁[n] for n=0:255 elements. Similarly set the histogram of source image 2 to be H₂[n]. Both distributions are shown in FIG. 2B.

2. Calculate the cumulative density functions C₁, C₂ by summing the elements of each memory array. Hence calculate

$C_1\left[n\right]=\sum _{m=1}^nH_1\left[m\right]$

and

$C_2\left[n\right]=\sum _{m=1}^nH_2\left[m\right]$

3. Create the required grayscale mapping by scanning the elements of C₁ and mapping the grayscale of each entry in C₁ to the grayscale value that indexes the value having the least difference in C₂

For each m=0:255

Find that value of n which makes |C₁[m]-C₂[n]| smallest. Define this value as T[m].

The mapping for m is then T[m]

The result is then an array T[m] which contains, for each element m=0:255, a corresponding number that is the grayscale to be assigned. The mapping is shown in FIG. 2B. To apply the transformation to source image then, each pixel in the source image is visited and its grayscale, I say, is used to look up in T[m] to yield the new grayscale. Hence for a pixel of grayscale I the new grayscale to be assigned is T[I].

There are variants on this process of course, but the above serves as a mechanism for achieving 1-D histogram matching. The resulting image should then have an intensity distributed similarly to the target image.

Unfortunately, the above method fails for color images because there are more dimensions involved. A naive approach based on the grayscale process above would be to separate each color component (red, green and blue) and match the histograms of each color plane individually. An example of doing this is shown in FIG. 1. The match is clearly poor.

However, it is an aspect of the current invention that the above approach can be successful by transforming (step 301) the color space before performing the histogram matching (steps 302 and 303), and then repeating this process iteratively 305, as shown in the diagram of FIG. 3A. To return to the original colorspace, the rotation is then undone in step 304 or reversed after each mapping. In other words by applying the matching operation for multiple directions, i.e., more that just the three conventional axes, in the color space until the changes to the pixels are no longer significant, even color images can be successfully matched.

The central part of the process outlined above is illustrated in greater detail in FIG. 3C. There it is assumed that, after some pre-processing to be described below when referring to FIG. 3B, the process has moved to an arbitrary iteration k. The image pixels involved are the current image I^(k) (which has been transformed up to iteration k) and the image data from G, the target image. G is the image containing the color distribution that is required. The first step 321 of the iteration is to change the coordinate system by rotating both the color of the samples of I^(k) and the samples of G.

Defining this rotation operation as R:

$\begin{array}{cc}\begin{array}{c}R=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(\alpha \right)& \sin \left(\alpha \right)\\ 0& -\sin \left(\alpha \right)& \cos \left(\alpha \right)\end{array}\right)\left(\begin{array}{ccc}\cos \left(\beta \right)& 0& -\sin \left(\beta \right)\\ 0& 1& 0\\ \sin \left(\beta \right)& 0& \cos \left(\beta \right)\end{array}\right)\\ \left(\begin{array}{ccc}\cos \left(\gamma \right)& \sin \left(\gamma \right)& 0\\ -\sin \left(\gamma \right)& \cos \left(\gamma \right)& 0\\ 0& 0& 1\end{array}\right)\end{array}& \left(1\right)\end{array}$

for some angles α,β,γ, then each pixel color vector I^(k)=(r, g, b) is transformed in I^(k)=(r′, g′, b′) as follows

$\begin{array}{cc}\left(\begin{array}{c}{r}^{\prime }\\ g^{\prime }\\ b^{\prime }\end{array}\right)=R\left(\begin{array}{c}r\\ g\\ b\end{array}\right)& \left(2\right)\end{array}$

This operation is can be also referred to as a color space conversion rotating one set of axes in the color space into different set of axes in color space.

In a second step 322, the samples of both distributions are projected onto the new axis by summing in the relevant direction. Hence summing along the transformed red and green color axes yields a projection along the transformed blue axis. This yields three histograms for the three color planes in each sequence. These histograms for the new axes are defined as H_I¹,H_I²,H_I³ and H_G¹,H_G²,H_G³ for the current and target image, respectively.

Then it is possible using the ID process above to find, in the next step 323, for each axis 1,2,3, the mappings t₁, t₂, t₃ that transform the histograms H_I¹,H_I²,H_I³ into H_G¹,H_G²,H_G³. The resulting transformation t maps a pixel of color I₁, I₂,I₃ onto t(I₁,I₂,I₃)=(t₁(I₁),t₂(I₂),t₃(I₃)).

The iteration is completed with a rotation in step 324 of the samples by R⁻¹ to return in the original coordinate system.

When no new iterations can modify the iterated samples of I^(k), the process is considered to have converged in step 325 to the final solution and the color statistics of I^(k) and G are identical.

The process described above is not restricted to 3 dimensions. It can be applied to vector spaces of any dimensionality. The basic steps described above remain the same, the only difference is that the rotation matrices fit the relevant dimension, i.e., for a vector space of 4 dimensions, the rotation matrix becomes a 4×4 matrix etc.

A full process according to an example of the invention includes an initialization step 311 and optional steps 312, 313 to accelerate the process as shown in FIG. 3B. This full process includes the following steps.

INTIALIZATION: Retrieve the pixels of the original image I and the target image G with G containing the desired color distribution. For example, the image pixel I=(I₁,I₂,I₃) encapsulates the red, green and blue components. Set k←0, I⁽⁰⁾←I.

REPEAT:

• • a. Take one rotation matrix R;
  • b. Rotate the pixels samples: I_r←RI^(k) and G_r←RG;
  • c. Project the samples on all 3 axis to get the histograms H_I¹,H_I²,H_I³ and H_G¹,H_G²,H_G³;
  • d. For each axis i, find the 1D transformation t_i that matches the histograms H_I¹,H_I²,H_I³ into H_G¹,H_G²,H_G³;
  • e. Remap the pixels colors I_r according to the 1D transformations. For example, a sample (I₁,I₂,I₃) is remapped into (t₁(I₁),t₂(I₂),t₃(I₃));
  • f. Rotate back the pixel colors to the original color space: I^(k+1)←R⁻¹I_r
  • g. k←k+1

2. UNTIL convergence on all histograms for every possible rotation

Although the above method can be applied to any kind of mapping, it might be desirable as a first step 312, 313 to use a different technique for registering global transformations between the two distributions. In particular, the novel method is not naturally designed to find rotation between the original and target dataset. This is not to be confused with the rotations used in the PDF transfer process of FIG. 3C. The rotation considered here is the rotation between the two distributions, whereas the rotations mentioned in the PDF transfer process above are applied to both distributions and only describe a change of coordinate system. The linear mapping that we are looking for in these optional steps of FIG. 3B is of the form G=AI+b.

A simple solution to this problem is to use the Principal Component Analysis 312 of each of the images.

(I−Ī)^t(I−Ī)=U_I^tD_IU_I   (3)

and

(G− G)^t(G− G)=U_G^tD_GU_G   (4)

where U_I=[u_I¹,u_I²,u_I³] and U_I=[u_G¹,u_G²,u_G³] are the 3×3 orthogonal matrices containing the eigenvectors of the covariance matrices of the pixel color vectors of I and G. The diagonal matrices D_I and D_G contain the eigenvalues corresponding to the eigenvectors in U_I and U_G.

Registering the pixel colors vectors of I⁽⁰⁾ to the target distribution G is then possible by employing the following linear transformation process 313:

I⁽⁰⁾←(U_GD_G^1/2)(U_ID_I^1/2)⁻¹(I−Ī)+ G  (5)

The color distribution of the image I⁽⁰⁾ shares the average value, variance and global color directions with the target image G. This operation cannot handle non-linear transformations but offers an accelerated process if it is used as a first step, i.e., before the actual pdf transfer process.

To obtain a correct transformation, the eigenvectors in U_I and U_G have to correspond. This is achieved by ordering them with respect to the magnitude of the corresponding eigenvalues and making sure that they do not point in opposite directions. i.e.

∀ i≦3, ^tu_Iⁱ·u_Gⁱ>0   (6)

The color mapping to the original picture transfers correctly the target color palette to the original picture but it might also produce some grain artifacts, as shown in FIG. 5. When the content differs or the dynamic range of both pictures is too different, the resulting mapping function can be stretched on some parts (see particularly FIG. 5), and thus enhances the noise level.

This effect can be understood by taking the simple example of a linear transformation t of the original picture X: t(X)→aX+b. The overall variance of the resulting picture is changed to var(t(X))=a var(X). Hence a greater stretching (a>1) produces greater noise.

A solution to this noise problem in accordance with an variant of the present invention is to reduce the grain by executing a post-processing process, as shown in FIG. 3D, that forces the level of noise to remain the same or approximately the same. In its most basic form, this variant includes the step 333 of adjusting the gradient field of the picture result so that it matches the original picture. If the gradient of both pictures are similar, the level of noise will be the same or at least approximately the same. Matching the gradient of a picture has been addressed in different applications such as image stitching [12] or high dynamic range compression [3], and it can be efficiently solved using a variational approach to be described in the following.

Denoting as I(x,y) the original picture and omitting the coordinates in the expressions to simplify the following discussion, i.e. I, J, ψ, φ, etc. actually refer to I(x,y), J(x,y), ψ(x,y) and φ(x,y), the color transformation is t:I→t(I). The problem posed by the above solution is to find a modified image J of the mapped picture t(I) that minimises on the whole picture range Ω the expression

$\begin{array}{cc}\underset{j}{\min }\int {\int }_{\Omega }\phi ·{\nabla J-\nabla I}^2+\psi ·{J-t\left(I\right)}^2& \left(7\right)\end{array}$

F=φ·∥∇J−∇I∥ ² +ψ·∥J−t(I)∥²

from which the following can be derived:

φ·J−div(ψ·∇J)=φ·t(I)−div(ψ·∇I)   (12)

This is an elliptic partial differential equation. The expression div(Ψ·∇I) can be approximated using standard finite differences as in [17] by: (13)

$\mathrm{div}\left(\psi ·\nabla I\right)\approx \sum _{l\in N_i}\frac{{\psi }_I+{\psi }_i}{2}\left(I_l-I_i\right)$

where N_i corresponds to the four neighboring pixels of i.

This yields the following linear system of step 332:

a(x,y)J(x,y−1)+b(x,y)J(x,y+1)+c(x,y)J(x−1,y)+d(x,y)J(x+1,y)+e(x,y)J(x,y)=f(x,y)   (14)

with with Neumann boundary condition ∇J|_∂Ω=∇I|_∂Ω to match the the gradient of J with the gradient of I at the picture border. The term ∥CJ−CI∥² forces the image gradient to be preserved. The term ∥J−t(I)∥² ensures that the colors remain close to the target picture. Without ∥J−t(I)∥², a solution of expression (7) will be actually the original picture I.

The weight fields φ and ψ affect the importance of both terms of expression (7). Numerous expressions for ψ and φ are possible. For example, the weight field φ can emphasize that only flat areas have to remain flat but that gradient can change at object border, in which case:

$\begin{array}{cc}\phi =\frac{{\phi }_0}{1+\nabla I}& \left(8\right)\end{array}$

(with Φ₀ constant). A possible choice 331 for the weight field ψ is to account for the possible stretching of the transformation t. Where ∇t is big, the grain becomes more visible:

$\begin{array}{cc}\psi =\frac{{\psi }_0}{1+\left(\nabla t\right)\left(I\right)}& \left(9\right)\end{array}$

where Ψ₀ is constant and (∇t)(I) refers to the gradient of t for the color I.

The minimization problem of expression (7) can be solved using the variational principle which states that the integral must satisfy the Euler-Lagrange equation:

$\begin{array}{cc}\frac{\partial F}{\partial J}-\frac{}{x}\frac{\partial F}{\partial J_x}-\frac{}{y}\frac{\partial F}{\partial J_y}=0\mathrm{where}& \left(10\right)\\ a\left(x,y\right)=-\frac{\psi \left(x,y-1\right)+\psi \left(x,y\right)}{2}b\left(x,y\right)=-\frac{\psi \left(x,y+1\right)+\psi \left(x,y\right)}{2}c\left(x,y\right)=-\frac{\psi \left(x+1,y\right)+\psi \left(x,y\right)}{2}d\left(x,y\right)=-\frac{\psi \left(x+1,y\right)+\psi \left(x,y\right)}{2}e\left(x,y\right)=\frac{1}{2}\left(4\psi \left(x,y\right)+\psi \left(x,y-1\right)+\psi \left(x,y+1\right)+\psi \left(x-1,y\right)+\psi \left(x+1,y\right)\right)+\phi \left(x,y\right)f\left(x,y\right)=\frac{1}{2}(\psi \left(x,y\right)+\psi \left(x,y-1\right)\left(I\left(x,y-1\right)-I\left(x,y\right)\right)+(\psi \left(x,y\right)+\psi \left(x,y+1\right)\left(I\left(x,y+1\right)-I\left(x,y\right)\right)+(\psi \left(x,y\right)+\psi \left(x-1,y\right)\left(I\left(x-1,y\right)-I\left(x,y\right)\right)+\left(\psi \left(x,y\right)+\psi \left(x+1,y\right)\left(I\left(x+1,y\right)-I\left(x,y\right)\right)\right)+\phi \left(x,y\right)I\left(x,y\right)& \left(11\right)\end{array}$

The final process is identical to a two-dimensional infinite impulse response (2D IIR) filter in which the coefficients a,b,c,d,e and f only dependent on the original picture I and the mapped picture t(I) and which thus can be computed beforehand. The unknown result picture J can be solved by standard iterative methods 333 like SOR, Gauss-Seidel with multigrid approach. Implementations of these numerical solvers are widely available and one can refer for instance to the respective numerical recipes [14].

FIG. 5 illustrates the result of noise smoothing on grayscale images where two consecutive archive frames (a) and (b) suffer from extreme brightness variation and the corresponding mapping transformation (e) is overstretched resulting in an increased level of noise on the mapped original frame (c) whilst the proposed grain artifact reducer is able to reproduce the noise level of the original picture, though leaving the top of the original picture saturated and not retrieved but succeeding in preserving the soft gradient in image (d).

FIG. 6 illustrates the result of noise smoothing for a color picture with the details of the picture being preserved, while the spurious graininess in sky has been washed out.

FIG. 7 is an example of employing the novel color grading process for matching lighting conditions by using the color properties of the sunset to synthesize the ‘evening’ scene depicted at sunset in accordance with an embodiment of the present invention.

The main step of these methods is to solve iteratively for J(x,y) as follows.

(15)

$J^{\left(k+1\right)}\left(x,y\right)=\frac{1}{e\left(x,y\right)}\left(f\left(x,y\right)-a\left(x,y\right)J^{\left(k\right)}\left(x,y-1\right)-b\left(x,y\right)J^{\left(k\right)}\left(x,y+1\right)-c\left(x,y\right)J^{\left(k\right)}\left(x-1,y\right)-d\left(x,y\right)J^{\left(k\right)}\left(x+1,y\right)\right)$

The main computational demands in the first step of the full process are the colorspace conversion and the look up table for mapping pixels. In the second noise smoothing step, the solution of the system requires a substantial amount of spatial processing. It is possible to propose the use of Graphics Hardware from companies such as NVIDIA and ATI to provide a hardware embodiment.

An efficient system can be created by using a general purpose CPU for control and data management while employing the Graphics subsystem for data processing. The main subsystems are proposed as below.

It is well known that Graphics Processing Units (GPU's) can perform multiply/add operations on images [6]. Colorspace conversion is in particular built in as a function on chips from NVIDIA and ATI. The calculation of RI in equation (2) can be performed on such chips directly.

The creation of the 1D mapping look up table can be performed using LUT operations on the GPU. The entire image can be mapped in a single operation using the GPU functionality.

It is well known that all the fundamental operations required for Gauss-Seidel Multigrid optimization can be performed on the GPU [1, 9]. A process description has been previously presented in [5]. That module can be used with minimal change for this application.

It is not required to present all the details of the GPU programming needed to achieve the above, as such programming is considered to be well within the scope of an ordinarily trained person in the relevant field when being taught the steps as described above. This programming implementation will change with generations of GPUs and generations of tool development kits. However, currently, the GPU represents a readily available cheap computational resource that can be taken advantage of for the purpose of implementation the present invention. With the processes as outlined above, frame rates at near real time are achieved with current hardware.

All embodiments of the present invention conveniently may be implemented using a conventional general purpose computer or micro-processor programmed according to the teachings of the present invention, as will be apparent to those skilled in the computer art. Appropriate software may readily be prepared by programmers of ordinary skill based on the teachings of the present disclosure, as will be apparent to those skilled in the software art.

A computer may implement the methods of the present invention, wherein the computer housing houses a motherboard which contains a CPU, memory (e.g., DRAM, ROM, EPROM, EEPROM, SRAM, SDRAM, and Flash RAM), and other optional special purpose logic devices (e.g., ASICS) or configurable logic devices (e.g., GAL and reprogrammable FPGA). The computer also includes plural input devices, (e.g., keyboard and mouse), and a display card for controlling a monitor. Additionally, the computer may include a floppy disk drive; other removable media devices (e.g. compact disc, tape, and removable magneto-optical media); and a hard disk or other fixed high density media drives, connected using an appropriate device bus (e.g., a SCSI bus, an Enhanced IDE bus, or an Ultra DMA bus). The computer may also include a compact disc reader, a compact disc reader/writer unit, or a compact disc jukebox, which may be connected to the same device bus or to another device bus. Examples of computer readable media associated with the present invention include compact discs, hard disks, floppy disks, tape, magneto-optical disks, PROMs (e.g., EPROM, EEPROM, Flash EPROM), DRAM, SRAM, SDRAM, etc. Stored on any one or on a combination of these computer readable media, the present invention includes software for controlling both the hardware of the computer and for enabling the computer to interact with a human user. Such software may include, but is not limited to, device drivers, operating systems and user applications, such as development tools. Computer program products of the present invention include any computer readable medium which stores computer program instructions (e.g., computer code devices) which when executed by a computer causes the computer to perform the method of the present invention. The computer code devices of the present invention may be any interpretable or executable code mechanism, including but not limited to, scripts, interpreters, dynamic link libraries, Java classes, and complete executable programs. Moreover, parts of the processing of the present invention may be distributed (e.g., between (1) multiple CPUs or (2) at least one CPU and at least one configurable logic device) for better performance, reliability, and/or cost.

While detailed descriptions of one or more embodiments of the invention have been given above, various alternatives, modifications, and equivalents will be apparent to those skilled in the art without varying from the spirit of the invention. Moreover, except where clearly inappropriate or otherwise expressly noted, it should be assumed that the features, devices and/or components of different embodiments can be substituted and/or combined. Thus, the above description should not be taken as limiting the scope of the invention, which is defined by the appended claims.

Claims

1. A method for transforming a first color distribution based on a second color distribution, comprising: determining, for each of the first and second color distributions, a one-dimensional histogram along a direction in a color space;matching the one-dimensional histogram determined for the first color distribution and the one-dimensional histogram determined for the second color distribution so as to generate a transform mapping;transforming the first color distribution based on the generated transform mapping; andrepeating the determining, matching, and transforming steps for other directions in the color space until the generated transform mapping converges.

2. The method of claim 1, wherein the repeating step is performed at least three times such that the matching step is performed for at least four directions in the color space.

3. The method of claim 1, wherein the repeating step comprises determining whether the generated transform mapping converges by comparing the one-dimensional histogram determined for the first color distribution and the one-dimensional histogram determined for the second color distribution.

4. The method of claim 3, further comprising: applying a transformation to the first and second color distributions between each iteration.

5. The method of claim 4, wherein the transformation is a rotation.

6. The method of claim 1, wherein the matching step comprises matching the one-dimensional histograms using a cumulative histogram.

7. The method of claim 1, further comprising: obtaining the first color distribution from a source image; andobtaining the second color distribution from a target image.

8. The method of claim 7, further comprising: noise smoothing the transformed source image.

9. The method of claim 8, further comprising: constraining gradients in the source image to be close to gradients in the target image.

10. The method of claim 8, further comprising: minimizing a weighted sum of a squared difference of image gradients between a desired noise reduced image and the target image, and a weighted sum of squared differences between said desired noise reduced image and the transformed source image.

11. The method of claim 7, further comprising: rendering the transformed source image using a graphics processing unit.

12. The method of claim 11, wherein the graphics processing unit is programmed to perform color space transformation, the method further comprising: mapping input pixels into output pixels using a look-up-table embedded in said graphics processing unit.

13. The method of claim 1, further comprising: using principal component analysis based on the first and second color distributions.

14. An apparatus for transforming a first color distribution based on a second color distribution, comprising: means for determining, for each of the first and second color distributions, a one-dimensional histogram along a direction in a color space;means for matching the one-dimensional histogram determined for the first color distribution and the one-dimensional histogram determined for the second color distribution so as to generate a transform mapping;means for transforming the first color distribution based on the generated transform mapping; andmeans for causing the repeated execution of the means for determining, means for matching, and means for transforming for other directions in the color space until the generated transform mapping converges.

15. The apparatus of claim 14, further comprising a graphics processing unit.

16. A computer program product embedded on a computer-readable medium, the computer program product including instructions which, when executed by a computer, cause the computer to transform a first color distribution based on a second color distribution by performing the steps of: determining, for each of the first and second color distributions, a one-dimensional histogram along a direction in a color space;matching the one-dimensional histogram determined for the first color distribution and the one-dimensional histogram determined for the second color distribution so as to generate a transform mapping;transforming the first color distribution based on the generated transform mapping; andrepeating the determining, matching, and transforming steps for other directions in the color space until the generated transform mapping converges.