Camera成像中的Demosaic

RAW

Raw指Camera sensor输出的未经处理的原生数据，也就是经过sensor的感光点对光进行采样后输出的数字信号。为了采集RGB，会将RGB的滤镜阵列排列在感光点上。
通常情况下，camera sensor成像输出bayer阵列，得到N*M大小的由RGB组成的bayer阵列。利用Bayer阵列，将其对应的RGB阵列分开，然后经过插值的方法得到RGB三通道的图像。

首先我们实现对Camera Sensor输出的bayer阵列即RAW进行读取。
Color Filter Array的阵列一般有：
Bayer_RGGB、Bayer_GRBG、Bayer_GBRG、Bayer_BGGR等
Bayer阵列即RAW的数据存储，可以根据sensor的datasheet确定，每个pixel的输出大小取决于sensor的ADC（模数转换器）的精度，一般包括：RAW8bit、RAW10bit、RAW12bit、RAW14bit、RAW16bit等格式，数值越大，表示精度越高，相应的sensor的颜色范围更广，对应的图像动态范围等会更好一些。目前手机camera sensor使用的是10/8bit的RAW。
下图表示8bit和10bit对应的一个颜色梯度

RAW10存储读取格式

为了节省空间，对于10bit的RAW的存储和读取，使用5个字节存储4个像素，存储的格式如下图所示，这样的存储即高效有不浪费空间：

RAW12存储读取格式

对于12bit的RAW的存储和读取，使用3个字节存储2个像素，如下图所示：

RAW16存储读取格式

对于16bit的RAW的存储和读取，使用2个字节存储1个像素，值得注意的是读取是需要明确16bit所表示的意义，一般也可能是10bit数据用16bit存储，如果是这样存储，需要知道使用的2bit所占的高位还是地位。

读取RAW10代码

