ISP--Demosaicking

前言

人眼之所以有能感受到自然界的颜色，是因为人眼的感光细胞中有三种锥体细胞对红绿蓝三种颜色敏感，所以我们就可以通过RGB三种颜色来表示一个颜色空间，通过这个颜色空间中的点就能表示自然界中所有的颜色。那么数码相机只要能类似人一样获取自然界中的这三个分量，那么就能复现人眼看到的颜色。相机系统用的感光器件只是一个光电转换器件，所以感光器件只对亮度分量敏感，无法感知颜色，所以需要通过滤光片将光线分解成RGB三个分量然后再用感光器件去接受。那么最直接的方式就是用三个滤光片分别过滤出RGB三个通道的分量，然后用三个感光器件去分别接受三个通道的强度，然后再将这三个通道的值叠加到一起就能复现出正常的颜色。这种设计称为3CCD，这种方式大概可以用如下简图表示。但是这种方式工艺复杂且成本较高，所以到目前位置在消费领域基本没有这种操作。

由于其缺陷，所以柯达公司的科学家Bryce Bayer（1929-2012）发明了一个突破性的解决方案。这个方案不需要使用昂贵的光学棱镜，也不需要使用3个CCD阵列，只需要在一个CCD阵列上制造三种不同的滤光膜，构成一个滤光膜阵列（Color Filter Array，CFA），就形成一个廉价而高效的解决方案。这种方式就是在感光器件上面通过交替的滤光透镜过滤出三中颜色分量形成如图的RGB三色交替的图像。因为自然界中绿色占比较大，所以给G分配了两个位置。这样每个像素点就只有一种颜色，那如何得到每个像素点的所有颜色呢？所以就有了去马赛克算法：后期再通过一定的算法通过周边的颜色恢复出确实的颜色，最终形成RGB的颜色，这种后期的处理方式就是本文讨论的重点，它实现了把图像从raw域转换到了rgb域。


参考：Demosiac 去马赛克 （CIP）

算法解释

简单的线性插值

如图是Bayer格式的raw图，RGB三种颜色交替覆盖，且绿色分量是RB分量的两倍。由于这种特殊的分布方式，所以可以通过最简单的线性插值的方式通过附近已知的颜色插值出同一通道缺失的分量。例如途中G13点为绿色，缺失了R和B，但是G13点为绿色，但是G13左右的R是已知的，那么就可以通过左右红色分量插值出该点缺失的红色，同理可以通过上下的蓝色分量插值出缺失的蓝色分量。对于R14缺失的G和B，但是周围相邻的有4个G是已知的，就可以通过这四个点插值出该点缺失的G，同理可以通过周围四个B插值出该点的B。

代码实现

%% --------------------------------
%% fuction: main file of Demosaic. The simple linear interpolation.
%% note: add RGGB format only, other formats will be added later
%% --------------------------------
clc;clear;close all;

%% ------------Raw Format----------------
filePath = 'images/kodim19_8bits_RGGB.raw';
bayerFormat = 'RGGB';
width = 512;
height= 768;
bits = 8;
%% --------------------------------------
bayerData = readRaw(filePath, bits, width, height);
figure();
imshow(bayerData);
title('raw image');

%% 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))
                    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'
                continue;
            case 'GBGR'
                continue;
            case 'BGGR'
                continue;
        end
    end
end
imDst = uint8(imDst(2:height+1,2:width+1,:));
figure,imshow(imDst);title('demosaic image');

orgImage = imread('images/kodim19.png');
figure, imshow(orgImage);title('org image');

左侧是实验原图，右侧是通过上述的简单的线性插值出来的效果，从图中可以看出简单的线性插值出来得到效果一般。存在画面整体清晰度变差，高频存在伪彩，边缘又伪像。

色差法和色比法

