github上找到的一个复现代码,这个只能参考,有些是错误的,我在文末会附上我们做的代码和视频 deblur

该论文的团队在2014年写了一篇报告,是对这篇论文的详细解释,对于理解这篇论文很有帮助,强烈建议看deconvmodelsolver_2014.pdf
这是我们报告的视频,最后有详细的讨论,视频中也出现了一些关键代码https://www.bilibili.com/video/BV1MK411p7cw


运动去模糊问题


论文讨论了在运动去模糊问题中,新的核估计方法对非盲卷积有着显著的意义。
作者发现了显著的边缘不总会提高核估计的效果,在确切的环境下还会使其边缘退化。


论文刚开始,提到了但是存在的两种去模糊的方法,一种叫稀疏先验,一种叫高斯先验,它们有自己的优势也有缺点。


稀疏先验
利用稀疏先验能够避免在迭代中陷于局部最优值。但稀疏先验会存在计算消耗过大的问题。


高斯先验
利用高斯先验能够缩短计算时间,但是高斯先验的问题在于会陷入噪声问题和对密集核估计不准确。


论文贡献


1,提出了一个新的2阶核估计算法,将计算消耗过大的非凸优化器从快速核初始化中分离,提出了一个更有效的核估计算法。
2,本文介绍了一种基于空间先验的能够保存潜在图像边缘信息的存储器。
3,通过迭代支持检测算法,可以自适应的加强空间约束核正确的保存参数。
4,最后通过tv-l1目标函数,使得核估计对噪声更健壮。


论文中相关原理介绍


理想的镜头成像时,一个对焦平面上的物点会投影为一个像点。
在这里插入图片描述
但事实上理想的镜头时不存在的,镜头总是存在一些缺陷会导致一个物点会投影为很多点。一个理想点经过相机后成的像由点扩散函数PSF(Point Spread Function)来描述。
而如果对此PSF做傅里叶变换,就可以得到这个镜头的光学传递函数OTF(Optical Transfer Function),它长这个样子:
在这里插入图片描述
假设有一个理想镜头,不受衍射的影响,用它所成的像为x, 而实际镜头的PSF为c, 实际镜头的成像是b,那么这三者之间的关系是一种典型的卷积关系:
在这里插入图片描述
已知模糊的图像b以及成像系统的PSF,恢复原始图像x的过程叫做非盲反卷积(Non-blind deconvolution),比如直接求卷积的维纳滤波(Wiener Filter)和通过迭代法求卷积的Richardson-Lucy方法。
在这里插入图片描述
非盲反卷积(Non-blind deconvolution)


非盲反卷积主要解决的问题一般包括以下几点:


1.减少出现在强边缘附近的振铃现象。
2.抑制噪音。
3.减少计算量。


而更多时候我们并不能提前测定PSF,而这种即在模糊核未知的情况下去卷积的方法叫做盲反卷积(Blind Deconvolution)。
在这里插入图片描述
对模糊核的估计我们就需要一些先验条件:


1.图像的梯度分布。清晰的自然图像的梯度符合一种叫做“Heavy-Tail”的分布形态。有很多平滑的区域且噪声较低的区域,其梯度也多为0。而图像清洗,其边缘会更明显,所以也有很多像素梯度较大的区域。
在这里插入图片描述
模糊的图像的边缘被糊掉了,所以更多的像素的梯度趋于0,因此其梯度直方图就会变化
在这里插入图片描述
2 - 模糊核的形态
对于运动模糊,其模糊核是稀疏的,有连续的轨迹,并且模糊核值都是非负数。
在这里插入图片描述


Kernel Estimation


接下来就是K的核估计
在这里插入图片描述
上图中需要注意:反卷积从最粗略的层开始,也就是说对图像金字塔进行迭代,先使用小尺寸的图片,最后使用大尺寸的图片。随着迭代的进行,慢慢会得到更准确的答案


在迭代(论文中提到的是盲去卷积过程中的迭代)中,使用冲击滤波器可以得到模糊图像的梯度边缘


接下来,通过梯度阈值去除小幅度边缘,只保留大边缘


