论文复现：Colorization Using Optimization

论文复现：Colorization Using Optimization"/>

论文复现：Colorization Using Optimization

Colorization Using Optimization这篇论文中介绍了一种简单而有效的灰度图像着色算法：在 Y U V YUV YUV色彩空间中，对灰度图像进行简单着色，再求解其他未知的像素点，填充到目标图像中得到彩色图像。

Y U V YUV YUV色彩空间中， Y Y Y表示图像的灰度， U 、 V U、V U、V表示图像的色度。灰度图像在 Y U V YUV YUV色彩空间中 Y Y Y值已知而 U 、 V U、V U、V值未知，经过简单着色后的图像 Y Y Y值已知且部分 U 、 V U、V U、V值未知，最后得到的目标图像 Y 、 U 、 V Y、U、V Y、U、V值均已知。

RGB格式转化为YUV格式：
Y = 0.299 × R + 0.587 × G + 0.114 × B Y = 0.299 \times R + 0.587 \times G + 0.114 \times B Y=0.299×R+0.587×G+0.114×B
U = − 0.147 × R − 0.289 × G + 0.436 × B U = -0.147 \times R - 0.289 \times G + 0.436 \times B U=−0.147×R−0.289×G+0.436×B
V = 0.615 × R − 0.515 × G − 0.100 × B V = 0.615 \times R - 0.515 \times G - 0.100 \times B V=0.615×R−0.515×G−0.100×B
YUV格式转化为RGB格式：
R = Y + 1.140 × V R = Y + 1.140 \times V R=Y+1.140×V
G = Y − 0.395 × U − 0.581 × V G = Y - 0.395 \times U - 0.581 \times V G=Y−0.395×U−0.581×V
B = Y + 2.032 × U B = Y + 2.032 \times U B=Y+2.032×U

算法的输入：灰度图像和简单着色后的图像
算法的输出：彩色图像

算法基于一个前提： Y U V YUV YUV色彩空间下， Y Y Y值相似的相邻像素点，其 U V UV UV值也应该相似。
如何来描述这种相似呢？论文中用到了成本函数 J ( U ) J(U) J(U)和 J ( V ) J(V) J(V)：
J ( U ) = ∑ r ( U r − ∑ s ∈ N ( r ) w r s U s ) 2 J(U)=\sum_{r}(U_{r}-\sum_{s\in N(r)}w_{rs}U_{s})^2 J(U)=r∑​(Ur​−s∈N(r)∑​wrs​Us​)2
J ( V ) = ∑ r ( V r − ∑ s ∈ N ( r ) w r s V s ) 2 J(V)=\sum_{r}(V_{r}-\sum_{s\in N(r)}w_{rs}V_{s})^2 J(V)=r∑​(Vr​−s∈N(r)∑​wrs​Vs​)2

这样就把问题转化为最小化成本函数的优化问题。直观来看，当括号中的内容等于 0 0 0时，成本函数能取得最小值 0 0 0。
U r − ∑ s ∈ N ( r ) w r s U s = 0 U_{r}- \sum_{s\in N(r)} w_{rs}U_{s}=0 Ur​−s∈N(r)∑​wrs​Us​=0
V r − ∑ s ∈ N ( r ) w r s V s = 0 V_{r}- \sum_{s\in N(r)} w_{rs}V_{s}=0 Vr​−s∈N(r)∑​wrs​Vs​=0

其中权重函数可以通过下式计算（每一个像素点与其邻域像素点的权重值需要有归一化的过程）。
W r s ∝ e − ( Y r − Y s ) 2 2 σ r 2 W_{rs}\propto e^{-\frac{(Y_{r}-Y_{s})^2}{2\sigma _{r}^2 } } Wrs​∝e−2σr2​(Yr​−Ys​)2​

用一个例子来描述算法的大致过程，如图所示是一个 3 ∗ 3 3*3 3∗3的图像， 2 、 6 、 7 2、6、7 2、6、7处已着色，我们要做的工作就是求解出其他每一个未知的 U V UV UV值。

