# Multivariate Histogram Matching for Image Color Matching

Suppose X is random variable that is spatial position in finite range, which its Cumulative Distribution Function (CDF) could be calculated easily, then if I have two random variable X and Y spatial position, can I make multivariate CDF by assuming those pair as one dimensional vector in the pair possible finite range?

## TL;DR

I am working on experiment to match 2 images color hues, inspired from histogram matching by using CDF of each image, and create lookup table (LUT) from it. The naive way in each image, one channel used as the random variable (X), when x is the color range [0..255], so CDF will satisfy f(x) = P(X <= x), and to match both function is to map f_src(x) to f_dst(x).

What I am trying to do is to transform source image distribution in the 2d hue space to match target image distribution, by taking Cb and Cr component of each image, treat it as position in CbCr space (p,q), getting its CDF, and do histogram matching. I attempted to treat (p,q) as flattened points (p+q*256), that is element of 'CrCb range' [0..65536], build CDF, and find the LUT.

source = np.asarray(Image.open('./misa_2.png').resize((640,480)).convert('YCbCr'))
target = np.asarray(Image.open('./misa_3.png').resize((640,480)).convert('YCbCr'))

# source CbCr as flattened pair of points
flatpt_src = source[:,:,1:].reshape(-1,source[:,:,1:].shape[-1]).swapaxes(1,0).copy()
flatpt_trg = np.concatenate((flatpt_trg, np.array([[0,0],[255,255]]).T), axis=1)
flatpt_src.sort(axis=1)

# Cb -> x, Cr -> y
x_src_flat = (flatpt_src[1,:] + flatpt_src[0,:]*256) # how [1,255] will be like if flattened
y_src_cdf = np.arange(1, x_src_flat.shape[-1] + 1) / x_src_flat.shape[-1] # 1/len ... len/len

# target CbCr as flattened pair of points
flatpt_trg = target[:,:,1:].reshape(-1,target[:,:,1:].shape[-1]).swapaxes(1,0).copy()
flatpt_trg.sort(axis=1)

# Cb -> x, Cr -> y
x_trg_flat = (flatpt_trg[1,:] + flatpt_trg[0,:]*256) # how [1,255] will be like if flattened
y_trg_cdf = np.arange(1, x_trg_flat.shape[-1] + 1) / x_trg_flat.shape[-1] # 1/len ... len/len

# Source CbCr CDF, Target CbCr CDF
plt.plot(x_src_flat,y_src_cdf) # blue
plt.plot(x_trg_flat,y_trg_cdf) # orange lookup table from projecting cdf source to cdf target

# check if lut plot shape is same as 'cdf_src(x) = cdf_trg(x)' mapping
plt.plot(flatpt_lut_index, lut2D) # blue
plt.plot(ka, el) # orange when LUT is applied to source, this is the result

result = np.dstack((source.copy()[:,:,0], CbCr_res))

axshows = (source[:,:,1], source[:,:,2], result[:,:,1], result[:,:,2])
axtitles = ('source Cb', 'source Cr', 'target Cb', 'target Cr')

f, axs = plt.subplots(2,2, figsize=(12.8,7.2), dpi=100)
f.tight_layout()
for ax, ashw, attl in zip(axs.flat, axshows, axtitles):
ax.imshow(ashw, cmap='gray')
ax.set_title(attl) the result is clearly what I didn't expecting, and at this point I am confused whether which is wrong, is it my 'multivariate' workaround or just the whole implementation.

(edit : 01/28 : update CDF plot and CbCr result-source comparison)