在迭代(论文中提到的是盲去卷积过程中的迭代)中,使用冲击滤波器可以得到模糊图像的梯度边缘


用阈值处理后的边缘图替换原来的图片,替换后的图片用于下一次的PSF估计


在迭代过程初期,阈值边缘图很粗糙,和输入的模糊图很不一样


随着迭代的增加,越来越多的细节被添加到边缘图中,PSF估计也进一步完善


金字塔图像的作用
1,首先是构造n个等级的图像金字塔


2,接下来用每一个图像金字塔先得到粗糙的图像结果,再把这粗糙的结果将作为初始化进行下一步的循环迭代优化。一直重复,知道所有的图像金字塔层处理完毕。


3,在每一层的图像金字塔中,都需要几次迭代才能选择边和初始PSF


4,最终得到的PSF可以提高模糊图像细节处的恢复


边缘选择
1,在这个阶段,通过显著的梯度图的构建和核估计来估计PSF


2,为了使该过程更快,采用粗糙图像恢复来快速获得模糊图像的估计。


3,先对模糊图像进行高斯模糊,然后使用2014年那篇论文中提到的公式Eq. (2.56) 进行冲击滤波器处理


4,如果模糊核检测的对象的比例太小了,得出边缘信息将对PSF估计不利


由于空域的卷积等效于频域的乘法,因此我们只需要在频域上做除法,就能很好的恢复出k了,我们把这个过程称为Deconvolution,去卷积:
在这里插入图片描述
模糊核核初始化
通过多尺寸设置估计模糊核。通过高斯先验进行估计。主要分为三步:
1.边缘锐利重建;2.核估计;3.粗糙图像的存储;


基于ISD的核精炼


通过ISD的方法可以保证去模糊图片的质量并去除噪声。
本文的一个核心贡献就是发现就是,显著的边缘并不会提高核估计,相反如果目标尺寸比核还小,边缘信息可能会摧毁核估计。
在这里插入图片描述
具有与模糊信号相同的高度的虚线是去模糊过程中的最佳解决方案。如a中的红色虚线,b中的绿色虚线。A中的绿色虚线不能保留总的变化


1,a和b中的绿色虚线,通过高斯模糊,变成了图中的实线部分


2,由于a中的潜在信号的宽度较小,所以进行高斯模糊处理的时候,高度会降低,这也会误导PSF的估计。


3,关于a图对PSF的误导估计,具体来说,就是矮一点的信号和高一点的信号,两者作为模糊信号都具有相同的总变化,因此在拉普拉斯正则化中会产生较小的能量(在参考论文中出现的)


4,相反,对于b图,它有较大的尺寸,比核估计大,并且保留的边缘的总变化,所以不会产生a图那样的歧义,


5,a和b作对比,可以知道,在进行模糊处理中,如果模糊图像的显著结构发生了变化,则冲击滤波器产生的边缘可能会对PSF产生误导


6,接下来我们就要解决这个问题


核心创新点:
在这里插入图片描述
1.0.5是为了避免在平坦区域中产生较大的r。
2.在有纹理的区域,梯度处理后,有纹理的地方,会有正负两种值,假如上面为正,下面为负,在使用Nh矩阵卷积求和的时候,正负就会出现抵消。相反不是纹理的区域,就不会出现正负抵消。


1,选择有用的边缘信息,消除具有精细结构的纹理边缘区域
2,公式中 B表示模糊图像 Nh(i)是以像素i为中心的h×h窗口 0.5是为了避免在平坦区域中产生较大的r ∇B(j)是梯度图
3,对于只要包含纹理图案的窗口, ∇B会在很大程度上抵消
在这里插入图片描述
4,对3进行细致的讲解,在有纹理的地方,梯度处理后,有纹理的地方,会有正负两种值,假设上面为正,下面为负,在使用Nh矩阵卷积求和的时候,正负就会出现抵消,相反不是纹理的区域,就不会出现正负抵消。如下图,在a图就会出现正负抵消,b图就不会。这样在有纹理的地方
在这里插入图片描述
这样在有纹理的地方在这里插入图片描述就会是个很小的数字。


