Edge Detection

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.

In [2]:
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

Diffusion and Convolution

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.

In [3]:
n = 256*2;
WARNING: Base.UTF8String is deprecated, use String
Out[3]:
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.

In [4]:
f0 = load_image("NtToolBox/src/data/hibiscus.png", n);
Out[4]:
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.

In [5]:
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). $$

In [6]:
cconv = (f, h) -> real(plan_ifft((plan_fft(f)*f).*(plan_fft(h)*h))*((plan_fft(f)*f).*(plan_fft(h)*h)));
Out[6]:
(::#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$.

In [10]:
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
Out[10]:
(::#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.

In [11]:
blur = (f, sigma) -> cconv(f, h(sigma));
Out[11]:
(::#23) (generic function with 1 method)

Exercise 1

Test blurring with several blurring size $\si$.

In [20]:
include("NtSolutions/segmentation_1_edge_detection/exo1.jl");
In [9]:
## Insert your code here.

Gradient Based Edge Detectiors

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.

In [22]:
s = [[n]; collect(1:n-1)]
nabla = f -> cat(3, f - f[s, :], f - f[:, s]);
Out[22]:
(::#25) (generic function with 1 method)

One thus has $ \nabla : \RR^N \mapsto \RR^{N \times 2}. $

In [23]:
v = nabla(f0);
Out[23]:
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.

In [24]:
figure(figsize = (10, 10))
imageplot(v[:,:,1], L"\frac{d}{dx}", [1,2,1])
imageplot(v[:,:,2], L"\frac{d}{dy}", [1,2,2])
Out[24]:
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$.

In [25]:
sigma = 1
d = sqrt(sum(nabla(blur(f0, sigma)).^2, 3));
Out[25]:
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.

In [26]:
figure(figsize=(5,5))
imageplot(d[:, :])