--- a
+++ b/src/hugin_base/algorithms/control_points/edgedetection.hxx
@@ -0,0 +1,2181 @@
+/************************************************************************/
+/*                                                                      */
+/*               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 {
+
+/** \addtogroup EdgeDetection Edge Detection
+    Edge detectors based on first and second derivatives,
+          and related post-processing.
+*/
+//@{
+
+/********************************************************/
+/*                                                      */
+/*           differenceOfExponentialEdgeImage           */
+/*                                                      */
+/********************************************************/
+
+/** \brief Detect and mark edges in an edge image using the Shen/Castan zero-crossing detector.
+
+    This operator applies an exponential filter to the source image
+    at the given <TT>scale</TT> and subtracts the result from the original image.
+    Zero crossings are detected in the resulting difference image. Whenever the
+    gradient at a zero crossing is greater than the given <TT>gradient_threshold</TT>,
+    an edge point is marked (using <TT>edge_marker</TT>) in the destination image on
+    the darker side of the zero crossing (note that zero crossings occur
+    <i>between</i> pixels). For example:
+
+    \code
+    sign of difference image     resulting edge points (*)
+
+        + - -                          * * .
+        + + -               =>         . * *
+        + + +                          . . .
+    \endcode
+
+    Non-edge pixels (<TT>.</TT>) remain untouched in the destination image.
+    The result can be improved by the post-processing operation \ref removeShortEdges().
+    A more accurate edge placement can be achieved with the function
+    \ref differenceOfExponentialCrackEdgeImage().
+
+    The source value type
+    (<TT>SrcAccessor::value_type</TT>) must be a linear algebra, i.e. addition,
+    subtraction and multiplication of the type with itself, and multiplication
+    with double and
+    \ref NumericTraits "NumericTraits" must
+    be defined. In addition, this type must be less-comparable.
+
+    <b> Declarations:</b>
+
+    pass arguments explicitly:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor,
+              class DestIterator, class DestAccessor,
+              class GradValue,
+              class DestValue = DestAccessor::value_type>
+        void differenceOfExponentialEdgeImage(
+               SrcIterator sul, SrcIterator slr, SrcAccessor sa,
+               DestIterator dul, DestAccessor da,
+               double scale, GradValue gradient_threshold,
+               DestValue edge_marker = NumericTraits<DestValue>::one())
+    }
+    \endcode
+
+    use argument objects in conjunction with \ref ArgumentObjectFactories:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor,
+              class DestIterator, class DestAccessor,
+              class GradValue,
+              class DestValue = DestAccessor::value_type>
+        inline
+        void differenceOfExponentialEdgeImage(
+               triple<SrcIterator, SrcIterator, SrcAccessor> src,
+               pair<DestIterator, DestAccessor> dest,
+               double scale, GradValue gradient_threshold,
+               DestValue edge_marker = NumericTraits<DestValue>::one())
+    }
+    \endcode
+
+    <b> Usage:</b>
+
+        <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
+    Namespace: vigra
+
+    \code
+    vigra::BImage src(w,h), edges(w,h);
+
+    // empty edge image
+    edges = 0;
+    ...
+
+    // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
+    vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges),
+                                     0.8, 4.0, 1);
+    \endcode
+
+    <b> Required Interface:</b>
+
+    \code
+    SrcImageIterator src_upperleft, src_lowerright;
+    DestImageIterator dest_upperleft;
+
+    SrcAccessor src_accessor;
+    DestAccessor dest_accessor;
+
+    SrcAccessor::value_type u = src_accessor(src_upperleft);
+    double d;
+    GradValue gradient_threshold;
+
+    u = u + u
+    u = u - u
+    u = u * u
+    u = d * u
+    u < gradient_threshold
+
+    DestValue edge_marker;
+    dest_accessor.set(edge_marker, dest_upperleft);
+    \endcode
+
+    <b> Preconditions:</b>
+
+    \code
+    scale > 0
+    gradient_threshold > 0
+    \endcode
+*/
+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         */
+/*                                                      */
+/********************************************************/
+
+/** \brief Detect and mark edges in a crack edge image using the Shen/Castan zero-crossing detector.
+
+    This operator applies an exponential filter to the source image
+    at the given <TT>scale</TT> and subtracts the result from the original image.
+    Zero crossings are detected in the resulting difference image. Whenever the
+    gradient at a zero crossing is greater than the given <TT>gradient_threshold</TT>,
+    an edge point is marked (using <TT>edge_marker</TT>) in the destination image
+    <i>between</i> the corresponding original pixels. Topologically, this means we
+    must insert additional pixels between the original ones to represent the
+    boundaries between the pixels (the so called zero- and one-cells, with the original
+    pixels being two-cells). Within VIGRA, such an image is called \ref CrackEdgeImage.
+    To allow insertion of the zero- and one-cells, the destination image must have twice the
+    size of the original (precisely, <TT>(2*w-1)</TT> by <TT>(2*h-1)</TT> pixels). Then the algorithm
+    proceeds as follows:
+
+    \code
+sign of difference image     insert zero- and one-cells     resulting edge points (*)
+
+                                     + . - . -                   . * . . .
+      + - -                          . . . . .                   . * * * .
+      + + -               =>         + . + . -           =>      . . . * .
+      + + +                          . . . . .                   . . . * *
+                                     + . + . +                   . . . . .
+    \endcode
+
+    Thus the edge points are marked where they actually are - in between the pixels.
+    An important property of the resulting edge image is that it conforms to the notion
+    of well-composedness as defined by Latecki et al., i.e. connected regions and edges
+    obtained by a subsequent \ref Labeling do not depend on
+    whether 4- or 8-connectivity is used.
+    The non-edge pixels (<TT>.</TT>) in the destination image remain unchanged.
+    The result conformes to the requirements of a \ref CrackEdgeImage. It can be further
+    improved by the post-processing operations \ref removeShortEdges() and
+    \ref closeGapsInCrackEdgeImage().
+
+    The source value type (<TT>SrcAccessor::value_type</TT>) must be a linear algebra, i.e. addition,
+    subtraction and multiplication of the type with itself, and multiplication
+    with double and
+    \ref NumericTraits "NumericTraits" must
+    be defined. In addition, this type must be less-comparable.
+
+    <b> Declarations:</b>
+
+    pass arguments explicitly:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor,
+              class DestIterator, class DestAccessor,
+              class GradValue,
+              class DestValue = DestAccessor::value_type>
+        void differenceOfExponentialCrackEdgeImage(
+               SrcIterator sul, SrcIterator slr, SrcAccessor sa,
+               DestIterator dul, DestAccessor da,
+               double scale, GradValue gradient_threshold,
+               DestValue edge_marker = NumericTraits<DestValue>::one())
+    }
+    \endcode
+
+    use argument objects in conjunction with \ref ArgumentObjectFactories:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor,
+              class DestIterator, class DestAccessor,
+              class GradValue,
+              class DestValue = DestAccessor::value_type>
+        inline
+        void differenceOfExponentialCrackEdgeImage(
+               triple<SrcIterator, SrcIterator, SrcAccessor> src,
+               pair<DestIterator, DestAccessor> dest,
+               double scale, GradValue gradient_threshold,
+               DestValue edge_marker = NumericTraits<DestValue>::one())
+    }
+    \endcode
+
+    <b> Usage:</b>
+
+        <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
+    Namespace: vigra
+
+    \code
+    vigra::BImage src(w,h), edges(2*w-1,2*h-1);
+
+    // empty edge image
+    edges = 0;
+    ...
+
+    // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
+    vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
+                                     0.8, 4.0, 1);
+    \endcode
+
+    <b> Required Interface:</b>
+
+    \code
+    SrcImageIterator src_upperleft, src_lowerright;
+    DestImageIterator dest_upperleft;
+
+    SrcAccessor src_accessor;
+    DestAccessor dest_accessor;
+
+    SrcAccessor::value_type u = src_accessor(src_upperleft);
+    double d;
+    GradValue gradient_threshold;
+
+    u = u + u
+    u = u - u
+    u = u * u
+    u = d * u
+    u < gradient_threshold
+
+    DestValue edge_marker;
+    dest_accessor.set(edge_marker, dest_upperleft);
+    \endcode
+
+    <b> Preconditions:</b>
+
+    \code
+    scale > 0
+    gradient_threshold > 0
+    \endcode
+
+    The destination image must have twice the size of the source:
+    \code
+    w_dest = 2 * w_src - 1
+    h_dest = 2 * h_src - 1
+    \endcode
+*/
+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                    */
+/*                                                      */
+/********************************************************/
+
+/** \brief Remove short edges from an edge image.
+
+    This algorithm can be applied as a post-processing operation of
+    \ref differenceOfExponentialEdgeImage() and \ref differenceOfExponentialCrackEdgeImage().
+    It removes all edges that are shorter than <TT>min_edge_length</TT>. The corresponding
+    pixels are set to the <TT>non_edge_marker</TT>. The idea behind this algorithms is
+    that very short edges are probably caused by noise and don't represent interesting
+    image structure. Technically, the algorithms executes a connected components labeling,
+    so the image's value type must be equality comparable.
+
+    If the source image fulfills the requirements of a \ref CrackEdgeImage,
+    it will still do so after application of this algorithm.
+
+    Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place,
+    i.e. on only one image. Also, the algorithm assumes that all non-edges pixels are already
+    marked with the given <TT>non_edge_marker</TT> value.
+
+    <b> Declarations:</b>
+
+    pass arguments explicitly:
+    \code
+    namespace vigra {
+        template <class Iterator, class Accessor, class SrcValue>
+        void removeShortEdges(
+               Iterator sul, Iterator slr, Accessor sa,
+               int min_edge_length, SrcValue non_edge_marker)
+    }
+    \endcode
+
+    use argument objects in conjunction with \ref ArgumentObjectFactories:
+    \code
+    namespace vigra {
+        template <class Iterator, class Accessor, class SrcValue>
+        inline
+        void removeShortEdges(
+               triple<Iterator, Iterator, Accessor> src,
+               int min_edge_length, SrcValue non_edge_marker)
+    }
+    \endcode
+
+    <b> Usage:</b>
+
+        <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
+    Namespace: vigra
+
+    \code
+    vigra::BImage src(w,h), edges(w,h);
+
+    // empty edge image
+    edges = 0;
+    ...
+
+    // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
+    vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges),
+                                     0.8, 4.0, 1);
+
+    // zero edges shorter than 10 pixels
+    vigra::removeShortEdges(srcImageRange(edges), 10, 0);
+    \endcode
+
+    <b> Required Interface:</b>
+
+    \code
+    SrcImageIterator src_upperleft, src_lowerright;
+    DestImageIterator dest_upperleft;
+
+    SrcAccessor src_accessor;
+    DestAccessor dest_accessor;
+
+    SrcAccessor::value_type u = src_accessor(src_upperleft);
+
+    u == u
+
+    SrcValue non_edge_marker;
+    src_accessor.set(non_edge_marker, src_upperleft);
+    \endcode
+*/
+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                */
+/*                                                      */
+/********************************************************/
+
+/** \brief Close one-pixel wide gaps in a cell grid edge image.
+
+    This algorithm is typically applied as a post-processing operation of
+    \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
+    the requirements of a \ref CrackEdgeImage, and will still do so after
+    application of this algorithm.
+
+    It closes one pixel wide gaps in the edges resulting from this algorithm.
+    Since these gaps are usually caused by zero crossing slightly below the gradient
+    threshold used in edge detection, this algorithms acts like a weak hysteresis
+    thresholding. The newly found edge pixels are marked with the given <TT>edge_marker</TT>.
+    The image's value type must be equality comparable.
+
+    Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place,
+    i.e. on only one image.
+
+    <b> Declarations:</b>
+
+    pass arguments explicitly:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor, class SrcValue>
+        void closeGapsInCrackEdgeImage(
+               SrcIterator sul, SrcIterator slr, SrcAccessor sa,
+               SrcValue edge_marker)
+    }
+    \endcode
+
+    use argument objects in conjunction with \ref ArgumentObjectFactories:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor, class SrcValue>
+        inline
+        void closeGapsInCrackEdgeImage(
+               triple<SrcIterator, SrcIterator, SrcAccessor> src,
+               SrcValue edge_marker)
+    }
+    \endcode
+
+    <b> Usage:</b>
+
+        <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
+    Namespace: vigra
+
+    \code
+    vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
+
+    // empty edge image
+    edges = 0;
+    ...
+
+    // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
+    vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
+                                         0.8, 4.0, 1);
+
+    // close gaps, mark with 1
+    vigra::closeGapsInCrackEdgeImage(srcImageRange(edges), 1);
+
+    // zero edges shorter than 20 pixels
+    vigra::removeShortEdges(srcImageRange(edges), 10, 0);
+    \endcode
+
+    <b> Required Interface:</b>
+
+    \code
+    SrcImageIterator src_upperleft, src_lowerright;
+
+    SrcAccessor src_accessor;
+    DestAccessor dest_accessor;
+
+    SrcAccessor::value_type u = src_accessor(src_upperleft);
+
+    u == u
+    u != u
+
+    SrcValue edge_marker;
+    src_accessor.set(edge_marker, src_upperleft);
+    \endcode
+*/
+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                  */
+/*                                                      */
+/********************************************************/
+
+/** \brief Beautify crack edge image for visualization.
+
+    This algorithm is applied as a post-processing operation of
+    \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
+    the requirements of a \ref CrackEdgeImage, but will <b> not</b> do so after
+    application of this algorithm. In particular, the algorithm removes zero-cells
+    marked as edges to avoid staircase effects on diagonal lines like this:
+
+    \code
+    original edge points (*)     resulting edge points
+
+          . * . . .                   . * . . .
+          . * * * .                   . . * . .
+          . . . * .           =>      . . . * .
+          . . . * *                   . . . . *
+          . . . . .                   . . . . .
+    \endcode
+
+    Therfore, this algorithm should only be applied as a vizualization aid, i.e.
+    for human inspection. The algorithm assumes that edges are marked with <TT>edge_marker</TT>,
+    and background pixels with <TT>background_marker</TT>. The image's value type must be
+    equality comparable.
+
+    Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place,
+    i.e. on only one image.
+
+    <b> Declarations:</b>
+
+    pass arguments explicitly:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor, class SrcValue>
+        void beautifyCrackEdgeImage(
+               SrcIterator sul, SrcIterator slr, SrcAccessor sa,
+               SrcValue edge_marker, SrcValue background_marker)
+    }
+    \endcode
+
+    use argument objects in conjunction with \ref ArgumentObjectFactories:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor, class SrcValue>
+        inline
+        void beautifyCrackEdgeImage(
+                   triple<SrcIterator, SrcIterator, SrcAccessor> src,
+               SrcValue edge_marker, SrcValue background_marker)
+    }
+    \endcode
+
+    <b> Usage:</b>
+
+        <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
+    Namespace: vigra
+
+    \code
+    vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
+
+    // empty edge image
+    edges = 0;
+    ...
+
+    // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
+    vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
+                                         0.8, 4.0, 1);
+
+    // beautify edge image for visualization
+    vigra::beautifyCrackEdgeImage(destImageRange(edges), 1, 0);
+
+    // show to the user
+    window.open(edges);
+    \endcode
+
+    <b> Required Interface:</b>
+
+    \code
+    SrcImageIterator src_upperleft, src_lowerright;
+
+    SrcAccessor src_accessor;
+    DestAccessor dest_accessor;
+
+    SrcAccessor::value_type u = src_accessor(src_upperleft);
+
+    u == u
+    u != u
+
+    SrcValue background_marker;
+    src_accessor.set(background_marker, src_upperleft);
+    \endcode
+*/
+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);
+}
+
+
+/** Helper class that stores edgel attributes.
+*/
+class Edgel
+{
+  public:
+        /** The edgel's sub-pixel x coordinate.
+        */
+    float x;
+
+        /** The edgel's sub-pixel y coordinate.
+        */
+    float y;
+
+        /** The edgel's strength (magnitude of the gradient vector).
+        */
+    float strength;
+
+        /**
+        The edgel's orientation. This is the angle
+        between the x-axis and the edge, so that the bright side of the
+        edge is on the right. The angle is measured
+        counter-clockwise in radians like this:
+
+
+        \code
+
+  edgel axis
+      \  (bright side)
+ (dark \
+ side)  \ /__
+         \\  \ orientation angle
+          \  |
+           +------------> x-axis
+           |
+           |
+           |
+           |
+    y-axis V
+        \endcode
+
+        So, for example a vertical edge with its dark side on the left
+        has orientation PI/2, and a horizontal edge with dark side on top
+        has orientation 0. Obviously, the edge's orientation changes
+        by PI if the contrast is reversed.
+
+        */
+    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                  */
+/*                                                      */
+/********************************************************/
+
+/** \brief Simple implementation of Canny's edge detector.
+
+    This operator first calculates the gradient vector for each
+    pixel of the image using first derivatives of a Gaussian at the
+    given scale. Then a very simple non-maxima supression is performed:
+    for each 3x3 neighborhood, it is determined whether the center pixel has
+    larger gradient magnitude than its two neighbors in gradient direction
+    (where the direction is rounded into octands). If this is the case,
+    a new \ref Edgel is appended to the given vector of <TT>edgels</TT>. The subpixel
+    edgel position is determined by fitting a parabola
+    to the three gradient magnitude values
+    mentioned above. The sub-pixel location of the parabola's tip
+    and the gradient magnitude and direction are written in the newly created edgel.
+
+    <b> Declarations:</b>
+
+    pass arguments explicitly:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor, class BackInsertable>
+        void cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
+                            BackInsertable & edgels, double scale);
+    }
+    \endcode
+
+    use argument objects in conjunction with \ref ArgumentObjectFactories:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor, class BackInsertable>
+        void
+        cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
+                       BackInsertable & edgels, double scale);
+    }
+    \endcode
+
+    <b> Usage:</b>
+
+    <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
+    Namespace: vigra
+
+    \code
+    vigra::BImage src(w,h);
+
+    // empty edgel list
+    std::vector<vigra::Edgel> edgels;
+    ...
+
+    // find edgels at scale 0.8
+    vigra::cannyEdgelList(srcImageRange(src), edgels, 0.8);
+    \endcode
+
+    <b> Required Interface:</b>
+
+    \code
+    SrcImageIterator src_upperleft;
+    SrcAccessor src_accessor;
+
+    src_accessor(src_upperleft);
+
+    BackInsertable edgels;
+    edgels.push_back(Edgel());
+    \endcode
+
+    SrcAccessor::value_type must be a type convertible to float
+
+    <b> Preconditions:</b>
+
+    \code
+    scale > 0
+    \endcode
+*/
+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                 */
+/*                                                      */
+/********************************************************/
+
+/** \brief Detect and mark edges in an edge image using Canny's algorithm.
+
+    This operator first calls \ref cannyEdgelList() to generate an
+    edgel list for the given image. Then it scans this list and selects edgels
+    whose strength is above the given <TT>gradient_threshold</TT>. For each of these
+    edgels, the edgel's location is rounded to the nearest pixel, and that
+    pixel marked with the given <TT>edge_marker</TT>.
+
+    <b> Declarations:</b>
+
+    pass arguments explicitly:
+    \code
+    namespace vigra {
+        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);
+    }
+    \endcode
+
+    use argument objects in conjunction with \ref ArgumentObjectFactories:
+    \code
+    namespace vigra {
+        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);
+    }
+    \endcode
+
+    <b> Usage:</b>
+
+    <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
+    Namespace: vigra
+
+    \code
+    vigra::BImage src(w,h), edges(w,h);
+
+    // empty edge image
+    edges = 0;
+    ...
+
+    // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
+    vigra::cannyEdgeImage(srcImageRange(src), destImage(edges),
+                                     0.8, 4.0, 1);
+    \endcode
+
+    <b> Required Interface:</b>
+
+    see also: \ref cannyEdgelList().
+
+    \code
+    DestImageIterator dest_upperleft;
+    DestAccessor dest_accessor;
+    DestValue edge_marker;
+
+    dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
+    \endcode
+
+    <b> Preconditions:</b>
+
+    \code
+    scale > 0
+    gradient_threshold > 0
+    \endcode
+*/
+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              */
+/*                                                      */
+/********************************************************/
+
+/** \brief Detect and mark edges in an edge image using Canny's algorithm.
+
+    The input pixels of this algorithms must be vectors of length 2 (see Required Interface below).
+    It first searches for all pixels whose gradient magnitude is larger
+    than the given <tt>gradient_threshold</tt> and larger than the magnitude of its two neighbors
+    in gradient direction (where these neighbors are determined by nearest neighbor
+    interpolation, i.e. according to the octant where the gradient points into).
+    The resulting edge pixel candidates are then subjected to topological thinning
+    so that the remaining edge pixels can be linked into edgel chains with a provable,
+    non-heuristic algorithm. Thinning is performed so that the pixels with highest gradient
+    magnitude survive. Optionally, the outermost pixels are marked as edge pixels
+    as well when <tt>addBorder</tt> is true. The remaining pixels will be marked in the destination
+    image with the value of <tt>edge_marker</tt> (all non-edge pixels remain untouched).
+
+    <b> Declarations:</b>
+
+    pass arguments explicitly:
+    \code
+    namespace vigra {
+        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 = true);
+    }
+    \endcode
+
+    use argument objects in conjunction with \ref ArgumentObjectFactories:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor,
+                  class DestIterator, class DestAccessor,
+                  class GradValue, class DestValue>
+        void cannyEdgeImageFromGradWithThinning(
+                   triple<SrcIterator, SrcIterator, SrcAccessor> src,
+                   pair<DestIterator, DestAccessor> dest,
+                   GradValue gradient_threshold,
+                   DestValue edge_marker, bool addBorder = true);
+    }
+    \endcode
+
+    <b> Usage:</b>
+
+    <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
+    Namespace: vigra
+
+    \code
+    vigra::BImage src(w,h), edges(w,h);
+
+    vigra::FVector2Image grad(w,h);
+    // compute the image gradient at scale 0.8
+    vigra::gaussianGradient(srcImageRange(src), destImage(grad), 0.8);
+
+    // empty edge image
+    edges = 0;
+    // find edges gradient larger than 4.0, mark with 1, and add border
+    vigra::cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destImage(edges),
+                                              4.0, 1, true);
+    \endcode
+
+    <b> Required Interface:</b>
+
+    \code
+    // the input pixel type must be a vector with two elements
+    SrcImageIterator src_upperleft;
+    SrcAccessor src_accessor;
+    typedef SrcAccessor::value_type SrcPixel;
+    typedef NormTraits<SrcPixel>::SquaredNormType SrcSquaredNormType;
+
+    SrcPixel g = src_accessor(src_upperleft);
+    SrcPixel::value_type g0 = g[0];
+    SrcSquaredNormType gn = squaredNorm(g);
+
+    DestImageIterator dest_upperleft;
+    DestAccessor dest_accessor;
+    DestValue edge_marker;
+
+    dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
+    \endcode
+
+    <b> Preconditions:</b>
+
+    \code
+    gradient_threshold > 0
+    \endcode
+*/
+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              */
+/*                                                      */
+/********************************************************/
+
+/** \brief Detect and mark edges in an edge image using Canny's algorithm.
+
+    This operator first calls \ref gaussianGradient() to compute the gradient of the input
+    image, ad then \ref cannyEdgeImageFromGradWithThinning() to generate an
+    edge image. See there for more detailed documentation.
+
+    <b> Declarations:</b>
+
+    pass arguments explicitly:
+    \code
+    namespace vigra {
+        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 = true);
+    }
+    \endcode
+
+    use argument objects in conjunction with \ref ArgumentObjectFactories:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor,
+                  class DestIterator, class DestAccessor,
+                  class GradValue, class DestValue>
+        void cannyEdgeImageWithThinning(
+                   triple<SrcIterator, SrcIterator, SrcAccessor> src,
+                   pair<DestIterator, DestAccessor> dest,
+                   double scale, GradValue gradient_threshold,
+                   DestValue edge_marker, bool addBorder = true);
+    }
+    \endcode
+
+    <b> Usage:</b>
+
+    <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
+    Namespace: vigra
+
+    \code
+    vigra::BImage src(w,h), edges(w,h);
+
+    // empty edge image
+    edges = 0;
+    ...
+
+    // find edges at scale 0.8 with gradient larger than 4.0, mark with 1, annd add border
+    vigra::cannyEdgeImageWithThinning(srcImageRange(src), destImage(edges),
+                                     0.8, 4.0, 1, true);
+    \endcode
+
+    <b> Required Interface:</b>
+
+    see also: \ref cannyEdgelList().
+
+    \code
+    DestImageIterator dest_upperleft;
+    DestAccessor dest_accessor;
+    DestValue edge_marker;
+
+    dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
+    \endcode
+
+    <b> Preconditions:</b>
+
+    \code
+    scale > 0
+    gradient_threshold > 0
+    \endcode
+*/
+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                  */
+/*                                                      */
+/********************************************************/
+
+/** \brief Improved implementation of Canny's edge detector.
+
+    This operator first computes pixels which are crossed by the edge using
+    cannyEdgeImageWithThinning(). The gradient magnitude in the 3x3 neighborhood of these
+    pixels are then projected onto the normal of the edge (as determined
+    by the gradient direction). The edgel's subpixel location is found by fitting a
+    parabola through the 9 gradient values and determining the parabola's tip.
+    A new \ref Edgel is appended to the given vector of <TT>edgels</TT>. Since the parabola
+    is fitted to 9 points rather than 3 points as in cannyEdgelList(), the accuracy is higher.
+
+    <b> Declarations:</b>
+
+    pass arguments explicitly:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor, class BackInsertable>
+        void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
+                               BackInsertable & edgels, double scale);
+    }
+    \endcode
+
+    use argument objects in conjunction with \ref ArgumentObjectFactories:
+    \code
+    namespace vigra {
+        template <class SrcIterator, class SrcAccessor, class BackInsertable>
+        void
+        cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
+                          BackInsertable & edgels, double scale);
+    }
+    \endcode
+
+    <b> Usage:</b>
+
+    <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
+    Namespace: vigra
+
+    \code
+    vigra::BImage src(w,h);
+
+    // empty edgel list
+    std::vector<vigra::Edgel> edgels;
+    ...
+
+    // find edgels at scale 0.8
+    vigra::cannyEdgelList3x3(srcImageRange(src), edgels, 0.8);
+    \endcode
+
+    <b> Required Interface:</b>
+
+    \code
+    SrcImageIterator src_upperleft;
+    SrcAccessor src_accessor;
+
+    src_accessor(src_upperleft);
+
+    BackInsertable edgels;
+    edgels.push_back(Edgel());
+    \endcode
+
+    SrcAccessor::value_type must be a type convertible to float
+
+    <b> Preconditions:</b>
+
+    \code
+    scale > 0
+    \endcode
+*/
+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);
+}
+
+
+
+//@}
+
+/** \page CrackEdgeImage Crack Edge Image
+
+Crack edges are marked <i>between</i> the pixels of an image.
+A Crack Edge Image is an image that represents these edges. In order
+to accomodate the cracks, the Crack Edge Image must be twice as large
+as the original image (precisely (2*w - 1) by (2*h - 1)). A Crack Edge Image
+can easily be derived from a binary image or from the signs of the
+response of a Laplacean filter. Consider the following sketch, where
+<TT>+</TT> encodes the foreground, <TT>-</TT> the background, and
+<TT>*</TT> the resulting crack edges.
+
+    \code
+sign of difference image         insert cracks         resulting CrackEdgeImage
+
+                                   + . - . -              . * . . .
+      + - -                        . . . . .              . * * * .
+      + + -               =>       + . + . -      =>      . . . * .
+      + + +                        . . . . .              . . . * *
+                                   + . + . +              . . . . .
+    \endcode
+
+Starting from the original binary image (left), we insert crack pixels
+to get to the double-sized image (center). Finally, we mark all
+crack pixels whose non-crack neighbors have different signs as
+crack edge points, while all other pixels (crack and non-crack) become
+region pixels.
+
+<b>Requirements on a Crack Edge Image:</b>
+
+<ul>
+    <li>Crack Edge Images have odd width and height.
+    <li>Crack pixels have at least one odd coordinate.
+    <li>Only crack pixels may be marked as edge points.
+    <li>Crack pixels with two odd coordinates must be marked as edge points
+        whenever any of their neighboring crack pixels was marked.
+</ul>
+
+The last two requirements ensure that both edges and regions are 4-connected.
+Thus, 4-connectivity and 8-connectivity yield identical connected
+components in a Crack Edge Image (so called <i>well-composedness</i>).
+This ensures that Crack Edge Images have nice topological properties
+(cf. L. J. Latecki: "Well-Composed Sets", Academic Press, 2000).
+*/
+
+
+} // namespace vigra
+
+#endif // VIGRA_EDGEDETECTION_HXX