立体匹配中的代价
AD (颜色差异的绝对值):
- 公式(1)(颜色差异的绝对值):
- AD变换反映像素点的灰度变化,在纹理丰富区域具有良好的匹配效果,是一种简单、易实现的代价衡量的方法。但是,基于单个像素点计算的匹配代价往往会受到图像噪声、光照不均等的影响,相似度可靠性不高。
- 公式(2)(梯度差异的绝对值):
- 参考论文:
- a local adaptive approach for dense stereo matching in architectural scene reconstruction(2013)
census变换
在1994年,Zabih和Woodfill 提出了Census transform,后来被广泛用于计算机视觉领域中的匹配代价计算。Census变换是在图像区域定义一个矩形窗口,用这个矩形窗口遍历整幅图像。选取中心像素作为参考像素,将矩形窗口中每个像素的灰度值与参考像素的灰度值进行比较,灰度值小于或等于参考值的像素标记为0,大于参考值的像素标记为1,最后再将它们按位连接,得到变换后的结果,变换后的结果是由0和1组成的二进制码流。如下图所示。
Census变换的匹配代价计算方法是计算左右影像对应的两个像素的Census变换值的汉明(Hamming)距离。
- 优点:Census变换对图像的明暗变化不敏感,能够容忍一定的噪声,因为比较的是相对灰度关系,所以即使左右影像亮度不一致,也能得到较好的匹配效果。而且,相比互信息具有并行度高的优点,Census变换值是局部窗口运算,每个像素可以独立运算,这个特性让其可以很好的设计多线程并行计算模型,无论是CPU并行还是GPU并行都能达到非常高的并行效率。
- 缺点:census变换在双目匹配中对于弱纹理区域和重复场景的匹配效果不好。
例如:灰度完全不同的两个邻域窗口,Census变换得到的代价相同。
- 参考文献:
- ZABIH R, WOODFILL J. Non-parametric local transforms for computing visual correspondence[M]. 1994: 151-158
- a local adaptive approach for dense stereo matching in architectural scene reconstruction(2013)
SAD-ALD:(颜色差异的绝对值之和+臂长差异)
这种代价衡量主要假设是:左右图像对应的像素具有相似的颜色信息,同时在垂直方向上支撑域具有相似的臂长。
- SAD:
- SAD+臂长:
- 优点:
在大多数区域可以减少错误,但是在颜色和形状重复的区域,这种优势会减弱。
- 参考文献:
- A Near Real-Time Color Stereo Matching Method for GPU
- SSD(Sum of Squared Difference):像素差的平方和
- 公式:
- 实现:
// Sum of Squred Difference Mat ssd(Mat left, Mat right, string type, int kernel_size=3,int max_offset = 79) { //图像的宽和高 int width = left.size().width; int height = left.size().height; Mat depth(height, width, 0); vector< vector<int> > min_ssd; // store min SSD values for (int i = 0; i < height; ++i) { vector<int> tmp(width, numeric_limits<int>::max()); min_ssd.push_back(tmp); } //生成灰度图像 Mat left_gray,right_gray; cvtColor(left,left_gray,CV_BGR2GRAY); cvtColor(right,right_gray,CV_BGR2GRAY); for (int offset = 0; offset <= max_offset; offset++) { Mat tmp(height, width, 0); // shift image depend on type to save calculation time //产生视差 if (type == "left") { for (int y = 0; y < height; y++) { for (int x = 0; x < offset; x++) { tmp.at<uchar>(y, x) = right.at<uchar>(y, x); } for (int x = offset; x < width; x++) { tmp.at<uchar>(y, x) = right.at<uchar>(y, x - offset); } } } else if (type == "right") { for (int y = 0; y < height; y++) { for (int x = 0; x < width - offset; x++) { tmp.at<uchar>(y, x) = left.at<uchar>(y, x + offset); } for (int x = width - offset; x < width; x++) { tmp.at<uchar>(y, x) = left.at<uchar>(y, x); } } } else { Mat tmp(0, 0, 0); return tmp; } // calculate each pixel's SSD value for (int y = 0; y < height; y++) { for (int x = 0; x < width; x++) { int start_x = max(0, x - kernel_size); int start_y = max(0, y - kernel_size); int end_x = min(width - 1, x + kernel_size); int end_y = min(height - 1, y + kernel_size); int sum_sd = 0; if (type == "left") { for (int i = start_y; i <= end_y; i++) { for (int j = start_x; j <= end_x; j++) { int delta = abs(left.at<uchar>(i, j) - tmp.at<uchar>(i, j)); sum_sd += delta * delta; } } } else { for (int i = start_y; i <= end_y; i++) { for (int j = start_x; j <= end_x; j++) { int delta = abs(right.at<uchar>(i, j) - tmp.at<uchar>(i, j)); sum_sd += delta * delta; } } } // smaller SSD value found if (sum_sd < min_ssd[y][x]) { min_ssd[y][x] = sum_sd; // for better visualization depth.at<uchar>(y, x) = (uchar)(offset * 3); } } } } return depth; }
- CC: 领域的互相关系数
- NCC :归一化互相关系数(归一化:消除对光照敏感的问题)
- 公式:
- 实现:
Mat ncc(Mat in1, Mat in2, string type, int max_disparity = 79,int kernel_size = 3) { int width = in1.size().width; int height = in1.size().height; Mat left,right; cvtColor(in1,left,CV_BGR2GRAY); cvtColor(in2,right,CV_BGR2GRAY); Mat depth(height, width, 0); vector< vector<double> > max_ncc; // store max NCC value for (int i = 0; i < height; ++i) { vector<double> tmp(width, -2); max_ncc.push_back(tmp); } for (int offset = 1; offset <= max_disparity; offset++) { Mat tmp(height, width, 0); // shift image depend on type to save calculation time if (type == "left") { for (int y = 0; y < height; y++) { for (int x = 0; x < offset; x++) { tmp.at<uchar>(y, x) = right.at<uchar>(y, x); } for (int x = offset; x < width; x++) { tmp.at<uchar>(y, x) = right.at<uchar>(y, x - offset); } } } else if (type == "right") { for (int y = 0; y < height; y++) { for (int x = 0; x < width - offset; x++) { tmp.at<uchar>(y, x) = left.at<uchar>(y, x + offset); } for (int x = width - offset; x < width; x++) { tmp.at<uchar>(y, x) = left.at<uchar>(y, x); } } } else { Mat tmp(0, 0, 0); return tmp; } // calculate each pixel's NCC value for (int y = 0; y < height; y++) { for (int x = 0; x < width; x++) { int start_x = max(0, x - kernel_size); int start_y = max(0, y - kernel_size); int end_x = min(width - 1, x + kernel_size); int end_y = min(height - 1, y + kernel_size); double n = (end_y - start_y) * (end_x - start_x); double res_ncc = 0; if (type == "left") { double left_mean = 0, right_mean = 0; double left_std = 0, right_std = 0; double numerator = 0; for (int i = start_y; i <= end_y; i++) { for (int j = start_x; j <= end_x; j++) { left_mean += left.at<uchar>(i, j); right_mean += tmp.at<uchar>(i, j); } } left_mean /= n; right_mean /= n; for (int i = start_y; i <= end_y; i++) { for (int j = start_x; j <= end_x; j++) { left_std += pow(left.at<uchar>(i, j) - left_mean, 2); right_std += pow(tmp.at<uchar>(i, j) - right_mean, 2); numerator += (left.at<uchar>(i, j) - left_mean) * (tmp.at<uchar>(i, j) - right_mean); } } numerator /= n; left_std /= n; right_std /= n; res_ncc = numerator / (sqrt(left_std) * sqrt(right_std)) / n; } else { double left_mean = 0, right_mean = 0; double left_std = 0, right_std = 0; double numerator = 0; for (int i = start_y; i <= end_y; i++) { for (int j = start_x; j <= end_x; j++) { left_mean += tmp.at<uchar>(i, j); right_mean += right.at<uchar>(i, j); } } left_mean /= n; right_mean /= n; for (int i = start_y; i <= end_y; i++) { for (int j = start_x; j <= end_x; j++) { left_std += pow(tmp.at<uchar>(i, j) - left_mean, 2); right_std += pow(right.at<uchar>(i, j) - right_mean, 2); numerator += (tmp.at<uchar>(i, j) - left_mean) * (right.at<uchar>(i, j) - right_mean); } } numerator /= n; left_std /= n; right_std /= n; res_ncc = numerator / (sqrt(left_std) * sqrt(right_std)) / n; } // greater NCC value found if (res_ncc > max_ncc[y][x]) { max_ncc[y][x] = res_ncc; // for better visualization depth.at<uchar>(y, x) = (uchar)(offset * 3); } } } } return depth; }
- ASW
- 公式:
- 实现:
Mat asw(Mat in1, Mat in2, string type, int max_disparity = 79,int kernel_size = 3) { int width = in1.size().width; int height = in1.size().height; double k = 3, gamma_c = 30, gamma_g = 20; // ASW parameters Mat depth(height, width, 0); vector< vector<double> > min_asw; // store min ASW value Mat left,right; cvtColor(in1,left,CV_BGR2GRAY); cvtColor(in2,right,CV_BGR2GRAY); for (int i = 0; i < height; ++i) { vector<double> tmp(width, numeric_limits<double>::max()); min_asw.push_back(tmp); } for (int offset = 1; offset <= max_disparity; offset++) { Mat tmp(height, width, 0); // shift image depend on type to save calculation time if (type == "left") { for (int y = 0; y < height; y++) { for (int x = 0; x < offset; x++) { tmp.at<uchar>(y, x) = right.at<uchar>(y, x); } for (int x = offset; x < width; x++) { tmp.at<uchar>(y, x) = right.at<uchar>(y, x - offset); } } } else if (type == "right") { for (int y = 0; y < height; y++) { for (int x = 0; x < width - offset; x++) { tmp.at<uchar>(y, x) = left.at<uchar>(y, x + offset); } for (int x = width - offset; x < width; x++) { tmp.at<uchar>(y, x) = left.at<uchar>(y, x); } } } else { Mat tmp(0, 0, 0); return tmp; } // calculate each pixel's ASW value for (int y = 0; y < height; y++) { for (int x = 0; x < width; x++) { int start_x = max(0, x - kernel_size); int start_y = max(0, y - kernel_size); int end_x = min(width - 1, x + kernel_size); int end_y = min(height - 1, y + kernel_size); double E = 0; if (type == "left") { double numerator = 0; double denominator = 0; for (int i = start_y; i <= end_y; i++) { for (int j = start_x; j <= end_x; j++) { double delta_c1 = fabs(left.at<uchar>(i, j) - left.at<uchar>(y, x)); double delta_c2 = fabs(tmp.at<uchar>(i, j) - tmp.at<uchar>(y, x)); double delta_g = sqrt((i - y) * (i - y) + (j - x) * (j - x)); double w1 = k * exp(-(delta_c1 / gamma_c + delta_g / gamma_g)); double w2 = k * exp(-(delta_c2 / gamma_c + delta_g / gamma_g)); numerator += w1 * w2 * fabs(left.at<uchar>(i, j) - tmp.at<uchar>(i, j)); denominator += w1 * w2; } } E = numerator / denominator; } else { double numerator = 0; double denominator = 0; for (int i = start_y; i <= end_y; i++) { for (int j = start_x; j <= end_x; j++) { double delta_c1 = fabs(right.at<uchar>(i, j) - right.at<uchar>(y, x)); double delta_c2 = fabs(tmp.at<uchar>(i, j) - tmp.at<uchar>(y, x)); double delta_g = sqrt((i - y) * (i - y) + (j - x) * (j - x)); double w1 = k * exp(-(delta_c1 / gamma_c + delta_g / gamma_g)); double w2 = k * exp(-(delta_c2 / gamma_c + delta_g / gamma_g)); numerator += w1 * w2 * fabs(right.at<uchar>(i, j) - tmp.at<uchar>(i, j)); denominator += w1 * w2; } } E = numerator / denominator; } // smaller ASW found if (E < min_asw[y][x]) { min_asw[y][x] = E; // for better visualization depth.at<uchar>(y, x) = (uchar)(offset * 3); } } } } return depth; }
结果: