离散拉普拉斯算子与LOG推导

离散拉普拉斯算子的推导

连续函数的拉普拉斯算子如下：
∇2f=∂2f∂x2+∂2f∂y2\nabla^{2} f=\frac{\partial^{2} f}{\partial x^{2}}+\frac{\partial^{2} f}{\partial y^{2}} ∇2f=∂x2∂2f+∂y2∂2f
在该式中，f(x,y)f(x,y)f(x,y)代表图像(x,y)(x,y)(x,y)位置的像素值。

利用导数定义，我们可以得到f(x,y)f(x,y)f(x,y)关于x,yx,yx,y的两侧偏导数如下：
∂f(x,y)∂x=limh→0f(x+h,y)−f(x,y)h∂f(x,y)∂x=limh→0f(x,y)−f(x−h,y)h∂f(x,y)∂y=limh→0f(x,y+h)−f(x,y)h∂f(x,y)∂y=limh→0f(x,y)−f(x,y−h)h∂f(x,y)∂x=limh→0f(x+h,y)−f(x−h,y)2h∂f(x,y)∂y=limh→0f(x,y+h)−f(x,y−h)2h\begin{aligned} \frac{\partial f(x,y)}{\partial x}&=\mathop{lim}\limits_{h \rightarrow 0} \frac{f(x+h, y)-f(x, y)}{h} \\ \frac{\partial f(x,y)}{\partial x}&=\mathop{lim}\limits_{h \rightarrow 0} \frac{f(x, y)-f(x-h, y)}{h} \\ \frac{\partial f(x,y)}{\partial y}&=\mathop{lim}\limits_{h \rightarrow 0}\frac{f(x, y+h)-f(x, y)}{h} \\ \frac{\partial f(x,y)}{\partial y}&=\mathop{lim}\limits_{h \rightarrow 0}\frac{f(x, y)-f(x, y-h)}{h} \\ \frac{\partial f(x,y)}{\partial x}&=\mathop{lim}\limits_{h \rightarrow 0} \frac{f(x+h, y)-f(x-h, y)}{2h} \\ \frac{\partial f(x,y)}{\partial y}&=\mathop{lim}\limits_{h \rightarrow 0}\frac{f(x, y+h)-f(x, y-h)}{2h} \\ \end{aligned} ∂x∂f(x,y)∂x∂f(x,y)∂y∂f(x,y)∂y∂f(x,y)∂x∂f(x,y)∂y∂f(x,y)=h→0limhf(x+h,y)−f(x,y)=h→0limhf(x,y)−f(x−h,y)=h→0limhf(x,y+h)−f(x,y)=h→0limhf(x,y)−f(x,y−h)=h→0lim2hf(x+h,y)−f(x−h,y)=h→0lim2hf(x,y+h)−f(x,y−h)
对于函数f(x,y)f(x,y)f(x,y)的二阶导数，利用导数定义以及上述得到的偏导数，可以得到以下式子。
∂2f∂x2=limh→0∂f(x+h,y)∂x−∂f(x,y)∂xh=limh→0f(x+h,y)−f(x,y)h−f(x,y)−f(x−h,y)hh=limh→0f(x+h,y)+f(x−h,y)−2f(x,y)h2∂2f∂y2=limh→0∂f(x,y+h)∂x−∂f(x,y)∂xh=limh→0f(x,y+h)−f(x,y)h−f(x,y)−f(x,y−h)hh=limh→0f(x,y+h)+f(x,y−h)−2f(x,y)h2\begin{aligned} \frac{\partial^{2} f}{\partial x^{2}} &=\mathop{lim}\limits_{h \rightarrow 0}\frac{\frac{\partial f(x+h, y)}{\partial x}-\frac{\partial f(x, y)}{\partial x}}{h}\\ &=\mathop{lim}\limits_{h \rightarrow 0}\frac{\frac{f(x+h, y)-f(x, y)}{h}-\frac{f(x, y)-f(x-h, y)}{h}}{h}\\ &=\mathop{lim}\limits_{h \rightarrow 0}\frac{f(x+h, y)+f(x-h, y)-2f(x,y)}{h^2}\\ \frac{\partial^{2} f}{\partial y^{2}} &=\mathop{lim}\limits_{h \rightarrow 0}\frac{\frac{\partial f(x, y+h)}{\partial x}-\frac{\partial f(x, y)}{\partial x}}{h}\\ &=\mathop{lim}\limits_{h \rightarrow 0}\frac{\frac{f(x, y+h)-f(x, y)}{h}-\frac{f(x, y)-f(x, y-h)}{h}}{h}\\ &=\mathop{lim}\limits_{h \rightarrow 0}\frac{f(x, y+h)+f(x, y-h)-2f(x,y)}{h^2}\\ \end{aligned} ∂x2∂2f∂y2∂2f=h→0limh∂x∂f(x+h,y)−∂x∂f(x,y)=h→0limhhf(x+h,y)−f(x,y)−hf(x,y)−f(x−h,y)=h→0limh2f(x+h,y)+f(x−h,y)−2f(x,y)=h→0limh∂x∂f(x,y+h)−∂x∂f(x,y)=h→0limhhf(x,y+h)−f(x,y)−hf(x,y)−f(x,y−h)=h→0limh2f(x,y+h)+f(x,y−h)−2f(x,y)

