• 中文
  • 3.0.31

    • 去马赛克示例介绍

    去马赛克是一种图像处理技术，用于将Bayer格式的图像数据转换为RGB格式的图像数据。去马赛克的原理是通过对Bayer格式的图像数据进行插值处理，生成RGB格式的图像数据，以便于后续的图像处理。 为了更好的展示如何使用DSP实现去马赛克的功能，本章从算法原理出发，展示了详细的去马赛克算法实现流程。

    原理介绍

    Bayer格式的图像数据是一种单通道的图像数据，其中每个像素点只包含一种颜色的信息， 常见的Bayer格式有RGGB、BGGR、GRBG、GBRG等。而RGB格式的图像数据是一种三通道的图像数据，其中每个像素点包含红、绿、蓝三种颜色的信息。 根据插值算法的不同，去马赛克算法可以分为最近邻插值、双线性插值、双三次插值等，在本示例中采用了一种利用色差恒定原理的自定义插值方式。

    插值算法分析

    在此示例中，RAW图片数据排布为RGGB、comp模式、每像素12bit。其具体示意图片如下：

    raw_rggb

    RGB格式的输出为R、G、B三通道，每通道16bit，每种通道信息独立存储。其具体示意图片如下：

    rgb_p

    下面以上图为基准，介绍插值计算的实现方法。 在插值计算中，Raw图片的每个位置都需要计算出其对应的RGB值，根据不同位置的计算方式将其分为以下几类：

    • Raw格式R计算R G B。
    • Raw格式G计算R G B。
    • Raw格式B计算R G B。

    Raw格式R计算R G B

    R通道的值为Raw的值，例如Raw图中的R00，R02点。

    G通道的值需要根据Raw的值和周围的值进行插值计算，将此插值计算方法命名为green_5x5。 以R24为例，其计算方式如下：

    • 首先判断插值方向,插值方向根据水平与垂直方向的一阶梯度变化加2倍二阶梯度变化的和作为判断依据:
      • 水平方向梯度变化：dh = abs((G23-G25)/2)+abs((R22+R26-2*R24)/2)。
      • 垂直方向梯度变化：dv = abs((G14-G34)/2)+abs((R04+R44-2*R24)/2)。
    • 当dh < dv，沿着水平方向插值：G = (G23+G25)/2。
    • 当dh > dv，沿着垂直方向插值：G = (G14+G34)/2。
    • 当dh == dv，G = (G23+G25+G14+G34)/4。
    • 插值完成，将计算出来的值clip到有效范围内，例如12bit数据，clip到0~4095范围。

    将上述插值用代码表示如下：

    // G[i][j] -- green_5x5[i][j]
dh = abs(G[i][j-1] - G[i][j+1]) + abs(R[i][j-2] - R[i][j] + R[i][j+2] - R[i][j]);
dv = abs(G[i-1][j] - G[i+1][j]) + abs(R[i-2][j] - R[i][j] + R[i+2][j] - R[i][j]);
if (dh < dv) {
    G[i][j] = (G[i][j-1] + G[i][j+1]) / 2;
} else if (dh > dv) {
    G[i][j] = (G[i-1][j] + G[i+1][j]) / 2;
} else {
    G[i][j] = (G[i][j-1] + G[i][j+1] + G[i-1][j] + G[i+1][j]) / 4;
}

    B通道的值需要根据Raw的值和周围的值进行插值计算，以R44为例，具体的计算方式为 B = G + (B33-G33+B35-G35+B53-G53+B55-G55)/4。

    其中G33、G35、G53、G55分别为对应位置G通道的值，计算方法命名为green_3x3。 以B33为例，green_3x3计算方式如下：

    • 计算水平方向梯度变化量绝对值，dh = abs(G32-G34)。
    • 计算垂直方向梯度变化量绝对值，dv = abs(G23-G43)。
    • 当dh < dv，沿着水平方向插值：G = (G32+G34)/2。
    • 当dh > dv，沿着垂直方向插值：G = (G23+G43)/2。
    • 当dh == dv，G = (G32+G34+G23+G43)/4。
    • 插值完成，将计算出来的值clip到有效范围内，例如12bit数据，clip到0~4095范围。

    将上述插值用代码表示如下：

    // G[i][j] -- green_3x3
