【去马赛克专题】demosaic算法之边缘判别插值

浩瀚之水_csdn

已于 2024-07-01 11:18:33 修改

阅读量310

点赞数

分类专栏：图像处理ISP专题文章标签：算法人工智能计算机视觉 ISP

于 2024-07-01 10:48:33 首次发布

原文链接：https://zhuanlan.zhihu.com/p/686883055

版权

图像处理ISP专题专栏收录该内容

35 篇文章 31 订阅

订阅专栏

1. 边缘判别插值简介

在bayer demosaic时，G通道的插值中边缘产生的错误导致最终估计的彩色图像出现伪彩色和锯齿。无论是之前提到的色差法还是色比法，由于G通道的插值通常都作为demosaic算法的第一步，因此G通道插值的准确性对后续R、B通道的重建精度有着至关重要的影响。

双线性插值中，G通道的锯齿导致图像拉链现象

从bayer到rgb：ISP中的demosaic技术

demosaic算法之基础插值 - 知乎 (zhihu.com)

边缘判别方法在基础插值方法的基础上不仅考虑了空间相关性和光谱相关性，还考虑了边缘效应。边缘判别法通过指定策略判断边缘的方向去控制G通道插值时的邻域采样选择使得插值沿着边缘方向进行，避免了边缘两侧成分的混叠。主要的插值策略依靠双线性插值和色差恒定法，对于伪彩色、拉链效应有着明显的改善。

2. 边缘判别插值算法详解

2.1 梯度矫正线性插值法

参考论文：Malvar-He-Cutler Linear Image Demosaicking

基础的色差和色比方法将G通道看作亮度成分， R、B通道看成色度成分。由于G通道的采样是R、B通道的两倍，所以G通道能代表图像中大部分高频信息。R、B通道作为色度成分变化比较缓慢，因此插值后的G分量图像能够更准确地反映图像中的轮廓。通常先在双线性插值的基础上对G通道进行重建，然后再利用色差或者色比恒定的假设重建R、B通道。

以一个R像素位置为例（图中标记+号的位置），可以看到当计算这个位置G像素的插值时，色差和色比法会丢弃这个位置的R值仅对邻域的G值做平均，但很明显它是有价值的信息。在梯度矫正线性插值法中，引入了R这一点的观测信息控制边缘处的叠加分量以进行梯度修正，将R值与R邻域的双线性插值估计进行比较。如果它与估计值不同，这可能意味着在这个像素处有一个剧烈的亮度变化（即处于边缘位置）。

具体来说，作者提出了一组矫正计算公式，假设在一个像素点上已知A通道需要插值C通道，那么插值结果为：

可以看到C通道的插值结果不仅由C通道邻域的加权平均表示，还由A通道在该点的边缘响应的增益分量控制，该分量为A通道与水平和垂直方向最近的4个A通道的差乘一个增益系数表示。类似于对于边缘部分的插值结果进行了一个锐化的过程。

设定三通道对应的增益系数为α、β、r，作者通过kodak数据集进行大量测试得到一组经验性系数α=0.5、β=0.625、r=0.75，最终表示为各通道之间的5X5插值滤波核如下：

目前大家熟知的图像画质评测软件Imatest中，demosaic算法就是用的这篇论文。

2.2 Adam-Hamilton梯度自适应插值

参考论文：Self-similarity driven color demosaicking

从2.2开始会介绍严格意义上通过边缘方向进行插值的方法，在介绍亚当-汉密尔顿插值算法前，首先介绍一阶微分和二阶微分的边缘插值。

一阶微分边缘敏感性插值算法通过亮度分量(G通道)确定各像素位置可能的边缘方向，通过计算出水平和竖直方向的梯度判断插值方向，具体流程如下：

(1) 在G通道缺失值的像素位置（i，j）处用G通道邻域像素计算水平和垂直方向梯度：

(2) 定义阈值T，表示水平和垂直梯度分量的均值，对G通道进行插值：

选择梯度小的方向进行插值

(3) 通过色差法重建R、B通道。

二阶微分的边缘插值基于色度分量确定边缘，此处使用拉普拉斯算子。图像经过二阶微分后在边缘产生一个陡峭的零交叉，根据这个零交叉可以来判断边缘。

对于边缘处由于一阶微分是一个极值点，二阶微分是零点，确定零点的位置比极值点更容易也更准确。具体流程如下：

(1) 以R通道中缺失的G值重建为例，在位置（i，j）处计算水平和垂直方向梯度，与一阶微分不同的是该算法利用了缺失像素位置的观测通道（异色的R通道）计算梯度：

(2) 定义阈值T为水平和垂直梯度分量的均值，插值G通道时当垂直梯度小就按竖直方向插值，如果梯度差不多就说明是平坦区域直接双线性插值：

(3) 通过色差法重建R、B通道。

亚当-汉密尔顿插值算法为了更准确判断边缘方向，同时考虑了同色通道的一阶梯度和异色通道的二阶差分，利用亮度和色彩的相似性去确定G通道插值的方向是选择水平方向像素加权还是垂直方向像素加权。

另外在色差法重建R、B通道时，对于R通道上的B缺失值（B通道上的R缺失值）的计算时，额外考虑到了对角邻域的方向判断，去确定按照45°方向像素加权还是135°方向的像素进行加权。该算法可以明显改善伪彩色和拉链效应，但缺陷是对于密集纹理的方向判断效果不佳，其次插值的时候本质是低通滤波器，那么插值得到的色彩的高频分量可能不准确。具体计算流程如下：

(1) 在G通道缺失值的位置（i，j）处计算水平和垂直方向梯度：

(2) 根据水平和垂直方向的大小关系判断插值方向，并对G通道进行插值：

可以看到G通道插值引入了R通道的信息，这就类似2.1提到的做法，也可以理解成是色差法的一种转换

(3) 通过色差法重建R和B通道，以R通道的重建为例，对于已知G通道值而缺失R值的位置插值，由于该位置水平或垂直方向邻域存在观测到的R像素，插值公式为：

而对于已知B通道缺失R通道值的像素，该位置只在对角邻域存在观测R值，此时引入额外判断对比45°和135°对角方向的梯度大小:

插值计算逻辑则为：

总结下来亚当-汉密尔顿算法以同色一阶微分加异色二阶差分算梯度，而且考虑了4个方向的梯度计算。插值时也会考虑通道相关性用同色和异色的像素进行加权，第三节的仿真可以发现高频区域有很明显的效果提升。

2.3 AHD (Adaptive homogeneity-directed)

参考论文：Adaptive homogeneity-directed demosaicing algorithm

采用边缘判别插值方法的demosaic算法存在两种类型的颜色伪影：第一种被称为误导色伪影，当插值方向选择错误时，就会产生误导。第二种类型的颜色伪影与插值中的限制有关：也就是说，即使有一个完美的方向选择器，插值由于频繁地从一个方向切换到另一个方向也会导致输出图像的不连续。

相较亚当-汉密尔顿插值法，AHD主要改进为提高边缘方向判断的精度和平滑度：算法根据图像的颜色同质性（color homogeneity）来自适应的在彩色伪影较少的方向上插值，图像的不同区域可能存在不同程度的颜色变化，算法分析相邻像素间的颜色差异和空间相关性进而决定最佳的插值方向，以解决颜色伪影和混叠。

