Task

Write a program that performs so-called canny edge detection on an image.



A possible algorithm consists of the following steps:

 Noise reduction.   May be performed by Gaussian filter. 
  
 Compute intensity gradient   (matrices 




G

x




{\displaystyle G_{x}}

 and 




G

y




{\displaystyle G_{y}}

)   and its magnitude   



G


{\displaystyle G}

:
           



G
=



G

x


2


+

G

y


2






{\displaystyle G={\sqrt {G_{x}^{2}+G_{y}^{2}}}}


May be performed by convolution of an image with Sobel operators. 
  
 Non-maximum suppression.   
For each pixel compute the orientation of intensity gradient vector:   



θ
=


a
t
a
n
2



(

G

y


,


G

x


)



{\displaystyle \theta ={\rm {atan2}}\left(G_{y},\,G_{x}\right)}

.     
Transform   angle 



θ


{\displaystyle \theta }

   to one of four directions:   0, 45, 90, 135 degrees.     
Compute new array   



N


{\displaystyle N}

:     if         



G

(

p

a


)

<
G

(
p
)

<
G

(

p

b


)



{\displaystyle G\left(p_{a}\right)<G\left(p\right)<G\left(p_{b}\right)}


where   



p


{\displaystyle p}

   is the current pixel,   




p

a




{\displaystyle p_{a}}

   and   




p

b




{\displaystyle p_{b}}

   are the two neighbour pixels in the direction of gradient,   
then     



N
(
p
)
=
G
(
p
)


{\displaystyle N(p)=G(p)}

,       otherwise   



N
(
p
)
=
0


{\displaystyle N(p)=0}

.   
Nonzero pixels in resulting array correspond to local maxima of   



G


{\displaystyle G}

   in direction   



θ
(
p
)


{\displaystyle \theta (p)}

. 
  
 Tracing edges with hysteresis.   
At this stage two thresholds for the values of   



G


{\displaystyle G}

   are introduced:   




T

m
i
n




{\displaystyle T_{min}}

   and   




T

m
a
x




{\displaystyle T_{max}}

.   
Starting from pixels with   



N
(
p
)
⩾

T

m
a
x




{\displaystyle N(p)\geqslant T_{max}}

,   
find all paths of pixels with   



N
(
p
)
⩾

T

m
i
n




{\displaystyle N(p)\geqslant T_{min}}

   and put them to the resulting image.


import core.stdc.stdio, std.math, std.typecons, std.string, std.conv,
       std.algorithm, std.ascii, std.array, bitmap, grayscale_image;
 
enum maxBrightness = 255;
 
alias Pixel = short;
alias IntT = typeof(size_t.init.signed);
 
// If normalize is true, map pixels to range 0...maxBrightness.
void convolution(bool normalize)(in Pixel[] inp, Pixel[] outp,
                                 in float[] kernel,
                                 in IntT nx, in IntT ny, in IntT kn)
pure nothrow @nogc in {
    assert(kernel.length == kn ^^ 2);
    assert(kn % 2 == 1);
    assert(nx > kn && ny > kn);
    assert(inp.length == outp.length);
} body {
    //immutable IntT kn = sqrti(kernel.length);
    immutable IntT khalf = kn / 2;
 
    static if (normalize) {
        float pMin = float.max, pMax = -float.max;
 
        foreach (immutable m; khalf .. nx - khalf) {
            foreach (immutable n; khalf .. ny - khalf) {
                float pixel = 0.0;
                size_t c;
                foreach (immutable j; -khalf .. khalf + 1) {
                    foreach (immutable i; -khalf .. khalf + 1) {
                        pixel += inp[(n - j) * nx + m - i] * kernel[c];
                        c++;
                    }
                }
 
                if (pixel < pMin) pMin = pixel;
                if (pixel > pMax) pMax = pixel;
            }
        }
    }
 
    foreach (immutable m; khalf .. nx - khalf) {
        foreach (immutable n; khalf .. ny - khalf) {
            float pixel = 0.0;
            size_t c;
            foreach (immutable j; -khalf .. khalf + 1) {
                foreach (immutable i; -khalf .. khalf + 1) {
                    pixel += inp[(n - j) * nx + m - i] * kernel[c];
                    c++;
                }
            }
 
            static if (normalize)
                pixel = maxBrightness * (pixel - pMin) / (pMax - pMin);
            outp[n * nx + m] = cast(Pixel)pixel;
        }
    }
}
 
 
void gaussianFilter(in Pixel[] inp, Pixel[] outp,
                    in IntT nx, in IntT ny, in float sigma)
