Geffray S., Klutchnikoff N., Vimond M. (2016), Illumination problems in digital images. A statistical point of view, Journal of Multivariate Analysis, 150: 191-213

Summary of the article: Interest is focused on a multi-dimensional signal \(R\) observed on a grid in the presence of both an illumination artifact and an additional additive bounded variance centered noise \(\varepsilon.\) The illumination artifact consists of colour or grey level intensity variations which are seen on the sampled image but which are not present in \(R.\) Such an assumption is classically modelled using a function \(L\) which acts multiplicatively on \(R.\) Our goal is to estimate \(R\) from observations of a random variable \(Y\) which obeys the regression model \(Y=RL+\varepsilon.\) We identify and propose a solution to an identifiability issue. We construct a consistent multi-step estimation procedure. We first make use of any consistent denoising method to estimate from the noisy data the gradient of the logarithm of the regression function. We assume that the artefact \(L\) consists of “smooth” variations. Then we project the denoised auxiliary estimate on a finite basis of “smooth” functions. We deduce the final estimator of \(R\) so that the proposed identifiability constraint is satisfied. We derive an upper bound for the sup-norm risk of our estimator. Applications to different images are presented.

Aim of this page: help you find your own way with our publicly available code!

Hardware requirement: to run our code on your data, you should be working with the latest free open source R software (www.r-project.org) or at least with R Version 3.3.1 installed on it. Moreover, your system should be run at least under windows seven or Ubuntu 10.4 or Mac OS X El Capitan 10.11.1.

If your are a beginner R user, then first go to our help page. You can skip it, if you already know how to use R.

Download our code: Download the compressed file in the prefered format and extract the archive into a new empty repertory (named as “GKV” throughout this tutorial):

• our R code in a .zip file
• our R code in a .tar file

Contents of the archive: Here is a schematic description of the contents of the file.

Our code is written to be used according to the following guidelines:

1. Set R working directory
2. Load additional packages
3. Load your data in R workspace. Images have to be in the .png format.
4. You can do some preliminary analysis (plot your data, plot graylevel histogram, check some pixel values, check data size)
5. Call the main function sample-analysis.R (and modify it if needed).
6. Plot the denoised image.
7. Plot the illumination estimation (actually you get a level plot of the illumination).
8. Plot the reflectance estimation (this is your processed image which is free of both noise and illumination).

Get R workspace ready

Set R working directory

Our code is written to work with “GKV” as R working directory.

# Specify the path to R working directory
workDir <- "/your-workDir-path-to/GKV"
setwd(workDir)

Load additional packages

The following packages are required: png, rasterVis. Check that these packages are installed on your system. If not, install them before loading them.

If you are working off-line, note that the rasterVis package depends on other packages that you also have to install: raster, lattice, sp, latticeExtra, RColorBrewer.

Load the technical functions

First, load the technical functions into R workspace as shown below. The script file fnct01-grayImage.R contains the code which reads, writes and plots a graylevel image. The script file fnct04-gkv-gray.R computes the illumination and the reflectance estimates. To this aim, it calls on the technical functions fnct02-gkv-utils.R and fnct03-gkv-grads.R.

# Run user-defined functions
source('./lib/fnct01-grayImage.R')
source('./lib/fnct04-gkv-gray.R')

Load and manipulate your image data in R workspace

Read, plot and write an image with default options

Let deep.png be the “.png” image we are using as an illustrated example.

To load an image into R workspace, use the function read.grayImage. The output of the function read.grayImage is a matrix containing the image pixel graylevel intensities. The intensities belong to the interval [0,1] where the value 0 stands for the black color and the value 1 stands for the white color. In our illustrated example, the matrix of pixel intensities is then stored in the object named as data.

# Run the read.grayImage function 
data <- read.grayImage('./data/deep.png')
# Run this if you want to check or need the image size
n1 <- dim(data)[1] # number of rows
n2 <- dim(data)[2] # number of columns
n1; n2

To plot this image (once it has been read), use the function plot.grayImage. By default, it is only necessary to specify the path to the image location from R working directory.

# Run the plot.grayImage function 
plot.grayImage(data)

Read, plot and write an image with user-specified options

To apply our method on this image, the image legend needs to be cut off. In fact, the region of interest of the image deep.png consists of all the pixel columns and the pixel rows from 1 to 422.

# Load data and overwrite the "data"" object
data <- read.grayImage('./data/deep.png', nrow=1:422, ncol=NULL)
# Run this if you want to check or need the size of the loaded part of the image 
n1 <- dim(data)[1] # number of rows
n2 <- dim(data)[2] # number of columns
n1; n2

# Run this if you want to plot the trimmed image
plot.grayImage(data)

More generally, a rectangular subregion of interest may be defined by the optional argument nrow= and ncol=. For instance, if the region of interest is a rectangle defined by the pixel rows from 10 to 30 and by the pixels columns from 50 to 70, then do:

# Run the read.grayImage function on a subimage
subdata <- read.grayImage('./data/deep.png', nrow=10:30, ncol=50:70)
plot.grayImage(subdata)

The object data may be written into a .png file with the function write.grayImage.

# Run this if you want to save the data matrix in the folder "cache" into a .png file
write.grayImage(data,'./cache/deep-orig.png')

The default setting of the functions plot.grayImage and write.grayImage is that the pixel intensities are encoded in 8 bits, which corresponds to 256 graylevels. Another number of graylevels may be specified through the optional argument cData=. The different possibilities are "bitmap" for binary images, "8bits" for 256 graylevel images (default), "10bits" for 1024 graylevel images and "12bits" for 4096 graylevel images.

