edgedetection.hxx

Go to the documentation of this file.

/************************************************************************/
/*                                                                      */
/*               Copyright 1998-2002 by Ullrich Koethe                  */
/*       Cognitive Systems Group, University of Hamburg, Germany        */
/*                                                                      */
/*    This file is part of the VIGRA computer vision library.           */
/*    ( Version 1.5.0, Dec 07 2006 )                                    */
/*    The VIGRA Website is                                              */
/*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
/*    Please direct questions, bug reports, and contributions to        */
/*        koethe@informatik.uni-hamburg.de          or                  */
/*        vigra@kogs1.informatik.uni-hamburg.de                         */
/*                                                                      */
/*    Permission is hereby granted, free of charge, to any person       */
/*    obtaining a copy of this software and associated documentation    */
/*    files (the "Software"), to deal in the Software without           */
/*    restriction, including without limitation the rights to use,      */
/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
/*    sell copies of the Software, and to permit persons to whom the    */
/*    Software is furnished to do so, subject to the following          */
/*    conditions:                                                       */
/*                                                                      */
/*    The above copyright notice and this permission notice shall be    */
/*    included in all copies or substantial portions of the             */
/*    Software.                                                         */
/*                                                                      */
/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
/*                                                                      */
/************************************************************************/


#ifndef VIGRA_EDGEDETECTION_HXX
#define VIGRA_EDGEDETECTION_HXX

#include <vector>
#include <queue>
#include <cmath>     // sqrt(), abs()
#include "vigra/utilities.hxx"
#include "vigra/numerictraits.hxx"
#include "vigra/stdimage.hxx"
#include "vigra/stdimagefunctions.hxx"
#include "vigra/recursiveconvolution.hxx"
#include "vigra/separableconvolution.hxx"
#include "vigra/labelimage.hxx"
#include "vigra/mathutil.hxx"
#include "vigra/pixelneighborhood.hxx"
#include "vigra/linear_solve.hxx"