```python
# -*- coding: utf-8 -*-
"""
Created on Thu Aug 27 21:06:48 2020
@author: PRCZ
"""
path = "1_0.RAW"   #GRBG
from time import perf_counter
import numpy as np
from skimage import io
from matplotlib import pyplot as plt
rawData = open(path, 'rb').read()
data_bytes = np.fromstring(rawData,'B')
def one_pixel(input):
    image = np.zeros((4))
    image[0] = (input[0] << 2) | (input[4] & 3)
    image[1] = (input[1] << 2 ) | (input[4]>>2 & 3)
    image[2] = (input[2] << 2 ) | (input[4]>>4 & 3)
    image[3] = (input[3] << 2 ) | (input[4]>>6 & 3)  
    return image
   
def read_raw(data):    
    img = []
    for i in range(0, 3000):
        data = data_bytes[i*5008:(i+1)*5008]
        for j in range(0, 5000, 5):
            input1 = data[j:j+5]
            out = one_pixel(input1)
            img.extend(out)
    img1 = np.matrix(img)
    return img1   
img1 = read_raw(data_bytes)
plt.imshow(img1)

对一个RAW10进行读取，可得下图：

Demosaic

在RAW的基础上进行Demosaic可以得到sensor对应处理后的RGB图，其实是对RGB三通道进行分离，然后对相邻的像素进行插值得到RGB三通的图像。

邻近插值

这里我使用了邻近插值的方法对RGB三通道进行插值。

我是用的RAW的格式是GRBG，可以根据不同的格式调整RGB的顺序。

# -*- coding: utf-8 -*-
"""
Created on Thu Aug 27 21:06:48 2020
@author: PRCZ
"""
path = "1_0.RAW"   #GRBG
from PIL import Image
import os,time,shutil,sys
from time import perf_counter
import numpy as np
from skimage import io
from matplotlib import pyplot as plt

rawData = open(path, 'rb').read()
data_bytes = np.fromstring(rawData,'B')

def G_channel(image, r, c):
    G = np.zeros((3000, 4000))
    G = image
    for i in range(c):
        for j in range(r):
            if (i==0) &(j%2==1) & (j < r-1):
                G[0, j] = 0.5*(image[0, j-1] + image[0, j+1])
            if (i==c-1) &(j%2==0) & (j > 0):
                G[i, j] = 0.5*(image[i, j-1] + image[0, j+1])
            if (j==0) &(i%2==1) & (i < c-1):
                G[i, 0] = 0.5*(image[i-1, 0] + image[i+1, 0])
            if (j==r-1) &(i%2==0) & (i < c-1) & (i >0):
                G[i, 0] = 0.5*(image[i-1, j] + image[i+1, j])
            if (i == c-1) & ( j==0):
                G[i,j] = 0.5*(image[i-1, 0]+image[i, 1])
            if (i == 0) & (j == r-1):
                G[i,j] = 0.5*(image[i, j-1] + image[i+1, j])
            if (i%2 ==1) &(j%2 == 0) &(j >0) &(j<r-1) & (i<c-1):
                G[i,j] = 1/4 * (image[i-1, j] +image[i+1, j]+ image[i, j-1] +image[i, j+1])
            if (i%2==0) & (j%2==1) & (i>0) & (i<c-1) &(j<r-1):
                G[i,j] = 1/4 * (image[i-1, j] +image[i+1, j]+ image[i, j-1] +image[i, j+1])
    return G

def R_channel(image, r, c):
    R = np.zeros((3000, 4000))
    R = image
    for i in range(c):
        for j in range(r):
            if (i%2==1) & (j%2==0) & (i< c-1) & (j>0):
                R[i,j] = 1/4 * (image[i-1, j-1]  + image[i-1, j+1] + image[i+1, j-1]+image[i+1, j+1])
            if (i%2==0) & (j%2==0) & (i>0) & (i < c-1):
                R[i,j] = 0.5*(image[i, j-1] + image[i, j+1])
            if (i%2==1) & (j%2==1) & (i < c-1):
                R[i, j] = 0.5*(image[i-1, j] + image[i+1, j])
    R[:, 0] = R[:, 1]
    R[c-1, :] = R[c-2, :]
    return R

def B_channel(image, r, c):
    B = np.zeros((3000, 4000))
    B = image
    for i in range(c):
        for j in range(r):
            if (i%2==0) & (j%2==1) & (i>0) & (j<r-1):
                B[i,j] = 1/4 * (image[i-1, j-1]  + image[i-1, j+1] + image[i+1, j-1]+image[i+1, j+1])
            if (i%2==1) & (j%2==1) & (j<r-1):
                B[i,j] = 0.5*(image[i, j-1] + image[i, j+1])
            if (i%2==0) & (j%2==0) & (i >0) & (j<r-1):
                B[i, j] = 0.5*(image[i-1, j] + image[i+1, j])
    B[0, :] = B[1, :]
    B[:, c-1] = B[:, c-2]
    return B
 
img1 = read_raw(data_bytes) - 64
img2 = img1.reshape((3000, 4000))
gmap = G_channel(img2, 4000, 3000)

img1 = read_raw(data_bytes) - 64
img3 = img1.reshape((3000, 4000))
rmap = R_channel(img3, 4000, 3000)

img1 = read_raw(data_bytes) - 64
img4 = img1.reshape((3000, 4000))
bmap = B_channel(img4, 4000, 3000)
    
Img = np.zeros((3000, 4000, 3))

Yuv = rmap*0.299+gmap*0.587+bmap*0.114
Img[:, :, 0] = rmap
Img[:, :, 1] = gmap
Img[:, :, 2] = bmap
Img = Img.astype(int)
plt.imshow(Img)

对一张RAW10图进行Demosaic，得到下图：

可以思考一下为什么demosaic得到的图是绿色的呢？
答：这个camera sensor烧录的AWB数据有关。（灰色世界法、模组一致性）

AWB+Demosaic

所以接下来，利用sensor eeprom中烧录的AWB数据对Demosaic后的图像进行一个校正，
我拿到的EEPROM中的AWB烧录数据是：r/g=0.599609， b/g=0.677734
将这两个数据作用的R、B通道上：

a = 1/0.599609  #AWB OTPdata rg
b = 1/0.677734  #AWB OTPdata bg
Img2 = Img1
Img2[:,:,0] = Img1[:,:,0]*a
Img2[:,:,2] = Img1[:,:,2]*b
Img2 = Img2.astype(int)
plt.imshow(Img2)

最后得到的图如下所示，这样看来这是一张较为正常的图像啦。