用 x 1 x_{1} x1​, x 2 x_{2} x2​, x 3 x_{3} x3​, x 4 x_{4} x4​, x 5 x_{5} x5​, x 6 x_{6} x6​, x 7 x_{7} x7​, x 8 x_{8} x8​, x 9 x_{9} x9​来表示每一处要求解的U值。
对于第一个像素点，遍历其邻域坐标计算出权重 w 12 w_{12} w12​, w 14 w_{14} w14​, w 15 w_{15} w15​，有： x 1 − w 12 x 2 − w 14 x 4 − w 15 x 5 = 0 x_{1}-w_{12}x_{2}-w_{14}x_{4}-w_{15}x_{5}=0 x1​−w12​x2​−w14​x4​−w15​x5​=0
对于第二个像素点，该像素点已知，则有： x 2 = u 2 x_{2}=u_{2} x2​=u2​
对于第三个像素点，遍历其邻域坐标计算出权重 w 32 w_{32} w32​, w 35 w_{35} w35​, w 36 w_{36} w36​，有： x 3 − w 32 x 2 − w 35 x 5 − w 36 x 6 = 0 x_{3}-w_{32}x_{2}-w_{35}x_{5}-w_{36}x_{6}=0 x3​−w32​x2​−w35​x5​−w36​x6​=0
对于第四个像素点，遍历其邻域坐标计算出权重 w 41 w_{41} w41​, w 42 w_{42} w42​, w 45 w_{45} w45​, w 47 w_{47} w47​, w 48 w_{48} w48​，有： x 4 − w 41 x 1 − w 42 x 2 − w 45 x 5 − w 47 x 7 − w 48 x 8 = 0 x_{4}-w_{41}x_{1}-w_{42}x_{2}-w_{45}x_{5}-w_{47}x_{7}-w_{48}x_{8}=0 x4​−w41​x1​−w42​x2​−w45​x5​−w47​x7​−w48​x8​=0
对于第五个像素点，遍历其邻域坐标计算出权重 w 51 w_{51} w51​, w 52 w_{52} w52​, w 53 w_{53} w53​, w 54 w_{54} w54​, w 56 w_{56} w56​, w 57 w_{57} w57​, w 58 w_{58} w58​, w 58 w_{58} w58​有： x 5 − w 51 x 1 − w 52 x 2 − w 53 x 3 − w 54 x 4 − w 56 x 6 − w 57 x 7 − w 58 x 8 − w 59 x 9 = 0 x_{5}-w_{51}x_{1}-w_{52}x_{2}-w_{53}x_{3}-w_{54}x_{4}-w_{56}x_{6}-w_{57}x_{7}-w_{58}x_{8}-w_{59}x_{9}=0 x5​−w51​x1​−w52​x2​−w53​x3​−w54​x4​−w56​x6​−w57​x7​−w58​x8​−w59​x9​=0
对于第六个像素点，该像素点已知，则有： x 6 = u 6 x_{6}=u_{6} x6​=u6​
对于第七个像素点，该像素点已知，则有： x 7 = u 7 x_{7}=u_{7} x7​=u7​
对于第八个像素点，遍历其邻域坐标计算出权重 w 84 w_{84} w84​, w 85 w_{85} w85​, w 86 w_{86} w86​, w 87 w_{87} w87​, w 89 w_{89} w89​，有： x 8 − w 84 x 4 − w 85 x 5 − w 86 x 6 − w 87 x 7 − w 89 x 9 = 0 x_{8}-w_{84}x_{4}-w_{85}x_{5}-w_{86}x_{6}-w_{87}x_{7}-w_{89}x_{9}=0 x8​−w84​x4​−w85​x5​−w86​x6​−w87​x7​−w89​x9​=0
对于第九个像素点，遍历其邻域坐标计算出权重 w 95 w_{95} w95​, w 96 w_{96} w96​, w 98 w_{98} w98​，有： x 9 − w 95 x 5 − w 96 x 6 − w 98 x 8 = 0 x_{9}-w_{95}x_{5}-w_{96}x_{6}-w_{98}x_{8}=0 x9​−w95​x5​−w96​x6​−w98​x8​=0
联立以上方程组，移项并写成矩阵形式：
( 1 0 0 − w 14 − w 15 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 − w 35 0 0 0 0 − w 41 0 0 1 − w 45 0 0 − w 48 0 − w 51 0 − w 53 − w 54 1 0 0 − w 58 − w 59 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 − w 84 − w 85 0 0 1 − w 89 0 0 0 0 − w 95 0 0 − w 98 1 ) ( x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 8 x 9 ) = ( w 12 u 2 u 2 w 32 u 2 + w 36 u 6 w 42 u 2 + w 47 u 7 w 52 u 2 + w 57 u 7 u 6 u 7 w 86 u 6 + w 87 u 7 w 9 6 u 6 ) \begin{pmatrix} 1 & 0 & 0 & -w_{14} & -w_{15} & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & -w_{35} & 0 & 0 & 0 & 0 \\ -w_{41} & 0 & 0 & 1 & -w_{45} & 0 & 0 & -w_{48} & 0 \\ -w_{51} & 0 & -w_{53} & -w_{54} & 1 & 0 & 0 & -w_{58} & -w_{59} \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & -w_{84} & -w_{85} & 0 & 0 & 1 & -w_{89} \\ 0 & 0 & 0 & 0 & -w_{95} & 0 & 0 & -w_{98} & 1 \end{pmatrix}\begin{pmatrix} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ x_{5} \\ x_{6} \\ x_{7} \\ x_{8} \\ x_{9} \end{pmatrix}=\begin{pmatrix} w_{12}u_{2} \\u_{2} \\w_{32}u_{2}+w_{36}u_{6} \\w_{42}u_{2}+w_{47}u_{7} \\w_{52}u_{2}+w_{57}u_{7} \\u_{6} \\u_{7} \\w_{86}u_{6}+w_{87}u_{7} \\w_{9{\tiny } 6}u_{6} \end{pmatrix} ⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛​100−w41​−w51​0000​010000000​0010−w53​0000​−w14​001−w54​00−w84​0​−w15​0−w35​−w45​100−w85​−w95​​000001000​000000100​000−w48​−w58​001−w98​​0000−w59​00−w89​1​⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞​⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛​x1​x2​x3​x4​x5​x6​x7​x8​x9​​⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞​=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛​w12​u2​u2​w32​u2​+w36​u6​w42​u2​+w47​u7​w52​u2​+w57​u7​u6​u7​w86​u6​+w87​u7​w96​u6​​⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞​

接下来的问题就是求解以上线性方程组，得到其他未知像素点的 U U U值（ V V V值的求解同理）。

使用VS+opencv+eigen来实现该算法。通过以上的例子可以知道，算法的主要工作是构建一个线性系统（ W x = u , W x = v Wx=u,Wx=v Wx=u,Wx=v）。

class Colorization {
private:Image gray_img;		//输入灰度图Image sketch_img;	//部分着色图Mat mask;     //差异图像Image dest_img;		//输出彩色图SparseMatrix<double> W; //稀疏矩阵 D-WVectorXd u;				//列向量BVectorXd v;Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> > solver; //线性方程组求解public:Colorization() {}//构造函数Colorization(Image& g_img, Image& s_img, Mat m_img, Image& d_img) {gray_img = g_img;sketch_img = s_img;mask = m_img;dest_img = d_img;}void GetColorfulImage() {W.resize(sketch_img.h * sketch_img.w, sketch_img.h * sketch_img.w);W.setIdentity();double sigma = 10000;//double mean = 0.0;//均值//double sum = 0.0;//for (int i = 0; i < sketch_img.h; i++) {//	for (int j = 0; j < sketch_img.w; j++) {//		double y = sketch_img.yuv_colors[i][j][0];//		sum += y;//	}//}//mean = sum / (sketch_img.h * sketch_img.w);//double sigma = 0.0;//方差//for (int i = 0; i < sketch_img.h; i++) {//	for (int j = 0; j < sketch_img.w; j++) {//		sigma = sigma + (sketch_img.yuv_colors[i][j][0] - mean) * (sketch_img.yuv_colors[i][j][0] - mean);//	}//}//sigma = sigma / (sketch_img.h * sketch_img.w);cout << "build laplacian matrix..." << endl; //计算权重矩阵（有颜色的像素点直接跳过不用计算它与邻域的权重）
//#pragma omp parallel forfor (int i = 0; i < gray_img.h; i++) {for (int j = 0; j < gray_img.w; j++) {int index = i * mask.cols + j;double data = mask.data[index];if (data == 255) continue;//如果像素点已知int center_index = i * gray_img.w + j;//像素点在矩阵中的行值double Y_ij = gray_img.yuv_colors[i][j][0];//像素点的Y值vector<int> neighbors = gray_img.GetNeighborsOf(i,j);//获取像素点邻域索引vector<double> weights;//记录权重clock_t t1 = clock();double w_sum = 0;//所有邻域像素点与中心像素点的权重值之和for (int k = 0; k < neighbors.size(); k++) {int neighbor_index = neighbors[k];int x = neighbor_index / gray_img.w;//计算邻域像素点在图像中的行、列值int y = neighbor_index % gray_img.w;double Y_ij_neighbor = gray_img.yuv_colors[x][y][0];//邻域像素点double w = exp(-pow(Y_ij - Y_ij_neighbor, 2) / (2 * sigma ));weights.push_back(w);w_sum += w;}clock_t t2 = clock();//权重值归一化for (int k = 0; k < neighbors.size(); k++) {int neighbor_index = neighbors[k];double neighbor_w = weights[k] / w_sum;W.coeffRef(center_index, neighbor_index) = - neighbor_w; //这里耗时间, 没什么好办法。。。那就等它执行？嗯，一行要1s}clock_t t3 = clock();//cout << i << "," << j << ", time 1: " << t2 - t1 << "; time 2: " << t3 - t2 << endl;}}cout << "build right hand matrix U and V" << endl;//构造U、VVectorXd u; u.resize(sketch_img.h * sketch_img.w); u.setZero();VectorXd v; v.resize(sketch_img.h * sketch_img.w); v.setZero();for (int i = 0; i < sketch_img.h; i++) {for (int j = 0; j < sketch_img.w; j++) {int index = i * mask.cols + j;double data = mask.data[index];if (data == 255) {u(i * sketch_img.w + j) = sketch_img.yuv_colors[i][j][1];v(i * sketch_img.w + j) = sketch_img.yuv_colors[i][j][2];}else {vector<int> ij_neighbor = sketch_img.GetNeighborsOf(i, j);for (int k = 0; k < ij_neighbor.size(); k++) {int ij_neighbor_index = ij_neighbor[k];int x = ij_neighbor_index / sketch_img.w;int y = ij_neighbor_index % sketch_img.w;double u_tmp = 0;double v_tmp = 0;double neighbor_data = mask.data[ij_neighbor_index];if (neighbor_data == 255) { //if neighbor is knownu_tmp = -W.coeffRef(i * sketch_img.w + j, ij_neighbor[k]) * sketch_img.yuv_colors[x][y][1];v_tmp = -W.coeffRef(i * sketch_img.w + j, ij_neighbor[k]) * sketch_img.yuv_colors[x][y][2];u(i * sketch_img.w + j) += u_tmp;v(i * sketch_img.w + j) += v_tmp;W.coeffRef(i * sketch_img.w + j, ij_neighbor[k]) = 0;}}}}}cout << "solve linear equations..." << endl;solverpute(W);VectorXd U = solver.solve(u);VectorXd V = solver.solve(v);//TODO:操作dest-img的YUV通道数据dest_img.ConvertRGB2YUV();for (int i = 0; i < sketch_img.h; i++)for (int j = 0; j < sketch_img.w; j++) {dest_img.yuv_colors[i][j][0] = gray_img.yuv_colors[i][j][0];dest_img.yuv_colors[i][j][1] = U(i * sketch_img.w + j);dest_img.yuv_colors[i][j][2] = V(i * sketch_img.w + j);}dest_img.ConvertYUV2RGB();}//保存结果void SaveResults(string path) {dest_img.Save(path);}
};

效果图

总结

• gray_img和sketch_img中 Y Y Y值是不同的，计算权重矩阵和对目标图像赋值时，应该使用gray_img.yuv_colors[i][j][0]
• 先想好算法思路再写代码，暴力解写出来的代码又乱又容易错。。。
• 使用灰度图和涂鸦图的差值图像来判定涂鸦图着色区域
• 算法中计算的U、V值转化为RGB值时，要保证RGB值在0~255，不然会出现很奇怪的颜色区域

//保证RGB值在0~255
double new_r = min(max(R, 0.0), 255.0);
double new_g = min(max(G, 0.0), 255.0);
double new_b = min(max(B, 0.0), 255.0);
rgb_colors[i][j] = Vec3d(new_r, new_g, new_b);