色差法和色比法其实是基于两个假设实现插值的。**色差恒定认为在一定范围内入射光和法线是比较接近的，因此可以认为在一定范围内的相邻点像素之间符合R-G或者B-G是相等的，由此产生了色差恒定的假设。**其中色比法假设在一个邻域范围内不同颜色通道的比值是固定的，简单来说就是相邻像素的R/G的值和B/G的值是一样的，那么设计算法时就可以利用这一点。一般情况下都会先插值出G的缺失值，然后通过与G的比值恒定插值出其他的缺失值。如上展示的RAW图，可以通过G9，G13， G15，G19做简单的线性插值恢复G14，然后通过R14/G14 = R13/G13的假设恢复出R13。同理可以恢复出其他缺失的颜色。
色差法和色比法类似，色差法假设在一个邻域内不同通道的颜色插值时恒定的。只是将色比法的比值转换为差值即可。

具体流程，这里我直接手绘了：

原理参考：https://zhuanlan.zhihu.com/p/592335382

代码实现：

%% --------------------------------
%% fuction: main file of Demosaic. The simple linear interpolation.
%% note: add RGGB format only, other formats will be added later
%% --------------------------------
clc;clear;close all;

%% ------------Raw Format----------------
filePath = 'images/kodim19_8bits_RGGB.raw';
bayerFormat = 'RGGB';
width = 512;
height= 768;
bits = 8;
%% --------------------------------------
bayerData = readRaw(filePath, bits, width, height);
figure();
imshow(bayerData);
title('raw image');

%% expand image inorder to make it easy to calculate edge pixels
bayerPadding = zeros(height + 4, width + 4);
bayerPadding(3:height+2, 3:width+2) = uint32(bayerData);

% 正确镜像填充（上下方向）
bayerPadding(1:2, 3:width+2) = bayerPadding(5:-1:4, 3:width+2); % 顶部镜像原图第3-4行
bayerPadding(height+3:height+4, 3:width+2) = bayerPadding(height+1:-1:height, 3:width+2); % 底部镜像

% 正确镜像填充（左右方向）
bayerPadding(3:height+2, 1:2) = bayerPadding(3:height+2, 5:-1:4); % 左侧镜像
bayerPadding(3:height+2, width+3:width+4) = bayerPadding(3:height+2, width+1:-1:width); % 右侧镜像

%% main code of imterpolation
imDst = zeros(height+4, width+4, 3); %创建三通道图像
for ver = 3:height + 2
    for hor = 3:width + 2
        Dh = bayerPadding(ver-1, hor)+bayerPadding(ver+1, hor);
        Dv = bayerPadding(ver, hor+1)+bayerPadding(ver, hor-1);
        switch bayerFormat
            case 'RGGB'
                % R
                if(1 == mod(ver, 2) && 1 == mod(hor, 2))   
                    imDst(ver, hor, 1) = bayerPadding(ver, hor); %直接填充红色通道
                    % G -> R 求均值，插值R上的G
                    imDst(ver, hor, 2) = (Dv+Dh)/4;
                % B通道
                elseif (0 == mod(ver, 2) && 0 == mod(hor, 2))    
                    imDst(ver, hor, 3) = bayerPadding(ver, hor);
                    % G -> B 插值B上的G
                    imDst(ver, hor, 2) = (Dv+Dh)/4;  
                elseif(1 == mod(ver, 2) && 0 == mod(hor, 2))
                    imDst(ver, hor, 2) = bayerPadding(ver, hor);   
                elseif(0 == mod(ver, 2) && 1 == mod(hor, 2))
                    imDst(ver, hor, 2) = bayerPadding(ver, hor);   
                end
        end
    end
