论文复现：Colorization Using Optimization_colorization using optimization中文翻译-CSDN博客

本文链接：https://blog.csdn.net/weixin_44178062/article/details/127833326

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

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

RGB格式转化为YUV格式：
$\times R + 0.587 \times G + 0.114 \times B$
$\times R - 0.289 \times G + 0.436 \times B$
$\times R - 0.515 \times G - 0.100 \times B$
YUV格式转化为RGB格式：
$\times V$
$\times U - 0.581 \times V$
$\times U$

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

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

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

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

用一个例子来描述算法的大致过程，如图所示是一个 $3 * 3$ 的图像， $2 、 6 、 7$ 处已着色，我们要做的工作就是求解出其他每一个未知的 $U V$ 值。
在这里插入图片描述
用 $x_{1}$ , $x_{2}$ , $x_{3}$ , $x_{4}$ , $x_{5}$ , $x_{6}$ , $x_{7}$ , $x_{8}$ , $x_{9}$ 来表示每一处要求解的U值。
对于第一个像素点，遍历其邻域坐标计算出权重 $w_{12}$ , $w_{14}$ , $w_{15}$ ，有： $x_{1}-w_{12}x_{2}-w_{14}x_{4}-w_{15}x_{5}=0$
对于第二个像素点，该像素点已知，则有： $x_{2}=u_{2}$
对于第三个像素点，遍历其邻域坐标计算出权重 $w_{32}$ , $w_{35}$ , $w_{36}$ ，有： $x_{3}-w_{32}x_{2}-w_{35}x_{5}-w_{36}x_{6}=0$
对于第四个像素点，遍历其邻域坐标计算出权重 $w_{41}$ , $w_{42}$ , $w_{45}$ , $w_{47}$ , $w_{48}$ ，有： $x_{4}-w_{41}x_{1}-w_{42}x_{2}-w_{45}x_{5}-w_{47}x_{7}-w_{48}x_{8}=0$
对于第五个像素点，遍历其邻域坐标计算出权重 $w_{51}$ , $w_{52}$ , $w_{53}$ , $w_{54}$ , $w_{56}$ , $w_{57}$ , $w_{58}$ , $w_{58}$ 有： $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_{6}=u_{6}$
对于第七个像素点，该像素点已知，则有： $x_{7}=u_{7}$
对于第八个像素点，遍历其邻域坐标计算出权重 $w_{84}$ , $w_{85}$ , $w_{86}$ , $w_{87}$ , $w_{89}$ ，有： $x_{8}-w_{84}x_{4}-w_{85}x_{5}-w_{86}x_{6}-w_{87}x_{7}-w_{89}x_{9}=0$
对于第九个像素点，遍历其邻域坐标计算出权重 $w_{95}$ , $w_{96}$ , $w_{98}$ ，有： $x_{9}-w_{95}x_{5}-w_{96}x_{6}-w_{98}x_{8}=0$
联立以上方程组，移项并写成矩阵形式：
$\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}$

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

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

class Colorization {
private:
	Image gray_img;		//输入灰度图
	Image sketch_img;	//部分着色图
	Mat mask;     //差异图像
	Image dest_img;		//输出彩色图
	SparseMatrix<double> W; //稀疏矩阵 D-W
	VectorXd u;				//列向量B
	VectorXd 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 for
		for (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、V
		VectorXd 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 known
							u_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;
		
		solver.compute(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$ 值是不同的，计算权重矩阵和对目标图像赋值时，应该使用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);