抱歉，您的浏览器无法访问本站
本页面需要浏览器支持（启用）JavaScript
了解详情 >

CodingStudio

努力进步

isp_pipeline练习实录

Interview

发布于：2023年4月7日

字数：7.4k 字

时长：38 分钟

次浏览

0 前言

  • 记录自己使用C++和Python练习ISP Pipeline开发
  • 设备: iphone 11
    • 软件: ProCam 8
  • 电脑环境: VS 2019, miniconda, opencv4

1 读取raw文件

1.1 Adobe DNG SDK

  • 使用Adobe DNG SDK处理DNG RAW数据
  • 准备工作
    1. 下载Adobe DNG SDK 下载直链： 跳转至下载页面
    2. 默认已经安装好VS2019，另外安装cmake(跳转至下载页面)
      • 注意在安装时将cmake添加进系统变量
    3. 将下载的Adobe DNG SDK解压缩到自己的dng_path
    4. 下载项目构建需要使用的libjpeg、zlib、expat
      1. 下载libjpeg(跳转至下载页面)，随后将压缩文件里的文件解压至dng_path\libjpeg， 将该文件夹内的jconfig.h删除, 并将jconfig.vc重命名为jconfig.h
      2. 下载zlib(跳转至下载页面)，随后将**压缩文件内根目录下的.c和.h**文件解压到dng_path\xmp\toolkit\third-party\zlib目录内
      3. 下载expat(跳转至下载页面, github链接)，随后将压缩文件内的lib文件夹解压到dng_path\xmp\toolkit\XMPCore\third-party\expat\public文件夹内
    5. 构建xmp
      1. 打开dng_path\xmp\toolkit\XMPCore\build\CMake64Static_VC16\XMPCore64.sln
      2. 选择Debug或Release模式下ARM64或x64进行生成
      3. 打开dng_path\xmp\toolkit\XMPFiles\build\CMake64Static_VC16\XMPFiles64.sln
      4. 选择Debug或Release模式下ARM64或x64进行生成
    6. 打开dng_path\dng_sdk\projects\win\dng_validate.sln, 选择合适的平台(ARM64或X64)即可生成
  • 为了方便调试,可将链接器的输出目录(OutputFile)和项目的输出路径(TargetPath)设置一致, 打开项目解决方案的属性，找到链接器->常规->输出文件, 将属性值设置为$(TargetPath)

1.2 libraw

  • 使用LibRaw处理DNG文件
  • 准备工作
    1. 下载LibRaw文件 跳转至下载页面
    2. 解压至LibRaw文件夹
    3. 在该文件夹下打开命令行，输入nmake -f Makefile.msvc进行编译
      • 注意已经安装好VS2019
    4. 编译得到的动态链接库文件在bin/libraw.dll，而链接库在lib/libraw.lib中，头文件在libraw文件夹中
      • 静态链接库为lib/libraw_static.lib
  • 推荐使用动态链接库进行项目开发

1.3 C++项目构建

  • 本教程使用LibRaw构建ISP-pipeline处理流程
    • 需配置好opencv4(C++)环境
  • 项目属性配置
    • 配置属性
      • 常规：输出目录$(SolutionDir)build
      • 常规：中间目录$(SolutionDir)build\dist\$(Configuration)\
    • C/C++
      • 常规：附加包含目录：C:\Program Files\opencv\build\include\opencv2,C:\Program Files\opencv\build\include,$(SolutionDir)include
    • 链接器
      • 常规：附加库目录：C:\Program Files\opencv\build\x64\vc15\lib,$(SolutionDir)lib
      • 输入：附加依赖项：libraw.lib,opencv_world460d.lib
    • 将libraw.dll文件放在输出目录中
  • 项目结构
1
2
3
4
5
6
7
8
9
10
11
12
├─ isp_pipeline
│  ├─build
│  │  ├─dist
│  │  ├─libraw.dll
│  ├─include
│  │  ├─libraw
│  │  ├─Raw.cpp
│  │  ├─Raw.h
│  ├─lib
│  │  ├─libraw.lib
│  ├─res
│  ├─demo.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
// main.cpp
#include <iostream>
#include <vector>
#include <ctime>

#include "libraw.h"
#include "module.h"

int main() {
    cv::utils::logging::setLogLevel(cv::utils::logging::LOG_LEVEL_ERROR);   // 只输出错误日志
    clock_t begin = clock();

    const char* file = "IMG_1.dng";
    // open dng file
    LibRaw* iProcessor = new LibRaw;    
    iProcessor->open_file(file);
    iProcessor->unpack();

    // raw2Mat and cut2roi
    cv::Mat raw = cv::Mat(iProcessor->imgdata.sizes.raw_height, iProcessor->imgdata.sizes.raw_width, CV_16UC1, iProcessor->imgdata.rawdata.raw_image);
    raw = raw(cv::Rect(iProcessor->imgdata.sizes.left_margin, iProcessor->imgdata.sizes.top_margin, iProcessor->imgdata.sizes.width, iProcessor->imgdata.sizes.height));

    clock_t end = clock();    
    std::cout << "-------------------------------------------" << std::endl;
    std::cout << "ISP Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC << " s." << std::endl;
    // recycle LibRaw and delete ptr
    iProcessor->recycle();
    delete iProcessor;
}
  • 注意事项：
    • 上述程序中获取raw数据在iProcessor->imgdata.rawdata.raw_image字段中，且在构建cv::Mat矩阵时，需要注意范围为raw_height到raw_width
      • 因为在CMOS中存在OB区，当获取的尺寸不对时，获取不到正确的数据。同时可以看出此时的图像存在黑边
      • 故需要使用cv::Rect进行截取，截取的参数为left_margin,top_margin,width,height
    • 在保存时，opencv会自动将CV_16UC1转换为CV_8UC1进行保存，而三通道的16位的数据(CV_16UC3)需要自己手动进行转换为(CV_8UC3)才可以正确保存

