Concise Radiometric Calibration Using The Power of Ranking
Han Gong Graham D. Finlayson Maryam M. Darrodi
cs.CV


[
    Jorge Drumond Silva
    December 30, 2023
=======================




Compared with raw images, the more common JPEG images are less useful for machine vision algorithms and professional photographers because JPEG-sRGB does not preserve a linear relation between pixel values and the light measured from the scene. A camera is said to be radiometrically calibrated if there is a computational model which can predict how the raw linear sensor image is mapped to the corresponding rendered image (JPEGs) and vice versa. This paper begins with the observation that the rank order of pixel values are mostly preserved post colour correction. We show that this observation is the key to solving for the whole camera pipeline (colour correction, tone and gamut mapping). Our rank-based calibration method is simpler than the prior art and so is parametrised by fewer variables which, concomitantly, can be solved for using less calibration data. Another advantage is that we can derive the camera pipeline from a single pair of raw-JPEG images. Experiments demonstrate that our method delivers state-of-the-art results (especially for the most interesting case of JPEG to raw).





§ INTRODUCTION

Many computer vision algorithms (photometric stereo <cit.>, photometric invariants <cit.>, shadow removal <cit.>, and colour constancy <cit.>) assume that the captured RGBs in images are linearly related to the actual scene radiance. However, the imaging pipeline in a digital camera is necessarily non-linear in order to produce perceptually-pleasing photos rather than its physically-meaningful counterparts. In this paper, we present a new rank-based radiometric calibration method which solves for the bi-directional mappings between the camera's RAW responses and the rendered RGBs produced by digital camera.

There is prior art in this field which models the pipeline with a large number of parameters (up to several thousand <cit.>) which both means a large corpus of data is required to uncover the pipeline and that there is at least tacitly the premise that the underlying pipeline is quite complex. The key insight in our approach is that post-colour correction (a 3 × 3 matrix correction) the linear corrected raw RGBs are to the greatest extent in the same rank order as the final rendered RGBs. Building on this insight, we develop a simple rank-based radiometric calibration model that “solves for” the camera pipeline with many fewer parameters and concomitantly needs much less training data.

In Fig. <ref>, we illustrate a conventional image reproduction pipeline that holds for many cameras <cit.>. An exemplar raw image, Fig. <ref>a, is mapped by a 3 × 3 colour correction matrix to give the image shown in Fig. <ref>b. The colour correction matrix implements several processing steps (illumination correction <cit.>, display RGB mapping <cit.>, and colour preference adjustments <cit.>). It is well-known that a display device cannot display all captured image colours that some RGBs will fall outside the RGB cube after mapping (the pixels marked in purple in Fig. <ref>b). We therefore need gamut mapping, <cit.>, to bring the colours back inside the cube as shown in Fig. <ref>c. Finally, the gamut mapped image is tone mapped to arrive at the final rendered output <cit.> shown in Fig. <ref>d. Tone mapping accounts for the display non-linearity <cit.>, dynamic range compression and some aspects of preference <cit.>.


The colour processing pipeline – for cameras, in general, can be written as Eqn. <ref>.

    [       P = f(Γ (Mρ))         =  Γ(f(Mρ))         ≈    LUT(ρ);                (1a)                (1b)                (1c) ]

Here and throughout this paper, ρ denotes a camera RAW and P refers to its rendered RGB counterpart. Respectively, the 3 × 3 correction matrix, gamut mapping and tone mapping are denoted by the matrix M and the functions Γ() and f(). The function f() can implement a single or three per-channel tone curves. Because gamut mapping only implements a small change in comparison with colour and tone mapping steps, the order of gamut mapping and tone mapping may be switched (Eqn. <ref>b & c), a property that we exploit in this paper. Equally, we can also roll all three processing steps into one and directly solve for a 3D LUT (Look-Up-Table) that maps RAW to rendered counterparts. This LUT function is denoted LUT() <cit.> in Eqn. <ref>c. Readers may refer to the top row of Fig. <ref> to link each mathematical function to our example processed image.

In radiometric calibration, given a set of ρ and P, one seeks to solve for the parametrised pipeline parts (M, Γ(), f() and LUT()). A disadvantage of the current best performing methods is that a great deal of data may be required to fit their assumed models. In Eqns. <ref>a and <ref>b, the gamut mapping step could be modelled by 1000s of Radial Basis functions <cit.> and in Eqn. <ref>c, the deployed LUT could also have several thousand control points. 


Our proposed method begins with the simple observation <cit.> that, assuming the gamut mapping step makes small changes to image colours, we expect the rank ordering of the rendered P to be the same as ρ multiplied by the correction matrix M (because a tone curves are always monotonically increasing). Suppose that two rendered (JPEG) responses – in the 1^st red colour channel – are denoted P^a_1 and P_1^b; and that P^a_1>P^b_1. The rank order of two corresponding raw red channel measurements post colour correction is written as M_1ρ^a>M_1ρ^b (where M_1 denotes the first row of M and ρ^a and ρ^b are a pair of raw RGBs). Rewriting: this implies that M_1(ρ^a-ρ^b)>0 which mathematically defines a half-space constraint. If we visualise the row vector M_1 as a point in 3-space then this inequality – which we call a ranking constraint – forces the point to be located in one half of 3-space but not the other. Because we have multiple pixels, each pair of pixels (2 raw and 2 JPEG RGBs) generates a half space constraint and intersecting all these constraints delimits the region in which M_1 must lie. Our experiments demonstrates that a small numbers of patches suffices to estimate M accurately. Once we have M we then find the best rank preserving tone curves f(). At this stage, only using M and f() we have a tolerable approximation of the pipeline. Indeed, we argue that our construction of M and f() also incorporates, to a first order, gamut mapping. Now we adopt (Eqn <ref>c) and find a 125-parameter LUT to “mop up” any remaining errors due to gamut mapping (higher order terms). 





§ RELATED WORK


 


Using the pipeline form Eqn <ref>b, Chakrabarti  <cit.> first solve for M and f() (in a least-squares sense) in iteration and then solve directly for Γ(). In their approach, f() is constrained to be a 7^th order monotonic polynomial. They model Γ() with the radial basis function (RBF) method of <cit.> where several thousands of RBFs are potentially used. A restriction of the above calibration is presented in <cit.> where  the gamut mapping Γ() is ignored. This less general model works tolerably well on many real pairs of raw and rendered images and this is a point we will return to later in this paper.
In either version (<cit.> or <cit.>), the coupled nature of the minimization means that a global minimum is not guaranteed to be found. Thus, a random start search is incorporated – multiple minimisations are carried out – in order to find their best parameters. 
Kim <cit.> solve for the pipeline in the form of Eqn. <ref>a and makes additional assumptions to decouple the optimization. They assume that images of the same scene are captured with respect to two or more exposures and their Γ() is a multi-thousand set of RBFs. Regarding solving for f(), Debevec  <cit.> showed how relating corresponding pixels under known exposure differences suffices to solve for f() (assuming there is no gamut mapping step). Importantly, in <cit.>, it was argued that for the set of desaturated pixels (RAWs far from the RGB cube boundary) the gamut mapping step has little or no effect and can be ignored.  Relative to this assumption, f() can be solved using the Debevec method. Given f() then the colour correction matrix M can be found (again using desaturated pixels). Though, for typical capture conditions, e.g. for most mobile phones,  multiple exposures are not available and so the requirement that multiple exposures are needed is a weakness in this method. Finally, in <cit.> a gamut mapping RBF network is “trained”.
Of course, if a large number of radial basis functions are used to model gamut mapping (as proposed in <cit.> or <cit.>) then solving for Γ() requires a large corpus of data. Further the application of gamut mapping is expensive and its parametrisation is large.

In <cit.> it was shown that is possible to ignore the underlying structure of the colour processing pipeline and directly solve for the best 3D surjective function – implemented as a LUT that maps the RAWs to rendered RGBs (Eqn. <ref>c).
Finally, in <cit.>, a method is presented for solving for f() by examining the edge distribution in an image. This method has the advantage that the method works for a single image (no need for multiple exposures) but the disadvantage that the method is sensitive to processing steps such as image sharpening which is used extensively in mobile phone image processing. 

 









§ THE RANK-BASED METHOD



As the reader shall see, to make the rank-based method work we need to assume that the gamut mapping step Γ() only makes small adjustments to colour. In fact our assumption is more nuanced. We assume that – to a first order – gamut mapping can mostly be implemented as an affine transform and that this affine transform can be folded into the colour correction  matrix M and the  monotonically increasing tone mapping functions f(). 



 §.§ Gamut Mapping as an Affine Transform


In Eqn. <ref>b, gamut mapping is applied when, after colour correction, colours are mapped outside the colour cube and become non-displayable. 

Let us use a Taylor expansion to model Γ() around a point a inside the gamut:

    Γ(Mρ)≈Γ(a)+J(a)(ρ-a)

where J is the 3 × 3 Jacobian (matrix of derivatives of Γ). Not only does Eqn. <ref> show that, to a first approximation, gamut mapping is an affine transform it is also one of the gamut mapping algorithms proposed in <cit.>. In particular, <cit.> solves, with good results, for the best affine transform that maps image colours inside the gamut and which are, simultaneously, close to the non-gamut mapped colours:

    min_T,oΣ_i ||TMρ_i+o-Mρ||^2  s.t. 0≤ TMρ_i+o≤1

In Eqn. <ref>, T and o are respectively a 3 × 3 matrix and 3 × 1 offset vector defining the affine gamut mapping algorithm. The 3-vectors of 0s and 1s are denoted 0 and 1. Eqn. <ref> is solved directly by Quadratic Programming <cit.>. The gamut mapping shown in Fig. <ref>c is the result of solving Eqn. <ref>.

Here, we make two important remarks about affine gamut mapping: 1) Gamut mapping and colour correction combined can be represented by the single affine transform: 3 × 3 matrix TM and offset o; 2) It follows that the rank-based method presented in the next section will actually solve for TM. The offset term can be incorporated directly in f() (since an offset does not change ranking).







 §.§ Rank-based estimation for colour correction

