Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This numerical tour explores local differential operators (grad, div, laplacian) and their use to perform edge detection.
using PyPlot
using NtToolBox
using Autoreload
arequire("NtToolBox")
WARNING: Base.UTF8String is deprecated, use String instead.
likely near C:\Users\Ayman\.julia\v0.5\Autoreload\src\constants.jl:1
WARNING: Base.UTF8String is deprecated, use String instead.
likely near C:\Users\Ayman\.julia\v0.5\Autoreload\src\constants.jl:1
WARNING: Base.UTF8String is deprecated, use String instead.
likely near C:\Users\Ayman\.julia\v0.5\Autoreload\src\constants.jl:7
WARNING: Base.UTF8String is deprecated, use String instead.
likely near C:\Users\Ayman\.julia\v0.5\Autoreload\src\constants.jl:7
WARNING: Base.UTF8String is deprecated, use String instead.
likely near C:\Users\Ayman\.julia\v0.5\Autoreload\src\constants.jl:12
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[2]:4
WARNING: readall is deprecated, use readstring instead.
in depwarn(::String, ::Symbol) at .\deprecated.jl:64
in readall(::IOStream, ::Vararg{IOStream,N}) at .\deprecated.jl:30
in #parse_file#6(::Array{Any,1}, ::Function, ::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\files.jl:63
in parse_file(::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\files.jl:61
in #arequire#10(::Symbol, ::Array{String,1}, ::Function, ::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:65
in arequire(::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:47
in include_string(::String, ::String) at .\loading.jl:441
in execute_request(::ZMQ.Socket, ::IJulia.Msg) at C:\Users\Ayman\.julia\v0.5\IJulia\src\execute_request.jl:157
in eventloop(::ZMQ.Socket) at C:\Users\Ayman\.julia\v0.5\IJulia\src\eventloop.jl:8
in (::IJulia.##13#19)() at .\task.jl:360
while loading In[2], in expression starting on line 4
WARNING: replacing module NtToolBox
To obtain robust edge detection method, it is required to first remove the noise and small scale features in the image. This can be achieved using a linear blurring kernel.
Size of the image.
n = 256*2;
WARNING: Base.UTF8String is deprecated, use String
512
instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31 WARNING: Base.UTF8String is deprecated, use String instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31 WARNING: Base.UTF8String is deprecated, use String instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31 WARNING: Base.UTF8String is deprecated, use String instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31 WARNING: Base.UTF8String is deprecated, use String instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31 WARNING: Base.UTF8String is deprecated, use String instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31
Load an image $f_0$ of $N=n \times n$ pixels.
f0 = load_image("NtToolBox/src/data/hibiscus.png", n);
512×512 Array{Float32,2}: 0.116279 0.117647 0.110807 0.109439 … 0.0287278 0.0287278 0.0314638 0.123119 0.117647 0.113543 0.113543 0.0246238 0.0273598 0.0314638 0.119015 0.109439 0.105335 0.117647 0.0191518 0.0246238 0.0328317 0.113543 0.108071 0.108071 0.127223 0.0218878 0.0259918 0.0328317 0.121751 0.121751 0.131327 0.150479 0.0218878 0.0273598 0.0383037 0.145007 0.155951 0.161423 0.173735 … 0.0287278 0.0328317 0.0451437 0.168263 0.183311 0.187414 0.187414 0.0328317 0.0424077 0.0437756 0.19699 0.205198 0.179207 0.172367 0.0314638 0.0396717 0.0396717 0.202462 0.192886 0.179207 0.172367 0.0396717 0.0396717 0.0410397 0.201094 0.213406 0.202462 0.206566 0.0437756 0.0410397 0.0424077 0.236662 0.229822 0.23119 0.218878 … 0.0451437 0.0396717 0.0451437 0.247606 0.242134 0.243502 0.28591 0.0465116 0.0451437 0.0396717 0.242134 0.261286 0.268126 0.321477 0.0519836 0.0519836 0.0410397 ⋮ ⋱ ⋮ 0.214774 0.229822 0.239398 0.236662 … 0.0971272 0.105335 0.110807 0.195622 0.209302 0.222982 0.23803 0.0930232 0.0943913 0.101231 0.184679 0.195622 0.21067 0.227086 0.0916553 0.0916553 0.0930232 0.21751 0.195622 0.19015 0.213406 0.0848153 0.0875513 0.0930232 0.336525 0.192886 0.186046 0.205198 0.0793434 0.0807114 0.0889193 0.478796 0.29275 0.166895 0.19015 … 0.0834473 0.0793434 0.0848153 0.53762 0.474692 0.235294 0.175103 0.0820793 0.0807114 0.0807114 0.534884 0.525308 0.403557 0.216142 0.0711354 0.0766074 0.0875513 0.518468 0.491108 0.478796 0.300958 0.0725034 0.0725034 0.0807114 0.481532 0.491108 0.480164 0.439124 0.0697675 0.0711354 0.0725034 0.417237 0.47606 0.474692 0.48974 … 0.0642955 0.0683994 0.0766074 0.341997 0.429549 0.467852 0.473324 0.0711354 0.0766074 0.0766074
Display it.
figure(figsize=(5,5))
imageplot(f0)
Blurring is achieved using convolution: $$ f \star h(x) = \sum_y f(y-x) h(x) $$ where we assume periodic boundary condition.
This can be computed in $O(N\log(N))$ operations using the FFT, since $$ g = f \star h \qarrq \forall \om, \quad \hat g(\om) = \hat f(\om) \hat h(\om). $$
cconv = (f, h) -> real(plan_ifft((plan_fft(f)*f).*(plan_fft(h)*h))*((plan_fft(f)*f).*(plan_fft(h)*h)));
(::#1) (generic function with 1 method)
Define a Gaussian blurring kernel of width $\si$: $$ h_\si(x) = \frac{1}{Z} e^{ -\frac{x_1^2+x_2^2}{2\si^2} }$$ where $Z$ ensure that $\hat h(0)=1$.
include("NtToolBox/src/ndgrid.jl")
t = [collect(0 : Base.div(n, 2)); collect(-Base.div(n, 2) + 1 : -1)]
(X2, X1) = meshgrid(t, t)
normalize = h -> h./sum(h)
h = sigma -> normalize(exp(-(X1.^2 + X2.^2)./(2*sigma^2)));
WARNING: Method definition ndgrid(AbstractArray{T<:Any, 1
(::#21) (generic function with 1 method)
}) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:3 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:3. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:6 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:6. WARNING: Method definition ndgrid_fill(Any, Any, Any, Any) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:13 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:13. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}...) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:19 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:19. WARNING: Method definition meshgrid(AbstractArray{T<:Any, 1}) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:33 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:33. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:36 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:36. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:44 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:44.
Define blurring operator.
blur = (f, sigma) -> cconv(f, h(sigma));
(::#23) (generic function with 1 method)
Exercise 1
Test blurring with several blurring size $\si$.
include("NtSolutions/segmentation_1_edge_detection/exo1.jl");
## Insert your code here.
The simplest edge detectors only make use of the first order derivatives.
For continuous functions, the gradient reads $$ \nabla f(x) = \pa{ \pd{f(x)}{x_1}, \pd{f(x)}{x_2} } \in \RR^2. $$
We discretize this differential operator using first order finite differences. $$ (\nabla f)_i = ( f_{i_1,i_2}-f_{i_1-1,i_2}, f_{i_1,i_2}-f_{i_1,i_2-1} ) \in \RR^2. $$ Note that for simplity we use periodic boundary conditions.
Compute its gradient, using (here decentered) finite differences.
s = [[n]; collect(1:n-1)]
nabla = f -> cat(3, f - f[s, :], f - f[:, s]);
(::#25) (generic function with 1 method)
One thus has $ \nabla : \RR^N \mapsto \RR^{N \times 2}. $
v = nabla(f0);
512×512×2 Array{Float32,3}: [:, :, 1] = -0.225718 -0.311902 -0.357045 … -0.0478796 -0.0451436 0.00683992 0.0 0.00273598 -0.00136801 0.0 -0.00410395 -0.00820793 -0.00820793 -0.00273598 0.00136797 -0.00547195 -0.00136796 0.00273598 0.00136801 0.0 0.00820793 0.0136799 0.0232558 0.00136797 0.00547196 0.0232558 0.0341998 0.0300958 … 0.00547195 0.00683997 0.0232558 0.0273598 0.0259918 0.00957594 -0.00136801 0.0287277 0.0218878 -0.00820793 -0.00273598 -0.00410394 0.00547196 -0.0123119 0.0 0.0 0.00136797 -0.00136797 0.0205198 0.0232558 0.00136797 0.00136801 0.0355677 0.0164159 0.0287278 … -0.00136797 0.00273598 0.0109439 0.0123119 0.0123119 0.00547196 -0.00547196 -0.00547194 0.0191518 0.0246238 0.00683992 0.00136797 ⋮ ⋱ ⋮ -0.0123119 -0.00410399 0.00547194 … -0.0123119 -0.00820793 -0.0191519 -0.0205198 -0.0164159 -0.0109439 -0.00957594 -0.0109439 -0.0136799 -0.0123119 -0.00273598 -0.00820793 0.0328317 0.0 -0.0205199 -0.00410399 0.0 0.119015 -0.00273597 -0.00410399 -0.00683992 -0.00410395 0.142271 0.0998632 -0.0191519 … -0.00136801 -0.00410398 0.0588235 0.181943 0.0683994 0.00136801 -0.00410395 -0.00273597 0.0506155 0.168263 -0.00410399 0.00683992 -0.0164158 -0.0341997 0.0752394 -0.00410394 -0.00683992 -0.0369358 0.0 0.00136799 -0.00136801 -0.00820793 -0.0642955 -0.0150478 -0.00547197 … -0.00273598 0.00410394 -0.0752393 -0.0465117 -0.0068399 0.00820793 0.0 [:, :, 2] = 0.0848153 0.00136796 -0.00683992 … 0.0 0.00273598 0.0916552 -0.00547196 -0.00410394 0.00273598 0.00410399 0.0861833 -0.00957594 -0.00410394 0.00547196 0.00820793 0.0807114 -0.00547195 0.0 0.00410399 0.00683992 0.0834473 0.0 0.0095759 0.00547196 0.0109439 0.0998632 0.0109439 0.00547194 … 0.00410394 0.0123119 0.124487 0.0150479 0.00410394 0.00957594 0.00136797 0.157319 0.00820798 -0.0259918 0.00820794 0.0 0.161423 -0.0095759 -0.0136799 0.0 0.00136797 0.158687 0.0123119 -0.0109439 -0.00273598 0.00136801 0.191518 -0.00683996 0.001368 … -0.00547196 0.00547196 0.207934 -0.00547194 0.001368 -0.00136797 -0.00547196 0.201094 0.0191518 0.00683993 0.0 -0.0109439 ⋮ ⋱ ⋮ 0.103967 0.0150479 0.00957593 … 0.00820793 0.00547195 0.0943912 0.0136799 0.0136799 0.00136801 0.00683992 0.0916553 0.0109439 0.0150479 0.0 0.00136797 0.124487 -0.0218878 -0.00547194 0.00273597 0.00547196 0.247606 -0.143639 -0.00683996 0.00136801 0.00820793 0.393981 -0.186046 -0.125855 … -0.00410399 0.00547196 0.456908 -0.0629274 -0.239398 -0.00136797 0.0 0.447332 -0.00957596 -0.121751 0.00547195 0.0109439 0.437757 -0.0273599 -0.0123119 0.0 0.00820793 0.409029 0.00957593 -0.0109439 0.00136797 0.00136801 0.340629 0.0588236 -0.00136805 … 0.00410394 0.00820793 0.26539 0.0875513 0.0383037 0.00547195 0.0
One can display each of its components.
figure(figsize = (10, 10))
imageplot(v[:,:,1], L"\frac{d}{dx}", [1,2,1])
imageplot(v[:,:,2], L"\frac{d}{dy}", [1,2,2])
PyObject <matplotlib.text.Text object at 0x00000000235E6940>
A simple edge detector is simply obtained by obtained the gradient magnitude of a smoothed image.
A very simple edge detector is obtained by simply thresholding the gradient magnitude above some $t>0$. The set $\Ee$ of edges is then $$ \Ee = \enscond{x}{ d_\si(x) \geq t } $$ where we have defined $$ d_\si(x) = \norm{\nabla f_\si(x)}, \qwhereq f_\si = f_0 \star h_\si. $$
Compute $d_\si$ for $\si=1$.
sigma = 1
d = sqrt(sum(nabla(blur(f0, sigma)).^2, 3));
512×512×1 Array{Float64,3}: [:, :, 1] = 0.111691 0.133643 0.137846 … 0.0173866 0.0249714 0.0607401 0.0649781 0.0764041 0.0818717 0.0127883 0.0166122 0.0374897 0.0369233 0.0253904 0.0178443 0.00522808 0.00905602 0.0251312 0.0351577 0.0206691 0.0103345 0.00329904 0.0087954 0.0247427 0.0400816 0.0284165 0.0204969 0.00440289 0.0103128 0.0274222 0.0487043 0.037103 0.0249514 … 0.00606915 0.0123566 0.0318627 0.0564674 0.0395644 0.0189339 0.00663194 0.0131735 0.0358444 0.0612456 0.0369361 0.00856977 0.00585356 0.0126213 0.0386501 0.0639991 0.0366253 0.00772051 0.00512322 0.0111313 0.0399362 0.0685126 0.0415887 0.0175301 0.00456429 0.00965071 0.0416953 0.0750434 0.0482735 0.028286 … 0.00301873 0.00934612 0.0448561 0.079609 0.053542 0.0342373 0.00276138 0.00959219 0.0464144 0.0809133 0.0559552 0.0276 0.00357066 0.00956376 0.0457979 ⋮ ⋱ ⋮ 0.0462413 0.0320244 0.0113311 … 0.00866148 0.0127822 0.0300774 0.0452542 0.0332004 0.0159553 0.00673444 0.011488 0.0288872 0.0446182 0.028477 0.0149304 0.00436103 0.00923803 0.0289221 0.0555034 0.0198506 0.00557253 0.0038092 0.0102026 0.039367 0.0861361 0.033497 0.0279929 0.00311363 0.0148083 0.0621332 0.124356 0.0624582 0.0662958 … 0.00282461 0.0203358 0.0862395 0.15224 0.0795974 0.0954839 0.00457114 0.024641 0.101889 0.164258 0.0790181 0.0930554 0.00622951 0.0271562 0.107805 0.166341 0.0824006 0.0629 0.0059643 0.0270328 0.106231 0.161628 0.0939792 0.0294855 0.00485015 0.0249172 0.0993019 0.148531 0.10324 0.0214316 … 0.0044168 0.0218981 0.0872197 0.130637 0.120792 0.0879376 0.00965963 0.0214207 0.0727358
Display it.
figure(figsize=(5,5))
imageplot(d[:, :])
Exercise 2
For $\si=1$, study the influence of the threshold value $t$.
include("NtSolutions/segmentation_1_edge_detection/exo2.jl")
## Insert your code here.
Exercise 3
Study the influence of $\si$.
include("NtSolutions/segmentation_1_edge_detection/exo3.jl")
## Insert your code here.
Defining a Laplacian requires to define a divergence operator. The divergence operator maps vector field to images. For continuous vector fields $v(x) \in \RR^2$, it is defined as $$ \text{div}(v)(x) = \pd{v_1(x)}{x_1} + \pd{v_2(x)}{x_2} \in \RR. $$ It is minus the adjoint of the gadient, i.e. $\text{div} = - \nabla^*$.
It is discretized, for $v=(v^1,v^2)$ as $$ \text{div}(v)_i = v^1_{i_1+1,i_2} + v^2_{i_1,i_2+1}. $$
function diiv(v)
v0 = v[:, :, 1]
v1 = v[:, :, 2]
t = [collect(2:n); [1]]
return v0[t, :] - v0 + v1[:,t] - v1
end;
diiv (generic function with 1 method)
The Laplacian operatore is defined as $\Delta=\text{div} \circ \nabla = -\nabla^* \circ \nabla$. It is thus a negative symmetric operator.
delta = f -> diiv(nabla(f));
(::#27) (generic function with 1 method)
Display $\Delta f_0$.
figure(figsize=(5,5))
imageplot(delta(f0))
Check that the relation $ \norm{\nabla f} = - \dotp{\Delta f}{f}. $
dotp = (a, b) -> sum(a.*b)
@sprintf("Should be 0: %.10f ", (dotp(nabla(f0), nabla(f0)) + dotp(delta(f0), f0)))
"Should be 0: 0.0000000000 "
The zero crossing of the Laplacian is a well known edge detector. This requires first blurring the image (which is equivalent to blurring the laplacian). The set $\Ee$ of edges is defined as: $$ \Ee = \enscond{x}{ \Delta f_\si(x) = 0 } \qwhereq f_\si = f_0 \star h_\si . $$
It was proposed by Marr and Hildreth:
Marr, D. and Hildreth, E., Theory of edge detection, In Proc. of the Royal Society London B, 207:187-217, 1980.
Display the zero crossing.
sigma = 4
figure(figsize = (5, 5))
NtToolBox.plot_levelset(delta(blur(f0, sigma)), 0, f0)
Exercise 4
Study the influence of $\si$.
include("NtSolutions/segmentation_1_edge_detection/exo4.jl")
## Insert your code here.
Zero-crossing of the Laplacian can be shown to return false edges corresponding to local minima of the gradient magnitude. Moreover, this operator gives poor localization at curved edges.
In order to improve over this basic detector, more advanced edge detectors make use of the second order derivatives. Several authors have advocated for this choice, in particular:
Haralick, R., Digital step edges from zero crossing of second directional derivatives, IEEE Trans. on Pattern Analysis and Machine Intelligence, 6(1):58-68, 1984.
Canny, J., A computational approach to edge detection, IEEE Trans. on PAMI, 8(6):679-698, 1986
Deriche, R., Using Canny's criteria to derive a recursively implemented optimal edge detector. International Journal of Computer Vision, 1:167-187, 1987.
They define the edge locations $\Ee$ as the zero-crossings of the second-order directional derivative in the gradient direction. $$ \Ee = \enscond{x}{ \dotp{ H(x) \times g_\si(x) }{ g_\si(x) } = 0 } \qwhereq g_\si = \nabla ( f_0 \star h_\si ) $$ where $\times$ is the matrix-vector multiplication.
Define centered first order derivatives.
dx = f -> (f[s, :] - f[t, :])./2
dy = f -> transpose(dx(transpose(f)));
(::#33) (generic function with 1 method)
Define second order derivatives.
s = [collect(2:n); [1]]
t = [[n]; collect(1:n-1)]
d2x = f -> f[s, :] + f[t, :] - 2.*f
d2y = f -> transpose(d2x(transpose(f)))
dxy = f -> dy(dx(f));
(::#39) (generic function with 1 method)
Define Hessian operator.
hessian = f -> cat(3, d2x(f), dxy(f), d2y(f));
(::#41) (generic function with 1 method)
Compute $g_\si = \nabla (f_0 \star h_\si). $
sigma = 6
g = grad(blur(f0, sigma));
512×512×2 Array{Float64,3}: [:, :, 1] = -0.00426809 -0.00417673 -0.00405771 … -0.00436316 -0.00432992 -0.00377098 -0.00364972 -0.0034979 -0.00391764 -0.00385979 -0.00309232 -0.00294029 -0.00275778 -0.0033008 -0.00321191 -0.00227817 -0.00209688 -0.00188768 -0.00255221 -0.00242939 -0.00138676 -0.0011803 -0.000950231 -0.00172285 -0.00156717 -0.000480565 -0.000255373 -1.18645e-5 … -0.000868113 -0.000684543 0.000381659 0.000617179 0.000865469 -4.08408e-5 0.000162355 0.00115129 0.00138743 0.00163104 0.000714417 0.000926722 0.00179386 0.0020203 0.00224945 0.00136512 0.00157495 0.00229017 0.00249664 0.00270163 0.00189198 0.00208779 0.00263473 0.0028115 0.00298302 … 0.00228784 0.00245904 0.0028322 0.00297036 0.00309972 0.00255488 0.00269227 0.00289296 0.00298447 0.00306362 0.00270087 0.00279671 ⋮ ⋱ ⋮ -0.00219075 -0.00219995 -0.00224713 … -0.00230116 -0.00222475 -0.00217272 -0.00218353 -0.00224072 -0.00230013 -0.00221236 -0.00223733 -0.00224835 -0.0023111 -0.00237554 -0.00228092 -0.00240131 -0.00241172 -0.00247558 -0.00254151 -0.00244609 -0.00266761 -0.00267681 -0.00273727 -0.00279993 -0.00271034 -0.00302239 -0.00302962 -0.00308213 … -0.00313789 -0.00306022 -0.00343455 -0.00343845 -0.00347837 -0.00352713 -0.00346573 -0.00385795 -0.00385627 -0.003879 -0.003926 -0.00388259 -0.0042366 -0.00422605 -0.00422711 -0.00428399 -0.00425698 -0.00451199 -0.00448835 -0.00446362 -0.00454811 -0.0045326 -0.00463151 -0.00458996 -0.00453592 … -0.00467036 -0.00465862 -0.00455644 -0.00449209 -0.0044062 -0.00461481 -0.0045973 [:, :, 2] = 0.0093996 0.00872447 0.00787371 … 0.00970781 0.00993398 0.0098217 0.00949097 0.00884348 0.00801579 0.00971585 0.00996722 0.00988353 0.00961222 0.0089953 0.00819297 0.00974713 0.0100251 0.00997234 0.00976426 0.00917782 0.0084003 0.00981018 0.010114 0.0100919 0.00994554 0.00938702 0.00863141 0.00990952 0.0102368 0.0102431 0.010152 0.00961708 0.00887876 … 0.0100451 0.0103925 0.0104235 0.0103772 0.00986059 0.00913394 0.0102122 0.010576 0.0106275 0.0106127 0.0101089 0.00938796 0.0104025 0.0107792 0.0108468 0.0108489 0.0103525 0.00963156 0.0106049 0.0109915 0.0110714 0.0110753 0.0105816 0.0098554 0.010807 0.0112014 0.0112903 0.0112818 0.0107866 0.0100503 … 0.0109967 0.0113972 0.0114927 0.0114585 0.0109582 0.0102075 0.0111629 0.0115684 0.0116684 0.0115967 0.0110875 0.0103186 0.0112962 0.0117058 0.0118083 ⋮ ⋱ ⋮ 0.00931962 0.00896744 0.00838988 … 0.00878848 0.00923592 0.00941433 0.00931042 0.00892026 0.00831605 0.0089013 0.00931233 0.00944834 0.00929961 0.00886308 0.00822311 0.00902954 0.0094001 0.00948797 0.00928859 0.00880032 0.00811825 0.00916624 0.00949471 0.00953157 0.00927819 0.00873646 0.00800932 0.00930266 0.00959013 0.00957634 0.00926899 0.008676 0.0079045 … 0.00942939 0.00967972 0.00961908 0.00926176 0.0086235 0.00781199 0.0095378 0.0097574 0.0096569 0.00925785 0.00858358 0.00773959 0.00962148 0.0098188 0.00968809 0.00925953 0.00856085 0.00769425 0.00967745 0.00986221 0.00971273 0.00927007 0.00855979 0.00768173 0.0097069 0.00988923 0.00973311 0.0092937 0.00858453 0.00770619 … 0.00971529 0.00990474 0.00975372 0.00933526 0.00863857 0.00777003 0.00971174 0.00991647 0.00978083
Compute $h_\si = H (f_0 \star h_\si). $
H = hessian(blur(f0, sigma));
512×512×3 Array{Float64,3}: [:, :, 1] = 0.000288345 0.000315367 0.000348489 … 0.000251649 0.000267383 0.000497115 0.000527007 0.000559813 0.00044552 0.000470137 0.000678653 0.000709429 0.00074012 0.000616846 0.000647873 0.000814154 0.000843406 0.000870095 0.000748584 0.000782526 0.000891407 0.000916587 0.000937452 0.000829357 0.00086222 0.000906198 0.000924924 0.000938367 … 0.000854741 0.000882623 0.000862224 0.000872551 0.000877333 0.000827273 0.000846898 0.000769632 0.00077025 0.000765574 0.000755258 0.000764367 0.000642571 0.000632871 0.000618407 0.000650703 0.000648225 0.000496311 0.000476342 0.000452177 0.000526855 0.000512846 0.000344557 0.000314858 0.000281397 … 0.000395866 0.000371251 0.000197469 0.000158858 0.000116693 0.000267044 0.000233228 6.07637e-5 1.41095e-5 -3.6097e-5 0.000145987 0.000104439 ⋮ ⋱ ⋮ 7.45577e-5 7.07651e-5 5.63327e-5 … 5.5045e-5 6.878e-5 1.80261e-5 1.64187e-5 6.40779e-6 1.03242e-6 1.23958e-5 -6.46055e-5 -6.48138e-5 -7.03816e-5 -7.54119e-5 -6.85644e-5 -0.000163983 -0.000163369 -0.000164476 -0.00016597 -0.000165167 -0.000266293 -0.000265091 -0.000261695 -0.000258424 -0.000264254 -0.000354787 -0.000352817 -0.000344855 … -0.000337962 -0.000349872 -0.000412156 -0.00040883 -0.000396248 -0.000389235 -0.000405518 -0.000423397 -0.000417817 -0.000400623 -0.000398874 -0.000416852 -0.000378649 -0.000369782 -0.000348113 -0.000357986 -0.000374392 -0.000275394 -0.000262301 -0.000236504 -0.000264121 -0.000275624 -0.000119521 -0.000101603 -7.23028e-5 … -0.000122248 -0.000126019 7.50735e-5 9.78648e-5 0.000129718 5.55476e-5 6.13156e-5 [:, :, 2] = 6.46021e-5 9.01545e-5 0.000112668 … 1.37154e-5 3.83611e-5 9.08156e-5 0.000120865 0.000147524 3.2604e-5 6.04338e-5 0.000120422 0.000151906 0.000179712 6.02662e-5 8.87841e-5 0.000151031 0.000181258 0.000207539 9.35248e-5 0.000120628 0.000179843 0.000206754 0.00022943 0.000128354 0.000152533 0.00020401 0.000226308 0.000244024 … 0.000160491 0.00018091 0.000220999 0.000238128 0.000250251 0.000186051 0.000202512 0.000228883 0.00024089 0.000247382 0.00020204 0.000214843 0.000226515 0.000233835 0.00023505 0.000206658 0.000216404 0.00021355 0.00021676 0.000213226 0.000199368 0.000206735 0.000190326 0.000189937 0.000182157 … 0.000180725 0.000186272 0.000157635 0.000153953 0.000142287 0.000152025 0.000156051 0.000116461 0.000109544 9.41806e-5 0.000114867 0.000117351 ⋮ ⋱ ⋮ 1.19047e-5 -2.36323e-5 -5.0532e-5 … 8.62517e-5 5.03263e-5 1.34067e-5 -3.10931e-5 -6.7784e-5 0.000101308 5.94529e-5 1.53501e-5 -3.54418e-5 -7.94337e-5 0.00011183 6.6403e-5 1.67372e-5 -3.70089e-5 -8.51002e-5 0.000115789 6.96013e-5 1.69772e-5 -3.59823e-5 -8.4519e-5 0.000112039 6.81307e-5 1.60317e-5 -3.23496e-5 -7.75732e-5 … 0.000100602 6.1957e-5 1.44676e-5 -2.58895e-5 -6.43336e-5 8.27919e-5 5.20203e-5 1.33983e-5 -1.6219e-5 -4.50963e-5 6.11173e-5 4.01593e-5 1.43094e-5 -2.89145e-6 -2.04111e-5 3.89635e-5 2.88628e-5 1.87926e-5 1.44648e-5 8.9047e-6 2.00893e-5 2.08788e-5 2.82274e-5 3.59917e-5 4.17701e-5 … 8.01967e-6 1.87422e-5 4.34688e-5 6.14573e-5 7.68642e-5 5.44131e-6 2.43054e-5 [:, :, 3] = -0.000422098 -0.000675133 -0.000850755 … 0.000226166 -0.00011228 -0.000392562 -0.000647484 -0.000827693 0.000251366 -8.369e-5 -0.000360113 -0.00061692 -0.000802333 0.000277943 -5.27378e-5 -0.000327668 -0.000586442 -0.000777511 0.00030378 -2.20316e-5 -0.0002976 -0.000558527 -0.000755602 0.000327258 6.36107e-6 -0.000271537 -0.000534928 -0.000738321 … 0.000347388 3.10766e-5 -0.000250323 -0.000516612 -0.000726651 0.000363812 5.14851e-5 -0.000234108 -0.000503842 -0.000720915 0.000376714 6.75939e-5 -0.00022254 -0.000496364 -0.000720934 0.000386653 7.98589e-5 -0.000215018 -0.000493651 -0.000726245 0.000394345 8.89487e-5 -0.000210932 -0.000495132 -0.000736303 … 0.000400467 9.55127e-5 -0.000209849 -0.000500376 -0.00075064 0.000405484 9.99967e-5 -0.000211619 -0.000509174 -0.000768926 0.000409546 0.000102538 ⋮ ⋱ ⋮ -9.47091e-5 -0.000352187 -0.000577555 … 0.000447447 0.000178409 -0.000137916 -0.000390159 -0.000604215 0.000411024 0.000136009 -0.00018836 -0.000436535 -0.000639971 0.000370561 8.78754e-5 -0.000242972 -0.00048827 -0.000682071 0.000328472 3.68534e-5 -0.000298154 -0.000541726 -0.000727145 0.000287472 -1.37879e-5 -0.000350093 -0.000592987 -0.000771501 … 0.000250332 -6.06397e-5 -0.000395147 -0.000638258 -0.000811508 0.000219601 -0.000100496 -0.000430235 -0.000674274 -0.000843991 0.000197315 -0.000130708 -0.000453198 -0.000698675 -0.0008666 0.000184761 -0.000149489 -0.000463038 -0.000710274 -0.000878069 0.000182322 -0.000156119 -0.000460016 -0.000709169 -0.000878346 … 0.000189446 -0.000151016 -0.000445574 -0.000696682 -0.000868548 0.000204736 -0.000135643
Compute $ a_\si(x) = h_\si(x) \times g_\si (x) $ (this is a matrix times vector operation).
a = H[:, :, 1:2].*cat(3, g[:, :, 1], g[:, :, 1]) + H[:, :, 2:3].*cat(3, g[:, :, 2], g[:, :, 2]);
512×512×2 Array{Float64,3}: [:, :, 1] = -6.23451e-7 -5.30653e-7 -5.26956e-7 … -9.61736e-7 -7.80978e-7 -1.01268e-6 -8.54558e-7 -7.75649e-7 -1.42042e-6 -1.21733e-6 -9.4109e-7 -7.19485e-7 -5.68711e-7 -1.43191e-6 -1.19553e-6 -3.80072e-7 -1.04973e-7 1.00927e-7 -9.64639e-7 -6.83686e-7 5.52468e-7 8.58962e-7 1.08951e-6 -1.14928e-7 2.11177e-7 1.63563e-6 1.94022e-6 2.1555e-6 … 9.25883e-7 1.28153e-6 2.62242e-6 2.8866e-6 3.04508e-6 1.9339e-6 2.2897e-6 3.31514e-6 3.5038e-6 3.57109e-6 2.7174e-6 3.03873e-6 3.61011e-6 3.69937e-6 3.65498e-6 3.15977e-6 3.41681e-6 3.50177e-6 3.48294e-6 3.32304e-6 3.22999e-6 3.40482e-6 3.05503e-6 2.934e-6 2.67015e-6 … 2.96544e-6 3.05369e-6 2.36554e-6 2.15891e-6 1.81411e-6 2.44095e-6 2.44877e-6 1.52634e-6 1.25668e-6 8.61223e-7 1.7389e-6 1.67781e-6 ⋮ ⋱ ⋮ -5.23903e-8 -3.67601e-7 -5.50544e-7 … 6.69947e-7 3.2077e-7 8.5656e-8 -3.1321e-7 -5.78053e-7 9.41037e-7 5.34307e-7 2.87293e-7 -1.68399e-7 -4.90533e-7 1.23035e-6 7.8642e-7 5.4924e-7 6.83103e-8 -2.83693e-7 1.5212e-6 1.06742e-6 8.67883e-7 3.9524e-7 3.93905e-8 1.79804e-6 1.36866e-6 1.2209e-6 7.88237e-7 4.4971e-7 … 2.03429e-6 1.66665e-6 1.54957e-6 1.18248e-6 8.75727e-7 2.18072e-6 1.90777e-6 1.75748e-6 1.472e-6 1.20499e-6 2.16608e-6 2.00753e-6 1.73668e-6 1.53797e-6 1.31446e-6 1.91787e-6 1.87412e-6 1.41678e-6 1.30112e-6 1.12407e-6 1.39992e-6 1.45251e-6 8.15899e-7 7.75324e-7 6.49848e-7 … 6.50372e-7 7.69881e-7 6.37247e-8 9.12861e-8 2.56726e-8 -2.02383e-7 -4.41598e-8 [:, :, 2] = -4.24328e-6 -6.26672e-6 … 5.48777e-6 2.18689e-6 -1.26888e-6 -4.06826e-6 -6.16713e-6 5.61763e-6 2.37769e-6 -1.06041e-6 -3.83387e-6 -5.99603e-6 5.74557e-6 2.58748e-6 -8.11086e-7 -3.54351e-6 -5.76233e-6 5.90249e-6 2.83372e-6 -5.15394e-7 -3.20919e-6 -5.48694e-6 6.1139e-6 3.12894e-6 -1.73888e-7 -2.85468e-6 -5.20224e-6 … 6.38781e-6 3.47089e-6 2.00088e-7 -2.51331e-6 -4.94713e-6 6.71075e-6 3.84009e-6 5.80038e-7 -2.22101e-6 -4.75906e-6 7.05264e-6 4.20503e-6 9.3228e-7 -2.00797e-6 -4.66619e-6 7.37729e-6 4.53202e-6 1.22497e-6 -1.89233e-6 -4.68246e-6 7.6532e-6 4.7944e-6 1.43588e-6 -1.87823e-6 -4.8068e-6 … 7.86028e-6 4.97767e-6 1.55575e-6 -1.95811e-6 -5.02591e-6 7.99107e-6 5.0792e-6 1.58693e-6 -2.11717e-6 -5.31855e-6 8.04733e-6 5.1043e-6 1.539e-6 ⋮ ⋱ ⋮ -9.08733e-7 -3.10622e-6 … 5.80442e-6 3.9341e-6 1.56764e-6 -1.31318e-6 -3.41243e-6 5.59616e-6 3.59457e-6 1.15353e-6 -1.78602e-6 -3.78936e-6 5.37593e-6 3.21765e-6 6.82299e-7 -2.29706e-6 -4.20768e-6 5.15768e-6 2.82447e-6 1.8102e-7 -2.81161e-6 -4.63645e-6 4.96115e-6 2.44319e-6 -3.16695e-7 -3.29346e-6 -5.04675e-6 … 4.81087e-6 2.10747e-6 -7.729e-7 -3.70944e-6 -5.415e-6 4.73069e-6 1.85071e-6 -1.15077e-6 -4.03474e-6 -5.72514e-6 4.73558e-6 1.69744e-6 -1.42223e-6 -4.25703e-6 -5.96903e-6 4.82396e-6 1.65523e-6 -1.57481e-6 -4.37719e-6 -6.14472e-6 4.97552e-6 1.71166e-6 -1.61416e-6 -4.40599e-6 -6.25308e-6 … 5.157e-6 1.83895e-6 -1.56028e-6 -4.35761e-6 -6.29441e-6 5.3346e-6 2.00515e-6 -1.43844e-6
Display the level set of $\dotp{a_\si(x)}{g_\si(x)}$.
figure(figsize=(5,5))
NtToolBox.plot_levelset(sum(a.*g, 3)[:, :], 0, f0);
Exercise 5
Study the influence of $\si$.
include("NtSolutions/segmentation_1_edge_detection/exo5.jl")
## Insert your code here.