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);

  } // namespace computerVision

} // namespace dlr


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

#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;

        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;
          }
        }

        size_t lastRow = gradImage.rows() - 1;
        size_t lastColumn = gradImage.columns() - 1;
        size_t columns = gradImage.columns();
        while(seedList.size() != 0) {
          Index2D seed = *seedList.begin();
          seedList.pop_front();
          size_t row = seed.getRow();
          size_t column = seed.getColumn();
          if(row == 0 || column == 0
             || row == lastRow || column == lastColumn) {
            continue;
          }
          index0 = row * columns + column;

          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 = index0 - columns;
          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 = 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 = 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 = index0 + columns;
          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 = 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));
          }
        }
        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)
    {
      // Argument checking.
      if(inputImage.rows() < gaussianSize + 3
         || inputImage.columns() < gaussianSize + 3) {
        DLR_THROW(ValueException, "applyCanny()",
                  "Argument inputImage has insufficient size.");
      }
      if(lowerThreshold > upperThreshold) {
        DLR_THROW(ValueException, "applyCanny()",
                  "Argument lowerThreshold must be less than or equal to "
                  "Arguments upperThreshold.");
      }

      // 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] + gradX[index0] * gradX[index0]);
        }

        // Pick edge thresholds.
        if(upperThreshold <= 0.0) {
          Float64 maxGrad = maximum(gradMagnitude.ravel());
          Float64 minGrad = minimum(gradMagnitude.ravel());
          upperThreshold = minGrad + 0.9 * (maxGrad - minGrad);
        }
        if(lowerThreshold <= 0.0) {
          lowerThreshold = 0.9 * upperThreshold;
        }
        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_ */

---