end
for ver = 3:height + 2
    for hor = 3:width+2
        Dh = bayerPadding(ver-1, hor)+bayerPadding(ver+1, hor);
        switch bayerFormat
            case 'RGGB'
                %R通道
                if(1==mod(ver,2)&&(1==mod(hor,2)))
                    %更新R上的G，插值R上的B：
                    imDst(ver,hor,2) = (2*bayerPadding(ver, hor)-bayerPadding(ver, hor-2)-bayerPadding(ver, hor+2))/4 + Dh/2;
                    imDst(ver,hor,3) = imDst(ver, hor, 2)  +( -imDst(ver-1, hor-1, 2) + bayerPadding(ver-1, hor-1) -imDst(ver-1, hor+1, 2)+ bayerPadding(ver-1, hor+1)-...
imDst(ver+1, hor-1, 2) +bayerPadding(ver+1, hor-1) -imDst(ver+1, hor+1, 2) +bayerPadding(ver+1, hor+1) )/4;
                %B通道
                elseif(0==mod(ver,2)&&0==mod(hor,2))
                    imDst(ver,hor,2) = (2*bayerPadding(ver, hor)-bayerPadding(ver, hor-2)-bayerPadding(ver, hor+2))/4 + Dh/2;
                    imDst(ver, hor, 1) = imDst(ver, hor, 2)  +( -imDst(ver-1, hor-1, 2) + bayerPadding(ver-1, hor-1) -imDst(ver-1, hor+1, 2)+ bayerPadding(ver-1, hor+1)-...
imDst(ver+1, hor-1, 2) +bayerPadding(ver+1, hor-1) -imDst(ver+1, hor+1, 2) +bayerPadding(ver+1, hor+1) )/4;
                elseif(1==mod(ver,2)&&0==mod(hor,2))
                    imDst(ver, hor, 1) = imDst(ver,hor,2)+ ( bayerPadding(ver, hor-1)- imDst(ver, hor-1, 2) - imDst(ver, hor+1, 2) + bayerPadding(ver, hor+1) )/2;
                    imDst(ver, hor, 3) = imDst(ver,hor,2) + ( bayerPadding(ver-1, hor)- imDst(ver-1, hor, 2) - imDst(ver+1, hor, 2) + bayerPadding(ver+1, hor) )/2;
                elseif(0==mod(ver,2)&&1==mod(hor,2))
                    imDst(ver, hor, 3) =imDst(ver,hor,2) + ( bayerPadding(ver, hor-1)- imDst(ver, hor-1, 2) - imDst(ver, hor+1, 2) + bayerPadding(ver, hor+1) )/2;
                    imDst(ver, hor, 1) = imDst(ver,hor,2) + ( bayerPadding(ver-1, hor)- imDst(ver-1, hor, 2) - imDst(ver+1, hor, 2) + bayerPadding(ver+1, hor) )/2;
                end
        end
    end

end
imDst = uint8(imDst(3:height+2, 3:width+2, :)); % 正确裁剪
figure,imshow(imDst);title('demosaic image');

orgImage = imread('images/kodim19.png');
figure, imshow(orgImage);title('org image');

效果：

伪彩还是存在的。G值的精确更新可以减少伪彩。

基于方向加权的方法

基本原理参考：小话demosaic （二）
通过分析边缘方向来选择插值方向，从而减少边缘的模糊，保留边缘信息。

RB缺失的G通道的插值

分别计算出各个方向的梯度值，计算公式如下：

例对于B44来说：对于B上的G插值

$K_n$是量化邻域像素间的颜色差异的系数。Kn值越小，表明该方向颜色变化平缓，插值可靠性高；反之则可能存在边缘，需降低权重。
各个梯度上的权重计算如下：
$$W_{n(i,j)}=\frac{1}{(1+I_{n(i,j)})}/\sum _{n=1}^{12}\frac{1}{(1+I_{n(i,j)})}$$
在边缘区域，沿边缘方向的Kn差异小，对应Wn更大；跨边缘方向Kn差异大，Wn趋近于零。
则插值的$G_{i,j}$：
$$G_{(i,j)}=B_{(i,j)}+\sum _{n=1}^{12}{W}_{n(i,j)}\ast K_{b,n}(i+v_n,j+h_n)$$
$$K_{b,n}(i+v_n,j+h_n)=G_{(}i+v_n,j+h_n)-B_{(}i+v_n,j+h_n)$$
即$G_{(i,j)}-B_{(i,j)}$=色差，B则是上下或者 左右的均值。

RB缺失的BR的插值

对于B上的R插值:
同理求四个方向的梯度值:

各个梯度上的权重计算如下：
$$W_{n(i,j)}=\frac{1}{1+I_{n(i,j)}}/\sum _{n=1}^4\frac{1}{1+I_{n(i,j)}}$$
则:
$$R_{(i,j)}=G_{(i,j)}-\sum _{n=1}^4{W}_{n(i,j)}\ast K_{rn}(i+v_n^{\prime },j+h_n^{\prime })$$
$$K_{rn}(i+v_n^{\prime },j+h_n^{\prime })=G_{(}i+v_n^{\prime },j+h_n^{\prime })-R_{(}i+v_n^{\prime },j+h_n^{\prime })$$
同理R则是上下或者 左右的均值。

G缺失的BR的插值

和BR缺失的G插值是一样的:

代码实现

%% --------------------------------
%% reference: "Directionally weighted color interpolation for digital cameras"
%% --------------------------------
clc;clear;close all;

%% ------------Raw Format----------------
filePath = 'images/kodim19_8bits_RGGB.raw';
bayerFormat = 'RGGB';
width = 512;
height= 768;
bits = 8;

%% ------------Global Value--------------
RC = 1;
GC = 2;
BC = 3;

%% --------------------------------------
orgImg = imread('images/kodim19.png');
figure();imshow(orgImg);title('org image');

bayerData = readRaw(filePath, bits, width, height);
figure();
imshow(bayerData);
title('raw image');

%% expand image inorder to make it easy to calculate edge pixels
addpath('../publicFunction');
bayerPadding = expandRaw(bayerData, 4);

imDst = zeros(height+8, width+8, 3); 

%% Interpolate the missing green value of blue/red samples 计算插值的 r和b差的g
for ver = 5: height + 4
    for hor = 5: width +4
        % R channal
        if(1 == mod(ver, 2) && 1 == mod(hor, 2))
            imDst(ver, hor, 1) = bayerPadding(ver, hor);
            neighborhoodData = bayerPadding(ver-4: ver+4, hor-4: hor+4);
            Wn = DW_Wn(neighborhoodData, 12);
            Kn = DW_Kn(neighborhoodData, 12);
            imDst(ver, hor, 2) = bayerPadding(ver, hor) + sum(Wn .* Kn);
        % B channal
        elseif (0 == mod(ver, 2) && 0 == mod(hor, 2))
            imDst(ver, hor, 3) = bayerPadding(ver, hor);
            neighborhoodData = bayerPadding(ver-4: ver+4, hor-4: hor+4);
            Wn = DW_Wn(neighborhoodData, 12);
            Kn = DW_Kn(neighborhoodData, 12);
            imDst(ver, hor, 2) = bayerPadding(ver, hor) + sum(Wn .* Kn);
        % Gr
        elseif (1 == mod(ver, 2) && 0 == mod(hor, 2))
            imDst(ver, hor, 2) = bayerPadding(ver, hor);
        % Gb
        elseif (0 == mod(ver, 2) && 1 == mod(hor, 2))
            imDst(ver, hor, 2) = bayerPadding(ver, hor);
        end
    end
end

% expand the imDst
imDst(:, 1: 4, :) = imDst(:, 5: 8, :);
imDst(:, width+5: width+8, :) = imDst(:, width+1: width+4, :);
imDst(1:4, : , :) = imDst(5: 8, :, :);
imDst(height+5: height+8, : , :) = imDst(height+1: height+4, :, :);

%% Interpolate the missing red/blue values of blue/red samples.
for ver = 5: height + 4
    for hor = 5: width +4
        % R channal
        if(1 == mod(ver, 2) && 1 == mod(hor, 2))
            neighborRaw = bayerPadding(ver-4: ver+4, hor-4: hor+4);
            Wn = DW_Wn(neighborRaw, 4);
            neighborhoodData = imDst(ver-4: ver+4, hor-4: hor+4, :);
            Kn = DW_Kn(neighborhoodData, 4, BC);
            imDst(ver, hor, 3) = imDst(ver, hor, 2) - sum(Wn .* Kn);
        % B channal
        elseif (0 == mod(ver, 2) && 0 == mod(hor, 2))
            neighborRaw = bayerPadding(ver-4: ver+4, hor-4: hor+4);
            Wn = DW_Wn(neighborRaw, 4);
            neighborhoodData = imDst(ver-4: ver+4, hor-4: hor+4, :);
            Kn = DW_Kn(neighborhoodData, 4, RC);
            imDst(ver, hor, 1) = imDst(ver, hor, 2) - sum(Wn .* Kn);
        else
            continue
        end
    end