5,
在这里插入图片描述
||∇B(j)|| 表示吧梯度图求绝对值,那么梯度图中的负数就变成正数,那么在梯度图和Nh矩阵中进行卷积求和的时候,在纹理区域的求和就不会出现0或者很小的数字。梯度值越大,导致在这里插入图片描述越大,从而导致分母越大


通过Mask过滤较小的r值。
在这里插入图片描述
该公式的意思就是在梯度小于某个值得时候,就不计算该梯度。梯度也就表示边缘得强度。


这个公式的作用是,删除属于较小r值窗口的像素


H(·)是Heaviside阶跃函数,对于负值和零值输出1,否则输出零。
τr是一个阈值


在这里插入图片描述
这个公式的意思:


在这里插入图片描述
在这里插入图片描述


最后的和估计效果对比


1, ( c)-(e)和(g)-(i)分别表示不使用边缘检测和使用边缘检测计算出的Is图的差别
2,( c)-(f)中包含更多边, (g)-(i)包含更少的边,这说明,包含更多的边不一定有益于核估计,所以适当的边缘选择是很重要的
3,为了最终推断出细微的结构,可以在迭代中减小τr和τs的值,以包含越来越多的边缘。
在这里插入图片描述


快速核估计


在这里插入图片描述


1,使用选定的边缘图,使用上面公式可以完成PSF的初始化


2,由于有效果的梯度图在这里插入图片描述,f被限制在一个简单的二次项中


3,最小化E(K)会导致PSF产生更多噪音


在这里插入图片描述
1,上面的公式是由Parseval 定理和Eq. (2.39) 推到出来的
根据帕塞瓦尔定理将上式推导出下式。
1,上面的公式是由Parseval 定理和Eq. (2.39) 推到出来的
在这里插入图片描述


ISD-Based Kernel Refinement


为了获得稀疏的PSFs,前期的方法大都忽视了内部模糊核的结构,进而降低了核的质量。如下图b。只保留了较大值的元素,则显然不能保存模糊核更精细的结构。s1-s3就是通过ISD方法迭代之后的模糊核。
在这里插入图片描述
我们将ki中的所有元素按照升序排列,并计算相邻两者
间的差异。然后我们从头开始检查这些差异,从d0开始并找到第一
个元素。以dj为例,它需要满足这个公式,然后我们将位置j的内核
值分配给€s
在这里插入图片描述
在这里插入图片描述


在这里插入图片描述


因为需要模糊核是最稀疏的,因此我们可以通过最小化下式来求得模糊核k。
在这里插入图片描述


而最小化上式,我们可以采用迭代加权最小二乘法,即IRLS。该方法的核心是,通过L0-norm松弛化,把高度不连续的问题松弛化为连续的甚至是光滑的逼近。并且可以通过2范数来模拟0范数。
在这里插入图片描述
ISD算法如下图:
在这里插入图片描述


Fast TV-l1 Deconvolution


假设数据适合高斯分布,对大多数场景并不是一个好的方法。为了达到一个更好的鲁棒性,提出了一个TV-l1模型用于去卷积。


在这里插入图片描述
公式推导:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述


其他文献中方法类比:


在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述


振铃效应(edge taper):Matlab (EdgeTaper(I , kernel))


在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述


我的代码,有问题,但是关键的地方可以参考


import cv2
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt
import math
from PIL import Image

np.set_printoptions(threshold=np.inf)
#公式1
‘’’
参考的gthub和博客
https://github.com/qq456cvb/deblur/blob/master/main.py
https://github.com/apvijay/shock_filter/blob/master/shock_filter.m
https://blog.csdn.net/kastolo/article/details/10512241
‘’’

#shock filter image