namespace vigra {


/********************************************************/
/*                                                      */
/*           differenceOfExponentialEdgeImage           */
/*                                                      */
/********************************************************/

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
void differenceOfExponentialEdgeImage(
               SrcIterator sul, SrcIterator slr, SrcAccessor sa,
           DestIterator dul, DestAccessor da,
           double scale, GradValue gradient_threshold, DestValue edge_marker)
{
    vigra_precondition(scale > 0,
                 "differenceOfExponentialEdgeImage(): scale > 0 required.");

    vigra_precondition(gradient_threshold > 0,
                 "differenceOfExponentialEdgeImage(): "
         "gradient_threshold > 0 required.");

    int w = slr.x - sul.x;
    int h = slr.y - sul.y;
    int x,y;

    typedef typename
        NumericTraits<typename SrcAccessor::value_type>::RealPromote
    TMPTYPE;
    typedef BasicImage<TMPTYPE> TMPIMG;

    TMPIMG tmp(w,h);
    TMPIMG smooth(w,h);

    recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0);
    recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0);

    recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale);
    recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale);

    typename TMPIMG::Iterator iy = smooth.upperLeft();
    typename TMPIMG::Iterator ty = tmp.upperLeft();
    DestIterator              dy = dul;

    static const Diff2D right(1, 0);
    static const Diff2D bottom(0, 1);


    TMPTYPE thresh = (gradient_threshold * gradient_threshold) *
                     NumericTraits<TMPTYPE>::one();
    TMPTYPE zero = NumericTraits<TMPTYPE>::zero();

    for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, ++dy.y)
    {
        typename TMPIMG::Iterator ix = iy;
        typename TMPIMG::Iterator tx = ty;
        DestIterator              dx = dy;

        for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x)
        {
            TMPTYPE diff = *tx - *ix;
            TMPTYPE gx = tx[right] - *tx;
            TMPTYPE gy = tx[bottom] - *tx;

            if((gx * gx > thresh) &&
                (diff * (tx[right] - ix[right]) < zero))
            {
                if(gx < zero)
                {
                    da.set(edge_marker, dx, right);
                }
                else
                {
                    da.set(edge_marker, dx);
                }
            }
            if(((gy * gy > thresh) &&
                (diff * (tx[bottom] - ix[bottom]) < zero)))
            {
                if(gy < zero)
                {
                    da.set(edge_marker, dx, bottom);
                }
                else
                {
                    da.set(edge_marker, dx);
                }
            }
        }
    }

    typename TMPIMG::Iterator ix = iy;
    typename TMPIMG::Iterator tx = ty;
    DestIterator              dx = dy;

    for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x)
    {
        TMPTYPE diff = *tx - *ix;
        TMPTYPE gx = tx[right] - *tx;

        if((gx * gx > thresh) &&
           (diff * (tx[right] - ix[right]) < zero))
        {
            if(gx < zero)
            {
                da.set(edge_marker, dx, right);
            }
            else
            {
                da.set(edge_marker, dx);
            }
        }
    }
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue>
inline
void differenceOfExponentialEdgeImage(
               SrcIterator sul, SrcIterator slr, SrcAccessor sa,
           DestIterator dul, DestAccessor da,
           double scale, GradValue gradient_threshold)
{
    differenceOfExponentialEdgeImage(sul, slr, sa, dul, da,
                                        scale, gradient_threshold, 1);
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
      class GradValue, class DestValue>
inline
void differenceOfExponentialEdgeImage(
           triple<SrcIterator, SrcIterator, SrcAccessor> src,
       pair<DestIterator, DestAccessor> dest,
       double scale, GradValue gradient_threshold,
       DestValue edge_marker)
{
    differenceOfExponentialEdgeImage(src.first, src.second, src.third,
                                        dest.first, dest.second,
                    scale, gradient_threshold,
                    edge_marker);
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue>
inline
void differenceOfExponentialEdgeImage(
           triple<SrcIterator, SrcIterator, SrcAccessor> src,
       pair<DestIterator, DestAccessor> dest,
       double scale, GradValue gradient_threshold)
{
    differenceOfExponentialEdgeImage(src.first, src.second, src.third,
                                        dest.first, dest.second,
                    scale, gradient_threshold, 1);
}

/********************************************************/
/*                                                      */
/*        differenceOfExponentialCrackEdgeImage         */
/*                                                      */
/********************************************************/

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
void differenceOfExponentialCrackEdgeImage(
               SrcIterator sul, SrcIterator slr, SrcAccessor sa,
           DestIterator dul, DestAccessor da,
           double scale, GradValue gradient_threshold,
           DestValue edge_marker)
{
    vigra_precondition(scale > 0,
                 "differenceOfExponentialCrackEdgeImage(): scale > 0 required.");

    vigra_precondition(gradient_threshold > 0,
                 "differenceOfExponentialCrackEdgeImage(): "
         "gradient_threshold > 0 required.");

    int w = slr.x - sul.x;
    int h = slr.y - sul.y;
    int x, y;

    typedef typename
        NumericTraits<typename SrcAccessor::value_type>::RealPromote
    TMPTYPE;
    typedef BasicImage<TMPTYPE> TMPIMG;

    TMPIMG tmp(w,h);
    TMPIMG smooth(w,h);

    TMPTYPE zero = NumericTraits<TMPTYPE>::zero();

    static const Diff2D right(1,0);
    static const Diff2D bottom(0,1);
    static const Diff2D left(-1,0);
    static const Diff2D top(0,-1);

    recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0);
    recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0);

    recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale);
    recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale);

    typename TMPIMG::Iterator iy = smooth.upperLeft();
    typename TMPIMG::Iterator ty = tmp.upperLeft();
    DestIterator              dy = dul;

    TMPTYPE thresh = (gradient_threshold * gradient_threshold) *
                     NumericTraits<TMPTYPE>::one();

    // find zero crossings above threshold
    for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, dy.y+=2)
    {
        typename TMPIMG::Iterator ix = iy;
        typename TMPIMG::Iterator tx = ty;
        DestIterator              dx = dy;

        for(int x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2)
        {
            TMPTYPE diff = *tx - *ix;
            TMPTYPE gx = tx[right] - *tx;
            TMPTYPE gy = tx[bottom] - *tx;

            if((gx * gx > thresh) &&
               (diff * (tx[right] - ix[right]) < zero))
            {
                da.set(edge_marker, dx, right);
            }
            if((gy * gy > thresh) &&
               (diff * (tx[bottom] - ix[bottom]) < zero))
            {
                da.set(edge_marker, dx, bottom);
            }
        }

        TMPTYPE diff = *tx - *ix;
        TMPTYPE gy = tx[bottom] - *tx;

        if((gy * gy > thresh) &&
           (diff * (tx[bottom] - ix[bottom]) < zero))
        {
            da.set(edge_marker, dx, bottom);
        }
    }

    typename TMPIMG::Iterator ix = iy;
    typename TMPIMG::Iterator tx = ty;
    DestIterator              dx = dy;

    for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2)
    {
        TMPTYPE diff = *tx - *ix;
        TMPTYPE gx = tx[right] - *tx;

        if((gx * gx > thresh) &&
           (diff * (tx[right] - ix[right]) < zero))
        {
            da.set(edge_marker, dx, right);
        }
    }

    iy = smooth.upperLeft() + Diff2D(0,1);
    ty = tmp.upperLeft() + Diff2D(0,1);
    dy = dul + Diff2D(1,2);

    static const Diff2D topleft(-1,-1);
    static const Diff2D topright(1,-1);
    static const Diff2D bottomleft(-1,1);
    static const Diff2D bottomright(1,1);

    // find missing 1-cells below threshold (x-direction)
    for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2)
    {
        typename TMPIMG::Iterator ix = iy;
        typename TMPIMG::Iterator tx = ty;
        DestIterator              dx = dy;

        for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2)
        {
            if(da(dx) == edge_marker) continue;

            TMPTYPE diff = *tx - *ix;

            if((diff * (tx[right] - ix[right]) < zero) &&
               (((da(dx, bottomright) == edge_marker) &&
                 (da(dx, topleft) == edge_marker)) ||
                ((da(dx, bottomleft) == edge_marker) &&
                 (da(dx, topright) == edge_marker))))

            {
                da.set(edge_marker, dx);
            }
        }
    }

    iy = smooth.upperLeft() + Diff2D(1,0);
    ty = tmp.upperLeft() + Diff2D(1,0);
    dy = dul + Diff2D(2,1);

    // find missing 1-cells below threshold (y-direction)
    for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2)
    {
        typename TMPIMG::Iterator ix = iy;
        typename TMPIMG::Iterator tx = ty;
        DestIterator              dx = dy;

        for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2)
        {
            if(da(dx) == edge_marker) continue;

            TMPTYPE diff = *tx - *ix;

            if((diff * (tx[bottom] - ix[bottom]) < zero) &&
               (((da(dx, bottomright) == edge_marker) &&
                 (da(dx, topleft) == edge_marker)) ||
                ((da(dx, bottomleft) == edge_marker) &&
                 (da(dx, topright) == edge_marker))))

            {
                da.set(edge_marker, dx);
            }
        }
    }

    dy = dul + Diff2D(1,1);

    // find missing 0-cells
    for(y=0; y<h-1; ++y, dy.y+=2)
    {
        DestIterator              dx = dy;

        for(int x=0; x<w-1; ++x, dx.x+=2)
        {
            static const Diff2D dist[] = {right, top, left, bottom };

            int i;
            for(i=0; i<4; ++i)
            {
            if(da(dx, dist[i]) == edge_marker) break;
            }

            if(i < 4) da.set(edge_marker, dx);
        }
    }
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
inline
void differenceOfExponentialCrackEdgeImage(
           triple<SrcIterator, SrcIterator, SrcAccessor> src,
       pair<DestIterator, DestAccessor> dest,
       double scale, GradValue gradient_threshold,
       DestValue edge_marker)
{
    differenceOfExponentialCrackEdgeImage(src.first, src.second, src.third,
                                        dest.first, dest.second,
                    scale, gradient_threshold,
                    edge_marker);
}

