RAW域算法处理之CFA Demosaic

主流的CMOS和CCD sensor几乎都是输出Bayer mosaic格式的RAW数据，这种数据格式是无法直接观看的，必须转换成常见的RGB或YUV格式才能被主流的图像处理软件支持。对于camera产品而言，一般还需要将RGB或YUV图像进一步转换成JPEG格式以方便进行存储。

上述图像处理过程统称图像信号处理（Image Signal Processing，ISP），广义的ISP包含了JPEG和H.264/265图像压缩处理，而狭义的ISP仅包括从RAW格式变换到RGB或YUV的处理过程。

一个典型的ISP流水线由一系列处理模块组成，这些模块首尾相连，在几百MHz的时钟驱动下同时高速运转，图像数据不断从一个模块转移至下一个模块，直到完成所有的算法处理，最终以YUV或RGB的形式从流水线的末级流出ISP。下图所示的是一个支持常见基本功能的ISP流水线。

从图中可以看到，图像数据在ISP内部经历了两次颜色空间变换，第一次变换发生在Demosaic模块，它把像素从即RAW域变换到RGB域，第二次变换发生在CSC模块，它把像素从RGB变到YUV域。下表对ISP 各模块的作用给予了简要说明。

ISP输入图像的格式

目前主流的CMOS sensor几乎都是输出Bayer mosaic格式的RAW数据。Bayer格式图片是伊士曼·柯达公司科学家Bryce Bayer（1929 –2012）发明的，拜耳阵列被广泛运用与数字图像处理领域。

常用的Bayer格式有RGGB,、GRBG、GBRG等多种，因此需要正确配置ISP以反应sensor的数据格式。

图像bayer格式介绍

bayer格式图片是伊士曼·柯达公司科学家Bryce Bayer发明的，Bryce Bayer所发明的拜耳阵列被广泛运用数字图像。

对于彩色图像，需要采集多种最基本的颜色，如rgb三种颜色，最简单的方法就是用滤镜的方法，红色的滤镜透过红色的波长，绿色的滤镜透过绿色的波长，蓝色的滤镜透过蓝色的波长。如果要采集rgb三个基本色，则需要三块滤镜，这样价格昂贵，且不好制造，因为三块滤镜都必须保证每一个像素点都对齐。当用bayer格式的时候，很好的解决了这个问题。bayer 格式图片在一块滤镜上设置的不同的颜色，通过分析人眼对颜色的感知发现，人眼对绿色比较敏感，所以一般bayer格式的图片绿色格式的像素是是r和g像素的和。

另外，Bayer格式是相机内部的原始图片, 一般后缀名为.raw。很多软件都可以查看, 比如PS。我们相机拍照下来存储在存储卡上的.jpeg或其它格式的图片, 都是从.raw格式转化过来的。如下图，为bayer色彩滤波阵列，由一半的G，1/4的R，1/4的B组成。

bayer格式图像传感器硬件

图像传感器的结构如下所示，每一个感光像素之间都有金属隔离层，光纤通过显微镜头，在色彩滤波器过滤之后，投射到相应的漏洞式硅的感光元件上。

当Image Sensor往外逐行输出数据时，像素的序列为GRGRGR.../BGBGBG...（顺序RGB）。这样阵列的Sensor设计，使得RGB传感器减少到了全色传感器的1/3，如下所示。

RAW数据格式

RAW数据的精度常见有8/10/12/14bit等规格。安防监控行业较多使用10/12bit精度的sensor，医疗行业则主要使用12bit以上精度sensor，单反和广电行业则主要使用14bit精度sensor。

当精度高于8bit时，就会自然带来如何存储和表示数据的问题。多数sensor支持固定用16bit数据表示10/12/14bit的有效数据，这种表示方法简单易行，但是在传输过程中会浪费一些带宽，在存储环节会浪费一些存储空间。