Let us denote the k^th row of M as M_k, let us assume that given two colour corrected RAWs, M_kρ^a and M_kρ^b that the rank order is the same as for the corresponding rendered RGBs:

    P^a_k>P^b_k ⇒ M_kρ^a>M_kρ^b ⇒  M_k(ρ^a-ρ^b)>0

Defining the difference vector d^j=ρ^a-ρ^b:

    M_kd^j>0

where it is understood the superscript ^j denotes the difference vector from the j^th of n2 pairs of n image pixel values. Suppose that we have a vector M_k where Eqn. <ref> holds, then the inequality cannot be true for -M_k. That is Eqn. <ref> defines a half plane constraint <cit.>. The vector d^j is perpendicular to the half-plane: any M_k less than 90 degrees to d^j is a possible solution. Given multiple difference vectors then we have multiple half-plane constraints which taken together delimit a region in 3-space where M_k must lie. Denoting the half-plane as H(d^j), M_k must satisfy:

    M_k∈⋂_j  H(d^j)



Let us visualise the computation of M_k using ranking. Without loss of generality let us assume that M_k,3=1. We rewrite Eqn. <ref> as

    M_k,1d_1^j+M_k,2d_2^j+d_3^j>0

If [a b c] is a solution to Eqn. <ref>, then [a/c b/c c/c] for Eqn. <ref> is also true since M_k,1=a/c and M_k,2=b/c. Solutions for [M_k,1,M_k,2] lie on one side of the line, the 3D half-space constraints maps directly to a 2D half-plane constraint. Or, if we consider the whole set of collations, the cone in 3D, defined by Eqn. <ref>, maps to a 2D convex region <cit.>. Denoting half-planes as P(d^j) we, equivalently, solve for

    [M_k,1,M_k,2]∈⋂_j  P(d^j)

