Visual Servoing Platform  version 3.6.1 under development (2024-02-13)
vpCannyEdgeDetection.cpp
1 /****************************************************************************
2  *
3  * ViSP, open source Visual Servoing Platform software.
4  * Copyright (C) 2005 - 2023 by Inria. All rights reserved.
5  *
6  * This software is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  * See the file LICENSE.txt at the root directory of this source
11  * distribution for additional information about the GNU GPL.
12  *
13  * For using ViSP with software that can not be combined with the GNU
14  * GPL, please contact Inria about acquiring a ViSP Professional
15  * Edition License.
16  *
17  * See https://visp.inria.fr for more information.
18  *
19  * This software was developed at:
20  * Inria Rennes - Bretagne Atlantique
21  * Campus Universitaire de Beaulieu
22  * 35042 Rennes Cedex
23  * France
24  *
25  * If you have questions regarding the use of this file, please contact
26  * Inria at visp@inria.fr
27  *
28  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
29  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
30  *
31 *****************************************************************************/
32 
33 #include <visp3/core/vpCannyEdgeDetection.h>
34 
35 #include <visp3/core/vpImageConvert.h>
36 
37 #if (VISP_CXX_STANDARD == VISP_CXX_STANDARD_98) // Check if cxx98
38 namespace
39 {
40 // Helper to apply the scale to the raw values of the filters
41 template <typename FilterType>
42 static void scaleFilter(vpArray2D<FilterType> &filter, const float &scale)
43 {
44  const unsigned int nbRows = filter.getRows();
45  const unsigned int nbCols = filter.getCols();
46  for (unsigned int r = 0; r < nbRows; ++r) {
47  for (unsigned int c = 0; c < nbCols; ++c) {
48  filter[r][c] = filter[r][c] * scale;
49  }
50  }
51 }
52 };
53 #endif
54 
55 // // Initialization methods
56 
58  : m_filteringAndGradientType(vpImageFilter::CANNY_GBLUR_SOBEL_FILTERING)
59  , m_gaussianKernelSize(3)
60  , m_gaussianStdev(1.f)
61  , m_areGradientAvailable(false)
62  , m_gradientFilterKernelSize(3)
63  , m_lowerThreshold(-1.f)
64  , m_lowerThresholdRatio(0.6f)
65  , m_upperThreshold(-1.f)
66  , m_upperThresholdRatio(0.8f)
67  , mp_mask(nullptr)
68 {
69  initGaussianFilters();
70  initGradientFilters();
71 }
72 
73 vpCannyEdgeDetection::vpCannyEdgeDetection(const int &gaussianKernelSize, const float &gaussianStdev
74  , const unsigned int &sobelAperture, const float &lowerThreshold, const float &upperThreshold
75  , const float &lowerThresholdRatio, const float &upperThresholdRatio
77 )
78  : m_filteringAndGradientType(filteringType)
79  , m_gaussianKernelSize(gaussianKernelSize)
80  , m_gaussianStdev(gaussianStdev)
81  , m_areGradientAvailable(false)
82  , m_gradientFilterKernelSize(sobelAperture)
83  , m_lowerThreshold(lowerThreshold)
84  , m_lowerThresholdRatio(lowerThresholdRatio)
85  , m_upperThreshold(upperThreshold)
86  , m_upperThresholdRatio(upperThresholdRatio)
87  , mp_mask(nullptr)
88 {
89  initGaussianFilters();
90  initGradientFilters();
91 }
92 
93 #ifdef VISP_HAVE_NLOHMANN_JSON
94 vpCannyEdgeDetection::vpCannyEdgeDetection(const std::string &jsonPath)
95 {
96  initFromJSON(jsonPath);
97 }
98 
99 void
100 vpCannyEdgeDetection::initFromJSON(const std::string &jsonPath)
101 {
102  std::ifstream file(jsonPath);
103  if (!file.good()) {
104  std::stringstream ss;
105  ss << "Problem opening file " << jsonPath << ". Make sure it exists and is readable" << std::endl;
106  throw vpException(vpException::ioError, ss.str());
107  }
108  json j;
109  try {
110  j = json::parse(file);
111  }
112  catch (json::parse_error &e) {
113  std::stringstream msg;
114  msg << "Could not parse JSON file : \n";
115  msg << e.what() << std::endl;
116  msg << "Byte position of error: " << e.byte;
117  throw vpException(vpException::ioError, msg.str());
118  }
119  from_json(j, *this);
120  file.close();
121  initGaussianFilters();
122  initGradientFilters();
123 }
124 #endif
125 
126 void
127 vpCannyEdgeDetection::initGaussianFilters()
128 {
129  if ((m_gaussianKernelSize % 2) == 0) {
130  throw(vpException(vpException::badValue, "The Gaussian kernel size should be odd"));
131  }
132  m_fg.resize(1, (m_gaussianKernelSize + 1) / 2);
133  vpImageFilter::getGaussianKernel(m_fg.data, m_gaussianKernelSize, m_gaussianStdev, true);
134 }
135 
136 void
137 vpCannyEdgeDetection::initGradientFilters()
138 {
139  if ((m_gradientFilterKernelSize % 2) != 1) {
140  throw vpException(vpException::badValue, "Gradient filters kernel size should be odd.");
141  }
142  m_gradientFilterX.resize(m_gradientFilterKernelSize, m_gradientFilterKernelSize);
143  m_gradientFilterY.resize(m_gradientFilterKernelSize, m_gradientFilterKernelSize);
144 
145 #if (VISP_CXX_STANDARD >= VISP_CXX_STANDARD_11)
146  auto scaleFilter = [](vpArray2D<float> &filter, const float &scale) {
147  for (unsigned int r = 0; r < filter.getRows(); r++) {
148  for (unsigned int c = 0; c < filter.getCols(); c++) {
149  filter[r][c] = filter[r][c] * scale;
150  }
151  }};
152 #endif
153 
154  float scaleX = 1.f;
155  float scaleY = 1.f;
156 
157  if (m_filteringAndGradientType == vpImageFilter::CANNY_GBLUR_SOBEL_FILTERING) {
158  scaleX = vpImageFilter::getSobelKernelX(m_gradientFilterX.data, (m_gradientFilterKernelSize - 1) / 2);
159  scaleY = vpImageFilter::getSobelKernelY(m_gradientFilterY.data, (m_gradientFilterKernelSize - 1) / 2);
160  }
161  else if (m_filteringAndGradientType == vpImageFilter::CANNY_GBLUR_SCHARR_FILTERING) {
162  // Compute the Scharr filters
163  scaleX = vpImageFilter::getScharrKernelX(m_gradientFilterX.data, (m_gradientFilterKernelSize - 1) / 2);
164  scaleY = vpImageFilter::getScharrKernelY(m_gradientFilterY.data, (m_gradientFilterKernelSize - 1) / 2);
165  }
166  else {
167  std::string errMsg = "[vpCannyEdgeDetection::initGradientFilters] Error: gradient filtering method \"";
168  errMsg += vpImageFilter::vpCannyFilteringAndGradientTypeToString(m_filteringAndGradientType);
169  errMsg += "\" has not been implemented yet\n";
171  }
172 
173  scaleFilter(m_gradientFilterX, scaleX);
174  scaleFilter(m_gradientFilterY, scaleY);
175 }
176 
177 // // Detection methods
178 #ifdef HAVE_OPENCV_CORE
180 vpCannyEdgeDetection::detect(const cv::Mat &cv_I)
181 {
182  vpImage<unsigned char> I_gray;
183  vpImageConvert::convert(cv_I, I_gray);
184  return detect(I_gray);
185 }
186 #endif
187 
190 {
191  vpImage<unsigned char> I_gray;
192  vpImageConvert::convert(I_color, I_gray);
193  return detect(I_gray);
194 }
195 
198 {
199  // // Clearing the previous results
200  m_edgeMap.resize(I.getHeight(), I.getWidth(), 0);
201  m_edgeCandidateAndGradient.clear();
202  m_edgePointsCandidates.clear();
203 
204  // // Step 1 and 2: filter the image and compute the gradient, if not given by the user
205  if (!m_areGradientAvailable) {
206  performFilteringAndGradientComputation(I);
207  }
208  m_areGradientAvailable = false; // Reset for next call
209 
210  // // Step 3: edge thining
211  float upperThreshold = m_upperThreshold;
212  float lowerThreshold = m_lowerThreshold;
213  if (upperThreshold < 0) {
214  upperThreshold = vpImageFilter::computeCannyThreshold(I, lowerThreshold, &m_dIx, &m_dIy, m_gaussianKernelSize,
215  m_gaussianStdev, m_gradientFilterKernelSize, m_lowerThresholdRatio,
216  m_upperThresholdRatio, m_filteringAndGradientType, mp_mask);
217  }
218  else if (m_lowerThreshold < 0) {
219  // Applying Canny recommendation to have the upper threshold 3 times greater than the lower threshold.
220  lowerThreshold = m_upperThreshold / 3.f;
221  }
222  // To ensure that if lowerThreshold = 0, we reject null gradient points
223  lowerThreshold = std::max<float>(lowerThreshold, std::numeric_limits<float>::epsilon());
224  performEdgeThinning(lowerThreshold);
225 
226  // // Step 4: hysteresis thresholding
227  performHysteresisThresholding(lowerThreshold, upperThreshold);
228 
229  // // Step 5: edge tracking
230  performEdgeTracking();
231  return m_edgeMap;
232 }
233 
234 void
235 vpCannyEdgeDetection::performFilteringAndGradientComputation(const vpImage<unsigned char> &I)
236 {
237  if (m_filteringAndGradientType == vpImageFilter::CANNY_GBLUR_SOBEL_FILTERING
238  || m_filteringAndGradientType == vpImageFilter::CANNY_GBLUR_SCHARR_FILTERING) {
239  // Computing the Gaussian blur
240  vpImage<float> Iblur;
241  vpImage<float> GIx;
242  vpImageFilter::filterX<unsigned char, float>(I, GIx, m_fg.data, m_gaussianKernelSize, mp_mask);
243  vpImageFilter::filterY<float, float>(GIx, Iblur, m_fg.data, m_gaussianKernelSize, mp_mask);
244 
245  // Computing the gradients
246  vpImageFilter::filter(Iblur, m_dIx, m_gradientFilterX, true, mp_mask);
247  vpImageFilter::filter(Iblur, m_dIy, m_gradientFilterY, true, mp_mask);
248  }
249  else {
250  std::string errmsg("Currently, the filtering operation \"");
251  errmsg += vpImageFilter::vpCannyFilteringAndGradientTypeToString(m_filteringAndGradientType);
252  errmsg += "\" is not handled.";
254  }
255 }
256 
269 void
270 getInterpolationWeightsAndOffsets(const float &gradientOrientation,
271  float &alpha, float &beta,
272  int &dRowGradAlpha, int &dRowGradBeta,
273  int &dColGradAlpha, int &dColGradBeta
274 )
275 {
276  float thetaMin = 0.f;
277  if (gradientOrientation < M_PI_4f) {
278  // Angles between 0 and 45 deg rely on the horizontal and diagonal points
279  dColGradAlpha = 1;
280  dColGradBeta = 1;
281  dRowGradAlpha = 0;
282  dRowGradBeta = -1;
283  }
284  else if (gradientOrientation >= M_PI_4f && gradientOrientation < M_PI_2f) {
285  // Angles between 45 and 90 deg rely on the diagonal and vertical points
286  thetaMin = M_PI_4f;
287  dColGradAlpha = 1;
288  dColGradBeta = 0;
289  dRowGradAlpha = -1;
290  dRowGradBeta = -1;
291  }
292  else if (gradientOrientation >= M_PI_2f && gradientOrientation < (3.f * M_PI_4f)) {
293  // Angles between 90 and 135 deg rely on the vertical and diagonal points
294  thetaMin = M_PI_2f;
295  dColGradAlpha = 0;
296  dColGradBeta = -1;
297  dRowGradAlpha = -1;
298  dRowGradBeta = -1;
299  }
300  else if (gradientOrientation >= (3.f * M_PI_4f) && gradientOrientation < M_PIf) {
301  // Angles between 135 and 180 deg rely on the vertical and diagonal points
302  thetaMin = 3.f * M_PI_4f;
303  dColGradAlpha = -1;
304  dColGradBeta = -1;
305  dRowGradAlpha = -1;
306  dRowGradBeta = 0;
307  }
308  beta = (gradientOrientation - thetaMin) / M_PI_4f;
309  alpha = 1.f - beta;
310 }
311 
321 float
322 getManhattanGradient(const vpImage<float> &dIx, const vpImage<float> &dIy, const int &row, const int &col)
323 {
324  float grad = 0.;
325  int nbRows = dIx.getRows();
326  int nbCols = dIx.getCols();
327  if (row >= 0
328  && row < nbRows
329  && col >= 0
330  && col < nbCols
331  ) {
332  float dx = dIx[row][col];
333  float dy = dIy[row][col];
334  grad = std::abs(dx) + std::abs(dy);
335  }
336  return grad;
337 }
338 
350 float
351 getGradientOrientation(const vpImage<float> &dIx, const vpImage<float> &dIy, const int &row, const int &col)
352 {
353  float gradientOrientation = 0.f;
354  float dx = dIx[row][col];
355  float dy = dIy[row][col];
356 
357  if (std::abs(dx) < std::numeric_limits<float>::epsilon()) {
358  gradientOrientation = M_PI_2f;
359  }
360  else {
361  // -dy because the y-axis of the image is oriented towards the bottom of the screen
362  // while we later work with a y-axis oriented towards the top when getting the theta quadrant.
363  gradientOrientation = static_cast<float>(std::atan2(-dy, dx));
364  if (gradientOrientation < 0.f) {
365  gradientOrientation += M_PIf; // + M_PI in order to be between 0 and M_PIf
366  }
367  }
368  return gradientOrientation;
369 }
370 
371 void
372 vpCannyEdgeDetection::performEdgeThinning(const float &lowerThreshold)
373 {
374  int nbRows = m_dIx.getRows();
375  int nbCols = m_dIx.getCols();
376 
377  for (int row = 0; row < nbRows; row++) {
378  for (int col = 0; col < nbCols; col++) {
379  if (mp_mask != nullptr) {
380  if (!(*mp_mask)[row][col]) {
381  // The mask tells us to ignore the current pixel
382  continue;
383  }
384  }
385 
386  // Computing the gradient orientation and magnitude
387  float grad = getManhattanGradient(m_dIx, m_dIy, row, col);
388 
389  if (grad < lowerThreshold) {
390  // The gradient is lower than minimum threshold => ignoring the point
391  continue;
392  }
393 
394  // Getting the offset along the horizontal and vertical axes
395  // depending on the gradient orientation
396  int dRowAlphaPlus = 0, dRowBetaPlus = 0;
397  int dColAphaPlus = 0, dColBetaPlus = 0;
398  float gradientOrientation = getGradientOrientation(m_dIx, m_dIy, row, col);
399  float alpha = 0.f, beta = 0.f;
400  getInterpolationWeightsAndOffsets(gradientOrientation, alpha, beta, dRowAlphaPlus, dRowBetaPlus, dColAphaPlus, dColBetaPlus);
401  int dRowAlphaMinus = -dRowAlphaPlus, dRowBetaMinus = -dRowBetaPlus;
402  int dColAphaMinus = -dColAphaPlus, dColBetaMinus = -dColBetaPlus;
403  float gradAlphaPlus = getManhattanGradient(m_dIx, m_dIy, row + dRowAlphaPlus, col + dColAphaPlus);
404  float gradBetaPlus = getManhattanGradient(m_dIx, m_dIy, row + dRowBetaPlus, col + dColBetaPlus);
405  float gradAlphaMinus = getManhattanGradient(m_dIx, m_dIy, row + dRowAlphaMinus, col + dColAphaMinus);
406  float gradBetaMinus = getManhattanGradient(m_dIx, m_dIy, row + dRowBetaMinus, col + dColBetaMinus);
407  float gradPlus = alpha * gradAlphaPlus + beta * gradBetaPlus;
408  float gradMinus = alpha * gradAlphaMinus + beta * gradBetaMinus;
409 
410  if (grad >= gradPlus && grad >= gradMinus) {
411  // Keeping the edge point that has the highest gradient
412  std::pair<unsigned int, unsigned int> bestPixel(row, col);
413  m_edgeCandidateAndGradient[bestPixel] = grad;
414  }
415  }
416  }
417 }
418 
419 void
420 vpCannyEdgeDetection::performHysteresisThresholding(const float &lowerThreshold, const float &upperThreshold)
421 {
422  std::map<std::pair<unsigned int, unsigned int>, float>::iterator it;
423  for (it = m_edgeCandidateAndGradient.begin(); it != m_edgeCandidateAndGradient.end(); it++) {
424  if (it->second >= upperThreshold) {
425  m_edgePointsCandidates[it->first] = STRONG_EDGE;
426  }
427  else if (it->second >= lowerThreshold && it->second < upperThreshold) {
428  m_edgePointsCandidates[it->first] = WEAK_EDGE;
429  }
430  }
431 }
432 
433 void
434 vpCannyEdgeDetection::performEdgeTracking()
435 {
436  std::map<std::pair<unsigned int, unsigned int>, EdgeType>::iterator it;
437  for (it = m_edgePointsCandidates.begin(); it != m_edgePointsCandidates.end(); it++) {
438  if (it->second == STRONG_EDGE) {
439  m_edgeMap[it->first.first][it->first.second] = 255;
440  }
441  else if (it->second == WEAK_EDGE) {
442  if (recursiveSearchForStrongEdge(it->first)) {
443  m_edgeMap[it->first.first][it->first.second] = 255;
444  }
445  }
446  }
447 }
448 
449 bool
450 vpCannyEdgeDetection::recursiveSearchForStrongEdge(const std::pair<unsigned int, unsigned int> &coordinates)
451 {
452  bool hasFoundStrongEdge = false;
453  int nbRows = m_dIx.getRows();
454  int nbCols = m_dIx.getCols();
455  m_edgePointsCandidates[coordinates] = ON_CHECK;
456  for (int dr = -1; dr <= 1 && !hasFoundStrongEdge; dr++) {
457  for (int dc = -1; dc <= 1 && !hasFoundStrongEdge; dc++) {
458  int idRow = dr + (int)coordinates.first;
459  idRow = std::max<int>(idRow, 0); // Avoid getting negative pixel ID
460  int idCol = dc + (int)coordinates.second;
461  idCol = std::max<int>(idCol, 0); // Avoid getting negative pixel ID
462 
463  // Checking if we are still looking for an edge in the limit of the image
464  if ((idRow < 0 || idRow >= nbRows)
465  || (idCol < 0 || idCol >= nbCols)
466  || (dr == 0 && dc == 0)
467  ) {
468  continue;
469  }
470 
471  try {
472  std::pair<unsigned int, unsigned int> key_candidate(idRow, idCol);
473  // Checking if the 8-neighbor point is in the list of edge candidates
474  EdgeType type_candidate = m_edgePointsCandidates.at(key_candidate);
475  if (type_candidate == STRONG_EDGE) {
476  // The 8-neighbor point is a strong edge => the weak edge becomes a strong edge
477  hasFoundStrongEdge = true;
478  }
479  else if (type_candidate == WEAK_EDGE) {
480  // Checking if the WEAK_EDGE neighbor has a STRONG_EDGE neighbor
481  hasFoundStrongEdge = recursiveSearchForStrongEdge(key_candidate);
482  }
483  }
484  catch (...) {
485  continue;
486  }
487  }
488  }
489  if (hasFoundStrongEdge) {
490  m_edgePointsCandidates[coordinates] = STRONG_EDGE;
491  m_edgeMap[coordinates.first][coordinates.second] = 255;
492  }
493  return hasFoundStrongEdge;
494 }
Implementation of a generic 2D array used as base class for matrices and vectors.
Definition: vpArray2D.h:125
unsigned int getCols() const
Definition: vpArray2D.h:274
Type * data
Address of the first element of the data array.
Definition: vpArray2D.h:138
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Definition: vpArray2D.h:299
unsigned int getRows() const
Definition: vpArray2D.h:284
vpImage< unsigned char > detect(const vpImage< vpRGBa > &I_color)
Detect the edges in an image. Convert the color image into a gray-scale image.
friend void from_json(const json &j, vpCannyEdgeDetection &detector)
Read the detector configuration from JSON. All values are optional and if an argument is not present,...
void initFromJSON(const std::string &jsonPath)
Initialize all the algorithm parameters using the JSON file whose path is jsonPath....
vpCannyEdgeDetection()
Default constructor of the vpCannyEdgeDetection class. The thresholds used during the hysteresis thre...
error that can be emitted by ViSP classes.
Definition: vpException.h:59
@ ioError
I/O error.
Definition: vpException.h:79
@ badValue
Used to indicate that a value is not in the allowed range.
Definition: vpException.h:85
@ notImplementedError
Not implemented.
Definition: vpException.h:81
static void convert(const vpImage< unsigned char > &src, vpImage< vpRGBa > &dest)
Various image filter, convolution, etc...
Definition: vpImageFilter.h:70
static FilterType getSobelKernelX(FilterType *filter, unsigned int size)
vpCannyFilteringAndGradientType
Canny filter and gradient operators to apply on the image before the edge detection stage.
@ CANNY_GBLUR_SOBEL_FILTERING
Apply Gaussian blur + Sobel operator on the input image.
@ CANNY_GBLUR_SCHARR_FILTERING
Apply Gaussian blur + Scharr operator on the input image.
static void filter(const vpImage< ImageType > &I, vpImage< FilterType > &If, const vpArray2D< FilterType > &M, bool convolve=false, const vpImage< bool > *p_mask=nullptr)
static std::string vpCannyFilteringAndGradientTypeToString(const vpCannyFilteringAndGradientType &type)
Cast a vpImageFilter::vpCannyFilteringAndGradientType into a string, to know its name.
static float computeCannyThreshold(const cv::Mat &cv_I, const cv::Mat *p_cv_dIx, const cv::Mat *p_cv_dIy, float &lowerThresh, const unsigned int &gaussianKernelSize=5, const float &gaussianStdev=2.f, const unsigned int &apertureGradient=3, const float &lowerThresholdRatio=0.6, const float &upperThresholdRatio=0.8, const vpCannyFilteringAndGradientType &filteringType=CANNY_GBLUR_SOBEL_FILTERING)
Compute the upper Canny edge filter threshold, using Gaussian blur + Sobel or + Scharr operators to c...
static void getGaussianKernel(FilterType *filter, unsigned int size, FilterType sigma=0., bool normalize=true)
static FilterType getScharrKernelY(FilterType *filter, unsigned int size)
static FilterType getSobelKernelY(FilterType *filter, unsigned int size)
static FilterType getScharrKernelX(FilterType *filter, unsigned int size)
unsigned int getWidth() const
Definition: vpImage.h:242
void resize(unsigned int h, unsigned int w)
resize the image : Image initialization
Definition: vpImage.h:780
unsigned int getCols() const
Definition: vpImage.h:175
unsigned int getHeight() const
Definition: vpImage.h:184
unsigned int getRows() const
Definition: vpImage.h:214