由于图像之间的间距插值为1，所以我们使用1来近似代替hhh，由此可以得到以下等式：
∂2f∂x2≈f(x+1,y)+f(x−1,y)−2f(x,y)∂2f∂y2≈f(x,y+1)+f(x,y−1)−2f(x,y)\frac{\partial^{2} f}{\partial x^{2}} \approx f(x+1,y)+f(x-1,y)-2f(x,y)\\ \frac{\partial^{2} f}{\partial y^{2}} \approx f(x,y+1)+f(x,y-1)-2f(x,y)\\ ∂x2∂2f≈f(x+1,y)+f(x−1,y)−2f(x,y)∂y2∂2f≈f(x,y+1)+f(x,y−1)−2f(x,y)
综合上述两式，可以得到离散拉普拉斯算子如下：
∇d2f=f(x+1,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)−4f(x,y)\nabla^{2}_d f=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y) ∇d2f=f(x+1,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)−4f(x,y)
将其表示为卷积核的形式，可以得到如下卷积核：
∇d2=[0101−41010]\nabla_{d}^{2}=\left[\begin{array}{ccc} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{array}\right] ∇d2=⎣⎡0101−41010⎦⎤
但是这样的卷积核不是代表旋转情况的卷积，为了表示旋转情况的卷积我们提出以下方式：

如图所示，我们使用新的坐标(u,v)(u,v)(u,v)，在坐标系中[11−11][xy]=[uv]\left[\begin{array}{cc} 1& 1\\ -1 & 1 \end{array}\right]\left[\begin{array}{c}x\\y\end{array}\right]=\left[\begin{array}{c}u\\v\end{array}\right][1−111][xy]=[uv]，这个变换相当于对原坐标轴逆时针旋转45o45^o45o角，并拉伸为原来的2倍。

便于计算方便，列出后面需要用到的变换。