end
% expand the imDst
imDst(:, 1: 4, :) = imDst(:, 5: 8, :);
imDst(:, width+5: width+8, :) = imDst(:, width+1: width+4, :);
imDst(1:4, : , :) = imDst(5: 8, :, :);
imDst(height+5: height+8, : , :) = imDst(height+1: height+4, :, :);

%% Interpolate missing red/blue values of green samples.
for ver = 5: height + 4
    for hor = 5: width +4
        neighborhoodData = imDst(ver-4: ver+4, hor-4: hor+4, :);
        % R channal
        if(1 == mod(ver, 2) && 1 == mod(hor, 2))
            continue
        % B channal
        elseif (0 == mod(ver, 2) && 0 == mod(hor, 2))
            continue
        % G
        else
            Wrn = DW_Wn(neighborhoodData, 12, GC, RC);
            Wbn = DW_Wn(neighborhoodData, 12, GC, BC);
            Krn = DW_Kn(neighborhoodData, 12, RC);
            Kbn = DW_Kn(neighborhoodData, 12, BC);
            imDst(ver, hor, 1) = imDst(ver, hor, 2) - sum(Wrn .* Krn);
            imDst(ver, hor, 3) = imDst(ver, hor, 2) - sum(Wbn .* Kbn);
        end
    end
end
% expand the imDst
imDst(:, 1: 4, :) = imDst(:, 5: 8, :);
imDst(:, width+5: width+8, :) = imDst(:, width+1: width+4, :);
imDst(1:4, : , :) = imDst(5: 8, :, :);
imDst(height+5: height+8, : , :) = imDst(height+1: height+4, :, :);
figure();imshow(uint8(imDst));title('now');

%% Adjust the estimated green values of red/blue samples
for ver = 5: height + 4
    for hor = 5: width +4
        neighborhoodData = imDst(ver-4: ver+4, hor-4: hor+4, :);
        % R channal
        if(1 == mod(ver, 2) && 1 == mod(hor, 2))
            Wrn = DW_Wn(neighborhoodData, 12, RC, GC);           
            Krn = DW_Kn(neighborhoodData, 12, RC);
            imDst(ver, hor, 2) = imDst(ver, hor, 1) + sum(Wrn .* Krn);
        % B channal
        elseif (0 == mod(ver, 2) && 0 == mod(hor, 2))
            Wbn = DW_Wn(neighborhoodData, 12, BC, GC);
            Kbn = DW_Kn(neighborhoodData, 12, BC);
            imDst(ver, hor, 2) = imDst(ver, hor, 3) + sum(Wbn .* Kbn);
        % G
        else
            continue
        end
    end
end
demosaicImg = imDst(5: height + 4, 5: width +4, :);
figure();imshow(uint8(demosaicImg));title('demosaicking image');