在CIELAB空间分别统计邻域的距离B、亮度差L、色度差C进行综合评价得到homogeneity map：

d是欧几里得距离

同质图的计算思路是在中心的一个邻域（集合B）内观察所有像素与中心像素的相似度，因此定义一个相似度函数H。一个像素与中心像素“相似“的标准是，CIELAB空间上的亮度误差小于一个阈值，色度误差也小于一个阈值，符合这个条件的待考察像素定义为相似，否则不相似。H函数的值定义为相似的数量与全部像素数量之比，该值越接近1，伪彩越轻微。

具体计算流程如下：

(1) 首先分别对三通道分别只进行水平和垂直单个方向的插值，得到Yh和Yv（跟亚当汉密尔顿的插值公式一样）。

以G通道为例，水平和垂直的插值公式和2.2相同

(2) 分别计算水平和垂直方向的同质图Hx和Hy，homo map两个方向中同质值大的方向就对应选择Yh和Yv中方向插值结果。

(3) 考虑到色差空间变化比较平缓，在邻域使用中值滤波对artifact进行m次的色差迭代修正：

AHD主体思想还是基于同一物体颜色相似的原理：如果插值方向是正确的，那么沿着正确边缘方向上的像素应属于同一物体，颜色的梯度应较小。如果插值方向是错误的，则在错误边缘上的像素将属于不同物体，颜色的梯度应较大。这种先插值再校验的方法称为“后验”方法，已经得到了大量的研究和应用。