int dh = abs(G[i][j-1] - G[i][j+1]);
int dv = abs(G[i-1][j] - G[i+1][j]);
if (dh < dv) {
    G[i][j] = (G[i][j-1] + G[i][j+1]) / 2;
} else if (dh > dv) {
    G[i][j] = (G[i-1][j] + G[i+1][j]) / 2;
} else {
    G[i][j] = (G[i][j-1] + G[i][j+1] + G[i-1][j] + G[i+1][j]) / 4;
}

// G[i-1][j-1] -- calculate by green_3x3 
// G[i-1][j+1] -- calculate by green_3x3 
// G[i+1][j-1] -- calculate by green_3x3 
// G[i+1][j+1] -- calculate by green_3x3 
B[i][j] = G[i][j] + 
          (B[i - 1][j - 1] - G[i - 1][j - 1] +
           B[i - 1][j + 1] - G[i - 1][j + 1] +
           B[i + 1][j - 1] - G[i + 1][j - 1] +
           B[i + 1][j + 1] - G[i + 1][j + 1]) / 4;

    Raw格式G计算R G B

    根据G的位置不同，将G分为G1和G2两种情况，G1为R同行的G通道，G2为B同列的G通道。

    R通道的值需要根据Raw的值和周围的值进行插值计算， 以G23为例，具体的计算方式如下：

    • 计算左侧位置的G通道值G22，计算方法为green_3x3。
    • 计算右侧位置的G通道值G24，计算方法为green_3x3。
    • R = G + (R22-G22+R24-G24)/2。

    将上述插值用代码表示如下：

    // G[i][j-1] -- calculate by green_3x3
// G[i][j+1] -- calculate by green_3x3
R[i][j] = G[i][j] + (R[i][j-1] - G[i][j-1] + R[i][j+1] - G[i][j+1]) / 2;

    以G32为例，具体的计算方式如下：

    • 计算上侧位置的G通道值G22，计算方法为green_3x3。
    • 计算下侧位置的G通道值G42，计算方法为green_3x3。
    • R = G + (R22-G22+R42-G42)/2。

    将上述插值用代码表示如下：

    // G[i-1][j] -- calculate by green_3x3
// G[i+1][j] -- calculate by green_3x3
R[i][j] = G[i][j] + (R[i-1][j] - G[i-1][j] + R[i+1][j] - G[i+1][j]) / 2;

    G通道的值为Raw的值，例如Raw图中的R01，R10点。

    B通道的值需要根据Raw的值和周围的值进行插值计算， 以G23为例，具体的计算方式如下：

    • 计算上侧位置的G通道值G13，计算方法为green_3x3。
    • 计算下侧位置的G通道值G33，计算方法为green_3x3。
    • B = G + (B13-G13+B33-G33)/2。

    将上述插值用代码表示如下：

    // G[i-1][j] -- calculate by green_3x3
// G[i+1][j] -- calculate by green_3x3
B[i][j] = G[i][j] + (B[i-1][j] - G[i-1][j] + B[i+1][j] - G[i+1][j]) / 2;

    以G32为例，具体的计算方式如下：

    • 计算左侧位置的G通道值G31，计算方法为green_3x3。
    • 计算右侧位置的G通道值G33，计算方法为green_3x3。
    • B = G + (B31-G31+B33-G33)/2。

    将上述插值用代码表示如下：

    // G[i][j-1] -- calculate by green_3x3
// G[i][j+1] -- calculate by green_3x3
B[i][j] = G[i][j] + (B[i][j-1] - G[i][j-1] + B[i][j+1] - G[i][j+1]) / 2;

    Raw格式B计算R G B

    G通道的值需要根据Raw的值和周围的值进行插值计算，计算方法为green_5x5。

    R通道的值需要根据Raw的值和周围的值进行插值计算，以B33为例，具体的计算方式为 R = G + (R22-G22+R24-G24+R42-G42+R44-G44)/4。

    其中G22、G24、G42、G44分别为对应位置G通道的值，计算方法为green_3x3。

    将上述插值用代码表示如下：

    // G[i-1][j-1] -- calculate by green_3x3 