/********************************************************/
/*                                                      */
/*                  removeShortEdges                    */
/*                                                      */
/********************************************************/

template <class Iterator, class Accessor, class Value>
void removeShortEdges(
               Iterator sul, Iterator slr, Accessor sa,
           unsigned int min_edge_length, Value non_edge_marker)
{
    int w = slr.x - sul.x;
    int h = slr.y - sul.y;
    int x,y;

    IImage labels(w, h);
    labels = 0;

    int number_of_regions =
                labelImageWithBackground(srcIterRange(sul,slr,sa),
                                     destImage(labels), true, non_edge_marker);

    ArrayOfRegionStatistics<FindROISize<int> >
                                         region_stats(number_of_regions);

    inspectTwoImages(srcImageRange(labels), srcImage(labels), region_stats);

    IImage::Iterator ly = labels.upperLeft();
    Iterator oy = sul;

    for(y=0; y<h; ++y, ++oy.y, ++ly.y)
    {
        Iterator ox(oy);
        IImage::Iterator lx(ly);

        for(x=0; x<w; ++x, ++ox.x, ++lx.x)
        {
            if(sa(ox) == non_edge_marker) continue;
            if((region_stats[*lx].count) < min_edge_length)
            {
                 sa.set(non_edge_marker, ox);
            }
        }
    }
}