pure nothrow in {
    assert(inp.length == outp.length);
} body {
    immutable IntT n = 2 * cast(IntT)(2 * sigma) + 3;
    immutable float mean = floor(n / 2.0);
    auto kernel = new float[n * n];
 
    debug fprintf(stderr,
                  "gaussianFilter: kernel size %d, sigma=%g\n",
                  n, sigma);
 
    size_t c;
    foreach (immutable i; 0 .. n) {
        foreach (immutable j; 0 .. n) {
            kernel[c] = exp(-0.5 * (((i - mean) / sigma) ^^ 2 +
                                    ((j - mean) / sigma) ^^ 2))
                        / (2 * PI * sigma * sigma);
            c++;
        }
    }
 
    convolution!true(inp, outp, kernel, nx, ny, n);
}
 
 
Image!Pixel cannyEdgeDetection(in Image!Pixel inp,
                               in IntT tMin, in IntT tMax,
                               in float sigma)
pure nothrow in {
    assert(inp !is null);
} body {
    immutable IntT nx = inp.nx.signed;
    immutable IntT ny = inp.ny.signed;
    auto outp = new Pixel[nx * ny];
 
    gaussianFilter(inp.image, outp, nx, ny, sigma);
 
    static immutable float[] Gx = [-1, 0, 1,
                                   -2, 0, 2,
                                   -1, 0, 1];
    auto after_Gx = new Pixel[nx * ny];
    convolution!false(outp, after_Gx, Gx, nx, ny, 3);
 
    static immutable float[] Gy = [ 1, 2, 1,
                                    0, 0, 0,
                                   -1,-2,-1];
    auto after_Gy = new Pixel[nx * ny];
    convolution!false(outp, after_Gy, Gy, nx, ny, 3);
 
    auto G = new Pixel[nx * ny];
    foreach (i; 1 .. nx - 1)
        foreach (j; 1 .. ny - 1) {
            immutable size_t c = i + nx * j;
            G[c] = cast(Pixel)hypot(after_Gx[c], after_Gy[c]);
        }
 
    // Non-maximum suppression, straightforward implementation.
    auto nms = new Pixel[nx * ny];
    foreach (immutable i; 1 .. nx - 1)
        foreach (immutable j; 1 .. ny - 1) {
            immutable IntT c = i + nx * j,
                           nn = c - nx,
                           ss = c + nx,
                           ww = c + 1,
                           ee = c - 1,
                           nw = nn + 1,
                           ne = nn - 1,
                           sw = ss + 1,
                           se = ss - 1;
 
            immutable aux = atan2(double(after_Gy[c]),
                                  double(after_Gx[c])) + PI;
            immutable float dir = float((aux % PI) / PI) * 8;
 
            if (((dir <= 1 || dir > 7) && G[c] > G[ee] &&
                 G[c] > G[ww]) || // 0 deg.
                ((dir > 1 && dir <= 3) && G[c] > G[nw] &&
                 G[c] > G[se]) || // 45 deg.
                ((dir > 3 && dir <= 5) && G[c] > G[nn] &&
                 G[c] > G[ss]) || // 90 deg.
                ((dir > 5 && dir <= 7) && G[c] > G[ne] &&
                 G[c] > G[sw]))   // 135 deg.
                nms[c] = G[c];
            else
                nms[c] = 0;
        }
 
    // Reuse array used as a stack. nx*ny/2 elements should be enough.
    IntT[] edges = (cast(IntT*)after_Gy.ptr)[0 .. after_Gy.length / 2];
    outp[] = Pixel.init;
    edges[] = 0;
 
    // Tracing edges with hysteresis. Non-recursive implementation.
    size_t c = 1;
    foreach (immutable j; 1 .. ny - 1) {
        foreach (immutable i; 1 .. nx - 1) {
            if (nms[c] >= tMax && outp[c] == 0) { // Trace edges.
                outp[c] = maxBrightness;
                IntT nedges = 1;
                edges[0] = c;
 
                do {
                    nedges--;
                    immutable IntT t = edges[nedges];
 
                    immutable IntT[8] neighbours = [
                        t - nx,      // nn
                        t + nx,      // ss
                        t + 1,       // ww
                        t - 1,       // ee
                        t - nx + 1,  // nw
                        t - nx - 1,  // ne
                        t + nx + 1,  // sw
                        t + nx - 1]; // se
 
                    foreach (immutable n; neighbours)
                        if (nms[n] >= tMin && outp[n] == 0) {
                            outp[n] = maxBrightness;
                            edges[nedges] = n;
                            nedges++;
                        }
                } while (nedges > 0);
            }
            c++;
        }
    }
 
    return Image!Pixel.fromData(outp, nx, ny);
}
 
 
void main(in string[] args) {
    immutable fileName = (args.length == 2) ? args[1] : "lena.pgm";
    Image!Pixel imIn;
    imIn = imIn.loadPGM(fileName);
    printf("Image size: %d x %d\n", imIn.nx, imIn.ny);
    imIn.cannyEdgeDetection(45, 50, 1.0f).savePGM("lena_canny.pgm");
}