def shock_filter(img,dt = 0.25):
img = cv2.blur(img, ksize=(3, 3))
#(5, 5)表示高斯矩阵的长与宽都是5,标准差取0
img = cv2.GaussianBlur(img,(5,5),0)
for i in range(4):
Ix = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=3)
Iy = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=3)
Ixx = cv2.Sobel(Ix, cv2.CV_64F, 1, 0, ksize=3)
Ixy = cv2.Sobel(Ix, cv2.CV_64F, 0, 1, ksize=3)
Iyy = cv2.Sobel(Iy, cv2.CV_64F, 0, 1, ksize=3)
DeltaI = Ix _ Ix _ Ixx + 2 _ Ix _ Iy _ Ixy + Iy _ Iy _ Iyy
gradI = np.linalg.norm(np.stack([Ix, Iy], -1), axis=-1)

IT = -np.sign(DeltaI)_gradI
img = img + (IT) _ dt;
return img

#公式2
def rEq2(img,h):
Bx = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=3)
By = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=3)

plt.imshow(Bx)
plt.title(“sobelX”)
plt.show()

plt.imshow(By)
plt.title(“sobelY”)
plt.show()

#定义窗口Nh 大小暂时定为3x3
Nh = np.ones((h,h));
h = math.ceil(h/2)

convX = signal.convolve2d(Bx, Nh,mode=‘same’)
convY = signal.convolve2d(By, Nh,mode=‘same’)

convXY = np.stack([convX, convY], -1)

#axis=0表示按列向量处理,求多个列向量的范数
up = np.linalg.norm(convXY, axis=-1)


Bxy = np.stack([Bx, By], -1)
down = signal.convolve2d(np.linalg.norm(Bxy, axis=-1), Nh,mode=‘same’) + 0.5

r = up/down

plt.imshow(r)
plt.title(“r”)
plt.show()


return r

#公式3
def Meq3(r,tau_r):
M = np.heaviside(r - tau_r, 0)
return M

#公式4
def edges(img,r,tau_r,tau_s):
L1 = shock_filter(img)
L2 = np.linalg.norm(L1,ord=2)
M = Meq3(r,tau_r)
L3 = M_L2 - tau_s
L4 = np.heaviside(L3,0)
Is = L1_L4

return Is

#公式6
def kEq6(B,Is,gamma):

Isx = cv2.Sobel(Is, cv2.CV_64F, 1, 0, ksize=3)
Isy = cv2.Sobel(Is, cv2.CV_64F, 0, 1, ksize=3)
IsxCon = np.conjugate(Isx)
IsyCon = np.conjugate(Isy)

Bx = cv2.Sobel(B, cv2.CV_64F, 1, 0, ksize=3)
By = cv2.Sobel(B, cv2.CV_64F, 0, 1, ksize=3)
fBx = np.fft.fft2(Bx)
fBy = np.fft.fft2(By)

fIsx = np.fft.fft2(Isx)
fIsy = np.fft.fft2(Isy)

up = IsxCon_fBx + IsyCon_fBy

down = fIsx_fIsx + fIsy_fIsy + gamma

k = np.fft.ifft2(up/down)

return k


#公式8

‘’’
python 自己编写实现Sobel算子
https://blog.csdn.net/GuaPiQ/article/details/99890040
‘’’


def Ieq8(B,k,LAMBDA,Is):

scharr_x, scharr_y = np.zeros_like(k), np.zeros_like(k)

scharr_x[0, 1] = 1
scharr_x[0, -1] = -1

scharr_y[1, 0] = 1
scharr_y[-1, 0] = -1

scharr_x /= 2.
scharr_y /= 2.

fftdx = np.fft.fft2(scharr_x)
fftdy = np.fft.fft2(scharr_y)

Isx = cv2.Sobel(Is, cv2.CV_64F, 1, 0, ksize=3)
Isy = cv2.Sobel(Is, cv2.CV_64F, 0, 1, ksize=3)

l1 = np.conjugate(k)_np.fft.fft2(B)
l2 = np.conjugate(fftdx)_np.fft.fft2(Isx)
l3 = np.conjugate(fftdy)_np.fft.fft2(Isy)
l4 = LAMBDA_(l2+l3)

up = l1+l4

l5 = np.conjugate(k)_np.fft.fft2(k)
l6 = np.conjugate(fftdx)_np.fft.fft2(fftdx)
l7 = np.conjugate(fftdy)_np.fft.fft2(fftdy)
l8 = LAMBDA_(l6+l7)