2.4 GBTF (Gradient Based Threshold Free）

参考论文：Gradient based threshold free color filter array interpolation

之前提到的方法大多根据梯度信息判断插值方向，GBTF一个优化是在色差空间进行梯度计算进行改进，这是因为G通道上包含了大量高频信息，错误的方向判断会导致边缘周围的像素混叠，而色差空间是相较平缓的，加权平均对边缘的影响较小。

另外的一个改进是对水平垂直两个方向解耦构建了N、S、W、E四个插值方向，通过各方向5X5邻域的梯度权重之和来判断各插值方向的权重。最后基于色差法对R和B通道进行重建。

GBTF具体的算法流程如下：

构建色差空间的梯度

(1) 分别对G通道进行水平和垂直方向的插值，得到Gh和Gv，插值公式与亚当-汉密尔顿法相同。

(2) 分别对R和B通道进行水平和垂直方向的插值，得到Rh、Bh和Rv、Gv，插值公式与亚当-汉密尔顿法相同。

(3) 计算水平方向和垂直方向的色差Kh和Kv。

(4)计算色差空间中水平和垂直方向的梯度gradh、gradv。

N、S、W、E四方向插值权重计算

如图所示，针对红色圆圈标记的G通道缺失值，在色差空间得到各点的梯度后，分别计算四个方向5X5区域的梯度和：

(1) N、S方向代表垂直方向，W、E代表水平方向，分别在色差空间计算gradv在N、S方向ROI区域的梯度和与gradh在W、E方向ROI区域的梯度和得到四个方向的权重weight。

(2) 色差梯度大的地方考虑为边缘部分，希望插值时避免沿着这个方向进行，因此最终取权重的倒数1/weight作为各方向的插值权重。

色差法对G通道插值

通过色差法对G通道缺失值进行加权插值重建，权重为上面得到的各方向权重。需要注意的是各方向的插值采样点选择从插值点开始朝该方向共计5个色差采样点，以R像素位置的G值重建为例：

在色差空间中，中心点的色差G-R等于四个方向各5个点的色差和与该方向权重加权的结果。

R和B通道插值

对于Gr、Gb处通过色差法重建R和B通道分量，对于R通道处B通道（B通道处R通道）分量会在色差空间进行7X7邻域的稀疏采样加权(一定程度的通道稀疏来抑制噪声)，然后再利用色差法得到对应的R和B值。

3. 代码仿真

一阶微分法

def derivative1(img,pattern):
    img = img.astype(np.float64)
    R_m, G_m, B_m = masks_Bayer(img, pattern)
    R_m = img * R_m
    G_m = img * G_m
    B_m = img * B_m
    h, w = img.shape
    R = np.zeros((h,w))
    G = np.zeros((h,w))
    B = np.zeros((h,w))

    for i in range(2,h-2):
        for j in range(2,w-2):
            gradv = np.abs(G_m[i - 1, j] - G_m[i + 1, j])
            gradh = np.abs(G_m[i,j-1]-G_m[i,j+1])
            T=(gradv+gradh)/2
            if ((i % 2 == 0) & (j % 2 == 0)):
                if((gradh>gradv)):
                #if ((gradh > T) & (gradv < T)):
                    G[i, j] = (G_m[i - 1, j] + G_m[i + 1, j])/2
              #  elif((gradv>T) & (gradh<T)):
                elif ((gradh < gradv)):
                    G[i, j] = (G_m[i, j-1] + G_m[i, j+1]) / 2
                else:
                    G[i, j] = (G_m[i - 1, j] + G_m[i + 1, j] + G_m[i, j - 1] + G_m[i, j + 1]) / 4
            elif ((i % 2 == 1) & (j % 2 == 1)):
                if((gradh>gradv)):
                    G[i, j] = (G_m[i - 1, j] + G_m[i + 1, j])/2
                elif ((gradh < gradv)):
                    G[i, j] = (G_m[i, j-1] + G_m[i, j+1]) / 2
                else:
                    G[i, j] = (G_m[i - 1, j] + G_m[i + 1, j] + G_m[i, j - 1] + G_m[i, j + 1]) / 4
            elif ((i % 2 == 0) & (j % 2 == 1)):
                G[i,j] = G_m[i,j]
            else:
                G[i,j] = G_m[i,j]

    Kr = np.zeros((h, w))
    Kb = np.zeros((h, w))
    for i in range(2,h-2):
        for j in range(2, w - 2):
            if ((i % 2 == 0) & (j % 2 == 0)):
                Kr[i, j] = G[i, j] - R_m[i, j]
                Kb[i, j] = G[i, j] - (B_m[i - 1, j - 1] + B_m[i - 1, j + 1] + B_m[i + 1, j - 1] + B_m[
                        i + 1, j + 1]) / 4
            if ((i % 2 == 1) & (j % 2 == 1)):
                Kb[i, j] = G[i, j] - B_m[i, j]
                Kr[i, j] = G[i, j] - (R_m[i - 1, j - 1] + R_m[i - 1, j + 1] + R_m[i + 1, j - 1] + R_m[
                        i + 1, j + 1]) / 4
            if ((i % 2 == 0) & (j % 2 == 1)):
                Kr[i, j] = G[i, j] - (R_m[i, j - 1] + R_m[i, j + 1]) / 2
                Kb[i, j] = G[i, j] - (B_m[i - 1, j] + B_m[i + 1, j]) / 2
            if ((i % 2 == 1) & (j % 2 == 0)):
                Kr[i, j] = G[i, j] - (R_m[i - 1, j] + R_m[i + 1, j]) / 2
                Kb[i, j] = G[i, j] - (B_m[i, j - 1] + B_m[i, j + 1]) / 2

    for i in range(2,h-2):
        for j in range(2,w-2):
            if ((i % 2 == 0) & (j % 2 == 0)):
                R[i, j] = R_m[i, j]
                G[i,j] = R_m[i,j]+(Kr[i-1,j]+Kr[i+1,j]+Kr[i,j-1]+Kr[i,j+1])/4
                B[i, j] = G[i, j] - (Kb[i-1,j]+Kb[i+1,j]+Kb[i,j-1]+Kb[i,j+1])/4
            if ((i % 2 == 1) & (j % 2 == 1)):
                B[i, j] = B_m[i, j]
                G[i,j] = B_m[i,j]+(Kb[i-1,j]+Kb[i+1,j]+Kb[i,j-1]+Kb[i,j+1])/4
                R[i, j] = G[i, j] - (Kr[i-1,j]+Kr[i+1,j]+Kr[i,j-1]+Kr[i,j+1])/4
            if ((i % 2 == 0) & (j % 2 == 1)):
                G[i, j] = G_m[i, j]
                R[i, j] = G[i, j] - (Kr[i - 1, j] + Kr[i + 1, j] + Kr[i, j - 1] + Kr[i, j + 1]) / 4
                B[i, j] = G[i, j] - (Kb[i - 1, j] + Kb[i + 1, j] + Kb[i, j - 1] + Kb[i, j + 1]) / 4
            if ((i % 2 == 1) & (j % 2 == 0)):
                G[i, j] = G_m[i, j]
                R[i, j] = G[i, j] - (Kr[i - 1, j] + Kr[i + 1, j] + Kr[i, j - 1] + Kr[i, j + 1]) / 4
                B[i, j] = G[i, j] - (Kb[i - 1, j] + Kb[i + 1, j] + Kb[i, j - 1] + Kb[i, j + 1]) / 4
    result_img = np.zeros((h, w, 3))
    result_img[:, :, 0] = R
    result_img[:, :, 1] = G
    result_img[:, :, 2] = B
    del R_m, G_m, B_m
    return result_img

二阶微分法

def derivative2(img,pattern):
    img = img.astype(np.float64)
    R_m, G_m, B_m = masks_Bayer(img, pattern)
    R_m = img * R_m
    G_m = img * G_m
    B_m = img * B_m
    h, w = img.shape
    R = np.zeros((h,w))
    G = np.zeros((h,w))
    B = np.zeros((h,w))

    for i in range(3,h-3):
        for j in range(3,w-3):
            gradvr = np.abs(R_m[i, j] - (R_m[i-2,j]+R_m[i+2,j])/2)
            gradhr = np.abs(R_m[i, j] - (R_m[i,j-2]+R_m[i,j+2])/2)
            gradvb = np.abs(B_m[i, j] - (B_m[i-2,j]+B_m[i+2,j])/2)
            gradhb = np.abs(B_m[i, j] - (B_m[i,j-2]+B_m[i,j+2])/2)
            Tr=(gradvr+gradhr)/2
            Tb = (gradvb + gradhb) / 2
            if ((i % 2 == 0) & (j % 2 == 0)):
                if((gradhr>Tr) & (gradvr<Tr)):
                    G[i, j] = (G_m[i - 1, j] + G_m[i + 1, j])/2
                elif((gradvr>Tr) & (gradhr<Tr)):
                    G[i, j] = (G_m[i, j-1] + G_m[i, j+1]) / 2
                else:
                    G[i, j] = (G_m[i - 1, j] + G_m[i + 1, j] + G_m[i, j - 1] + G_m[i, j + 1]) / 4
            elif ((i % 2 == 1) & (j % 2 == 1)):
                if((gradhb>Tb) & (gradvb<Tb)):
                    G[i, j] = (G_m[i - 1, j] + G_m[i + 1, j])/2
                elif((gradvb>Tb) & (gradhb<Tb)):
                    G[i, j] = (G_m[i, j-1] + G_m[i, j+1]) / 2
                else:
                    G[i, j] = (G_m[i - 1, j] + G_m[i + 1, j] + G_m[i, j - 1] + G_m[i, j + 1]) / 4
            elif ((i % 2 == 0) & (j % 2 == 1)):
                G[i,j] = G_m[i,j]
            else:
                G[i,j] = G_m[i,j]

    Kr = np.zeros((h, w))
    Kb = np.zeros((h, w))
    for i in range(2,h-2):
        for j in range(2, w - 2):
            if ((i % 2 == 0) & (j % 2 == 0)):
                Kr[i, j] = G[i, j] - R_m[i, j]
                Kb[i, j] = G[i, j] - (B_m[i - 1, j - 1] + B_m[i - 1, j + 1] + B_m[i + 1, j - 1] + B_m[
                        i + 1, j + 1]) / 4
            if ((i % 2 == 1) & (j % 2 == 1)):
                Kb[i, j] = G[i, j] - B_m[i, j]
                Kr[i, j] = G[i, j] - (R_m[i - 1, j - 1] + R_m[i - 1, j + 1] + R_m[i + 1, j - 1] + R_m[
                        i + 1, j + 1]) / 4
            if ((i % 2 == 0) & (j % 2 == 1)):
                Kr[i, j] = G[i, j] - (R_m[i, j - 1] + R_m[i, j + 1]) / 2
                Kb[i, j] = G[i, j] - (B_m[i - 1, j] + B_m[i + 1, j]) / 2
            if ((i % 2 == 1) & (j % 2 == 0)):
                Kr[i, j] = G[i, j] - (R_m[i - 1, j] + R_m[i + 1, j]) / 2
                Kb[i, j] = G[i, j] - (B_m[i, j - 1] + B_m[i, j + 1]) / 2

    for i in range(2,h-2):
        for j in range(2,w-2):
            if ((i % 2 == 0) & (j % 2 == 0)):
                R[i, j] = R_m[i, j]
                G[i,j] = R_m[i,j]+(Kr[i-1,j]+Kr[i+1,j]+Kr[i,j-1]+Kr[i,j+1])/4
                B[i, j] = G[i, j] - (Kb[i-1,j]+Kb[i+1,j]+Kb[i,j-1]+Kb[i,j+1])/4
            if ((i % 2 == 1) & (j % 2 == 1)):
                B[i, j] = B_m[i, j]
                G[i,j] = B_m[i,j]+(Kb[i-1,j]+Kb[i+1,j]+Kb[i,j-1]+Kb[i,j+1])/4
                R[i, j] = G[i, j] - (Kr[i-1,j]+Kr[i+1,j]+Kr[i,j-1]+Kr[i,j+1])/4
            if ((i % 2 == 0) & (j % 2 == 1)):
                G[i, j] = G_m[i, j]
                R[i, j] = G[i, j] - (Kr[i - 1, j] + Kr[i + 1, j] + Kr[i, j - 1] + Kr[i, j + 1]) / 4
                B[i, j] = G[i, j] - (Kb[i - 1, j] + Kb[i + 1, j] + Kb[i, j - 1] + Kb[i, j + 1]) / 4
            if ((i % 2 == 1) & (j % 2 == 0)):
                G[i, j] = G_m[i, j]
                R[i, j] = G[i, j] - (Kr[i - 1, j] + Kr[i + 1, j] + Kr[i, j - 1] + Kr[i, j + 1]) / 4
                B[i, j] = G[i, j] - (Kb[i - 1, j] + Kb[i + 1, j] + Kb[i, j - 1] + Kb[i, j + 1]) / 4
    result_img = np.zeros((h, w, 3))
    result_img[:, :, 0] = 0
    result_img[:, :, 1] = 0
    result_img[:, :, 2] = B
    del R_m, G_m, B_m
    return result_img

梯度矫正线性插值法

def highquality(img,pattern):
    img = img.astype(np.float64)
    R_m, G_m, B_m = masks_Bayer(img, pattern)
    R_m = img * R_m
    G_m = img * G_m
    B_m = img * B_m
    h, w = img.shape
    R = np.zeros((h,w))
    G = np.zeros((h,w))
    B = np.zeros((h,w))
    G_RB = np.array([[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
    RB_BR = np.array([[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
    RB_Gr = np.array([[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
    RB_Gb = np.array([[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]])/8

    for i in range(3,h-3):
        for j in range(3,w-3):
            if ((i % 2 == 0) & (j % 2 == 0)):
                G[i,j] = np.sum(img[i-2:i+3,j-2:j+3]*G_RB)
                B[i,j] = np.sum(img[i-2:i+3,j-2:j+3]*RB_BR)
                R[i,j] = R_m[i,j]
            elif ((i % 2 == 1) & (j % 2 == 1)):
                G[i,j] =np.sum(img[i-2:i+3,j-2:j+3]*G_RB)
                R[i,j] = np.sum(img[i-2:i+3,j-2:j+3]*RB_BR)
                B[i,j] = B_m[i,j]
            elif ((i % 2 == 0) & (j % 2 == 1)):
                R[i,j] = np.sum(img[i-2:i+3,j-2:j+3]*RB_Gr)
                B[i,j] = np.sum(img[i-2:i+3,j-2:j+3]*RB_Gb)
                G[i,j] = G_m[i,j]
            else:
                R[i,j] = np.sum(img[i-2:i+3,j-2:j+3]*RB_Gb)
                B[i,j] = np.sum(img[i-2:i+3,j-2:j+3]*RB_Gr)
                G[i,j] = G_m[i,j]

    result_img = np.zeros((h, w, 3))

    result_img[:, :, 0] = R
    result_img[:, :, 1] = G
    result_img[:, :, 2] = B

    del R_m, G_m, B_m
    return result_img

Adam-Hamilton梯度自适应插值

def AH_interpolate(img, pattern, gamma, max_value):
    X = img
    Rm, Gm, Bm = masks_Bayer(img, pattern)

    # green
    Hg1 = np.array([   0, 1/2,   0, 1/2,    0])
    Hg2 = np.array([-1/4,   0, 1/2,   0, -1/4])
    Hg = Hg1 + Hg2 * gamma
    Hg = Hg.reshape(1, -1)
    G = Gm * X + (Rm + Bm) * signal.convolve(X, Hg, 'same')

    # red/blue
    Hr = [[1/4, 1/2, 1/4],
          [1/2,   1, 1/2],
          [1/4, 1/2, 1/4]]
    R = G + signal.convolve(Rm*(X-G), Hr, 'same')
    B = G + signal.convolve(Bm*(X-G), Hr, 'same')

    R = np.clip(R, 0, max_value)
    G = np.clip(G, 0, max_value)
    B = np.clip(B, 0, max_value)

    return R, G, B
def AH_interpolateX(img, pattern, gamma, max_value):
    h, w = img.shape
    Y = np.zeros((h, w, 3))

    R, G, B = AH_interpolate(img, pattern, gamma, max_value)

    Y[:, :, 0] = R
    Y[:, :, 1] = G
    Y[:, :, 2] = B

    return Y
def AH_interpolateY(img, pattern, gamma, max_value):
    if (pattern == "RGGB"):
        new_pattern="RGGB"
    elif (pattern == "GRBG"):
        new_pattern = "GBRG"
    elif (pattern == "GBRG"):
        new_pattern = "GRBG"
    elif (pattern == "BGGR"):
        new_pattern = "BGGR"
    else:
        print("pattern must be one of :  RGGB, GBRG, GBRG, or BGGR")
        return

    h, w = img.shape
    Y = np.zeros((h, w, 3))

    new_img = img.T
    R, G, B = AH_interpolate(new_img, new_pattern, gamma, max_value)

    Y[:, :, 0] = R.T
    Y[:, :, 1] = G.T
    Y[:, :, 2] = B.T

    return Y
def AH_gradient(img, pattern):
    X = img
    Rm, Gm, Bm = masks_Bayer(img, pattern)

    # green
    Hg1 = np.array([0,1,0,-1,0])
    Hg2 = np.array([-1,0,2,0,-1])

    Hg1 = Hg1.reshape(1,-1)
    Hg2 = Hg2.reshape(1,-1)
    Ga = (Rm + Bm) * (np.abs(signal.convolve(X, Hg1, 'same')) + np.abs(signal.convolve(X, Hg2, 'same')))

    return Ga
def AH_gradientX(img, pattern):
    Ga = AH_gradient(img, pattern)
    return Ga
def AH_gradientY(img,pattern):
    if (pattern == "RGGB"):
        new_pattern = "RGGB"
    elif (pattern == "GRBG"):
        new_pattern = "GBRG"
    elif (pattern == "GBRG"):
        new_pattern = "GRBG"
    elif (pattern == "BGGR"):
        new_pattern = "BGGR"
    else:
        print("pattern must be one of :  RGGB, GBRG, GBRG, or BGGR")
        return

    new_img = img.T
    Ga = AH_gradient(new_img, new_pattern)
    new_Ga = Ga.T

    return new_Ga
def AH(img, pattern):
    img = img.astype(np.float64)
    R_m, G_m, B_m = masks_Bayer(img, pattern)
    R_m = img * R_m
    G_m = img * G_m
    B_m = img * B_m
    h, w = img.shape
    R = np.zeros((h, w))
    G = np.zeros((h, w))
    B = np.zeros((h, w))

    for i in range(2, h - 2):
        for j in range(2, w - 2):
            gradv = np.abs(G_m[i - 1, j] - G_m[i + 1, j])
            gradh = np.abs(G_m[i, j - 1] - G_m[i, j + 1])
            gradvr = np.abs(2*R_m[i, j] - (R_m[i-2,j]+R_m[i+2,j]))
            gradhr = np.abs(2*R_m[i, j] - (R_m[i,j-2]+R_m[i,j+2]))
            gradvb = np.abs(2*B_m[i, j] - (B_m[i-2,j]+B_m[i+2,j]))
            gradhb = np.abs(2*B_m[i, j] - (B_m[i,j-2]+B_m[i,j+2]))
            if ((i % 2 == 0) & (j % 2 == 0)):
                if ((gradh - gradv+gradhr-gradvr)>0):
                    G[i, j] = ((G_m[i - 1, j] + G_m[i + 1, j]) / 2)+(2*R_m[i, j] - (R_m[i-2,j]+R_m[i+2,j]))/4
                elif ((gradh - gradv+gradhr-gradvr)<0):
                    G[i, j] = ((G_m[i, j - 1] + G_m[i, j + 1]) / 2)+(2*R_m[i, j] - (R_m[i,j-2]+R_m[i,j+2]))/4
                else:
                    G[i, j] = ((G_m[i - 1, j] + G_m[i + 1, j] + G_m[i, j - 1] + G_m[i, j + 1]) / 4)+((2*R_m[i, j] - (R_m[i-2,j]+R_m[i+2,j]))+(2*R_m[i, j] - (R_m[i,j-2]+R_m[i,j+2])))/8
            elif ((i % 2 == 1) & (j % 2 == 1)):
                if ((gradh - gradv+gradhb-gradvb)>0):
                    G[i, j] = ((G_m[i - 1, j] + G_m[i + 1, j]) / 2)+(2*B_m[i, j] - (B_m[i-2,j]+B_m[i+2,j]))/4
                elif ((gradh - gradv+gradhr-gradvr)<0):
                    G[i, j] = ((G_m[i, j - 1] + G_m[i, j + 1]) / 2)+(2*B_m[i, j] - (B_m[i,j-2]+B_m[i,j+2]))/4
                else:
                    G[i, j] = ((G_m[i - 1, j] + G_m[i + 1, j] + G_m[i, j - 1] + G_m[i, j + 1]) / 4)+((2*B_m[i, j] - (B_m[i-2,j]+B_m[i+2,j]))+(2*B_m[i, j] - (B_m[i,j-2]+B_m[i,j+2])))/8
            elif ((i % 2 == 0) & (j % 2 == 1)):
                G[i, j] = G_m[i, j]
            else:
                G[i, j] = G_m[i, j]

    Kr = np.zeros((h, w))
    Kb = np.zeros((h, w))
    for i in range(2,h-2):
        for j in range(2,w-2):
            if ((i % 2 == 0) & (j % 2 == 0)):
                Kr[i,j] = G[i,j]-R_m[i,j]
                Kb[i, j] = G[i,j]-(B_m[i-1,j-1]+B_m[i-1,j+1]+B_m[i+1,j-1]+B_m[i+1,j+1])/4
            if ((i % 2 == 1) & (j % 2 == 1)):
                Kb[i,j] = G[i,j]-B_m[i,j]
                Kr[i, j] = G[i,j]-(R_m[i-1,j-1]+R_m[i-1,j+1]+R_m[i+1,j-1]+R_m[i+1,j+1])/4
            if ((i % 2 == 0) & (j % 2 == 1)):
                Kr[i, j] = G[i, j] - (R_m[i, j-1]+R_m[i, j+1])/2
                Kb[i, j] = G[i, j] - (B_m[i-1, j]+B_m[i+1, j])/2
            if ((i % 2 == 1) & (j % 2 == 0)):
                Kr[i, j] = G[i, j] - (R_m[i-1, j]+R_m[i+1, j])/2
                Kb[i, j] = G[i, j] - (B_m[i, j-1]+B_m[i, j+1])/2

    for i in range(2,h-2):
        for j in range(2,w-2):
            grad45r = np.abs(R_m[i - 1, j+1] - R_m[i + 1, j-1])
            grad135r = np.abs(R_m[i-1, j - 1] - R_m[i+1, j + 1])
            grad45b = np.abs(B_m[i - 1, j+1] - B_m[i + 1, j-1])
            grad135b = np.abs(B_m[i-1, j - 1] - B_m[i+1, j + 1])
            grad45g = np.abs(2*G[i, j] - (G[i-1,j+1]+G[i+1,j-1]))
            grad135g = np.abs(2*G[i, j] - (G[i-1,j-1]+G[i+1,j+1]))

            if ((i % 2 == 0) & (j % 2 == 0)):
                R[i, j] = R_m[i, j]
                if((grad45b+grad45g-grad135b-grad135g)>0):
                    B[i, j] =(B_m[i-1, j - 1] + B_m[i+1, j + 1])/2+(2*G[i, j] - (G[i-1,j-1]+G[i+1,j+1]))/4
                elif((grad45b + grad45g - grad135b - grad135g)<0):
                    B[i, j] = (B_m[i - 1, j+1] + B_m[i + 1, j-1]) / 2 + (2*G[i, j] - (G[i-1,j+1]+G[i+1,j-1])) / 4
                else:
                    B[i, j] = (B_m[i-1, j - 1] + B_m[i+1, j + 1]+B_m[i - 1, j + 1] + B_m[i + 1, j - 1]) / 4 + ((2*G[i, j] - (G[i-1,j+1]+G[i+1,j-1]))+(2*G[i, j] - (G[i-1,j-1]+G[i+1,j+1]))) / 8
            if ((i % 2 == 1) & (j % 2 == 1)):
                B[i, j] = B_m[i, j]
                if((grad45r+grad45g-grad135r-grad135g)>0):
                    R[i, j] =(R_m[i-1, j - 1] + R_m[i+1, j + 1])/2+(2*G[i, j] - (G[i-1,j-1]+G[i+1,j+1]))/4
                elif((grad45b + grad45g - grad135b - grad135g)<0):
                    R[i, j] = (R_m[i - 1, j+1] + R_m[i + 1, j-1]) / 2 + (2*G[i, j] - (G[i-1,j+1]+G[i+1,j-1])) / 4
                else:
                    R[i, j] = (R_m[i-1, j - 1] + R_m[i+1, j + 1]+R_m[i - 1, j + 1] + R_m[i + 1, j - 1]) / 4 + ((2*G[i, j] - (G[i-1,j+1]+G[i+1,j-1]))+(2*G[i, j] - (G[i-1,j-1]+G[i+1,j+1]))) / 8

            if ((i % 2 == 0) & (j % 2 == 1)):
                R[i, j] = G[i, j] - ( Kr[i, j - 1] + Kr[i, j + 1]) / 2
                B[i, j] = G[i, j] - (Kb[i - 1, j] + Kb[i + 1, j] ) / 2
            if ((i % 2 == 1) & (j % 2 == 0)):
                R[i, j] = G[i, j] - (Kr[i - 1, j] + Kr[i + 1, j] ) / 2
                B[i, j] = G[i, j] - (Kb[i, j - 1] + Kb[i, j + 1]) / 2
    result_img = np.zeros((h, w, 3))
    #R, G, B = MNartifact(R, G, B, 2)
    result_img[:, :, 0] = R
    result_img[:, :, 1] = G
    result_img[:, :, 2] = B
    del R_m, G_m, B_m
    return result_img
def Adamhamilton(img, pattern):
    imgh, imgw = img.shape
    imgs = 10
    gamma = 1
    max_value = 255
    # X,Y方向插值
    # 扩展大小
    # Y = [X(n + 1:-1: 2, n + 1: -1:2,:)       X(n + 1: -1:2,:,:)        X(n + 1: -1:2, end - 1: -1:end - n,:)
    # X(:, n + 1: -1:2,:)               X                                X(:, end - 1: -1:end - n,:)
    # X(end - 1: -1:end - n, n + 1: -1:2,:)  X(end - 1: -1:end - n,:,:)  X(end - 1: -1:end - n, end - 1: -1:end - n,:)];
    f = np.pad(img, ((imgs, imgs), (imgs, imgs)), 'reflect')

    Yx = AH_interpolateX(f, pattern, gamma, max_value)
    Yy = AH_interpolateY(f, pattern, gamma, max_value)

    Hx = AH_gradientX(f, pattern)
    Hy = AH_gradientY(f, pattern)

    # set output to Yy if Hy <= Hx
    index = np.where(Hy <= Hx)
    R = Yx[:, :, 0]
    G = Yx[:, :, 1]
    B = Yx[:, :, 2]
    Ry = Yy[:, :, 0]
    Gy = Yy[:, :, 1]
    By = Yy[:, :, 2]

    Rs = R
    Gs = G
    Bs = B
    Rs[index] = Ry[index]
    Gs[index] = Gy[index]
    Bs[index] = By[index]

    h, w = Rs.shape
    Y = np.zeros((h, w, 3))
    Y[:, :, 0] = Rs
    Y[:, :, 1] = Gs
    Y[:, :, 2] = Bs
    # 调整size和值的范畴
    Y = np.clip(Y, 0, max_value)
    resultY = Y[imgs:imgs + imgh, imgs:imgs + imgw, :]
    return resultY

AHD

def MNparamA(YxLAB, YyLAB):
    # epsilon = min(max(left, right), max(top, bottm))
    X = YxLAB
    Y = YyLAB

    kernel_H1 = np.array([1, -1, 0])
    kernel_H1 = kernel_H1.reshape(1, -1)
    kernel_H2 = np.array([0, -1, 1])
    kernel_H2 = kernel_H2.reshape(1, -1)

    kernel_V1 = kernel_H1.reshape(1, -1).T
    kernel_V2 = kernel_H2.reshape(1, -1).T

    eLM1 = np.maximum(np.abs(signal.convolve(X[:,:, 0], kernel_H1, 'same')), np.abs(signal.convolve(X[:,:, 0], kernel_H2, 'same')))
    eLM2 = np.maximum(np.abs(signal.convolve(Y[:,:, 0], kernel_V1, 'same')), np.abs(signal.convolve(Y[:,:, 0], kernel_V2, 'same')))
    eL = np.minimum(eLM1, eLM2)

    eCx = np.maximum(signal.convolve(X[:,:, 1], kernel_H1, 'same')**2 + signal.convolve(X[:,:, 2], kernel_H1, 'same')**2, signal.convolve(X[:,:, 1], kernel_H2, 'same')**2 + signal.convolve(X[:,:, 2], kernel_H2, 'same')**2)
    eCy = np.maximum(signal.convolve(Y[:,:, 1], kernel_V1, 'same')**2 + signal.convolve(Y[:,:, 2], kernel_V1, 'same')**2, signal.convolve(Y[:,:, 1], kernel_V2, 'same')**2 + signal.convolve(Y[:,:, 2], kernel_V2, 'same')**2)
    eC = np.minimum(eCx, eCy)
    eC = eC**0.5

    return eL, eC
def MNballset(delta):
    # MNballset returns a set of convolution filters describing
    # the relative locations in the elements of the ball set.
    # This algorithm was developed according to Hirakawa's master's
    # thesis.
    index = delta
    H = np.zeros((index*2+1, index*2+1, (index*2+1)**2))  # initialize

    k = 0
    for i in range(-index, index+1):
        for j in range(-index, index+1):
            if(np.linalg.norm([i,j]) <= delta):
                p = np.linalg.norm([i,j])
                H[index+i, index+j, k] = 1  # included
                k = k + 1

    H = H[:, :, 0:k]
    return H
    # 计算相似度ƒ
def MNhomogeneity(LAB_image, delta, epsilonL, epsilonC):
    X = LAB_image

    epsilonC_sq = epsilonC**2

    H = MNballset(delta)

    h, w, c = LAB_image.shape
    K = np.zeros((h,w))

    kh, kw, kc = H.shape

    # 注意浮点数精度可能会有影响
    for i in range(kc):
        L = np.abs(signal.convolve(X[:,:, 0], H[:,:, i], 'same')-X[:,:, 0]) <= epsilonL  # level set
        C = ((signal.convolve(X[:,:, 1], H[:,:, i], 'same')-X[:,:, 1])**2 + (signal.convolve(X[:,:, 2], H[:,:, i], 'same')-X[:,:, 2])** 2) <= epsilonC_sq  # color set
        U = C & L  # metric neighborhood
        K = K + U  # homogeneity
    return K
    # 去artifact
def MNartifact(R, G, B, iterations):
    h, w = R.shape
    Rt = np.zeros((h, w, 8))
    Bt = np.zeros((h, w, 8))
    Grt = np.zeros((h, w, 4))
    Gbt = np.zeros((h, w, 4))

    kernel_1 = np.array([[1, 0, 0], [0, 0, 0], [0, 0, 0]])
    kernel_2 = np.array([[0, 1, 0], [0, 0, 0], [0, 0, 0]])
    kernel_3 = np.array([[0, 0, 1], [0, 0, 0], [0, 0, 0]])
    kernel_4 = np.array([[0, 0, 0], [1, 0, 0], [0, 0, 0]])
    kernel_5 = np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]])
    kernel_6 = np.array([[0, 0, 0], [0, 0, 1], [0, 0, 0]])
    kernel_7 = np.array([[0, 0, 0], [0, 0, 0], [1, 0, 0]])
    kernel_8 = np.array([[0, 0, 0], [0, 0, 0], [0, 1, 0]])
    kernel_9 = np.array([[0, 0, 0], [0, 0, 0], [0, 0, 1]])

    for i in range(iterations):
        Rt[:, :, 0] = signal.convolve(R - G, kernel_1, 'same')
        Rt[:, :, 1] = signal.convolve(R - G, kernel_2, 'same')
        Rt[:, :, 2] = signal.convolve(R - G, kernel_3, 'same')
        Rt[:, :, 3] = signal.convolve(R - G, kernel_4, 'same')
        Rt[:, :, 4] = signal.convolve(R - G, kernel_6, 'same')
        Rt[:, :, 5] = signal.convolve(R - G, kernel_7, 'same')
        Rt[:, :, 6] = signal.convolve(R - G, kernel_8, 'same')
        Rt[:, :, 7] = signal.convolve(R - G, kernel_9, 'same')
        Rm = np.median(Rt, axis=2)
        R = G + Rm
        Bt[:, :, 0] = signal.convolve(B - G, kernel_1, 'same')
        Bt[:, :, 1] = signal.convolve(B - G, kernel_2, 'same')
        Bt[:, :, 2] = signal.convolve(B - G, kernel_3, 'same')
        Bt[:, :, 3] = signal.convolve(B - G, kernel_4, 'same')
        Bt[:, :, 4] = signal.convolve(B - G, kernel_6, 'same')
        Bt[:, :, 5] = signal.convolve(B - G, kernel_7, 'same')
        Bt[:, :, 6] = signal.convolve(B - G, kernel_8, 'same')
        Bt[:, :, 7] = signal.convolve(B - G, kernel_9, 'same')
        Bm = np.median(Bt, axis=2)
        B = G + Bm

        Grt[:, :, 0] = signal.convolve(G - R, kernel_2, 'same')
        Grt[:, :, 1] = signal.convolve(G - R, kernel_4, 'same')
        Grt[:, :, 2] = signal.convolve(G - R, kernel_6, 'same')
        Grt[:, :, 3] = signal.convolve(G - R, kernel_8, 'same')
        Grm = np.median(Grt, axis=2)
        Gr = R + Grm

        Gbt[:, :, 0] = signal.convolve(G - B, kernel_2, 'same')
        Gbt[:, :, 1] = signal.convolve(G - B, kernel_4, 'same')
        Gbt[:, :, 2] = signal.convolve(G - B, kernel_6, 'same')
        Gbt[:, :, 3] = signal.convolve(G - B, kernel_8, 'same')
        Gbm = np.median(Gbt, axis=2)
        Gb = B + Gbm
        G = (Gr + Gb) / 2

    return R, G, B
    # internal function
def labf(t):
    d = t**(1/3)
    index = np.where(t<=0.008856)
    d[index] = 7.787 * t[index] + 16 / 116
    return d
def RGB2LAB(X):
    a = np.array([
        [3.40479,  -1.537150, -0.498535],
        [-0.969256, 1.875992, 0.041556],
        [0.055648, -0.204043, 1.057311]])
    ai = np.linalg.inv(a)

    h, w, c = X.shape

    R = X[:, :, 0]
    G = X[:, :, 1]
    B = X[:, :, 2]
    planed_R = R.flatten()
    planed_G = G.flatten()
    planed_B = B.flatten()

    planed_image = np.zeros((c, h*w))
    planed_image[0, :] = planed_R
    planed_image[1, :] = planed_G
    planed_image[2, :] = planed_B

    planed_lab = np.dot(ai, planed_image)
    planed_1 = planed_lab[0, :]
    planed_2 = planed_lab[1, :]
    planed_3 = planed_lab[2, :]
    L1 = np.reshape(planed_1, (h,w))
    L2 = np.reshape(planed_2, (h,w))
    L3 = np.reshape(planed_3, (h,w))

    # color  space conversion  into LAB
    result_lab = np.zeros((h, w, c))
    result_lab[:,:, 0] = 116 * labf(L2/255)-16
    result_lab[:,:, 1] = 500 * (labf(L1/255) - labf(L2/255))
    result_lab[:,:, 2] = 200 * (labf(L2/255) - labf(L3/255))

    return result_lab
#AHD https://blog.csdn.net/liuxingwan/article/details/108842247
def AHD(img, pattern):
    delta = 2
    gamma = 1
    maxvalue = 4095
    iterations = 2
    imgh, imgw = img.shape
    imgs = 10

    # X,Y方向插值
    # 扩展大小
    # Y = [X(n + 1:-1: 2, n + 1: -1 :2,:)       X(n + 1: -1:2,:,:)        X(n + 1: -1:2, end - 1: -1:end - n,:)
    # X(:, n + 1: -1:2,:)               X                                X(:, end - 1: -1:end - n,:)
    # X(end - 1: -1:end - n, n + 1: -1:2,:)  X(end - 1: -1:end - n,:,:)  X(end - 1: -1:end - n, end - 1: -1:end - n,:)];
    f = np.pad(img, ((imgs,imgs),(imgs,imgs)), 'reflect')

    Yx = AH_interpolateX(f, pattern, gamma, maxvalue)
    Yy = AH_interpolateY(f, pattern, gamma, maxvalue)

    # 转LAB
    YxLAB = RGB2LAB(Yx)
    YyLAB = RGB2LAB(Yy)

    # 色彩差异的运算
    epsilonL, epsilonC = MNparamA(YxLAB, YyLAB)
    Hx = MNhomogeneity(YxLAB, delta, epsilonL, epsilonC)
    Hy = MNhomogeneity(YyLAB, delta, epsilonL, epsilonC)
    f_kernel = np.ones((3, 3))
    Hx = signal.convolve(Hx, f_kernel, 'same')
    Hy = signal.convolve(Hy, f_kernel, 'same')

    # 选择X,Y
    # set output initially to Yx
    R = Yx[:,:, 0]
    G = Yx[:,:, 1]
    B = Yx[:,:, 2]
    Ry = Yy[:,:, 0]
    Gy = Yy[:,:, 1]
    By = Yy[:,:, 2]
    # color_show(Yx, 255)
    # color_show(Yy, 255)

    # set output to Yy if Hy >= Hx
    # 所有的都找到
    bigger_index = np.where(Hy >= Hx)
    Rs = R
    Gs = G
    Bs = B
    Rs[bigger_index] = Ry[bigger_index]
    Gs[bigger_index] = Gy[bigger_index]
    Bs[bigger_index] = By[bigger_index]

    h, w = Rs.shape
    YT = np.zeros((h, w, 3))
    YT[:, :, 0] = Rs
    YT[:, :, 1] = Gs
    YT[:, :, 2] = Bs
    # color_show(YT, 255)

    # 去掉artifact
    Rsa, Gsa, Bsa = MNartifact(Rs, Gs, Bs, iterations)  # find
    # R and B
    # values

    Y = np.zeros((h, w, 3))
    Y[:,:, 0] = Rsa
    Y[:,:, 1] = Gsa
    Y[:,:, 2] = Bsa

    # 调整size和值的范畴
    Y = np.clip(Y, 0, maxvalue)
    # color_show(Y, 255)
    resultY = Y[imgs:imgs+imgh, imgs:imgs+imgw, :]
    return resultY

GBTF

def GBTF(img, pattern):
    img = img.astype(np.float64)
    R_m, G_m, B_m = masks_Bayer(img, pattern)
    R_m = img * R_m
    G_m = img * G_m
    B_m = img * B_m
    h, w = img.shape
    R = np.zeros((h, w))
    G = np.zeros((h, w))
    B = np.zeros((h, w))
    Rh = np.zeros((h, w))
    Gh = np.zeros((h, w))
    Bh = np.zeros((h, w))
    Rv = np.zeros((h, w))
    Gv = np.zeros((h, w))
    Bv = np.zeros((h, w))
    diffh = np.zeros((h, w))
    diffv = np.zeros((h, w))
    Grdiff = np.zeros((h, w))
    Gbdiff = np.zeros((h, w))
    for i in range(2, h - 2):
        for j in range(2, w - 2):
            if ((i % 2 == 0) & (j % 2 == 0)):
                Gv[i, j] = ((G_m[i - 1, j] + G_m[i + 1, j]) / 2)+(2*R_m[i, j] - (R_m[i-2,j]+R_m[i+2,j]))/4
                Gh[i, j] = ((G_m[i, j - 1] + G_m[i, j + 1]) / 2)+(2*R_m[i, j] - (R_m[i,j-2]+R_m[i,j+2]))/4
                diffh[i,j] = Gh[i, j]-R_m[i, j]
                diffv[i, j] = Gv[i, j]- R_m[i, j]
            elif ((i % 2 == 1) & (j % 2 == 1)):
                Gv[i, j] = ((G_m[i - 1, j] + G_m[i + 1, j]) / 2)+(2*B_m[i, j] - (B_m[i-2,j]+B_m[i+2,j]))/4
                Gh[i, j] = ((G_m[i, j - 1] + G_m[i, j + 1]) / 2)+(2*B_m[i, j] - (B_m[i,j-2]+B_m[i,j+2]))/4
                diffh[i,j] = Gh[i, j]-B_m[i, j]
                diffv[i, j] = Gv[i, j]- B_m[i, j]
            elif ((i % 2 == 0) & (j % 2 == 1)):
                Rh[i, j] = ((R_m[i, j - 1] + R_m[i, j + 1]) / 2) + (2 * G_m[i, j] - (G_m[i, j - 2] + G_m[i, j + 2])) / 4
                Bv[i, j] = ((B_m[i - 1, j] + B_m[i + 1, j]) / 2) + (2 * G_m[i, j] - (G_m[i - 2, j] + G_m[i + 2, j])) / 4
                diffh[i,j] = G_m[i, j]-Rh[i, j]
                diffv[i, j] = G_m[i, j]- Bv[i, j]
            else:
                Rv[i, j] = ((R_m[i - 1, j] + R_m[i + 1, j]) / 2) + (2 * G_m[i, j] - (G_m[i - 2, j] + G_m[i + 2, j])) / 4
                Bh[i, j] = ((B_m[i, j - 1] + B_m[i, j + 1]) / 2) + (2 * G_m[i, j] - (G_m[i, j - 2] + G_m[i, j + 2])) / 4
                diffh[i, j] = G_m[i, j] - Bh[i, j]
                diffv[i, j] = G_m[i, j] - Rv[i, j]
#G-R/G-B的任意一个代表梯度
    graddiffh = np.zeros((h, w))
    graddiffv = np.zeros((h, w))
    for i in range(2, h - 2):
        for j in range(2, w - 2):
            graddiffh[i,j] = np.abs(diffh[i, j+1]-diffh[i, j-1])
            graddiffv[i, j] = np.abs(diffv[i+1, j] - diffv[i-1,j])
#将水平和垂直方向色差梯度和倒数作为权值
    for i in range(4, h - 4):
        for j in range(4, w - 4):
            weight=np.array([0,0,0,0])
            weight = weight.astype(np.float64)
            sum=0
            startpoint=np.array([[-4,-2],[0,-2],[-2,-4],[-2,0]])
            if (((i % 2 == 0) & (j % 2 == 0))|((i % 2 == 1) & (j % 2 == 1))):
                for dir in range(4):
                    startx=startpoint[dir,0]
                    starty = startpoint[dir, 1]
                    for x in range(5):
                        for y in range(5):
                            if(dir<2):
                                weight[dir]+=graddiffv[i+x+startx, j+y+starty]
                            else:
                                weight[dir] += graddiffh[i+x + startx, j+y + starty]
                    weight[dir]*=weight[dir]
                    weight[dir] = 1.0 / weight[dir]
                    sum += weight[dir]
                #G通道插值
                    if ((i % 2 == 0) & (j % 2 == 0)):
                        Grdiff[i, j]=(weight[0]*(diffv[i-4,j]+diffv[i,j]+diffv[i-3,j]+diffv[i-2,j]+diffv[i-1,j])+
                                                       weight[1] * (
                                                       diffv[i, j] + diffv[i+1, j] + diffv[i+2, j] + diffv[i+3, j] + diffv[
                                                           i+4, j])+
                                                       weight[2] * (
                                                       diffh[i, j-4] + diffh[i, j-3] + diffh[i, j-2] + diffh[i, j-1] + diffh[
                                                           i, j])+
                                                       weight[3] * (
                                                       diffh[i, j] + diffh[i, j+1] + diffh[i, j+2] + diffh[i, j+3] + diffh[
                                                           i, j+4]))/(5*sum)
                        G[i, j] = R_m[i, j]+Grdiff[i, j]
                    elif ((i % 2 == 1) & (j % 2 == 1)):
                        Gbdiff[i,j]=(weight[0]*(diffv[i-4,j]+diffv[i,j]+diffv[i-3,j]+diffv[i-2,j]+diffv[i-1,j])+
                                                       weight[1] * (
                                                       diffv[i, j] + diffv[i+1, j] + diffv[i+2, j] + diffv[i+3, j] + diffv[
                                                           i+4, j])+
                                                       weight[2] * (
                                                       diffh[i, j-4] + diffh[i, j-3] + diffh[i, j-2] + diffh[i, j-1] + diffh[
                                                           i, j])+
                                                       weight[3] * (
                                                       diffh[i, j] + diffh[i, j+1] + diffh[i, j+2] + diffh[i, j+3] + diffh[
                                                           i, j+4]))/(5*sum)
                        G[i, j] = B_m[i, j] +Gbdiff[i,j]
#R/B B/R插值
    G=G+G_m
    PrbData= np.array(
        [[0, 0, -0.03125, 0, -0.03125, 0, 0],
        [0, 0, 0, 0, 0, 0, 0],
        [-0.03125, 0, 0.3125, 0, 0.3125, 0, -0.03125],
        [0, 0, 0, 0, 0, 0, 0],
        [-0.03125, 0, 0.3125, 0, 0.3125, 0, -0.03125],
        [0, 0, 0, 0, 0, 0, 0],
        [0, 0, -0.03125, 0, -0.03125, 0, 0]]
    )
    for i in range(4, h - 4):
        for j in range(4, w - 4):
            if ((i % 2 == 0) & (j % 2 == 0)):
                B[i,j] = G[i,j]-np.sum(Gbdiff[i-3:i+4,j-3:j+4]*PrbData)
                R[i,j] = R_m[i,j]
            elif ((i % 2 == 1) & (j % 2 == 1)):
                R[i,j] = G[i,j]-np.sum(Grdiff[i-3:i+4,j-3:j+4]*PrbData)
                B[i,j] = B_m[i,j]

    for i in range(4, h - 4):
        for j in range(4, w - 4):
            if ((i % 2 == 0) & (j % 2 == 1)):
                B[i,j]=G[i,j]+(B[i-1,j]-G[i-1,j]+B[i+1,j]-G[i+1,j]+B[i,j-1]-G[i,j-1]+B[i,j+1]-G[i,j+1])/4
                R[i, j] = G[i, j] + (R[i - 1, j] - G[i - 1, j] + R[i + 1, j] - G[i + 1, j] + R[i, j - 1] - G[i, j - 1] +
                                     R[i, j + 1] - G[i, j + 1]) / 4
            if ((i % 2 == 1) & (j % 2 == 0)):
                B[i,j]=G[i,j]+(B[i-1,j]-G[i-1,j]+B[i+1,j]-G[i+1,j]+B[i,j-1]-G[i,j-1]+B[i,j+1]-G[i,j+1])/4
                R[i, j] = G[i, j] + (R[i - 1, j] - G[i - 1, j] + R[i + 1, j] - G[i + 1, j] + R[i, j - 1] - G[i, j - 1] +
                                     R[i, j + 1] - G[i, j + 1]) / 4

    result_img = np.zeros((h, w, 3))
    # R, G, B = MNartifact(R, G, B, 2)
    result_img[:, :, 0] = R
    result_img[:, :, 1] = G
    result_img[:, :, 2] = B
    del R_m, G_m, B_m
    return result_img

结果对比

测试图：Kodak18	CPSNR	MAE
色差恒定法(Kr-Kb)	32.92	3.10
一阶微分边缘插值	32.47	3.22
二阶微分边缘插值	32.52	3.21
梯度矫正线性插值	33.20	3.16
Adam-Hamilton	34.71	2.77
AHD	37.16	2.39
GBTF	39.51	1.91

对比客观数据指标可以看到一阶微分、二阶微分这类在色差法上简单增加边缘判断的算法以及梯度矫正线性插值这种单纯经验性滤波核的方法并没有带来很大的效果提升，与基础插值法里PSNR最高的色差恒定法比起来也只增加了1-2dB。而Adam-Hamilton、AHD、GBTF这种结合多种基础算法思想于一体的边缘处理方法则有着明显客观指标提升。