Color transfer between images

Beware of the fact that the function mean is different under Matlab and Scilab. It computes the total mean of a matrix under Scilab, whereas it returns the mean value along the first array dimension of the matrix under Matlab.

Write a function affine_transfer which apply an affine transform to \(u\) such that it has the same mean and the same standard deviation as \(v\) on each channel.

function w = affine_transfer(u,v)
  for i = 1:3
      INSERT YOUR CODE 
  end
endfunction

Display and comment the result.

w = affine_transfer(u,v);
w = min(max(w,0),255);             
ccview(w);                       // Matlab : imshow

Another solution consists in applying histogram specifications separately on each channel of the images \(u\) and \(v\).

function w = specification_channels(u,v)
  for i = 1:3
    uch = u(:,:,i);
    vch = v(:,:,i);
    INSERT YOUR CODE 
  end
endfunction

Display and comment the result.

w = specification_channels(u,v);
w = min(max(w,0),255);
ccview(w);                     // imshow

The RGB color space is not very well suited for these operations since the channels R, G and B are strongly correlated. Instead, we can first convert the RGB images to another color space with less correlation between channels. For instance, we can use a color opponent space with an achromatic channel and two chromatic channels representing an approximate red-green opponency and yellow- blue opponency :

$$ O_1 = \frac 1 {\sqrt{2}} (R-G),\; O_2 = \frac 1 {\sqrt{6}} (R+G-2B),\; O_3 = \frac 1 {\sqrt{3}} (R+G+B)$$

Apply the separable color transfers (affine and specification) to \(u\) and \(v\) in this opponent color space. Is the result more satisfying ?