down = l5+l8

I = np.fft.ifft2(up/down)


return I

image = cv2.imread(‘2.png’)
gray_img=cv2.cvtColor(image,cv2.COLOR_BGR2GRAY)


# 核的大小设置
ksize = 25;

# 定义全为0的核
kernel = np.zeros((ksize,ksize));

# n代表图像金字塔的高度或者有多少层
n = 3

# 定义窗口Nh的大小
h = 18

tau_r = 0.1
tau_s = 0.05

# 在粗核估计阶段,设置λ= 2e-3和γ= 10来抵抗噪声。 在核细化中,将γ= 1设置为1。将最终图像去卷积中的λ设置为2e-2。
gamma = 10

LAMBDA = 2_math.exp(-3)

for i in range(n):

gray_img = cv2.pyrDown(gray_img) # 上采样操作
plt.imshow(gray_img)
plt.title(“pyrDown”)
plt.show()

#公式2
r = rEq2(gray_img,h)


for it in range(10):
#公式4
Edge = edges(gray_img,r,tau_r,tau_s)
plt.imshow(Edge)
plt.title(“edges”)
plt.show()

#公式6
k = np.real(kEq6(gray_img,Edge,gamma))
# k = np.maximum(k, 0)

plt.imshow(k)
plt.title(“K”)
plt.show()

#公式8
I = np.real(Ieq8(gray_img,k,LAMBDA,Edge))
# print(I)
plt.imshow(I)
plt.title(“I”)
plt.show()

#为了在内核优化过程中推断出微妙的结构,我们在迭代中减小了τr和τs的值(每遍均除以1.1)
tau_s /= 1.1
tau_r /= 1.1

gray_img = I








下面是工作刚开始的时候写的总结,也可以参看参考




两阶段稀疏核估计


模糊过程建模为:
B=I ⊗ k + ε,
其中 I 是潜像,也就是 不是模糊状态的原图,k是模糊核,ε是图像噪声,⊗表示卷积,B是观察到的模糊图像


_论文中有许多关于PSF(点扩展估计函数)描述,PSF的介绍链接_


这一部分会介绍一种用于PSF估计的两阶段方法
第一阶段:计算内核的粗略版本,其中无需执行太多稀疏性。
第二阶段:采用了非凸优化,从阶段1传播了初始内核估计,不需要大量计算即可得出最终结果


以上两阶段,现在看不懂没关系,这是对整个方法的概括,需要看完后面的,再返回来看,才能明白


1 第一阶段:内核初始化


第一步,在多尺度设置下估计模糊核。使用存在封闭形式解的高斯先验时,可以产生高效率。
高斯先验的相关资料


算法:
三个主要步骤:1,锋利的边缘构建,2,核估计,3,粗略图像恢复。




算法1:核初始化




输入:模糊图像B和全0核(大小为hxh)
用等级索引{1,2,…,n}来构建一个图像金字塔
for




l



l


l
= 1 to n do
      计算所有像素的梯度置信度r(等式2)
      for i = 1 to m (m是迭代次数)
            (a)根据置信度r(等式4)选择边缘








I


s




∇I^{s}


Is

            (b)用高斯先验估计核(等式6)
            (c)用空间先验(等式8)估计潜像





I


l




I^{l}


Il
,并更新
      end for
      高档图像(这里应该指的是恢复的好的图)





I



l


+


1








I


l







I^{l+1} ← I^{l}↑


Il+1Il

end for
输出:核估计





k


0




k^{0}


k0
和尖锐边缘梯度图








I


s




∇I^{s}


Is




这个算法中的一些相关资料:
高斯金字塔
高斯金字塔与拉普拉斯金字塔
python实现opencv学习十五:高斯金字塔和拉普拉斯金字塔


在这里插入图片描述


图1,运动模糊中去模糊


图中的两个绿色虚线表示:两个潜在信号。
图中两个蓝色线表示:绿色虚线(潜在信号)使用相同的高斯核后的曲线(变模糊)
在(a)中,模糊信号没有保留全部变化,使得核估计不明确。
实际上,在常见的优化过程中,红色曲线比绿色曲线更可能是潜在信号。
底部的橙色线指示输入的核宽度。