template <class Iterator, class Accessor, class Value>
inline
void removeShortEdges(
           triple<Iterator, Iterator, Accessor> src,
       unsigned int min_edge_length, Value non_edge_marker)
{
    removeShortEdges(src.first, src.second, src.third,
                     min_edge_length, non_edge_marker);
}

/********************************************************/
/*                                                      */
/*             closeGapsInCrackEdgeImage                */
/*                                                      */
/********************************************************/

template <class SrcIterator, class SrcAccessor, class SrcValue>
void closeGapsInCrackEdgeImage(
               SrcIterator sul, SrcIterator slr, SrcAccessor sa,
           SrcValue edge_marker)
{
    int w = (slr.x - sul.x) / 2;
    int h = (slr.y - sul.y) / 2;
    int x, y;

    int count1, count2, count3;

    static const Diff2D right(1,0);
    static const Diff2D bottom(0,1);
    static const Diff2D left(-1,0);
    static const Diff2D top(0,-1);

    static const Diff2D leftdist[] = {
        Diff2D(0, 0), Diff2D(-1, 1), Diff2D(-2, 0), Diff2D(-1, -1)};
    static const Diff2D rightdist[] = {
        Diff2D(2, 0), Diff2D(1, 1), Diff2D(0, 0), Diff2D(1, -1)};
    static const Diff2D topdist[] = {
        Diff2D(1, -1), Diff2D(0, 0), Diff2D(-1, -1), Diff2D(0, -2)};
    static const Diff2D bottomdist[] = {
        Diff2D(1, 1), Diff2D(0, 2), Diff2D(-1, 1), Diff2D(0, 0)};

    int i;

    SrcIterator sy = sul + Diff2D(0,1);
    SrcIterator sx;

    // close 1-pixel wide gaps (x-direction)
    for(y=0; y<h; ++y, sy.y+=2)
    {
        sx = sy + Diff2D(2,0);

        for(x=2; x<w; ++x, sx.x+=2)
        {
            if(sa(sx) == edge_marker) continue;

            if(sa(sx, left) != edge_marker) continue;
            if(sa(sx, right) != edge_marker) continue;

            count1 = 0;
            count2 = 0;
            count3 = 0;

            for(i=0; i<4; ++i)
            {
                if(sa(sx, leftdist[i]) == edge_marker)
                {
                    ++count1;
                    count3 ^= 1 << i;
                }
                if(sa(sx, rightdist[i]) == edge_marker)
                {
                    ++count2;
                    count3 ^= 1 << i;
                }
            }

            if(count1 <= 1 || count2 <= 1 || count3 == 15)
            {
                sa.set(edge_marker, sx);
            }
        }
   }

    sy = sul + Diff2D(1,2);

    // close 1-pixel wide gaps (y-direction)
    for(y=2; y<h; ++y, sy.y+=2)
    {
        sx = sy;

        for(x=0; x<w; ++x, sx.x+=2)
        {
            if(sa(sx) == edge_marker) continue;

            if(sa(sx, top) != edge_marker) continue;
            if(sa(sx, bottom) != edge_marker) continue;

            count1 = 0;
            count2 = 0;
            count3 = 0;

            for(i=0; i<4; ++i)
            {
                if(sa(sx, topdist[i]) == edge_marker)
                {
                    ++count1;
                    count3 ^= 1 << i;
                }
                if(sa(sx, bottomdist[i]) == edge_marker)
                {
                    ++count2;
                    count3 ^= 1 << i;
                }
            }

            if(count1 <= 1 || count2 <= 1 || count3 == 15)
            {
                sa.set(edge_marker, sx);
            }
        }
    }
}

template <class SrcIterator, class SrcAccessor, class SrcValue>
inline
void closeGapsInCrackEdgeImage(
           triple<SrcIterator, SrcIterator, SrcAccessor> src,
       SrcValue edge_marker)
{
    closeGapsInCrackEdgeImage(src.first, src.second, src.third,
                    edge_marker);
}

/********************************************************/
/*                                                      */
/*              beautifyCrackEdgeImage                  */
/*                                                      */
/********************************************************/