// G[i-1][j+1] -- calculate by green_3x3 
// G[i+1][j-1] -- calculate by green_3x3 
// G[i+1][j+1] -- calculate by green_3x3 
B[i][j] = G[i][j] + 
          (B[i - 1][j - 1] - G[i - 1][j - 1] +
           B[i - 1][j + 1] - G[i - 1][j + 1] +
           B[i + 1][j - 1] - G[i + 1][j - 1] +
           B[i + 1][j + 1] - G[i + 1][j + 1]) / 4;

    B通道的值为Raw的值，例如Raw图中的B11、B13点。

    DSP 实现流程

    DSP的实现主要围绕着以下几点进行优化：

    1. 通过使用tile_manager框架，将数据分块处理，并将分块数据的搬运和计算构造pingpong流水线。
    2. 通过SIMD方法实现green_3x3、green_5x5插值方法。
    3. 使用数据分类等方法，将相同数据的计算过程进行合并，充分使用宽位寄存器的批量处理能力。

    tile manager框架的使用

    tile manager是cadence提供的一种算子开发框架，通过将整体数据分块处理，将数据的搬运和计算构造pingpong流水线，提高数据处理效率。 其中frame是框架中对完整输入或者输出的描述，tile是框架中对frame进行分块的描述，是框架中的计算处理单元。 算子开发中通常按照如下步骤使用：

    • 创建算子参数结构体
    • 实现tile分块函数，将frame数据按照顺序进行分块处理
    • 实现每个tile的计算函数，完成单个tile的处理
    • 通过tile manager相关接口，构造pingpong流水线，实现完整算子功能

    在Demosaicing算子中，Frame结构体中包含一个输入frame和三个输出frame，Tile结构体种包含一个输入tile和三个输出tile。 这样就满足了在每个输入tile中都可以计算出对应的三个RGB tile的目的。 具体参数结构体实现如下：

    typedef struct {
  xi_frame input_frame;
  xi_frame output_frame_r;
  xi_frame output_frame_g;
  xi_frame output_frame_b;
} hbCommonArgDemosaic;

typedef struct {
  xi_tile src_tile;
  xi_tile dst_tile_r;
  xi_tile dst_tile_g;
  xi_tile dst_tile_b;
} hbTileArgDemosaic;

    下一步要实现的是tile的分块函数，按照一定逻辑将frame数据分成逐个处理的块。在Demosaicing算子中，输入输出的每个像素的位置映射是一对一的，位置也不会发生变化。 例如，Raw图片中的G21点，在RGB图片中的每个通道中坐标也是(2,1)。因此，输入tile的坐标和设置和输出一致。 在Demosaicing计算中，Raw图片的边界值计算需要使用的数据超出了原图大小，因此需要填充边界值。 在上文描述的插值算法中，最多使用的是5x5的范围，因此，输入的边界值大小设置为2，输出的边界值大小设置为0。

    具体实现如下：

    static XI_ERR_TYPE setup_demosaic_tiles(xi_pTile src, xi_pTile dst_r,
                                        xi_pTile dst_g, xi_pTile dst_b,
                                        int update_only_flag) {
  XI_TILE_SET_X_COORD(src, XI_TILE_GET_X_COORD(dst_r));
  XI_TILE_SET_Y_COORD(src, XI_TILE_GET_Y_COORD(dst_r));
  XI_TILE_SET_WIDTH(src, XI_TILE_GET_WIDTH(dst_r));
  XI_TILE_SET_HEIGHT(src, XI_TILE_GET_HEIGHT(dst_r));

  if (update_only_flag == 0) {
    XI_TILE_SET_TYPE(dst_r, XI_S16);
    XI_TILE_SET_TYPE(dst_g, XI_S16);
    XI_TILE_SET_TYPE(dst_b, XI_S16);
    XI_TILE_SET_TYPE(src, XI_S16);

    XI_TILE_SET_EDGE_WIDTH(src, 2);
    XI_TILE_SET_EDGE_HEIGHT(src, 2);
  }
  return XI_ERR_OK;
}