2 DPC

  • 原因:
    • CMOS器件自身工艺瑕疵造成
    • 光线采集器件存在缺陷
  • 解决方法
    • 简单方法: 在5×5邻域内同一颜色通道相对中心像素有8个邻近像素, 计算中心像素与临近像素的差值, 若差值有正有负则正常; 否则利用阈值进行判断, 差值的绝对值都超过阈值则为坏点则用中位数替换
    • 梯度方法:
      • 计算四个方向的梯度: δ=(xi,0+xi,2−2∗xi,1)\delta = (x_{i,0}+x_{i,2}-2*x{i,1}), 得到12个梯度值
      • 求各个方向梯度的绝对值的中值, 利用四个中值的最小值作为边缘方向
      • 若边缘为水平(竖直方向)
        • 中心点的梯度大于同方向梯度绝对值和的4倍，则中心点为坏点
        • 若∣P4−Pc∣<∣Pc−P5∣|P_4-P_c| < |P_c -P_5|, 则中心像素值更靠近P4P_4, 则Pc=P4+(P2+P7−P1−P6)/2P_c = P_4+(P_2+P_7-P_1-P_6)/2
      • 若边缘为45度(135度方向)
        • 计算135度方向梯度两两之差的绝对值的和, 若和大于100且45度中心梯度大于两边梯度的3倍且135度中心梯度大于两边梯度的3倍, 则为坏点; 若和小于100且45度中心梯度大于两边梯度的3倍, 则为坏点
        • 若∣P3−Pc∣<∣Pc−P6∣|P_3-P_c| < |P_c -P_6|, 则中心像素值更靠近P3P_3, 则Pc=P3+(P4+P7−P2−P5)/2P_c = P_3+(P_4+P_7-P_2-P_5)/2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
// 使用梯度法进行校正
int median(int a, int b, int c) {
	return a >= b ? (b >= c ? b : (a >= c ? c : a)) : (a >= c ? a : (b >= c ? c : b));
}

cv::Mat padding(cv::Mat& img, uchar* pads)
{
	cv::Mat img_pad = cv::Mat(cv::Size(img.cols + pads[2] + pads[3], img.rows + pads[0] + pads[1]), img.type());
	cv::copyMakeBorder(img, img_pad, pads[0], pads[1], pads[2], pads[3], cv::BORDER_REFLECT);
	return img_pad;
}