template <class SrcIterator, class SrcAccessor, class SrcValue>
void beautifyCrackEdgeImage(
               SrcIterator sul, SrcIterator slr, SrcAccessor sa,
           SrcValue edge_marker, SrcValue background_marker)
{
    int w = (slr.x - sul.x) / 2;
    int h = (slr.y - sul.y) / 2;
    int x, y;

    SrcIterator sy = sul + Diff2D(1,1);
    SrcIterator sx;

    static const Diff2D right(1,0);
    static const Diff2D bottom(0,1);
    static const Diff2D left(-1,0);
    static const Diff2D top(0,-1);

    //  delete 0-cells at corners
    for(y=0; y<h; ++y, sy.y+=2)
    {
        sx = sy;

        for(x=0; x<w; ++x, sx.x+=2)
        {
            if(sa(sx) != edge_marker) continue;

            if(sa(sx, right) == edge_marker && sa(sx, left) == edge_marker) continue;
            if(sa(sx, bottom) == edge_marker && sa(sx, top) == edge_marker) continue;

            sa.set(background_marker, sx);
        }
    }
}

template <class SrcIterator, class SrcAccessor, class SrcValue>
inline
void beautifyCrackEdgeImage(
           triple<SrcIterator, SrcIterator, SrcAccessor> src,
           SrcValue edge_marker, SrcValue background_marker)
{
    beautifyCrackEdgeImage(src.first, src.second, src.third,
                    edge_marker, background_marker);
}


class Edgel
{
  public:
    float x;

    float y;

    float strength;

    float orientation;

    Edgel()
    : x(0.0f), y(0.0f), strength(0.0f), orientation(0.0f)
    {}

    Edgel(float ix, float iy, float is, float io)
    : x(ix), y(iy), strength(is), orientation(io)
    {}
};

template <class Image1, class Image2, class BackInsertable>
void internalCannyFindEdgels(Image1 const & gx,
                             Image1 const & gy,
                             Image2 const & magnitude,
                             BackInsertable & edgels, std::vector<int> p)
{

    typedef typename Image1::value_type PixelType;
    double t = 0.5 / VIGRA_CSTD::sin(M_PI/8.0);

    //last element in edgel list is edgel that holds orientation
    //of interest point

    //orientation assignment
    std::vector<int > point= p;

    PixelType gradx = gx(p[0],p[1]);
    PixelType grady = gy(p[0],p[1]);

    double orientation = VIGRA_CSTD::atan2(-grady, gradx) - M_PI * 1.5;
    if(orientation < 0.0)
        orientation += 2.0*M_PI;
    Edgel edgel1;
    edgel1.orientation=orientation;
    edgels.push_back(edgel1);
    //EOF orientation assignment

    for(int y=1; y<gx.height()-1; ++y)
    {
        for(int x=1; x<gx.width()-1; ++x)
        {
            gradx = gx(x,y);
            grady = gy(x,y);
            double mag = magnitude(x, y);

            int dx = (int)VIGRA_CSTD::floor(gradx*t/mag + 0.5);
            int dy = (int)VIGRA_CSTD::floor(grady*t/mag + 0.5);

            int x1 = x - dx,
                x2 = x + dx,
                y1 = y - dy,
                y2 = y + dy;

            PixelType m1 = magnitude(x1, y1);
            PixelType m3 = magnitude(x2, y2);

            if(m1 < mag && m3 <= mag)
            {
                Edgel edgel;

                // local maximum => quadratic interpolation of sub-pixel location
                PixelType del = (m1 - m3) / 2.0 / (m1 + m3 - 2.0*mag);
                edgel.x = x + dx*del;
                edgel.y = y + dy*del;
                edgel.strength = mag;
                orientation = VIGRA_CSTD::atan2(-grady, gradx) - M_PI * 1.5;
                if(orientation < 0.0)
                    orientation += 2.0*M_PI;
                edgel.orientation = orientation;
                edgels.push_back(edgel);
            }
        }
    }
}

/********************************************************/
/*                                                      */
/*                      cannyEdgelList                  */
/*                                                      */
/********************************************************/

template <class SrcIterator, class SrcAccessor, class BackInsertable>
void cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
                        BackInsertable & edgels, double scale, std::vector<int> p)
{
    int w = lr.x - ul.x;
    int h = lr.y - ul.y;

    // calculate image gradients
    typedef typename
        NumericTraits<typename SrcAccessor::value_type>::RealPromote
        TmpType;

    BasicImage<TmpType> tmp(w,h), dx(w,h), dy(w,h);

    Kernel1D<double> smooth, grad;

    smooth.initGaussian(scale);
    grad.initGaussianDerivative(scale, 1);

    separableConvolveX(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad));
    separableConvolveY(srcImageRange(tmp), destImage(dx), kernel1d(smooth));

    separableConvolveY(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad));
    separableConvolveX(srcImageRange(tmp), destImage(dy), kernel1d(smooth));

    combineTwoImages(srcImageRange(dx), srcImage(dy), destImage(tmp),
                     MagnitudeFunctor<TmpType>());


    // find edgels
    internalCannyFindEdgels(dx, dy, tmp, edgels, p);
}

