图像超分辨率

作者

又是一篇简单的记录。图像超分辨率大致分为以下几类：单张LR->单张HR，这类基本依靠先验或者训练学习，实际运用上等同于猜。多张LR->多张HR，基本上是视频的超分辨率，多张LR->单张HR，这种就是典型的堆栈超分。堆栈超分是建立在亚像素配准和多帧融合的策略上，通过多帧图像采集不同的相位。亚像素配准又可以使用块匹配、光流、sift，融合策略五花八门，Google HDR做了不少这方面的工作，单帧超分的RAISR、多帧超分的Handheld Multi-frame SR等，可多多研究。

Raisr是典型的单帧超分，其将图像分成很多11×11的图像块进行训练，映射到特征空间（梯度、强度、相似度），生成离散特征空间的卷积核。在超分时，在11×11的图像块上，计算特征，映射的特征空间查找卷积核，对其进行超分。

多帧超分，可以通过sift特征匹配、光流、金字塔块匹配等计算亚像素，然后在亚像素水平上插值，做mask，然后stacking，进行多帧超分。

Hand-held Multi-frame super resolution实现记录

大致实现了Google Research 2019 paper中的内容，简单记录下。多帧超分的第一步是对齐，块匹配+光流的方式进行图像对齐，要求图像不能有太大的明暗变化。金字塔块匹配用了4层金字塔，块的大小用了16 32 和 64，这个根据噪声水平来选取。每一层的搜索范围设置了2pixel，没有设太大（感觉也可以把金字塔层数设成2层，把搜索范围改大）。光流用LK光流迭代就可以，要带上块匹配时候的初始值。

Motion Robustness计算的时候需要用到噪声模型，有些DNG图像里面存了噪声模型，直接拿来用就可以  N=SI+O 就是噪声方差，作为论文中的sigma_md，然后根据这个噪声方差来仿真得到dist_md。同时还可以根据噪声方差计算得到图像的SNR。

SNR=10log（I2/(SI+O)）。

在motion robustness计算里面还有一个高频拒绝要注意，对于一些实验室场景，可能出现密集的空间高频，这个时候就需要把高频区给找到（找的方法就是通过一个低通滤波器，然后看方差变化了多大），这些区域在最后stacking的时候不做融合。

Kernel Reconstruction实现的时候细节比较多，比如要根据snr做线性分段来选参数，算特征值的时候，如果要参数能和论文对上，就要把图像归一化到[0,1]，公式里面的 I 也是要区域求和的，不过最后实现出来有点儿问题，直接在RAW上做有很多artifacts，后来没有解决，并且跳出来kernel的形状感觉不太对。

简单写了一个Python的，Python运行实在太慢了，没写完，转投matlab搞，这个python的里面gamma之类的都没调，大致放个框架吧，这个里面还没有超分，要做超分的话增加采样频率就可以，这个里面robust mask没有算，也没有金字塔光流，只有一个简单的lk，姑且参考：

import os
import glob
import rawpy
import imageio
import cv2
from matplotlib import pyplot as plt
import numpy as np
import numpy.matlib
from scipy.interpolate import RectBivariateSpline
from scipy.interpolate import interp2d

result_dir = 'result/'
raw_dir = 'D:/data/bursts/33TJ_20150810_140637_256/'
raw_ext = '.DNG'
save_ext = '.bmp'

raw_files = sorted(glob.glob(raw_dir + '*' + raw_ext))

#process first frame
print('Processing ' + raw_files[0])
_, filename = os.path.split(raw_files[0])
filename, _ = os.path.splitext(filename)

#load raw image
with rawpy.imread(raw_files[0]) as raw:
    rgb = raw.postprocess(gamma=(1,1), no_auto_bright=True,
    output_bps=8)
    sz = np.shape(rgb)
    color = raw.raw_colors.copy()
    color = color[:sz[0],:sz[1]]
    bayer = np.empty([sz[0],sz[1]])
    bayer[color==0] = rgb[color==0,0]
    bayer[color==1] = rgb[color==1,1]
    bayer[color==2] = rgb[color==2,2]

cv2.imwrite(result_dir + filename + save_ext, rgb)