通过滤波对图像进行预平滑处理,来解决震荡的滤波PDE问题(这里没懂,有相关论文可以参考, Osher, S., Rudin, L.: Feature-oriented image enhancement using shock filters.
SIAM Journal on Numerical Analysis 27, 919–940 (1990)










I


/





t


=





s


i


g


n


(


Δ


I


)











I








,



I


0



=



G


σ






I


i


n


p


u


t


,


(


1


)



∂I/∂t = −sign(ΔI)||∇I||, I_0 = G_σ ⊗ Iinput, (1)


I/t=sign(ΔI)I,I0=GσIinput,(1)

其中







I


=


(


I


x


,


I


y



)












Δ


I


=



I


x


2




I



x


x




+


2



I


x




I


y




I



x


y




+



I


y


2




I



y


y





∇I=(Ix,Iy)^{‘},ΔI= I^{2}_ xI_{xx} + 2I_xI_yI_{xy} + I^{2} _yI_{yy}


I=(Ix,Iy)ΔI=Ix2Ixx+2IxIyIxy+Iy2Iyy
分别是一阶和二阶空间导数。





I


0




I_0


I0
表示高斯平滑输入图像,它用作迭代更新







I


/





t



∂I /∂t


I/t
的初始输入


公式1的代码:


def shock_filter(img):
#(3, 3)表示高斯矩阵的长与宽都是3,标准差取0
img = cv2.GaussianBlur(img,(3,3),0)

#it表示迭代次数,这里表示迭代10次
it = 10
dt = 0.01
#迭代10次
for i in range(it):
Ix = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=3)
Iy = cv2.Sobel(img, cv2.CV_64F, 0, 1, ksize=3)
Ixx = cv2.Sobel(Ix, cv2.CV_64F, 1, 0, ksize=3)
Ixy = cv2.Sobel(Ix, cv2.CV_64F, 0, 1, ksize=3)
Iyy = cv2.Sobel(Iy, cv2.CV_64F, 0, 1, ksize=3)
img -= dt _ np.sign(Ix _ Ix _ Ixx + 2 _ Ix _ Iy _ Ixy + Iy _ Iy _ Iyy) * np.linalg.norm(np.stack([Ix, Iy], -1), axis=-1)
plt.imshow(img)
plt.show()
return img

关于_cv2.GaussianBlur_,这里查了一些资料,高斯模糊和高斯平滑,网上的资料讲的比较模糊,没有说明两者的区别,大多使用GaussianBlur实现
我这里也就使用GaussianBlur来实现文中所说的高斯平滑
github上使用的cv2.blur,我这里就修改为 cv2.GaussianBlur
相关参考opencv高斯滤波cv2.GaussianBlur


代码中也出现了sobel的用法,网上有篇博客也写得不错,可以参考:
Python+opencv利用sobel进行边缘检测(细节讲解)


sign(x)函数其实就是一个分段函数,x>0,取值为1;x=0,取值为0;x<0,取值为-1.
在这里插入图片描述
相关资料:【Python】Numpy库之符号函数sign()的介绍及用法


在公式中,













I









||∇I||


I
代表求







I



∇I


I
的范数,其对应代码中为:np.linalg.norm
参考资料:np.linalg.norm(求范数)










I


=


(


I


x


,


I


y



)










∇I=(Ix,Iy)^{‘}


I=(Ix,Iy)
在代码中为:np.stack([Ix, Iy], -1) ,stack表示数组的堆叠,axis=-1也就是在最后一维后增加一维,numpy中stack()的方法的作用
但是,这个代码是否完成了公式的运算,公式代表什么,代码表示出了公式么?我还没弄懂
同样 Ix = cv2.Sobel(img, cv2.CV_64F, 1, 0, ksize=3) 是否复现了公式





I


x




I_x


Ix
我也不知道


在这里插入图片描述
再来深入分析一下这张图