static void setup_updates_demosaic_tiles(xvTileManager *pxv_tm,
                                         RefTile *ref_tile, void *com_args,
                                         void *tiles_args,
                                         int update_only_flag) {
  hbTileArgDemosaic *tile_args = (hbTileArgDemosaic *)tiles_args;
  hbCommonArgDemosaic *common_args = (hbCommonArgDemosaic *)com_args;
  tile_args->dst_tile_r.x = ref_tile->x;
  tile_args->dst_tile_r.y = ref_tile->y;
  tile_args->dst_tile_r.width = ref_tile->tileWidth;
  tile_args->dst_tile_r.height = ref_tile->tileHeight;
  tile_args->dst_tile_r.edgeWidth = 0;
  tile_args->dst_tile_r.edgeHeight = 0;

  tile_args->dst_tile_g.x = ref_tile->x;
  tile_args->dst_tile_g.y = ref_tile->y;
  tile_args->dst_tile_g.width = ref_tile->tileWidth;
  tile_args->dst_tile_g.height = ref_tile->tileHeight;
  tile_args->dst_tile_g.edgeWidth = 0;
  tile_args->dst_tile_g.edgeHeight = 0;

  tile_args->dst_tile_b.x = ref_tile->x;
  tile_args->dst_tile_b.y = ref_tile->y;
  tile_args->dst_tile_b.width = ref_tile->tileWidth;
  tile_args->dst_tile_b.height = ref_tile->tileHeight;
  tile_args->dst_tile_b.edgeWidth = 0;
  tile_args->dst_tile_b.edgeHeight = 0;

  setup_demosaic_tiles(&tile_args->src_tile, &tile_args->dst_tile_r,
                       &tile_args->dst_tile_g, &tile_args->dst_tile_b,
                       update_only_flag);
  if (update_only_flag == 0) {
    xvRegisterTile(pxv_tm, (xvTile *)&(tile_args->src_tile), NULL,
                   (xvFrame *)&(common_args->input_frame), XV_INPUT_TILE);
    xvRegisterTile(pxv_tm, (xvTile *)&(tile_args->dst_tile_r), NULL,
                   (xvFrame *)&(common_args->output_frame_r), XV_OUTPUT_TILE);
    xvRegisterTile(pxv_tm, (xvTile *)&(tile_args->dst_tile_g), NULL,
                   (xvFrame *)&(common_args->output_frame_g), XV_OUTPUT_TILE);
    xvRegisterTile(pxv_tm, (xvTile *)&(tile_args->dst_tile_b), NULL,
                   (xvFrame *)&(common_args->output_frame_b), XV_OUTPUT_TILE);
  }
}

    第三步要实现的是每个tile的计算函数，是算子的核心处理逻辑，对应源码中的demosaic_kernel函数。

    最后，调用框架相关接口，实现多个tile的流水线处理，对应源码中的demosaic_rggb_comp函数。 代码核心如下：

      ret = xvExecuteFullImageKernel(pxv_tm, &refTile, sizeof(hbTileArgDemosaic),
                                 (void *)com_args, setup_updates_demosaic_tiles,
                                 demosaic_kernel, 0, 1, XV_AUTOSIZE_FLAG);

    实现插值算法green_3x3、green_5x5

    在上文对插值算法的分析中，可以提取出核心的插值逻辑 green_3x3 和 green_5x5。通过SIMD优化，实现代码如下：

    int32_t green_3x3(xb_vecNx16 &val, xb_vecNx16 &l, xb_vecNx16 &t, xb_vecNx16 &r,
                  xb_vecNx16 &b) {
  /* 3x3 kernel
    *     top     *
    left  mid     right
    *     bottom  *
    */
  // abs(*left - *right)
  val = IVP_ANDNX16(val, 0);
  xb_vecNx16 lr = IVP_ABSNX16(IVP_SUBNX16(l, r));
  xb_vecNx16 tb = IVP_ABSNX16(IVP_SUBNX16(t, b));
  vboolN vlt = IVP_LTNX16(lr, tb);
  vboolN vgt = IVP_GTNX16(lr, tb);
  vboolN veq = (lr == tb);  // IVP_EQNX16(lr, tb);

  // (*left + *right + *top + *bottom) / 4;
  IVP_ADDNX16T(val, l, r, veq);
  IVP_ADDNX16T(val, val, b, veq);
  IVP_ADDNX16T(val, val, t, veq);
  val = val >> 1;

  // (*left + *right) / 2;
  // (*top + *bottom) / 2;
  IVP_ADDNX16T(val, l, r, vlt);
  IVP_ADDNX16T(val, t, b, vgt);
  val = val >> 1;
  return XVTM_SUCCESS;
}