O = [1/sqrt(3), 1/sqrt(3), 1/sqrt(3) ; 1/sqrt(2) -1/sqrt(2) 0; 1/sqrt(6) 1/sqrt(6) -2/sqrt(6)];
u_opp = matrix(matrix(u,nr*nc,3)*O',nr,nc,3);      // Matlab : replace 'matrix'  by 'reshape'
v_opp = matrix(matrix(v,nr*nc,3)*O',nr,nc,3);
//w_opp = affine_transfer(u_opp,v_opp);
w_opp = specification_channels(u_opp,v_opp);
w_opp = matrix(matrix(w_opp,nr*nc,3)*inv(O)',nr,nc,3);
w_opp = min(max(w_opp,0),255);
ccview(w_opp);                                     // Matlab  : imshow

As adviced in the works of Reinhard et al., these separable transformations can also be applied in the \(L\alpha\beta\) color space (see the following pdf).

LMS = [0.3811 0.5783 0.0402 ; 0.1967 0.7244 0.0782 ; 0.0241 0.1288 0.8444];
u_lms = log10(0.1+matrix(matrix(u,nr*nc,3)*LMS',nr,nc,3));
v_lms = log10(0.1+matrix(matrix(v,nr*nc,3)*LMS',nr,nc,3));
u_lab = matrix(matrix(u_lms,nr*nc,3)*O',nr,nc,3);
v_lab = matrix(matrix(v_lms,nr*nc,3)*O',nr,nc,3);
w_lab = affine_transfer(u_lab,v_lab);
//w_lab = specification_channels(u_lab,v_lab);
w_lms = matrix(matrix(w_lab,nr*nc,3)*inv(O)',nr,nc,3);
w_lms = matrix(matrix(exp(w_lms*log(10))-0.1,nr*nc,3)*inv(LMS)',nr,nc,3);

Display and comment the results.

w_lms = min(max(w_lms,0),255);
ccview(w_lms);

In order to transport the full color distribution of \(u\) on the one of \(v\), we propose to write these distributions as 3D point clouds \begin{equation} X = \begin{pmatrix} X_1 \\ \vdots\\ X_n \end{pmatrix}, \;\;\text{ and }\;\; Y = \begin{pmatrix} Y_1 \\ \vdots\\ Y_n \end{pmatrix}, \end{equation} with \(n\) the number of pixels in \(u\) and \(v\), \(X_i\) the color \((R_i,G_i,B_i)\) of the \(i\)-th pixel in \(u\) and \(Y_i\) the color of the \(i\)-th pixel in \(v\) (\(X\) and \(Y\) are \(n\times 3\) matrices). We look for the permutation \(\sigma\) of \(\{1,\dots, n\}\) which minimizes the quantity \begin{equation} \label{eq:transport_appariement} \sum_{i=1}^n \|X_i - Y_{\sigma(i)}\|^2_2. \end{equation} The assignment \(\sigma\) defines a mapping between the set of colors \(\{X_i\}_{i=1,\dots n}\) and the set of colors \(\{Y_j\}_{j=1,\dots n}\).

function [sX,sY] = transport1D(X,Y)
    [Xs,sX] = gsort(X,'g','i');   // We have Xs = X(sX)  // Matlab : sort(X)
    [Ys,sY] = gsort(Y,'g','i');
endfunction

function Z=transport3D(X,Y,N,e)  // X, Y and Z are  n x 3 matrices 
  Z = X;  
  for k=1:N

    // orthogonal basis, uniform distribution on the sphere
    u = rand(3,3,'normal');
    q = qr(u);

    // projection
    for j=1:3
      Zt = Z*q(:,j);
      Yt = Y*q(:,j);
     INSERT YOUR CODE
    end
  end
endfunction

X = rand(20,3);
Y = rand(20,3);
Z = transport3D(X,Y,5000,0.01);
param3d1(X(:,1),X(:,2),list(X(:,3),-1));  // points of X are represented as +
param3d1(Y(:,1),Y(:,2),list(Y(:,3),-2));  // points of Y are represented as x
param3d1(Z(:,1),Z(:,2),list(Z(:,3),-3));   // points of Y are represented as o 

Z=X;
figure;param3d1(X(:,1),X(:,2),list(X(:,3),-1));
param3d1(Y(:,1),Y(:,2),list(Y(:,3),-2));
param3d1(Z(:,1),Z(:,2),list(Z(:,3),-3));
f=gca(); // handle current graphic
sleep(1000); // pause
for i=1:60
   Z = transport3D(Z,Y,100,0.01);
   f.children(3).data=Z; // modify the coordinates of the third point cloud
   sleep(100);
end

u = u(1:4:nr,1:4:nc,:);
v = v(1:4:nr,1:4:nc,:);
[nr,nc,nch] = size(u);

X = matrix(u,nr*nc,3);
Y = matrix(v,nr*nc,3);

Z = transport3D(X,Y,1000,0.01);
w = matrix(Z,nr,nc,3);
w = min(max(w,0),255);
ccview(w);

Knowing the optimal (or an approximation) assignment between the point clouds \(X\) and \(Y\), you can compute barycenters between these discrete distributions. Use the previous results to compute the midway color distribution between \(u\) and \(v\) and the corresponding midway images.

A common drawback of classical methods aiming at color and contrast modifications is the revealing of artefacts (JPEG blocs, color inconsistancies, noise enhancement) or the attenuation of details and textures (see for instance the following web page). Let \(u\) be an image and \(g(u)\) the same image after color or contrast modification, we write \(\mathcal{M}(u) = g(u) - u\). All artefacts observable in \(g(u)\) can be seen as irregularities in the map \(\mathcal{M}(u)\). In order to reduce these artefacts, we propose to filter this difference map thanks to an operator \(Y_u\) and to reconstruct the image :

$$T(g(u)) = u + Y_u(g(u)-u).$$

Let us begin with a very simple averaging filter for \(Y_u\).

K = ones(4,4); K = K/sum(K);
for i=1:3
  wreg(:,:,i) = INSERT YOUR CODE HERE
end

Display the result

wreg = min(max(wreg,0),255);
ccview(wreg);

A more interesting result can be obtained by filtering the difference map \(\mathcal{M}(u)\) in order to follow the regularity of \(u\). To this aim, we define the operator \(Y_u\), a weighted average with weights depending on the similarity of pixels in the original image \(u\). The effect of this operator on an image \(u_0\) is defined as $$Y_u (u_0)(x) = \frac{1}{c(x)} \sum_{y \in \mathcal{V}(x)} u_0(y) w_u (x,y),$$ with \(\mathcal{V}(x)\) a square neighborhood of \(x\) and \(c(x)\) the normalization constant \(c(x) = \sum_{y \in \mathcal{V}(x)} w_u (x,y)\). The weights \(w_u (x, y)\) are defined as $$w_u (x, y) = e^{-\| u(x)−u(y) \|^2 / h^2 },$$ with \(\|.\|\) the Euclidean distance in \(\mathbb{R}^3\). and \(\mathcal{V}(x)\) a square neighborhood of \(x\). Observe that if we apply \(Y_u\) to the image \(u\), we obtain the Yaroslavsky filter.

• Write a function crossed_filter which computes the previously described filter.

function b=crossed_filter(u,u0,f,sigma)
    // the goal is to regularize u0 following the regularity of u 
    [nrow,ncol]=size(u);

    for i=1:nrow
        for j=1:ncol
         // Extract a neighborhood of (i,j) 
         iMin = max(i-f,1);
         iMax = min(i+f,nrow);
         jMin = max(j-f,1);
         jMax = min(j+f,ncol);
         v = u(iMin:iMax,jMin:jMax,:);
         v0 = u0(iMin:iMax,jMin:jMax,:);
         // Compute gaussian weights 
         INSERT YOUR CODE 
        // Compute the filtered image at (i,j) 
         INSERT YOUR CODE 
        end
    end
endfunction

Apply the previous filter to \(\mathcal{M}(u)\) and display the result \(T(g(u))\). In order to remove the aforementioned artefacts, you can iterate the filter.