当ϕ(u,v)=f(x,y)\phi(u,v)=f(x,y)ϕ(u,v)=f(x,y)时，以下变换成立（以下证明都建立在这个基础上）：
ϕ(u+1,v)=f(x+1,y+1)ϕ(u,v+1)=f(x−1,y+1)ϕ(u−1,v)=f(x−1,y−1)ϕ(u,v−1)=f(x+1,y−1)\begin{aligned} &\phi(u+1,v)=f(x+1,y+1)\\ &\phi(u,v+1)=f(x-1,y+1)\\ &\phi(u-1,v)=f(x-1,y-1)\\ &\phi(u,v-1)=f(x+1,y-1)\\ \end{aligned} ϕ(u+1,v)=f(x+1,y+1)ϕ(u,v+1)=f(x−1,y+1)ϕ(u−1,v)=f(x−1,y−1)ϕ(u,v−1)=f(x+1,y−1)
由此我们可以得到新的拉普拉斯算子如下：
∇×2ϕ=∂2ϕ∂u2+∂2ϕ∂v2\nabla_{\times}^{2} \phi=\frac{\partial^{2} \phi}{\partial u^{2}}+\frac{\partial^{2} \phi}{\partial v^{2}} ∇×2ϕ=∂u2∂2ϕ+∂v2∂2ϕ
利用导数定义分别对∂2ϕ∂u2,∂2ϕ∂v2\frac{\partial^{2} \phi}{\partial u^{2}},\frac{\partial^{2} \phi}{\partial v^{2}}∂u2∂2ϕ,∂v2∂2ϕ进行化简可以得到：
∂2ϕ∂u2=limh→0∂ϕ(u+h,v)∂u−∂ϕ(u−h,v)∂u2h=limh→0ϕ(u+h,v)−ϕ(u,v)h−ϕ(u,v)−ϕ(u−h,v)h2h=limh→0ϕ(u+h,v)+ϕ(u−h,v)−2ϕ(u,v)2h2∂2ϕ∂v2=limh→0∂ϕ(u,v+h)∂v−∂ϕ(u,v−h)∂v2h=limh→0ϕ(u,v+h)−ϕ(u,v)h−ϕ(u,v)−ϕ(u,v−h)h2h=limh→0ϕ(u,v+h)+ϕ(u,v−h)−2ϕ(u,v)2h2\begin{aligned} \frac{\partial^{2} \phi}{\partial u^{2}} &=\mathop{lim}\limits_{h \rightarrow 0}\frac{\frac{\partial \phi(u+h, v)}{\partial u}-\frac{\partial \phi(u-h, v)}{\partial u}}{2h}\\ &=\mathop{lim}\limits_{h \rightarrow 0}\frac{\frac{\phi(u+h, v)-\phi(u, v)}{h}-\frac{\phi(u, v)-\phi(u-h, v)}{h}}{2h}\\ &=\mathop{lim}\limits_{h \rightarrow 0}\frac{\phi(u+h, v)+\phi(u-h, v)-2\phi(u,v)}{2h^2}\\ \frac{\partial^{2} \phi}{\partial v^{2}} &=\mathop{lim}\limits_{h \rightarrow 0}\frac{\frac{\partial \phi(u, v+h)}{\partial v}-\frac{\partial \phi(u, v-h)}{\partial v}}{2h}\\ &=\mathop{lim}\limits_{h \rightarrow 0}\frac{\frac{\phi(u, v+h)-\phi(u, v)}{h}-\frac{\phi(u, v)-\phi(u, v-h)}{h}}{2h}\\ &=\mathop{lim}\limits_{h \rightarrow 0}\frac{\phi(u, v+h)+\phi(u, v-h)-2\phi(u,v)}{2h^2}\\ \end{aligned} ∂u2∂2ϕ∂v2∂2ϕ=h→0lim2h∂u∂ϕ(u+h,v)−∂u∂ϕ(u−h,v)=h→0lim2hhϕ(u+h,v)−ϕ(u,v)−hϕ(u,v)−ϕ(u−h,v)=h→0lim2h2ϕ(u+h,v)+ϕ(u−h,v)−2ϕ(u,v)=h→0lim2h∂v∂ϕ(u,v+h)−∂v∂ϕ(u,v−h)=h→0lim2hhϕ(u,v+h)−ϕ(u,v)−hϕ(u,v)−ϕ(u,v−h)=h→0lim2h2ϕ(u,v+h)+ϕ(u,v−h)−2ϕ(u,v)
令h=2h=\sqrt{2}h=2可以得到以下形式。
∇×2≈ϕ(u+1,v)+ϕ(u−1,v)+ϕ(u,v+1)+ϕ(u,v−1)4−ϕ(u,v)=f(x+1,y+1)+f(x−1,y−1)+f(x−1,y+1)+f(x+1,y−1)2−2f(x,y)\begin{aligned} \nabla_{\times}^{2} &\approx \frac{\phi(u+1, v)+\phi(u-1, v)+\phi(u, v+1)+\phi(u, v-1)}{4}- \phi(u, v) \\ &=\frac{f(x+1, y+1)+f(x-1, y-1)+f(x-1, y+1)+f(x+1, y-1)}{2}-2 f(x, y) \end{aligned} ∇×2≈4ϕ(u+1,v)+ϕ(u−1,v)+ϕ(u,v+1)+ϕ(u,v−1)−ϕ(u,v)=2f(x+1,y+1)+f(x−1,y−1)+f(x−1,y+1)+f(x+1,y−1)−2f(x,y)
将其表示为卷积核的形式，可以得到：
∇×2=[0.500.50−200.500.5]\nabla_{\times}^{2}=\left[\begin{array}{ccc} 0.5 & 0 & 0.5 \\ 0 & -2 & 0 \\ 0.5 & 0 & 0.5 \end{array}\right] ∇×2=⎣⎡0.500.50−200.500.5⎦⎤
进行组合可以得到：
∇γ2=(1−γ)∇d2+γ∇×2=(1−γ)[0101−41010]+γ[0.500.50−200.500.5]\nabla_{\gamma}^{2}=(1-\gamma) \nabla_{d}^{2}+\gamma \nabla_{\times}^{2}=(1-\gamma)\left[\begin{array}{ccc} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{array}\right]+\gamma\left[\begin{array}{ccc} 0.5 & 0 & 0.5 \\ 0 & -2 & 0 \\ 0.5 & 0 & 0.5 \end{array}\right] ∇γ2=(1−γ)∇d2+γ∇×2=(1−γ)⎣⎡0101−41010⎦⎤+γ⎣⎡0.500.50−200.500.5⎦⎤
其中λ∈[0,1]\lambda \in [0,1]λ∈[0,1]