int32_t green_5x5(xb_vecNx16 &val, xb_vecNx16 &mid, xb_vecNx16 &l,
                  xb_vecNx16 &l2, xb_vecNx16 &t, xb_vecNx16 &t2, xb_vecNx16 &r,
                  xb_vecNx16 &r2, xb_vecNx16 &b, xb_vecNx16 &b2) {
  /* 5x5 kernel
    *      *     top2     *      *
    *      *     top      *      *
    left2  left  mid      right  right2
    *      *     bottom   *      *
    *      *     bottom2  *      *
    */
  // abs(*left - *right) + abs(*left2 + *right2 - 2 * *ptr)
  xb_vecNx16 lr =
      IVP_ADDNX16(IVP_ABSNX16(IVP_SUBNX16(l, r)),
                  IVP_ABSNX16(IVP_SUBNX16(IVP_ADDNX16(l2, r2), (mid) << 1)));
  xb_vecNx16 tb =
      IVP_ADDNX16(IVP_ABSNX16(IVP_SUBNX16(t, b)),
                  IVP_ABSNX16(IVP_SUBNX16(IVP_ADDNX16(t2, b2), (mid) << 1)));
  vboolN vlt = IVP_LTNX16(lr, tb);
  vboolN vgt = IVP_GTNX16(lr, tb);
  vboolN veq = (lr == tb);  // IVP_EQNX16(lr, tb);

  // (*left + *right + *top + *bottom) / 4;
  IVP_ADDNX16T(val, l, r, veq);
  IVP_ADDNX16T(val, val, b, veq);
  IVP_ADDNX16T(val, val, t, veq);
  val = val >> 1;

  // (*left + *right) / 2;
  // (*top + *bottom) / 2;
  IVP_ADDNX16T(val, l, r, vlt);
  IVP_ADDNX16T(val, t, b, vgt);
  val = val >> 1;
  return XVTM_SUCCESS;
}

    其中，函数的输入参数为计算插值所需要的不同位置的数组，每个元素计算所需的值根据数组下标一一对应。

    主体计算逻辑合并和SIMD优化

    根据分析，在计算插值过程中，可以将奇数和偶数位置的数据分开处理，比如偶数行元素为RG，可以通过奇偶性将R和G分开进行计算，最后再将计算结果合并。 在插值计算中，每个位置的计算最多涉及到周围5x5的内核大小，因此定义5x6的混合内核如下：

      1  2  3  4  5  6
  7  8  9  10 11 12
  13 14 15 16 17 18
  19 20 21 22 23 24
  25 26 27 28 29 30

    算子的处理按照DSP中int16的最大处理长度进行，因此第一层循环按照宽度，步长为1024/16 = 64，第二层循环按照高度进行，步长为1。 可以看出，在内层的处理逻辑中，混合内核的前4行数据可以复用，可以节省重复的计算逻辑，减少耗时。 流程大致如下：

    1. 读取数据，并按照奇偶位置将数据进行拆分。
    2. 将拆分后的数据进行判断，执行不同的去马赛克逻辑。
    3. 计算结果按照奇偶位置进行合并，并保存到输出数据中。
    4. 将插值内核中的数据进行更新，充分利用本次计算结果，进行下一轮循环。

    第一步，读取数据，并按照奇偶位置将数据进行拆分。DSP实现中将读取数据并按照奇偶分开的过程封装为了LOAD_AND_SELECT_ODD_EVEN_WITH_OFFSET宏函数， 该函数的功能为读取最多64*3=192个数据，并按照奇偶位置拆分，将拆分后的数据按照混合内核的位置进行选择，保存到对应的寄存器中。 其定义及使用代码如下：

    #define LOAD_AND_SELECT_ODD_EVEN_WITH_OFFSET(v0, v1, v2, v3, v4, v5, load_ptr) \
  a_load_src = IVP_LANX16_PP(load_ptr);                                        \
  IVP_LAVNX16_XP(load_v1, a_load_src, load_ptr,                                \
                 2 * ((width - x) + 2 * src_edge_w));                          \
  IVP_LAVNX16_XP(load_v2, a_load_src, load_ptr,                                \
                 2 * ((width - x) - XCHAL_IVPN_SIMD_WIDTH + 2 * src_edge_w));  \
  IVP_LAVNX16_XP(                                                              \
      load_tail, a_load_src, load_ptr,                                         \
      2 * ((width - x) - 2 * XCHAL_IVPN_SIMD_WIDTH + 2 * src_edge_w));         \
  IVP_DSELNX16I(v1, v0, load_v2, load_v1, IVP_DSELI_DEINTERLEAVE_1);           \
  IVP_DSELNX16_OFFSET(v2, v3, load_v1, load_v2, load_tail, 2);                 \
  IVP_DSELNX16_OFFSET(v4, v5, load_v1, load_v2, load_tail, 4);


