成像芯片是数码相机的关键部分，生产厂家通常采用单片CCD或CMOS图像传感器以降低数码相机的生产成本和体积，并在其表面覆盖一层彩色滤波阵列CFA，其中Bayer CFA应用最为广泛。CFA使每个成像点只能获得物理三原色（红、绿、蓝）中的一种颜色分量。因此要获得全彩色图像，必须进行插值获得丢失的其余两个分量，这个过程被称为彩色插值或去马赛克DEMOSAICING。论文主要是对DEMOSAICING过程做研究，提出了双线性法和梯度法等四种不同的算法来实现从灰色图像到真彩色图像的还原，然后分别在MATLAB上进行实现。论文分别得出了四种不同算法的彩色还原图像，并通过对原始图像中有代表性的一小块区域进行放大处理，比较和找出了其中的优势算法。结构上，论文首先介绍数字图像处理的一些基础知识，然后回顾了部分典型的插值算法，并选择了在数字图像处理上独具优势的MATLAB软件，详细记述了运用四种不同的算法来实现从灰色图像到真彩色图像还原的全过程，并对图像放大，对比后得出结论。
由于 Bayer 图像是经过CFA采样得到的，需要进行插值才能恢复全彩色图像。然而，高质量的实时插值算法却一直是日本、美国、韩国等厂商产品的高附加值技术。国外对插值算法的研究日趋成熟，然而国内对这方面的研究相对较少。
就目前现状来说，核心问题是如何在保证图像质量的前提下，尽可能地降低插值算法的复杂度。本文主要对插值算法进行研究，力求在保证高质量插值效果的同时控制算法保持在较低的复杂度，以尽可能地提高算法的实用性。
图像传感器利用CFA来检测不同颜色的光的光强，常见的彩色滤波阵列如图1-3所示。
在大多数数码相机里，运用的较多的都是Bayer CFA。其中，R、G、B分别代表红色Red、绿色Green、蓝色Blue。另外在一些比较高端数码相机里，会使用CMY CFA，其中C、Y、M分别代表青色Cyan、黄色Yellow、紫色Magenta。
Bayer CFA根据绿色在可见光谱中的位置最宽和绿色能体现最多的细节这两个特点提出的，它符合人眼对绿色最敏感的这一视觉特性，因此应用最广。关于Bayer CFA的插值算法的研究也最多。Bayer CFA由一组红色和绿色的滤镜与一组绿色和蓝色的滤镜交叉循环组成，采样得到总像素中有一半的是绿，有四分之一是红色，另外四分之一是蓝色。
用Bayer CFA制作的传感器获得的Bayer图像比较昏暗，色彩也并不鲜明，它在每个像素点只有红、绿或蓝一种颜色分量，实际效果如图1-4所示。其中（a）图为传感器获取的理想的Bayer格式的马赛克图片，已经可以从该图片上判断出图像的轮廓外形，只要靠得足够近，就能够比较清楚的看到一个个红、绿、蓝色的点，验证了每个像素只有一种基色；（b）图为（a）中红色线框内的图像放大，可以清晰的观察到Bayer格式的排列规律。
2.1.1 三色原理
在人的视觉系统中，存在着两种感光细胞：杆状细胞和锥状细胞。其中杆状细胞为暗视器官，它们在光照弱时对物体较为敏感，将物体的总体形状形成于视野内。锥状细胞为明视器官，在光照强时起作用。人眼对红、绿、蓝最为敏感，在人眼的大约七百万个锥状细胞中，对红光敏感的占65%，对绿光敏感的占33%，对蓝光敏感的占2%。任何颜色都可通过红、绿、蓝这三种颜色按不同的比例混合而成，同样，任何颜色都可以分解成红、绿、蓝三种颜色，这就是三色原理，也称三原色原理。需要说明的是红绿蓝作为三原色，它们之间是相互独立的，任何一种颜色都不能由其余的两种颜色组成。混合色的饱和度由三原色的比例来决定，其亮度是三种颜色的亮度之和。
2.1.2 颜色的三个属性
颜色是外界光刺激作用于人的视觉器官而产生的一种主观感觉，而人眼的视觉器官只能分辨颜色的色调、色饱和度和亮度三种特征，因此称这三个特征为颜色的三个基本属性。HSI彩色模型正是基于颜色的这三个基本属性提出的。其中色调H(Hue)与光的波长有关，它表示人的感官对不同颜色的感受；饱和度S(Saturation)表示颜色的纯度，饱和度越大，颜色看起来就会越鲜艳；强度I(Intensity)对应成像亮度和图像灰度，是颜色的明亮程度。
常规DEMOSAICING方法简介
在现有文献中提出了大量的彩色插值算法，本文在此介绍两种典型而有效的算法，为第五章MATLAB实现图像彩色还原作准备。
4.1双线性插值法
双线性插值算法（bilinear[8]）是基于Bayer CFA提出来的，它先找出最接近像素的四个图素，然后在它们之间作差补效果。此法属于单通道独立插值算法，简单、易懂，它是设计新型插值算法的基础。
图4-1 Bayer采样阵列
在图4-1中，在每个像素位置只有一种相应的颜色分量，以为例，此处只有蓝色信息，为了得到全彩信息，必须利用相邻位置的红色和绿色分量来恢复出(4,4)位置的红色和绿色分量。双线性插值算法采用3×3滤波器，并取其相邻位置同色分量的平均值，因此可以得到：
这个程序主要包含五大功能模块：一是设置图像的大小,二是打开原始图像，三是对图像边界进行处理，四是对图像主要部分进行处理，五是将图像显示出来。
1）设置图像的大小。已知原始图像的宽度width=1600，高度height=1200，因此在进行图像处理之前，我们先设定图像的大小：w=1600；h=1200。本程序基于Bayer格式，因此，先预设四个1600×1200大小的矩阵：I_data=zeros(h,w)；R=zeros(h,w)；G=zeros(h,w)；B=zeros(h,w)。其中I_data用于存放原始图像的数据，R、G、B分别用于存放Bayer格式图像所需的R、G、B的数据。程序如下：
w=1600;
h=1200;
I_data=zeros(h,w);
R=zeros(h,w);
G=zeros(h,w);
B=zeros(h,w);
2）打开原始图像。由下面小段程序实现：
fid=fopen('00157.raw','r');%%指针指向文件00157.raw，且只读
[A1,cc]=fread(fid,w*h,'uint8');%%打开文件00157.raw并以8bit方式采集数据
A2=uint16(A1);
for i=1:h
for j=1:w
I_data(i,j)=A2((i-1)*w+j);%%将数据赋值给矩阵I
end
end
由于原始图像的格式并不是MATLAB本身能直接处理的数据格式，因此采用了指针fid指向名为00157.raw的文件（这儿就是为何前面要预设路径set path的原因，否则软件没法正确找到所要处理的文件），并用fread()函数将图像的1600×1200个数据以8bit形式存入数列A1中。此时并不符合MATLAB显示图像的格式标准，因此采用了for循环将数据对应的存入I_data=zeros(h,w)矩阵中，此时再用cat()、figure和imshow()函数即可将原始图像显示出来。原始图像如图5-8所示。
3）图像边界处理。由下面小段程序实现：
for i=1:h
for j=1:w
bit_i=bitget(i,1); %%以二进制方式取个数的低位的第一位
bit_j=bitget(j,1); %%用于判断不在边界时点的位置
bit_ij=bit_i*2+ bit_j;
if (i==1)|(i==h)|(j==1)|(j==w)%%此时在图像边界
R(i,j)=A2((i-1)*w+j);%%取边界本身位置的值
G(i,j)=A2((i-1)*w+j);
B(i,j)=A2((i-1)*w+j);
else%%此时是图像主要部分
用(i==1)|(i==h)|(j==1)|(j==w)判断图像的边界。因为已经在图像的边界，外侧没有数据，此时仍然采用3×3的方式，MATLAB软件将会报错，没法处理，因此采用了for循环将位于图像边界的数据直接赋值到R、G、B矩阵中。这并不意味着没有对图像边界上的图像进行处理，而是在处理过程中将位于图像边界上的数据以加权系数为1的方式处理的。
4）图像主要部分的处理。由下面小段程序实现：
else%%此时是图像主要部分
switch(bit_ij)%%判断所需处理的像素所在的位置
case 0%%此时点位于典型位置B22
R(i,j)=(A2((i-2)*w+j-1)+ A2((i-2)*w+j+1)+ A2(i*w+j-1)+ A2(i*w+j+1))/4;
G(i,j)=(A2((i-1)*w+j-1)+ A2((i-1)*w+j+1)+ A2((i-2)*w+j)+ A2(i*w+j))/4;
B(i,j)=A2((i-1)*w+j);
case 1%%此时点位于典型位置G23
R(i,j)=(A2((i-2)*w+j)+ A2(i*w+j))/2;
G(i,j)=A2((i-1)*w+j);
B(i,j)= (A2((i-1)*w+j-1)+ A2((i-1)*w+j+1))/2;
case 2%%此时点位于典型位置G32
R(i,j)=(A2((i-1)*w+j-1)+ A2((i-1)*w+j+1))/2;
G(i,j)=A2((i-1)*w+j);
B(i,j)= (A2((i-2)*w+j)+ A2(i*w+j))/2;
case 3%%此时点位于典型位置R33
R(i,j)=A2((i-1)*w+j);
G(i,j)= (A2((i-1)*w+j-1)+ A2((i-1)*w+j+1)+ A2((i-2)*w+j)+ A2(i*w+j))/4;
B(i,j)= (A2((i-2)*w+j-1)+ A2((i-2)*w+j+1)+ A2(i*w+j-1)+ A2(i*w+j+1))/4;
otherwis
R(i,j)=255;
end
end
end
end
用bit_ij=bit_i*2+ bit_j计算数据的临时值，用switch()函数判断数据在图像中的具体位置，有B22、G23、G32和R33四个典型位置，此时严格按照第4.1节介绍的双线性插值算法处理数据，将得到的数据对应赋值到R、G、B矩阵中。
5）显示图像。由下面小段程序实现：
ff=cat(3,uint8(I_data),uint8(I_data),uint8(I_data));
figure,imshow(ff);%%显示原始图像
ff1=cat(3,uint8(R),uint8(G),uint8(B));
figure,imshow(ff1);%%显示还原得到的彩色图像
经过3×3双线性插值算法还原后的图像如图5-9所示。
5.3.2 3×3梯度法
这个程序也包含五大功能模块。由于是针对同一原始图像做处理，因此其中设置图像的大小、打开原始图像、对图像边界进行处理和显示图像四大模块与3×3双线性法相同，只在对图像主要部分进行处理的部分有区别。
图像主要部分的处理,由下面小段程序实现：
else%%此时是图像主要部分
a(i-1,j-1)=abs(A2((i-1)*w+j-1)-A2((i-1)*w+j+1));%%水平方向梯度
b(i-1,j-1)=abs(A2((i-2)*w+j)-A2(i*w+j));%%垂直方向梯度
switch(bit_ij)%%先插值恢复出G分量
case 1%%G23，此时和case 2处理方法相同
case 2%%G32
G(i,j)=A2((i-1)*w+j);
case 0%%B22，此时和case 3处理方法相同
case 3%%R55
if (a(i-1,j-1)>b(i-1,j-1)) %%判断进行插值的方向
G(i,j)=(A2((i-2)*w+j)+A2(i*w+j))/2;%%垂直方向上进行插值
else if (a(i-1,j-1)<b(i-1,j-1))
G(i,j)=(A2((i-1)*w+j-1)+A2((i-1)*w+j+1))/2; %%水平方向上进行插值
else%%此时水平方向和垂直方向等效
G(i,j)=(A2((i-1)*w+j-1)+ A2((i-1)*w+j+1)+ A2((i-2)*w+j)+ A2(i*w+j))/4;
end
end
end
switch(bit_ij) %%利用G分量恢复出R分量和B分量
case 0%%判断所需处理的点所在的位置，此时位于B22
R(i,j)=(A2((i-2)*w+j-1)+A2((i-2)*w+j+1)+A2(i*w+j-1)+A2(i*w+j+1))/4+(G(i-1,j-1)+G(i-1,j+1)+G(i+1,j-1)+G(i+1,j+1))/4-G(i,j);
B(i,j)=A2((i-1)*w+j);
case 1%%G23
R(i,j)=(A2((i-2)*w+j)+ A2(i*w+j))/2+(G(i-1,j)+G(i+1,j))/2-G(i,j);
B(i,j)=(A2((i-1)*w+j-1)+A2((i-1)*w+j+1))/2+(G(i,j-1)+G(i,j+1))/2-G(i,j);
case 2%%G32
R(i,j)=(A2((i-1)*w+j-1)+A2((i-1)*w+j+1))/2+(G(i,j-1)+G(i,j+1))/2-G(i,j);
B(i,j)=(A2((i-2)*w+j)+A2(i*w+j))/2+(G(i-1,j)+G(i+1,j))/2-G(i,j);
case 3%%R33
R(i,j)=A2((i-1)*w+j);
B(i,j)=(A2((i-2)*w+j-1)+A2((i-2)*w+j+1)+A2(i*w+j-1)+A2(i*w+j+1))/4+(G(i-1,j-1)+G(i-1,j+1)+G(i+1,j-1)+G(i+1,j+1))/4-G(i,j);
otherwise
R(i,j)=255;
end
end
end
end
此时严格按照第4.2节介绍的基于梯度的插值算法处理数据，将得到的数据对应赋值到R、G、B矩阵中。经过插值还原后的图像如图5-10所示。
5.3.3 5×5双线性法
这个程序也包含五大功能模块。其中设置图像的大小、打开原始图像和显示图像三大模块与3×3双线性法相同，只在对图像边界和主要部分进行处理的部分有区别。
1）图像边界处理。由下面小段程序实现：
if (i==1)|(i==h)|(j==1)|(j==w) %%此时在图像5×5的边界
R(i,j)=A2((i-1)*w+j); %%取边界本身位置的值
G(i,j)=A2((i-1)*w+j);
B(i,j)=A2((i-1)*w+j);
else if (((i==2)|(i==h-1))&((j~=1)|(j~=w)))|(((j==2)|(j==w-1))&((i~=1)|(i~=h))) %%此时点在图像3×3的边界
switch(bit_ij) %%判断所需处理的像素所在的位置
case 0%%B22%%按3×3方法取值
R(i,j)=(A2((i-2)*w+j-1)+ A2((i-2)*w+j+1)+ A2(i*w+j-1)+ A2(i*w+j+1))/4;
G(i,j)=(A2((i-1)*w+j-1)+ A2((i-1)*w+j+1)+ A2((i-2)*w+j)+ A2(i*w+j))/4;
B(i,j)=A2((i-1)*w+j);
case 1%%G23
R(i,j)=(A2((i-2)*w+j)+ A2(i*w+j))/2;
G(i,j)=(A2((i-1)*w+j)+(A2((i-2)*w+j-1)+A2((i-2)*w+j+1)+A2(i*w+j-1)+A2(i*w+j+1))/4)/2;
B(i,j)=(A2((i-1)*w+j-1)+ A2((i-1)*w+j+1))/2;
case 2%%G32
R(i,j)=(A2((i-1)*w+j-1)+ A2((i-1)*w+j+1))/2;
G(i,j)=(A2((i-1)*w+j)+(A2((i-2)*w+j-1)+A2((i-2)*w+j+1)+A2(i*w+j-1)+A2(i*w+j+1))/4)/2;
B(i,j)=(A2((i-2)*w+j)+A2(i*w+j))/2;
case 3%%R55
R(i,j)=A2((i-1)*w+j);
G(i,j)=(A2((i-1)*w+j-1)+ A2((i-1)*w+j+1)+ A2((i-2)*w+j)+ A2(i*w+j))/4;
B(i,j)=(A2((i-2)*w+j-1)+ A2((i-2)*w+j+1)+ A2(i*w+j-1)+ A2(i*w+j+1))/4;
end
else%%此时是图像主要部分
由于此时采用5×5双线性法，因此图像将会在5×5和3×3位置上出现两个边界。在5×5位置取边界本身位置的值，在3×3位置按3×3方法取值。
2）图像主要部分的处理。由下面小段程序实现：
else%%此时是图像主要部分
switch(bit_ij) %%判断所需处理的像素所在的位置
case 0%%B44，按5×5方法取值
R(i,j)=(A2((i-2)*w+j-1)+ A2((i-2)*w+j+1)+ A2(i*w+j-1)+ A2(i*w+j+1))/4;
G(i,j)=(A2((i-1)*w+j-1)+A2((i-1)*w+j+1)+A2((i-2)*w+j)+A2(i*w+j))/8+(A2((i-3)*w+j-1)+A2((i-3)*w+j+1)+A2((i-2)*w+j-2)+A2((i-2)*w+j+2)+A2(i*w+j-2)+A2(i*w+j+2)+A2((i+1)*w+j-1)+A2((i+1)*w+j+1))/16;
B(i,j)=A2((i-1)*w+j)/2+(A2((i-3)*w+j)+A2((i-1)*w+j-2)+A2((i-1)*w+j+2)+A2((i+1)*w+j)+A2((i-3)*w+j-2)+A2((i-3)*w+j+2)+A2((i+1)*w+j-2)+A2((i+1)*w+j+2))/16;
case 1%%(G43)
R(i,j)=(A2((i-2)*w+j)+A2(i*w+j))/4+(A2((i-2)*w+j-2)+A2((i-2)*w+j+2)+A2(i*w+j-2)+A2(i*w+j+2))/8;
G(i,j)=(A2((i-1)*w+j)+(A2((i-2)*w+j-1)+A2((i-2)*w+j+1)+A2(i*w+j-1)+A2(i*w+j+1))/4+(A2((i-3)*w+j-2)+A2((i-3)*w+j)+A2((i-3)*w+j+2)+A2((i-1)*w+j-2)+A2((i-1)*w+j+2)+A2((i+1)*w+j-2)+A2((i+1)*w+j)+A2((i+1)*w+j+2))/8)/3;
B(i,j)=((A2((i-1)*w+j-1)+A2((i-1)*w+j+1))/2+(A2((i-3)*w+j-1)+A2((i-3)*w+j+1)+A2((i+1)*w+j-1)+A2((i+1)*w+j+1))/4)/2;
case 2%%(G34)
R(i,j)=((A2((i-1)*w+j-1)+A2((i-1)*w+j+1))/2+(A2((i-3)*w+j-1)+A2((i-3)*w+j+1)+A2((i+1)*w+j-1)+A2((i+1)*w+j+1))/4)/2;
G(i,j)=(A2((i-1)*w+j)+(A2((i-2)*w+j-1)+A2((i-2)*w+j+1)+A2(i*w+j-1)+A2(i*w+j+1))/4+(A2((i-3)*w+j-2)+A2((i-3)*w+j)+A2((i-3)*w+j+2)+A2((i-1)*w+j-2)+A2((i-1)*w+j+2)+A2((i+1)*w+j-2)+A2((i+1)*w+j)+A2((i+1)*w+j+2))/8)/3;
B(i,j)=((A2((i-2)*w+j)+A2(i*w+j))/2+(A2((i-2)*w+j-2)+A2(i*w+j-2)+A2((i-2)*w+j+2)+A2(i*w+j+2))/4)/2;
case 3%%R33
R(i,j)=(A2((i-1)*w+j)+(A2((i-3)*w+j-2)+A2((i-3)*w+j)+A2((i-3)*w+j+2)+A2((i-1)*w+j-2)+A2((i-1)*w+j+2)+A2((i+1)*w+j-2)+A2((i+1)*w+j)+A2((i+1)*w+j+2))/8)/2;
G(i,j)=((A2((i-1)*w+j-1)+A2((i-1)*w+j+1)+A2((i-2)*w+j)+A2(i*w+j))/4+(A2((i-3)*w+j-1)+A2((i-3)*w+j+1)+A2((i-2)*w+j-2)+A2((i-2)*w+j+2)+A2(i*w+j-2)+A2(i*w+j+2)+A2((i+1)*w+j-1)+A2((i+1)*w+j+1))/8)/2;
B(i,j)=(A2((i-2)*w+j-1)+ A2((i-2)*w+j+1)+ A2(i*w+j-1)+ A2(i*w+j+1))/4;
otherwise
R(i,j)=255;
end
end
end
end
end
此时基于第4.1节介绍的双线性插值算法，并将其由3×3扩张到5×5进行处理数据，将得到的数据对应赋值到R、G、B矩阵中。经过插值还原后的图像如图5-11所示。
5.3.4 5×5不同系数双线性法
这个程序也包含五大功能模块。其中设置图像的大小、打开原始图像、对图像边界处理和显示图像四大模块与5×5双线性法相同，只在对图像主要部分进行处理有区别。
图像主要部分的处理,由下面小段程序实现：
else%%此时是图像主要部分
switch(bit_ij) %%判断所需处理的像素所在的位置
case 0%%B44，按5×5方法取值，只是加权系数与它不同
R(i,j)=(A2((i-2)*w+j-1)+ A2((i-2)*w+j+1)+ A2(i*w+j-1)+ A2(i*w+j+1))/4;
G(i,j)=(A2((i-1)*w+j-1)+A2((i-1)*w+j+1)+A2((i-2)*w+j)+A2(i*w+j)+A2((i-3)*w+j-1)+A2((i-3)*w+j+1)+A2((i-2)*w+j-2)+A2((i-2)*w+j+2)+A2(i*w+j-2)+A2(i*w+j+2)+A2((i+1)*w+j-1)+A2((i+1)*w+j+1))/12;
B(i,j)=(A2((i-1)*w+j)+A2((i-3)*w+j)+A2((i-1)*w+j-2)+A2((i-1)*w+j+2)+A2((i+1)*w+j)+A2((i-3)*w+j-2)+A2((i-3)*w+j+2)+A2((i+1)*w+j-2)+A2((i+1)*w+j+2))/9;
case 1%%(G43)
R(i,j)=(A2((i-2)*w+j)+A2(i*w+j)+A2((i-2)*w+j-2)+A2((i-2)*w+j+2)+A2(i*w+j-2)+A2(i*w+j+2))/6;
G(i,j)=(A2((i-1)*w+j)+A2((i-2)*w+j-1)+A2((i-2)*w+j+1)+A2(i*w+j-1)+A2(i*w+j+1)+A2((i-3)*w+j-2)+A2((i-3)*w+j)+A2((i-3)*w+j+2)+A2((i-1)*w+j-2)+A2((i-1)*w+j+2)+A2((i+1)*w+j-2)+A2((i+1)*w+j)+A2((i+1)*w+j+2))/13;
B(i,j)=(A2((i-1)*w+j-1)+A2((i-1)*w+j+1)+A2((i-3)*w+j-1)+A2((i-3)*w+j+1)+A2((i+1)*w+j-1)+A2((i+1)*w+j+1))/6;
case 2%%(G34)
R(i,j)=(A2((i-1)*w+j-1)+A2((i-1)*w+j+1)+A2((i-3)*w+j-1)+A2((i-3)*w+j+1)+A2((i+1)*w+j-1)+A2((i+1)*w+j+1))/6;
G(i,j)=(A2((i-1)*w+j)+A2((i-2)*w+j-1)+A2((i-2)*w+j+1)+A2(i*w+j-1)+A2(i*w+j+1)+A2((i-3)*w+j-2)+A2((i-3)*w+j)+A2((i-3)*w+j+2)+A2((i-1)*w+j-2)+A2((i-1)*w+j+2)+A2((i+1)*w+j-2)+A2((i+1)*w+j)+A2((i+1)*w+j+2))/13;
B(i,j)=(A2((i-2)*w+j)+A2(i*w+j)+A2((i-2)*w+j-2)+A2(i*w+j-2)+A2((i-2)*w+j+2)+A2(i*w+j+2))/6;
case 3%%R33
R(i,j)=(A2((i-1)*w+j)+A2((i-3)*w+j-2)+A2((i-3)*w+j)+A2((i-3)*w+j+2)+A2((i-1)*w+j-2)+A2((i-1)*w+j+2)+A2((i+1)*w+j-2)+A2((i+1)*w+j)+A2((i+1)*w+j+2))/9;
G(i,j)=(A2((i-1)*w+j-1)+A2((i-1)*w+j+1)+A2((i-2)*w+j)+A2(i*w+j)+A2((i-3)*w+j-1)+A2((i-3)*w+j+1)+A2((i-2)*w+j-2)+A2((i-2)*w+j+2)+A2(i*w+j-2)+A2(i*w+j+2)+A2((i+1)*w+j-1)+A2((i+1)*w+j+1))/9;
B(i,j)=(A2((i-2)*w+j-1)+ A2((i-2)*w+j+1)+ A2(i*w+j-1)+ A2(i*w+j+1))/4;
otherwise
R(i,j)=255;
end
end
end
end
end
此时处理数据的方法与“5.3.3 5×5双线性法”完全相同，只是所选取的系数不同，将得到的数据对应赋值到R、G、B矩阵中。经过插值还原后的图像如图5-12所示。