template <class SrcIterator, class SrcAccessor, class BackInsertable>
inline void
cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
               BackInsertable & edgels, double scale, std::vector<int> p)
{
    cannyEdgelList(src.first, src.second, src.third, edgels, scale,p);
}

/********************************************************/
/*                                                      */
/*                       cannyEdgeImage                 */
/*                                                      */
/********************************************************/

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
void cannyEdgeImage(
           SrcIterator sul, SrcIterator slr, SrcAccessor sa,
           DestIterator dul, DestAccessor da,
           double scale, GradValue gradient_threshold, DestValue edge_marker)
{
    std::vector<Edgel> edgels;

    cannyEdgelList(sul, slr, sa, edgels, scale);

    for(unsigned int i=0; i<edgels.size(); ++i)
    {
        if(gradient_threshold < edgels[i].strength)
        {
            Diff2D pix((int)(edgels[i].x + 0.5), (int)(edgels[i].y + 0.5));

            da.set(edge_marker, dul, pix);
        }
    }
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
inline void cannyEdgeImage(
           triple<SrcIterator, SrcIterator, SrcAccessor> src,
           pair<DestIterator, DestAccessor> dest,
           double scale, GradValue gradient_threshold, DestValue edge_marker)
{
    cannyEdgeImage(src.first, src.second, src.third,
                   dest.first, dest.second,
                   scale, gradient_threshold, edge_marker);
}

/********************************************************/

namespace detail {

template <class DestIterator>
int neighborhoodConfiguration(DestIterator dul)
{
    int v = 0;
    NeighborhoodCirculator<DestIterator, EightNeighborCode> c(dul, EightNeighborCode::SouthEast);
    for(int i=0; i<8; ++i, --c)
    {
        v = (v << 1) | ((*c != 0) ? 1 : 0);
    }

    return v;
}

template <class GradValue>
struct SimplePoint
{
    Diff2D point;
    GradValue grad;

    SimplePoint(Diff2D const & p, GradValue g)
    : point(p), grad(g)
    {}

    bool operator<(SimplePoint const & o) const
    {
        return grad < o.grad;
    }

    bool operator>(SimplePoint const & o) const
    {
        return grad > o.grad;
    }
};

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
void cannyEdgeImageFromGrad(
           SrcIterator sul, SrcIterator slr, SrcAccessor grad,
           DestIterator dul, DestAccessor da,
           GradValue gradient_threshold, DestValue edge_marker)
{
    typedef typename SrcAccessor::value_type PixelType;
    typedef typename NormTraits<PixelType>::SquaredNormType NormType;

    NormType zero = NumericTraits<NormType>::zero();
    double tan22_5 = M_SQRT2 - 1.0;
    typename NormTraits<GradValue>::SquaredNormType g2thresh = squaredNorm(gradient_threshold);

    int w = slr.x - sul.x;
    int h = slr.y - sul.y;

    sul += Diff2D(1,1);
    dul += Diff2D(1,1);
    Diff2D p(0,0);

    for(int y = 1; y < h-1; ++y, ++sul.y, ++dul.y)
    {
        SrcIterator sx = sul;
        DestIterator dx = dul;
        for(int x = 1; x < w-1; ++x, ++sx.x, ++dx.x)
        {
            PixelType g = grad(sx);
            NormType g2n = squaredNorm(g);
            if(g2n < g2thresh)
                continue;

            NormType g2n1, g2n3;
            // find out quadrant
            if(abs(g[1]) < tan22_5*abs(g[0]))
            {
                // north-south edge
                g2n1 = squaredNorm(grad(sx, Diff2D(-1, 0)));
                g2n3 = squaredNorm(grad(sx, Diff2D(1, 0)));
            }
            else if(abs(g[0]) < tan22_5*abs(g[1]))
            {
                // west-east edge
                g2n1 = squaredNorm(grad(sx, Diff2D(0, -1)));
                g2n3 = squaredNorm(grad(sx, Diff2D(0, 1)));
            }
            else if(g[0]*g[1] < zero)
            {
                // north-west-south-east edge
                g2n1 = squaredNorm(grad(sx, Diff2D(1, -1)));
                g2n3 = squaredNorm(grad(sx, Diff2D(-1, 1)));
            }
            else
            {
                // north-east-south-west edge
                g2n1 = squaredNorm(grad(sx, Diff2D(-1, -1)));
                g2n3 = squaredNorm(grad(sx, Diff2D(1, 1)));
            }

            if(g2n1 < g2n && g2n3 <= g2n)
            {
                da.set(edge_marker, dx);
            }
        }
    }
}

} // namespace detail

/********************************************************/
/*                                                      */
/*              cannyEdgeImageWithThinning              */
/*                                                      */
/********************************************************/

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
void cannyEdgeImageFromGradWithThinning(
           SrcIterator sul, SrcIterator slr, SrcAccessor sa,
           DestIterator dul, DestAccessor da,
           GradValue gradient_threshold,
           DestValue edge_marker, bool addBorder)
{
    int w = slr.x - sul.x;
    int h = slr.y - sul.y;

    BImage edgeImage(w, h, BImage::value_type(0));
    BImage::traverser eul = edgeImage.upperLeft();
    BImage::Accessor ea = edgeImage.accessor();
    if(addBorder)
        initImageBorder(destImageRange(edgeImage), 1, 1);
    detail::cannyEdgeImageFromGrad(sul, slr, sa, eul, ea, gradient_threshold, 1);

    static bool isSimplePoint[256] = {
        0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,
        0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
        1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
        0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1,
        0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0,
        0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0,
        1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1,
        0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
        0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0,
        1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0,
        0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1,
        1, 0, 1, 0 };

    eul += Diff2D(1,1);
    sul += Diff2D(1,1);
    int w2 = w-2;
    int h2 = h-2;

    typedef detail::SimplePoint<GradValue> SP;
    // use std::greater becaus we need the smallest gradients at the top of the queue
    std::priority_queue<SP, std::vector<SP>, std::greater<SP> >  pqueue;

    Diff2D p(0,0);
    for(; p.y < h2; ++p.y)
    {
        for(p.x = 0; p.x < w2; ++p.x)
        {
            BImage::traverser e = eul + p;
            if(*e == 0)
                continue;
            int v = detail::neighborhoodConfiguration(e);
            if(isSimplePoint[v])
            {
                pqueue.push(SP(p, norm(sa(sul+p))));
                *e = 2; // remember that it is already in queue
            }
        }
    }

    static const Diff2D dist[] = { Diff2D(-1,0), Diff2D(0,-1),
                                   Diff2D(1,0),  Diff2D(0,1) };

    while(pqueue.size())
    {
        p = pqueue.top().point;
        pqueue.pop();

        BImage::traverser e = eul + p;
        int v = detail::neighborhoodConfiguration(e);
        if(!isSimplePoint[v])
            continue; // point may no longer be simple because its neighbors changed

        *e = 0; // delete simple point

        for(int i=0; i<4; ++i)
        {
            Diff2D pneu = p + dist[i];
            if(pneu.x == -1 || pneu.y == -1 || pneu.x == w2 || pneu.y == h2)
                continue; // do not remove points at the border

            BImage::traverser eneu = eul + pneu;
            if(*eneu == 1) // point is boundary and not yet in the queue
            {
                int v = detail::neighborhoodConfiguration(eneu);
                if(isSimplePoint[v])
                {
                    pqueue.push(SP(pneu, norm(sa(sul+pneu))));
                    *eneu = 2; // remember that it is already in queue
                }
            }
        }
    }

    initImageIf(destIterRange(dul, dul+Diff2D(w,h), da),
                maskImage(edgeImage), edge_marker);
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
inline void cannyEdgeImageFromGradWithThinning(
           triple<SrcIterator, SrcIterator, SrcAccessor> src,
           pair<DestIterator, DestAccessor> dest,
           GradValue gradient_threshold,
           DestValue edge_marker, bool addBorder)
{
    cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third,
                               dest.first, dest.second,
                               gradient_threshold, edge_marker, addBorder);
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
inline void cannyEdgeImageFromGradWithThinning(
           SrcIterator sul, SrcIterator slr, SrcAccessor sa,
           DestIterator dul, DestAccessor da,
           GradValue gradient_threshold, DestValue edge_marker)
{
    cannyEdgeImageFromGradWithThinning(sul, slr, sa,
                               dul, da,
                               gradient_threshold, edge_marker, true);
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
inline void cannyEdgeImageFromGradWithThinning(
           triple<SrcIterator, SrcIterator, SrcAccessor> src,
           pair<DestIterator, DestAccessor> dest,
           GradValue gradient_threshold, DestValue edge_marker)
{
    cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third,
                               dest.first, dest.second,
                               gradient_threshold, edge_marker, true);
}

/********************************************************/
/*                                                      */
/*              cannyEdgeImageWithThinning              */
/*                                                      */
/********************************************************/

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
void cannyEdgeImageWithThinning(
           SrcIterator sul, SrcIterator slr, SrcAccessor sa,
           DestIterator dul, DestAccessor da,
           double scale, GradValue gradient_threshold,
           DestValue edge_marker, bool addBorder)
{
    // mark pixels that are higher than their neighbors in gradient direction
    typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    BasicImage<TinyVector<TmpType, 2> > grad(slr-sul);
    gaussianGradient(srcIterRange(sul, slr, sa), destImage(grad), scale);
    cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destIter(dul, da),
                               gradient_threshold, edge_marker, addBorder);
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
inline void cannyEdgeImageWithThinning(
           triple<SrcIterator, SrcIterator, SrcAccessor> src,
           pair<DestIterator, DestAccessor> dest,
           double scale, GradValue gradient_threshold,
           DestValue edge_marker, bool addBorder)
{
    cannyEdgeImageWithThinning(src.first, src.second, src.third,
                               dest.first, dest.second,
                               scale, gradient_threshold, edge_marker, addBorder);
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
inline void cannyEdgeImageWithThinning(
           SrcIterator sul, SrcIterator slr, SrcAccessor sa,
           DestIterator dul, DestAccessor da,
           double scale, GradValue gradient_threshold, DestValue edge_marker)
{
    cannyEdgeImageWithThinning(sul, slr, sa,
                               dul, da,
                               scale, gradient_threshold, edge_marker, true);
}

template <class SrcIterator, class SrcAccessor,
          class DestIterator, class DestAccessor,
          class GradValue, class DestValue>
inline void cannyEdgeImageWithThinning(
           triple<SrcIterator, SrcIterator, SrcAccessor> src,
           pair<DestIterator, DestAccessor> dest,
           double scale, GradValue gradient_threshold, DestValue edge_marker)
{
    cannyEdgeImageWithThinning(src.first, src.second, src.third,
                               dest.first, dest.second,
                               scale, gradient_threshold, edge_marker, true);
}

/********************************************************/

template <class Image1, class Image2, class BackInsertable>
void internalCannyFindEdgels3x3(Image1 const & grad,
                                Image2 const & mask,
                                BackInsertable & edgels)
{
    typedef typename Image1::value_type PixelType;
    typedef typename PixelType::value_type ValueType;

    for(int y=1; y<grad.height()-1; ++y)
    {
        for(int x=1; x<grad.width()-1; ++x)
        {
            if(!mask(x,y))
                continue;

            ValueType gradx = grad(x,y)[0];
            ValueType grady = grad(x,y)[1];
            double mag = hypot(gradx, grady),
                   c = gradx / mag,
                   s = grady / mag;

            Matrix<double> ml(3,3), mr(3,1), l(3,1), r(3,1);
            l(0,0) = 1.0;

            for(int yy = -1; yy <= 1; ++yy)
            {
                for(int xx = -1; xx <= 1; ++xx)
                {
                    double u = c*xx + s*yy;
                    double v = norm(grad(x+xx, y+yy));
                    l(1,0) = u;
                    l(2,0) = u*u;
                    ml += outer(l);
                    mr += v*l;
                }
            }

            linearSolve(ml, mr, r);

            Edgel edgel;

            // local maximum => quadratic interpolation of sub-pixel location
            ValueType del = -r(1,0) / 2.0 / r(2,0);
            edgel.x = x + c*del;
            edgel.y = y + s*del;
            edgel.strength = mag;
            double orientation = VIGRA_CSTD::atan2(-grady, gradx) - M_PI * 1.5;
            if(orientation < 0.0)
                orientation += 2.0*M_PI;
            edgel.orientation = orientation;
            edgels.push_back(edgel);
        }
    }
}


/********************************************************/
/*                                                      */
/*                   cannyEdgelList3x3                  */
/*                                                      */
/********************************************************/

template <class SrcIterator, class SrcAccessor, class BackInsertable>
void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
                        BackInsertable & edgels, double scale)
{
    int w = lr.x - ul.x;
    int h = lr.y - ul.y;

    typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
    gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);

    UInt8Image edges(lr-ul);
    cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destImage(edges),
                                       0.0, 1, false);

    // find edgels
    internalCannyFindEdgels3x3(grad, edges, edgels);
}

template <class SrcIterator, class SrcAccessor, class BackInsertable>
inline void
cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
               BackInsertable & edgels, double scale)
{
    cannyEdgelList3x3(src.first, src.second, src.third, edgels, scale);
}




} // namespace vigra

#endif // VIGRA_EDGEDETECTION_HXX

---