// load 5x6 kernel row 1
src_vptr = OFFSET_PTR_NX16(psrc_data, -src_edge_h, sstride, -src_edge_w + x);
LOAD_AND_SELECT_ODD_EVEN_WITH_OFFSET(p_temp, p2, p3, p4, p5, p_temp, src_vptr);

    第二步，将拆分后的数据进行判断，执行不同的去马赛克逻辑。在计算插值过程中，相同颜色的计算逻辑是一致的，因此将相同颜色的计算逻辑合并。 在具体实现中，行的奇偶判断取y起始坐标加y偏移坐标的最后一位，即(y+ys)%2。列的奇偶判断取x偏移的最后一位，即(x+xs)%2。因为列的处理步长一定是偶数， 所以x偏移不会影响奇偶性。具体实现代码如下：

    
uint8_t y_odd_flag = (y ^ ys) & 1;  // (y + ys) % 2;
xb_vecNx16 dst_vr1, dst_vg1, dst_vb1, dst_vr2, dst_vg2, dst_vb2;
if (y_odd_flag) {
  if (xs & 1) {
    // x position start from odd position, data is b g b g b g
    // calc rgb from b
    green_5x5(dst_vg1, p15, p14, p13, p9, p3, p16, p17, p21, p27);
    CAL_ROUND_BY_GREEN(dst_vr1, dst_vg1, p8, p10, p20, p22, g33_8, g33_10,
                        p2, p7);
    dst_vb1 = p15;

    // calc rgb from g
    CAL_TB_BY_GREEN(dst_vr2, p16, p10, p22, g33_10, p7);
    dst_vg2 = p16;
    CAL_LR_BY_GREEN(dst_vb2, p16, p15, p17, g33_15, g33_17);
  } else {
    // x position start from even position, data is g b g b g b
    // calc rgb from g
    CAL_TB_BY_GREEN(dst_vr1, p15, p9, p21, g33_9, p5);
    dst_vg1 = p15;
    CAL_LR_BY_GREEN(dst_vb1, p15, p14, p16, g33_14, g33_16);

    // calc rgb from b
    green_5x5(dst_vg2, p16, p15, p14, p10, p4, p17, p18, p22, p28);
    CAL_ROUND_BY_GREEN(dst_vr2, dst_vg2, p9, p11, p21, p23, g33_9, g33_11,
                        p5, p12);
    dst_vb2 = p16;
  }
} else {
  if (xs & 1) {
    // x position start from odd position, data is g r g r g r
    // calc rgb from g
    CAL_LR_BY_GREEN(dst_vr1, p15, p14, p16, g33_14, g33_16);
    dst_vg1 = p15;
    CAL_TB_BY_GREEN(dst_vb1, p15, p9, p21, g33_9, p5);

    // calc rgb from r
    dst_vr2 = p16;
    green_5x5(dst_vg2, p16, p15, p14, p10, p4, p17, p18, p22, p28);
    CAL_ROUND_BY_GREEN(dst_vb2, dst_vg2, p9, p11, p21, p23, g33_9, g33_11,
                        p5, p12);
  } else {
    // x position start from even position, data is r g r g r g
    // calc rgb from r
    dst_vr1 = p15;
    green_5x5(dst_vg1, p15, p14, p13, p9, p3, p16, p17, p21, p27);
    CAL_ROUND_BY_GREEN(dst_vb1, dst_vg1, p8, p10, p20, p22, g33_8, g33_10,
                        p2, p7);

    // calc rgb from g
    CAL_LR_BY_GREEN(dst_vr2, p16, p15, p17, g33_15, g33_17);
    dst_vg2 = p16;
    CAL_TB_BY_GREEN(dst_vb2, p16, p10, p22, g33_10, p7);
  }
}

    第三步，将计算结果按照奇偶位置进行合并，并保存到输出数据中。主要实现方法为通过双选操作，将数据合并。 还需要注意的是，在处理结果之前，要将结果限制到合法的输出范围。具体实现代码如下：

    // clamp data from 16bit to 12bit