The intersection problem of Eqn. <ref> is easily visualised. In Fig. <ref>a we show the intersection of 4 half plane constraints and indicate the solution set where M_k must lie. 


We solve for M one sensor channel at a time. Empirically, we have to be careful not to generate too many half planes. In our experiment, we generate half-planes from all pairs of up to 50 randomly selected unique RAWs, generating 2450 half-planes. Due to noise or small deviations in real camera data, it is likely that no common intersection can be found that satisfy every half-planes constraint. To solve this problem, we generate 100,000 unit length vectors that are uniformly distributed on the surface of the unit sphere <cit.>, which is visualised in Fig. <ref>b. With respect to this sampling, the furthest distance between any point and its nearest neighbour is less than 1.15 degrees. So, the orientation of the rows of M are found to this accuracy. For each point on the sphere (a possible row of M_k), we count how many half-space constraints are satisfied. The point on the unit sphere that satisfies most half-plane constraints – or the median of multiple points if there is a tie – defines M_k.
To make our approach robust, we find randomly select 50 colours 25 times and for each trial find the best M_k. Overall, we find the M that places all the corresponding raw and rendered image RGBs in the most similar rank order. That is, if we plot the mapped raw red responses, for example, against the rendered red JPEG corresponding values then the graph should be a monotonically increasing function. How well a monotonically increasing function fits our data can be used to judge the efficacy of each M.

