1. 简介
微软在2004年发表的这篇High-Quality Linear Interpolation for Demosaicing of Bayer-Patterned Color Imagesdemosaic文章,也是matlab自带函数demosaic()使用的方法。该文描述方法在线性解马赛克方法中性能优异,能与一些非线性方法媲美,却有更少的资源消耗,也利于硬件实现。下面的内容中,首先介绍了bilinear demosaic方法,然后介绍Gradient Based Demosaic方法,最后对两者的结果进行了比较。
2.Bilinear Demosaic
目前市场上几乎所有的相机都是单图像传感器的相机,而传感器上的每个像元只能感应R/G/B中的一种颜色,传感器的能感应的颜色大多按bayer pattern拍了(如图1)。为了获得每个像元的R/G/B分量,需要通过demosaic方法插值出该像元不能感应的颜色。
如图1中R像元的G分量,就需要插值获得。而双线性插值,就是通过临近像元取均值获得对应的G分量,公式如下:
其中, (m,n)=(0,−1),(0,1),(−1,0),(1,0) ( m , n ) = ( 0 , − 1 ) , ( 0 , 1 ) , ( − 1 , 0 ) , ( 1 , 0 ) 。同理,每个B像元的G分量,也用相同的方法获得。对于R/B分量,也采用类似的方法进行计算,只是其插值元素为对角线上的4个像元。或者在G像元,R/B分量采用临近的2个像元进行插值。Bilinear Demosaic的Matlab代码为:
function [ demosaiced_img ] = demosaic_bilinear( cfa_image, cfa_pattern )
%DEMOSAIC_BILINEAR Summary of this function goes here
% Detailed explanation goes here
demosaiced_img = zeros(size(cfa_image,1), size(cfa_image, 2), 3);
demosaiced_img(:,:,2) = calc_g_channel(cfa_image,cfa_pattern);
[demosaiced_img(:,:,1),demosaiced_img(:,:,3)] = ...
calc_rb_channel(cfa_image, cfa_pattern);
end
function [ g_channel ] = calc_g_channel(cfa_image, cfa_pattern)
kernel = [0 1/4 0; 1/4 0 1/4; 0 1/4 0];
padded_cfa_image = zeros(size(cfa_image,1)+2,size(cfa_image,2)+2);
padded_cfa_image(2:end-1,2:end-1) = cfa_image;
padded_cfa_image(2:end-1,1) = padded_cfa_image(2:end-1,3);
padded_cfa_image(2:end-1,end) = padded_cfa_image(2:end-1,end-2);
padded_cfa_image(1,:) = padded_cfa_image(3,:);
padded_cfa_image(end,:) = padded_cfa_image(end-2,:);
filtered_cfa = filter2(kernel, padded_cfa_image, 'valid');
g_channel = cfa_image;
if cfa_pattern(1) == 2 && cfa_pattern(4) == 2
g_channel(1:2:end,2:2:end) = filtered_cfa(1:2:end,2:2:end);
g_channel(2:2:end,1:2:end) = filtered_cfa(2:2:end,1:2:end);
elseif cfa_pattern(2) == 2 && cfa_pattern(3) == 2
g_channel(1:2:end,1:2:end) = filtered_cfa(1:2:end,1:2:end);
g_channel(2:2:end,2:2:end) = filtered_cfa(2:2:end,2:2:end);
end
end
function [ r_channel, b_channel ] = calc_rb_channel(cfa_image, cfa_pattern)
kernel_v = [0 1/2 0; 0 0 0; 0 1/2 0];
kernel_h = [0 0 0; 1/2 0 1/2; 0 0 0];
kernel_d = [1/4 0 1/4; 0 0 0; 1/4 0 1/4];
padded_cfa_image = zeros(size(cfa_image,1)+2,size(cfa_image,2)+2);
padded_cfa_image(2:end-1,2:end-1) = cfa_image;
padded_cfa_image(2:end-1,1) = padded_cfa_image(2:end-1,3);
padded_cfa_image(2:end-1,end) = padded_cfa_image(2:end-1,end-2);
padded_cfa_image(1,:) = padded_cfa_image(3,:);
padded_cfa_image(end,:) = padded_cfa_image(end-2,:);
filtered_v = filter2(kernel_v, padded_cfa_image, 'valid');
filtered_h = filter2(kernel_h, padded_cfa_image, 'valid');
filtered_d = filter2(kernel_d, padded_cfa_image, 'valid');
r_channel = cfa_image;
b_channel = r_channel;
if cfa_pattern(1) == 2 && cfa_pattern(4) == 2
if cfa_pattern(2) == 3 && cfa_pattern(3) == 1
b_channel(1:2:end,1:2:end) = filtered_h(1:2:end,1:2:end);
b_channel(2:2:end,2:2:end) = filtered_v(2:2:end,2:2:end);
b_channel(2:2:end,1:2:end) = filtered_d(2:2:end,1:2:end);
r_channel(1:2:end,1:2:end) = filtered_v(1:2:end,1:2:end);
r_channel(1:2:end,2:2:end) = filtered_d(1:2:end,2:2:end);
r_channel(2:2:end,2:2:end) = filtered_h(2:2:end,2:2:end);
else
r_channel(1:2:end,1:2:end) = filtered_h(1:2:end,1:2:end);
r_channel(2:2:end,2:2:end) = filtered_v(2:2:end,2:2:end);
r_channel(2:2:end,1:2:end) = filtered_d(2:2:end,1:2:end);
b_channel(1:2:end,1:2:end) = filtered_v(1:2:end,1:2:end);
b_channel(1:2:end,2:2:end) = filtered_d(1:2:end,2:2:end);
b_channel(2:2:end,2:2:end) = filtered_h(2:2:end,2:2:end);
end
elseif cfa_pattern(2) == 2 && cfa_pattern(3) == 2
if cfa_pattern(1) == 1 && cfa_pattern(4) == 3
r_channel(1:2:end,2:2:end) = filtered_h(1:2:end,2:2:end);
r_channel(2:2:end,1:2:end) = filtered_v(2:2:end,1:2:end);
r_channel(2:2:end,2:2:end) = filtered_d(2:2:end,2:2:end);
b_channel(1:2:end,1:2:end) = filtered_d(1:2:end,1:2:end);
b_channel(1:2:end,2:2:end) = filtered_v(1:2:end,2:2:end);
b_channel(2:2:end,1:2:end) = filtered_h(2:2:end,1:2:end);
else
b_channel(1:2:end,2:2:end) = filtered_h(1:2:end,2:2:end);
b_channel(2:2:end,1:2:end) = filtered_v(2:2:end,1:2:end);
b_channel(2:2:end,2:2:end) = filtered_d(2:2:end,2:2:end);
r_channel(1:2:end,1:2:end) = filtered_d(1:2:end,1:2:end);
r_channel(1:2:end,2:2:end) = filtered_v(1:2:end,2:2:end);
r_channel(2:2:end,1:2:end) = filtered_h(2:2:end,1:2:end);
end
end
end
3. Gradient-based Bilinear Demosaic
相对于纯粹的双线性插值,基于梯度的双线性插值方法不仅采用单通道平面内插值方式,而且考虑待插值像素本身通道的梯度值,将该梯度值以一定比例加到双线性插值的结果上。例如,为了计算R像素处的G通道分量,采用如下公式:
其中, g^(i,j) g ^ ( i , j ) 是待估算的G通道值, g^B(i,j) g ^ B ( i , j ) 是采用双线性插值获得的G通道值, α α 是控制梯度增强的定系数, ΔR(i,j) Δ R ( i , j ) 是R通道在该处的梯度,计算方式为:
其中, (m,n)={(0,−1),(0,2),(−2,0),(2,0)} ( m , n ) = { ( 0 , − 1 ) , ( 0 , 2 ) , ( − 2 , 0 ) , ( 2 , 0 ) } 。同理,用 ΔB(i,j) Δ B ( i , j ) 代替 ΔR(i,j) Δ R ( i , j ) ,即可获得G通道在B像素处的值。
而在计算G通道处的R分量是,采用如下公式:
ΔG(i,j) Δ G ( i , j ) 是由9个点作为支持区域,获得的梯度。对于B像素值所在的位置,R分量的计算方式为:
ΔB(i,j) Δ B ( i , j ) 是由5个点作为支持区域,获得的梯度。根据对称性,B通道的计算方式与R通道类似。
原文中指出, α,β,γ α , β , γ 分别取值 12,58,34 1 2 , 5 8 , 3 4 时,效果最好。根据以上讨论,可以将本文提出的方法转为滤波方式来完成插值,具体分为6个kernel,然后每个kernel乘以 18 1 8 ,滤波完成即可获得每个像素的RGB分量。这6个kernel如下:
R位置的G分量kernel
0 | 0 | -1 | 0 | 0 |
---|---|---|---|---|
0 | 0 | 2 | 0 | 0 |
-1 | 2 | 4 | 2 | -1 |
0 | 0 | 2 | 0 | 0 |
0 | 0 | -1 | 0 | 0 |
B位置的G分量kernel
0 | 0 | -1 | 0 | 0 |
---|---|---|---|---|
0 | 0 | 2 | 0 | 0 |
-1 | 2 | 4 | 2 | -1 |
0 | 0 | 2 | 0 | 0 |
0 | 0 | -1 | 0 | 0 |
G在R行B列的位置,R分量的kernel
0 | 0 | +1/2 | 0 | 0 |
---|---|---|---|---|
0 | -1 | 0 | -1 | 0 |
-1 | 4 | 5 | 4 | -1 |
0 | -1 | 0 | -1 | 0 |
0 | 0 | +1/2 | 0 | 0 |
G在B行R列位置,R分量的kernel
0 | 0 | -1 | 0 | 0 |
---|---|---|---|---|
0 | -1 | 4 | -1 | 0 |
+1/2 | 0 | 5 | 0 | +1/2 |
0 | -1 | 4 | -1 | 0 |
0 | 0 | -1 | 0 | 0 |
B在B行B列,R分量的kernel
0 | 0 | -3/2 | 0 | 0 |
---|---|---|---|---|
0 | 2 | 0 | 2 | 0 |
-3/2 | 0 | 6 | 0 | -3/2 |
0 | 2 | 0 | 2 | 0 |
0 | 0 | -3/2 | 0 | 0 |
G在B行,R列,B分量的kernel
0 | 0 | +1/2 | 0 | 0 |
---|---|---|---|---|
0 | -1 | 0 | -1 | 0 |
-1 | 4 | 5 | 4 | -1 |
0 | -1 | 0 | -1 | 0 |
0 | 0 | +1/2 | 0 | 0 |
G在R行,B列,B分量的kernel
0 | 0 | -1 | 0 | 0 |
---|---|---|---|---|
0 | -1 | 4 | -1 | 0 |
+1/2 | 0 | 5 | 0 | +1/2 |
0 | -1 | 4 | -1 | 0 |
0 | 0 | -1 | 0 | 0 |
R在R行,R列,B分量的kernel
0 | 0 | -3/2 | 0 | 0 |
---|---|---|---|---|
0 | 2 | 0 | 2 | 0 |
-3/2 | 0 | 6 | 0 | -3/2 |
0 | 2 | 0 | 2 | 0 |
0 | 0 | -3/2 | 0 | 0 |
具体的实现代码为:
function [ rgb_image ] = grad_bilinear_demosaic( cfa_image, cfa_pattern )
%GRAD_BILINEAR_DEMOSAIC Summary of this function goes here
% Detailed explanation goes here
if isequal(cfa_pattern,[1 2 2 3])
rgb_image = grad_bilinear_kernel_rggb(cfa_image);
else
end
end
function [ rgb_image ] = grad_bilinear_kernel_rggb(cfa_image)
% 1. Filter.
kernel_g_rb = [0 0 -1 0 0;
0 0 2 0 0;
-1 2 4 2 -1;
0 0 2 0 0;
0 0 -1 0 0] / 8;
kernel_rb_g_0 = [ 0 0 1/2 0 0;
0 -1 0 -1 0;
-1 4 5 4 -1;
0 -1 0 -1 0;
0 0 1/2 0 0] / 8;
kernel_rb_g_1 = kernel_rb_g_0';
kernel_rb_br = [ 0 0 -3/2 0 0;
0 2 0 2 0;
-3/2 0 6 0 -3/2;
0 2 0 2 0;
0 0 -3/2 0 0] / 8;
padded_cfa = paddarray_symmetric(cfa_image, 2);
filter_cfa_g_rb = filter2(kernel_g_rb, padded_cfa, 'valid');
filter_cfa_rb_g_0 = filter2(kernel_rb_g_0, padded_cfa, 'valid');
filter_cfa_rb_g_1 = filter2(kernel_rb_g_1, padded_cfa, 'valid');
filter_cfa_rb_br = filter2(kernel_rb_br, padded_cfa, 'valid');
% 2. Get the rgb image.
rgb_image = zeros(size(cfa_image,1),size(cfa_image,2),3);
% g channel
rgb_image(:,:,2) = cfa_image;
rgb_image(1:2:end,1:2:end,2) = filter_cfa_g_rb(1:2:end,1:2:end);
rgb_image(2:2:end,2:2:end,2) = filter_cfa_g_rb(2:2:end,2:2:end);
% r channel
rgb_image(:,:,1) = cfa_image;
rgb_image(1:2:end,2:2:end,1) = filter_cfa_rb_g_0(1:2:end,2:2:end);
rgb_image(2:2:end,1:2:end,1) = filter_cfa_rb_g_1(2:2:end,1:2:end);
rgb_image(2:2:end,2:2:end,1) = filter_cfa_rb_br(2:2:end,2:2:end);
% b channel
rgb_image(:,:,3) = cfa_image;
rgb_image(1:2:end,2:2:end,3) = filter_cfa_rb_g_1(1:2:end,2:2:end);
rgb_image(2:2:end,1:2:end,3) = filter_cfa_rb_g_0(2:2:end,1:2:end);
rgb_image(1:2:end,1:2:end,3) = filter_cfa_rb_br(1:2:end,1:2:end);
end
function [ padded_array ] = paddarray_symmetric( input_array, pad_sz )
%PADDARRAY_SYMMETRIC Pad the array using symmetric method.
% Using the board row/column as the symmetry axises.
% Detailed explanation goes here
org_sz = size(input_array);
padded_sz = org_sz + pad_sz + pad_sz;
padded_array = zeros(padded_sz);
padded_array(pad_sz+1:pad_sz+org_sz(1),pad_sz+1:pad_sz+org_sz(2)) = input_array;
% add left
padded_array(:,1:pad_sz) = flip(padded_array(:,pad_sz+2:pad_sz+pad_sz+1),2);
% add right
padded_array(:,end-pad_sz+1:end) = flip(padded_array(:,end-pad_sz-pad_sz:end-pad_sz-1),2);
% add top
padded_array(1:pad_sz,:) = flip(padded_array(pad_sz+2:pad_sz+pad_sz+1,:),1);
% add bottom
padded_array(end-pad_sz+1:end,:) = flip(padded_array(end-pad_sz-pad_sz:end-pad_sz-1,:),1);
end
4. 实验结果
对以上结果的验证,采用如下步骤:
1.将rgb图像读入,按照Bayer Pattern的方式,将图像处理成Bayer Pattern的马赛克图像。
2.用Bilinear Demosaic算法和Gradient-based Bilinear Demosaic算法,将图像插值成rgb图像。
3.将原图,两种算法的结果进行比较。
从实验结果可以看出,Grad-based demosaic相对有更高的对比度和解析度,和Bilinear算法相比,有较大提升。而且Grad-based demosaic可用滤波的方式实现,这使得硬件实现更加容易和高效。