#bayer 2x2 downsample
bayer = bayer.reshape(sz[0]//2,2,sz[1]//2,2)
bayer = bayer.mean(axis=(1,3))
gray_ref = bayer.copy()
szh = np.shape(gray_ref)

#generate bayer pattern
cfa_pattern = np.zeros([sz[0], sz[1], 3], dtype=np.float32)
cfa_pattern[:,:,0] = np.matlib.repmat([[1,0],[0,0]],szh[0],szh[1])
cfa_pattern[:,:,1] = np.matlib.repmat([[0,1],[1,0]],szh[0],szh[1])
cfa_pattern[:,:,2] = np.matlib.repmat([[0,0],[0,1]],szh[0],szh[1])

#accumulate weight
accum_value = np.zeros([sz[0], sz[1], 3], dtype=np.float32)
accum_weight = np.zeros([sz[0], sz[1], 3], dtype=np.float32)
cfa_weight = cfa_pattern.copy()

accum_value = accum_value + rgb * cfa_weight
accum_weight = accum_weight + cfa_weight


for k in range(1,len(raw_files)):
    print('Processing ' + raw_files[k])
    _, filename = os.path.split(raw_files[k])
    filename, _ = os.path.splitext(filename)

    #load raw image
    with rawpy.imread(raw_files[k]) as raw:
        rgb = raw.postprocess(gamma=(1,1), no_auto_bright=True,
        output_bps=8)
        color = raw.raw_colors.copy()
        color = color[:np.shape(rgb)[0],:np.shape(rgb)[1]]
        #bayer = raw.raw_image.copy()
        sz = np.shape(rgb)
        bayer = np.empty([sz[0],sz[1]])
        bayer[color==0] = rgb[color==0,0]
        bayer[color==1] = rgb[color==1,1]
        bayer[color==2] = rgb[color==2,2]

    cv2.imwrite(result_dir + filename + save_ext, rgb)

    #bayer pattern 2x2 downsample
    bayer = bayer.reshape(sz[0]//2,2,sz[1]//2,2)
    bayer = bayer.mean(axis=(1,3))
    gray_cur = bayer

    #estimate optical flow
    pts = np.empty([np.size(gray_cur),1,2],dtype=np.float32)
    xv, yv = np.meshgrid(range(szh[1]), range(szh[0]))
    xv = xv.reshape(np.size(xv),order='C')
    yv = yv.reshape(np.size(yv),order='C')
    pts[:,0,0] = xv.astype(np.float32)
    pts[:,0,1] = yv.astype(np.float32)

    lk_params = dict(winSize = (15,15), maxLevel = 2,
        criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))
    p1, st, err = cv2.calcOpticalFlowPyrLK(gray_ref.astype(np.uint8), gray_cur.astype(np.uint8), pts, None, **lk_params)

    dstxh = p1[:,0,0]
    dstyh = p1[:,0,1]
    dstxh = dstxh.reshape(szh,order='C')
    dstyh = dstyh.reshape(szh,order='C')
    dstx_interp = RectBivariateSpline(range(szh[0]),range(szh[1]),dstxh,kx=1,ky=1)
    dstx = 2*dstx_interp(np.arange(0,szh[0],0.5),np.arange(0,szh[1],0.5))
    dsty_interp = RectBivariateSpline(range(szh[0]),range(szh[1]),dstyh,kx=1,ky=1)
    dsty = 2*dsty_interp(np.arange(0,szh[0],0.5),np.arange(0,szh[1],0.5))

    #warp grayscale image
    gray_cur_interp2 = RectBivariateSpline(range(szh[0]),range(szh[1]),gray_cur)
    gray_cur_warp = gray_cur_interp2(dstyh,dstxh,grid=False)
    cv2.imwrite(result_dir + filename + '_warp' + save_ext,gray_cur_warp.astype(np.uint8))

    #calculate robust mask

    #fusion
    for chns in range(3):
        rgb_interp2 = RectBivariateSpline(range(sz[0]),range(sz[1]),rgb[:,:,chns])
        rgb[:,:,chns] = rgb_interp2(dsty,dstx,grid=False)
        cfa_weight_interp2 = RectBivariateSpline(range(sz[0]),range(sz[1]),cfa_pattern[:,:,chns],kx=1,ky=1)
        cfa_weight[:,:,chns] = cfa_weight_interp2(dsty,dstx,grid=False)

    cv2.imwrite(result_dir + filename + '_warp_color' + save_ext,rgb.astype(np.uint8))

    accum_value = accum_value + rgb * cfa_weight
    accum_weight = accum_weight + cfa_weight

    #recover stacking image
    rgb_stack = accum_value / accum_weight
    cv2.imwrite(result_dir + filename + '_stack' + save_ext,rgb_stack.astype(np.uint8))
    print('Frame ' + filename + ' Done.')

#finish
#rgb_stack = accum_value / accum_weight
cv2.imwrite(result_dir + 'stack_rgb' + save_ext,rgb_stack.astype(np.uint8))

print(raw_files)

 

谢谢！

分享到微信

打开微信“扫一扫”，打开网页后点击屏幕右上角分享按钮
分享到微博

发表评论