Ranking can only estimate M up to an unknown scaling of its rows. Suppose for a rendered achromatic RGB P^A=[0.5 0.5 0.5]^⊺ and the corresponding raw ρ^A=[a b c]^⊺, we apply: 0.5diag(Mρ^A)^-1Mρ^A=[0.5 0.5 0.5]^⊺ where diag() places a vector along the diagonal of a diagonal matrix. After this step, M← 0.5diag(Mρ^A)^-1M maps achromatic colours correctly.

Because M← DM (where in the example, D=0.5diag(Mρ^A)^-1) we might also solve for D in a least-squares sense by including all colours indexed ^i close to the achromatic axis: min_D ∑_i
||DMρ^i-P^i|| (our experiment does not include this additional step).




 §.§ Rank-preserving optimization of tone curves

We now solve for the optimal per-channel tone curves which map colour corrected RAWs to corresponding rendered RGBs. Let us denote the i^th colour corrected RAW and rendered RGB pixel pairs for the k^th channel as (ρ_k,i,P_k,i). Then, the k^th-channel rank-preserving tone curve f_k() is optimised as a 7^th order monotonic and smooth polynomial function as follows:

    min_f_k()Σ_i ||f_k (ρ_k,i)-P_k,i||^2 + λ∫_t||f_k”(t)||^2 dt   s.t.  f_k'()≥0
    .

where the first term is for data fitness, the second term is for curve smoothness and λ is a small weight (10^-5). This polynomial fitting is solved by Quadratic Programming <cit.>. 
In this paper, we further denote the combination of all 3-channel mappings f_1-3() as f().



 §.§ Gamut correction step

As argued previously, we propose that f(Mρ) has the expressive power to implement colour correction, tone correction and gamut mapping (to the first order in a Taylor expansion). However, we wish to add a further gamut mapping step for the higher order terms. But, since our hypothesis is that much of the gamut mapping will have been accounted for we are going to adopt a simple small parameter solution. Further, this additional correction is going to be carried out at the end of the process, we adopt Eqn. <ref>b.
Specifically, we find a 5×5×5 LUT by using lattice regression <cit.> that minimises min_LUT()Σ_i ||LUT(f(Mρ_i))-P_i||^2.




 §.§ Rank-based recovery of raw

Suppose we wish to map rendered RGBs to RAWs. Using the method presented in Section 3.2, M has already been solved in the RAW-to-JPEG forward estimation phrase. Now, in a least-squares optimal way, we use the same polynomial fitting method (Eqn. <ref>) to find f^-1 by optimising min_f^-1()Σ_i ||f^-1(P_i)-Mρ_i||. Finally, we solve for the backward LUT() by optimising min_LUT()Σ_i ||LUT(M^-1f^-1(P_i))-ρ_i || where the LUT is fitted by a 5×5×5 lattice regression <cit.>.



 §.§ Parameter counting