在模糊图像中,我们看到的模糊信号就是以蓝色的线显示
由于左信号在水平方向上较窄,因此在(a)中模糊处理会降低其高度,从而在潜在信号恢复中产生歧义。


具体来说,将运动稀疏先验强加于潜像的梯度图 上的运动模糊方法将在计算未模糊信号时偏向于红色虚线,因为该版本的梯度幅值较小。(_这里的梯度图就是指算法中的图像金字塔的图_)


此外,红色信号比绿色信号更好地保留了总变化。


所以,对于使用锐边预测的方法组,是一种更合适的解决方案。


此示例表明,如果图像结构的大小在模糊后发生显着变化,则相应的边缘信息可能会误认为核估计。


相比之下,(b)所示的较大尺度的对象可以产生稳定的核估计,因为它比核宽,可以保留潜在信号沿其边缘的总变化。




上面的这几句话估计在看的你也没懂,其实我也模模糊糊的,这是直接翻译论文中的原话。
在这里插入图片描述
图2.图像结构对核估计的影响。
(a)图像模糊。(b)Fergus等人的结果。 [1],(c)Shan等人的结果。 [3]。(d)r映射(由等式(2))。(e)-(g)PoI s映射,使用Poisson重建在第1次,第2次和第7次迭代中可视化,而不考虑r。(h)不使用r映射对结果进行模糊处理。(i)-(k)∇Is映射根据等式计算。 (4)。(l)我们的最终结果。 模糊PSF的尺寸为45×45。


图2中,模糊的输入(如(a)所示)包含沿许多小比例对象的丰富边缘信息。(c)是通过大量手动调整参数来计算的。 但是,仍然由于主要的上述小结构问题而仍然找不到正确的核估计。


我们提出了一种新的准则,用于选择用于核估计的信息性边缘。
测量梯度有用性的新指标定义为:
在这里插入图片描述
B表示模糊图像,





N



h


(


x


)





N_{h(x)}


Nh(x)
是以像素x为中心的




h


×


h



h×h


h×h
窗口。0.5是为了防止在平坦区域产生较大的r。
狭窄物体(尖峰)的带符号







B


(


y


)



∇B(y)


B(y)
大部分会在







y






N



h


(


x


)







B


(


y


)



\sum y∈N_{h(x)}∇B(y)


yNh(x)B(y)
中抵消。
小r意味着涉及尖峰或平坦区域,这会导致中和许多梯度分量。
以上的三句官方论文解释,我没太懂,只能翻译过来。


先上公式2的代码:


ksize = 25
image = cv2.imread(‘2.PNG’,)
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY).astype(np.float32) / 255.

# 这里金字塔设置为3层
n_pyramids = 3
pyramids = [image]
ksizes = [ksize]

