Revision: 2947
http://hugin.svn.sourceforge.net/hugin/?rev=2947&view=rev
Author: dangelo
Date: 2008-03-14 01:50:14 -0700 (Fri, 14 Mar 2008)
Log Message:
-----------
Moved matchpoint into main hugin tree. Fixed crash when computing the descriptors.
Modified Paths:
--------------
hugin/trunk/src/CMakeLists.txt
hugin/trunk/src/matchpoint/CMakeLists.txt
hugin/trunk/src/matchpoint/Descriptor.cpp
hugin/trunk/src/matchpoint/MatchPoint.cpp
Added Paths:
-----------
hugin/trunk/src/matchpoint/
hugin/trunk/src/matchpoint/edgedetection.hxx
Removed Paths:
-------------
gsoc07_featuredetection/
Modified: hugin/trunk/src/CMakeLists.txt
===================================================================
--- hugin/trunk/src/CMakeLists.txt 2008-03-13 18:24:34 UTC (rev 2946)
+++ hugin/trunk/src/CMakeLists.txt 2008-03-14 08:50:14 UTC (rev 2947)
@@ -23,6 +23,7 @@
add_subdirectory(foreign)
add_subdirectory(hugin_base)
add_subdirectory(tools)
+add_subdirectory(matchpoint)
add_subdirectory(deghosting)
# build vips stuff, only if vips was found
Copied: hugin/trunk/src/matchpoint (from rev 2945, gsoc07_featuredetection)
Modified: hugin/trunk/src/matchpoint/CMakeLists.txt
===================================================================
--- gsoc07_featuredetection/CMakeLists.txt 2008-03-12 18:43:05 UTC (rev 2945)
+++ hugin/trunk/src/matchpoint/CMakeLists.txt 2008-03-14 08:50:14 UTC (rev 2947)
@@ -1,25 +1,10 @@
-PROJECT(MatchPoint)
-
-#if you don't want the full compiler output, remove the following line
-#SET(CMAKE_VERBOSE_MAKEFILE ON)
-
-# check for opencv
-#FIND_PATH(OPENCV_INCLUDE_DIR cv.h PATH_SUFFIXES opencv)
-#INCLUDE_DIRECTORIES(${OPENCV_INCLUDE_DIR})
-
-#add definitions, compiler switches, etc.
-ADD_DEFINITIONS(-Wall -O2 -g)
-
#list all source files here
-ADD_EXECUTABLE(MatchPoint MatchPoint.cpp APImage.cpp HessianDetector.cpp Descriptor.cpp edgedetection.hxx)
+ADD_EXECUTABLE(matchpoint MatchPoint.cpp APImage.cpp HessianDetector.cpp Descriptor.cpp)
#need to link to some other libraries ? just add them here
-TARGET_LINK_LIBRARIES(MatchPoint /usr/lib/libvigraimpex.so png jpeg tiff)
-
-#list all source files here
-#ADD_EXECUTABLE(MatchPoint.cpp APImage.cpp HessianDetector.cpp Descriptor.cpp FeatureMatcher.cpp)
-#ADD_EXECUTABLE(matchpoint APImage.cpp HessianDetector.cpp Descriptor.cpp)
-#need to link to some other libraries ? just add them here
-#TARGET_LINK_LIBRARIES(FeatureMatching highgui cv cxcore png jpeg)
-#TARGET_LINK_LIBRARIES(/usr/lib/libvigraimpex.so png jpeg tiff)
+TARGET_LINK_LIBRARIES(matchpoint ${image_libs})
+# matchpoint is not ready for general use yet
+#
+#install(TARGETS matchpoint
+# DESTINATION ${BINDIR})
Modified: hugin/trunk/src/matchpoint/Descriptor.cpp
===================================================================
--- gsoc07_featuredetection/Descriptor.cpp 2008-03-12 18:43:05 UTC (rev 2945)
+++ hugin/trunk/src/matchpoint/Descriptor.cpp 2008-03-14 08:50:14 UTC (rev 2947)
@@ -653,6 +653,12 @@
int yStart=(interestPoint[1]-regionSize);
int yEnd=(interestPoint[1]+regionSize);
+ // only analyse pixels in image
+ int xStartImg = xStart < 0 ? 0 : xStart;
+ int xEndImg = xEnd >= this->image->getWidthBW() ? this->image->getWidthBW() -1 : xEnd;
+ int yStartImg = yStart < 0 ? 0 : yStart;
+ int yEndImg = yEnd >= this->image->getHeightBW() ? this->image->getHeightBW() -1 : yEnd;
+
/*
double result;
result = atan (gy/gx) * 180 / PI;
@@ -676,8 +682,8 @@
iPointRelative[1]=regionSize;
// find edgels at scale of the interest point
- vigra::cannyEdgelList(srcIterRange( this->image->imgBW->upperLeft() + vigra::Diff2D(xStart, yStart),
- this->image->imgBW->upperLeft() + vigra::Diff2D(xEnd, yEnd)),
+ vigra::cannyEdgelList(srcIterRange( this->image->imgBW->upperLeft() + vigra::Diff2D(xStartImg, yStartImg),
+ this->image->imgBW->upperLeft() + vigra::Diff2D(xEndImg, yEndImg)),
edgels, this->_getMaxima(interestPoint[0], interestPoint[1]),iPointRelative);
/*this->image->_cannyEdgelList1(srcIterRange( this->image->imgBW->upperLeft() + vigra::Diff2D(xStart, yStart),
this->image->imgBW->upperLeft() + vigra::Diff2D(xEnd, yEnd)),
@@ -703,7 +709,8 @@
edgePoint=*iter2;
//TODO discard edgels that are not in circle
- double distance=_euclidianDistance(interestPoint[0],interestPoint[1],edgePoint.x+xStart, edgePoint.y+yStart);
+ double distance=_euclidianDistance(interestPoint[0],interestPoint[1],edgePoint.x+xStartImg, edgePoint.y+yStartImg);
+
if(distance>regionSize) {
iter2++; continue;
}
Modified: hugin/trunk/src/matchpoint/MatchPoint.cpp
===================================================================
--- gsoc07_featuredetection/MatchPoint.cpp 2008-03-12 18:43:05 UTC (rev 2945)
+++ hugin/trunk/src/matchpoint/MatchPoint.cpp 2008-03-14 08:50:14 UTC (rev 2947)
@@ -83,7 +83,7 @@
verbose = true;
break;
case 't':
- testFileOutput=true;
+ testFileOutput=true;
break;
/*case 'o':
output = optarg;
@@ -141,11 +141,12 @@
d1.createDescriptors();
if(testFileOutput) {
- string testOutputPath = input1.append(".key");
- if (verbose) cerr << "Generating output file for matlab test suite: " <<testOutputPath << endl;
- d1.printDescriptors(testOutputPath); //for matlab test suite
+// string testOutputPath = input1.append(".key");
+ if (verbose) cerr << "Generating output file for matlab test suite: " << output << endl;
+ d1.printDescriptors(output); //for matlab test suite
+ } else {
+ d1.generateAutopanoXML(output);
}
- d1.generateAutopanoXML(output);
return EXIT_SUCCESS;
}
Copied: hugin/trunk/src/matchpoint/edgedetection.hxx (from rev 2946, gsoc07_featuredetection/edgedetection.hxx)
===================================================================
--- hugin/trunk/src/matchpoint/edgedetection.hxx (rev 0)
+++ hugin/trunk/src/matchpoint/edgedetection.hxx 2008-03-14 08:50:14 UTC (rev 2947)
@@ -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@... or */
+/* vigra@... */
+/* */
+/* 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
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|