计算机视觉中的变分方法-扩散(Diffusion)
最近在看一个计算机视觉中的变分方法系列的视频,是德国慕尼黑工大出的,讲课老师是LSD-SLAM的作者Daniel Cremers,老师讲得很清楚,看了还是很有收获的。我已经变成Cremers大神的脑残粉了,有兴趣看视频的戳这里Variational Methods for Computer Vision
Diffusion equation:
扩散是一种物理过程,是让空间中的物质的浓度分布u(x,t)u(x,t)更加均匀一些。这个过程可以用两个基础的等式来描述:
Fick′slawFick′slaw : 空间物质的浓度的差别导致在浓度的负梯度方向上会有流jj。 这个也很好理解,意思就是说浓度高出的物质会往浓度低处扩散:
j=−g∇u
j=−g∇u
其中, gg是扩散系数(diffusivity),表示扩散过程的快慢
continuitycontinuity equationequation :
∂u∂t=−divj
∂u∂t=−divj
这里,divj=∇⋅j=∂jx∂x+∂jy∂ydivj=∇·j=∂jx∂x+∂jy∂y 称为散度。关于散度,其实在高等数学中有过介绍,通俗来讲,对于空间场中一点,如果该点散度大于00,则表示该点向外扩散物质(好比是该点是水龙头,向外流水);如果该点散度等于00, 那就是扩散保持平衡,进多少出多少;如果散度小于00,那么就说明该点在吸收物质(就像黑洞一样吸收空间场中该点附近的物质)。
关于散度的更多资料,可以参见知乎上这个回答在图像处理中,散度 div 具体的作用是什么
由上面两个基本的等式,联合起来就得到了今天要讲的扩散方程(DiffusionDiffusion equationequation)
∂u∂t=div(g∇u)
∂u∂t=div(g∇u)
Example: Linear Diffusion Equation
下面以一维线性扩散方程为例来说明。
对于线性情况,g=1g=1:
∂tu=∂2tu
∂tu=∂t2u
初始条件:
u(x,t=0)=f(x)
u(x,t=0)=f(x)
这个方程有唯一解:
u(x,t)=(G2t√)∗f(x)=∫G2t√(x−x′)f(x′)dx′
u(x,t)=(G2t)∗f(x)=∫G2t(x−x′)f(x′)dx′
其中,Gσ=12πσ√exp(−x22σ2)Gσ=12πσexp(−x22σ2),是一个高斯核,σ=2t−−√σ=2t
可以看到,高斯滤波其实是扩散的一种特例。但是我们都知道,用高斯滤波器对一个图像进行平滑滤波,由于高斯滤波器的各向同性,会使图像的边缘细节都变模糊,有时候这不是我们想要的结果。
Nonlinear and Anisotropic Diffusion
一般形式下的扩散方程:
∂tu=div(g∇u)
∂tu=div(g∇u)
对图像滤波时,要想保持图像的边缘细节,可以在图像边缘信息强的地方少扩散一些,那么怎么做呢?
我们用梯度的模来作为检测边缘的算子|∇u|=u2x+u2y−−−−−−√|∇u|=ux2+uy2,那么在边缘处|∇u||∇u|的值就会比较大 ,然后再这些地方让扩散速率变小,可以构造这样的 gg:
g(|∇u|)=11+|∇u|2/λ2−−−−−−−−−−−√
g(|∇u|)=11+|∇u|2/λ2
其中,λ>0λ>0,称为对比参数,在|∇u|>>λ|∇u|>>λ的区域,扩散速度接近于00,不受扩散的影响,所以可以保持该区域的细节。
关于这部分的详细资料,可以参考图像处理的经典论文11Scale-space and edge detection using anisotropic diffusion
有限差分的实现
上面讲了各向异性的扩散方程,现在就来说明一下如何编程实现。这部分内参考的是 Weickert, J Anisotropic diffusion in image processing里的内容。
非线性扩散方程:
∂tu=∂x(g|∇u|∂xu)+∂y(g|∇u|∂yu)
∂tu=∂x(g|∇u|∂xu)+∂y(g|∇u|∂yu)
用差分来代替微分:∂tu=ut+1ij−utijτ∂tu=uijt+1−uijtτ
非线性扩散方程右边第一项就可以表示为:
∂x(g∂xu)=((g∂xu)ti+1/2,j−(g∂xu)ti−1/2,j)∂x(g∂xu)=((g∂xu)i+1/2,jt−(g∂xu)i−1/2,jt)
=(gti+1/2,j(uti+1,j−uti,j)−gti−1/2,j(uti,j−uti−1,j))=(gi+1/2,jt(ui+1,jt−ui,jt)−gi−1/2,jt(ui,jt−ui−1,jt))
其中,gti+1/2,j=gi+1,jgi,j−−−−−−−√gi+1/2,jt=gi+1,jgi,j,说明只要这两个像素点处有一个的扩散速率gg为00,那么插值得到的gti+1/2,jgi+1/2,jt就会为00,而不是去这两者的平均值。
关于这段差分实现的公式部分,需要说明的是扩散方程中对x,yx,y是进行了二阶差分,注意在上面公式中,第一次对xx方向差分选择的两个点是(i,j)(i,j)旁边的两个点(i+1/2,j)和(i−1/2,j)(i+1/2,j)和(i−1/2,j),然后又进行了一次差分,得到的结果中,用到的像素点位置只有(i,j),(i+1,j),(i−1,j)(i,j),(i+1,j),(i−1,j),这样还是在一个3x33x3的窗口操作的,如果按照以前的第一次差分是右边的(i+1,j)(i+1,j)减左边的(i−1,j)(i−1,j),那么结果就会出现(i+2,j),(i−2,j)(i+2,j),(i−2,j)项,最后就是相当于用了5x55x5的窗口,大的窗口对于细节的保持是不利的。
Anisotropic Diffusion Matlab代码示例关于代码实现的这部分内容,可以进一步参考这里。
使用示例:
function diff_im = anisodiff2D(im, num_iter, delta_t, kappa, option)
%ANISODIFF2D Conventional anisotropic diffusion
% DIFF_IM = ANISODIFF2D(IM, NUM_ITER, DELTA_T, KAPPA, OPTION) perfoms
% conventional anisotropic diffusion (Perona & Malik) upon a gray scale
% image. A 2D network structure of 8 neighboring nodes is considered for
% diffusion conduction.
%
% ARGUMENT DESCRIPTION:
% IM - gray scale image (MxN).
% NUM_ITER - number of iterations.
% DELTA_T - integration constant (0 <= delta_t <= 1/7).
% Usually, due to numerical stability this
% parameter is set to its maximum value.
% KAPPA - gradient modulus threshold that controls the conduction.
% OPTION - conduction coefficient functions proposed by Perona & Malik:
% 1 - c(x,y,t) = exp(-(nablaI/kappa).^2),
% privileges high-contrast edges over low-contrast ones.
% 2 - c(x,y,t) = 1./(1 + (nablaI/kappa).^2),
% privileges wide regions over smaller ones.
%
% OUTPUT DESCRIPTION:
% DIFF_IM - (diffused) image with the largest scale-space parameter.
%
% Example
% -------------
% s = phantom(512) + randn(512);
% num_iter = 15;
% delta_t = 1/7;
% kappa = 30;
% option = 2;
% ad = anisodiff2D(s,num_iter,delta_t,kappa,option);
% figure, subplot 121, imshow(s,[]), subplot 122, imshow(ad,[])
%
% See also anisodiff1D, anisodiff3D.
% References:
% P. Perona and J. Malik.
% Scale-Space and Edge Detection Using Anisotropic Diffusion.
% IEEE Transactions on Pattern Analysis and Machine Intelligence,
% 12(7):629-639, July 1990.
%
% G. Grieg, O. Kubler, R. Kikinis, and F. A. Jolesz.
% Nonlinear Anisotropic Filtering of MRI Data.
% IEEE Transactions on Medical Imaging,
% 11(2):221-232, June 1992.
%
% MATLAB implementation based on Peter Kovesi's anisodiff(.):
% P. D. Kovesi. MATLAB and Octave Functions for Computer Vision and Image Processing.
% School of Computer Science & Software Engineering,
% The University of Western Australia. Available from:
% <http://www.csse.uwa.edu.au/~pk/research/matlabfns/>.
%
% Credits:
% Daniel Simoes Lopes
% ICIST
% Instituto Superior Tecnico - Universidade Tecnica de Lisboa
% danlopes (at) civil ist utl pt
% http://www.civil.ist.utl.pt/~danlopes
%
% May 2007 original version.
% Convert input image to double.
im = double(im);
% PDE (partial differential equation) initial condition.
diff_im = im;
% Center pixel distances.
dx = 1;
dy = 1;
dd = sqrt(2);
% 2D convolution masks - finite differences.
hN = [0 1 0; 0 -1 0; 0 0 0];
hS = [0 0 0; 0 -1 0; 0 1 0];
hE = [0 0 0; 0 -1 1; 0 0 0];
hW = [0 0 0; 1 -1 0; 0 0 0];
hNE = [0 0 1; 0 -1 0; 0 0 0];
hSE = [0 0 0; 0 -1 0; 0 0 1];
hSW = [0 0 0; 0 -1 0; 1 0 0];
hNW = [1 0 0; 0 -1 0; 0 0 0];
% Anisotropic diffusion.
for t = 1:num_iter
% Finite differences. [imfilter(.,.,'conv') can be replaced by conv2(.,.,'same')]
nablaN = imfilter(diff_im,hN,'conv');
nablaS = imfilter(diff_im,hS,'conv');
nablaW = imfilter(diff_im,hW,'conv');
nablaE = imfilter(diff_im,hE,'conv');
nablaNE = imfilter(diff_im,hNE,'conv');
nablaSE = imfilter(diff_im,hSE,'conv');
nablaSW = imfilter(diff_im,hSW,'conv');
nablaNW = imfilter(diff_im,hNW,'conv');
% Diffusion function.
if option == 1
cN = exp(-(nablaN/kappa).^2);
cS = exp(-(nablaS/kappa).^2);
cW = exp(-(nablaW/kappa).^2);
cE = exp(-(nablaE/kappa).^2);
cNE = exp(-(nablaNE/kappa).^2);
cSE = exp(-(nablaSE/kappa).^2);
cSW = exp(-(nablaSW/kappa).^2);
cNW = exp(-(nablaNW/kappa).^2);
elseif option == 2
cN = 1./(1 + (nablaN/kappa).^2);
cS = 1./(1 + (nablaS/kappa).^2);
cW = 1./(1 + (nablaW/kappa).^2);
cE = 1./(1 + (nablaE/kappa).^2);
cNE = 1./(1 + (nablaNE/kappa).^2);
cSE = 1./(1 + (nablaSE/kappa).^2);
cSW = 1./(1 + (nablaSW/kappa).^2);
cNW = 1./(1 + (nablaNW/kappa).^2);
end
% Discrete PDE solution.
diff_im = diff_im + ...
delta_t*(...
(1/(dy^2))*cN.*nablaN + (1/(dy^2))*cS.*nablaS + ...
(1/(dx^2))*cW.*nablaW + (1/(dx^2))*cE.*nablaE + ...
(1/(dd^2))*cNE.*nablaNE + (1/(dd^2))*cSE.*nablaSE + ...
(1/(dd^2))*cSW.*nablaSW + (1/(dd^2))*cNW.*nablaNW );
% Iteration warning.
fprintf('\rIteration %d\n',t);
end
---------------------
作者:蜗牛一步一步往上爬
来源:CSDN
原文:https://blog.csdn.net/yc461515457/article/details/50847526
版权声明:本文为博主原创文章,转载请附上博文链接!
左边是原图,右边是Anisotropic Diffusion结果图
计算机视觉中的变分方法-扩散(Diffusion)相关推荐
- MICCAI 2022中的医学扩散模型
来源:知乎-JunMa 地址:https://zhuanlan.zhihu.com/p/569555913 没有意外,diffusion model的热浪也涌入了医学图像,这是今年MICCAI看到的几 ...
- 图像处理和计算机视觉中的经典论文
转自:http://www.cnblogs.com/moondark/archive/2012/04/20/2459594.html 感谢水木上同领域的同学分享,有了他的整理,让我很方便的获得了CV方 ...
- 图像处理和计算机视觉中的经典论文(部分)
自己视野狭小,不敢说全部,只是把自己熟悉的方向中的部分经典文章列出来了.经典的论文,读得怎么透都不过分.有人说关于配准的文章太多了,其实我也不太关注这方面,不过由于它们引用率都比较高,就都列出来了,不 ...
- Talk | 阿姆斯特丹大学博士生胡涛:计算机视觉中的标签效率学习
本期为TechBeat人工智能社区第509期线上Talk! 北京时间6月29日(周四)20:00,阿姆斯特丹大学博士生-胡涛的Talk将准时在TechBeat人工智能社区开播! 他与大家分享的主题是: ...
- 深度学习在计算机视觉中的应用长篇综述
深度学习在计算机视觉中的应用长篇综述 前言 2012年ImageNet比赛,使深度学习在计算机视觉领域在全世界名声大震,由此人工智能的全球大爆发.第一个研究CNN的专家使Yann LeCun,现就职于 ...
- 使用Python,EoN模拟网络中的疾病扩散模型,并结合matplotlib绘图
使用Python,EoN模拟网络中的疾病扩散模型,并结合matplotlib绘图 1. EoN是什么 2. 安装 3. 效果图 4. 源代码 4.1 源码 4.2 源码 参考 写这篇博客源于博友的提问 ...
- 如何通过深度学习,完成计算机视觉中的所有工作?
Mask-RCNN做对象检测和实例分割 作者 | George Seif 译者 | 天道酬勤,责编 | Carol 出品 | AI科技大本营(ID:rgznai100) Mask-RCNN做对象检测和 ...
- 计算机视觉中的自注意力
点击上方"视学算法",选择加"星标"或"置顶" 重磅干货,第一时间送达 作者丨Branislav Holländer 来源丨CV技术指南 编 ...
- 理解计算机视觉中的损失函数
点击上方"小白学视觉",选择加"星标"或"置顶" 重磅干货,第一时间送达 导读 损失函数在模型的性能中起着关键作用.选择正确的损失函数可以帮 ...
最新文章
- c语言hook内核作用,hook内核函数的基本方法(附代码,经测试不蓝屏)
- 生成树的计数 Matrix-Tree(矩阵树)定理
- 别再说你不会ElasticSearch,都给你整理好了
- C# -- 文件的压缩与解压(GZipStream)
- 泛型集合 有序泛型 c#
- Android getRunningTasks和getRunningAppProcesses失效
- 【LeetCode】【字符串】题号:*151. 翻转字符串里的单词
- 关于验证码错误的解决办法
- 戴钊《自我教练:迈向自我实现之路》读书笔记
- 有哪些写项目策划书注意事项
- Unity-Photon Pun2个人总结
- 算法学习(2)----丢番图方程
- 学习笔记三:MLP基本原理、矩阵求导术推反向传播、激活函数、Xavier
- obs源码分析【五】:音频采集线程
- AES128加密算法与实现1
- Android Sensor感应器介绍(三)获取用户移动方向,指南针原理
- unity的异步加载
- INSERT 语句 增加条件,条件成立插入,否则不插入
- DM8060驱动器的使用
- JDBC技术预研选型依据【转】