The default setting of these functions also forces the pixel intensities to be in [0,1] in the following way. If some pixel intensities are negative, they are set to 0 (black). If some pixel intensities are larger than 1, they are set to 1 (white). An optional argument map="scaled" may be chosen. It forces the data to be in [0,1] by standardizing the graylevel intensity matrix.

Caution: R is not good at image displaying (as you can see on this page). That is why a recommended practice is that you write your image into a .png file and that you then open it with a suitable image viewer.

Plot the histogram of the image pixel intensities

The function hist.grayImage plots the histogram of the image pixel intensities. The default setting of this function is that the pixel intensities are encoded in 8 bits (that means 256 gray levels). Another number of graylevels may be specified through the optional argument cData=. The default setting of this function also forces the pixel intensities to be in [0,1] in the following way. If some pixel intensities are negative, they are set to 0 (black). If some pixel intensities are larger than 1, they are set to 1 (white). An optional argument map="scaled" may be chosen. It forces the data to be in [0,1] by standardizing the graylevel intensity matrix.

## Run this to plot the histogram of pixel intensities
hist.grayImage(data)

Reflectance estimation

To run our method and estimate an image reflectance, use the function illumination. The first input argument (data) is to be set as the matrix of pixel intensities. The second input argument (degree=) is the polynomial degree of the log-illumination function. The third input argument (method=) is the method chosen to estimate the gradient of the image. For the moment there are two choices for this latter argument. If the image is noise free, you can use the double finite difference method (methode="direct"). Otherwise you can use a local polynomial procedure (methode="gaussian") with a truncated Gaussian kernel.
The output of the function is a list consisting of different objects which are: the denoised image F, the reflectance estimation R, the illumination estimation L, the estimation of the polynomial coefficients of the log-illumination gamma, the parameters of the procedure pars, and the size of the image data.size.

# Run GKV estimation method
out <- illumination(data, degree=c(5,6), method='gaussian')
str(out)

Here the degree of the polynomial is set as (5,6). The parameter gamma is the matrix of polynomial coefficients of the log-illumination.

# Coefficients of the log-illumination polynome
round(out$gamma,2)

The gradient of the image is estimated by a local polynomial procedure (out$pars$gradient) with a truncated Gaussian kernel. The support of the truncated Gaussian kernel is of size 9 by 9 pixels. The bandwidth parameter is automatically chosen by cross-validation.

# Bandwidth parameter of local polynomials procedure which had been automatically chosen by the cross-validation procedure
out$pars$gradient$bandwidth

The three output images, that is the denoised image out$F, the estimated reflectance out$R and the estimated illumination out$L, can then be plotted.

Caution: do NOT use the optional argument map="scaled" when plotting the reflectance estimation. Indeed, our procedure is written to provide a consistent reflectance estimate the values of which are in [0,1]. If this is not the case, it indicates a convergence problem and and acts as a warning against non-reliable results.

# Plot the results:
# Plot the denoised image (but still with its illumination bias)
plot.grayImage(out$F)
# Plot the reflectance estimation
plot.grayImage(out$R)

Strictly speaking, the illumination image is not an image but a multiplicative factor of illumination effect. If the illumination at a given pixel is larger than 1, the ground truth intensity at this pixel is whiten by the illumination effect. If the illumination at a given pixel is less than 1, the ground truth intensity at this pixel is darken by the illumination effect. That is why we propose a level plot as illumination plot to get a visual impression of the repartition of the illumination effect.

# Plot the illumination estimation
plot.grayIllumination(out$L)

# Run this if you want to save the matrix out$R in the folder "figure" into a .png file
write.grayImage(out$R,'./figure/deep-R.png')

# Run this if you want to save the matrix out$L in the folder "figure" into a .png file
png(file="./figure/deep-L.png", width=ncol(out$L), height=(nrow(out$L)+50), res= 72)
plot.grayIllumination(out$L)
dev.off()

How to set the method parameters?

Estimation of the image gradient method=. If the image is noise free, both methods ("direct" and "gaussian") work and produce similar outputs. Note that the "gaussian" is more time-consuming because the smoothing parameter is chosen by cross-validation. If the image is noisy, the "gaussian" method provides results which are more consistant with the ground truth. The user may specify a pair of positive values for the smoothing parameter of the local polynomial procedure through the optional argument bandwidth= of the function illumination (in which case, no cross-validation procedure is done). The trick is that what we call the bandwidth is actually the inverse of bandwidth. Whatever, we advise the unexperienced user against the use of this optional argument.

# Run the GKV estimation method
out_direct <- illumination(data, degree=c(5,6), method='direct')
str(out_direct$pars)
plot.grayImage(out_direct$R)
plot.grayIllumination(out_direct$L)

Degree of the log-illumination polynome degree=. As of now, we do not propose any automatic procedure to select this argument. The user has to test several plausible values of pairs of positive integers. To this aim, we formulate some guidelines according to our experience. First, we highlight that the choice for this parameter is to be made within a (small) finite set. Indeed, if the degree of the polynome is chosen too large, the linear system that computes the parameter gamma of our procedure is singular. Secondly, an important point is the fact that a significant number of pixel intensity larger than 1 (sum(out$R>1)) Indeed, the reflectance to be recovered is supposed to be in [0,1], the estimate is proved to be consistent and no post-treatment is applied on the reflectance estimate out$R. Thus, if the estimated reflectance matrix gets a significant number of pixel intensity larger than 1 (sum(out$R>1)), this means that the choice of the argument degree is not well suited to the illumination artefact at hand.

# this is an example of the GKV estimation method run with an inconvenient degree choice
out_bad_degree <- illumination(data, degree=c(5,8), method='direct')

Our system is set in french hence the language of the error message. This error message says that the linear system to be solved is computationally singular and that its reciprocal condition number is far too small.

---

Authors: Geffray S., Klutchnikoff N., Vimond M.

Last modification: 09-09-2016