IVP_CLAMP_12BIT(dst_vr1);
IVP_CLAMP_12BIT(dst_vg1);
IVP_CLAMP_12BIT(dst_vb1);
IVP_CLAMP_12BIT(dst_vr2);
IVP_CLAMP_12BIT(dst_vg2);
IVP_CLAMP_12BIT(dst_vb2);

// merge data and store
IVP_DSELNX16I(dst_vr2, dst_vr1, dst_vr2, dst_vr1, IVP_DSELI_INTERLEAVE_1);
IVP_DSELNX16I(dst_vg2, dst_vg1, dst_vg2, dst_vg1, IVP_DSELI_INTERLEAVE_1);
IVP_DSELNX16I(dst_vb2, dst_vb1, dst_vb2, dst_vb1, IVP_DSELI_INTERLEAVE_1);

dst_vptr_r = OFFSET_PTR_NX16(pdst_r, y, dstride, x);
IVP_SAVNX16_XP(dst_vr1, a_store_dst, dst_vptr_r, 2 * (width - x));
IVP_SAVNX16_XP(dst_vr2, a_store_dst, dst_vptr_r,
                2 * ((width - x) - XCHAL_IVPN_SIMD_WIDTH));
IVP_SAPOSNX16_FP(a_store_dst, dst_vptr_r);
dst_vptr_g = OFFSET_PTR_NX16(pdst_g, y, dstride, x);
IVP_SAVNX16_XP(dst_vg1, a_store_dst, dst_vptr_g, 2 * (width - x));
IVP_SAVNX16_XP(dst_vg2, a_store_dst, dst_vptr_g,
                2 * ((width - x) - XCHAL_IVPN_SIMD_WIDTH));
IVP_SAPOSNX16_FP(a_store_dst, dst_vptr_g);
dst_vptr_b = OFFSET_PTR_NX16(pdst_b, y, dstride, x);
IVP_SAVNX16_XP(dst_vb1, a_store_dst, dst_vptr_b, 2 * (width - x));
IVP_SAVNX16_XP(dst_vb2, a_store_dst, dst_vptr_b,
                2 * ((width - x) - XCHAL_IVPN_SIMD_WIDTH));
IVP_SAPOSNX16_FP(a_store_dst, dst_vptr_b);

    第四步，将插值内核中的数据进行更新，充分利用本次计算结果，进行下一轮循环。在内层循环中，很多参数是可以复用的， 5x6的混合内核每次循环向下移动一行，因此只需要更新最后一行的数据，寄存器数组、green_3x3的临时结果都可以使用旧值。具体实现代码如下：

    // update green
g33_8 = g33_14;
g33_9 = g33_15;
g33_10 = g33_16;
g33_11 = g33_17;
g33_14 = p2;
g33_15 = p5;
g33_16 = p7;
g33_17 = p12;
// update p data of row 1~4
p3 = p9;
p4 = p10;

p8 = p14;
p9 = p15;
p10 = p16;
p11 = p17;

p13 = p19;
p14 = p20;
p15 = p21;
p16 = p22;
p17 = p23;
p18 = p24;

p19 = p_temp;
p20 = p26;
p21 = p27;
p22 = p28;
p23 = p29;
p24 = p_temp2;

    至此，整个Demosaicing算法的实现流程已经完成。此示例仅供参考，在实际应用场景中，可以根据此示例进行修改和优化，或许可以进一步提升算子效果和性能。