#这里的for循环没有执行,n_pyramids - 1为0
for i in range(n_pyramids - 1):
pyramids.append(cv2.pyrDown(pyramids[-1]))
ksizes.append(ksize // 2)
if ksizes[-1] % 2 == 0:
ksizes[-1] += 1
plt.imshow(pyramids[-1])
plt.show()

tau_r = 0.1
tau_s = 0.05

for i in range(n_pyramids):
print(“this is “ + str(i) + “ pyramids”)
#一个个读取img,图片由大到小,从图片金字塔底部读取
img = pyramids[i]

n_x = np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
n_y = n_x.transpose()
Bx = signal.convolve2d(img, n_x , mode=‘same’)
By = signal.convolve2d(img, n_y, mode=‘same’)
convx = signal.convolve2d(Bx, np.ones((ksizes[i], ksizes[i])), mode=‘same’)
convy = signal.convolve2d(By, np.ones((ksizes[i], ksizes[i])), mode=‘same’)

#r为梯度置信度
r = np.linalg.norm(np.stack([convx, convy], -1), axis=-1) / \
(signal.convolve2d(np.linalg.norm(np.stack([Bx, By], -1), axis=-1), np.ones((ksizes[i], ksizes[i])),
mode=‘same’) + 0.5 / 255.)
plt.imshow(r)
plt.show()

首先是算法1 中提到的图像金字塔,n_pyramids = 3表示把图像金字塔分为3层,底层为原图,越往上图片像素越低,图片尺寸越小。


经过3次循环,使用cv2.pyrDown(pyramids[-1]),将图片一次次变小,然后将3张图像放入pyramids数组中,形成图像金字塔


这里设置的核大小为25,ksize=25


接下来看这个公式的代码实现:








y







N



h


(


x


)







B


(


y


)



\sum_y∈N_{h(x)} ∇B(y)


yNh(x)B(y)







N



h


(


x


)





N_{h(x)}


Nh(x)
窗口是一个




h


×


h



h×h


h×h
矩阵
本文定义为 n_x = np.array([[-1,0,1],[-2,0,2],[-1,0,1]]) 为什么会使用这个矩阵,我也不知道,网上找到的,在github的代码里使用的是1x3的矩阵,论文中要求使用




h





h



h_h


hh
的矩阵,所以这里我做了修改.


signal.convolve2d(img, n_x , mode=‘same’)表示:图像平滑二维离散卷积,参数same表示卷积结果和原图像的高、宽相等
参考资料:OpenCV–Python 图像平滑之二维离散卷积


经过运算后,运行结果:
图像金字塔的图:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
我是使用 jupyter notebook 运行的,虽然大小看起来一样,但是看它左边核下边的尺度,会发现越来越小,图像也越来越模糊。


然后就是用这三张图做循环,得到公式2的结果:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述


结果就是这样的,我们在和原论文做一下对比
在这里插入图片描述
感觉上差不多,颜色上的原因是jupyter notebook的原因。


然后,我们使用遮罩排除属于较小r值窗口的像素





M


=


H


(


r






τ


r



)


,


(


3


)



M = H(r − τ_r), (3)


M=H(rτr),(3)

这个公式实现起来很简单,一行代码就搞定了


np.heaviside(r - tau_r, 0)
  • 1







τ


r




τ_r


τr
是一个阈值
np.heaviside是单位阶跃函数,heaviside和上文写的sign(x)有相似的地方
在这里插入图片描述
在这里插入图片描述
参考资料:heaviside()


下一部分,用于核估计的最终选择的边缘被确定为:
在这里插入图片描述
其中,







I



〜I


I
(这里不知道怎么打出




I



I


I
上面一个〜)表示经冲击滤波的图像,并且





τ


s




τ_s


τs
是梯度幅度的阈值。
先看看代码吧,论文读到这里,我都不知道在讲什么了


    for it in range(10):
print(“this is “ + str(it+1) + “ iteration”)
Ix = signal.convolve2d(I, np.array(n_x), mode=‘same’)
Iy = signal.convolve2d(I, np.array(n_y), mode=‘same’)
shocked = shock_filter(I)

Ix_shocked = signal.convolve2d(shocked, np.array([[-1, 0, 1]]) / 2, mode=‘same’)
Iy_shocked = signal.convolve2d(shocked, np.array([[-1], [0], [1]]) / 2, mode=‘same’)

print(“Eq(4):”)
edges = np.stack([Ix_shocked, Iy_shocked], -1) _ np.heaviside(np.heaviside(r - tau_r, 0) * np.linalg.norm(np.stack([Ix_shocked, Iy_shocked], -1), axis=-1) - tau_s, 0)[..., None]
plt.imshow( edges[..., 0])
plt.show()

    这里的迭代此时定为10次,运行结果为:
    在这里插入图片描述
    在这里插入图片描述
    第一张图是第一次循环的结果,第二章图是最后一次循环的结果,做到这里,我就感觉不对了,和论文中的结果对不上了。


    我还做了公式6和公式8的结果,代码就不看了,直接上结果,结果应该也有问题
    在这里插入图片描述
    在这里插入图片描述
    公式6:第一张图是第一次循环的结果,第二章图是最后一次循环的结果


    在这里插入图片描述
    在这里插入图片描述
    公式8:第一张图是第一次循环的结果,第二章图是最后一次循环的结果


    最终的结论是:复现失败,望后人继续努力


     
    转载自:https://blog.csdn.net/WhiffeYF/article/details/105737964?spm=1001.2014.3001.5502