Wn计算：这里的1-4方向kn=0.5，up主是根据插值点的距离来选择的，1-4方向近，权重越大，In越小。(存疑？

function Wn = DW_Wn(neighborhoodData, directionNum, channelDeal, channelAdded)
% DW_Wn.m    get rawData from HiRawImage
%   Input:
%       neighborhoodData    the data of neighborhood range
%       directionNum        the num of direction
%       channalAdded        the channel need to be added
%   Output:
%       Wn                  The weignt of each direction
% Note: 
if nargin < 3
    channelDeal = 1;
    channelAdded = 1;
end
In = zeros(directionNum, 1);
Wn = zeros(directionNum, 1);
[h, w, c] = size(neighborhoodData);
centerH = round(h/2);
centerW = round(w/2);
sumIn = 0;
k=4;
switch directionNum
    case 12
        In(1) = 0.5*(abs(neighborhoodData(centerH, centerW-1, channelAdded) - neighborhoodData(centerH, centerW+1, channelAdded)) + ...
                abs(neighborhoodData(centerH, centerW-2, channelDeal) - neighborhoodData(centerH, centerW, channelDeal)));

        In(2) =  0.5*(abs(neighborhoodData(centerH-1, centerW, channelAdded) - neighborhoodData(centerH+1, centerW, channelAdded)) + ...
                abs(neighborhoodData(centerH-2, centerW, channelDeal) - neighborhoodData(centerH, centerW, channelDeal)));

        In(3) =  0.5*(abs(neighborhoodData(centerH, centerW+1, channelAdded) - neighborhoodData(centerH, centerW-1, channelAdded)) + ...
                abs(neighborhoodData(centerH, centerW+2, channelDeal) - neighborhoodData(centerH, centerW, channelDeal)));

        In(4) =  0.5*(abs(neighborhoodData(centerH+1, centerW, channelAdded) - neighborhoodData(centerH-1, centerW, channelAdded)) + ...
                abs(neighborhoodData(centerH+2, centerW) - neighborhoodData(centerH, centerW)));    

        In(5) = (abs(neighborhoodData(centerH-1, centerW-2, channelAdded) - neighborhoodData(centerH+1, centerW+2, channelAdded)) + ...
                     abs(neighborhoodData(centerH-2, centerW-4, channelDeal) - neighborhoodData(centerH, centerW, channelDeal)));

        In(6) = (abs(neighborhoodData(centerH-2, centerW-1, channelAdded) - neighborhoodData(centerH+2, centerW+1, channelAdded)) + ...
                     abs(neighborhoodData(centerH-4, centerW-2, channelDeal) - neighborhoodData(centerH, centerW, channelDeal)));

        In(7) = (abs(neighborhoodData(centerH-2, centerW+1, channelAdded) - neighborhoodData(centerH+2, centerW-1, channelAdded)) + ...
                     abs(neighborhoodData(centerH-4, centerW+2, channelDeal) - neighborhoodData(centerH, centerW, channelDeal)));

        In(8) = (abs(neighborhoodData(centerH-1, centerW+2, channelAdded) - neighborhoodData(centerH+1, centerW-2, channelAdded)) + ...
                     abs(neighborhoodData(centerH-2, centerW+4, channelDeal) - neighborhoodData(centerH, centerW, channelDeal)));

        In(9) = (abs(neighborhoodData(centerH+1, centerW+2, channelAdded) - neighborhoodData(centerH-1, centerW-2, channelAdded)) + ...
                     abs(neighborhoodData(centerH+2, centerW+4, channelDeal) - neighborhoodData(centerH, centerW, channelDeal))); 

        In(10) = (abs(neighborhoodData(centerH+2, centerW+1, channelAdded) - neighborhoodData(centerH-2, centerW-1, channelAdded)) + ...
                      abs(neighborhoodData(centerH+4, centerW+2, channelDeal) - neighborhoodData(centerH, centerW, channelDeal)));

        In(11) = (abs(neighborhoodData(centerH+2, centerW-1, channelAdded) - neighborhoodData(centerH-2, centerW+1, channelAdded)) + ...
                      abs(neighborhoodData(centerH+4, centerW-2, channelDeal) - neighborhoodData(centerH, centerW, channelDeal)));

        In(12) = (abs(neighborhoodData(centerH+1, centerW-2, channelAdded) - neighborhoodData(centerH-1, centerW+2, channelAdded)) + ...
                      abs(neighborhoodData(centerH+2, centerW-4, channelDeal) - neighborhoodData(centerH, centerW, channelDeal)));  
        for i =1 :12
                sumIn = sumIn + (1/(1+In(i)));
        end
        for i =1 :12
            Wn(i) = (1/(1+In(i)))/sumIn;
        end
    case 4
        In(1) = abs(neighborhoodData(centerH-1, centerW-1) - neighborhoodData(centerH+1, centerW+1)) + ...
                abs(neighborhoodData(centerH-2, centerW-2) - neighborhoodData(centerH, centerW));

        In(2) = abs(neighborhoodData(centerH-1, centerW+1) - neighborhoodData(centerH+1, centerW-1)) + ...
                abs(neighborhoodData(centerH-2, centerW+2) - neighborhoodData(centerH, centerW));

        In(3) = abs(neighborhoodData(centerH+1, centerW+1) - neighborhoodData(centerH-1, centerW-1)) + ...
                abs(neighborhoodData(centerH+2, centerW+2) - neighborhoodData(centerH, centerW));

        In(4) = abs(neighborhoodData(centerH+1, centerW-1) - neighborhoodData(centerH-1, centerW+1)) + ...
                abs(neighborhoodData(centerH+2, centerW-2) - neighborhoodData(centerH, centerW));

        for i =1 :4
            sumIn = sumIn + (1/(1+In(i)));
        end
        for i =1 :4
            Wn(i) = (1/(1+In(i)))/sumIn;
        end
end   

$K_{b,n}$或$K_{r,n}$计算：

function Kn = DW_Kn(neighborhoodData, directionNum, channelAdded)
% DW_Kn.m    get rawData from HiRawImage
%   Input:
%       neighborhoodData    the data of neighborhood range
%       directionNum        the num of direction
%       channelAdded        the channal to be interpolated
%   Output:
%       Kn                  The color defference of each direction
% Note: 

% The default value of channalAdded is 1 if there is no input of it
if nargin < 3
    channelAdded = 1;
end
Kn = zeros(directionNum, 1);
[h, w, c] = size(neighborhoodData);
centerH = round(h/2);
centerW = round(w/2);

if c == 1
    % Take the average of the two adjacent samples of the desired color in the
    % Bayer array if the neighborhoodData is Raw
    Kn(1) = neighborhoodData(centerH, centerW-1) - (neighborhoodData(centerH, centerW) + neighborhoodData(centerH, centerW-2)) / 2;
    Kn(2) = neighborhoodData(centerH-1, centerW) - (neighborhoodData(centerH, centerW) + neighborhoodData(centerH-2, centerW)) / 2;
    Kn(3) = neighborhoodData(centerH, centerW+1) - (neighborhoodData(centerH, centerW) + neighborhoodData(centerH, centerW+2)) / 2;
    Kn(4) = neighborhoodData(centerH+1, centerW) - (neighborhoodData(centerH, centerW) + neighborhoodData(centerH+2, centerW)) / 2;

    Kn(5) = neighborhoodData(centerH-1, centerW-2) - (neighborhoodData(centerH-2, centerW-2) + neighborhoodData(centerH, centerW-2)) / 2;
    Kn(6) = neighborhoodData(centerH-2, centerW-1) - (neighborhoodData(centerH-2, centerW-2) + neighborhoodData(centerH-2, centerW)) / 2;
    Kn(7) = neighborhoodData(centerH-2, centerW+1) - (neighborhoodData(centerH-2, centerW) + neighborhoodData(centerH-2, centerW+2)) / 2;
    Kn(8) = neighborhoodData(centerH-1, centerW+2) - (neighborhoodData(centerH-2, centerW+2) + neighborhoodData(centerH, centerW+2)) / 2;
    Kn(9) = neighborhoodData(centerH+1, centerW+2) - (neighborhoodData(centerH+2, centerW+2) + neighborhoodData(centerH, centerW+2)) / 2;
    Kn(10) = neighborhoodData(centerH+2, centerW+1) - (neighborhoodData(centerH+2, centerW+2) + neighborhoodData(centerH+2, centerW)) / 2;
    Kn(11) = neighborhoodData(centerH+2, centerW-1) - (neighborhoodData(centerH+2, centerW-2) + neighborhoodData(centerH+2, centerW)) / 2;
    Kn(12) = neighborhoodData(centerH+1, centerW-2) - (neighborhoodData(centerH+2, centerW-2) + neighborhoodData(centerH, centerW-2)) / 2;
else
    switch directionNum
    case 12
        Kn(1) = neighborhoodData(centerH, centerW-1, 2) - neighborhoodData(centerH, centerW-1, channelAdded);
        Kn(2) = neighborhoodData(centerH-1, centerW, 2) - neighborhoodData(centerH-1, centerW, channelAdded );
        Kn(3) = neighborhoodData(centerH, centerW+1, 2) - neighborhoodData(centerH, centerW+1, channelAdded);
        Kn(4) = neighborhoodData(centerH+1, centerW, 2) - neighborhoodData(centerH+1, centerW, channelAdded);
        
        Kn(5) = neighborhoodData(centerH-1, centerW-2, 2) - neighborhoodData(centerH-1, centerW-2, channelAdded);
        Kn(6) = neighborhoodData(centerH-2, centerW-1, 2) - neighborhoodData(centerH-2, centerW-1, channelAdded);
        Kn(7) = neighborhoodData(centerH-2, centerW+1, 2) - neighborhoodData(centerH-2, centerW+1, channelAdded);
        Kn(8) = neighborhoodData(centerH-1, centerW+2, 2) - neighborhoodData(centerH-1, centerW+2, channelAdded);
        Kn(9) = neighborhoodData(centerH+1, centerW+2, 2) - neighborhoodData(centerH+1, centerW+2, channelAdded);
        Kn(10) = neighborhoodData(centerH+2, centerW+1, 2) - neighborhoodData(centerH+2, centerW+1, channelAdded);
        Kn(11) = neighborhoodData(centerH+2, centerW-1, 2) - neighborhoodData(centerH+2, centerW-1, channelAdded);
        Kn(12) = neighborhoodData(centerH+1, centerW-2, 2) - neighborhoodData(centerH+1, centerW-2, channelAdded);
    case 4
        Kn(1) = neighborhoodData(centerH-1, centerW-1, 2) - neighborhoodData(centerH-1, centerW-1, channelAdded);
        Kn(2) = neighborhoodData(centerH-1, centerW+1, 2) - neighborhoodData(centerH-1, centerW+1, channelAdded);
        Kn(3) = neighborhoodData(centerH+1, centerW+1, 2) - neighborhoodData(centerH+1, centerW+1, channelAdded);
        Kn(4) = neighborhoodData(centerH+1, centerW-1, 2) - neighborhoodData(centerH+1, centerW-1, channelAdded);
    end
end

对于Kn，应该是保留梯度较小的方向的高权重，抑制梯度较大的方向。这里做了个小修改：

        [minValues, minIndices] = mink(In, k);
        sumIn = 0;
        for i =1 :12
                isPresent = ismember(i, minIndices);
            if  isPresent %是最小值
                sumIn = sumIn + (1/(1+0.5*In(i)));
            else
                sumIn = sumIn + (1/(1+1.25*In(i)));
            end 
        end
        for i =1 :12
            isPresent = ismember(i, minIndices);
            if  isPresent %是最小值 高权重
                 Wn(i) = (1/(1+0.5*In(i)))/sumIn;
            else %梯度大，边缘 低权重
                Wn(i) = (1/(1+1.25*In(i)))/sumIn;
            end 
            
        end

效果：和以上没进行边缘检测的方法相比，质量是显著提升的。

基于边缘检测的方法

1）先计算水平梯度和垂直梯度
$$e_h=(|C_{i-1,j-1}-C_{i-1,j+1}|+2|C_{i,j-1}-C_{i,j+1}|+|C_{i+1,j-1}-C_{i+1,j+1}|)/4;$$
$$e_v=(|C_{i-1,j-1}-C_{i+1,j-1}|+2|C_{i-1,j}-C_{i+1,j}|+|C_{i-1,j+1}-C_{i+1,j+1}|)/4;$$

2）计算两个方向上的权重
$$W_h=\frac{1}{e_h}/\frac{1}{e_h+e_v}$$
$$W_v=\frac{1}{e_v}/\frac{1}{e_h+e_v}$$

计算缺失的G

最后

计算缺失的RB值/计算缺失的G值

跟色差法同理

每天补一补坑。