Added EDColor support to EdgeDrawing class.

Fixed a bug where params.PFmode = true caused incorrect edge detection results.
Updated edge_drawing.py to support the new EDColor feature.
Created a new edge_drawing_demo.cpp for demonstrating the EDColor feature.
Updated ximgproc.bib to include references for EDColor.
Added tests for PFmode and EDColor in test_fld.cpp.
Added CV_WRAP to Params::read and Params::write functions to expose them to Python bindings.
the class NFALUT uses the code from original ED code
pull/3764/head
Suleyman TURKMEN 9 months ago
parent 438e67f472
commit b8d6f196d3
  1. 9
      modules/ximgproc/doc/ximgproc.bib
  2. 48
      modules/ximgproc/include/opencv2/ximgproc.hpp
  3. 10
      modules/ximgproc/include/opencv2/ximgproc/edge_drawing.hpp
  4. 149
      modules/ximgproc/samples/edge_drawing.py
  5. 185
      modules/ximgproc/samples/edge_drawing_demo.cpp
  6. 62
      modules/ximgproc/samples/fld_lines.cpp
  7. 673
      modules/ximgproc/src/edge_drawing.cpp
  8. 210
      modules/ximgproc/src/edge_drawing_common.hpp
  9. 46
      modules/ximgproc/test/test_fld.cpp

