Converter: added new color transfer modes: mkl, mkl-m, idt, idt-m

imagelib/color_transfer.py

@@ -1,9 +1,83 @@
 import cv2
 import numpy as np
 
+import scipy as sp
 import scipy.sparse
 from scipy.sparse.linalg import spsolve
 
+def color_transfer_mkl(x0, x1):
+    eps = np.finfo(float).eps
+    
+    h,w,c = x0.shape
+    h1,w1,c1 = x1.shape
+    
+    x0 = x0.reshape ( (h*w,c) )
+    x1 = x1.reshape ( (h1*w1,c1) )
+    
+    a = np.cov(x0.T)
+    b = np.cov(x1.T)
+
+    Da2, Ua = np.linalg.eig(a)
+    Da = np.diag(np.sqrt(Da2.clip(eps, None))) 
+
+    C = np.dot(np.dot(np.dot(np.dot(Da, Ua.T), b), Ua), Da)
+
+    Dc2, Uc = np.linalg.eig(C)
+    Dc = np.diag(np.sqrt(Dc2.clip(eps, None))) 
+
+    Da_inv = np.diag(1./(np.diag(Da)))
+
+    t = np.dot(np.dot(np.dot(np.dot(np.dot(np.dot(Ua, Da_inv), Uc), Dc), Uc.T), Da_inv), Ua.T) 
+
+    mx0 = np.mean(x0, axis=0)
+    mx1 = np.mean(x1, axis=0)
+
+    result = np.dot(x0-mx0, t) + mx1
+    return np.clip ( result.reshape ( (h,w,c) ), 0, 1)
+    
+def color_transfer_idt(i0, i1, bins=256, n_rot=20):
+    relaxation = 1 / n_rot
+    h,w,c = i0.shape
+    h1,w1,c1 = i1.shape
+    
+    i0 = i0.reshape ( (h*w,c) )
+    i1 = i1.reshape ( (h1*w1,c1) )
+    
+    n_dims = c
+    
+    d0 = i0.T
+    d1 = i1.T
+    
+    for i in range(n_rot):
+        
+        r = sp.stats.special_ortho_group.rvs(n_dims).astype(np.float32)
+        
+        d0r = np.dot(r, d0)
+        d1r = np.dot(r, d1)
+        d_r = np.empty_like(d0)
+        
+        for j in range(n_dims):
+            
+            lo = min(d0r[j].min(), d1r[j].min())
+            hi = max(d0r[j].max(), d1r[j].max())
+            
+            p0r, edges = np.histogram(d0r[j], bins=bins, range=[lo, hi])
+            p1r, _     = np.histogram(d1r[j], bins=bins, range=[lo, hi])
+
+            cp0r = p0r.cumsum().astype(np.float32)
+            cp0r /= cp0r[-1]
+
+            cp1r = p1r.cumsum().astype(np.float32)
+            cp1r /= cp1r[-1]
+            
+            f = np.interp(cp0r, cp1r, edges[1:])
+            
+            d_r[j] = np.interp(d0r[j], edges[1:], f, left=0, right=bins)
+        
+        d0 = relaxation * np.linalg.solve(r, (d_r - d0r)) + d0
+
+    return np.clip ( d0.T.reshape ( (h,w,c) ), 0, 1)
+
 def laplacian_matrix(n, m):
     mat_D = scipy.sparse.lil_matrix((m, m))
     mat_D.setdiag(-1, -1)