Dear Johannes,
Am 31.07.2017 um 04:18 schrieb johannes hanika:
oh, nice, thanks for looking into this!
i can't parse your diff, it has so many white space changes. could you
maybe fix it or explain what it does? is it just a special case for
the zero-patch case?
yeah the solver sucks. unfortunately eigen is written in c++, so i'm
hesitant to pull it in as a dependency.
in the color look up module we have to solve a linear system of size
(N+4)x(N+4) where N is the number of color patches. The problem is that
for N<=4 the coefficient matrix A has rank N, not the full rank N+4.
This means that the system has non-unique solutions. Thus we have to
pick a reasonable solution among the infinite set of solutions for N<=4.
I think the most reasonable solution is to set as many as possible
coefficients to zero. For the remaining coefficients, one can derive a
smaller full-rank linear system. This approach is numerically very
stable and explained in more detail in the python prototype attached to
this e-mail.
Assuming that all source colors are different, the (N+4)x(N+4)
coefficient matrix has full rank and can be solved by various methods.
LU-decomposition with partial pivoting but also numerical solution via
singular value decomposition as implemented in the python script work
quite well. I also tried various variants of Newton iteration, which
often also gives accurate results. In the Newton iteration approach,
one minimizes |A*x - b|, where x is the vector of unknowns and b the
right-hand-side of the linear system. For the Newton iteration, one has
to find an approximate solution to a linear system of the form 2*A^2 * y
= c. I tried various methods for solving the latter equation.
The modified version of colochecker.c in my github repository deals only
with the N=0 case. The diffs may show an unreasonable large change set
because my editor has automatically removed white spaces at the end of
lines.
Heiko
--
-- Number Crunch Blog @ https://www.numbercrunch.de
-- Cluster Computing @ http://www.clustercomputing.de
-- Professional @ https://www.mpi-hd.mpg.de/personalhomes/bauke
-- Social Networking @ https://www.researchgate.net/profile/Heiko_Bauke
___________________________________________________________________________
darktable developer mailing list
to unsubscribe send a mail to [email protected]#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# calculate interpolation coefficients as needed in
# the darktable color look up module
# import math and plotting modules
import numpy as np
from pylab import *
from scipy.optimize import minimize, curve_fit
from scipy.linalg import *
# the source colors
S=array([37.99, 13.56, 14.06, # dark skin
65.71, 18.13, 17.81, # light skin
49.93, -4.88, -21.93, # blue sky
43.14, -13.10, 21.91, # foliage
55.11, 8.84, -25.40, # blue flower
70.72, -33.40, -0.20, # bluish green
62.66, 36.07, 57.10, # orange
40.02, 10.41, -45.96, # purple red
51.12, 48.24, 16.25, # moderate red
30.33, 22.98, -21.59, # purple
72.53, -23.71, 57.26, # yellow green
71.94, 19.36, 67.86, # orange yellow
28.78, 14.18, -50.30, # blue
55.26, -38.34, 31.37, # green
42.10, 53.38, 28.19, # red
81.73, 4.04, 79.82, # yellow
51.94, 49.99, -14.57, # magenta
51.04, -28.63, -28.64, # cyan
96.54, -0.43, 1.19, # white
81.26, -0.64, -0.34, # neutral 8
66.77, -0.73, -0.50, # neutral 65
50.87, -0.15, -0.27, # neutral 5
35.66, -0.42, -1.23, # neutral 35
20.46, -0.08, -0.97 # black
])
S=reshape(S, (24, 3))
channel=1 # select a channel, (0, 1, 2) for (L, a, b)
# reduce number of color mappings for testing
# try N>4 for the nummerically instable cases
N=5
q=array(range(0, 24))
np.random.shuffle(q)
S=S[q[0:N], :] # random selection of original colors
# construct coefficient matrix
# the radial basis function
def kernel(a, b):
r2=dot(a-b, a-b)
return r2*log(r2) if r2>1e-12 else 0
A=zeros((N+4, N+4))
for i in range(0, N):
for j in range(0, N):
A[i, j]=kernel(S[i, :], S[j, :])
A[i, N+0]=A[N+0, i]=1
A[i, N+1]=A[N+1, i]=S[i, 0]
A[i, N+2]=A[N+2, i]=S[i, 1]
A[i, N+3]=A[N+3, i]=S[i, 2]
T=zeros((N+4))
T[0:N]=S[:, channel] + 10*randn(N) # target color channel, start from original color and modify
# quadratic cost function to be minimized (A must be symetric)
def N2(x):
return norm(dot(A, x)-T)
# gradient of cost function
def dN2(x):
return 2*dot(A, dot(A, x)-T)
# jacobian of cost function
def ddN2():
return 2*dot(A, A)
# solve the linear system, treat special cases N=0, N=1, N=2, N=3, N=4
# extra, solve cases N>4 by different methods
# for N<=4 the linear system has rank N, i.e. less than the dimension
# of the coefficient matrix the linear system, which is N+4
# we have to pick a reasonable solution among the infinite set of solutions for N<=4
if N==0:
# set all coefficents of the interpolating function to zero
x=zeros((N+4))
# set the coefficient of the linear function correponding to the target channel to one
x[N+1+channel]=1
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
elif N==1:
# set all coefficents of the interpolating function to zero
x=zeros((N+4))
# set the coefficient of the linear function correponding to the target channel to the ratio of target and source
x[N+1+channel]=T[0]/S[0, channel]
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
elif N==2:
# set all coefficents of the interpolating function to zero
x=zeros((N+4))
# determine the coefficient of the linear function correponding to the target channel and the coefficinet of the constant offset
x[N]=(T[0]*S[1, channel]-T[1]*S[0, channel])/(S[1, channel]-S[0, channel])
x[N+1+channel]=(T[1]-T[0])/(S[1, channel]-S[0, channel])
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
elif N==3:
# set all coefficents of the interpolating function to zero
x=zeros((N+4))
# determine the coefficient of the linear function correponding to
# the target channel, the coefficinet of the constant offset and
# the coefficients of the linear functions correponding to the
# other two color channels, the latter coefficinets are forced to
# be equal to render the solution unique
c1, c2={ 0, 1, 2 } - { channel }
M=zeros((3, 3))
M[0, 0]=1
M[1, 0]=1
M[2, 0]=1
M[0, 1]=S[0, channel]
M[1, 1]=S[1, channel]
M[2, 1]=S[2, channel]
M[0, 2]=S[0, c1]+S[0, c2]
M[1, 2]=S[1, c1]+S[1, c2]
M[2, 2]=S[2, c1]+S[2, c2]
a=solve(M, T[0:3])
x[3]=a[0]
x[4+channel]=a[1]
x[4+c1]=a[2]
x[4+c2]=a[2]
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
elif N==4:
# set all coefficents of the interpolating function to zero
x=zeros((N+4))
# determine the coefficient of the linear function correponding to
# the target channel, the coefficinet of the constant offset and
# the coefficients of the linear functions correponding to the
# other two color channels, the latter coefficinets may be different
M=zeros((4, 4))
M[0, 0]=1
M[1, 0]=1
M[2, 0]=1
M[3, 0]=1
M[:, 1:]=S[:, :]
a=solve(M, T[0:4])
x[4:]=a
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
else:
# solve directly
x=zeros((N+4))
try:
print("Direct solution via LU decomposition and partial pivoting:")
x=solve(A, T)
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
except:
print("Direct solution via LU cdecomposition failed.\n")
X=x
# solve by diagonalization of A
# is surprisingly inaccurate
x=zeros((N+4))
try:
print("Solution via diagonalization of A:")
Lambda, U=eigh(A)
x=dot(U, (1/Lambda)*dot(U.T, T))
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
except:
print("Solution via diagonalization of A failed.\n")
X=vstack((X, x))
# solve by sigular value decomposition of A
x=zeros((N+4))
try:
print("Solution via sigular value decomposition of A:")
u, s, v=svd(A, full_matrices=True)
x=dot(v.T, (1/s)*dot(u.T, T))
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
except:
print("Solution via sigular value decomposition of A failed.\n")
X=vstack((X, x))
# solve by miminizing quadratic cost function via Newton iteration
# solve 2*A^2 * y = b via LU factorization of 2*A^2 in each Newton iteration step
# in rare cases, this method gives results that differ from those as obtained
# by the other methods
x=zeros((N+4))
try:
print("Solution via Newton iteration and LU factorization of 2*A^2:")
a=polyfit(S[:, channel], T[0:N], 1)
x[N]=a[1]
x[N+1+channel]=a[0]
x[N]=1
x[N+1+channel]=0
p_L_U=lu_factor(ddN2())
for i in range(0, 128):
if N2(x)<1e-9:
break
y=dN2(x)
x-=lu_solve(p_L_U, y)
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
except:
print("Solution via Newton iteration and LU factorization of 2*A^2 failed.\n")
X=vstack((X, x))
# solve by miminizing quadratic cost function via Newton iteration
# solve 2*A^2 * y = b via diagonalization of 2*A^2 in each Newton iteration step
# in rare cases, this method gives results that differ from those as obtained
# by the other methods
x=zeros((N+4))
try:
print("Solution via Newton iteration and diagonalization 2*A^2:")
a=polyfit(S[:, channel], T[0:N], 1)
x[N]=a[1]
x[N+1+channel]=a[0]
x[N]=1
x[N+1+channel]=0
Lambda, U=eigh(ddN2())
for i in range(0, 128):
if N2(x)<1e-9:
break
y=dN2(x)
x-=dot(U, (1/Lambda)*dot(U.T, y))
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
except:
print("Solution via Newton iteration and diagonalization 2*A^2 failed.\n")
X=vstack((X, x))
# solve by miminizing quadratic cost function via Newton iteration
# solve 2*A^2 * y = b via singular value decomposition of 2*A^2 in each Newton iteration step
x=zeros((N+4))
try:
print("Solution via Newton iteration and singular value decomposition 2*A^2:")
a=polyfit(S[:, channel], T[0:N], 1)
x[N]=a[1]
x[N+1+channel]=a[0]
x[N]=1
x[N+1+channel]=0
u, s, v=svd(ddN2(), full_matrices=True)
for i in range(0, 128):
if N2(x)<1e-9:
break
y=dN2(x)
x-=dot(v.T, (1/s)*dot(u.T, y))
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
except:
print("Solution via Newton iteration and singular value decomposition 2*A^2 failed.\n")
X=vstack((X, x))
# solve by miminizing quadratic cost function via Newton iteration
# solve 2*A^2 * y = b via QR factorization of 2*A^2 in each Newton iteration step
x=zeros((N+4))
try:
print("Solution via Newton iteration and QR factorization of 2*A^2:")
a=polyfit(S[:, channel], T[0:N], 1)
x[N]=a[1]
x[N+1+channel]=a[0]
x[N]=1
x[N+1+channel]=0
q, r=qr(ddN2())
for i in range(0, 128):
if N2(x)<1e-9:
break
y=dN2(x)
x-=solve_triangular(r, dot(q.T, y))
print('residual: {0}\nsolution: {1}\n'.format(N2(x), x))
except:
print("Solution via Newton iteration and QR factorization of 2*A^2 failed.\n")
X=vstack((X, x))
# compare results as obtained by the differen methods
for j in range(0, X.shape[0]):
for i in range(0, j):
print("{} {} {}".format(i, j, norm(X[i, :]-X[j, :])))