将两个矩阵相加可以得到：
∇γ2=[0.5λ1−λ0.5λ1−λ2λ−41−λ0.5λ1−λ0.5λ]\nabla^{2}_{\gamma}=\left[\begin{array}{ccc} 0.5 \lambda & 1-\lambda & 0.5 \lambda \\ 1-\lambda & 2 \lambda-4 & 1-\lambda \\ 0.5 \lambda & 1-\lambda & 0.5 \lambda \end{array}\right] ∇γ2=⎣⎡0.5λ1−λ0.5λ1−λ2λ−41−λ0.5λ1−λ0.5λ⎦⎤
令λ=2a1+a\lambda=\frac{2a}{1+a}λ=1+a2a，我们可以对上述式子进行化简，得到：
∇2=11+a[a1−aa1−a−41−aa1−aa]\nabla^{2}=\frac{1}{1+a}\left[\begin{array}{ccc} a & 1-a & a \\ 1-a & -4 & 1-a \\ a & 1-a & a \end{array}\right] ∇2=1+a1⎣⎡a1−aa1−a−41−aa1−aa⎦⎤

高斯拉普拉斯算子的推导

高斯函数如下：
Gσ(x,y)=12πσ2exp⁡(−x2+y22σ2)G_{\sigma}(x, y)=\frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \left(-\frac{x^{2}+y^{2}}{2 \sigma^{2}}\right) Gσ(x,y)=2πσ21exp(−2σ2x2+y2)