void dpc(cv::Mat& src, uchar thres)
{
	clock_t begin = clock();
	uchar pads[4] = { 2, 2, 2, 2 };
	cv::Mat src_p = padding(src, pads);
	ushort* p0, * p1, * p2, * p3, * p4, * p5, * p6, * p7, * p8, * p;
	int grad_h1 = 0, grad_h2 = 0, grad_h3 = 0, grad_v1 = 0, grad_v2 = 0, grad_v3 = 0;
	int grad_45_1 = 0, grad_45_2 = 0, grad_45_3 = 0, grad_135_1 = 0, grad_135_2 = 0, grad_135_3 = 0;
	int grad_h = 0, grad_v = 0, grad_45 = 0, grad_135 = 0;
	std::vector<int> gradient(4, 0);
	int grad_sum = 0;
	for (int i = 0; i < src_p.rows - 4; ++i) {
		std::cout << "\r" << "DPC: ";
		std::cout << std::setw(6) << std::fixed << std::setprecision(2) << (float)i / (src_p.rows - 2) * 100 << "%";
		p0 = src_p.ptr<ushort>(i + 2);
		p1 = src_p.ptr<ushort>(i);
		p2 = src_p.ptr<ushort>(i);
		p3 = src_p.ptr<ushort>(i);
		p4 = src_p.ptr<ushort>(i + 2);
		p5 = src_p.ptr<ushort>(i + 2);
		p6 = src_p.ptr<ushort>(i + 4);
		p7 = src_p.ptr<ushort>(i + 4);
		p8 = src_p.ptr<ushort>(i + 4);
		p = src.ptr<ushort>(i);
		for (int j = 0; j < src_p.cols - 4; j++) {
			grad_h1 = abs(p1[j] + p3[j + 4] - 2 * p2[j + 2]);
			grad_h2 = abs(p4[j] + p5[j + 4] - 2 * p0[j + 2]);
			grad_h3 = abs(p6[j] + p8[j + 4] - 2 * p7[j + 2]);
			grad_v1 = abs(p1[j] + p4[j] - 2 * p6[j]);
			grad_v2 = abs(p2[j + 2] + p7[j + 2] - 2 * p0[j + 2]);
			grad_v3 = abs(p3[j + 4] + p8[j + 4] - 2 * p5[j + 4]);
			grad_45_1 = 2 * abs(p2[j + 2] - p4[j]);
			grad_45_2 = abs(p3[j + 4] + p6[j] - 2 * p0[j + 2]);
			grad_45_3 = 2 * abs(p5[j + 4] - p7[j + 2]);
			grad_135_1 = 2 * abs(p2[j + 2] -p5[j + 4]);
			grad_135_2 = abs(p1[j] + p8[j + 4] - 2 * p0[j + 2]);
			grad_135_3 = 2 * abs(p4[j] - p7[j + 2]);

			grad_h = median(grad_h1, grad_h2, grad_h3);
			grad_v = median(grad_v1, grad_v2, grad_v3);
			grad_45 = median(grad_45_1, grad_45_2, grad_45_3);
			grad_135 = median(grad_135_1, grad_135_2, grad_135_3);
			gradient = { grad_h, grad_v, grad_45, grad_135 };
			auto minPosition = std::min_element(gradient.begin(), gradient.end());
			if (minPosition == gradient.begin() && grad_h2 > 4*(grad_h1+grad_h3))
			{
				if (abs(p4[j] - p0[j + 2]) < abs(p0[j + 2] - p5[j + 4])) {
					p[j] = p4[j] + (p2[j + 2] + p7[j + 2] - p1[j] - p6[j]) / 2;
				}
				else {
					p[j] = p5[j + 4] + (p2[j + 2] + p7[j + 2] - p3[j + 4] - p8[j + 4]) / 2;
				}
			}
			if (minPosition == gradient.begin()+1 && grad_v2 > 4 * (grad_v1 + grad_v3))
			{
				if (abs(p2[j + 2] - p0[j + 2]) < abs(p0[j + 2] - p7[j + 2])) {
					p[j] = p2[j + 2] + (p4[j] + p5[j + 4] - p1[j] - p3[j + 4]) / 2;
				}
				else {
					p[j] = p7[j + 2] + (p4[j] + p5[j + 4] - p6[j] - p8[j + 4]) / 2;
				}
			}
			if (minPosition == gradient.begin()+2)
			{
				grad_sum = abs(grad_135_1 - grad_135_2) + abs(grad_135_1 - grad_135_3) + abs(grad_135_2 - grad_135_3);
				if (grad_sum > 100) {
					if (grad_45_2 > 3 * (grad_45_1 + grad_45_3) && grad_135_2 > 3 * (grad_135_1 + grad_135_3)) {
						if (abs(p3[j + 4] - p0[j + 2]) < abs(p0[j + 2] - p6[j])) {
							p[j] = p3[j + 4] + (p4[j] + p7[j + 2] - p2[j + 2] - p5[j + 4]) / 2;
						}
						else {
							p[j] = p6[j] - (p4[j] + p7[j + 2] - p2[j + 2] - p5[j + 4]) / 2;
						}
					}
				}
				else {
					if (grad_45_2 > 3 * (grad_45_1 + grad_45_3)) {
						if (abs(p3[j + 4] - p0[j + 2]) < abs(p0[j + 2] - p6[j])) {
							p[j] = p3[j + 4] + (p4[j] + p7[j + 2] - p2[j + 2] - p5[j + 4]) / 2;
						}
						else {
							p[j] = p6[j] - (p4[j] + p7[j + 2] - p2[j + 2] - p5[j + 4]) / 2;
						}
					}
				}
			}
			if (minPosition == gradient.begin()+3)
			{
				grad_sum = abs(grad_45_1 - grad_45_2) + abs(grad_45_1 - grad_45_3) + abs(grad_45_2 - grad_45_3);
				if (grad_sum > 100) {
					if (grad_135_2 > 3 * (grad_135_1 + grad_135_3) && grad_45_2 > 3 * (grad_45_1 + grad_45_3)) {
						if (abs(p1[j] - p0[j + 2]) < abs(p0[j + 2] - p8[j + 4])) {
							p[j] = p1[j] + (p5[j + 4] + p6[j] - p2[j + 2] - p4[j]) / 2;
						}
						else {
							p[j] = p8[j + 4] - (p5[j + 4] + p6[j] - p2[j + 2] - p4[j]) / 2;
						}
					}
				}
				else {
					if (grad_135_2 > 3 * (grad_135_1 + grad_135_3)) {
						if (abs(p1[j] - p0[j + 2]) < abs(p0[j + 2] - p8[j + 4])) {
							p[j] = p1[j] + (p5[j + 4] + p6[j] - p2[j + 2] - p4[j]) / 2;
						}
						else {
							p[j] = p8[j + 4] - (p5[j + 4] + p6[j] - p2[j + 2] - p4[j]) / 2;
						}
					}
				}
			}

		}
	}
	clock_t end = clock();
	std::cout << "\r" << "DPC Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}


uchar dpc_thres = 30;
dpc(raw, dpc_thres);                        //DPC，坏点矫正，使用梯度进行解决

3 BLC

  • 原因：
    1. 常用的光电器件在没有光照时仍然有电压输出
      • CMOS内部为PN结, 工作在反向电压下, 没有光照时依然有微小电流为暗电流
    2. sensor到isp中有AD转换，但AD存在灵敏度问题, 通常认为添加一个固定值, 使低于AD阈值时也可以进行转换, 保留更多暗部细节
  • 校正方法
    • sensor端算法: 利用OB区的像素值(取平均, 曲线拟合)进行较准
    • isp端算法:
      1. 扣除暗电流: 利用黑帧RAW图每个通道的平均值进行校正, 并对GrG_r和GbG_b进行归一化(R和B通道不进行归一化, 后续使用AWB进行校正)
      2. ISO联动: 暗电流与Gain值和温度相关, 简历ISO与暗电流之间的关系, 每次查表校正
      3. 曲线拟合: 在黑帧中选取部分位置和校正值存储, 后续其他像素校正值通过线性插值进行精准校正
  • 练习实录
    • 在dng文件中包含该文件的暗电平校正值和白电平值, 包含在imgdata.color.dng_levels.dng_black与imgdata.color.dng_levels.dng_whitelevel[0]字段中
    • 利用上述两个字段, 使用简单的扣除暗电流方法即可校正暗电流
    • 注:在该方案中, 并不对GG通道进行归一化，发现对该通道归一化, 会导致偏色, 原因未知
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
std::vector<cv::Mat> split_bayer(cv::Mat& bayer_array, std::string& bayer_pattern) {
	cv::Mat R = cv::Mat(bayer_array.size() / 2, bayer_array.type());
	cv::Mat Gr = cv::Mat(bayer_array.size() / 2, bayer_array.type());
	cv::Mat Gb = cv::Mat(bayer_array.size() / 2, bayer_array.type());
	cv::Mat B = cv::Mat(bayer_array.size() / 2, bayer_array.type());
	ushort* p_R, * p_Gr, * p_Gb, * p_B, * p_bayer1, * p_bayer2;
	for (int i = 0; i < bayer_array.rows / 2; ++i) {
		p_R = R.ptr<ushort>(i);
		p_Gr = Gr.ptr<ushort>(i);
		p_Gb = Gb.ptr<ushort>(i);
		p_B = B.ptr<ushort>(i);
		p_bayer1 = bayer_array.ptr<ushort>(2 * i);
		p_bayer2 = bayer_array.ptr<ushort>(2 * i + 1);
		for (int j = 0; j < bayer_array.cols / 2; ++j) {
			if (bayer_pattern == "gbrg") {
				p_R[j] = p_bayer1[2 * j + 1];
				p_Gr[j] = p_bayer2[2 * j + 1];
				p_Gb[j] = p_bayer1[2 * j];
				p_B[j] = p_bayer2[2 * j];
			} 
			else if (bayer_pattern == "rggb") {
				p_R[j] = p_bayer1[2 * j];
				p_Gr[j] = p_bayer2[2 * j];
				p_Gb[j] = p_bayer1[2 * j + 1];
				p_B[j] = p_bayer2[2 * j + 1];
			}
			else if (bayer_pattern == "bggr") {
				p_R[j] = p_bayer2[2 * j + 1];
				p_Gr[j] = p_bayer1[2 * j + 1];
				p_Gb[j] = p_bayer2[2 * j];
				p_B[j] = p_bayer1[2 * j];
			}
			else if (bayer_pattern == "grbg") {
				p_R[j] = p_bayer2[2 * j];
				p_Gr[j] = p_bayer1[2 * j];
				p_Gb[j] = p_bayer2[2 * j + 1];
				p_B[j] = p_bayer1[2 * j + 1];
			}
			else {
				throw "bayer pattern is not declared!\n";
			}
		}
	}
	return std::vector<cv::Mat> { R,Gr,Gb,B };
}

cv::Mat reconstruct_bayer(std::vector<cv::Mat>& sub_arrays, std::string& bayer_pattern) {
	cv::Mat bayer_array = cv::Mat(sub_arrays[0].size() * 2, sub_arrays[0].type());
	ushort* p_R, * p_Gr, * p_Gb, * p_B, * p_bayer1, * p_bayer2;
	for (int i = 0; i < sub_arrays[0].rows; ++i) {
		p_R = sub_arrays[0].ptr<ushort>(i);
		p_Gr = sub_arrays[1].ptr<ushort>(i);
		p_Gb = sub_arrays[2].ptr<ushort>(i);
		p_B = sub_arrays[3].ptr<ushort>(i);
		p_bayer1 = bayer_array.ptr<ushort>(2 * i);
		p_bayer2 = bayer_array.ptr<ushort>(2 * i + 1);
		for (int j = 0; j < sub_arrays[0].cols; ++j) {
			if (bayer_pattern == "gbrg") {
				p_bayer1[2 * j + 1] = p_R[j];
				p_bayer2[2 * j + 1] = p_Gr[j];
				p_bayer1[2 * j] = p_Gb[j];
				p_bayer2[2 * j] = p_B[j];
			}
			else if (bayer_pattern == "rggb") {
				p_bayer1[2 * j] = p_R[j];
				p_bayer2[2 * j] = p_Gr[j];
				p_bayer1[2 * j + 1] = p_Gb[j];
				p_bayer2[2 * j + 1] = p_B[j];
			}
			else if (bayer_pattern == "bggr") {
				p_bayer2[2 * j + 1] = p_R[j];
				p_bayer1[2 * j + 1] = p_Gr[j];
				p_bayer2[2 * j] = p_Gb[j];
				p_bayer1[2 * j] = p_B[j];
			}
			else if (bayer_pattern == "grbg") {
				p_bayer2[2 * j] = p_R[j];
				p_bayer1[2 * j] = p_Gr[j];
				p_bayer2[2 * j + 1] = p_Gb[j];
				p_bayer1[2 * j + 1] = p_B[j];
			}
			else {
				throw "bayer pattern is not declared!\n";
			}
		}
	}
	return bayer_array;
}

void blc(cv::Mat& src, std::string bayer_pattern, float black_level, float white_level)
{
	clock_t begin = clock();
	std::vector<cv::Mat> bayer4 = split_bayer(src, bayer_pattern);
	int pixel_R = 0, pixel_Gr = 0, pixel_Gb = 0, pixel_B = 0;
	ushort* p_R, * p_Gr, * p_Gb, * p_B;
	for (int i = 0; i < bayer4[0].rows; ++i) {
		p_R = bayer4[0].ptr<ushort>(i);
		p_Gr = bayer4[1].ptr<ushort>(i);
		p_Gb = bayer4[2].ptr<ushort>(i);
		p_B = bayer4[3].ptr<ushort>(i);
		for (int j = 0; j < bayer4[0].cols; ++j) {
			pixel_R = p_R[j] - black_level;
			pixel_B = p_B[j] - black_level;
			pixel_Gr = p_Gr[j] - black_level;
			pixel_Gb = p_Gb[j] - black_level;

			p_R[j] = ((pixel_R < 0) ? 0 : (pixel_R > white_level ? white_level : pixel_R));
			p_Gr[j] = ((pixel_Gr < 0) ? 0 : (pixel_Gr > white_level ? white_level : pixel_Gr));
			p_Gb[j] = ((pixel_Gb < 0) ? 0 : (pixel_Gb > white_level ? white_level : pixel_Gb));
			p_B[j] = ((pixel_B < 0) ? 0 : (pixel_B > white_level ? white_level : pixel_B));
		}
	}
	src = reconstruct_bayer(bayer4, bayer_pattern);
	clock_t end = clock();
	std::cout << "BLC Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}

std::string BAYER = "rggb";
float black = iProcessor->imgdata.color.dng_levels.dng_black;
float white = iProcessor->imgdata.color.dng_levels.dng_whitelevel[0];

blc(raw, BAYER, black, white);              //BLC

4 AAF(抗混叠滤波)

  • 混叠产生的原因
    • 在后续的Demosaics会进行插值操作, 该操作可能会导致混叠
    • 由于无法对一个函数无限地取样, 因此在数字图像中会出现混叠现象(空间混叠和时间混叠)
    • 空间混叠, 由于欠采样导致, 会引入伪影
      • 对图像放大视为过取样, 混叠不严重; 而对图像缩小视为欠取样, 混叠较为严重
    • 根本原因：高频信息在采样过程中重叠到低频段
  • 解决方法：
    • 使用抗混叠滤波器进行滤波;(本质上使用低通滤波器进行滤波, 衰减高频分量, 使重新采样后的分辨率可以表示高频信息)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
void aaf(cv::Mat& src, cv::Mat& kernel)
{
	clock_t begin = clock();
	cv::Mat origin_img = src.clone();
	cv::filter2D(origin_img, src, CV_16UC1, kernel);
	clock_t end = clock();
	std::cout << "AAF Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms."<< std::endl;
}

cv::Mat aaf_kernel = (cv::Mat_<float>(5, 5) <<
        1, 0, 1, 0, 1,
        0, 0, 0, 0, 0,
        1, 0, 8, 0, 1,
        0, 0, 0, 0, 0,
        1, 0, 1, 0, 1) / 16;

aaf(raw, aaf_kernel);                       //AAF, 抗混叠滤波，抑制频谱混叠

5 WBGC(白平衡校正)

  • 产生原因: 传感器不具有人眼的不同光照色温下的色彩恒常性, 白平衡模块需要将人眼看来白色的物体进行色彩的还原, 使其在照片上呈现白色
  • 校正方法:
    • 灰度世界法: 任何一幅图像具有足够的色彩变化, 则RGB分量的均值会趋于0(一般选取G通道作为参考值计算RB通道的Gain值进行校正)
    • 完美反射法: 图像中最亮的点就是白色或者镜面反射出来的, 最亮的点就是光源的属性, 但该点本身应该是白色, 以此为基础计算Gain值进行校正(计算图像RGB分别的最大值, 用三个最大值的最大值除以每个通道的最大值即可得到每个通道的Gain值)
    • 基于色温的方法
      1. 不同色温环境下灰卡得到Gain值关系曲线与R/G与色温关系曲线
      2. 获取了一张图像的拍摄色温, 通过第二张图获取R/G的值随后可根据图一计算出B/G从而计算R和B的Gain值
  • 练习实录
    • 在dng文件中具有AWB矩阵, 该矩阵存放在imgdata.color.dng_levels.asshotneutral字段中
    • 利用上述字段分别对图像的各个通道进行校正
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
void wbgc(cv::Mat& src, std::string bayer_pattern, std::vector<float> gain, ushort clip)
{
	gain = { 1 / gain[0] * 1024, 1 / gain[1] * 1024, 1 / gain[2] * 1024 };
	clock_t begin = clock();
	std::vector<cv::Mat> rggb = split_bayer(src, bayer_pattern);
	rggb[0] *= (gain[0] / 1024);
	rggb[1] *= (gain[1] / 1024);
	rggb[2] *= (gain[1] / 1024);
	rggb[3] *= (gain[2] / 1024);
	src = reconstruct_bayer(rggb, bayer_pattern);
	ushort* p_img;
	for (int i = 0; i < src.rows; ++i) {
		p_img = src.ptr<ushort>(i);
		for (int j = 0; j < src.cols; ++j) {
			p_img[j] = (p_img[j] < 0) ? 0 : (p_img[j] > clip ? clip : p_img[j]);
		}
	}
	clock_t end = clock();
	std::cout << "AWB Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}

std::vector<float> parameter(iProcessor->imgdata.color.dng_levels.asshotneutral, iProcessor->imgdata.color.dng_levels.asshotneutral + 3);
ushort CLIP = 4095;

wbgc(raw, BAYER, parameter, CLIP);          //WBGC，白平衡较准

6 BNR

  • 在RAW域进行降噪处理效果会更好, 此时的噪声没有经过各种模块的调节而导致变大
  • 常用方法:
    • 中值滤波
    • 双边滤波
    • BM3D降噪
  • 练习实录
    • 在Raw数据分为4个通道, 在每个通道上应用中值滤波进行处理
1
2
3
4
5
6
7
8
9
10
11
12
13
14
void bnr(cv::Mat& src, std::string bayer_pattern, uchar ksize) {
	clock_t begin = clock();
	std::vector<cv::Mat> rggb = split_bayer(src, bayer_pattern);
	medianBlur(rggb[0], rggb[0], ksize);
	medianBlur(rggb[1], rggb[1], ksize);
	medianBlur(rggb[2], rggb[2], ksize);
	medianBlur(rggb[3], rggb[3], ksize);
	src = reconstruct_bayer(rggb, bayer_pattern);
	clock_t end = clock();
	std::cout  << "BNF Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}

uchar bnr_size = 3;
bnr(raw, BAYER, bnr_size);

7 Demosaics

  • 原理：一般的CMOS为前照式CMOS, 使用CFA阵列将光线分解为RGB三个分量用感光器件接收
  • 解决方法:
    1. 简单线性插值: 效果一般, 清晰度降低, 高频产生伪彩, 边缘产生伪像
    2. 色比法: 一个邻域内不同颜色通道的比值固定, 首先插值计算出G的缺失值, 随后根据比值恒定计算出其他的缺失值(清晰度变差, 产生伪彩和伪像)
  • 练习实录
    • 使用opencv自身的算法进行实现
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
void cfa(cv::Mat& src, std::string bayer_pattern)
{
	clock_t begin = clock();
	if (bayer_pattern == "gbrg") {
		cv::cvtColor(src, src, cv::COLOR_BayerGB2RGB);
	}
	else if (bayer_pattern == "rggb") {
		cv::cvtColor(src, src, cv::COLOR_BayerRG2RGB);
	}
	else if (bayer_pattern == "bggr") {
		cv::cvtColor(src, src, cv::COLOR_BayerBG2RGB);
	}
	else if (bayer_pattern == "grbg") {
		cv::cvtColor(src, src, cv::COLOR_BayerGR2RGB);
	}
	else {
		throw "bayer pattern is not declared!\n";
	}
	clock_t end = clock();
	std::cout << "CFA Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}


cfa(raw, BAYER);                            //CFA

8 CCM

  • 原因: sensor和人眼对颜色的感知存在差异, 且不同sensor之间的颜色感知也存在差异; 因此需要将sensor感光数据转换至人眼感光数据
  • 解决方法
    • 3D-LUT查表法
    • 多项式拟合方法: 拟合颜色校正矩阵(CCM矩阵, 该矩阵的每个行相等,即r=g=br=g=b)
    • 使用24色标准色卡(常规做法): 通常在CIE_LAB空间具有ΔE=ΔL2+ΔA2+ΔB2,L为亮度,AB为色相\Delta E = \sqrt{\Delta L^2+ \Delta A^2 + \Delta B^2}, L为亮度, AB为色相参数, 更精确表示人眼的特性
    • 通常是求解最优点, 而不是直接计算出CCM矩阵
  • 经过校正后, sensorRGB贴近sRGB, 且具有了负响应
  • 练习实录
    • 在dng文件中具有CCM矩阵,该字段在imgdata.color.cmatrix字段中
    • 获取该字段后,使用矩阵相乘即可实现CCM校正
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
void ccm(cv::Mat& src, cv::Mat& ccm)
{
	clock_t begin = clock();
	cv::Mat img0 = src.clone();
	cv::Vec3s* p, * p0;
	for (int i = 0; i < src.rows; ++i) {
		p0 = img0.ptr<cv::Vec3s>(i);
		p = src.ptr<cv::Vec3s>(i);
		for (int j = 0; j < src.cols; ++j) {
			for (int k = 0; k < 3; ++k) {
				p[j][k] = ccm.ptr<float>(k)[0] * p0[j][0] + ccm.ptr<float>(k)[1] * p0[j][1] + ccm.ptr<float>(k)[2] * p0[j][2];
			}
		}
	}
	clock_t end = clock();
	std::cout << "CCM Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}

cv::Mat ccmatrix = cv::Mat(3, 4, CV_32F, iProcessor->imgdata.color.cmatrix)(cv::Rect(0, 0, 3, 3));

ccm(raw, ccmatrix);                         //CCM, sensor RGB转sRGB，转换到人眼相适应的颜色空间，颜色较准

9 Gamma校正

  • 产生原因: 人眼视觉特征和显示器的视觉特性不一致导致，为了使人眼看起来更加舒服，对暗区加强以提高动态范围和暗区细节
  • 解决方法
    • 查表法
    • 线性插值法: 在Gamma曲线上取一些采样点进行存储, 校正时插值获取校正值
  • 练习实录
    • 根据Gamma值建立表, 对每个像素进行查表解决
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
void gc(cv::Mat& src, float gamma, ushort clip)
{
	clock_t begin = clock();
	std::vector<int> LUT{};
	for (int i = 0; i < clip + 1; ++i) {
		float lx = std::pow(((float)i / clip), (float)gamma) * (float)sdr_max_value;
		LUT.push_back((int)lx);
	}

	cv::Vec3s* p;
	for (int i = 0; i < src.rows; ++i) {
		p = src.ptr<cv::Vec3s>(i);
		for (int j = 0; j < src.cols; ++j) {
			for (int k = 0; k < 3; ++k) {
				p[j][k] = ((p[j][k] < 0) ? ushort(LUT[0]) : (p[j][k] > clip ? ushort(LUT[clip]) : ushort(LUT[p[j][k]])));
			}
		}
	}
	clock_t end = clock();
	std::cout << "GaC Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}


float gamma = 0.42;

gc(raw, gamma, CLIP);                       //GC

10 CSC

  • 原因: 将颜色从RGB空间转换至其他空间进行处理, 通常转换至YUV空间, 以便进一步对色度和饱和度进行处理
  • 练习实录:
    • 将颜色转换至YUV420颜色空间, 该颜色空间下YUV三个通道相互分离
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
void csc(cv::Mat& src, cv::Mat& y, cv::Mat& u, cv::Mat& v, float YUV420[3][4])
{
	clock_t begin = clock();
	cv::Vec3s* p;
	float* p_YUV420;
	ushort* p_y, * p_u, * p_v;
	for (int i = 0; i < src.rows; ++i) {
		p = src.ptr<cv::Vec3s>(i);
		p_y = y.ptr<ushort>(i);
		p_u = u.ptr<ushort>(i);
		p_v = v.ptr<ushort>(i);
		for (int j = 0; j < src.cols; ++j) {
			p_y[j] = (p[j][0] * YUV420[0][0] + p[j][1] * YUV420[0][1] + p[j][2] * YUV420[0][2] + YUV420[0][3]);
			p_u[j] = (p[j][0] * YUV420[1][0] + p[j][1] * YUV420[1][1] + p[j][2] * YUV420[1][2] + YUV420[1][3]);
			p_v[j] = (p[j][0] * YUV420[2][0] + p[j][1] * YUV420[2][1] + p[j][2] * YUV420[2][2] + YUV420[2][3]);
		}
	}

	clock_t end = clock();
	std::cout << "CSC Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}

float RGB2YUV420[3][4] = {
        {0.257, 0.504, 0.098, 16},
        {-.148, -.291, 0.439, 128 },
        {0.439, -.368, -.071, 128},
    };

cv::Mat y = cv::Mat(cv::Size(raw.cols, raw.rows), CV_16U);
cv::Mat u = cv::Mat(cv::Size(raw.cols, raw.rows), CV_16U);
cv::Mat v = cv::Mat(cv::Size(raw.cols, raw.rows), CV_16U);
csc(raw, y, u, v, RGB2YUV420);              //CSC

11 NLM

  • 对Y通道进行保边降噪
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
void nlm(cv::Mat& src, uchar ds, uchar Ds, uchar h, uchar clip)
{
	clock_t begin = clock();
	ushort center_y = 0, center_x = 0, start_y = 0, start_x = 0, pixel_pad = 0;
	double sweight = .0, average = .0, wmax = .0, w = .0, dist = .0, pixel = 0;
	ushort* p_src, *p_p, *p_p0, * p_neighbor, * p_center_w;

	uchar pads[4] = { Ds, Ds, Ds, Ds };
	cv::Mat src_p = padding(src, pads);
	for (int y = 0; y < src_p.rows - 2 * Ds; ++y) {
		std::cout << "\r" << "NLM: ";
		std::cout << std::setw(6) << std::fixed << std::setprecision(2) << (float)y / (src_p.rows - 2) * 100 << "%";
		center_y = y + Ds;
		p_src = src.ptr<ushort>(y);
		p_p0 = src_p.ptr<ushort>(center_y);
		for (int x = 0; x < src_p.cols - 2 * Ds; ++x) {
			center_x = x + Ds;
			//calWeights
			sweight = .0, average = .0, wmax = .0, dist = .0;
			for (int m = 0; m < 2 * Ds + 1 - 2 * ds - 1; ++m) {
				start_y = center_y - Ds + ds + m;
				p_p = src_p.ptr<ushort>(start_y);
				for (int n = 0; n < 2 * Ds + 1 - 2 * ds - 1; ++n) {
					start_x = center_x - Ds + ds + n;
					pixel_pad = p_p[start_x];
					if (m != center_y || n != center_x) {
						for (int i = 0; i < 2 * ds + 1; ++i) {
							p_neighbor = src_p.ptr<ushort>(start_y - ds + i);
							p_center_w = src_p.ptr<ushort>(center_y - ds + i);
							for (int j = 0; j < 2 * ds + 1; ++j) {
								dist += (((p_neighbor[start_x - ds + j] - p_center_w[center_x - ds + j]) ^ 2) / (2 * ds + 1) ^ 2);
							}
						}
						w = exp(-dist / (h ^ 2));
						if (w > wmax) wmax = w;
						sweight += w;
						average += (w * pixel_pad);
					}
				}
			}
			average += (wmax * p_p0[center_x]);
			sweight += wmax;
			pixel = average / sweight;
			pixel = (pixel <= 0) ? 0 : (pixel >= clip ? clip : pixel);
			p_src[x] = (ushort)pixel;
		}
	}
	clock_t end = clock();
	std::cout << "\r"  << "NLM Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}

uchar nlm_dw = 1;
uchar nlm_Dw = 3;
uchar nlm_thres = 10;
float bnf_dw[5][5] = {
	{0.5, 0.75,  2, 0.75, 0.5 },
	{0.75,   4,  8,    4, 0.75},
	{2,      4, 64,     4,    2},
	{0.75,   4,  8,    4, 0.75},
	{0.5, 0.75,  2, 0.75, 0.5 },
};
nlm(y, nlm_dw, nlm_Dw, nlm_thres, CLIP);    //NLM

12 BNF

  • 双边滤波：非线性滤波，结合图像的空间邻近度和像素相似度的折中处理，同时考虑空间与信息和灰度相似度，达到保边去噪的目的；
    • 方法：使用两个滤波器（一个函数由几何空间距离决定滤波器系数，另一个由像素差值决定滤波器系数）
      • g(i,j)=∑k,lf(k,l)w(i,j,k,l)∑k,lw(i,j,k,l)g(i, j)=\frac{\sum_{k, l} f(k, l) w(i, j, k, l)}{\sum_{k, l} w(i, j, k, l)}
      • w(i,j,k,l)=exp⁡(−(i−k)2+(j−l)22σd2−∥f(i,j)−f(k,l)∥22σr2)w(i, j, k, l)=\exp \left(-\frac{(i-k)^2+(j-l)^2}{2 \sigma_d^2}-\frac{\|f(i, j)-f(k, l)\|^2}{2 \sigma_r^2}\right).
      • 简单来说双边滤波模板组合由两个模板生成：第一个为高斯模板，第二个为以灰度级的插值作为函数系数生成的模板，两模板点乘即可得到最终滤波模板
    • 可以很好的保留图像边缘细节而滤除掉低频分量的噪声，但效率较低，花费时间较长
  • 练习实录
    • 使用双边滤波进一步保边降噪
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

void bnf(cv::Mat& src, uchar sigmoid_s, uchar sigmoid_p, ushort clip)
{
	clock_t begin = clock();

	uchar pads[4] = { 2, 2, 2, 2 };
	cv::Mat src_p = padding(src, pads);
	int sum_weight = 0, sum_imgA = 0;
	ushort pixel_center = 0;
	double dive = .0, w = .0;
	ushort* p, * p_p;

	for (int y = 0; y < src_p.rows - 4; ++y) {
		std::cout << "\r" << "BNF: ";
		std::cout << std::setw(6) << std::fixed << std::setprecision(2) << (float)y / (src_p.rows - 4) * 100 << "%";
		p = src.ptr<ushort>(y);
		for (int x = 0; x < src_p.cols - 4; ++x) {
			pixel_center = src_p.ptr<ushort>(y + 2)[x + 2];
			sum_imgA = 0;
			sum_weight = 0;
			dive = .0;
			for (int i = 0; i < 5; ++i) {
				p_p = src_p.ptr<ushort>(y + i);
				for (int j = 0; j < 5; ++j) {
					w = exp(-((abs(p_p[x + j] - pixel_center)) ^ 2) / 2 / sigmoid_p ^ 2) * exp(-((2 - i) ^ 2 + (2 - j) ^ 2) / 2 / sigmoid_s);
					sum_imgA += p_p[x + j] * w;
					sum_weight += w;
				}
			}
			dive = sum_imgA / sum_weight;
			dive = (dive < 0) ? 0 : (dive > clip ? clip : dive);
			p[x] = (ushort)dive;
		}
	}
	clock_t end = clock();
	std::cout << "\r" << "BNF Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}

uchar bnf_s_p = 5;
uchar bnf_s_s = 5;
bnf(y, bnf_s_s, bnf_s_p, CLIP);   //BNF

13 EE

  • 练习实录
    • 使用边缘增强增强边缘
    • 通过使用拉普拉斯算子提取出边缘信息, 使用Gain值将边缘增加到原图上
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
void ee(cv::Mat& src, cv::Mat& edgemap, char edge_filter[3][5], ushort clip)
{
	 clock_t begin = clock();

	uchar pads[4] = { 1,1,2,2 };
	cv::Mat src_p = padding(src, pads);
	double em_img, ee_img, tmp_em_img;
	ushort* p_src, * p_edgemap;
	for (int y = 0; y < src_p.rows - 2; ++y) {
		std::cout << "\r" << "EEH: ";
		std::cout << std::setw(6) << std::fixed << std::setprecision(2) << (float)y / (src_p.rows - 4) * 100 << "%";
		p_edgemap = edgemap.ptr<ushort>(y);
		p_src = src.ptr<ushort>(y);
		for (int x = 0; x < src_p.cols - 4; ++x) {
			em_img = 0.0;
			for (int i = 0; i < 3; ++i) {
				for (int j = 0; j < 5; ++j) {
					em_img += src_p.ptr<ushort>(y + i)[x + j] * edge_filter[i][j];
				}
			}
			em_img = em_img / 8;
			ee_img = src_p.ptr<ushort>(y + 1)[x + 2] + em_img;
			ee_img = (ee_img < 0) ? 0 : (ee_img > clip ? clip : ee_img);
			p_edgemap[x] = (ushort)em_img;
			p_src[x] = (ushort)ee_img;
		}
	}
	clock_t end = clock();
	std::cout << "\r" << "EEH Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}

char edge_filter[3][5] = {
	{-1, 0, -1, 0, -1},
	{-1, 0, 8,  0, -1},
	{-1, 0, -1, 0, -1} };

cv::Mat edge_map = cv::Mat(cv::Size(y.cols, y.rows), CV_16U);
ee(y, edge_map, edge_filter, CLIP);           //EE

14 BCC

  • 练习实录
    • 调节y通道的亮度和对比度
  • 算法流程:
    • 对每个通道的值, 分别加上亮度增益
    • 对比度调节时, 将像素值减一个中间值(255/2=127255/2=127), 随后使用原像素值加乘以对比度的差值pixel+(pixel−127)∗contrastpixel+(pixel-127)*contrast
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
void bcc(cv::Mat& src, uchar brightness, uchar contrast, ushort bcc_clip)
{
	clock_t begin = clock();
	contrast = contrast / (2 ^ 5);
	ushort* p;
	int pixel;
	for (int y = 0; y < src.rows; ++y) {
		p = src.ptr<ushort>(y);
		for (int x = 0; x < src.cols; ++x) {
			pixel = p[x];
			pixel = pixel + brightness + (pixel - 127) * contrast;
			pixel = (pixel < 0) ? 0 : (pixel > bcc_clip ? bcc_clip : pixel);
			p[x] = (ushort)pixel;
		}
	}
	clock_t end = clock();
	std::cout << "BCC Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}

uchar brightness = 10;
uchar contrast = 10;
bcc(y, brightness, contrast, CLIP);                                         //BCC

15 转换并存储

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void yuv2rgb(cv::Mat& src, cv::Mat& y, cv::Mat& u, cv::Mat& v)
{
	clock_t begin = clock();
	ushort* p_y, * p_u, * p_v;
	cv::Vec3s* p;
	for (int i = 0; i < src.rows; ++i){
		p = src.ptr<cv::Vec3s>(i);
		p_y = y.ptr<ushort>(i);
		p_u = u.ptr<ushort>(i);
		p_v = v.ptr<ushort>(i);
		for (int j = 0; j < src.cols; ++j){
			p[j][0] = cv::saturate_cast<ushort>(p_y[j] + 1.402 * (p_v[j] - 128));									// R
			p[j][1] = cv::saturate_cast<ushort>(p_y[j] - 0.34413 * (p_u[j] - 128) - 0.71414 * (p_v[j] - 128));		// G
			p[j][2] = cv::saturate_cast<ushort>(p_y[j] + 1.772 * (p_u[j] - 128));									// B
		}
	}
	clock_t end = clock();
	std::cout << "Y2R Done! Elapsed " << double(end - begin) / CLOCKS_PER_SEC * 1000 << " ms." << std::endl;
}

yuv2rgb(raw, y, u, v);
raw.convertTo(raw, CV_8UC3);
cv::imwrite("result.png", raw);

更新于：2023年10月16日

Interview

评论