对与一些带宽和存储资源特别紧张的场合，有些sensor会支持压缩表示以节约带宽，但是这就需要ISP能够支持相应的压缩格式，否则就需要增加格式适配处理环节，无论是采用硬件还是软件方式实现都会增加系统的复杂度和成本。

另外，很多相机、摄像机产品都会在镜头光路中插入一个红外截止滤光片，其作用是允许波长小于截止频率的光进入系统，而阻止波长大于截止频率的光。常用的截止频率有630nm和650nm两种，它们允许不同程度的红光进入成像系统，因此会对白平衡和颜色校正算法产生影响，所以ISP的参数配置必须要与实际选用的截止频率相配置。

bayer格式插值红蓝算法实现

　　每一个像素仅仅包括了光谱的一部分，必须通过插值来实现每个像素的RGB值。为了从Bayer格式得到每个像素的RGB格式，我们需要通过插值填补缺失的2个色彩。插值的方法有很多（包括领域、线性、3*3等），速度与质量权衡，最好的线性插值补偿算法。其中算法如下：

　　R和B通过线性领域插值，但这有四种不同的分布，如下图所示：

　(a) (b)

　　在（a）与（b）中，R和B分别取邻域的平均值。

　　在（c）与（d）中，取领域的4个B或R的均值作为中间像素的B值。

bayer格式插值绿算法实现

　　由于人眼对绿光反应最敏感，对紫光和红光则反应较弱，因此为了达到更好的画质，需要对G特殊照顾。在上述（c）与（d）中，扩展开来就是上图的（e）与（f）中间像素G的取值，者也有一定的算法要求，不同的算法效果上会有差异。经过相关的研究，

　　（e）中间像素G值的算法如下：

　　（f）中间像素G值的算法如下：

　　CMOS摄像头这部分转换是在内部用ADC或者ISP完成的，生产商为了降低成本必然会使得图像失真。当然用外部处理器来实现转换，如果处理器的速度足够NB，能够胜任像素的操作，用上面的算法来进行转换，皆大欢喜。不过上述算法将直接成倍提高了算法的复杂度，速度上将会有所限制。因此为了速度的提成，可以直接通过来4领域G取均值来中间像素的G值，将会降低一倍的速率，而在性能上差之甚微，算法如下：

　　如果能够通过损失图像的额质量，来达到更快的速度，还可以取G1、G2的均值来实现，但是这样的做法会导致边沿以及跳变部分的失真。

clc;clear;close all;

%% ------------Raw Format----------------
filePath = '../images/kodim19_8bits_RGGB.raw';
bayerFormat = 'RGGB';
width = 512;
% height= 768;
height= 512;
bits = 8;
%% --------------------------------------
file=fopen(filePath, 'r');
format = sprintf('uint8=>uint8');
bayerData = fread(file, [width height], format);
bayerData=bayerData';
%% expand image inorder to make it easy to calculate edge pixels
bayerPadding = zeros(height + 2,width+2);
bayerPadding(2:height+1,2:width+1) = uint32(bayerData);
bayerPadding(1,:) = bayerPadding(3,:);
bayerPadding(height+2,:) = bayerPadding(height,:);
bayerPadding(:,1) = bayerPadding(:,3);
bayerPadding(:,width+2) = bayerPadding(:,width);