@ -388,6 +388,15 @@ pages={617--632},
publisher={Elsevier}
}
@article{akinlar201782,
title = {ColorED: Color edge and segment detection by Edge Drawing (ED)},
author = {Cuneyt Akinlar and Cihan Topal},
journal = {Journal of Visual Communication and Image Representation},
volume = {44},
pages = {82-94},
year = {2017},
publisher={Academic Press}
@article{loke2021accelerated,
title={Accelerated superpixel image segmentation with a parallelized DBSCAN algorithm},
author={Loke, Seng Cheong and MacDonald, Bruce A and Parsons, Matthew and W{\"u}nsche, Burkhard Claus},

@ -83,18 +83,42 @@
@defgroup ximgproc_fast_line_detector Fast line detector
@defgroup ximgproc_edge_drawing EdgeDrawing
EDGE DRAWING LIBRARY FOR GEOMETRIC FEATURE EXTRACTION AND VALIDATION
Edge Drawing (ED) algorithm is an proactive approach on edge detection problem. In contrast to many other existing edge detection algorithms which follow a subtractive
approach (i.e. after applying gradient filters onto an image eliminating pixels w.r.t. several rules, e.g. non-maximal suppression and hysteresis in Canny), ED algorithm
works via an additive strategy, i.e. it picks edge pixels one by one, hence the name Edge Drawing. Then we process those random shaped edge segments to extract higher level
edge features, i.e. lines, circles, ellipses, etc. The popular method of extraction edge pixels from the thresholded gradient magnitudes is non-maximal supression that tests
every pixel whether it has the maximum gradient response along its gradient direction and eliminates if it does not. However, this method does not check status of the
neighboring pixels, and therefore might result low quality (in terms of edge continuity, smoothness, thinness, localization) edge segments. Instead of non-maximal supression,
ED points a set of edge pixels and join them by maximizing the total gradient response of edge segments. Therefore it can extract high quality edge segments without need for
an additional hysteresis step.
@defgroup ximgproc_edge_drawing Edge Drawing
Edge Drawing (ED) algorithm for geometric feature extraction and validation.
The Edge Drawing (ED) algorithm is a proactive approach to the edge detection problem.
In contrast to many existing edge detection algorithms, which follow a subtractive
approach (i.e., applying gradient filters and eliminating pixels based on several rules,
such as non-maximal suppression and hysteresis in the Canny Edge Detector), the ED algorithm
operates via an additive strategy. It selects edge pixels one by one and connects them,
hence the name Edge Drawing.
ED offers several key advantages:
1. **Additive Strategy**: Instead of eliminating non-edge pixels after gradient filtering,
ED incrementally builds up edge segments by selecting and connecting pixels based on
their gradient response. This differs from traditional methods, which rely on
non-maximal suppression and hysteresis to filter out non-edge pixels.
2. **Edge Pixel Selection**: ED selects edge pixels by analyzing their local gradient
response, while also considering neighboring pixels. This results in smoother and
more continuous edge segments, as ED aims to maximize the overall gradient strength
along the edge segment.
3. **Edge Segment Formation**: Traditional methods, such as non-maximal suppression,
check whether a pixel has the maximum gradient response along its gradient direction,
eliminating it otherwise. However, this approach doesn't consider neighboring pixels,
often resulting in lower-quality edge segments. ED, on the other hand, joins a set of
edge pixels together by maximizing the total gradient response of the segment, leading
to high-quality, well-localized edges.
4. **Higher-Level Feature Extraction**: After forming edge segments, ED enables the
extraction of higher-level geometric features such as lines, circles, ellipses, and
other shapes, making it useful for tasks involving geometric feature extraction and validation.
The ED algorithm produces continuous, smooth, and localized edge segments, making it ideal
for applications requiring precise edge detection and geometric shape analysis.
@defgroup ximgproc_fourier Fourier descriptors

@ -15,7 +15,7 @@ namespace ximgproc
//! @addtogroup ximgproc_edge_drawing
//! @{
/** @brief Class implementing the ED (EdgeDrawing) @cite topal2012edge, EDLines @cite akinlar2011edlines, EDPF @cite akinlar2012edpf and EDCircles @cite akinlar2013edcircles algorithms
/** @brief Class implementing the ED (EdgeDrawing) @cite topal2012edge, EDLines @cite akinlar2011edlines, EDPF @cite akinlar2012edpf, EDCircles @cite akinlar2013edcircles and ColorED @cite akinlar201782 algorithms.
*/
class CV_EXPORTS_W EdgeDrawing : public Algorithm
@ -66,13 +66,13 @@ public:
//! Default value is 1.3
CV_PROP_RW double MaxErrorThreshold;
void read(const FileNode& fn);
void write(FileStorage& fs) const;
CV_WRAP void read(const FileNode& fn);
CV_WRAP void write(FileStorage& fs) const;
};
/** @brief Detects edges in a grayscale image and prepares them to detect lines and ellipses.
/** @brief Detects edges in a grayscale or color image and prepares them to detect lines and ellipses.
@param src 8-bit, single-channel, grayscale input image.
@param src 8-bit, single-channel (CV_8UC1) or color (CV_8UC3, CV_8UC4) input image.
*/
CV_WRAP virtual void detectEdges(InputArray src) = 0;

@ -1,7 +1,24 @@
#!/usr/bin/python
'''
This example illustrates how to use cv.ximgproc.EdgeDrawing class.
This example script illustrates how to use cv.ximgproc.EdgeDrawing class.
It uses the OpenCV library to load an image, and then use the EdgeDrawing class
to detect edges, lines, and ellipses. The detected features are then drawn and displayed.
The main loop allows the user changing parameters of EdgeDrawing by pressing following keys:
to toggle the grayscale conversion press 'space' key
to increase MinPathLength value press '/' key
to decrease MinPathLength value press '*' key
to increase MinLineLength value press '+' key
to decrease MinLineLength value press '-' key
to toggle NFAValidation value press 'n' key
to toggle PFmode value press 'p' key
to save parameters to file press 's' key
to load parameters from file press 'l' key
The program exits when the Esc key is pressed.
Usage:
ed.py [<image_name>]
@ -16,71 +33,117 @@ import cv2 as cv
import random as rng
import sys
rng.seed(12345)
def main():
try:
fn = sys.argv[1]
except IndexError:
fn = 'board.jpg'
src = cv.imread(cv.samples.findFile(fn))
gray = cv.cvtColor(src, cv.COLOR_BGR2GRAY)
cv.imshow("source", src)
ssrc = src.copy()*0
def EdgeDrawingDemo(src, ed, EDParams, convert_to_gray):
rng.seed(12345)
ssrc = np.zeros_like(src)
lsrc = src.copy()
esrc = src.copy()
ed = cv.ximgproc.createEdgeDrawing()
img_to_detect = cv.cvtColor(src, cv.COLOR_BGR2GRAY) if convert_to_gray else src
# you can change parameters (refer the documentation to see all parameters)
EDParams = cv.ximgproc_EdgeDrawing_Params()
EDParams.MinPathLength = 50 # try changing this value between 5 to 1000
EDParams.PFmode = False # defaut value try to swich it to True
EDParams.MinLineLength = 10 # try changing this value between 5 to 100
EDParams.NFAValidation = True # defaut value try to swich it to False
cv.imshow("source image", img_to_detect)
ed.setParams(EDParams)
print("")
print("convert_to_gray:", convert_to_gray)
print("MinPathLength:", EDParams.MinPathLength)
print("MinLineLength:", EDParams.MinLineLength)
print("PFmode:", EDParams.PFmode)
print("NFAValidation:", EDParams.NFAValidation)
tm = cv.TickMeter()
tm.start()
# Detect edges
# you should call this before detectLines() and detectEllipses()
ed.detectEdges(gray)
ed.detectEdges(img_to_detect)
segments = ed.getSegments()
lines = ed.detectLines()
ellipses = ed.detectEllipses()
#Draw detected edge segments
for i in range(len(segments)):
color = (rng.randint(0,256), rng.randint(0,256), rng.randint(0,256))
cv.polylines(ssrc, [segments[i]], False, color, 1, cv.LINE_8)
tm.stop()
print("Detection time : {:.2f} ms. using the parameters above".format(tm.getTimeMilli()))
# Draw detected edge segments
for segment in segments:
color = (rng.randint(0, 256), rng.randint(0, 256), rng.randint(0, 256))
cv.polylines(ssrc, [segment], False, color, 1, cv.LINE_8)
cv.imshow("detected edge segments", ssrc)
#Draw detected lines
if lines is not None: # Check if the lines have been found and only then iterate over these and add them to the image
# Draw detected lines
if lines is not None: # Check if the lines have been found and only then iterate over these and add them to the image
lines = np.uint16(np.around(lines))
for i in range(len(lines)):
cv.line(lsrc, (lines[i][0][0], lines[i][0][1]), (lines[i][0][2], lines[i][0][3]), (0, 0, 255), 1, cv.LINE_AA)
for line in lines:
cv.line(lsrc, (line[0][0], line[0][1]), (line[0][2], line[0][3]), (0, 0, 255), 1, cv.LINE_AA)
cv.imshow("detected lines", lsrc)
#Draw detected circles and ellipses
if ellipses is not None: # Check if circles and ellipses have been found and only then iterate over these and add them to the image
for i in range(len(ellipses)):
center = (int(ellipses[i][0][0]), int(ellipses[i][0][1]))
axes = (int(ellipses[i][0][2])+int(ellipses[i][0][3]),int(ellipses[i][0][2])+int(ellipses[i][0][4]))
angle = ellipses[i][0][5]
color = (0, 0, 255)
if ellipses[i][0][2] == 0:
color = (0, 255, 0)
cv.ellipse(esrc, center, axes, angle,0, 360, color, 2, cv.LINE_AA)
# Draw detected circles and ellipses
if ellipses is not None: # Check if circles and ellipses have been found and only then iterate over these and add them to the image
for ellipse in ellipses:
center = (int(ellipse[0][0]), int(ellipse[0][1]))
axes = (int(ellipse[0][2] + ellipse[0][3]), int(ellipse[0][2] + ellipse[0][4]))
angle = ellipse[0][5]
color = (0, 255, 0) if ellipse[0][2] == 0 else (0, 0, 255)
cv.ellipse(esrc, center, axes, angle, 0, 360, color, 2, cv.LINE_AA)
cv.imshow("detected circles and ellipses", esrc)
cv.waitKey(0)
print('Done')
def main():
try:
fn = sys.argv[1]
except IndexError:
fn = 'board.jpg'
src = cv.imread(cv.samples.findFile(fn))
if src is None:
print("Error loading image")
return
ed = cv.ximgproc.createEdgeDrawing()
# Set parameters (refer to the documentation for all parameters)
EDParams = cv.ximgproc_EdgeDrawing_Params()
EDParams.MinPathLength = 10 # try changing this value by pressing '/' and '*' keys
EDParams.MinLineLength = 10 # try changing this value by pressing '+' and '-' keys
EDParams.PFmode = False # default value is False, try switching by pressing 'p' key
EDParams.NFAValidation = True # default value is True, try switching by pressing 'n' key
convert_to_gray = True
key = 0
while key != 27:
ed.setParams(EDParams)
EdgeDrawingDemo(src, ed, EDParams, convert_to_gray)
key = cv.waitKey()
if key == 32: # space key
convert_to_gray = not convert_to_gray
if key == 112: # 'p' key
EDParams.PFmode = not EDParams.PFmode
if key == 110: # 'n' key
EDParams.NFAValidation = not EDParams.NFAValidation
if key == 43: # '+' key
EDParams.MinLineLength = EDParams.MinLineLength + 5
if key == 45: # '-' key
EDParams.MinLineLength = max(0, EDParams.MinLineLength - 5)
if key == 47: # '/' key
EDParams.MinPathLength = EDParams.MinPathLength + 20
if key == 42: # '*' key
EDParams.MinPathLength = max(0, EDParams.MinPathLength - 20)
if key == 115: # 's' key
fs = cv.FileStorage("ed-params.xml",cv.FileStorage_WRITE)
EDParams.write(fs)
fs.release()
print("parameters saved to ed-params.xml")
if key == 108: # 'l' key
fs = cv.FileStorage("ed-params.xml",cv.FileStorage_READ)
if fs.isOpened():
EDParams.read(fs.root())
fs.release()
print("parameters loaded from ed-params.xml")
if __name__ == '__main__':
print(__doc__)

@ -0,0 +1,185 @@
/* edge_drawing.cpp
This example illustrates how to use cv.ximgproc.EdgeDrawing class.
It uses the OpenCV library to load an image, and then use the EdgeDrawing class
to detect edges, lines, and ellipses. The detected features are then drawn and displayed.
The main loop allows the user changing parameters of EdgeDrawing by pressing following keys:
to toggle the grayscale conversion press 'space' key
to increase MinPathLength value press '/' key
to decrease MinPathLength value press '*' key
to increase MinLineLength value press '+' key
to decrease MinLineLength value press '-' key
to toggle NFAValidation value press 'n' key
to toggle PFmode value press 'p' key
to save parameters to file press 's' key
to load parameters from file press 'l' key
The program exits when the Esc key is pressed.
*/
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/ximgproc.hpp>
#include <iostream>
void EdgeDrawingDemo(const cv::Mat src, cv::Ptr<cv::ximgproc::EdgeDrawing> ed, bool convert_to_gray);
void EdgeDrawingDemo(const cv::Mat src, cv::Ptr<cv::ximgproc::EdgeDrawing> ed, bool convert_to_gray)
{
cv::Mat ssrc = cv::Mat::zeros(src.size(), src.type());
cv::Mat lsrc = src.clone();
cv::Mat esrc = src.clone();
std::cout << std::endl << "convert_to_gray: " << convert_to_gray << std::endl;
std::cout << "MinPathLength: " << ed->params.MinPathLength << std::endl;
std::cout << "MinLineLength: " << ed->params.MinLineLength << std::endl;
std::cout << "PFmode: " << ed->params.PFmode << std::endl;
std::cout << "NFAValidation: " << ed->params.NFAValidation << std::endl;
cv::TickMeter tm;
tm.start();
cv::Mat img_to_detect;
if (convert_to_gray)
{
cv::cvtColor(src, img_to_detect, cv::COLOR_BGR2GRAY);
}
else
{
img_to_detect = src;
}
cv::imshow("source image", img_to_detect);
tm.start();
// Detect edges
ed->detectEdges(img_to_detect);
std::vector<std::vector<cv::Point>> segments = ed->getSegments();
std::vector<cv::Vec4f> lines;
ed->detectLines(lines);
std::vector<cv::Vec6d> ellipses;
ed->detectEllipses(ellipses);
tm.stop();
cv::RNG& rng = cv::theRNG();
cv::setRNGSeed(0);
// Draw detected edge segments
for (const auto& segment : segments)
{
cv::Scalar color(rng.uniform(0, 256), rng.uniform(0, 256), rng.uniform(0, 256));
cv::polylines(ssrc, segment, false, color, 1, cv::LINE_8);
}
cv::imshow("detected edge segments", ssrc);
// Draw detected lines
if (!lines.empty()) // Check if the lines have been found and only then iterate over these and add them to the image
{
for (size_t i = 0; i < lines.size(); i++)
{
cv::line(lsrc, cv::Point2d(lines[i][0], lines[i][1]), cv::Point2d(lines[i][2], lines[i][3]), cv::Scalar(0, 0, 255), 1, cv::LINE_AA);
}
}
cv::imshow("detected lines", lsrc);
// Draw detected circles and ellipses
if (!ellipses.empty()) // Check if circles and ellipses have been found and only then iterate over these and add them to the image
{
for (const auto& ellipse : ellipses)
{
cv::Point center((int)ellipse[0], (int)ellipse[1]);
cv::Size axes((int)ellipse[2] + (int)ellipse[3], (int)ellipse[2] + (int)ellipse[4]);
double angle(ellipse[5]);
cv::Scalar color = (ellipse[2] == 0) ? cv::Scalar(0, 255, 0) : cv::Scalar(0, 0, 255);
cv::ellipse(esrc, center, axes, angle, 0, 360, color, 1, cv::LINE_AA);
}
}
cv::imshow("detected circles and ellipses", esrc);
std::cout << "Total Detection Time : " << tm.getTimeMilli() << "ms." << std::endl;
}
int main(int argc, char** argv)
{
std::string filename = (argc > 1) ? argv[1] : "board.jpg";
cv::Mat src = cv::imread(cv::samples::findFile(filename));
if (src.empty())
{
std::cerr << "Error: Could not open or find the image!" << std::endl;
return -1;
}
cv::Ptr<cv::ximgproc::EdgeDrawing> ed = cv::ximgproc::createEdgeDrawing();
// Set parameters (refer to the documentation for all parameters)
ed->params.MinPathLength = 10; // try changing this value by pressing '/' and '*' keys
ed->params.MinLineLength = 10; // try changing this value by pressing '+' and '-' keys
ed->params.PFmode = false; // default value is false, try switching by pressing 'p' key
ed->params.NFAValidation = true; // default value is true, try switching by pressing 'n' key
bool convert_to_gray = true;
int key = 0;
while (key != 27)
{
EdgeDrawingDemo(src, ed, convert_to_gray);
key = cv::waitKey(0);
switch (key)
{
case 32: // space key
convert_to_gray = !convert_to_gray;
break;
case 'p': // 'p' key
ed->params.PFmode = !ed->params.PFmode;
break;
case 'n': // 'n' key
ed->params.NFAValidation = !ed->params.NFAValidation;
break;
case '+': // '+' key
ed->params.MinLineLength = std::max(0, ed->params.MinLineLength + 5);
break;
case '-': // '-' key
ed->params.MinLineLength = std::max(0, ed->params.MinLineLength - 5);
break;
case '/': // '/' key
ed->params.MinPathLength += 20;
break;
case '*': // '*' key
ed->params.MinPathLength = std::max(0, ed->params.MinPathLength - 20);
break;
case 's': // 's' key
{
cv::FileStorage fs("ed-params.xml", cv::FileStorage::WRITE);
ed->params.write(fs);
fs.release();
std::cout << "Parameters saved to ed-params.xml" << std::endl;
}
break;
case 'l': // 'l' key
{
cv::FileStorage fs("ed-params.xml", cv::FileStorage::READ);
if (fs.isOpened())
{
ed->params.read(fs.root());
fs.release();
std::cout << "Parameters loaded from ed-params.xml" << std::endl;
}
}
break;
default:
break;
}
}
return 0;
}

@ -24,6 +24,7 @@ int main(int argc, char** argv)
if( image.empty() )
{
parser.printMessage();
return -1;
}
@ -54,7 +55,7 @@ int main(int argc, char** argv)
vector<Vec4f> lines;
// Because of some CPU's power strategy, it seems that the first running of
// an algorithm takes much longer. So here we run the algorithm 10 times
// an algorithm takes much longer. So here we run the algorithm 5 times
// to see the algorithm's processing time with sufficiently warmed-up
// CPU performance.
for (int run_count = 0; run_count < 5; run_count++) {
@ -71,64 +72,7 @@ int main(int argc, char** argv)
Mat line_image_fld(image);
fld->drawSegments(line_image_fld, lines);
imshow("FLD result", line_image_fld);
waitKey(1);
Ptr<EdgeDrawing> ed = createEdgeDrawing();
ed->params.EdgeDetectionOperator = EdgeDrawing::SOBEL;
ed->params.GradientThresholdValue = 38;
ed->params.AnchorThresholdValue = 8;
vector<Vec6d> ellipses;
for (int run_count = 0; run_count < 5; run_count++) {
double freq = getTickFrequency();
lines.clear();
int64 start = getTickCount();
// Detect edges
//you should call this before detectLines() and detectEllipses()
ed->detectEdges(image);
// Detect lines
ed->detectLines(lines);
double duration_ms = double(getTickCount() - start) * 1000 / freq;
cout << "Elapsed time for EdgeDrawing detectLines " << duration_ms << " ms." << endl;
start = getTickCount();
// Detect circles and ellipses
ed->detectEllipses(ellipses);
duration_ms = double(getTickCount() - start) * 1000 / freq;
cout << "Elapsed time for EdgeDrawing detectEllipses " << duration_ms << " ms." << endl;
}
Mat edge_image_ed = Mat::zeros(image.size(), CV_8UC3);
vector<vector<Point> > segments = ed->getSegments();
for (size_t i = 0; i < segments.size(); i++)
{
const Point* pts = &segments[i][0];
int n = (int)segments[i].size();
polylines(edge_image_ed, &pts, &n, 1, false, Scalar((rand() & 255), (rand() & 255), (rand() & 255)), 1);
}
imshow("EdgeDrawing detected edges", edge_image_ed);
Mat line_image_ed(image);
fld->drawSegments(line_image_ed, lines);
// Draw circles and ellipses
for (size_t i = 0; i < ellipses.size(); i++)
{
Point center((int)ellipses[i][0], (int)ellipses[i][1]);
Size axes((int)ellipses[i][2] + (int)ellipses[i][3], (int)ellipses[i][2] + (int)ellipses[i][4]);
double angle(ellipses[i][5]);
Scalar color = ellipses[i][2] == 0 ? Scalar(255, 255, 0) : Scalar(0, 255, 0);
ellipse(line_image_ed, center, axes, angle, 0, 360, color, 2, LINE_AA);
}
imshow("EdgeDrawing result", line_image_ed);
waitKey();
return 0;
}

@ -22,8 +22,6 @@ mutable Mat_<uchar> dirImage;
int gradThresh;
int op;
bool SumFlag;
int* grads;
bool PFmode;
};
void ComputeGradientBody::operator() (const Range& range) const
@ -78,9 +76,6 @@ void ComputeGradientBody::operator() (const Range& range) const
gradRow[x] = (ushort)sum;
if (PFmode)
grads[sum]++;
if (sum >= gradThresh)
{
if (gx >= gy)
@ -92,6 +87,9 @@ void ComputeGradientBody::operator() (const Range& range) const
}
}
// Look up table size for fast color space conversion
#define LUT_SIZE (1024*4096)
class EdgeDrawingImpl : public EdgeDrawing
{
public:
@ -126,15 +124,38 @@ protected:
Mat smoothImage;
uchar *edgeImg; // pointer to edge image data
uchar *smoothImg; // pointer to smoothed image data
int segmentNos;
Mat srcImage;
int op; // edge detection operator
int gradThresh; // gradient threshold
int anchorThresh; // anchor point threshold
double divForTestSegment;
double* dH;
int* grads;
int np;
private:
Mat L_Img;
Mat a_Img;
Mat b_Img;
uchar* smooth_L;
uchar* smooth_a;
uchar* smooth_b;
std::vector<std::vector<cv::Point>> segments;
static double LUT1[LUT_SIZE + 1];
static double LUT2[LUT_SIZE + 1];
static bool LUT_Initialized;
void MyRGB2LabFast();
void ComputeGradientMapByDiZenzo();
void smoothChannel(const Mat& src, Mat& dst, double sigma);
static void fixEdgeSegments(std::vector<std::vector<cv::Point>> map);
static void InitColorEDLib();
void ComputeGradient();
void ComputeAnchorPoints();
void JoinAnchorPointsUsingSortedAnchors();
@ -153,10 +174,7 @@ private:
uchar *dirImg; // pointer to direction image data
ushort *gradImg; // pointer to gradient image data
int op; // edge detection operator
int gradThresh; // gradient threshold
int anchorThresh; // anchor point threshold
int segmentNos;
std::vector<EDLineSegment> lines;
int linesNo;
int min_line_len;
@ -337,9 +355,9 @@ void EdgeDrawingImpl::write(cv::FileStorage& fs) const
EdgeDrawingImpl::EdgeDrawingImpl()
{
params = EdgeDrawing::Params();
nfa = new NFALUT(1, 1/2, 1, 1);
dH = new double[MAX_GRAD_VALUE];
grads = new int[MAX_GRAD_VALUE];
nfa = new NFALUT(1, 1/8, 0.5);
dH = new double[ED_MAX_GRAD_VALUE];
grads = new int[ED_MAX_GRAD_VALUE];
}
EdgeDrawingImpl::~EdgeDrawingImpl()
@ -351,7 +369,7 @@ EdgeDrawingImpl::~EdgeDrawingImpl()
void EdgeDrawingImpl::detectEdges(InputArray src)
{
CV_Assert(!src.empty() && src.type() == CV_8UC1);
CV_Assert(!src.empty() && (src.type() == CV_8UC1 || src.type() == CV_8UC3 || src.type() == CV_8UC4));
gradThresh = params.GradientThresholdValue;
anchorThresh = params.AnchorThresholdValue;
op = params.EdgeDetectionOperator;
@ -366,6 +384,9 @@ void EdgeDrawingImpl::detectEdges(InputArray src)
if (anchorThresh < 0)
anchorThresh = 0;
if (params.MinPathLength < 3)
params.MinPathLength = 3;
segmentNos = 0;
anchorNos = 0;
anchorPoints.clear();
@ -377,58 +398,157 @@ void EdgeDrawingImpl::detectEdges(InputArray src)
height = srcImage.rows;
width = srcImage.cols;
edgeImage = Mat(height, width, CV_8UC1, Scalar(0)); // initialize edge Image
gradImage = Mat(height, width, CV_16UC1); // gradImage contains short values
edgeImage = Mat::zeros(height, width, CV_8UC1);
gradImage = Mat::zeros(height, width, CV_16UC1); // gradImage contains short values
dirImage = Mat(height, width, CV_8UC1);
if (params.Sigma < 1.0)
smoothImage = srcImage;
else if (params.Sigma == 1.0)
GaussianBlur(srcImage, smoothImage, Size(5, 5), params.Sigma);
else
GaussianBlur(srcImage, smoothImage, Size(), params.Sigma); // calculate kernel from sigma
// Assign Pointers from Mat's data
smoothImg = smoothImage.data;
gradImg = (ushort*)gradImage.data;
edgeImg = edgeImage.data;
gradImg = (ushort*)gradImage.data;
dirImg = dirImage.data;
if (params.PFmode)
if (srcImage.type() == CV_8UC1)
{
memset(dH, 0, sizeof(double) * MAX_GRAD_VALUE);
memset(grads, 0, sizeof(int) * MAX_GRAD_VALUE);
}
if (params.Sigma < 1.0)
smoothImage = srcImage;
else if (params.Sigma == 1.0)
GaussianBlur(srcImage, smoothImage, Size(5, 5), params.Sigma);
else
GaussianBlur(srcImage, smoothImage, Size(), params.Sigma); // calculate kernel from sigma
smoothImg = smoothImage.data;
ComputeGradient(); // COMPUTE GRADIENT & EDGE DIRECTION MAPS
ComputeAnchorPoints(); // COMPUTE ANCHORS
JoinAnchorPointsUsingSortedAnchors(); // JOIN ANCHORS
if (params.PFmode)
{
memset(gradImg, 0, sizeof(short) * width * height);
memset(grads, 0, sizeof(int) * ED_MAX_GRAD_VALUE);
memset(dH, 0, sizeof(double) * ED_MAX_GRAD_VALUE);
memset(edgeImg, 0, width * height); // clear edge image
ComputeGradient(); // COMPUTE GRADIENT & EDGE DIRECTION MAPS
ComputeAnchorPoints(); // COMPUTE ANCHORS
JoinAnchorPointsUsingSortedAnchors(); // JOIN ANCHORS
GaussianBlur(srcImage, smoothImage, Size(), 0.4);
if (params.PFmode)
for (int i = 1; i < height - 1; i++) {
for (int j = 1; j < width - 1; j++) {
int com1 = smoothImg[(i + 1) * width + j + 1] - smoothImg[(i - 1) * width + j - 1];
int com2 = smoothImg[(i - 1) * width + j + 1] - smoothImg[(i + 1) * width + j - 1];
int gx = abs(com1 + com2 + (smoothImg[i * width + j + 1] - smoothImg[i * width + j - 1]));
int gy = abs(com1 - com2 + (smoothImg[(i + 1) * width + j] - smoothImg[(i - 1) * width + j]));
int g = gx + gy;
gradImg[i * width + j] = (ushort)g;
grads[g]++;
}
}
// Compute probability function H
int size = (width - 2) * (height - 2);
for (int i = ED_MAX_GRAD_VALUE - 1; i > 0; i--)
grads[i - 1] += grads[i];
for (int i = 0; i < ED_MAX_GRAD_VALUE; i++)
dH[i] = (double)grads[i] / ((double)size);
divForTestSegment = 2.25; // Some magic number :-)
np = 0;
for (int i = 0; i < segmentNos; i++)
{
int len = (int)segmentPoints[i].size();
np += (len * (len - 1)) / 2;
}
// Validate segments
for (int i = 0; i < segmentNos; i++)
TestSegment(i, 0, (int)segmentPoints[i].size() - 1);
ExtractNewSegments();
}
}
else
{
// Compute probability function H
int size = (width - 2) * (height - 2);
// Convert RGB2Lab
MyRGB2LabFast();
// Smooth Channels
smoothChannel(L_Img, L_Img, params.Sigma);
smoothChannel(a_Img, a_Img, params.Sigma);
smoothChannel(b_Img, b_Img, params.Sigma);
smooth_L = L_Img.data;
smooth_a = a_Img.data;
smooth_b = b_Img.data;
ComputeGradientMapByDiZenzo(); // Compute Gradient & Edge Direction Maps
// Compute anchors with the user supplied parameters
anchorThresh = 0; // anchorThresh used as zero while computing anchor points if selectStableAnchors set.
// Finding higher number of anchors is OK, because we have following validation steps in selectStableAnchors.
ComputeAnchorPoints(); // COMPUTE ANCHORS
anchorThresh = params.AnchorThresholdValue; // set it to its initial argument value for further anchor validation.
anchorPoints.clear(); // considering validation step below, it should constructed again.
for (int i = 1; i < height - 1; i++) {
for (int j = 1; j < width - 1; j++) {
if (edgeImg[i * width + j] != ANCHOR_PIXEL) continue;
// Take only "stable" anchors
// 0 degree edge
if (edgeImg[i * width + j - 1] && edgeImg[i * width + j + 1]) {
int diff1 = gradImg[i * width + j] - gradImg[(i - 1) * width + j];
int diff2 = gradImg[i * width + j] - gradImg[(i + 1) * width + j];
if (diff1 >= anchorThresh && diff2 >= anchorThresh) edgeImg[i * width + j] = 255;
for (int i = MAX_GRAD_VALUE - 1; i > 0; i--)
grads[i - 1] += grads[i];
continue;
}
for (int i = 0; i < MAX_GRAD_VALUE; i++)
dH[i] = (double)grads[i] / ((double)size);
// 90 degree edge
if (edgeImg[(i - 1) * width + j] && edgeImg[(i + 1) * width + j]) {
int diff1 = gradImg[i * width + j] - gradImg[i * width + j - 1];
int diff2 = gradImg[i * width + j] - gradImg[i * width + j + 1];
if (diff1 >= anchorThresh && diff2 >= anchorThresh) edgeImg[i * width + j] = 255;
divForTestSegment = 2.25; // Some magic number :-)
memset(edgeImg, 0, width * height); // clear edge image
np = 0;
for (int i = 0; i < segmentNos; i++)
{
int len = (int)segmentPoints[i].size();
np += (len * (len - 1)) / 2;
continue;
}
// 135 degree diagonal
if (edgeImg[(i - 1) * width + j - 1] && edgeImg[(i + 1) * width + j + 1]) {
int diff1 = gradImg[i * width + j] - gradImg[(i - 1) * width + j + 1];
int diff2 = gradImg[i * width + j] - gradImg[(i + 1) * width + j - 1];
if (diff1 >= anchorThresh && diff2 >= anchorThresh) edgeImg[i * width + j] = 255;
continue;
}
// 45 degree diagonal
if (edgeImg[(i - 1) * width + j + 1] && edgeImg[(i + 1) * width + j - 1]) {
int diff1 = gradImg[i * width + j] - gradImg[(i - 1) * width + j - 1];
int diff2 = gradImg[i * width + j] - gradImg[(i + 1) * width + j + 1];
if (diff1 >= anchorThresh && diff2 >= anchorThresh) edgeImg[i * width + j] = 255;
}
}
}
// Validate segments
for (int i = 0; i < segmentNos; i++)
TestSegment(i, 0, (int)segmentPoints[i].size() - 1);
for (int i = 0; i < width * height; i++)
if (edgeImg[i] == ANCHOR_PIXEL)
edgeImg[i] = 0;
else if (edgeImg[i] == 255) {
edgeImg[i] = ANCHOR_PIXEL;
int y = i / width;
int x = i % width;
anchorPoints.push_back(Point(x, y)); // push validated anchor point to vector
}
anchorNos = (int)anchorPoints.size(); // get # of anchor pixels
JoinAnchorPointsUsingSortedAnchors(); // JOIN ANCHORS
ExtractNewSegments();
// Fix 1 pixel errors in the edge map
fixEdgeSegments(segments);
}
}
@ -473,8 +593,6 @@ void EdgeDrawingImpl::ComputeGradient()
body.gradThresh = gradThresh;
body.SumFlag = params.SumFlag;
body.op = op;
body.grads = grads;
body.PFmode = params.PFmode;
parallel_for_(Range(1, smoothImage.rows - 1), body);
}
@ -562,24 +680,24 @@ void EdgeDrawingImpl::JoinAnchorPointsUsingSortedAnchors()
{
stack[++top].r = i;
stack[top].c = j;
stack[top].dir = DOWN;
stack[top].dir = ED_DOWN;
stack[top].parent = 0;
stack[++top].r = i;
stack[top].c = j;
stack[top].dir = UP;
stack[top].dir = ED_UP;
stack[top].parent = 0;
}
else
{
stack[++top].r = i;
stack[top].c = j;
stack[top].dir = RIGHT;
stack[top].dir = ED_RIGHT;
stack[top].parent = 0;
stack[++top].r = i;
stack[top].c = j;
stack[top].dir = LEFT;
stack[top].dir = ED_LEFT;
stack[top].parent = 0;
}
@ -609,7 +727,7 @@ StartOfWhile:
len++;
chainLen++;
if (dir == LEFT)
if (dir == ED_LEFT)
{
while (dirImg[r * width + c] == EDGE_HORIZONTAL)
{
@ -680,12 +798,12 @@ StartOfWhile:
stack[++top].r = r;
stack[top].c = c;
stack[top].dir = DOWN;
stack[top].dir = ED_DOWN;
stack[top].parent = noChains;
stack[++top].r = r;
stack[top].c = c;
stack[top].dir = UP;
stack[top].dir = ED_UP;
stack[top].parent = noChains;
len--;
@ -695,7 +813,7 @@ StartOfWhile:
chains[parent].children[0] = noChains;
noChains++;
}
else if (dir == RIGHT)
else if (dir == ED_RIGHT)
{
while (dirImg[r * width + c] == EDGE_HORIZONTAL)
{
@ -766,12 +884,12 @@ StartOfWhile:
stack[++top].r = r;
stack[top].c = c;
stack[top].dir = DOWN; // Go down
stack[top].dir = ED_DOWN; // Go down
stack[top].parent = noChains;
stack[++top].r = r;
stack[top].c = c;
stack[top].dir = UP; // Go up
stack[top].dir = ED_UP; // Go up
stack[top].parent = noChains;
len--;
@ -782,7 +900,7 @@ StartOfWhile:
noChains++;
}
else if (dir == UP)
else if (dir == ED_UP)
{
while (dirImg[r * width + c] == EDGE_VERTICAL)
{
@ -853,12 +971,12 @@ StartOfWhile:
stack[++top].r = r;
stack[top].c = c;
stack[top].dir = RIGHT;
stack[top].dir = ED_RIGHT;
stack[top].parent = noChains;
stack[++top].r = r;
stack[top].c = c;
stack[top].dir = LEFT;
stack[top].dir = ED_LEFT;
stack[top].parent = noChains;
len--;
@ -939,12 +1057,12 @@ StartOfWhile:
stack[++top].r = r;
stack[top].c = c;
stack[top].dir = RIGHT;
stack[top].dir = ED_RIGHT;
stack[top].parent = noChains;
stack[++top].r = r;
stack[top].c = c;
stack[top].dir = LEFT;
stack[top].dir = ED_LEFT;
stack[top].parent = noChains;
len--;
@ -1187,17 +1305,13 @@ StartOfWhile:
delete[] pixels;
}
int* EdgeDrawingImpl::sortAnchorsByGradValue1()
{
int SIZE = 128 * 256;
int* C = new int[SIZE];
memset(C, 0, sizeof(int) * SIZE);
int* EdgeDrawingImpl::sortAnchorsByGradValue1() {
const int SIZE = 128 * 256;
std::vector<int> C(SIZE, 0);
// Count the number of grad values
for (int i = 1; i < height - 1; i++)
{
for (int j = 1; j < width - 1; j++)
{
for (int i = 1; i < height - 1; i++) {
for (int j = 1; j < width - 1; j++) {
if (edgeImg[i * width + j] != ANCHOR_PIXEL)
continue;
@ -1207,16 +1321,15 @@ int* EdgeDrawingImpl::sortAnchorsByGradValue1()
}
// Compute indices
for (int i = 1; i < SIZE; i++)
for (int i = 1; i < SIZE; i++) {
C[i] += C[i - 1];
}
int noAnchors = C[SIZE - 1];
int* A = new int[noAnchors];
for (int i = 1; i < height - 1; i++)
{
for (int j = 1; j < width - 1; j++)
{
for (int i = 1; i < height - 1; i++) {
for (int j = 1; j < width - 1; j++) {
if (edgeImg[i * width + j] != ANCHOR_PIXEL)
continue;
@ -1226,52 +1339,35 @@ int* EdgeDrawingImpl::sortAnchorsByGradValue1()
}
}
delete[] C;
return A;
}
int EdgeDrawingImpl::LongestChain(Chain* chains, int root)
{
if (root == -1 || chains[root].len == 0)
int EdgeDrawingImpl::LongestChain(Chain* chains, int root) {
if (root == -1 || chains[root].len == 0) {
return 0;
}
int len0 = 0;
if (chains[root].children[0] != -1)
len0 = LongestChain(chains, chains[root].children[0]);
int len1 = 0;
if (chains[root].children[1] != -1)
len1 = LongestChain(chains, chains[root].children[1]);
int max = 0;
int leftLength = LongestChain(chains, chains[root].children[0]);
int rightLength = LongestChain(chains, chains[root].children[1]);
if (len0 >= len1)
{
max = len0;
chains[root].children[1] = -1;
if (leftLength >= rightLength) {
chains[root].children[1] = -1; // Invalidate the right child
return chains[root].len + leftLength;
}
else
{
max = len1;
chains[root].children[0] = -1;
else {
chains[root].children[0] = -1; // Invalidate the left child
return chains[root].len + rightLength;
}
return chains[root].len + max;
}
int EdgeDrawingImpl::RetrieveChainNos(Chain* chains, int root, int chainNos[])
{
int EdgeDrawingImpl::RetrieveChainNos(Chain* chains, int root, int chainNos[]) {
int count = 0;
while (root != -1)
{
chainNos[count] = root;
count++;
while (root != -1) {
chainNos[count++] = root;
if (chains[root].children[0] != -1)
root = chains[root].children[0];
else
root = chains[root].children[1];
// Move to the next child in the chain
root = (chains[root].children[0] != -1) ? chains[root].children[0] : chains[root].children[1];
}
return count;
@ -1291,7 +1387,7 @@ void EdgeDrawingImpl::detectLines(OutputArray _lines)
max_distance_between_two_lines = params.MaxDistanceBetweenTwoLines;
max_error = params.MaxErrorThreshold;
if (min_line_len == -1) // If no initial value given, compute it
if (min_line_len < 0) // If no initial value given, compute it
min_line_len = ComputeMinLineLength();
if (min_line_len < 9) // avoids small line segments in the result. Might be deleted!
@ -1318,7 +1414,7 @@ void EdgeDrawingImpl::detectLines(OutputArray _lines)
JoinCollinearLines();
if (params.NFAValidation)
if (params.NFAValidation && srcImage.channels() < 3)
ValidateLineSegments();
// Delete redundant space from lines
@ -1440,7 +1536,7 @@ void EdgeDrawingImpl::SplitSegment2Lines(double* x, double* y, int noPixels, int
index--;
ComputeClosestPoint(x[index], y[index], lastA, lastB, lastInvert, ex, ey);
if ((sx == ex) & (sy == ey))
if ((sx == ex) && (sy == ey))
break;
// Add the line segment to lines
@ -1517,7 +1613,8 @@ void EdgeDrawingImpl::ValidateLineSegments()
{
int lutSize = (width + height) / 8;
double prob = 1.0 / 8; // probability of alignment
nfa = new NFALUT(lutSize, prob, width, height);
double logNT = 2.0 * (log10((double)width) + log10((double)height));
nfa = new NFALUT(lutSize, prob, logNT);
}
int* x = new int[(width + height) * 4];
@ -2037,7 +2134,7 @@ double EdgeDrawingImpl::ComputeMinDistanceBetweenTwoLines(EDLineSegment* ls1, ED
double dy = ls1->sy - ls2->sy;
double d = sqrt(dx * dx + dy * dy);
double min = d;
int which = SOUTH_SOUTH;
int which = ED_SOUTH_SOUTH;
dx = ls1->sx - ls2->ex;
dy = ls1->sy - ls2->ey;
@ -2045,7 +2142,7 @@ double EdgeDrawingImpl::ComputeMinDistanceBetweenTwoLines(EDLineSegment* ls1, ED
if (d < min)
{
min = d;
which = SOUTH_EAST;
which = ED_SOUTH_EAST;
}
dx = ls1->ex - ls2->sx;
@ -2054,7 +2151,7 @@ double EdgeDrawingImpl::ComputeMinDistanceBetweenTwoLines(EDLineSegment* ls1, ED
if (d < min)
{
min = d;
which = EAST_SOUTH;
which = ED_EAST_SOUTH;
}
dx = ls1->ex - ls2->ex;
@ -2063,7 +2160,7 @@ double EdgeDrawingImpl::ComputeMinDistanceBetweenTwoLines(EDLineSegment* ls1, ED
if (d < min)
{
min = d;
which = EAST_EAST;
which = ED_EAST_EAST;
}
if (pwhich)
@ -2370,51 +2467,53 @@ void EdgeDrawingImpl::TestSegment(int i, int index1, int index2)
TestSegment(i, start, index2);
}
//----------------------------------------------------------------------------------------------
// After the validation of the edge segments, extracts the valid ones
// In other words, updates the valid segments' pixel arrays and their lengths
// After validating the edge segments, this function extracts the valid ones.
// Specifically, it updates the valid segments' pixel arrays and their lengths based on validation.
void EdgeDrawingImpl::ExtractNewSegments()
{
vector< vector<Point> > validSegments;
int noSegments = 0;
for (int i = 0; i < segmentNos; i++) {
int start = 0;
while (start < (int)segmentPoints[i].size()) {
// Find the first valid start point
while (start < (int)segmentPoints[i].size()) {
int r = segmentPoints[i][start].y;
int c = segmentPoints[i][start].x;
if (edgeImg[r * width + c]) break;
if (edgeImg[r * width + c]) break; // Found valid point
start++;
}
int end = start + 1;
// Find the end of the valid segment
while (end < (int)segmentPoints[i].size()) {
int r = segmentPoints[i][end].y;
int c = segmentPoints[i][end].x;
if (edgeImg[r * width + c] == 0) break;
if (edgeImg[r * width + c] == 0) break; // End of segment
end++;
}
// Length of the segment
int len = end - start;
// Only accept segments of length >= 10
if (len >= 10) {
// A new segment. Accepted only only long enough (whatever that means)
//segments[noSegments].pixels = &map->segments[i].pixels[start];
//segments[noSegments].noPixels = len;
validSegments.push_back(vector<Point>());
vector<Point> subVec(&segmentPoints[i][start], &segmentPoints[i][end - 1]);
validSegments[noSegments] = subVec;
noSegments++;
vector<Point> subVec(segmentPoints[i].begin() + start, segmentPoints[i].begin() + end); // Safer bounds
validSegments.push_back(subVec); // Add to valid segments
}
// Move start to next unprocessed point
start = end + 1;
}
}
// Replace the old segments with valid segments
segmentPoints = validSegments;
segmentNos = noSegments;
segmentNos = (int)validSegments.size(); // Update number of segments
}
double EdgeDrawingImpl::NFA(double prob, int len)
@ -3119,7 +3218,8 @@ void EdgeDrawingImpl::ValidateCircles(bool validate)
{
int lutSize = (width + height) / 8;
double prob = 1.0 / 8; // probability of alignment
nfa = new NFALUT(lutSize, prob, width, height); // create look up table
double logNT = 2.0 * (log10((double)width) + log10((double)height));
nfa = new NFALUT(lutSize, prob, logNT); // create look up table
}
// Validate circles & ellipses
@ -3297,12 +3397,24 @@ void EdgeDrawingImpl::ValidateCircles(bool validate)
// This produces less false positives, but occationally misses on some valid circles
}
out:
// compute gx & gy
int com1 = smoothImg[(r + 1) * width + c + 1] - smoothImg[(r - 1) * width + c - 1];
int com2 = smoothImg[(r - 1) * width + c + 1] - smoothImg[(r + 1) * width + c - 1];
int com1, com2, gx, gy;
if (srcImage.channels() > 1)
{
com1 = smooth_L[(r + 1) * width + c + 1] - smooth_L[(r - 1) * width + c - 1];
com2 = smooth_L[(r - 1) * width + c + 1] - smooth_L[(r + 1) * width + c - 1];
gx = com1 + com2 + smooth_L[r * width + c + 1] - smooth_L[r * width + c - 1];
gy = com1 - com2 + smooth_L[(r + 1) * width + c] - smooth_L[(r - 1) * width + c];
}
else
{
com1 = smoothImg[(r + 1) * width + c + 1] - smoothImg[(r - 1) * width + c - 1];
com2 = smoothImg[(r - 1) * width + c + 1] - smoothImg[(r + 1) * width + c - 1];
gx = com1 + com2 + smoothImg[r * width + c + 1] - smoothImg[r * width + c - 1];
gy = com1 - com2 + smoothImg[(r + 1) * width + c] - smoothImg[(r - 1) * width + c];
}
int gx = com1 + com2 + smoothImg[r * width + c + 1] - smoothImg[r * width + c - 1];
int gy = com1 - com2 + smoothImg[(r + 1) * width + c] - smoothImg[(r - 1) * width + c];
double pixelAngle = nfa->myAtan2((double)gx, (double)-gy);
double derivX, derivY;
@ -3778,7 +3890,7 @@ void EdgeDrawingImpl::JoinArcs1()
EX = arcs[CandidateArcNo].sx;
EY = arcs[CandidateArcNo].sy;
break;
} //end-switch
}
break; // Do not look at the other candidates
}
@ -5745,5 +5857,276 @@ void EdgeDrawingImpl::ROTATE(double** a, int i, int j, int k, int l, double tau,
a[k][l] = h + s * (g - h * tau);
}
void EdgeDrawingImpl::MyRGB2LabFast()
{
const uchar* blueImg;
const uchar* greenImg;
const uchar* redImg;
// Split channels (OpenCV uses BGR)
vector<Mat> bgr(3);
split(srcImage, bgr);
blueImg = bgr[0].data;
greenImg = bgr[1].data;
redImg = bgr[2].data;
// Initialize LUTs if necessary
if (!LUT_Initialized)
InitColorEDLib();
L_Img.create(height, width, CV_8U);
a_Img.create(height, width, CV_8U);
b_Img.create(height, width, CV_8U);
std::vector<double> L(width * height);
std::vector<double> a(width * height);
std::vector<double> b(width * height);
// Conversion from RGB to Lab
for (int i = 0; i < width * height; i++)
{
double red = LUT1[static_cast<int>(redImg[i] / 255.0 * LUT_SIZE + 0.5)] * 100;
double green = LUT1[static_cast<int>(greenImg[i] / 255.0 * LUT_SIZE + 0.5)] * 100;
double blue = LUT1[static_cast<int>(blueImg[i] / 255.0 * LUT_SIZE + 0.5)] * 100;
double x = red * 0.4124564 + green * 0.3575761 + blue * 0.1804375;
double y = red * 0.2126729 + green * 0.7151522 + blue * 0.0721750;
double z = red * 0.0193339 + green * 0.1191920 + blue * 0.9503041;
x /= 95.047;
y /= 100.0;
z /= 108.883;
x = LUT2[static_cast<int>(x * LUT_SIZE + 0.5)];
y = LUT2[static_cast<int>(y * LUT_SIZE + 0.5)];
z = LUT2[static_cast<int>(z * LUT_SIZE + 0.5)];
L[i] = (116.0 * y) - 16;
a[i] = 500.0 * (x / y);
b[i] = 200.0 * (y - z);
}
// Scaling to [0, 255]
auto scale_to_uchar = [](std::vector<double>& src, uchar* dst, int size) {
double minVal = *std::min_element(src.begin(), src.end());
double maxVal = *std::max_element(src.begin(), src.end());
double scale = 255.0 / (maxVal - minVal);
for (int i = 0; i < size; i++) {
dst[i] = static_cast<uchar>((src[i] - minVal) * scale);
}
};
scale_to_uchar(L, L_Img.data, width * height);
scale_to_uchar(a, a_Img.data, width * height);
scale_to_uchar(b, b_Img.data, width * height);
}
void EdgeDrawingImpl::ComputeGradientMapByDiZenzo()
{
int max = 0;
for (int i = 1; i < height - 1; i++) {
for (int j = 1; j < width - 1; j++) {
int com1, com2, gxCh1, gxCh2, gxCh3, gyCh1, gyCh2, gyCh3;
if(params.EdgeDetectionOperator == PREWITT)
{
// Prewitt for channel1
com1 = smooth_L[(i + 1) * width + j + 1] - smooth_L[(i - 1) * width + j - 1];
com2 = smooth_L[(i - 1) * width + j + 1] - smooth_L[(i + 1) * width + j - 1];
gxCh1 = com1 + com2 + (smooth_L[i * width + j + 1] - smooth_L[i * width + j - 1]);
gyCh1 = com1 - com2 + (smooth_L[(i + 1) * width + j] - smooth_L[(i - 1) * width + j]);
// Prewitt for channel2
com1 = smooth_a[(i + 1) * width + j + 1] - smooth_a[(i - 1) * width + j - 1];
com2 = smooth_a[(i - 1) * width + j + 1] - smooth_a[(i + 1) * width + j - 1];
gxCh2 = com1 + com2 + (smooth_a[i * width + j + 1] - smooth_a[i * width + j - 1]);
gyCh2 = com1 - com2 + (smooth_a[(i + 1) * width + j] - smooth_a[(i - 1) * width + j]);
// Prewitt for channel3
com1 = smooth_b[(i + 1) * width + j + 1] - smooth_b[(i - 1) * width + j - 1];
com2 = smooth_b[(i - 1) * width + j + 1] - smooth_b[(i + 1) * width + j - 1];
gxCh3 = com1 + com2 + (smooth_b[i * width + j + 1] - smooth_b[i * width + j - 1]);
gyCh3 = com1 - com2 + (smooth_b[(i + 1) * width + j] - smooth_b[(i - 1) * width + j]);
}
else
{
// Sobel for channel1
com1 = smooth_L[(i + 1) * width + j + 1] - smooth_L[(i - 1) * width + j - 1];
com2 = smooth_L[(i - 1) * width + j + 1] - smooth_L[(i + 1) * width + j - 1];
gxCh1 = com1 + com2 + 2 * (smooth_L[i * width + j + 1] - smooth_L[i * width + j - 1]);
gyCh1 = com1 - com2 + 2 * (smooth_L[(i + 1) * width + j] - smooth_L[(i - 1) * width + j]);
// Sobel for channel2
com1 = smooth_a[(i + 1) * width + j + 1] - smooth_a[(i - 1) * width + j - 1];
com2 = smooth_a[(i - 1) * width + j + 1] - smooth_a[(i + 1) * width + j - 1];
gxCh2 = com1 + com2 + 2 * (smooth_a[i * width + j + 1] - smooth_a[i * width + j - 1]);
gyCh2 = com1 - com2 + 2 * (smooth_a[(i + 1) * width + j] - smooth_a[(i - 1) * width + j]);
// Sobel for channel3
com1 = smooth_b[(i + 1) * width + j + 1] - smooth_b[(i - 1) * width + j - 1];
com2 = smooth_b[(i - 1) * width + j + 1] - smooth_b[(i + 1) * width + j - 1];
gxCh3 = com1 + com2 + 2 * (smooth_b[i * width + j + 1] - smooth_b[i * width + j - 1]);
gyCh3 = com1 - com2 + 2 * (smooth_b[(i + 1) * width + j] - smooth_b[(i - 1) * width + j]);
}
int gxx = gxCh1 * gxCh1 + gxCh2 * gxCh2 + gxCh3 * gxCh3;
int gyy = gyCh1 * gyCh1 + gyCh2 * gyCh2 + gyCh3 * gyCh3;
int gxy = gxCh1 * gyCh1 + gxCh2 * gyCh2 + gxCh3 * gyCh3;
#if 1
// Di Zenzo's formulas from Gonzales & Woods - Page 337
double theta = atan2(2.0 * gxy, (double)(gxx - gyy)) / 2; // Gradient Direction
int grad = (int)(sqrt(((gxx + gyy) + (gxx - gyy) * cos(2 * theta) + 2 * gxy * sin(2 * theta)) / 2.0) + 0.5); // Gradient Magnitude
#else
// Koschan & Abidi - 2005 - Signal Processing Magazine
double theta = atan2(2.0 * gxy, (double)(gxx - gyy)) / 2; // Gradient Direction
double cosTheta = cos(theta);
double sinTheta = sin(theta);
int grad = (int)(sqrt(gxx * cosTheta * cosTheta + 2 * gxy * sinTheta * cosTheta + gyy * sinTheta * sinTheta) + 0.5); // Gradient Magnitude
#endif
// Gradient is perpendicular to the edge passing through the pixel
if (theta >= -3.14159 / 4 && theta <= 3.14159 / 4)
dirImg[i * width + j] = EDGE_VERTICAL;
else
dirImg[i * width + j] = EDGE_HORIZONTAL;
gradImg[i * width + j] = (ushort)grad;
if (grad > max)
max = grad;
}
}
// Scale the gradient values to 0-255
double scale = 255.0 / max;
for (int i = 0; i < width * height; i++)
gradImg[i] = (ushort)(gradImg[i] * scale);
}
void EdgeDrawingImpl::smoothChannel(const Mat& src, Mat& dst, double sigma)
{
if (sigma == 1.0)
GaussianBlur(src, dst, Size(5, 5), 1);
else if (sigma == 1.5)
GaussianBlur(src, dst, Size(7, 7), 1.5); // seems to be better?
else
GaussianBlur(src, dst, Size(), sigma);
}
//---------------------------------------------------------
// Fix edge segments having one or two pixel fluctuations
// An example one pixel problem getting fixed:
// x
// x x --> xxx
//
// An example two pixel problem getting fixed:
// xx
// x x --> xxxx
//
void EdgeDrawingImpl::fixEdgeSegments(std::vector<std::vector<cv::Point>> map)
{
/// First fix one pixel problems: There are four cases
for (int i = 0; i < (int)map.size(); i++) {
int cp = (int)map[i].size() - 2; // Current pixel index
int n2 = 0; // next next pixel index
while (n2 < (int)map[i].size()) {
int n1 = cp + 1; // next pixel
cp = cp % map[i].size(); // Roll back to the beginning
n1 = n1 % map[i].size(); // Roll back to the beginning
int r = map[i][cp].y;
int c = map[i][cp].x;
int r1 = map[i][n1].y;
int c1 = map[i][n1].x;
int r2 = map[i][n2].y;
int c2 = map[i][n2].x;
// 4 cases to fix
if (r2 == r - 2 && c2 == c) {
if (c1 != c) {
map[i][n1].x = c;
}
cp = n2;
n2 += 2;
}
else if (r2 == r + 2 && c2 == c) {
if (c1 != c) {
map[i][n1].x = c;
}
cp = n2;
n2 += 2;
}
else if (r2 == r && c2 == c - 2) {
if (r1 != r) {
map[i][n1].y = r;
}
cp = n2;
n2 += 2;
}
else if (r2 == r && c2 == c + 2) {
if (r1 != r) {
map[i][n1].y = r;
}
cp = n2;
n2 += 2;
}
else {
cp++;
n2++;
}
}
}
}
void EdgeDrawingImpl::InitColorEDLib()
{
if (LUT_Initialized)
return;
double inc = 1.0 / LUT_SIZE;
for (int i = 0; i <= LUT_SIZE; i++) {
double d = i * inc;
if (d >= 0.04045) LUT1[i] = pow(((d + 0.055) / 1.055), 2.4);
else LUT1[i] = d / 12.92;
}
inc = 1.0 / LUT_SIZE;
for (int i = 0; i <= LUT_SIZE; i++) {
double d = i * inc;
if (d > 0.008856) LUT2[i] = pow(d, 1.0 / 3.0);
else LUT2[i] = (7.787 * d) + (16.0 / 116.0);
}
LUT_Initialized = true;
}
bool EdgeDrawingImpl::LUT_Initialized = false;
double EdgeDrawingImpl::LUT1[LUT_SIZE + 1] = { 0 };
double EdgeDrawingImpl::LUT2[LUT_SIZE + 1] = { 0 };
} // namespace cv
} // namespace ximgproc

@ -13,15 +13,15 @@
#define ANCHOR_PIXEL 254
#define EDGE_PIXEL 255
#define LEFT 1
#define RIGHT 2
#define UP 3
#define DOWN 4
#define ED_LEFT 1
#define ED_RIGHT 2
#define ED_UP 3
#define ED_DOWN 4
#define SOUTH_SOUTH 0
#define SOUTH_EAST 1
#define EAST_SOUTH 2
#define EAST_EAST 3
#define ED_SOUTH_SOUTH 0
#define ED_SOUTH_EAST 1
#define ED_EAST_SOUTH 2
#define ED_EAST_EAST 3
// Circular arc, circle thresholds
#define VERY_SHORT_ARC_ERROR 0.40 // Used for very short arcs (>= CANDIDATE_CIRCLE_RATIO1 && < CANDIDATE_CIRCLE_RATIO2)
@ -37,39 +37,45 @@
// Ellipse thresholds
#define CANDIDATE_ELLIPSE_RATIO 0.50 // 50% -- If 50% of the ellipse is detected, it may be candidate for validation
#define ELLIPSE_ERROR 1.50 // Used for ellipses. (used to be 1.65 for what reason?)
#define MAX_GRAD_VALUE 128*256
#define ED_MAX_GRAD_VALUE 128*256
using namespace std;
using namespace cv;
#define TABSIZE 100000
#define MAX_LUT_SIZE 1024
#define RELATIVE_ERROR_FACTOR 100.0
class NFALUT
{
public:
NFALUT(int size, double _prob, int _w, int _h);
NFALUT(int size, double _prob, double _logNT);
~NFALUT();
int* LUT; // look up table
int LUTSize;
double prob;
int w, h;
double logNT;
bool checkValidationByNFA(int n, int k);
static double myAtan2(double yy, double xx);
private:
double nfa(int n, int k);
static double Comb(double n, double k);
static double log_gamma_lanczos(double x);
static double log_gamma_windschitl(double x);
static double log_gamma(double x);
static int double_equal(double a, double b);
};
NFALUT::NFALUT(int size, double _prob, int _w, int _h)
NFALUT::NFALUT(int size, double _prob, double _logNT)
{
LUTSize = size > 60 ? 60 : size;
LUT = new int[LUTSize];
w = _w;
h = _h;
prob = _prob;
logNT = _logNT;
LUT[0] = 1;
int j = 1;
@ -77,24 +83,23 @@ NFALUT::NFALUT(int size, double _prob, int _w, int _h)
{
LUT[i] = LUTSize + 1;
double ret = nfa(i, j);
if (ret >= 1.0)
if (ret < 0)
{
while (j < i)
{
j++;
ret = nfa(i, j);
if (ret <= 1.0)
if (ret >= 0)
break;
}
if (ret >= 1.0)
if (ret < 0)
continue;
}
LUT[i] = j;
}
}
NFALUT::~NFALUT()
{
delete[] LUT;
@ -103,7 +108,7 @@ NFALUT::~NFALUT()
bool NFALUT::checkValidationByNFA(int n, int k)
{
if (n >= LUTSize)
return nfa(n, k) <= 1.0;
return nfa(n, k) >= 0.0;
else
return k >= LUT[n];
}
@ -118,28 +123,140 @@ double NFALUT::myAtan2(double yy, double xx)
return angle / 180 * CV_PI;
}
double NFALUT::Comb(double n, double k) //fast combination computation
double NFALUT::nfa(int n, int k)
{
if (k > n)
return 0;
static double inv[TABSIZE]; /* table to keep computed inverse values */
double tolerance = 0.1; /* an error of 10% in the result is accepted */
double log1term, term, bin_term, mult_term, bin_tail, err, p_term;
int i;
/* check parameters */
if (n < 0 || k<0 || k>n || prob <= 0.0 || prob >= 1.0) return -1.0;
/* trivial cases */
if (n == 0 || k == 0) return -logNT;
if (n == k) return -logNT - (double)n * log10(prob);
/* probability term */
p_term = prob / (1.0 - prob);
/* compute the first term of the series */
/*
binomial_tail(n,k,p) = sum_{i=k}^n bincoef(n,i) * p^i * (1-p)^{n-i}
where bincoef(n,i) are the binomial coefficients.
But
bincoef(n,k) = gamma(n+1) / ( gamma(k+1) * gamma(n-k+1) ).
We use this to compute the first term. Actually the log of it.
*/
log1term = log_gamma((double)n + 1.0) - log_gamma((double)k + 1.0)
- log_gamma((double)(n - k) + 1.0)
+ (double)k * log(prob) + (double)(n - k) * log(1.0 - prob);
term = exp(log1term);
/* in some cases no more computations are needed */
if (double_equal(term, 0.0)) { /* the first term is almost zero */
if ((double)k > (double)n * prob) /* at begin or end of the tail? */
return -log1term / M_LN10 - logNT; /* end: use just the first term */
else
return -logNT; /* begin: the tail is roughly 1 */
}
/* compute more terms if needed */
bin_tail = term;
for (i = k + 1; i <= n; i++) {
/*
As
term_i = bincoef(n,i) * p^i * (1-p)^(n-i)
and
bincoef(n,i)/bincoef(n,i-1) = n-1+1 / i,
then,
term_i / term_i-1 = (n-i+1)/i * p/(1-p)
and
term_i = term_i-1 * (n-i+1)/i * p/(1-p).
1/i is stored in a table as they are computed,
because divisions are expensive.
p/(1-p) is computed only once and stored in 'p_term'.
*/
bin_term = (double)(n - i + 1) * (i < TABSIZE ?
(inv[i] != 0.0 ? inv[i] : (inv[i] = 1.0 / (double)i)) :
1.0 / (double)i);
mult_term = bin_term * p_term;
term *= mult_term;
bin_tail += term;
if (bin_term < 1.0) {
/* When bin_term<1 then mult_term_j<mult_term_i for j>i.
Then, the error on the binomial tail when truncated at
the i term can be bounded by a geometric series of form
term_i * sum mult_term_i^j. */
err = term * ((1.0 - pow(mult_term, (double)(n - i + 1))) /
(1.0 - mult_term) - 1.0);
/* One wants an error at most of tolerance*final_result, or:
tolerance * abs(-log10(bin_tail)-logNT).
Now, the error that can be accepted on bin_tail is
given by tolerance*final_result divided by the derivative
of -log10(x) when x=bin_tail. that is:
tolerance * abs(-log10(bin_tail)-logNT) / (1/bin_tail)
Finally, we truncate the tail if the error is less than:
tolerance * abs(-log10(bin_tail)-logNT) * bin_tail */
if (err < tolerance * fabs(-log10(bin_tail) - logNT) * bin_tail) break;
}
}
return -log10(bin_tail) - logNT;
}
double r = 1;
for (double d = 1; d <= k; ++d)
double NFALUT::log_gamma_lanczos(double x)
{
static double q[7] = { 75122.6331530, 80916.6278952, 36308.2951477,
8687.24529705, 1168.92649479, 83.8676043424,
2.50662827511 };
double a = (x + 0.5) * log(x + 5.5) - (x + 5.5);
double b = 0.0;
int n;
for (n = 0; n < 7; n++)
{
r *= n--;
r /= d;
a -= log(x + (double)n);
b += q[n] * pow(x, (double)n);
}
return r;
return a + log(b);
}
double NFALUT::nfa(int n, int k)
double NFALUT::log_gamma_windschitl(double x)
{
double sum = 0;
double p = 0.125;
for (int i = k; i <= n; i++)
sum += Comb(n, i) * pow(p, i) * pow(1 - p, n - i);
return 0.918938533204673 + (x - 0.5) * log(x) - x
+ 0.5 * x * log(x * sinh(1 / x) + 1 / (810.0 * pow(x, 6.0)));
}
double NFALUT::log_gamma(double x)
{
return x > 15 ? log_gamma_windschitl(x) : log_gamma_lanczos(x);
}
int NFALUT::double_equal(double a, double b)
{
double abs_diff, aa, bb, abs_max;
return sum * w * w * h * h;
/* trivial case */
if (a == b) return true;
abs_diff = fabs(a - b);
aa = fabs(a);
bb = fabs(b);
abs_max = aa > bb ? aa : bb;
/* DBL_MIN is the smallest normalized number, thus, the smallest
number whose relative error is bounded by DBL_EPSILON. For
smaller numbers, the same quantization steps as for DBL_MIN
are used. Then, for smaller numbers, a meaningful "relative"
error should be computed by dividing the difference by DBL_MIN. */
if (abs_max < DBL_MIN) abs_max = DBL_MIN;
/* equal if relative error <= factor x eps */
return (abs_diff / abs_max) <= (RELATIVE_ERROR_FACTOR * DBL_EPSILON);
}
struct StackNode
@ -220,18 +337,19 @@ struct mEllipse {
// Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
//
struct EllipseEquation {
double coeff[7]; // coeff[1] = A
EllipseEquation() {
for (int i = 0; i<7; i++) coeff[i] = 0;
}
double A() { return coeff[1]; }
double B() { return coeff[2]; }
double C() { return coeff[3]; }
double D() { return coeff[4]; }
double E() { return coeff[5]; }
double F() { return coeff[6]; }
static constexpr int COEFF_SIZE = 7;
double coeff[COEFF_SIZE]; // coeff[1] = A
// Constructor using an initialization list
EllipseEquation() : coeff{} {}
// Accessor functions marked as const
double A() const { return coeff[1]; }
double B() const { return coeff[2]; }
double C() const { return coeff[3]; }
double D() const { return coeff[4]; }
double E() const { return coeff[5]; }
double F() const { return coeff[6]; }
};
// ================================ CIRCLES ================================

@ -7,7 +7,7 @@ namespace opencv_test { namespace {
const Size img_size(320, 240);
const int FLD_TEST_SEED = 0x134679;
const int EPOCHS = 10;
const int EPOCHS = 5;
class FLDBase : public testing::Test
{
@ -44,6 +44,7 @@ class ximgproc_ED: public FLDBase
{
detector = createEdgeDrawing();
}
string filename = cvtest::TS::ptr()->get_data_path() + "cv/imgproc/beads.jpg";
protected:
Ptr<EdgeDrawing> detector;
@ -275,16 +276,18 @@ TEST_F(ximgproc_ED, rotatedRect)
ASSERT_EQ(EPOCHS, passedtests);
}
TEST_F(ximgproc_ED, ManySmallCircles)
TEST_F(ximgproc_ED, detectLinesAndEllipses)
{
string picture_name = "cv/imgproc/beads.jpg";
Mat gray_image;
vector<Vec6d> ellipses;
string filename = cvtest::TS::ptr()->get_data_path() + picture_name;
test_image = imread(filename, IMREAD_GRAYSCALE);
test_image = imread(filename);
EXPECT_FALSE(test_image.empty()) << "Invalid test image: " << filename;
vector<Vec6d> ellipses;
detector->detectEdges(test_image);
cvtColor(test_image, test_image, COLOR_BGR2BGRA);
cvtColor(test_image, gray_image, COLOR_BGR2GRAY);
detector->detectEdges(gray_image);
detector->detectEllipses(ellipses);
detector->detectLines(lines);
@ -295,5 +298,34 @@ TEST_F(ximgproc_ED, ManySmallCircles)
EXPECT_GE(lines.size(), lines_size);
EXPECT_LE(lines.size(), lines_size + 2);
EXPECT_EQ(ellipses.size(), ellipses_size);
detector->params.PFmode = true;
detector->detectEdges(gray_image);
detector->detectEllipses(ellipses);
detector->detectLines(lines);
segments_size = 2717;
lines_size = 6197;
ellipses_size = 2446;
EXPECT_EQ(detector->getSegments().size(), segments_size);
EXPECT_GE(lines.size(), lines_size);
EXPECT_LE(lines.size(), lines_size + 2);
EXPECT_EQ(ellipses.size(), ellipses_size);
detector->params.MinLineLength = 10;
detector->detectEdges(test_image);
detector->detectEllipses(ellipses);
detector->detectLines(lines);
detector->detectEllipses(ellipses);
segments_size = 6230;
lines_size = 11133;
ellipses_size = 2431;
EXPECT_EQ(detector->getSegments().size(), segments_size);
EXPECT_GE(lines.size(), lines_size);
EXPECT_LE(lines.size(), lines_size + 2);
EXPECT_GE(ellipses.size(), ellipses_size);
EXPECT_LE(ellipses.size(), ellipses_size + 2);
}
}} // namespace

Loading…
Cancel
Save