由于原图像可能存在噪声，所以需要先对图像进行高斯卷积再进行拉普拉斯变换：
∇2[Gσ(x,y)∗f(x,y)]=[∇2Gσ(x,y)]∗f(x,y)=LoG∗f(x,y)\nabla^{2}\left[G_{\sigma}(x, y) * f(x, y)\right]=\left[\nabla^{2} G_{\sigma}(x, y)\right] * f(x, y)=L o G * f(x, y) ∇2[Gσ(x,y)∗f(x,y)]=[∇2Gσ(x,y)]∗f(x,y)=LoG∗f(x,y)
由于卷积操作满足以下性质：
ddt[h(t)∗f(t)]=ddt∫f(τ)h(t−τ)dτ=∫f(τ)ddth(t−τ)dτ=f(t)∗ddth(t)\frac{d}{d t}[h(t) * f(t)]=\frac{d}{d t} \int f(\tau) h(t-\tau) d \tau=\int f(\tau) \frac{d}{d t} h(t-\tau) d \tau=f(t) * \frac{d}{d t} h(t) dtd[h(t)∗f(t)]=dtd∫f(τ)h(t−τ)dτ=∫f(τ)dtdh(t−τ)dτ=f(t)∗dtdh(t)
因此，LogLogLog可以先对高斯函数进行偏导操作，然后进行卷积求解，公式表示如下：
∂∂xGσ(x,y)=∂∂xe−(x2+y2)/2σ2=−xσ2e−(x2+y2)/2σ2\frac{\partial}{\partial x} G_{\sigma}(x, y)=\frac{\partial}{\partial x} e^{-\left(x^{2}+y^{2}\right) / 2 \sigma^{2}}=-\frac{x}{\sigma^{2}} e^{-\left(x^{2}+y^{2}\right) / 2 \sigma^{2}} ∂x∂Gσ(x,y)=∂x∂e−(x2+y2)/2σ2=−σ2xe−(x2+y2)/2σ2
求二阶导数可以得到：
∂2∂2xGσ(x,y)=x2σ4e−(x2+y2)/2σ2−1σ2e−(x2+y2)/2σ2=x2−σ2σ4e−(x2+y2)/2σ2\frac{\partial^{2}}{\partial^{2} x} G_{\sigma}(x, y)=\frac{x^{2}}{\sigma^{4}} e^{-\left(x^{2}+y^{2}\right) / 2 \sigma^{2}}-\frac{1}{\sigma^{2}} e^{-\left(x^{2}+y^{2}\right) / 2 \sigma^{2}}=\frac{x^{2}-\sigma^{2}}{\sigma^{4}} e^{-\left(x^{2}+y^{2}\right) / 2 \sigma^{2}} ∂2x∂2Gσ(x,y)=σ4x2e−(x2+y2)/2σ2−σ21e−(x2+y2)/2σ2=σ4x2−σ2e−(x2+y2)/2σ2
同理可得：
∂2∂2yGσ(x,y)=y2−σ2σ4e−(x2+y2)/2σ2\frac{\partial^{2}}{\partial^{2} y} G_{\sigma}(x, y)=\frac{y^{2}-\sigma^{2}}{\sigma^{4}} e^{-\left(x^{2}+y^{2}\right) / 2 \sigma^{2}} ∂2y∂2Gσ(x,y)=σ4y2−σ2e−(x2+y2)/2σ2
综上所述，可以得到LoGLoGLoG的定义如下：
LoG≜∇2Gσ(x,y)=∂2∂x2Gσ(x,y)+∂2∂y2Gσ(x,y)=x2+y2−2σ2σ4e−(x2+y2)/2σ2Lo G \triangleq \nabla^{2} G_{\sigma}(x, y)=\frac{\partial^{2}}{\partial x^{2}} G_{\sigma}(x, y)+\frac{\partial^{2}}{\partial y^{2}} G_{\sigma}(x, y)=\frac{x^{2}+y^{2}-2 \sigma^{2}}{\sigma^{4}} e^{-\left(x^{2}+y^{2}\right) / 2 \sigma^{2}} LoG≜∇2Gσ(x,y)=∂x2∂2Gσ(x,y)+∂y2∂2Gσ(x,y)=σ4x2+y2−2σ2e−(x2+y2)/2σ2
2维高斯拉普拉斯算子可以通过任何一个方形核进行逼近，只要保证该核的所有元素的和或均值为0，如下一个5×5的核进行逼近：
[001000121012−16210121000100]\left[\begin{array}{ccccc} 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 2 & 1 & 0 \\ 1 & 2 & -16 & 2 & 1 \\ 0 & 1 & 2 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 \end{array}\right] ⎣⎢⎢⎢⎢⎡001000121012−16210121000100⎦⎥⎥⎥⎥⎤
图像中的边缘检测可以使用以下步骤：

将 LoG 应用于图像
图像中过零点的检测
阈值的过零点只保留那些强的(正极大值与负极小值之间差别很大)