%% main code of imterpolation
imDst = zeros(height+2, width+2, 3);
for ver = 2:height + 1
    for hor = 2:width + 1
        switch bayerFormat
            case 'RGGB'
                % G B -> R
                if(0 == mod(ver, 2) && 0 == mod(hor, 2))
                    imDst(ver, hor, 1) = bayerPadding(ver, hor);
                    % G -> R
                    imDst(ver, hor, 2) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor) +...
                                         bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/4;
                    % B -> R
                    imDst(ver, hor, 3) = (bayerPadding(ver-1, hor-1) + bayerPadding(ver-1, hor+1) + ...
                                         bayerPadding(ver+1, hor-1) + bayerPadding(ver+1, hor+1))/4; 
                % G R -> B
                elseif (1 == mod(ver, 2) && 1 == mod(hor, 2))    
                    imDst(ver, hor, 3) = bayerPadding(ver, hor);
                    % G -> B
                    imDst(ver, hor, 2) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor) +...
                                         bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/4;
                    % R -> B
                    imDst(ver, hor, 1) = (bayerPadding(ver-1, hor-1) + bayerPadding(ver-1, hor+1) + ...
                                         bayerPadding(ver+1, hor-1) + bayerPadding(ver+1, hor+1))/4; 
                elseif(0 == mod(ver, 2) && 1 == mod(hor, 2)) %计算G的R和G的B
                    imDst(ver, hor, 2) = bayerPadding(ver, hor);
                    % R -> Gr
                    imDst(ver, hor, 1) = (bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/2;
                    % B -> Gr
                    imDst(ver, hor, 3) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor))/2;
                elseif(1 == mod(ver, 2) && 0 == mod(hor, 2))
                    imDst(ver, hor, 2) = bayerPadding(ver, hor);
                    % B -> Gb
                    imDst(ver, hor, 3) = (bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/2;
                    % R -> Gb
                    imDst(ver, hor, 1) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor))/2;
                end
            case 'GRBG'
                % G B -> R
                if(0 == mod(ver, 2) && 1 == mod(hor, 2))
                    imDst(ver, hor, 1) = bayerPadding(ver, hor);
                    % G -> R
                    imDst(ver, hor, 2) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor) +...
                                         bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/4;
                    % B -> R
                    imDst(ver, hor, 3) = (bayerPadding(ver-1, hor-1) + bayerPadding(ver-1, hor+1) + ...
                                         bayerPadding(ver+1, hor-1) + bayerPadding(ver+1, hor+1))/4; 
                % G R -> B
                elseif (1 == mod(ver, 2) && 0 == mod(hor, 2))    
                    imDst(ver, hor, 3) = bayerPadding(ver, hor);
                    % G -> B
                    imDst(ver, hor, 2) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor) +...
                                         bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/4;
                    % R -> B
                    imDst(ver, hor, 1) = (bayerPadding(ver-1, hor-1) + bayerPadding(ver-1, hor+1) + ...
                                         bayerPadding(ver+1, hor-1) + bayerPadding(ver+1, hor+1))/4; 
                elseif(1 == mod(ver, 2) && 1 == mod(hor, 2))
                    imDst(ver, hor, 2) = bayerPadding(ver, hor);
                    % R -> Gr
                    imDst(ver, hor, 3) = (bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/2;
                    % B -> Gr
                    imDst(ver, hor, 1) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor))/2;
                elseif(0 == mod(ver, 2) && 0 == mod(hor, 2))
                    imDst(ver, hor, 2) = bayerPadding(ver, hor);
                    % B -> Gb
                    imDst(ver, hor, 1) = (bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/2;
                    % R -> Gb
                    imDst(ver, hor, 3) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor))/2;
                end
            case 'GBRG'
                % G B -> R
                if(1 == mod(ver, 2) && 0 == mod(hor, 2))
                    imDst(ver, hor, 1) = bayerPadding(ver, hor);
                    % G -> R
                    imDst(ver, hor, 2) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor) +...
                                         bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/4;
                    % B -> R
                    imDst(ver, hor, 3) = (bayerPadding(ver-1, hor-1) + bayerPadding(ver-1, hor+1) + ...
                                         bayerPadding(ver+1, hor-1) + bayerPadding(ver+1, hor+1))/4; 
                % G R -> B
                elseif (0 == mod(ver, 2) && 1 == mod(hor, 2))    
                    imDst(ver, hor, 3) = bayerPadding(ver, hor);
                    % G -> B
                    imDst(ver, hor, 2) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor) +...
                                         bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/4;
                    % R -> B
                    imDst(ver, hor, 1) = (bayerPadding(ver-1, hor-1) + bayerPadding(ver-1, hor+1) + ...
                                         bayerPadding(ver+1, hor-1) + bayerPadding(ver+1, hor+1))/4; 
                elseif(1 == mod(ver, 2) && 1 == mod(hor, 2))
                    imDst(ver, hor, 2) = bayerPadding(ver, hor);
                    % R -> Gr
                    imDst(ver, hor, 1) = (bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/2;
                    % B -> Gr
                    imDst(ver, hor, 3) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor))/2;
                elseif(0 == mod(ver, 2) && 0 == mod(hor, 2))
                    imDst(ver, hor, 2) = bayerPadding(ver, hor);
                    % B -> Gb
                    imDst(ver, hor, 3) = (bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/2;
                    % R -> Gb
                    imDst(ver, hor, 1) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor))/2;
                end
            case 'BGGR'
                % G R -> B
                if(0 == mod(ver, 2) && 0 == mod(hor, 2))
                    imDst(ver, hor, 3) = bayerPadding(ver, hor);
                    % G -> B
                    imDst(ver, hor, 2) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor) +...
                                         bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/4;
                    % R -> B
                    imDst(ver, hor, 1) = (bayerPadding(ver-1, hor-1) + bayerPadding(ver-1, hor+1) + ...
                                         bayerPadding(ver+1, hor-1) + bayerPadding(ver+1, hor+1))/4; 
                % G B -> R
                elseif (1 == mod(ver, 2) && 1 == mod(hor, 2))    
                    imDst(ver, hor, 1) = bayerPadding(ver, hor);
                    % G -> B
                    imDst(ver, hor, 2) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor) +...
                                         bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/4;
                    % R -> B
                    imDst(ver, hor, 3) = (bayerPadding(ver-1, hor-1) + bayerPadding(ver-1, hor+1) + ...
                                         bayerPadding(ver+1, hor-1) + bayerPadding(ver+1, hor+1))/4; 
                elseif(0 == mod(ver, 2) && 1 == mod(hor, 2))
                    imDst(ver, hor, 2) = bayerPadding(ver, hor);
                    % R -> Gr
                    imDst(ver, hor, 3) = (bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/2;
                    % B -> Gr
                    imDst(ver, hor, 1) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor))/2;
                elseif(1 == mod(ver, 2) && 0 == mod(hor, 2))
                    imDst(ver, hor, 2) = bayerPadding(ver, hor);
                    % B -> Gb
                    imDst(ver, hor, 1) = (bayerPadding(ver, hor-1) + bayerPadding(ver, hor+1))/2;
                    % R -> Gb
                    imDst(ver, hor, 3) = (bayerPadding(ver-1, hor) + bayerPadding(ver+1, hor))/2;
                end
        end
    end
end
imDst = uint8(imDst(2:height+1,2:width+1,:));
imshow(imDst);
figure,imshow(imDst);
title('demosaic image');
figure();imshow(uint8(bayerData));
title('raw image');

RAW域的最后一步处理是Demosaic，将像素从RAW域变换到RGB域进行下一阶段的处理。Demosaic 算法的主要难点在于，RAW域的任何一个像点（photosite）只包含一个真实的采样值，而构成像素（R,G,B）的其它两个值需要从周围像点中预测得到。既然是预测，就一定会发生预测不准的情况，这是不可避免的，而预测不准会带来多种负面影响，包括拉链效应（zipper artifacts），边缘模糊，颜色误差等。

Demosaic 拉链效应和伪彩

Demosaic 伪彩

Demosaic 拉链效应，边缘模糊，伪彩。所以Demosaic 算法的主要挑战就是尽量提高算法的准确性，减少图像边缘损失和颜色误差。

posted @ 2023-05-25 11:26 luckylan 阅读(1083) 评论(0) 编辑收藏举报

刷新页面返回顶部

lanlancky