Assuming we solve for 3 independent tone curves then our method requires 9 (for M) + 8 × 3 (for f()) + 125× 3 (for Γ()) = 408 parameters which is significantly less (even an order of magnitude less) than <cit.>.



§ EVALUATION




Our evaluation is based on the most challenging dataset that we have encountered:  <cit.> which contains the RAW/JPEG intensity pairs of 140 colour checker patches viewed under multiple viewing conditions. Specifically, the colour chart is captured by 8 cameras (3 of which are JPEG-only) and under 16 illuminants across many different exposures. 

Below, we carried out the same experiment described in <cit.>. We are interested in validating whether our method, with much reduced number of parameters can produce, similar or even better results compared with <cit.>. We evaluate both RAW-to-JPEG and JPEG-to-RAW.
The dataset <cit.> captures a sort of “worst-case” viewing conditios. Normally, when we capture a picture there is a single prevailing illuminant colour. In the dataset of Chakrabarti , all camera processing parameters are turned off and then the same reflectances are viewed under multiple coloured lights. As Forsyth observed <cit.>, the reddest red camera response cannot be observed under a blue light. And, then he exploited this observation to solve for the colour of the light. Yet, in this dataset the reddest red, the greenest green and the bluest blue can all appear simultaneously in the same image. Practically, we believe the need to divine a pipeline for the all lights all surfaces case means the prior art pipelines are probably more complex than they need to be. 

As described in <cit.>, for each camera, we estimate the parameters of a calibration model using different subsets of available RAW-JPEG pairs. For each subset and a selected camera, the root mean-squared error (RMSE)  between the prediction and ground truth is validated by using all available RAW-JPEG pairs. Table <ref>a shows the RAW-to-JPEG mapping error (where pixel intensities are coded as integers in the interval [0,255]. In the table, Prob denotes the Chakrabarti method (with several thousands parameters) and RB the rank-based method with 408 parameters.
We found that our forward errors are close to the results of <cit.>, especially for the condition of less than 3 illuminants which are more likely to occur in the real world. Evidently, for the many illuminant case the prior art has a small advantage. Remembering that JPEGs are coded as integers in [0,255] the RMSE is typically 1 or less (for RB compared to Prob). Practically, when the “fits” are viewed visually (by looking at images) it has hard to see the difference between the two methods.




For computer vision, we are more interested in the performance of JPEG-to-RAW mapping which is shown in Table <ref>b and Table <ref>c. In <cit.>, a probabilistic framework for mapping rendered RGB to raw was presented. Here we take their mean estimates as the most likely RAW predictions. We found that our method generally reduces the errors of <cit.> by ∼ 34%. Our supplementary material also includes the additional experiment results compared with “<cit.> + our LUT” for interested readers.

The reader might be interested why our simple method seems to work so well going from rendered to raw (better than <cit.>) but not quite as well as the prior art in the forward direction (albeit visually almost indistinguishable). Our hypothesis here is that the LUT in the forward direction is applied post the tone curve. This curve (at least for dark values) has a very high slope and, consequently, the coarsely quantised 5×5×5 LUT cannot capture gamut mapping well. Yet, in the reverse direction (JPEG to RAW) the LUT is applied in linear raw where a course uniform quantisation is more justified.











§ CALIBRATION WITH SMALL NUMBERS OF PARAMETERS






We wished to visually validate our claim that we can calibrate with few parameters. We took 4 RAW+JPEG pairs (for different cameras) from <cit.>. We then uniformly selected 140 corresponding pixels from the RAW and JPEG. We solved for all parameters in our rank-based method. We then applied our model to the rest of the image. The result of this experiment for 4 images (JPEG-to-RAW) is shown in Fig. <ref>.






§ CONCLUSION


In this paper we have shown how the rank order of image responses is a powerful tool for solving for the individual steps in a camera processing pipeline (colour correction, gamut and tone mapping). A simple ranking argument, relating colour corrected RAWs to corresponding rendered RGBs suffices to solve for the colour correction matrix. Then, the rank-preserving tone map is found and, finally, a simple gamut correction step is derived. Compared with the prior art, our rank-based method requires the fewest assumptions and delivers state-of-the-art radiometric calibration results. 


