canny.h

Go to the documentation of this file.

#ifndef _DLRCOMPUTERVISION_CANNY_H_
#define _DLRCOMPUTERVISION_CANNY_H_

#include <dlrComputerVision/image.h>

namespace dlr {

  namespace computerVision {
    

    template <ImageFormat FORMAT>
    Image<GRAY1>
    applyCanny(const Image<FORMAT>& inputImage,
               size_t gaussianSize = 5,
               double upperThreshold = 0.0,
               double lowerThreshold = 0.0,
               double autoUpperThresholdFactor = 3.0,
               double autoLowerThresholdFactor = 0.0);

  } // namespace computerVision

} // namespace dlr


/* =============== Implementation follows =============== */

#include <limits>
#include <list>
#include <dlrComputerVision/filter.h>
#include <dlrComputerVision/kernels.h>
#include <dlrComputerVision/nonMaximumSuppress.h>
#include <dlrComputerVision/sobel.h>
#include <dlrComputerVision/utilities.h>
#include <dlrNumeric/index2D.h>
#include <dlrNumeric/utilities.h>

namespace dlr {

  namespace computerVision {

    namespace privateCode {

      Image<GRAY1>
      traceEdges(const Image<GRAY_FLOAT64>& gradImage, double threshold)
      {
        Image<GRAY1> edgeImage(gradImage.rows(), gradImage.columns());
        edgeImage = false;

        // Make a pass through the image pushing all the seed points
        // onto our list of edges to trace.
        std::list<Index2D> seedList;
        size_t index0 = 0;
        for(size_t row = 0; row < gradImage.rows(); ++row) {
          for(size_t column = 0; column < gradImage.columns(); ++column) {
            if(gradImage[index0] > threshold) {
              edgeImage[index0] = true;
              seedList.push_front(Index2D(static_cast<int>(row),
                                          static_cast<int>(column)));
            }
            ++index0;
          }
        }

        // Work our way through all edges that must be traced.
        size_t lastRow = gradImage.rows() - 1;
        size_t lastColumn = gradImage.columns() - 1;
        size_t columns = gradImage.columns();
        while(seedList.size() != 0) {

          // Get the next edge to trace.
          Index2D seed = *seedList.begin();
          seedList.pop_front();

          // Don't trace edges on the very outside border of the image.
          size_t row = seed.getRow();
          size_t column = seed.getColumn();
          if(row == 0 || column == 0
             || row == lastRow || column == lastColumn) {
            continue;
          }

          // Use single indexing on the assumption that it's faster.
          index0 = row * columns + column;

          // Inspect each neighbor in turn, marking as appropriate,
          // and pushing any newly marked neighbor pixels onto the
          // list so that they will themselves be traced.
          size_t neighborIndex = index0 - columns - 1;
          if(gradImage(neighborIndex) != 0.0
             && edgeImage(neighborIndex) == false) {
            edgeImage(neighborIndex) = true;
            seedList.push_front(Index2D(static_cast<int>(row) - 1,
                                        static_cast<int>(column) - 1));
          }
          ++neighborIndex;
          if(gradImage(neighborIndex) != 0.0
             && edgeImage(neighborIndex) == false) {
            edgeImage(neighborIndex) = true;
            seedList.push_front(Index2D(static_cast<int>(row) - 1,
                                        static_cast<int>(column)));
          }
          ++neighborIndex;
          if(gradImage(neighborIndex) != 0.0
             && edgeImage(neighborIndex) == false) {
            edgeImage(neighborIndex) = true;
            seedList.push_front(Index2D(static_cast<int>(row) - 1,
                                        static_cast<int>(column) + 1));
          }
          neighborIndex = index0 - 1;
          if(gradImage(neighborIndex) != 0.0
             && edgeImage(neighborIndex) == false) {
            edgeImage(neighborIndex) = true;
            seedList.push_front(Index2D(static_cast<int>(row),
                                        static_cast<int>(column) - 1));
          }
          neighborIndex = index0 + 1;
          if(gradImage(neighborIndex) != 0.0
             && edgeImage(neighborIndex) == false) {
            edgeImage(neighborIndex) = true;
            seedList.push_front(Index2D(static_cast<int>(row),
                                        static_cast<int>(column) + 1));
          }
          neighborIndex = index0 + columns - 1;
          if(gradImage(neighborIndex) != 0.0
             && edgeImage(neighborIndex) == false) {
            edgeImage(neighborIndex) = true;
            seedList.push_front(Index2D(static_cast<int>(row) + 1,
                                        static_cast<int>(column) - 1));
          }
          ++neighborIndex;
          if(gradImage(neighborIndex) != 0.0
             && edgeImage(neighborIndex) == false) {
            edgeImage(neighborIndex) = true;
            seedList.push_front(Index2D(static_cast<int>(row) + 1,
                                        static_cast<int>(column)));
          }
          ++neighborIndex;
          if(gradImage(neighborIndex) != 0.0
             && edgeImage(neighborIndex) == false) {
            edgeImage(neighborIndex) = true;
            seedList.push_front(Index2D(static_cast<int>(row) + 1,
                                        static_cast<int>(column) + 1));
          }
        }
        return edgeImage;
      }
      
    } // namespace privateCode

    
    // This function applies the canny edge operator in the X
    // direction.
    template <ImageFormat FORMAT>
    Image<GRAY1>
    applyCanny(const Image<FORMAT>& inputImage,
               size_t gaussianSize,
               double upperThreshold,
               double lowerThreshold,
               double autoUpperThresholdFactor,
               double autoLowerThresholdFactor)
    {
      // Argument checking.
      if(inputImage.rows() < gaussianSize + 3
         || inputImage.columns() < gaussianSize + 3) {
        DLR_THROW(ValueException, "applyCanny()",
                  "Argument inputImage has insufficient size, or argument "
                  "gaussianSize is too large.");
      }
      if(lowerThreshold > upperThreshold) {
        DLR_THROW(ValueException, "applyCanny()",
                  "Argument lowerThreshold must be less than or equal to "
                  "Arguments upperThreshold.");
      }
      autoLowerThresholdFactor =
        std::min(autoLowerThresholdFactor, autoUpperThresholdFactor);

      // We use non-normalized convolution kernels for the gradient,
      // which means that our user-specified thresholds don't match
      // the magnitude of our gradients.  Each gradient component is
      // 8 times as large as it should be, so the magnitude of the
      // gradient is 8*sqrt(2) times as large.  We solve this by
      // scaling the thresholds here.
      lowerThreshold *= std::sqrt(2.0) * 8.0;
      upperThreshold *= std::sqrt(2.0) * 8.0;
      
      // Step 1: Blur with a gaussian kernel to reduce noise.
      Image<GRAY_FLOAT64> blurredImage;
      if(gaussianSize == 0) {
        blurredImage = convertColorspace<GRAY_FLOAT64>(inputImage);
      } else {
        Kernel<double> gaussian =
          getGaussianKernel<double>(gaussianSize, gaussianSize);
        blurredImage = filter2D<GRAY_FLOAT64, GRAY_FLOAT64>(
          gaussian, inputImage, 0.0);
      }

      // Step 2: Compute derivatives of the blurred image, and discard
      // any which are less than the lower threshold.
      Image<GRAY_FLOAT64> gradX = applySobelX(blurredImage);
      Image<GRAY_FLOAT64> gradY = applySobelY(blurredImage);
      Image<GRAY_FLOAT64> gradMagnitude(gradX.rows(), gradX.columns());

      if(lowerThreshold > 0.0 && upperThreshold > 0.0) {
        // Discard values less than the lower threshold.
        for(size_t index0 = 0; index0 < gradX.size(); ++index0) {
          double tmpVal = std::sqrt(
            gradX[index0] * gradX[index0] + gradY[index0] * gradY[index0]);
          gradMagnitude[index0] = (tmpVal > lowerThreshold) ? tmpVal : 0.0;
        }
      } else {
        // Temporarily retain all gradient values.
        for(size_t index0 = 0; index0 < gradX.size(); ++index0) {
          gradMagnitude[index0] = std::sqrt(
            gradX[index0] * gradX[index0] + gradY[index0] * gradY[index0]);
        }

        // Pick edge thresholds.
        size_t startRow = (gaussianSize + 1) / 2;
        size_t endRow = gradMagnitude.rows() - startRow;
        size_t startColumn = startRow;
        size_t endColumn = gradMagnitude.columns() - startColumn;

        // Paranoid check should never fail.
        if((startRow >= endRow) || (startColumn >= endColumn)) {
          DLR_THROW(ValueException, "applyCanny()",
                    "Filter kernel is too large for image.");
        }
        
        // Hack(xxx): Zero out borders of images to avoid spurious edges.
        for(size_t row = 0; row < startRow; ++row) {
          for(size_t column = 0; column < gradMagnitude.columns(); ++column) {
            gradMagnitude(row, column) = 0.0;
          }
        }
        for(size_t row = endRow; row < gradMagnitude.rows(); ++row) {
          for(size_t column = 0; column < gradMagnitude.columns(); ++column) {
            gradMagnitude(row, column) = 0.0;
          }
        }
        for(size_t row = 0; row < gradMagnitude.rows(); ++row) {
          for(size_t column = 0; column < startColumn; ++column) {
            gradMagnitude(row, column) = 0.0;
          }
          for(size_t column = endColumn; column < gradMagnitude.columns();
              ++column) {
            gradMagnitude(row, column) = 0.0;
          }
        }

        // Compute mean and variance of gradient values.
        size_t numberOfPixels =
          (endRow - startRow) * (endColumn - startColumn);
        Float64 sumOfGradient = 0.0;
        Float64 sumOfGradientSquared = 0.0;
        for(size_t row = startRow; row < endRow; ++row) {
          Float64 subSum = 0.0;
          Float64 subSumOfSquares = 0.0;
          for(size_t column = startColumn; column < endColumn; ++column) {
            Float64 testValue = gradMagnitude(row, column);
            // Changing how we compute threshold...
            //
            // if(testValue < minGrad) {minGrad = testValue;}
            // if(testValue > maxGrad) {maxGrad = testValue;}
            subSum += testValue;
            subSumOfSquares += testValue * testValue;
          }
          sumOfGradient += subSum;
          sumOfGradientSquared += subSumOfSquares;
        }
        double gradientMean = sumOfGradient / numberOfPixels;
        double gradientVariance =
          sumOfGradientSquared / numberOfPixels - gradientMean * gradientMean;
        double gradientSigma = std::sqrt(gradientVariance);

        if(upperThreshold <= 0.0) {
          upperThreshold =
            gradientMean + autoUpperThresholdFactor * gradientSigma;
          upperThreshold = std::max(upperThreshold, 0.0);
        }
        if(lowerThreshold <= 0.0) {
          lowerThreshold = 
            gradientMean + autoLowerThresholdFactor * gradientSigma;
          lowerThreshold = std::min(lowerThreshold, upperThreshold);
          lowerThreshold = std::max(lowerThreshold, 0.0);
        }

        // Now zero out gradients that for sure can never be edges.
        for(size_t index0 = 0; index0 < gradX.size(); ++index0) {
          double tmpVal = gradMagnitude[index0];
          gradMagnitude[index0] = (tmpVal > lowerThreshold) ? tmpVal : 0.0;
        }
      }
          
      // Step 3: Non-maximum suppression.
      Image<GRAY_FLOAT64> edgeCandidates =
        nonMaximumSuppress(gradMagnitude, gradX, gradY);

      // Step 4: Threshold with hysteresis.
      Image<GRAY1> edgeImage =
        privateCode::traceEdges(edgeCandidates, upperThreshold);
      return edgeImage;
    }

  } // namespace computerVision
  
} // namespace dlr

#endif /* #ifndef _DLRCOMPUTERVISION_KERNEL_H_ */