diff --git a/ED.cpp b/ED.cpp new file mode 100644 index 0000000..d2284c1 --- /dev/null +++ b/ED.cpp @@ -0,0 +1,1086 @@ +#include "ED.h" +#include "EDColor.h" +#include + +using namespace cv; +using namespace std; + +ED::ED(Mat _srcImage, GradientOperator _op, int _gradThresh, int _anchorThresh,int _scanInterval, int _minPathLen ,double _sigma, bool _sumFlag) +{ + // Check parameters for sanity + if (_gradThresh < 1) _gradThresh = 1; + if (_anchorThresh < 0) _anchorThresh = 0; + if (_sigma < 1.0) _sigma = 1.0; + + srcImage = _srcImage; + + height = srcImage.rows; + width = srcImage.cols; + + op = _op; + gradThresh = _gradThresh; + anchorThresh = _anchorThresh; + scanInterval = _scanInterval; + minPathLen = _minPathLen; + sigma = _sigma; + sumFlag = _sumFlag; + + segmentNos = 0; + segmentPoints.push_back(vector()); // create empty vector of points for segments + + edgeImage = Mat(height, width, CV_8UC1, Scalar(0)); // initialize edge Image + smoothImage = Mat(height, width, CV_8UC1); + gradImage = Mat(height, width, CV_16SC1); // gradImage contains short values + + srcImg = srcImage.data; + + //// Detect Edges By Edge Drawing Algorithm //// + + /*------------ SMOOTH THE IMAGE BY A GAUSSIAN KERNEL -------------------*/ + if (sigma == 1.0) + GaussianBlur(srcImage, smoothImage, Size(5, 5), sigma); + else + GaussianBlur(srcImage, smoothImage, Size(), sigma); // calculate kernel from sigma + + // Assign Pointers from Mat's data + smoothImg = smoothImage.data; + gradImg = (short*)gradImage.data; + edgeImg = edgeImage.data; + + dirImg = new unsigned char[width*height]; + + /*------------ COMPUTE GRADIENT & EDGE DIRECTION MAPS -------------------*/ + ComputeGradient(); + + /*------------ COMPUTE ANCHORS -------------------*/ + ComputeAnchorPoints(); + + /*------------ JOIN ANCHORS -------------------*/ + JoinAnchorPointsUsingSortedAnchors(); + + delete[] dirImg; +} + +// This constructor for use of EDLines and EDCircle with ED given as constructor argument +// only the necessary attributes are coppied +ED::ED(const ED & cpyObj) +{ + height = cpyObj.height; + width = cpyObj.width; + + srcImage = cpyObj.srcImage.clone(); + + op = cpyObj.op; + gradThresh = cpyObj.gradThresh; + anchorThresh = cpyObj.anchorThresh; + scanInterval = cpyObj.scanInterval; + minPathLen = cpyObj.minPathLen; + sigma = cpyObj.sigma; + sumFlag = cpyObj.sumFlag; + + edgeImage = cpyObj.edgeImage.clone(); + smoothImage = cpyObj.smoothImage.clone(); + gradImage = cpyObj.gradImage.clone(); + + srcImg = srcImage.data; + + smoothImg = smoothImage.data; + gradImg = (short*)gradImage.data; + edgeImg = edgeImage.data; + + segmentPoints = cpyObj.segmentPoints; + segmentNos = cpyObj.segmentNos; +} + +// This constructor for use of EDColor with use of direction and gradient image +// It finds edge image for given gradient and direction image +ED::ED(short *_gradImg, uchar *_dirImg, int _width, int _height, int _gradThresh, int _anchorThresh, int _scanInterval, int _minPathLen, bool selectStableAnchors) +{ + height = _height; + width = _width; + + gradThresh = _gradThresh; + anchorThresh = _anchorThresh; + scanInterval = _scanInterval; + minPathLen = _minPathLen; + + gradImg = _gradImg; + dirImg = _dirImg; + + edgeImage = Mat(height, width, CV_8UC1, Scalar(0)); // initialize edge Image + + edgeImg = edgeImage.data; + + if (selectStableAnchors) { + + // 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(); + anchorThresh = _anchorThresh; // 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= anchorThresh && diff2 >= anchorThresh) edgeImg[i*width + j] = 255; + + continue; + } //end-if + + // 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; + + continue; + } //end-if + + // 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; + } //end-if + + // 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; + } //end-if + + } //end-for + } //end-for + + for (int i = 0; i()); // create empty vector of points for segments + + JoinAnchorPointsUsingSortedAnchors(); +} + +ED::ED(EDColor &obj) +{ + width = obj.getWidth(); + height = obj.getHeight(); + segmentPoints = obj.getSegments(); + segmentNos = obj.getSegmentNo(); +} + +ED::ED() +{ + // +} + + +Mat ED::getEdgeImage() +{ + return edgeImage; +} + +Mat ED::getAnchorImage() +{ + Mat anchorImage = Mat(edgeImage.size(), edgeImage.type(), Scalar(0)); + + std::vector::iterator it; + + for (it = anchorPoints.begin(); it != anchorPoints.end(); it++) + anchorImage.at(*it) = 255; + + return anchorImage; +} + +Mat ED::getSmoothImage() +{ + return smoothImage; +} + +Mat ED::getGradImage() +{ + Mat result8UC1; + convertScaleAbs(gradImage, result8UC1); + + return result8UC1; +} + +int ED::getSegmentNo() +{ + return segmentNos; +} + +int ED::getAnchorNo() +{ + return anchorNos; +} + +std::vector ED::getAnchorPoints() +{ + return anchorPoints; +} + +std::vector> ED::getSegments() +{ + return segmentPoints; +} + +std::vector> ED::getSortedSegments() +{ + // sort segments from largest to smallest + std::sort(segmentPoints.begin(), segmentPoints.end(), [](const std::vector & a, const std::vector & b) { return a.size() > b.size(); }); + + return segmentPoints; +} + +Mat ED::drawParticularSegments(std::vector list) +{ + Mat segmentsImage = Mat(edgeImage.size(), edgeImage.type(), Scalar(0)); + + std::vector::iterator it; + std::vector::iterator itInt; + + for (itInt = list.begin(); itInt != list.end(); itInt++) + for (it = segmentPoints[*itInt].begin(); it != segmentPoints[*itInt].end(); it++) + segmentsImage.at(*it) = 255; + + return segmentsImage; +} + + +void ED::ComputeGradient() +{ + // Initialize gradient image for row = 0, row = height-1, column=0, column=width-1 + for (int j = 0; j= gradThresh) { + if (gx >= gy) dirImg[index] = EDGE_VERTICAL; + else dirImg[index] = EDGE_HORIZONTAL; + } //end-if + } // end-for + } // end-for +} + +void ED::ComputeAnchorPoints() +{ + //memset(edgeImg, 0, width*height); + for (int i = 2; i= anchorThresh && diff2 >= anchorThresh) { + edgeImg[i*width + j] = ANCHOR_PIXEL; + anchorPoints.push_back(Point(j, i)); + } + + } + else { + // horizontal edge + 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] = ANCHOR_PIXEL; + anchorPoints.push_back(Point(j, i)); + } + } // end-else + } //end-for-inner + } //end-for-outer + + anchorNos = anchorPoints.size(); // get the total number of anchor points +} + +void ED::JoinAnchorPointsUsingSortedAnchors() +{ + int *chainNos = new int[(width + height) * 8]; + + Point *pixels = new Point[width*height]; + StackNode *stack = new StackNode[width*height]; + Chain *chains = new Chain[width*height]; + + // sort the anchor points by their gradient value in decreasing order + int *A = sortAnchorsByGradValue1(); + + // Now join the anchors starting with the anchor having the greatest gradient value + int totalPixels = 0; + + for (int k = anchorNos - 1; k >= 0; k--) { + int pixelOffset = A[k]; + + int i = pixelOffset / width; + int j = pixelOffset % width; + + //int i = anchorPoints[k].y; + //int j = anchorPoints[k].x; + + if (edgeImg[i*width + j] != ANCHOR_PIXEL) continue; + + chains[0].len = 0; + chains[0].parent = -1; + chains[0].dir = 0; + chains[0].children[0] = chains[0].children[1] = -1; + chains[0].pixels = NULL; + + + int noChains = 1; + int len = 0; + int duplicatePixelCount = 0; + int top = -1; // top of the stack + + if (dirImg[i*width + j] == EDGE_VERTICAL) { + stack[++top].r = i; + stack[top].c = j; + stack[top].dir = DOWN; + stack[top].parent = 0; + + stack[++top].r = i; + stack[top].c = j; + stack[top].dir = UP; + stack[top].parent = 0; + + } + else { + stack[++top].r = i; + stack[top].c = j; + stack[top].dir = RIGHT; + stack[top].parent = 0; + + stack[++top].r = i; + stack[top].c = j; + stack[top].dir = LEFT; + stack[top].parent = 0; + } //end-else + + // While the stack is not empty + StartOfWhile: + while (top >= 0) { + int r = stack[top].r; + int c = stack[top].c; + int dir = stack[top].dir; + int parent = stack[top].parent; + top--; + + if (edgeImg[r*width + c] != EDGE_PIXEL) duplicatePixelCount++; + + chains[noChains].dir = dir; // traversal direction + chains[noChains].parent = parent; + chains[noChains].children[0] = chains[noChains].children[1] = -1; + + + int chainLen = 0; + + chains[noChains].pixels = &pixels[len]; + + pixels[len].y = r; + pixels[len].x = c; + len++; + chainLen++; + + if (dir == LEFT) { + while (dirImg[r*width + c] == EDGE_HORIZONTAL) { + edgeImg[r*width + c] = EDGE_PIXEL; + + // The edge is horizontal. Look LEFT + // + // A + // B x + // C + // + // cleanup up & down pixels + if (edgeImg[(r - 1)*width + c] == ANCHOR_PIXEL) edgeImg[(r - 1)*width + c] = 0; + if (edgeImg[(r + 1)*width + c] == ANCHOR_PIXEL) edgeImg[(r + 1)*width + c] = 0; + + // Look if there is an edge pixel in the neighbors + if (edgeImg[r*width + c - 1] >= ANCHOR_PIXEL) { c--; } + else if (edgeImg[(r - 1)*width + c - 1] >= ANCHOR_PIXEL) { r--; c--; } + else if (edgeImg[(r + 1)*width + c - 1] >= ANCHOR_PIXEL) { r++; c--; } + else { + // else -- follow max. pixel to the LEFT + int A = gradImg[(r - 1)*width + c - 1]; + int B = gradImg[r*width + c - 1]; + int C = gradImg[(r + 1)*width + c - 1]; + + if (A > B) { + if (A > C) r--; + else r++; + } + else if (C > B) r++; + c--; + } //end-else + + if (edgeImg[r*width + c] == EDGE_PIXEL || gradImg[r*width + c] < gradThresh) { + if (chainLen > 0) { + chains[noChains].len = chainLen; + chains[parent].children[0] = noChains; + noChains++; + } // end-if + goto StartOfWhile; + } //end-else + + + pixels[len].y = r; + pixels[len].x = c; + len++; + chainLen++; + } //end-while + + stack[++top].r = r; + stack[top].c = c; + stack[top].dir = DOWN; + stack[top].parent = noChains; + + stack[++top].r = r; + stack[top].c = c; + stack[top].dir = UP; + stack[top].parent = noChains; + + len--; + chainLen--; + + chains[noChains].len = chainLen; + chains[parent].children[0] = noChains; + noChains++; + + } + else if (dir == RIGHT) { + while (dirImg[r*width + c] == EDGE_HORIZONTAL) { + edgeImg[r*width + c] = EDGE_PIXEL; + + // The edge is horizontal. Look RIGHT + // + // A + // x B + // C + // + // cleanup up&down pixels + if (edgeImg[(r + 1)*width + c] == ANCHOR_PIXEL) edgeImg[(r + 1)*width + c] = 0; + if (edgeImg[(r - 1)*width + c] == ANCHOR_PIXEL) edgeImg[(r - 1)*width + c] = 0; + + // Look if there is an edge pixel in the neighbors + if (edgeImg[r*width + c + 1] >= ANCHOR_PIXEL) { c++; } + else if (edgeImg[(r + 1)*width + c + 1] >= ANCHOR_PIXEL) { r++; c++; } + else if (edgeImg[(r - 1)*width + c + 1] >= ANCHOR_PIXEL) { r--; c++; } + else { + // else -- follow max. pixel to the RIGHT + int A = gradImg[(r - 1)*width + c + 1]; + int B = gradImg[r*width + c + 1]; + int C = gradImg[(r + 1)*width + c + 1]; + + if (A > B) { + if (A > C) r--; // A + else r++; // C + } + else if (C > B) r++; // C + c++; + } //end-else + + if (edgeImg[r*width + c] == EDGE_PIXEL || gradImg[r*width + c] < gradThresh) { + if (chainLen > 0) { + chains[noChains].len = chainLen; + chains[parent].children[1] = noChains; + noChains++; + } // end-if + goto StartOfWhile; + } //end-else + + + pixels[len].y = r; + pixels[len].x = c; + len++; + chainLen++; + } //end-while + + stack[++top].r = r; + stack[top].c = c; + stack[top].dir = DOWN; // Go down + stack[top].parent = noChains; + + stack[++top].r = r; + stack[top].c = c; + stack[top].dir = UP; // Go up + stack[top].parent = noChains; + + len--; + chainLen--; + + chains[noChains].len = chainLen; + chains[parent].children[1] = noChains; + noChains++; + + } + else if (dir == UP) { + while (dirImg[r*width + c] == EDGE_VERTICAL) { + edgeImg[r*width + c] = EDGE_PIXEL; + + // The edge is vertical. Look UP + // + // A B C + // x + // + // Cleanup left & right pixels + if (edgeImg[r*width + c - 1] == ANCHOR_PIXEL) edgeImg[r*width + c - 1] = 0; + if (edgeImg[r*width + c + 1] == ANCHOR_PIXEL) edgeImg[r*width + c + 1] = 0; + + // Look if there is an edge pixel in the neighbors + if (edgeImg[(r - 1)*width + c] >= ANCHOR_PIXEL) { r--; } + else if (edgeImg[(r - 1)*width + c - 1] >= ANCHOR_PIXEL) { r--; c--; } + else if (edgeImg[(r - 1)*width + c + 1] >= ANCHOR_PIXEL) { r--; c++; } + else { + // else -- follow the max. pixel UP + int A = gradImg[(r - 1)*width + c - 1]; + int B = gradImg[(r - 1)*width + c]; + int C = gradImg[(r - 1)*width + c + 1]; + + if (A > B) { + if (A > C) c--; + else c++; + } + else if (C > B) c++; + r--; + } //end-else + + if (edgeImg[r*width + c] == EDGE_PIXEL || gradImg[r*width + c] < gradThresh) { + if (chainLen > 0) { + chains[noChains].len = chainLen; + chains[parent].children[0] = noChains; + noChains++; + } // end-if + goto StartOfWhile; + } //end-else + + + pixels[len].y = r; + pixels[len].x = c; + + len++; + chainLen++; + } //end-while + + stack[++top].r = r; + stack[top].c = c; + stack[top].dir = RIGHT; + stack[top].parent = noChains; + + stack[++top].r = r; + stack[top].c = c; + stack[top].dir = LEFT; + stack[top].parent = noChains; + + len--; + chainLen--; + + chains[noChains].len = chainLen; + chains[parent].children[0] = noChains; + noChains++; + + } + else { // dir == DOWN + while (dirImg[r*width + c] == EDGE_VERTICAL) { + edgeImg[r*width + c] = EDGE_PIXEL; + + // The edge is vertical + // + // x + // A B C + // + // cleanup side pixels + if (edgeImg[r*width + c + 1] == ANCHOR_PIXEL) edgeImg[r*width + c + 1] = 0; + if (edgeImg[r*width + c - 1] == ANCHOR_PIXEL) edgeImg[r*width + c - 1] = 0; + + // Look if there is an edge pixel in the neighbors + if (edgeImg[(r + 1)*width + c] >= ANCHOR_PIXEL) { r++; } + else if (edgeImg[(r + 1)*width + c + 1] >= ANCHOR_PIXEL) { r++; c++; } + else if (edgeImg[(r + 1)*width + c - 1] >= ANCHOR_PIXEL) { r++; c--; } + else { + // else -- follow the max. pixel DOWN + int A = gradImg[(r + 1)*width + c - 1]; + int B = gradImg[(r + 1)*width + c]; + int C = gradImg[(r + 1)*width + c + 1]; + + if (A > B) { + if (A > C) c--; // A + else c++; // C + } + else if (C > B) c++; // C + r++; + } //end-else + + if (edgeImg[r*width + c] == EDGE_PIXEL || gradImg[r*width + c] < gradThresh) { + if (chainLen > 0) { + chains[noChains].len = chainLen; + chains[parent].children[1] = noChains; + noChains++; + } // end-if + goto StartOfWhile; + } //end-else + + pixels[len].y = r; + pixels[len].x = c; + + len++; + chainLen++; + } //end-while + + stack[++top].r = r; + stack[top].c = c; + stack[top].dir = RIGHT; + stack[top].parent = noChains; + + stack[++top].r = r; + stack[top].c = c; + stack[top].dir = LEFT; + stack[top].parent = noChains; + + len--; + chainLen--; + + chains[noChains].len = chainLen; + chains[parent].children[1] = noChains; + noChains++; + } //end-else + + } //end-while + + + if (len - duplicatePixelCount < minPathLen) { + for (int k = 0; k 0) { + // Retrieve the chainNos + int count = RetrieveChainNos(chains, chains[0].children[1], chainNos); + + // Copy these pixels in the reverse order + for (int k = count - 1; k >= 0; k--) { + int chainNo = chainNos[k]; + +#if 1 + /* See if we can erase some pixels from the last chain. This is for cleanup */ + + int fr = chains[chainNo].pixels[chains[chainNo].len - 1].y; + int fc = chains[chainNo].pixels[chains[chainNo].len - 1].x; + + int index = noSegmentPixels - 2; + while (index >= 0) { + int dr = abs(fr - segmentPoints[segmentNos][index].y); + int dc = abs(fc - segmentPoints[segmentNos][index].x); + + if (dr <= 1 && dc <= 1) { + // neighbors. Erase last pixel + segmentPoints[segmentNos].pop_back(); + noSegmentPixels--; + index--; + } + else break; + } //end-while + + if (chains[chainNo].len > 1 && noSegmentPixels > 0) { + fr = chains[chainNo].pixels[chains[chainNo].len - 2].y; + fc = chains[chainNo].pixels[chains[chainNo].len - 2].x; + + int dr = abs(fr - segmentPoints[segmentNos][noSegmentPixels - 1].y); + int dc = abs(fc - segmentPoints[segmentNos][noSegmentPixels - 1].x); + + if (dr <= 1 && dc <= 1) chains[chainNo].len--; + } //end-if +#endif + + for (int l = chains[chainNo].len - 1; l >= 0; l--) { + segmentPoints[segmentNos].push_back(chains[chainNo].pixels[l]); + noSegmentPixels++; + } //end-for + + chains[chainNo].len = 0; // Mark as copied + } //end-for + } //end-if + + totalLen = LongestChain(chains, chains[0].children[0]); + if (totalLen > 1) { + // Retrieve the chainNos + int count = RetrieveChainNos(chains, chains[0].children[0], chainNos); + + // Copy these chains in the forward direction. Skip the first pixel of the first chain + // due to repetition with the last pixel of the previous chain + int lastChainNo = chainNos[0]; + chains[lastChainNo].pixels++; + chains[lastChainNo].len--; + + for (int k = 0; k= 0) { + int dr = abs(fr - segmentPoints[segmentNos][index].y); + int dc = abs(fc - segmentPoints[segmentNos][index].x); + + if (dr <= 1 && dc <= 1) { + // neighbors. Erase last pixel + segmentPoints[segmentNos].pop_back(); + noSegmentPixels--; + index--; + } + else break; + } //end-while + + int startIndex = 0; + int chainLen = chains[chainNo].len; + if (chainLen > 1 && noSegmentPixels > 0) { + int fr = chains[chainNo].pixels[1].y; + int fc = chains[chainNo].pixels[1].x; + + int dr = abs(fr - segmentPoints[segmentNos][noSegmentPixels - 1].y); + int dc = abs(fc - segmentPoints[segmentNos][noSegmentPixels - 1].x); + + if (dr <= 1 && dc <= 1) { startIndex = 1; } + } //end-if +#endif + + /* Start a new chain & copy pixels from the new chain */ + for (int l = startIndex; l()); // create empty vector of points for segments + + // Copy the rest of the long chains here + for (int k = 2; k= 10) { + + // Retrieve the chainNos + int count = RetrieveChainNos(chains, k, chainNos); + + // Copy the pixels + noSegmentPixels = 0; + for (int k = 0; k= 0) { + int dr = abs(fr - segmentPoints[segmentNos][index].y); + int dc = abs(fc - segmentPoints[segmentNos][index].x); + + if (dr <= 1 && dc <= 1) { + // neighbors. Erase last pixel + segmentPoints[segmentNos].pop_back(); + noSegmentPixels--; + index--; + } + else break; + } //end-while + + int startIndex = 0; + int chainLen = chains[chainNo].len; + if (chainLen > 1 && noSegmentPixels > 0) { + int fr = chains[chainNo].pixels[1].y; + int fc = chains[chainNo].pixels[1].x; + + int dr = abs(fr - segmentPoints[segmentNos][noSegmentPixels - 1].y); + int dc = abs(fc - segmentPoints[segmentNos][noSegmentPixels - 1].x); + + if (dr <= 1 && dc <= 1) { startIndex = 1; } + } //end-if +#endif + /* Start a new chain & copy pixels from the new chain */ + for (int l = startIndex; l()); // create empty vector of points for segments + segmentNos++; + } //end-if + } //end-for + + } //end-else + + } //end-for-outer + + // pop back last segment from vector + // because of one preallocation in the beginning, it will always empty + segmentPoints.pop_back(); + + // Clean up + delete[] A; + delete[] chains; + delete[] stack; + delete[] chainNos; + delete[] pixels; +} + +void ED::sortAnchorsByGradValue() +{ + auto sortFunc = [&](const Point &a, const Point &b) + { + return gradImg[a.y * width + a.x] > gradImg[b.y * width + b.x]; + }; + + std::sort(anchorPoints.begin(), anchorPoints.end(), sortFunc); + + /* + ofstream myFile; + myFile.open("anchorsNew.txt"); + for (int i = 0; i < anchorPoints.size(); i++) { + int x = anchorPoints[i].x; + int y = anchorPoints[i].y; + + myFile << i << ". value: " << gradImg[y*width + x] << " Cord: (" << x << "," << y << ")" << endl; + } + myFile.close(); + + + vector temp(anchorNos); + + int x, y, i = 0; + char c; + std::ifstream infile("cords.txt"); + while (infile >> x >> c >> y && c == ',') { + temp[i] = Point(x, y); + i++; + } + + anchorPoints = temp; + */ +} + +int * ED::sortAnchorsByGradValue1() +{ + int SIZE = 128 * 256; + int *C = new int[SIZE]; + memset(C, 0, sizeof(int)*SIZE); + + // Count the number of grad values + for (int i = 1; i= len1) { + max = len0; + chains[root].children[1] = -1; + + } + else { + max = len1; + chains[root].children[0] = -1; + } //end-else + + return chains[root].len + max; +} //end-LongestChain + +int ED::RetrieveChainNos(Chain * chains, int root, int chainNos[]) +{ + int count = 0; + + while (root != -1) { + chainNos[count] = root; + count++; + + if (chains[root].children[0] != -1) root = chains[root].children[0]; + else root = chains[root].children[1]; + } //end-while + + return count; +} diff --git a/ED.h b/ED.h new file mode 100644 index 0000000..53c40c1 --- /dev/null +++ b/ED.h @@ -0,0 +1,129 @@ +/************************************************************************************************************** +* Edge Drawing (ED) and Edge Drawing Parameter Free (EDPF) source codes. +* Copyright (C) 2016, Cuneyt Akinlar & Cihan Topal +* E-mails of the authors: cuneytakinlar@gmail.com, cihant@anadolu.edu.tr +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. + +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. + +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . + +* By using this library you are implicitly assumed to have accepted all of the above statements, +* and accept to cite the following papers: +* +* [1] C. Topal and C. Akinlar, “Edge Drawing: A Combined Real-Time Edge and Segment Detector,” +* Journal of Visual Communication and Image Representation, 23(6), 862-872, DOI: 10.1016/j.jvcir.2012.05.004 (2012). +* +* [2] C. Akinlar and C. Topal, “EDPF: A Real-time Parameter-free Edge Segment Detector with a False Detection Control,” +* International Journal of Pattern Recognition and Artificial Intelligence, 26(1), DOI: 10.1142/S0218001412550026 (2012). +**************************************************************************************************************/ + +#ifndef _ED_ +#define _ED_ + +#include +#include "EDColor.h" + +/// Special defines +#define EDGE_VERTICAL 1 +#define EDGE_HORIZONTAL 2 + +#define ANCHOR_PIXEL 254 +#define EDGE_PIXEL 255 + +#define LEFT 1 +#define RIGHT 2 +#define UP 3 +#define DOWN 4 + +enum GradientOperator { PREWITT_OPERATOR = 101, SOBEL_OPERATOR = 102, SCHARR_OPERATOR = 103, LSD_OPERATOR = 104 }; + +struct StackNode { + int r, c; // starting pixel + int parent; // parent chain (-1 if no parent) + int dir; // direction where you are supposed to go +}; + +// Used during Edge Linking +struct Chain { + + int dir; // Direction of the chain + int len; // # of pixels in the chain + int parent; // Parent of this node (-1 if no parent) + int children[2]; // Children of this node (-1 if no children) + cv::Point *pixels; // Pointer to the beginning of the pixels array +}; + +class ED { + +public: + ED(cv::Mat _srcImage, GradientOperator _op = PREWITT_OPERATOR, int _gradThresh = 20, int _anchorThresh = 0, int _scanInterval = 1, int _minPathLen = 10, double _sigma = 1.0, bool _sumFlag = true); + ED(const ED &cpyObj); + ED(short* gradImg, uchar *dirImg, int _width, int _height, int _gradThresh, int _anchorThresh, int _scanInterval = 1, int _minPathLen = 10, bool selectStableAnchors = true); + ED(EDColor &cpyObj); + ED(); + + cv::Mat getEdgeImage(); + cv::Mat getAnchorImage(); + cv::Mat getSmoothImage(); + cv::Mat getGradImage(); + + int getSegmentNo(); + int getAnchorNo(); + + std::vector getAnchorPoints(); + std::vector> getSegments(); + std::vector> getSortedSegments(); + + cv::Mat drawParticularSegments(std::vector list); + +protected: + int width; // width of source image + int height; // height of source image + uchar *srcImg; + std::vector > segmentPoints; + double sigma; // Gaussian sigma + cv::Mat smoothImage; + uchar *edgeImg; // pointer to edge image data + uchar *smoothImg; // pointer to smoothed image data + int segmentNos; + int minPathLen; + cv::Mat srcImage; + +private: + void ComputeGradient(); + void ComputeAnchorPoints(); + void JoinAnchorPointsUsingSortedAnchors(); + void sortAnchorsByGradValue(); + int* sortAnchorsByGradValue1(); + + static int LongestChain(Chain *chains, int root); + static int RetrieveChainNos(Chain *chains, int root, int chainNos[]); + + int anchorNos; + std::vector anchorPoints; + std::vector edgePoints; + + cv::Mat edgeImage; + cv::Mat gradImage; + + uchar *dirImg; // pointer to direction image data + short *gradImg; // pointer to gradient image data + + GradientOperator op; // operation used in gradient calculation + int gradThresh; // gradient threshold + int anchorThresh; // anchor point threshold + int scanInterval; + bool sumFlag; +}; + + +#endif \ No newline at end of file diff --git a/EDCircles.cpp b/EDCircles.cpp new file mode 100644 index 0000000..a0a6160 --- /dev/null +++ b/EDCircles.cpp @@ -0,0 +1,3659 @@ +#include "EDCircles.h" + +using namespace cv; +using namespace std; + +EDCircles::EDCircles(Mat srcImage) + : EDPF(srcImage) +{ + // Arcs & circles to be detected + // If the end-points of the segment is very close to each other, + // then directly fit a circle/ellipse instread of line fitting + noCircles1 = 0; + circles1 = new Circle[(width + height) * 8]; + + // ----------------------------------- DETECT LINES --------------------------------- + int bufferSize = 0; + for (int i = 0; i < segmentPoints.size(); i++) + bufferSize += segmentPoints[i].size(); + + // Compute the starting line number for each segment + segmentStartLines = new int[segmentNos + 1]; + + bm = new BufferManager(bufferSize * 8); + vector lines; + + +#define CIRCLE_MIN_LINE_LEN 6 + + for (int i = 0; i < segmentNos; i++) { + + // Make note of the starting line number for this segment + segmentStartLines[i] = lines.size(); + + int noPixels = segmentPoints[i].size(); + + if (noPixels < 2 * CIRCLE_MIN_LINE_LEN) + continue; + + double *x = bm->getX(); + double *y = bm->getY(); + + for (int j = 0; j= 4 * CIRCLE_MIN_LINE_LEN) { + // If the end-points of the segment is close to each other, then assume a circular/elliptic structure + double dx = x[0] - x[noPixels - 1]; + double dy = y[0] - y[noPixels - 1]; + double d = sqrt(dx*dx + dy*dy); + double r = noPixels / TWOPI; // Assume a complete circle + + double maxDistanceBetweenEndPoints = MAX(3, r / 4); + + // If almost closed loop, then try to fit a circle/ellipse + if (d <= maxDistanceBetweenEndPoints) { + double xc, yc, r, circleFitError = 1e10; + + CircleFit(x, y, noPixels, &xc, &yc, &r, &circleFitError); + + EllipseEquation eq; + double ellipseFitError = 1e10; + + if (circleFitError > LONG_ARC_ERROR) { + // Try fitting an ellipse + if (EllipseFit(x, y, noPixels, &eq)) + ellipseFitError = ComputeEllipseError(&eq, x, y, noPixels); + } //end-if + + if (circleFitError <= LONG_ARC_ERROR) { + addCircle(circles1, noCircles1, xc ,yc, r, circleFitError, x, y, noPixels); + bm->move(noPixels); + continue; + + } + else if (ellipseFitError <= ELLIPSE_ERROR) { + double major, minor; + ComputeEllipseCenterAndAxisLengths(&eq, &xc, &yc, &major, &minor); + + // Assume major is longer. Otherwise, swap + if (minor > major) { double tmp = major; major = minor; minor = tmp; } + + if (major < 8 * minor) { + addCircle(circles1, noCircles1,xc, yc, r, circleFitError, &eq, ellipseFitError, x, y, noPixels); + bm->move(noPixels); + } //end-if + + continue; + } //end-else + } //end-if + } //end-if + + // Otherwise, split to lines + EDLines::SplitSegment2Lines(x, y, noPixels, i, lines); + } + + segmentStartLines[segmentNos] = lines.size(); + + // ------------------------------- DETECT ARCS --------------------------------- + + info = new Info[lines.size()]; + + // Compute the angle information for each line segment + for (int i = 0; iex - l2->sx; + double dy = l1->ey - l2->sy; + double d = sqrt(dx*dx + dy*dy); + if (d >= 15) { info[j].angle = 10; info[j].sign = 2; info[j].taken = false; continue; } +#endif + + // Compute the angle between the lines & their turn direction + double v1x = l1->ex - l1->sx; + double v1y = l1->ey - l1->sy; + double v1Len = sqrt(v1x*v1x + v1y*v1y); + + double v2x = l2->ex - l2->sx; + double v2y = l2->ey - l2->sy; + double v2Len = sqrt(v2x*v2x + v2y*v2y); + + double dotProduct = (v1x*v2x + v1y*v2y) / (v1Len*v2Len); + if (dotProduct > 1.0) dotProduct = 1.0; + else if (dotProduct < -1.0) dotProduct = -1.0; + + info[j].angle = acos(dotProduct); + info[j].sign = (v1x*v2y - v2x*v1y) >= 0 ? 1 : -1; // compute cross product + info[j].taken = false; + } //end-for + } //end-for + + // This is how much space we will allocate for circles buffers + int maxNoOfCircles = lines.size() / 3 + noCircles1 * 2; + + edarcs1 = new EDArcs(maxNoOfCircles); + DetectArcs(lines); // Detect all arcs + + // Try to join arcs that are almost perfectly circular. + // Use the distance between the arc end-points as a metric in choosing in choosing arcs to join + edarcs2 = new EDArcs(maxNoOfCircles); + JoinArcs1(); + + // Try to join arcs that belong to the same segment + edarcs3 = new EDArcs(maxNoOfCircles); + JoinArcs2(); + + // Try to combine arcs that belong to different segments + edarcs4 = new EDArcs(maxNoOfCircles); // The remaining arcs + JoinArcs3(); + + // Finally, go over the arcs & circles, and generate candidate circles + GenerateCandidateCircles(); + + //----------------------------- VALIDATE CIRCLES -------------------------- + noCircles2 = 0; + circles2 = new Circle[maxNoOfCircles]; + GaussianBlur(srcImage, smoothImage, Size(), 0.50); // calculate kernel from sigma; + + ValidateCircles(); + + //----------------------------- JOIN CIRCLES -------------------------- + noCircles3 = 0; + circles3 = new Circle[maxNoOfCircles]; + JoinCircles(); + + noCircles = 0; + noEllipses = 0; + for (int i = 0; i < noCircles3; i++) { + if (circles3[i].isEllipse) + noEllipses++; + else + noCircles++; + } + + for (int i = 0; i lines; + + +#define CIRCLE_MIN_LINE_LEN 6 + + for (int i = 0; i < segmentNos; i++) { + + // Make note of the starting line number for this segment + segmentStartLines[i] = lines.size(); + + int noPixels = segmentPoints[i].size(); + + if (noPixels < 2 * CIRCLE_MIN_LINE_LEN) + continue; + + double *x = bm->getX(); + double *y = bm->getY(); + + for (int j = 0; j= 4 * CIRCLE_MIN_LINE_LEN) { + // If the end-points of the segment is close to each other, then assume a circular/elliptic structure + double dx = x[0] - x[noPixels - 1]; + double dy = y[0] - y[noPixels - 1]; + double d = sqrt(dx*dx + dy*dy); + double r = noPixels / TWOPI; // Assume a complete circle + + double maxDistanceBetweenEndPoints = MAX(3, r / 4); + + // If almost closed loop, then try to fit a circle/ellipse + if (d <= maxDistanceBetweenEndPoints) { + double xc, yc, r, circleFitError = 1e10; + + CircleFit(x, y, noPixels, &xc, &yc, &r, &circleFitError); + + EllipseEquation eq; + double ellipseFitError = 1e10; + + if (circleFitError > LONG_ARC_ERROR) { + // Try fitting an ellipse + if (EllipseFit(x, y, noPixels, &eq)) + ellipseFitError = ComputeEllipseError(&eq, x, y, noPixels); + } //end-if + + if (circleFitError <= LONG_ARC_ERROR) { + addCircle(circles1, noCircles1, xc, yc, r, circleFitError, x, y, noPixels); + bm->move(noPixels); + continue; + + } + else if (ellipseFitError <= ELLIPSE_ERROR) { + double major, minor; + ComputeEllipseCenterAndAxisLengths(&eq, &xc, &yc, &major, &minor); + + // Assume major is longer. Otherwise, swap + if (minor > major) { double tmp = major; major = minor; minor = tmp; } + + if (major < 8 * minor) { + addCircle(circles1, noCircles1, xc, yc, r, circleFitError, &eq, ellipseFitError, x, y, noPixels); + bm->move(noPixels); + } //end-if + + continue; + } //end-else + } //end-if + } //end-if + + // Otherwise, split to lines + EDLines::SplitSegment2Lines(x, y, noPixels, i, lines); + } + + segmentStartLines[segmentNos] = lines.size(); + + // ------------------------------- DETECT ARCS --------------------------------- + + info = new Info[lines.size()]; + + // Compute the angle information for each line segment + for (int i = 0; iex - l2->sx; + double dy = l1->ey - l2->sy; + double d = sqrt(dx*dx + dy*dy); + if (d >= 15) { info[j].angle = 10; info[j].sign = 2; info[j].taken = false; continue; } +#endif + + // Compute the angle between the lines & their turn direction + double v1x = l1->ex - l1->sx; + double v1y = l1->ey - l1->sy; + double v1Len = sqrt(v1x*v1x + v1y*v1y); + + double v2x = l2->ex - l2->sx; + double v2y = l2->ey - l2->sy; + double v2Len = sqrt(v2x*v2x + v2y*v2y); + + double dotProduct = (v1x*v2x + v1y*v2y) / (v1Len*v2Len); + if (dotProduct > 1.0) dotProduct = 1.0; + else if (dotProduct < -1.0) dotProduct = -1.0; + + info[j].angle = acos(dotProduct); + info[j].sign = (v1x*v2y - v2x*v1y) >= 0 ? 1 : -1; // compute cross product + info[j].taken = false; + } //end-for + } //end-for + + // This is how much space we will allocate for circles buffers + int maxNoOfCircles = lines.size() / 3 + noCircles1 * 2; + + edarcs1 = new EDArcs(maxNoOfCircles); + DetectArcs(lines); // Detect all arcs + + // Try to join arcs that are almost perfectly circular. + // Use the distance between the arc end-points as a metric in choosing in choosing arcs to join + edarcs2 = new EDArcs(maxNoOfCircles); + JoinArcs1(); + + // Try to join arcs that belong to the same segment + edarcs3 = new EDArcs(maxNoOfCircles); + JoinArcs2(); + + // Try to combine arcs that belong to different segments + edarcs4 = new EDArcs(maxNoOfCircles); // The remaining arcs + JoinArcs3(); + + // Finally, go over the arcs & circles, and generate candidate circles + GenerateCandidateCircles(); + + //----------------------------- VALIDATE CIRCLES -------------------------- + noCircles2 = 0; + circles2 = new Circle[maxNoOfCircles]; + GaussianBlur(srcImage, smoothImage, Size(), 0.50); // calculate kernel from sigma; + + ValidateCircles(); + + //----------------------------- JOIN CIRCLES -------------------------- + noCircles3 = 0; + circles3 = new Circle[maxNoOfCircles]; + JoinCircles(); + + noCircles = 0; + noEllipses = 0; + for (int i = 0; i < noCircles3; i++) { + if (circles3[i].isEllipse) + noEllipses++; + else + noCircles++; + } + + for (int i = 0; i lines; + + +#define CIRCLE_MIN_LINE_LEN 6 + + for (int i = 0; i < segmentNos; i++) { + + // Make note of the starting line number for this segment + segmentStartLines[i] = lines.size(); + + int noPixels = segmentPoints[i].size(); + + if (noPixels < 2 * CIRCLE_MIN_LINE_LEN) + continue; + + double *x = bm->getX(); + double *y = bm->getY(); + + for (int j = 0; j= 4 * CIRCLE_MIN_LINE_LEN) { + // If the end-points of the segment is close to each other, then assume a circular/elliptic structure + double dx = x[0] - x[noPixels - 1]; + double dy = y[0] - y[noPixels - 1]; + double d = sqrt(dx*dx + dy*dy); + double r = noPixels / TWOPI; // Assume a complete circle + + double maxDistanceBetweenEndPoints = MAX(3, r / 4); + + // If almost closed loop, then try to fit a circle/ellipse + if (d <= maxDistanceBetweenEndPoints) { + double xc, yc, r, circleFitError = 1e10; + + CircleFit(x, y, noPixels, &xc, &yc, &r, &circleFitError); + + EllipseEquation eq; + double ellipseFitError = 1e10; + + if (circleFitError > LONG_ARC_ERROR) { + // Try fitting an ellipse + if (EllipseFit(x, y, noPixels, &eq)) + ellipseFitError = ComputeEllipseError(&eq, x, y, noPixels); + } //end-if + + if (circleFitError <= LONG_ARC_ERROR) { + addCircle(circles1, noCircles1, xc, yc, r, circleFitError, x, y, noPixels); + bm->move(noPixels); + continue; + + } + else if (ellipseFitError <= ELLIPSE_ERROR) { + double major, minor; + ComputeEllipseCenterAndAxisLengths(&eq, &xc, &yc, &major, &minor); + + // Assume major is longer. Otherwise, swap + if (minor > major) { double tmp = major; major = minor; minor = tmp; } + + if (major < 8 * minor) { + addCircle(circles1, noCircles1, xc, yc, r, circleFitError, &eq, ellipseFitError, x, y, noPixels); + bm->move(noPixels); + } //end-if + + continue; + } //end-else + } //end-if + } //end-if + + // Otherwise, split to lines + EDLines::SplitSegment2Lines(x, y, noPixels, i, lines); + } + + segmentStartLines[segmentNos] = lines.size(); + + // ------------------------------- DETECT ARCS --------------------------------- + + info = new Info[lines.size()]; + + // Compute the angle information for each line segment + for (int i = 0; iex - l2->sx; + double dy = l1->ey - l2->sy; + double d = sqrt(dx*dx + dy*dy); + if (d >= 15) { info[j].angle = 10; info[j].sign = 2; info[j].taken = false; continue; } +#endif + + // Compute the angle between the lines & their turn direction + double v1x = l1->ex - l1->sx; + double v1y = l1->ey - l1->sy; + double v1Len = sqrt(v1x*v1x + v1y*v1y); + + double v2x = l2->ex - l2->sx; + double v2y = l2->ey - l2->sy; + double v2Len = sqrt(v2x*v2x + v2y*v2y); + + double dotProduct = (v1x*v2x + v1y*v2y) / (v1Len*v2Len); + if (dotProduct > 1.0) dotProduct = 1.0; + else if (dotProduct < -1.0) dotProduct = -1.0; + + info[j].angle = acos(dotProduct); + info[j].sign = (v1x*v2y - v2x*v1y) >= 0 ? 1 : -1; // compute cross product + info[j].taken = false; + } //end-for + } //end-for + + // This is how much space we will allocate for circles buffers + int maxNoOfCircles = lines.size() / 3 + noCircles1 * 2; + + edarcs1 = new EDArcs(maxNoOfCircles); + DetectArcs(lines); // Detect all arcs + + // Try to join arcs that are almost perfectly circular. + // Use the distance between the arc end-points as a metric in choosing in choosing arcs to join + edarcs2 = new EDArcs(maxNoOfCircles); + JoinArcs1(); + + // Try to join arcs that belong to the same segment + edarcs3 = new EDArcs(maxNoOfCircles); + JoinArcs2(); + + // Try to combine arcs that belong to different segments + edarcs4 = new EDArcs(maxNoOfCircles); // The remaining arcs + JoinArcs3(); + + // Finally, go over the arcs & circles, and generate candidate circles + GenerateCandidateCircles(); + + // EDCircles does not use validation when constructed wia EDColor object. + // TODO :: apply validation to color image + + //----------------------------- VALIDATE CIRCLES -------------------------- + //noCircles2 = 0; + //circles2 = new Circle[maxNoOfCircles]; + //GaussianBlur(srcImage, smoothImage, Size(), 0.50); // calculate kernel from sigma; + + //ValidateCircles(); + + //----------------------------- JOIN CIRCLES -------------------------- + //noCircles3 = 0; + //circles3 = new Circle[maxNoOfCircles]; + //JoinCircles(); + + + noCircles = 0; + noEllipses = 0; + for (int i = 0; i < noCircles1; i++) { + if (circles1[i].isEllipse) + noEllipses++; + else + noCircles++; + } + + for (int i = 0; i EDCircles::getCircles() +{ + return circles; +} + +vector EDCircles::getEllipses() +{ + return ellipses; +} + +int EDCircles::getCirclesNo() +{ + return noCircles; +} + +int EDCircles::getEllipsesNo() +{ + return noEllipses; +} + +void EDCircles::GenerateCandidateCircles() +{ + // Now, go over the circular arcs & add them to circles1 + MyArc *arcs = edarcs4->arcs; + for (int i = 0; inoArcs; i++) { + if (arcs[i].isEllipse) { + // Ellipse + if (arcs[i].coverRatio >= CANDIDATE_ELLIPSE_RATIO && arcs[i].ellipseFitError <= ELLIPSE_ERROR) { + addCircle(circles1, noCircles1, arcs[i].xc, arcs[i].yc, arcs[i].r, arcs[i].circleFitError, &arcs[i].eq, arcs[i].ellipseFitError, + arcs[i].x, arcs[i].y, arcs[i].noPixels); + + } + else { + // double coverRatio = arcs[i].coverRatio; + double coverRatio = MAX(ArcLength(arcs[i].sTheta, arcs[i].eTheta) / TWOPI, arcs[i].coverRatio); + if ((coverRatio >= FULL_CIRCLE_RATIO && arcs[i].circleFitError <= LONG_ARC_ERROR) || + (coverRatio >= HALF_CIRCLE_RATIO && arcs[i].circleFitError <= HALF_ARC_ERROR) || + (coverRatio >= CANDIDATE_CIRCLE_RATIO2 && arcs[i].circleFitError <= SHORT_ARC_ERROR)) { + addCircle(circles1, noCircles1, arcs[i].xc, arcs[i].yc, arcs[i].r, arcs[i].circleFitError, arcs[i].x, arcs[i].y, arcs[i].noPixels); + } //end-if + } //end-else + + } + else { + // If a very short arc, ignore + if (arcs[i].coverRatio < CANDIDATE_CIRCLE_RATIO1) continue; + + // If the arc is long enough and the circleFitError is small enough, assume a circle + if ((arcs[i].coverRatio >= FULL_CIRCLE_RATIO && arcs[i].circleFitError <= LONG_ARC_ERROR) || + (arcs[i].coverRatio >= HALF_CIRCLE_RATIO && arcs[i].circleFitError <= HALF_ARC_ERROR) || + (arcs[i].coverRatio >= CANDIDATE_CIRCLE_RATIO2 && arcs[i].circleFitError <= SHORT_ARC_ERROR)){ + + addCircle(circles1, noCircles1, arcs[i].xc, arcs[i].yc, arcs[i].r, arcs[i].circleFitError, arcs[i].x, arcs[i].y, arcs[i].noPixels); + + continue; + } //end-if + + if (arcs[i].coverRatio < CANDIDATE_CIRCLE_RATIO2) continue; + + // Circle is not possible. Try an ellipse + EllipseEquation eq; + double ellipseFitError = 1e10; + double coverRatio; + + int noPixels = arcs[i].noPixels; + if (EllipseFit(arcs[i].x, arcs[i].y, noPixels, &eq)) { + ellipseFitError = ComputeEllipseError(&eq, arcs[i].x, arcs[i].y, noPixels); + coverRatio = noPixels / computeEllipsePerimeter(&eq); + } //end-if + + if (arcs[i].coverRatio > coverRatio) coverRatio = arcs[i].coverRatio; + + if (coverRatio >= CANDIDATE_ELLIPSE_RATIO && ellipseFitError <= ELLIPSE_ERROR) { + addCircle(circles1, noCircles1, arcs[i].xc, arcs[i].yc, arcs[i].r, arcs[i].circleFitError, &eq, ellipseFitError, arcs[i].x, arcs[i].y, arcs[i].noPixels); + } //end-if + } //end-else + } //end-for +} + +void EDCircles::DetectArcs(vector lines) +{ + + double maxLineLengthThreshold = MAX(width, height) / 5; + + double MIN_ANGLE = PI / 30; // 6 degrees + double MAX_ANGLE = PI / 3; // 60 degrees + + for (int iter = 1; iter <= 2; iter++) { + if (iter == 2) MAX_ANGLE = PI / 1.9; // 95 degrees + // if (iter == 2) MAX_ANGLE = PI/2.25; // 80 degrees + + for (int curSegmentNo = 0; curSegmentNo= maxLineLengthThreshold) { firstLine++; continue; } + + // Skip lines that cannot be part of an arc + if (info[firstLine].angle < MIN_ANGLE || info[firstLine].angle > MAX_ANGLE) { firstLine++; continue; } + + // Find a group of lines (at least 3) with the same sign & angle < MAX_ANGLE degrees + int lastLine = firstLine + 1; + while (lastLine < stopLine - 1) { + if (info[lastLine].taken) break; + if (info[lastLine].sign != info[firstLine].sign) break; + + if (lines[lastLine].len >= maxLineLengthThreshold) break; // very long lines cannot be part of an arc + if (info[lastLine].angle < MIN_ANGLE) break; + if (info[lastLine].angle > MAX_ANGLE) break; + + lastLine++; + } //end-while + + bool specialCase = false; + int wrapCase = -1; // 1: wrap the first two lines with the last line, 2: wrap the last two lines with the first line +#if 0 + // If we do not have 3 lines, then continue + if (lastLine - firstLine <= 1) { firstLine = lastLine; continue; } +#else + // if (lastLine-firstLine == 0){firstLine=lastLine; continue;} + if (lastLine - firstLine == 1) { + // Just 2 lines. If long enough, then try to combine. Angle between 15 & 45 degrees. Min. length = 40 + int totalLineLength = lines[firstLine].len + lines[firstLine + 1].len; + int shorterLen = lines[firstLine].len; + int longerLen = lines[firstLine + 1].len; + + if (lines[firstLine + 1].len < shorterLen) { + shorterLen = lines[firstLine + 1].len; + longerLen = lines[firstLine].len; + } //end-if + + if (info[firstLine].angle >= PI / 12 && info[firstLine].angle <= PI / 4 && totalLineLength >= 40 && shorterLen * 2 >= longerLen) { + specialCase = true; + } //end-if + + // If the two lines do not make up for arc generation, then try to wrap the lines to the first OR last line. + // There are two wrapper cases: + if (specialCase == false) { + // Case 1: Combine the first two lines with the last line of the segment + if (firstLine == segmentStartLines[curSegmentNo] && info[stopLine - 1].angle >= MIN_ANGLE && info[stopLine - 1].angle <= MAX_ANGLE) { + wrapCase = 1; + specialCase = true; + } //end-if + + // Case 2: Combine the last two lines with the first line of the segment + else if (lastLine == stopLine - 1 && info[lastLine].angle >= MIN_ANGLE && info[lastLine].angle <= MAX_ANGLE) { + wrapCase = 2; + specialCase = true; + } //end-if + } // end-if + + // If still not enough for arc generation, then skip + if (specialCase == false) { + firstLine = lastLine; + continue; + } //end-else + } //end-if +#endif + + // Copy the pixels of this segment to an array + int noPixels = 0; + double *x = bm->getX(); + double *y = bm->getY(); + + // wrapCase 1: Combine the first two lines with the last line of the segment + if (wrapCase == 1) { + int index = lines[stopLine - 1].firstPixelIndex; + + for (int n = 0; nmove(noPixels); + + // Try to fit a circle to the entire arc of lines + double xc, yc, radius, circleFitError; + CircleFit(x, y, noPixels, &xc, &yc, &radius, &circleFitError); + + double coverage = noPixels / (TWOPI*radius); + + // In the case of the special case, the arc must cover at least 22.5 degrees + if (specialCase && coverage < 1.0 / 16) { info[firstLine].taken = true; firstLine = lastLine; continue; } + + // If only 3 lines, use the SHORT_ARC_ERROR + double MYERROR = SHORT_ARC_ERROR; + if (lastLine - firstLine >= 3) MYERROR = LONG_ARC_ERROR; + if (circleFitError <= MYERROR) { + // Add this to the list of arcs + if (wrapCase == 1) { + x += lines[stopLine - 1].len; + y += lines[stopLine - 1].len; + noPixels -= lines[stopLine - 1].len; + + } + else if (wrapCase == 2) { + noPixels -= lines[segmentStartLines[curSegmentNo]].len; + } //end-else + + if ((coverage >= FULL_CIRCLE_RATIO && circleFitError <= LONG_ARC_ERROR)) { + addCircle(circles1, noCircles1, xc, yc, radius, circleFitError, x, y, noPixels); + + } + else { + double sTheta, eTheta; + ComputeStartAndEndAngles(xc, yc, radius, x, y, noPixels, &sTheta, &eTheta); + + addArc(edarcs1->arcs, edarcs1->noArcs, xc, yc, radius, circleFitError, sTheta, eTheta, info[firstLine].sign, curSegmentNo, + (int)x[0], (int)y[0], (int)x[noPixels - 1], (int)y[noPixels - 1], x, y, noPixels); + } //end-else + + for (int m = firstLine; m= FULL_CIRCLE_RATIO); + if (isAlmostClosedLoop || (iter == 1 && coverage >= 0.25)) { // an arc covering at least 90 degrees + EllipseEquation eq; + double ellipseFitError = 1e10; + + bool valid = EllipseFit(x, y, noPixels, &eq); + if (valid) ellipseFitError = ComputeEllipseError(&eq, x, y, noPixels); + + MYERROR = ELLIPSE_ERROR; + if (isAlmostClosedLoop == false) MYERROR = 0.75; + + if (ellipseFitError <= MYERROR) { + // Add this to the list of arcs + if (wrapCase == 1) { + x += lines[stopLine - 1].len; + y += lines[stopLine - 1].len; + noPixels -= lines[stopLine - 1].len; + + } + else if (wrapCase == 2) { + noPixels -= lines[segmentStartLines[curSegmentNo]].len; + } //end-else + + if (isAlmostClosedLoop) { + addCircle(circles1, noCircles1, xc, yc, radius, circleFitError, &eq, ellipseFitError, x, y, noPixels); // Add an ellipse for validation + + } + else { + double sTheta, eTheta; + ComputeStartAndEndAngles(xc, yc, radius, x, y, noPixels, &sTheta, &eTheta); + + addArc(edarcs1->arcs, edarcs1->noArcs, xc, yc, radius, circleFitError, sTheta, eTheta, info[firstLine].sign, curSegmentNo, &eq, ellipseFitError, + (int)x[0], (int)y[0], (int)x[noPixels - 1], (int)y[noPixels - 1], x, y, noPixels); + } //end-else + + for (int m = firstLine; m LONG_ARC_ERROR) { noPixels = noPixelsSave; break; } // Adding this line made the error big. So, we do not use this line + + // OK. Longer arc + XC = xc; + YC = yc; + R = r; + Error = error; + + info[curLine].taken = true; + curLine++; + } //end-while + + double coverage = noPixels / (TWOPI*radius); + if ((coverage >= FULL_CIRCLE_RATIO && circleFitError <= LONG_ARC_ERROR)) { + addCircle(circles1, noCircles1, XC, YC, R, Error, x, y, noPixels); + + } + else { + // Add this to the list of arcs + double sTheta, eTheta; + ComputeStartAndEndAngles(XC, YC, R, x, y, noPixels, &sTheta, &eTheta); + + addArc(edarcs1->arcs, edarcs1->noArcs, XC, YC, R, Error, sTheta, eTheta, info[firstLine].sign, curSegmentNo, + (int)x[0], (int)y[0], (int)x[noPixels - 1], (int)y[noPixels - 1], x, y, noPixels); + } //end-else + + x += noPixels; + y += noPixels; + + firstLine = curLine; + info[curLine].taken = false; // may reuse the last line? + } //end-while-current-arc-of-lines + + firstLine = lastLine; + } //end-while-entire-edge-segment + } //end-for + } //end-for-iter +} + +//----------------------------------------------------------------- +// Go over all circles & ellipses and validate them +// The idea here is to look at all pixels of a circle/ellipse +// rather than only the pixels of the lines making up the circle/ellipse +// +void EDCircles::ValidateCircles() +{ + double prec = PI / 16; // Alignment precision + double prob = 1.0 / 8; // probability of alignment + + double max = width; if (height > max) max = height; + double min = width; if (height < min) min = height; + + double *px = new double[8 * (width + height)]; + double *py = new double[8 * (width + height)]; + + // logNT & LUT for NFA computation + double logNT = 2 * log10((double)(width*height)) + log10((double)(width + height)); + + int lutSize = (width + height) / 8; + nfa = new NFALUT(lutSize, prob, logNT); // create look up table + + // Validate circles & ellipses + bool validateAgain; + int count = 0; + for (int i = 0; ixc; + double yc = circle->yc; + double radius = circle->r; + + // Skip potential invalid circles (sometimes these kinds of candidates get generated!) + if (radius > MAX(width, height)) { i++; continue; } + + validateAgain = false; + + int noPoints = 0; + + if (circle->isEllipse) { + noPoints = (int)(computeEllipsePerimeter(&circle->eq)); + + if (noPoints % 2) noPoints--; + ComputeEllipsePoints(circle->eq.coeff, px, py, noPoints); + + } + else { + ComputeCirclePoints(xc, yc, radius, px, py, &noPoints); + } //end-else + + int pr = -1; // previous row + int pc = -1; // previous column + + int tr = -100; + int tc = -100; + int tcount = 0; + + int noPeripheryPixels = 0; + int noEdgePixels = 0; + int aligned = 0; + for (int j = 0; j= height - 1) continue; + if (c <= 0 || c >= width - 1) continue; + + pr = r; + pc = c; + + int dr = abs(r - tr); + int dc = abs(c - tc); + if (dr + dc >= 2) { + tr = r; + tc = c; + tcount++; + } //end-if + +#if 1 + // + // See if there is an edge pixel within 1 pixel vicinity + // + if (edgeImg[r*width + c] != 255) { + // y-cy=-x-cx y-cy=x-cx + // \ / + // \ IV. / + // \ / + // \ / + // III. + I. quadrant + // / \ + // / \ + // / II. \ + // / \ + // + // (x, y)-->(x-cx, y-cy) + // + + int x = c; + int y = r; + + int diff1 = (int)(y - yc - x + xc); + int diff2 = (int)(y - yc + x - xc); + + if (diff1 < 0) { + if (diff2 > 0) { + // I. quadrant + c = x - 1; + if (c >= 1 && edgeImg[r*width + c] == 255) goto out; + c = x + 1; + if (c= 2 && edgeImg[r*width + c] == 255) goto out; + c = x + 2; + if (c= 1 && edgeImg[r*width + c] == 255) goto out; + r = y + 1; + if (r= 2 && edgeImg[r*width + c] == 255) goto out; + r = y + 2; + if (r 0) { + // II. quadrant + r = y - 1; + if (r >= 1 && edgeImg[r*width + c] == 255) goto out; + r = y + 1; + if (r= 2 && edgeImg[r*width + c] == 255) goto out; + r = y + 2; + if (r= 1 && edgeImg[r*width + c] == 255) goto out; + c = x + 1; + if (c= 2 && edgeImg[r*width + c] == 255) goto out; + c = x + 2; + if (cedgeImg[r*width+c] = 254; +#endif + + // 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 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; + if (circle->isEllipse) { + // Ellipse + derivX = 2 * circle->eq.A()*c + circle->eq.B()*r + circle->eq.D(); + derivY = circle->eq.B()*c + 2 * circle->eq.C()*r + circle->eq.E(); + + } + else { + // circle + derivX = c - xc; + derivY = r - yc; + } //end-else + + double idealPixelAngle = nfa->myAtan2(derivX, -derivY); + double diff = fabs(pixelAngle - idealPixelAngle); + if (diff <= prec || diff >= PI - prec) aligned++; + } //end-for + + double circumference; + if (circle->isEllipse) circumference = computeEllipsePerimeter(&circle->eq); + else circumference = TWOPI*radius; + + // Validate by NFA + bool isValid = nfa->checkValidationByNFA(noPeripheryPixels, aligned); + + if (isValid) { + circles2[count++] = circles1[i]; + + } + else if (circle->isEllipse == false && circle->coverRatio >= CANDIDATE_ELLIPSE_RATIO) { + // Fit an ellipse to this circle, and try to revalidate + double ellipseFitError = 1e10; + EllipseEquation eq; + + if (EllipseFit(circle->x, circle->y, circle->noPixels, &eq)) { + ellipseFitError = ComputeEllipseError(&eq, circle->x, circle->y, circle->noPixels); + } //end-if + + if (ellipseFitError <= ELLIPSE_ERROR) { + circle->isEllipse = true; + circle->ellipseFitError = ellipseFitError; + circle->eq = eq; + + validateAgain = true; + } //end-if + } //end-else + + if (validateAgain == false) i++; + } //end-for + + noCircles2 = count; + + delete px; + delete py; + delete nfa; +} + +void EDCircles::JoinCircles() +{ + // Sort the circles wrt their radius + sortCircle(circles2, noCircles2); + + int noCircles = noCircles2; + Circle *circles = circles2; + + for (int i = 0; i CENTER_DISTANCE_THRESHOLD) continue; + + double diff1, diff2; + if (circles[j].isEllipse) { + diff1 = fabs(majorAxisLength - circles[j].majorAxisLength); + diff2 = fabs(minorAxisLength - circles[j].minorAxisLength); + + } + else { + diff1 = fabs(majorAxisLength - circles[j].r); + diff2 = fabs(minorAxisLength - circles[j].r); + } //end-else + + if (diff1 > AXIS_LENGTH_DIFF_THRESHOLD) continue; + if (diff2 > AXIS_LENGTH_DIFF_THRESHOLD) continue; + + // Add to candidates + candidateCircles[noCandidateCircles] = j; + noCandidateCircles++; + } //end-for + + // Try to join the current arc with the candidate arc (if there is one) + double XC = circles[i].xc; + double YC = circles[i].yc; + double R = circles[i].r; + + double CircleFitError = circles[i].circleFitError; + bool CircleFitValid = false; + + EllipseEquation Eq; + double EllipseFitError; + bool EllipseFitValid = false; + + if (noCandidateCircles > 0) { + int noPixels = circles[i].noPixels; + double *x = bm->getX(); + double *y = bm->getY(); + memcpy(x, circles[i].x, noPixels * sizeof(double)); + memcpy(y, circles[i].y, noPixels * sizeof(double)); + + for (int j = 0; jarcs, edarcs1->noArcs); + + int noArcs = edarcs1->noArcs; + MyArc *arcs = edarcs1->arcs; + + bool *taken = new bool[noArcs]; + for (int i = 0; iarcs[edarcs2->noArcs++] = arcs[i]; continue; } + + // Current arc + bool CircleEqValid = false; + double XC = arcs[i].xc; + double YC = arcs[i].yc; + double R = arcs[i].r; + double CircleFitError = arcs[i].circleFitError; + int Turn = arcs[i].turn; + int NoPixels = arcs[i].noPixels; + + int SX = arcs[i].sx; + int SY = arcs[i].sy; + int EX = arcs[i].ex; + int EY = arcs[i].ey; + + // Take the pixels making up this arc + int noPixels = arcs[i].noPixels; + + double *x = bm->getX(); + double *y = bm->getY(); + memcpy(x, arcs[i].x, noPixels * sizeof(double)); + memcpy(y, arcs[i].y, noPixels * sizeof(double)); + + angles.clear(); + angles.set(arcs[i].sTheta, arcs[i].eTheta); + + while (1) { + bool extendedArc = false; + + // Find other arcs to join with + noCandidateArcs = 0; + + for (int j = i + 1; j radiusDiffThreshold) continue; + + // If 50% of the current arc overlaps with the existing arc, then ignore this arc + if (angles.overlap(arcs[j].sTheta, arcs[j].eTheta) >= 0.50) continue; + + // Compute the distances + // 1: (SX, SY)-(sx, sy) + double dx = SX - arcs[j].sx; + double dy = SY - arcs[j].sy; + double d = sqrt(dx*dx + dy*dy); + int which = 1; + + // 2: (SX, SY)-(ex, ey) + dx = SX - arcs[j].ex; + dy = SY - arcs[j].ey; + double d2 = sqrt(dx*dx + dy*dy); + if (d2 < d) { d = d2; which = 2; } + + // 3: (EX, EY)-(sx, sy) + dx = EX - arcs[j].sx; + dy = EY - arcs[j].sy; + d2 = sqrt(dx*dx + dy*dy); + if (d2 < d) { d = d2; which = 3; } + + // 4: (EX, EY)-(ex, ey) + dx = EX - arcs[j].ex; + dy = EY - arcs[j].ey; + d2 = sqrt(dx*dx + dy*dy); + if (d2 < d) { d = d2; which = 4; } + + // Endpoints must be very close + double maxDistanceBetweenEndpoints = minR*1.75; //1.5; + if (d > maxDistanceBetweenEndpoints) continue; + + // This is to give precedence to better matching arc + d += diff; + + // They have to turn in the same direction + if (which == 2 || which == 3) { + if (Turn != arcs[j].turn) continue; + } + else { + if (Turn == arcs[j].turn) continue; + } //end-else + + // Add to candidate arcs in sorted order. User insertion sort + int index = noCandidateArcs - 1; + while (index >= 0) { + if (candidateArcs[index].dist < d) break; + + candidateArcs[index + 1] = candidateArcs[index]; + index--; + } //end-while + + // Add the new candidate arc to the candidate list + index++; + candidateArcs[index].arcNo = j; + candidateArcs[index].which = which; + candidateArcs[index].dist = d; + noCandidateArcs++; + } //end-for + + // Try to join the current arc with the candidate arc (if there is one) + if (noCandidateArcs > 0) { + for (int j = 0; j LONG_ARC_ERROR) { + // No match. Continue with the next candidate + noPixels = noPixelsSave; + + } + else { + // Match. Take it + extendedArc = true; + CircleEqValid = true; + XC = xc; + YC = yc; + R = r; + CircleFitError = circleFitError; + NoPixels = noPixels; + + taken[CandidateArcNo] = true; + taken[i] = true; + + angles.set(arcs[CandidateArcNo].sTheta, arcs[CandidateArcNo].eTheta); + + // Update the end points of the new arc + switch (Which) { + // (SX, SY)-(sy, sy) + case 1: SX = EX, SY = EY; EX = arcs[CandidateArcNo].ex; EY = arcs[CandidateArcNo].ey; + if (Turn == 1) Turn = -1; else Turn = 1; // reverse the turn direction + break; + + // (SX, SY)-(ex, ey) + case 2: SX = EX, SY = EY; EX = arcs[CandidateArcNo].sx; EY = arcs[CandidateArcNo].sy; + if (Turn == 1) Turn = -1; else Turn = 1; // reverse the turn direction + break; + + // (EX, EY)-(sx, sy) + case 3: EX = arcs[CandidateArcNo].ex; EY = arcs[CandidateArcNo].ey; break; + + // (EX, EY)-(ex, ey) + case 4: EX = arcs[CandidateArcNo].sx; EY = arcs[CandidateArcNo].sy; break; + } //end-switch + + break; // Do not look at the other candidates + } //end-if + } //end-for + } //end-if + + if (extendedArc == false) break; + } //end-while + + if (CircleEqValid == false) { + // Add to arcs + edarcs2->arcs[edarcs2->noArcs++] = arcs[i]; + + } + else { + // Add the current OR the extended arc to the new arcs + double sTheta, eTheta; + angles.computeStartEndTheta(sTheta, eTheta); + + double coverage = ArcLength(sTheta, eTheta) / TWOPI; + if ((coverage >= FULL_CIRCLE_RATIO && CircleFitError <= LONG_ARC_ERROR)) + addCircle(circles1, noCircles1, XC, YC, R, CircleFitError, x, y, NoPixels); + else + addArc(edarcs2->arcs, edarcs2->noArcs,XC, YC, R, CircleFitError, sTheta, eTheta, Turn, arcs[i].segmentNo, SX, SY, EX, EY, x, y, NoPixels, angles.overlapRatio()); + + bm->move(NoPixels); + } //end-if + } //end-for + + delete taken; + delete candidateArcs; +} + +void EDCircles::JoinArcs2() +{ + AngleSet angles; + + // Sort the arcs with respect to their length so that longer arcs are at the beginning + sortArc(edarcs2->arcs, edarcs2->noArcs); + + int noArcs = edarcs2->noArcs; + MyArc *arcs = edarcs2->arcs; + + bool *taken = new bool[noArcs]; + for (int i = 0; igetX(); + double *y = bm->getY(); + memcpy(x, arcs[i].x, noPixels * sizeof(double)); + memcpy(y, arcs[i].y, noPixels * sizeof(double)); + + angles.clear(); + angles.set(arcs[i].sTheta, arcs[i].eTheta); + + while (1) { + bool extendedArc = false; + + // Find other arcs to join with + noCandidateArcs = 0; + + for (int j = i + 1; j radiusDiffThreshold) continue; + + // If 75% of the current arc overlaps with the existing arc, then ignore this arc + if (angles.overlap(arcs[j].sTheta, arcs[j].eTheta) >= 0.75) continue; + + // Compute the distances + // 1: (SX, SY)-(sx, sy) + double dx = SX - arcs[j].sx; + double dy = SY - arcs[j].sy; + double d = sqrt(dx*dx + dy*dy); + int which = 1; + + // 2: (SX, SY)-(ex, ey) + dx = SX - arcs[j].ex; + dy = SY - arcs[j].ey; + double d2 = sqrt(dx*dx + dy*dy); + if (d2 < d) { d = d2; which = 2; } + + // 3: (EX, EY)-(sx, sy) + dx = EX - arcs[j].sx; + dy = EY - arcs[j].sy; + d2 = sqrt(dx*dx + dy*dy); + if (d2 < d) { d = d2; which = 3; } + + // 4: (EX, EY)-(ex, ey) + dx = EX - arcs[j].ex; + dy = EY - arcs[j].ey; + d2 = sqrt(dx*dx + dy*dy); + if (d2 < d) { d = d2; which = 4; } + + // Endpoints must be very close + double maxDistanceBetweenEndpoints = 5; //10; + if (d > maxDistanceBetweenEndpoints) continue; + + // Add to candidate arcs in sorted order. User insertion sort + int index = noCandidateArcs - 1; + while (index >= 0) { + if (candidateArcs[index].dist < d) break; + + candidateArcs[index + 1] = candidateArcs[index]; + index--; + } //end-while + + // Add the new candidate arc to the candidate list + index++; + candidateArcs[index].arcNo = j; + candidateArcs[index].which = which; + candidateArcs[index].dist = d; + noCandidateArcs++; + } //end-for + + // Try to join the current arc with the candidate arc (if there is one) + if (noCandidateArcs > 0) { + for (int j = 0; j ELLIPSE_ERROR) { + // No match. Continue with the next candidate + noPixels = noPixelsSave; + + } + else { + // Match. Take it + extendedArc = true; + EllipseEqValid = true; + Eq = eq; + EllipseFitError = ellipseFitError; + NoPixels = noPixels; + + taken[CandidateArcNo] = true; + taken[i] = true; + + R = (R + arcs[CandidateArcNo].r) / 2.0; + + angles.set(arcs[CandidateArcNo].sTheta, arcs[CandidateArcNo].eTheta); + + // Update the end points of the new arc + switch (Which) { + // (SX, SY)-(sy, sy) + case 1: SX = EX, SY = EY; EX = arcs[CandidateArcNo].ex; EY = arcs[CandidateArcNo].ey; + if (Turn == 1) Turn = -1; else Turn = 1; // reverse the turn direction + break; + + // (SX, SY)-(ex, ey) + case 2: SX = EX, SY = EY; EX = arcs[CandidateArcNo].sx; EY = arcs[CandidateArcNo].sy; + if (Turn == 1) Turn = -1; else Turn = 1; // reverse the turn direction + break; + + // (EX, EY)-(sx, sy) + case 3: EX = arcs[CandidateArcNo].ex; EY = arcs[CandidateArcNo].ey; break; + + // (EX, EY)-(ex, ey) + case 4: EX = arcs[CandidateArcNo].sx; EY = arcs[CandidateArcNo].sy; break; + } //end-switch + + break; // Do not look at the other candidates + } //end-if + } //end-for + } //end-if + + if (extendedArc == false) break; + } //end-while + + if (EllipseEqValid == false) { + // Add to arcs + edarcs3->arcs[edarcs3->noArcs++] = arcs[i]; + + } + else { + // Add the current OR the extended arc to the new arcs + double sTheta, eTheta; + angles.computeStartEndTheta(sTheta, eTheta); + + double XC, YC, R, CircleFitError; + CircleFit(x, y, NoPixels, &XC, &YC, &R, &CircleFitError); + + double coverage = ArcLength(sTheta, eTheta) / TWOPI; + if ((coverage >= FULL_CIRCLE_RATIO && CircleFitError <= LONG_ARC_ERROR)) + addCircle(circles1, noCircles1, XC, YC, R, CircleFitError, x, y, NoPixels); + else + addArc(edarcs3->arcs, edarcs3->noArcs, XC, YC, R, CircleFitError, sTheta, eTheta, Turn, arcs[i].segmentNo, &Eq, EllipseFitError, SX, SY, EX, EY, x, y, NoPixels, angles.overlapRatio()); + + // Move buffer pointers + bm->move(NoPixels); + } //end-if + } //end-for + + delete taken; + delete candidateArcs; +} + +void EDCircles::JoinArcs3() +{ + AngleSet angles; + + // Sort the arcs with respect to their length so that longer arcs are at the beginning + sortArc(edarcs3->arcs, edarcs3->noArcs); + + int noArcs = edarcs3->noArcs; + MyArc *arcs = edarcs3->arcs; + + bool *taken = new bool[noArcs]; + for (int i = 0; igetX(); + double *y = bm->getY(); + memcpy(x, arcs[i].x, noPixels * sizeof(double)); + memcpy(y, arcs[i].y, noPixels * sizeof(double)); + + angles.clear(); + angles.set(arcs[i].sTheta, arcs[i].eTheta); + + while (1) { + bool extendedArc = false; + + // Find other arcs to join with + noCandidateArcs = 0; + + for (int j = i + 1; j minR) continue; + + // If 50% of the current arc overlaps with the existing arc, then ignore this arc + if (angles.overlap(arcs[j].sTheta, arcs[j].eTheta) >= 0.50) continue; + + // Compute the distances + // 1: (SX, SY)-(sx, sy) + double dx = SX - arcs[j].sx; + double dy = SY - arcs[j].sy; + double d = sqrt(dx*dx + dy*dy); + int which = 1; + + // 2: (SX, SY)-(ex, ey) + dx = SX - arcs[j].ex; + dy = SY - arcs[j].ey; + double d2 = sqrt(dx*dx + dy*dy); + if (d2 < d) { d = d2; which = 2; } + + // 3: (EX, EY)-(sx, sy) + dx = EX - arcs[j].sx; + dy = EY - arcs[j].sy; + d2 = sqrt(dx*dx + dy*dy); + if (d2 < d) { d = d2; which = 3; } + + // 4: (EX, EY)-(ex, ey) + dx = EX - arcs[j].ex; + dy = EY - arcs[j].ey; + d2 = sqrt(dx*dx + dy*dy); + if (d2 < d) { d = d2; which = 4; } + + // Endpoints must be very close + if (diff <= 0.50*minR) { if (d > minR*0.75) continue; } + else if (diff <= 0.75*minR) { if (d > minR*0.50) continue; } + else if (diff <= 1.00*minR) { if (d > minR*0.25) continue; } + else continue; + + // This is to allow more circular arcs a precedence + d += diff; + + // They have to turn in the same direction + if (which == 2 || which == 3) { + if (Turn != arcs[j].turn) continue; + } + else { + if (Turn == arcs[j].turn) continue; + } //end-else + + // Add to candidate arcs in sorted order. User insertion sort + int index = noCandidateArcs - 1; + while (index >= 0) { + if (candidateArcs[index].dist < d) break; + + candidateArcs[index + 1] = candidateArcs[index]; + index--; + } //end-while + + // Add the new candidate arc to the candidate list + index++; + candidateArcs[index].arcNo = j; + candidateArcs[index].which = which; + candidateArcs[index].dist = d; + noCandidateArcs++; + } //end-for + + // Try to join the current arc with the candidate arc (if there is one) + if (noCandidateArcs > 0) { + for (int j = 0; j ELLIPSE_ERROR) { + // No match. Continue with the next candidate + noPixels = noPixelsSave; + + } + else { + // Match. Take it + extendedArc = true; + EllipseEqValid = true; + Eq = eq; + EllipseFitError = ellipseFitError; + NoPixels = noPixels; + + taken[CandidateArcNo] = true; + taken[i] = true; + + R = (R + arcs[CandidateArcNo].r) / 2.0; + + angles.set(arcs[CandidateArcNo].sTheta, arcs[CandidateArcNo].eTheta); + + // Update the end points of the new arc + switch (Which) { + // (SX, SY)-(sy, sy) + case 1: SX = EX, SY = EY; EX = arcs[CandidateArcNo].ex; EY = arcs[CandidateArcNo].ey; + if (Turn == 1) Turn = -1; else Turn = 1; // reverse the turn direction + break; + + // (SX, SY)-(ex, ey) + case 2: SX = EX, SY = EY; EX = arcs[CandidateArcNo].sx; EY = arcs[CandidateArcNo].sy; + if (Turn == 1) Turn = -1; else Turn = 1; // reverse the turn direction + break; + + // (EX, EY)-(sx, sy) + case 3: EX = arcs[CandidateArcNo].ex; EY = arcs[CandidateArcNo].ey; break; + + // (EX, EY)-(ex, ey) + case 4: EX = arcs[CandidateArcNo].sx; EY = arcs[CandidateArcNo].sy; break; + } //end-switch + + break; // Do not look at the other candidates + } //end-if + } //end-for + } //end-if + + if (extendedArc == false) break; + } //end-while + + if (EllipseEqValid == false) { + // Add to arcs + edarcs4->arcs[edarcs4->noArcs++] = arcs[i]; + + } + else { + // Add the current OR the extended arc to the new arcs + double sTheta, eTheta; + angles.computeStartEndTheta(sTheta, eTheta); + + double XC, YC, R, CircleFitError; + CircleFit(x, y, NoPixels, &XC, &YC, &R, &CircleFitError); + + double coverage = ArcLength(sTheta, eTheta) / TWOPI; + if ((coverage >= FULL_CIRCLE_RATIO && CircleFitError <= LONG_ARC_ERROR)) + addCircle(circles1, noCircles1, XC, YC, R, CircleFitError, x, y, NoPixels); + else + addArc(edarcs4->arcs, edarcs4->noArcs, XC, YC, R, CircleFitError, sTheta, eTheta, Turn, arcs[i].segmentNo, &Eq, EllipseFitError, SX, SY, EX, EY, x, y, NoPixels, angles.overlapRatio()); + + bm->move(NoPixels); + } //end-if + } //end-for + + delete taken; + delete candidateArcs; +} + +Circle * EDCircles::addCircle(Circle *circles, int &noCircles,double xc, double yc, double r, double circleFitError, double * x, double * y, int noPixels) +{ + circles[noCircles].xc = xc; + circles[noCircles].yc = yc; + circles[noCircles].r = r; + circles[noCircles].circleFitError = circleFitError; + circles[noCircles].coverRatio = noPixels / (TWOPI*r); + + circles[noCircles].x = x; + circles[noCircles].y = y; + circles[noCircles].noPixels = noPixels; + + circles[noCircles].isEllipse = false; + + noCircles++; + + return &circles[noCircles - 1]; +} + +Circle * EDCircles::addCircle(Circle *circles, int &noCircles,double xc, double yc, double r, double circleFitError, EllipseEquation * pEq, double ellipseFitError, double * x, double * y, int noPixels) +{ + circles[noCircles].xc = xc; + circles[noCircles].yc = yc; + circles[noCircles].r = r; + circles[noCircles].circleFitError = circleFitError; + circles[noCircles].coverRatio = noPixels / computeEllipsePerimeter(pEq); + + circles[noCircles].x = x; + circles[noCircles].y = y; + circles[noCircles].noPixels = noPixels; + + circles[noCircles].eq = *pEq; + circles[noCircles].ellipseFitError = ellipseFitError; + circles[noCircles].isEllipse = true; + + noCircles++; + + return &circles[noCircles - 1]; +} + +void EDCircles::sortCircles(Circle *circles, int noCircles) +{ + for (int i = 0; i circles[max].r) max = j; + } //end-for + + if (max != i) { + Circle t = circles[i]; + circles[i] = circles[max]; + circles[max] = t; + } // end-if + } //end-for +} //end-sort + + + +// --------------------------------------------------------------------------- +// Given an ellipse equation, computes the length of the perimeter of the ellipse +// Calculates the ellipse perimeter wrt the Ramajunan II formula +// +double EDCircles::computeEllipsePerimeter(EllipseEquation * eq) +{ + double mult = 1; + + double A = eq->A()*mult; + double B = eq->B()*mult; + double C = eq->C()*mult; + double D = eq->D()*mult; + double E = eq->E()*mult; + double F = eq->F()*mult; + + double A2, B2, C2, D2, E2, F2, theta; //rotated coefficients + double D3, E3, F3; //ellipse form coefficients + double cX, cY, a, b; //(cX,cY) center, a & b: semimajor & semiminor axes + double h; //h = (a-b)^2 / (a+b)^2 + bool rotation = false; + +#define pi 3.14159265 + + //Normalize coefficients + B /= A; C /= A; D /= A; E /= A; F /= A; A /= A; + + if (B == 0) //Then not need to rotate the axes + { + A2 = A; + B2 = B; + C2 = C; + D2 = D; + E2 = E; + F2 = F; + + } + + else if (B != 0) //Rotate the axes + { + rotation = true; + + //Determine the rotation angle (in radians) + theta = atan(B / (A - C)) / 2; + + //Compute the coefficients wrt the new coordinate system + A2 = 0.5 * (A * (1 + cos(2 * theta) + B * sin(2 * theta) + C * (1 - cos(2 * theta)))); + + B2 = (C - A) * sin(2 * theta) + B * cos(2 * theta); //B2 should turn to be zero? + + C2 = 0.5 * (A * (1 - cos(2 * theta) - B * sin(2 * theta) + C * (1 + cos(2 * theta)))); + + D2 = D * cos(theta) + E * sin(theta); + + E2 = -D * sin(theta) + E * cos(theta); + + F2 = F; + } + + //Transform the conic equation into the ellipse form + D3 = D2 / A2; //normalize x term's coef + //A3 = 1; //A2 / A2 + + E3 = E2 / C2; //normalize y term's coef + //C3 = 1; //C2 / C2 + + cX = -(D3 / 2); //center X + cY = -(E3 / 2); //center Y + + F3 = A2 * pow(cX, 2.0) + C2 * pow(cY, 2.0) - F2; + + //semimajor axis + a = sqrt(F3 / A2); + //semiminor axis + b = sqrt(F3 / C2); + + //Center coordinates have to be re-transformed if rotation is applied! + if (rotation) + { + double tmpX = cX, tmpY = cY; + cX = tmpX * cos(theta) - tmpY * sin(theta); + cY = tmpX * sin(theta) + tmpY * cos(theta); + } + + //Perimeter Computation(s) + h = pow((a - b), 2.0) / pow((a + b), 2.0); + + //Ramajunan I + //double P1 = pi * (a + b) * (3 - sqrt(4 - h)); + ///printf("Perimeter of the ellipse is %.5f (Ramajunan I)\n", P1); + + //Ramajunan II + double P2 = pi * (a + b) * (1 + 3 * h / (10 + sqrt(4 - 3 * h))); + // printf("Perimeter of the ellipse is %.5f (Ramajunan II)\n", P2); + + // High-school formula + // double P3 = 2 * pi * sqrt(0.5 * (a*a + b*b)); + // printf("Perimeter of the ellipse is %.5f (simple formula)\n", P3); + + return P2; +#undef pi +} + +double EDCircles::ComputeEllipseError(EllipseEquation * eq, double * px, double * py, int noPoints) +{ + double error = 0; + + double A = eq->A(); + double B = eq->B(); + double C = eq->C(); + double D = eq->D(); + double E = eq->E(); + double F = eq->F(); + + double xc, yc, major, minor; + ComputeEllipseCenterAndAxisLengths(eq, &xc, &yc, &major, &minor); + + for (int i = 0; i fabs(dy)) { + // The line equation is of the form: y = mx+n + double m = dy / dx; + double n = yc - m*xc; + + // a*x^2 + b*x + c + double a = A + B*m + C*m*m; + double b = B*n + 2 * C*m*n + D + E*m; + double c = C*n*n + E*n + F; + double det = b*b - 4 * a*c; + if (det<0) det = 0; + double x1 = -(b + sqrt(det)) / (2 * a); + double x2 = -(b - sqrt(det)) / (2 * a); + + double y1 = m*x1 + n; + double y2 = m*x2 + n; + + dx = px[i] - x1; + dy = py[i] - y1; + double d1 = dx*dx + dy*dy; + + dx = px[i] - x2; + dy = py[i] - y2; + double d2 = dx*dx + dy*dy; + + if (d1A()*mult; + double B = eq->B()*mult; + double C = eq->C()*mult; + double D = eq->D()*mult; + double E = eq->E()*mult; + double F = eq->F()*mult; + + double A2, B2, C2, D2, E2, F2, theta; //rotated coefficients + double D3, E3, F3; //ellipse form coefficients + double cX, cY, a, b; //(cX,cY) center, a & b: semimajor & semiminor axes + bool rotation = false; + +#define pi 3.14159265 + + //Normalize coefficients + B /= A; C /= A; D /= A; E /= A; F /= A; A /= A; + + if (B == 0) //Then not need to rotate the axes + { + A2 = A; + B2 = B; + C2 = C; + D2 = D; + E2 = E; + F2 = F; + + // printf("No Rotation is applied\n\n"); + } + else if (B != 0) //Rotate the axes + { + rotation = true; + + //Determine the rotation angle (in radians) + theta = atan(B / (A - C)) / 2; + + //Compute the coefficients wrt the new coordinate system + A2 = 0.5 * (A * (1 + cos(2 * theta) + B * sin(2 * theta) + C * (1 - cos(2 * theta)))); + + B2 = (C - A) * sin(2 * theta) + B * cos(2 * theta); //B2 should turn to be zero? + + C2 = 0.5 * (A * (1 - cos(2 * theta) - B * sin(2 * theta) + C * (1 + cos(2 * theta)))); + + D2 = D * cos(theta) + E * sin(theta); + + E2 = -D * sin(theta) + E * cos(theta); + + F2 = F; + + // printf("Rotation in degrees: %.2f\n\n", theta * 180 / pi); + } + + //Transform the conic equation into the ellipse form + D3 = D2 / A2; //normalize x term's coef + //A3 = 1; //A2 / A2 + + E3 = E2 / C2; //normalize y term's coef + //C3 = 1; //C2 / C2 + + cX = -(D3 / 2); //center X + cY = -(E3 / 2); //center Y + + F3 = A2 * pow(cX, 2.0) + C2 * pow(cY, 2.0) - F2; + + //semimajor axis + a = sqrt(F3 / A2); + //semiminor axis + b = sqrt(F3 / C2); + + //Center coordinates have to be re-transformed if rotation is applied! + if (rotation) + { + double tmpX = cX, tmpY = cY; + cX = tmpX * cos(theta) - tmpY * sin(theta); + cY = tmpX * sin(theta) + tmpY * cos(theta); + } + + *pxc = cX; + *pyc = cY; + + *pmajorAxisLength = a; + *pminorAxisLength = b; + + //if (a > b) { + // *pmajorAxisLength = a; + // *pminorAxisLength = b; + //}q + //else { + // *pmajorAxisLength = b; + // *pminorAxisLength = a; + //} //end-else + + return theta; +#undef pi +} + +// --------------------------------------------------------------------------- +// Given an ellipse equation, computes "noPoints" many consecutive points +// on the ellipse periferi. These points can be used to draw the ellipse +// noPoints must be an even number. +// +void EDCircles::ComputeEllipsePoints(double * pvec, double * px, double * py, int noPoints) +{ + if (noPoints % 2) noPoints--; + int npts = noPoints / 2; + + double **u = AllocateMatrix(3, npts + 1); + double **Aiu = AllocateMatrix(3, npts + 1); + double **L = AllocateMatrix(3, npts + 1); + double **B = AllocateMatrix(3, npts + 1); + double **Xpos = AllocateMatrix(3, npts + 1); + double **Xneg = AllocateMatrix(3, npts + 1); + double **ss1 = AllocateMatrix(3, npts + 1); + double **ss2 = AllocateMatrix(3, npts + 1); + double *lambda = new double[npts + 1]; + double **uAiu = AllocateMatrix(3, npts + 1); + double **A = AllocateMatrix(3, 3); + double **Ai = AllocateMatrix(3, 3); + double **Aib = AllocateMatrix(3, 2); + double **b = AllocateMatrix(3, 2); + double **r1 = AllocateMatrix(2, 2); + double Ao, Ax, Ay, Axx, Ayy, Axy; + + double pi = 3.14781; + double theta; + int i; + int j; + double kk; + + memset(lambda, 0, sizeof(double)*(npts + 1)); + + Ao = pvec[6]; + Ax = pvec[4]; + Ay = pvec[5]; + Axx = pvec[1]; + Ayy = pvec[3]; + Axy = pvec[2]; + + A[1][1] = Axx; A[1][2] = Axy / 2; + A[2][1] = Axy / 2; A[2][2] = Ayy; + b[1][1] = Ax; b[2][1] = Ay; + + // Generate normals linspace + for (i = 1, theta = 0.0; i <= npts; i++, theta += (pi / npts)) + { + u[1][i] = cos(theta); + u[2][i] = sin(theta); + } + + inverse(A, Ai, 2); + + AperB(Ai, b, Aib, 2, 2, 2, 1); + A_TperB(b, Aib, r1, 2, 1, 2, 1); + r1[1][1] = r1[1][1] - 4 * Ao; + + AperB(Ai, u, Aiu, 2, 2, 2, npts); + for (i = 1; i <= 2; i++) + for (j = 1; j <= npts; j++) + uAiu[i][j] = u[i][j] * Aiu[i][j]; + + for (j = 1; j <= npts; j++) + { + if ((kk = (r1[1][1] / (uAiu[1][j] + uAiu[2][j]))) >= 0.0) + lambda[j] = sqrt(kk); + else + lambda[j] = -1.0; + } + + // Builds up B and L + for (j = 1; j <= npts; j++) + L[1][j] = L[2][j] = lambda[j]; + for (j = 1; j <= npts; j++) + { + B[1][j] = b[1][1]; + B[2][j] = b[2][1]; + } + + for (j = 1; j <= npts; j++) + { + ss1[1][j] = 0.5 * (L[1][j] * u[1][j] - B[1][j]); + ss1[2][j] = 0.5 * (L[2][j] * u[2][j] - B[2][j]); + ss2[1][j] = 0.5 * (-L[1][j] * u[1][j] - B[1][j]); + ss2[2][j] = 0.5 * (-L[2][j] * u[2][j] - B[2][j]); + } + + AperB(Ai, ss1, Xpos, 2, 2, 2, npts); + AperB(Ai, ss2, Xneg, 2, 2, 2, npts); + + for (j = 1; j <= npts; j++) + { + if (lambda[j] == -1.0) + { + px[j - 1] = -1; + py[j - 1] = -1; + px[j - 1 + npts] = -1; + py[j - 1 + npts] = -1; + } + else + { + px[j - 1] = Xpos[1][j]; + py[j - 1] = Xpos[2][j]; + px[j - 1 + npts] = Xneg[1][j]; + py[j - 1 + npts] = Xneg[2][j]; + } + } + + DeallocateMatrix(u, 3); + DeallocateMatrix(Aiu, 3); + DeallocateMatrix(L, 3); + DeallocateMatrix(B, 3); + DeallocateMatrix(Xpos, 3); + DeallocateMatrix(Xneg, 3); + DeallocateMatrix(ss1, 3); + DeallocateMatrix(ss2, 3); + delete lambda; + DeallocateMatrix(uAiu, 3); + DeallocateMatrix(A, 3); + DeallocateMatrix(Ai, 3); + DeallocateMatrix(Aib, 3); + DeallocateMatrix(b, 3); + DeallocateMatrix(r1, 2); +} + + +// Tries to join the last two arcs if their end-points are very close to each other +// and if they are part of the same segment. This is useful in cases where an arc on a segment +// is broken due to a noisy patch along the arc, and the long arc is broken into two or more arcs. +// This function will join such broken arcs +// +void EDCircles::joinLastTwoArcs(MyArc *arcs, int &noArcs) +{ + if (noArcs < 2) return; + + int prev = noArcs - 2; + int last = noArcs - 1; + + if (arcs[prev].segmentNo != arcs[last].segmentNo) return; + if (arcs[prev].turn != arcs[last].turn) return; + if (arcs[prev].isEllipse || arcs[last].isEllipse) return; + + // The radius difference between the arcs must be very small + double minR = MIN(arcs[prev].r, arcs[last].r); + double radiusDiffThreshold = minR*0.25; + + double diff = fabs(arcs[prev].r - arcs[last].r); + if (diff > radiusDiffThreshold) return; + + // End-point distance + double dx = arcs[prev].ex - arcs[last].sx; + double dy = arcs[prev].ey - arcs[last].sy; + double d = sqrt(dx*dx + dy*dy); + + double endPointDiffThreshold = 10; + if (d > endPointDiffThreshold) return; + + // Try join + int noPixels = arcs[prev].noPixels + arcs[last].noPixels; + + double xc, yc, r, circleFitError; + CircleFit(arcs[prev].x, arcs[prev].y, noPixels, &xc, &yc, &r, &circleFitError); + + if (circleFitError <= LONG_ARC_ERROR) { + arcs[prev].noPixels = noPixels; + arcs[prev].circleFitError = circleFitError; + + arcs[prev].xc = xc; + arcs[prev].yc = yc; + arcs[prev].r = r; + arcs[prev].ex = arcs[last].ex; + arcs[prev].ey = arcs[last].ey; + // arcs[prev].eTheta = arcs[last].eTheta; -- Fails in a very nasty way in a very special case (recall circles9 with cvsmooth(7x7)!) So, do not use. + + AngleSet angles; + angles.set(arcs[prev].sTheta, arcs[prev].eTheta); + angles.set(arcs[last].sTheta, arcs[last].eTheta); + angles.computeStartEndTheta(arcs[prev].sTheta, arcs[prev].eTheta); + + // arcs[prev].coverRatio = noPixels/(TWOPI*r); + arcs[prev].coverRatio = ArcLength(arcs[prev].sTheta, arcs[prev].eTheta) / (TWOPI); + + noArcs--; + } //end-if +} + + +//----------------------------------------------------------------------- +// Add a new arc to arcs +// +void EDCircles::addArc(MyArc *arcs, int &noArcs, double xc, double yc, double r, double circleFitError, double sTheta, double eTheta, int turn, int segmentNo, int sx, int sy, int ex, int ey, double * x, double * y, int noPixels, double overlapRatio) +{ + arcs[noArcs].xc = xc; + arcs[noArcs].yc = yc; + arcs[noArcs].r = r; + arcs[noArcs].circleFitError = circleFitError; + + arcs[noArcs].sTheta = sTheta; + arcs[noArcs].eTheta = eTheta; + arcs[noArcs].coverRatio = ArcLength(sTheta, eTheta) / (TWOPI); + + arcs[noArcs].turn = turn; + + arcs[noArcs].segmentNo = segmentNo; + + arcs[noArcs].isEllipse = false; + + arcs[noArcs].sx = sx; + arcs[noArcs].sy = sy; + arcs[noArcs].ex = ex; + arcs[noArcs].ey = ey; + + arcs[noArcs].x = x; + arcs[noArcs].y = y; + arcs[noArcs].noPixels = noPixels; + + noArcs++; + + // See if you can join the last two arcs + joinLastTwoArcs(arcs, noArcs); +} + +//------------------------------------------------------------------------- +// Add an elliptic arc to the list of arcs +// +void EDCircles::addArc(MyArc * arcs, int & noArcs, double xc, double yc, double r, double circleFitError, double sTheta, double eTheta, int turn, int segmentNo, EllipseEquation * pEq, double ellipseFitError, int sx, int sy, int ex, int ey, double * x, double * y, int noPixels, double overlapRatio) +{ + arcs[noArcs].xc = xc; + arcs[noArcs].yc = yc; + arcs[noArcs].r = r; + arcs[noArcs].circleFitError = circleFitError; + + arcs[noArcs].sTheta = sTheta; + arcs[noArcs].eTheta = eTheta; + arcs[noArcs].coverRatio = (double)((1.0 - overlapRatio)*noPixels) / computeEllipsePerimeter(pEq); + // arcs[noArcs].coverRatio = noPixels/ComputeEllipsePerimeter(pEq); + // arcs[noArcs].coverRatio = ArcLength(sTheta, eTheta)/(TWOPI); + + arcs[noArcs].turn = turn; + + arcs[noArcs].segmentNo = segmentNo; + + arcs[noArcs].isEllipse = true; + arcs[noArcs].eq = *pEq; + arcs[noArcs].ellipseFitError = ellipseFitError; + + arcs[noArcs].sx = sx; + arcs[noArcs].sy = sy; + arcs[noArcs].ex = ex; + arcs[noArcs].ey = ey; + + arcs[noArcs].x = x; + arcs[noArcs].y = y; + arcs[noArcs].noPixels = noPixels; + + noArcs++; +} + +//-------------------------------------------------------------- +// Given a circular arc, computes the start & end angles of the arc in radians +// +void EDCircles::ComputeStartAndEndAngles(double xc, double yc, double r, double * x, double * y, int len, double * psTheta, double * peTheta) +{ + double sx = x[0]; + double sy = y[0]; + double ex = x[len - 1]; + double ey = y[len - 1]; + double mx = x[len / 2]; + double my = y[len / 2]; + + double d = (sx - xc) / r; + if (d > 1.0) d = 1.0; else if (d < -1.0) d = -1.0; + double theta1 = acos(d); + + double sTheta; + if (sx >= xc) { + if (sy >= yc) { + // I. quadrant + sTheta = theta1; + } + else { + // IV. quadrant + sTheta = TWOPI - theta1; + } //end-else + } + else { + if (sy >= yc) { + // II. quadrant + sTheta = theta1; + } + else { + // III. quadrant + sTheta = TWOPI - theta1; + } //end-else + } //end-else + + d = (ex - xc) / r; + if (d > 1.0) d = 1.0; else if (d < -1.0) d = -1.0; + theta1 = acos(d); + + double eTheta; + if (ex >= xc) { + if (ey >= yc) { + // I. quadrant + eTheta = theta1; + } + else { + // IV. quadrant + eTheta = TWOPI - theta1; + } //end-else + } + else { + if (ey >= yc) { + // II. quadrant + eTheta = theta1; + } + else { + // III. quadrant + eTheta = TWOPI - theta1; + } //end-else + } //end-else + + // Determine whether the arc is clockwise (CW) or counter-clockwise (CCW) + double circumference = TWOPI*r; + double ratio = len / circumference; + + if (ratio <= 0.25 || ratio >= 0.75) { + double angle1, angle2; + + if (eTheta > sTheta) { + angle1 = eTheta - sTheta; + angle2 = TWOPI - eTheta + sTheta; + + } + else { + angle1 = sTheta - eTheta; + angle2 = TWOPI - sTheta + eTheta; + } //end-else + + angle1 = angle1 / TWOPI; + angle2 = angle2 / TWOPI; + + double diff1 = fabs(ratio - angle1); + double diff2 = fabs(ratio - angle2); + + if (diff1 < diff2) { + // angle1 is correct + if (eTheta > sTheta) { + ; + } + else { + double tmp = sTheta; + sTheta = eTheta; + eTheta = tmp; + } // end-else + + } + else { + // angle2 is correct + if (eTheta > sTheta) { + double tmp = sTheta; + sTheta = eTheta; + eTheta = tmp; + + } + else { + ; + } // end-else + } //end-else + + } + else { + double v1x = mx - sx; + double v1y = my - sy; + double v2x = ex - mx; + double v2y = ey - my; + + // cross product + double cross = v1x*v2y - v1y*v2x; + if (cross < 0) { + // swap sTheta & eTheta + double tmp = sTheta; + sTheta = eTheta; + eTheta = tmp; + } //end-if + } //end-else + + double diff = fabs(sTheta - eTheta); + if (diff < (TWOPI / 120)) { + sTheta = 0; + eTheta = 6.26; // 359 degrees + } //end-if + + // Round the start & etheta to 0 if very close to 6.28 or 0 + if (sTheta >= 6.26) sTheta = 0; + if (eTheta < 1.0 / TWOPI) eTheta = 6.28; // if less than 1 degrees, then round to 6.28 + + *psTheta = sTheta; + *peTheta = eTheta; +} + +void EDCircles::sortArc(MyArc * arcs, int noArcs) +{ + for (int i = 0; i arcs[max].coverRatio) max = j; + } //end-for + + if (max != i) { + MyArc t = arcs[i]; + arcs[i] = arcs[max]; + arcs[max] = t; + } // end-if + } //end-for +} + + +//--------------------------------------------------------------------- +// Fits a circle to a given set of points. There must be at least 2 points +// The circle equation is of the form: (x-xc)^2 + (y-yc)^2 = r^2 +// Returns true if there is a fit, false in case no circles can be fit +// +bool EDCircles::CircleFit(double * x, double * y, int N, double * pxc, double * pyc, double * pr, double * pe) +{ + *pe = 1e20; + if (N < 3) return false; + + double xAvg = 0; + double yAvg = 0; + + for (int i = 0; i circles[max].r) max = j; + } //end-for + + if (max != i) { + Circle t = circles[i]; + circles[i] = circles[max]; + circles[max] = t; + } // end-if + } //end-for +} + +bool EDCircles::EllipseFit(double * x, double * y, int noPoints, EllipseEquation * pResult, int mode) +{ + double **D = AllocateMatrix(noPoints + 1, 7); + double **S = AllocateMatrix(7, 7); + double **Const = AllocateMatrix(7, 7); + double **temp = AllocateMatrix(7, 7); + double **L = AllocateMatrix(7, 7); + double **C = AllocateMatrix(7, 7); + + double **invL = AllocateMatrix(7, 7); + double *d = new double[7]; + double **V = AllocateMatrix(7, 7); + double **sol = AllocateMatrix(7, 7); + double tx, ty; + int nrot = 0; + + memset(d, 0, sizeof(double) * 7); + + + switch (mode) + { + case (FPF): + //fprintf(stderr, "EllipseFit: FPF mode"); + Const[1][3] = -2; + Const[2][2] = 1; + Const[3][1] = -2; + break; + case (BOOKSTEIN): + //fprintf(stderr, "EllipseFit: BOOKSTEIN mode"); + Const[1][1] = 2; + Const[2][2] = 1; + Const[3][3] = 2; + } //end-switch + + if (noPoints < 6) + return false; + + // Now first fill design matrix + for (int i = 1; i <= noPoints; i++) { + tx = x[i - 1]; + ty = y[i - 1]; + + D[i][1] = tx * tx; + D[i][2] = tx * ty; + D[i][3] = ty * ty; + D[i][4] = tx; + D[i][5] = ty; + D[i][6] = 1.0; + } //end-for + + //pm(Const,"Constraint"); + // Now compute scatter matrix S + A_TperB(D, D, S, noPoints, 6, noPoints, 6); + //pm(S,"Scatter"); + + choldc(S, 6, L); + //pm(L,"Cholesky"); + + inverse(L, invL, 6); + //pm(invL,"inverse"); + + AperB_T(Const, invL, temp, 6, 6, 6, 6); + AperB(invL, temp, C, 6, 6, 6, 6); + //pm(C,"The C matrix"); + + jacobi(C, 6, d, V, nrot); + //pm(V,"The Eigenvectors"); /* OK */ + //pv(d,"The eigevalues"); + + A_TperB(invL, V, sol, 6, 6, 6, 6); + //pm(sol,"The GEV solution unnormalized"); /* SOl */ + + // Now normalize them + for (int j = 1; j <= 6; j++) /* Scan columns */ + { + double mod = 0.0; + for (int i = 1; i <= 6; i++) + mod += sol[i][j] * sol[i][j]; + for (int i = 1; i <= 6; i++) + sol[i][j] /= sqrt(mod); + } + + //pm(sol,"The GEV solution"); /* SOl */ + + double zero = 10e-20; + double minev = 10e+20; + int solind = 0; + int i; + switch (mode) + { + case (BOOKSTEIN): // smallest eigenvalue + for (i = 1; i <= 6; i++) + if (d[i] < minev && fabs(d[i]) > zero) + solind = i; + break; + case (FPF): + for (i = 1; i <= 6; i++) + if (d[i] < 0 && fabs(d[i]) > zero) + solind = i; + } + +#if 0 + fprintf(stderr, "SOLUTIONS: Selected: %d\n", solind); + for (i = 1; i <= 6; i++) { + fprintf(stderr, "d[i]: %e, a: %.7lf, b: %.7lf, c: %.7lf, d: %.7lf, e: %.7lf, f: %.7lf\n", + d[i], sol[1][i], sol[2][i], sol[3][i], sol[4][i], sol[5][i], sol[6][i]); + } //end-for + +#endif + + bool valid = true; + if (solind == 0) valid = false; + + if (valid) { + // Now fetch the right solution + for (int j = 1; j <= 6; j++) { + pResult->coeff[j] = sol[j][solind]; + } //end-for + } //end-if + + DeallocateMatrix(D, noPoints + 1); + DeallocateMatrix(S, 7); + DeallocateMatrix(Const, 7); + DeallocateMatrix(temp, 7); + DeallocateMatrix(L, 7); + DeallocateMatrix(C, 7); + DeallocateMatrix(invL, 7); + delete d; + DeallocateMatrix(V, 7); + DeallocateMatrix(sol, 7); + + if (valid) { + int len = (int)computeEllipsePerimeter(pResult); + if (len <= 0 || len > 50000) valid = false; + } //end-else + + return valid; +} + +double ** EDCircles::AllocateMatrix(int noRows, int noColumns) +{ + double **m = new double *[noRows]; + + for (int i = 0; i= 1; k--) sum -= a[i][k] * a[j][k]; + if (i == j) + { + if (sum <= 0.0) + // printf("\nA is not poitive definite!"); + { + } + else + p[i] = sqrt(sum); + } + else + { + a[j][i] = sum / p[i]; + } + } + } + for (i = 1; i <= n; i++) { + for (j = i; j <= n; j++) { + if (i == j) + l[i][i] = p[i]; + else + { + l[j][i] = a[j][i]; + l[i][j] = 0.0; + } + } //end-for-inner + } // end-for-outer + + delete p; +} + +int EDCircles::inverse(double ** TB, double ** InvB, int N) +{ + int k, i, j, p, q; + double mult; + double D, temp; + double maxpivot; + int npivot; + double **B = AllocateMatrix(N + 1, N + 2); + double **A = AllocateMatrix(N + 1, 2 * N + 2); + double **C = AllocateMatrix(N + 1, N + 1); + double eps = 10e-20; + + for (k = 1; k <= N; k++) + for (j = 1; j <= N; j++) + B[k][j] = TB[k][j]; + + for (k = 1; k <= N; k++) + { + for (j = 1; j <= N + 1; j++) + A[k][j] = B[k][j]; + for (j = N + 2; j <= 2 * N + 1; j++) + A[k][j] = (double)0; + A[k][k - 1 + N + 2] = (double)1; + } + for (k = 1; k <= N; k++) + { + maxpivot = fabs((double)A[k][k]); + npivot = k; + for (i = k; i <= N; i++) + if (maxpivot < fabs((double)A[i][k])) + { + maxpivot = fabs((double)A[i][k]); + npivot = i; + } + if (maxpivot >= eps) + { + if (npivot != k) + for (j = k; j <= 2 * N + 1; j++) + { + temp = A[npivot][j]; + A[npivot][j] = A[k][j]; + A[k][j] = temp; + }; + D = A[k][k]; + for (j = 2 * N + 1; j >= k; j--) + A[k][j] = A[k][j] / D; + for (i = 1; i <= N; i++) + { + if (i != k) + { + mult = A[i][k]; + for (j = 2 * N + 1; j >= k; j--) + A[i][j] = A[i][j] - mult * A[k][j]; + } + } + } + else + { // printf("\n The matrix may be singular !!") ; + + DeallocateMatrix(B, N + 1); + DeallocateMatrix(A, N + 1); + DeallocateMatrix(C, N + 1); + + return (-1); + } //end-else + } + /** Copia il risultato nella matrice InvB ***/ + for (k = 1, p = 1; k <= N; k++, p++) + for (j = N + 2, q = 1; j <= 2 * N + 1; j++, q++) + InvB[p][q] = A[k][j]; + + DeallocateMatrix(B, N + 1); + DeallocateMatrix(A, N + 1); + DeallocateMatrix(C, N + 1); + + return (0); +} + +void EDCircles::DeallocateMatrix(double ** m, int noRows) +{ + for (int i = 0; i 4 && fabs(d[ip]) + g == fabs(d[ip]) + // && fabs(d[iq]) + g == fabs(d[iq])) + + if (i > 4 && g == 0.0) + a[ip][iq] = 0.0; + else if (fabs(a[ip][iq]) > tresh) + { + h = d[iq] - d[ip]; + if (g == 0.0) + t = (a[ip][iq]) / h; + else + { + theta = 0.5 * h / (a[ip][iq]); + t = 1.0 / (fabs(theta) + sqrt(1.0 + theta * theta)); + if (theta < 0.0) t = -t; + } + c = 1.0 / sqrt(1 + t * t); + s = t * c; + tau = s / (1.0 + c); + h = t * a[ip][iq]; + z[ip] -= h; + z[iq] += h; + d[ip] -= h; + d[iq] += h; + a[ip][iq] = 0.0; + for (j = 1; j <= ip - 1; j++) + { + ROTATE(a, j, ip, j, iq, tau, s); + } + for (j = ip + 1; j <= iq - 1; j++) + { + ROTATE(a, ip, j, j, iq, tau, s); + } + for (j = iq + 1; j <= n; j++) + { + ROTATE(a, ip, j, iq, j, tau, s); + } + for (j = 1; j <= n; j++) + { + ROTATE(v, j, ip, j, iq, tau, s); + } + ++nrot; + } + } + } + for (ip = 1; ip <= n; ip++) + { + b[ip] += z[ip]; + d[ip] = b[ip]; + z[ip] = 0.0; + } + } + //printf("Too many iterations in routine JACOBI"); + delete b; + delete z; +} + +void EDCircles::ROTATE(double ** a, int i, int j, int k, int l, double tau, double s) +{ + double g, h; + g = a[i][j]; h = a[k][l]; a[i][j] = g - s * (h + g * tau); + a[k][l] = h + s * (g - h * tau); +} + +void AngleSet::_set(double sTheta, double eTheta) { + int arc = next++; + + angles[arc].sTheta = sTheta; + angles[arc].eTheta = eTheta; + angles[arc].next = -1; + + // Add the current arc to the linked list + int prev = -1; + int current = head; + while (1) { + // Empty list? + if (head < 0) { head = arc; break; } + + // End of the list. Add to the end + if (current < 0) { + angles[prev].next = arc; + break; + } //end-if + + if (angles[arc].eTheta <= angles[current].sTheta) { + // Add before current + if (prev < 0) { + angles[arc].next = current; + head = arc; + + } + else { + angles[arc].next = current; + angles[prev].next = arc; + } //end-else + + break; + + } + else if (angles[arc].sTheta >= angles[current].eTheta) { + // continue + prev = current; + current = angles[current].next; + + // End of the list? + if (current < 0) { + angles[prev].next = arc; + break; + } //end-if + + } + else { + // overlaps with current. Join + // First delete current from the list + if (prev < 0) head = angles[head].next; + else angles[prev].next = angles[current].next; + + // Update overlap amount. + if (angles[arc].eTheta < angles[current].eTheta) { + overlapAmount += angles[arc].eTheta - angles[current].sTheta; + } + else { + overlapAmount += angles[current].eTheta - angles[arc].sTheta; + } //end-else + + // Now join current with arc + if (angles[current].sTheta < angles[arc].sTheta) angles[arc].sTheta = angles[current].sTheta; + if (angles[current].eTheta > angles[arc].eTheta) angles[arc].eTheta = angles[current].eTheta; + current = angles[current].next; + } //end-else + } //end-while +} + +void AngleSet::set(double sTheta, double eTheta) +{ + if (eTheta > sTheta) { + _set(sTheta, eTheta); + + } + else { + _set(sTheta, TWOPI); + _set(0, eTheta); + } //end-else +} + +double AngleSet::_overlap(double sTheta, double eTheta) +{ + double o = 0; + + int current = head; + while (current >= 0) { + if (sTheta > angles[current].eTheta) { current = angles[current].next; continue; } + else if (eTheta < angles[current].sTheta) break; + + // 3 cases. + if (sTheta < angles[current].sTheta && eTheta > angles[current].eTheta) { + o += angles[current].eTheta - angles[current].sTheta; + + } + else if (sTheta < angles[current].sTheta) { + o += eTheta - angles[current].sTheta; + + } + else { + o += angles[current].eTheta - sTheta; + } //end-else + + current = angles[current].next; + } //end-while + + return o; +} + +double AngleSet::overlap(double sTheta, double eTheta) +{ + double o; + + if (eTheta > sTheta) { + o = _overlap(sTheta, eTheta); + + } + else { + o = _overlap(sTheta, TWOPI); + o += _overlap(0, eTheta); + } //end-else + + return o / ArcLength(sTheta, eTheta); +} + +void AngleSet::computeStartEndTheta(double & sTheta, double & eTheta) +{ + // Special case: Just one arc + if (angles[head].next < 0) { + sTheta = angles[head].sTheta; + eTheta = angles[head].eTheta; + + return; + } //end-if + + // OK. More than one arc. Find the biggest gap + int current = head; + int nextArc = angles[current].next; + + double biggestGapSTheta = angles[current].eTheta; + double biggestGapEtheta = angles[nextArc].sTheta; + double biggestGapLength = biggestGapEtheta - biggestGapSTheta; + + double start, end, len; + while (1) { + current = nextArc; + nextArc = angles[nextArc].next; + if (nextArc < 0) break; + + start = angles[current].eTheta; + end = angles[nextArc].sTheta; + len = end - start; + + if (len > biggestGapLength) { + biggestGapSTheta = start; + biggestGapEtheta = end; + biggestGapLength = len; + } //end-if + + } //end-while + + // Compute the gap between the last arc & the first arc + start = angles[current].eTheta; + end = angles[head].sTheta; + len = TWOPI - start + end; + if (len > biggestGapLength) { + biggestGapSTheta = start; + biggestGapEtheta = end; + } //end-if + + sTheta = biggestGapEtheta; + eTheta = biggestGapSTheta; +} + +double AngleSet::coverRatio() +{ + int current = head; + + double total = 0; + while (current >= 0) { + total += angles[current].eTheta - angles[current].sTheta; + current = angles[current].next; + } //end-while + + return total / (TWOPI); +} + + diff --git a/EDCircles.h b/EDCircles.h new file mode 100644 index 0000000..a1b52e9 --- /dev/null +++ b/EDCircles.h @@ -0,0 +1,303 @@ +/************************************************************************************************************** +* EDCircles source codes. +* Copyright (C) 2016, Cuneyt Akinlar & Cihan Topal +* E-mails of the authors: cuneytakinlar@gmail.com, cihant@anadolu.edu.tr +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. + +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. + +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +* +* By using this library you are implicitly assumed to have accepted all of the above statements, +* and accept to cite the following papers: +* +* [1] C. Akinlar and C. Topal, “EDCircles: A Real-time Circle Detector with a False Detection Control,” +* Pattern Recognition, 46(3), 725-740, 2013. +* +* [2] C. Akinlar and C. Topal, “EDCircles: Realtime Circle Detection by Edge Drawing (ED),” +* IEEE Int'l Conf. on Acoustics, Speech, and Signal Processing (ICASSP)}, March 2012. +**************************************************************************************************************/ + +#ifndef _EDCIRCLES_ +#define _EDCIRCLES_ + +#include "EDPF.h" +#include "EDLines.h" + +#define PI 3.141592653589793238462 +#define TWOPI (2*PI) + +// Circular arc, circle thresholds +#define VERY_SHORT_ARC_ERROR 0.40 // Used for very short arcs (>= CANDIDATE_CIRCLE_RATIO1 && < CANDIDATE_CIRCLE_RATIO2) +#define SHORT_ARC_ERROR 1.00 // Used for short arcs (>= CANDIDATE_CIRCLE_RATIO2 && < HALF_CIRCLE_RATIO) +#define HALF_ARC_ERROR 1.25 // Used for arcs with length (>=HALF_CIRCLE_RATIO && < FULL_CIRCLE_RATIO) +#define LONG_ARC_ERROR 1.50 // Used for long arcs (>= FULL_CIRCLE_RATIO) + +#define CANDIDATE_CIRCLE_RATIO1 0.25 // 25% -- If only 25% of the circle is detected, it may be a candidate for validation +#define CANDIDATE_CIRCLE_RATIO2 0.33 // 33% -- If only 33% of the circle is detected, it may be a candidate for validation +#define HALF_CIRCLE_RATIO 0.50 // 50% -- If 50% of a circle is detected at any point during joins, we immediately make it a candidate +#define FULL_CIRCLE_RATIO 0.67 // 67% -- If 67% of the circle is detected, we assume that it is fully covered + +// 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 BOOKSTEIN 0 // method1 for ellipse fit +#define FPF 1 // method2 for ellipse fit + +enum ImageStyle{NONE=0, CIRCLES, ELLIPSES, BOTH}; + +// Circle equation: (x-xc)^2 + (y-yc)^2 = r^2 +struct mCircle { + cv::Point2d center; + double r; + mCircle(cv::Point2d _center, double _r) { center = _center; r = _r; } +}; + +// Ellipse equation: Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0 +struct mEllipse { + cv::Point2d center; + cv::Size axes; + double theta; + mEllipse(cv::Point2d _center, cv::Size _axes, double _theta) { center = _center; axes = _axes; theta = _theta; } +}; + +//---------------------------------------------------------- +// Ellipse Equation is of the form: +// 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; + } //end-EllipseEquation + + 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]; } +}; + +// ================================ CIRCLES ================================ +struct Circle { + double xc, yc, r; // Center (xc, yc) & radius. + double circleFitError; // circle fit error + double coverRatio; // Percentage of the circle covered by the arcs making up this circle [0-1] + + double *x, *y; // Pointers to buffers containing the pixels making up this circle + int noPixels; // # of pixels making up this circle + + // If this circle is better approximated by an ellipse, we set isEllipse to true & eq contains the ellipse's equation + EllipseEquation eq; + double ellipseFitError; // ellipse fit error + bool isEllipse; + double majorAxisLength; // Length of the major axis + double minorAxisLength; // Length of the minor axis +}; + +// ------------------------------------------- ARCS ---------------------------------------------------- +struct MyArc { + double xc, yc, r; // center x, y and radius + double circleFitError; // Error during circle fit + + double sTheta, eTheta; // Start & end angle in radius + double coverRatio; // Ratio of the pixels covered on the covering circle [0-1] (noPixels/circumference) + + int turn; // Turn direction: 1 or -1 + + int segmentNo; // SegmentNo where this arc belongs + + int sx, sy; // Start (x, y) coordinate + int ex, ey; // End (x, y) coordinate of the arc + + double *x, *y; // Pointer to buffer containing the pixels making up this arc + int noPixels; // # of pixels making up the arc + + bool isEllipse; // Did we fit an ellipse to this arc? + EllipseEquation eq; // If an ellipse, then the ellipse's equation + double ellipseFitError; // Error during ellipse fit +}; + +// =============================== AngleSet ================================== + +//------------------------------------------------------------------------- +// add a circular arc to the list of arcs +// +inline double ArcLength(double sTheta, double eTheta) { + if (eTheta > sTheta) return eTheta - sTheta; + else return TWOPI - sTheta + eTheta; +} // end-ArcLength + +// A fast implementation of the AngleSet class. The slow implementation is really bad. About 10 times slower than this! +struct AngleSetArc { + double sTheta; + double eTheta; + int next; // Next AngleSetArc in the linked list +}; + +struct AngleSet { + AngleSetArc angles[360]; + int head; + int next; // Next AngleSetArc to be allocated + double overlapAmount; // Total overlap of the arcs in angleSet. Computed during set() function + + AngleSet() { clear(); } //end-AngleSet + void clear() { head = -1; next = 0; overlapAmount = 0; } + double overlapRatio() { return overlapAmount / (TWOPI); } + + void _set(double sTheta, double eTheta); + void set(double sTheta, double eTheta); + + double _overlap(double sTheta, double eTheta); + double overlap(double sTheta, double eTheta); + + void computeStartEndTheta(double &sTheta, double &eTheta); + double coverRatio(); +}; + + +struct EDArcs { + MyArc *arcs; + int noArcs; + +public: + EDArcs(int size = 10000) { + arcs = new MyArc[size]; + noArcs = 0; + } //end-EDArcs + + ~EDArcs() { + delete arcs; + } //end-~EDArcs +}; + + +//----------------------------------------------------------------- +// Buffer manager +struct BufferManager { + double *x, *y; + int index; + + BufferManager(int maxSize) { + x = new double[maxSize]; + y = new double[maxSize]; + index = 0; + } //end-BufferManager + + ~BufferManager() { + delete x; + delete y; + } //end-~BufferManager + + double *getX() { return &x[index]; } + double *getY() { return &y[index]; } + void move(int size) { index += size; } +}; + +struct Info { + int sign; // -1 or 1: sign of the cross product + double angle; // angle with the next line (in radians) + bool taken; // Is this line taken during arc detection +}; + + +class EDCircles: public EDPF { +public: + EDCircles(cv::Mat srcImage); + EDCircles(ED obj); + EDCircles(EDColor obj); + + cv::Mat drawResult(bool, ImageStyle); + + std::vector getCircles(); + std::vector getEllipses(); + int getCirclesNo(); + int getEllipsesNo(); + +private: + int noEllipses; + int noCircles; + std::vector circles; + std::vector ellipses; + + Circle *circles1; + Circle *circles2; + Circle *circles3; + int noCircles1; + int noCircles2; + int noCircles3; + + EDArcs *edarcs1; + EDArcs *edarcs2; + EDArcs *edarcs3; + EDArcs *edarcs4; + + int *segmentStartLines; + BufferManager *bm; + Info *info; + NFALUT *nfa; + + void GenerateCandidateCircles(); + void DetectArcs(std::vector lines); + void ValidateCircles(); + void JoinCircles(); + void JoinArcs1(); + void JoinArcs2(); + void JoinArcs3(); + + // circle utility functions + static Circle *addCircle(Circle *circles, int &noCircles,double xc, double yc, double r, double circleFitError, double *x, double *y, int noPixels); + static Circle *addCircle(Circle *circles, int &noCircles,double xc, double yc, double r, double circleFitError, EllipseEquation *pEq, double ellipseFitError, double *x, double *y, int noPixels); + static void sortCircles(Circle *circles, int noCircles); + static bool CircleFit(double *x, double *y, int N, double *pxc, double *pyc, double *pr, double *pe); + static void ComputeCirclePoints(double xc, double yc, double r, double *px, double *py, int *noPoints); + static void sortCircle(Circle *circles, int noCircles); + + // ellipse utility functions + static bool EllipseFit(double *x, double *y, int noPoints, EllipseEquation *pResult, int mode=FPF); + static double **AllocateMatrix(int noRows, int noColumns); + static void A_TperB(double **_A, double **_B, double **_res, int _righA, int _colA, int _righB, int _colB); + static void choldc(double **a, int n, double **l); + static int inverse(double **TB, double **InvB, int N); + static void DeallocateMatrix(double **m, int noRows); + static void AperB_T(double **_A, double **_B, double **_res, int _righA, int _colA, int _righB, int _colB); + static void AperB(double **_A, double **_B, double **_res, int _righA, int _colA, int _righB, int _colB); + static void jacobi(double **a, int n, double d[], double **v, int nrot); + static void ROTATE(double **a, int i, int j, int k, int l, double tau, double s); + static double computeEllipsePerimeter(EllipseEquation *eq); + static double ComputeEllipseError(EllipseEquation *eq, double *px, double *py, int noPoints); + static double ComputeEllipseCenterAndAxisLengths(EllipseEquation *eq, double *pxc, double *pyc, double *pmajorAxisLength, double *pminorAxisLength); + static void ComputeEllipsePoints(double *pvec, double *px, double *py, int noPoints); + + // arc utility functions + static void joinLastTwoArcs(MyArc *arcs, int &noArcs); + static void addArc(MyArc *arcs, int &noArchs,double xc, double yc, double r, double circleFitError, // Circular arc + double sTheta, double eTheta, int turn, int segmentNo, + int sx, int sy, int ex, int ey, + double *x, double *y, int noPixels, double overlapRatio = 0.0); + static void addArc(MyArc *arcs, int &noArchs, double xc, double yc, double r, double circleFitError, // Elliptic arc + double sTheta, double eTheta, int turn, int segmentNo, + EllipseEquation *pEq, double ellipseFitError, + int sx, int sy, int ex, int ey, + double *x, double *y, int noPixels, double overlapRatio = 0.0); + + static void ComputeStartAndEndAngles(double xc, double yc, double r, + double *x, double *y, int len, + double *psTheta, double *peTheta); + + static void sortArc(MyArc *arcs, int noArcs); +}; + +#endif // ! _EDCIRCLES_ + diff --git a/EDColor.cpp b/EDColor.cpp new file mode 100644 index 0000000..99dd58a --- /dev/null +++ b/EDColor.cpp @@ -0,0 +1,624 @@ +#include "EDColor.h" +#include "ED.h" + +using namespace cv; +using namespace std; + +EDColor::EDColor(Mat srcImage, int gradThresh, int anchor_thresh , double sigma, bool validateSegments) +{ + inputImage = srcImage.clone(); + + // check parameters for sanity + if (sigma < 1) sigma = 1; + if (gradThresh < 1) gradThresh = 1; + if (anchor_thresh < 0) anchor_thresh = 0; + + if (validateSegments) { // setup for validation + anchor_thresh = 0; + divForTestSegment = 2.25; + } + + // split channels (OpenCV uses BGR) + Mat bgr[3]; + split(srcImage, bgr); + blueImg = bgr[0].data; + greenImg = bgr[1].data; + redImg = bgr[2].data; + + height = srcImage.rows; + width = srcImage.cols; + + // Allocate space for L*a*b color space + L_Img = new uchar[width*height]; + a_Img = new uchar[width*height]; + b_Img = new uchar[width*height]; + + // Convert RGB2Lab + MyRGB2LabFast(); + + // Allocate space for smooth channels + smooth_L = new uchar[width*height]; + smooth_a = new uchar[width*height]; + smooth_b = new uchar[width*height]; + + // Smooth Channels + smoothChannel(L_Img, smooth_L, sigma); + smoothChannel(a_Img, smooth_a, sigma); + smoothChannel(b_Img, smooth_b, sigma); + + // Allocate space for direction and gradient images + dirImg = new uchar[width*height]; + gradImg = new short[width*height]; + + // Compute Gradient & Edge Direction Maps + ComputeGradientMapByDiZenzo(); + + + // Validate edge segments if the flag is set + if (validateSegments) { + // Get Edge Image using ED + ED edgeObj = ED(gradImg, dirImg, width, height, gradThresh, anchor_thresh, 1, 10, false); + segments = edgeObj.getSegments(); + edgeImage = edgeObj.getEdgeImage(); + + sigma /= 2.5; + smoothChannel(L_Img, smooth_L, sigma); + smoothChannel(a_Img, smooth_a, sigma); + smoothChannel(b_Img, smooth_b, sigma); + + edgeImg = edgeImage.data; // validation steps uses pointer to edgeImage + + validateEdgeSegments(); + + // Extract the new edge segments after validation + extractNewSegments(); + } + + else { + ED edgeObj = ED(gradImg, dirImg, width, height, gradThresh, anchor_thresh); + segments = edgeObj.getSegments(); + edgeImage = edgeObj.getEdgeImage(); + segmentNo = edgeObj.getSegmentNo(); + } + + // Fix 1 pixel errors in the edge map + fixEdgeSegments(segments, 1); + + // clean space + delete[] L_Img; + delete[] a_Img; + delete[] b_Img; + + delete[] smooth_L; + delete[] smooth_a; + delete[] smooth_b; + + delete[] gradImg; + delete[] dirImg; +} + +cv::Mat EDColor::getEdgeImage() +{ + return edgeImage; +} + + +std::vector> EDColor::getSegments() +{ + return segments; +} + +int EDColor::getSegmentNo() +{ + return segmentNo; +} + +int EDColor::getWidth() +{ + return width; +} + +int EDColor::getHeight() +{ + return height; +} + +void EDColor::MyRGB2LabFast() +{ + // Inialize LUTs if necessary + if (!LUT_Initialized) + InitColorEDLib(); + + // First RGB 2 XYZ + double red, green, blue; + double x, y, z; + + // Space for temp. allocation + double *L = new double[width*height]; + double *a = new double[width*height]; + double *b = new double[width*height]; + + for (int i = 0; imax) max = L[i]; + } //end-for + + double scale = 255.0 / (max - min); + for (int i = 0; imax) max = a[i]; + } //end-for + + scale = 255.0 / (max - min); + for (int i = 0; imax) max = b[i]; + } //end-for + + scale = 255.0 / (max - min); + for (int i = 0; i= -3.14159 / 4 && theta <= 3.14159 / 4) + dirImg[i*width + j] = EDGE_VERTICAL; + else + dirImg[i*width + j] = EDGE_HORIZONTAL; + + gradImg[i*width + j] = grad; + if (grad > max) max = grad; + + } + } // end outer for + + // Scale the gradient values to 0-255 + double scale = 255.0 / max; + for (int i = 0; i0; i--) + grads[i - 1] += grads[i]; + + for (int i = 0; i index1) { + int r = segments[i][end].y; + int c = segments[i][end].x; + + if (gradImg[r*width + c] <= minGrad) end--; + else break; + } //end-while + + int start = minGradIndex + 1; + while (start < index2) { + int r = segments[i][start].y; + int c = segments[i][start].x; + + if (gradImg[r*width + c] <= minGrad) start++; + else break; + } //end-while + + testSegment(i, index1, end); + 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 +// +void EDColor::extractNewSegments() +{ + vector< vector > validSegments; + int noSegments = 0; + + for (int i = 0; i < segments.size(); i++) { + int start = 0; + while (start < segments[i].size()) { + + while (start < segments[i].size()) { + int r = segments[i][start].y; + int c = segments[i][start].x; + + if (edgeImg[r*width + c]) break; + start++; + } //end-while + + int end = start + 1; + while (end < segments[i].size()) { + int r = segments[i][end].y; + int c = segments[i][end].x; + + if (edgeImg[r*width + c] == 0) break; + end++; + } //end-while + + int len = end - start; + 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()); + vector subVec(&segments[i][start], &segments[i][end - 1]); + validSegments[noSegments] = subVec; + noSegments++; + } //end-else + + start = end + 1; + } //end-while + } //end-for + + // Update + segments = validSegments; + segmentNo = noSegments; // = validSegments.size() + +} + + +double EDColor::NFA(double prob, int len) +{ + double nfa = np; + for (int i = 0; i EPSILON; i++) + nfa *= prob; + + return nfa; +} + + +//--------------------------------------------------------- +// 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 EDColor::fixEdgeSegments(std::vector> map, int noPixels) +{ + /// First fix one pixel problems: There are four cases + for (int i = 0; i < map.size(); i++) { + int cp = map[i].size() - 2; // Current pixel index + int n2 = 0; // next next pixel index + + while (n2 < 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; + } //end-if + + cp = n2; + n2 += 2; + + } + else if (r2 == r + 2 && c2 == c) { + if (c1 != c) { + map[i][n1].x = c; + } //end-if + + cp = n2; + n2 += 2; + + } + else if (r2 == r && c2 == c - 2) { + if (r1 != r) { + map[i][n1].y = r; + } //end-if + + cp = n2; + n2 += 2; + + } + else if (r2 == r && c2 == c + 2) { + if (r1 != r) { + map[i][n1].y = r; + } //end-if + + cp = n2; + n2 += 2; + + } + else { + cp++; + n2++; + } //end-else + } //end-while + } // end-for +} + +void EDColor::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; + } //end-for + + 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); + } //end-for + + LUT_Initialized = true; +} + +bool EDColor::LUT_Initialized = false; +double EDColor::LUT1[LUT_SIZE + 1] = {0}; +double EDColor::LUT2[LUT_SIZE + 1] = {0}; \ No newline at end of file diff --git a/EDColor.h b/EDColor.h new file mode 100644 index 0000000..aa333b0 --- /dev/null +++ b/EDColor.h @@ -0,0 +1,108 @@ +/************************************************************************************************************** +* Color Edge Drawing (ED) and Color Edge Drawing Parameter Free (EDPF) source codes. +* Copyright (c) 2017, Cuneyt Akinlar & Cihan Topal +* E-mails of the authors: cuneytakinlar@gmail.com, cihant@anadolu.edu.tr +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. + +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. + +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . + +* By using this library you are implicitly assumed to have accepted all of the above statements, +* and accept to cite the following papers: +* +* [1] C. Topal and C. Akinlar, “Edge Drawing: A Combined Real-Time Edge and Segment Detector,” +* Journal of Visual Communication and Image Representation, 23(6), 862-872, doi:10.1016/j.jvcir.2012.05.004 ( +* +* [2] C. Akinlar and C. Topal, “EDPF: A Real-time Parameter-free Edge Segment Detector with a False Detection Con +* International Journal of Pattern Recognition and Artificial Intelligence, 26(1), doi:10.1142/S0218001412550 +* +* [3] C. Akinlar, C. Topal, "ColorED: Color Edge and Segment Detection by Edge Drawing (ED)," +* submitted to the Journal of Visual Communication and Image Representation (2017). +**************************************************************************************************************/ + +#ifndef _EDColor_ +#define _EDColor_ + +#include + +// Look up table size for fast color space conversion +#define LUT_SIZE (1024*4096) + +// Special defines +#define EDGE_VERTICAL 1 +#define EDGE_HORIZONTAL 2 +#define EDGE_45 3 +#define EDGE_135 4 + +#define MAX_GRAD_VALUE 128*256 +#define EPSILON 1.0 +#define MIN_PATH_LEN 10 + + +class EDColor { +public: + EDColor(cv::Mat srcImage, int gradThresh = 20, int anchor_thresh = 4, double sigma = 1.5, bool validateSegments=false); + cv::Mat getEdgeImage(); + std::vector> getSegments(); + int getSegmentNo(); + + int getWidth(); + int getHeight(); + + cv::Mat inputImage; +private: + uchar *L_Img; + uchar *a_Img; + uchar *b_Img; + + uchar *smooth_L; + uchar *smooth_a; + uchar *smooth_b; + + uchar *dirImg; + short *gradImg; + + cv::Mat edgeImage; + uchar *edgeImg; + + const uchar *blueImg; + const uchar *greenImg; + const uchar *redImg; + + int width; + int height; + + double divForTestSegment; + double *H; + int np; + int segmentNo; + + std::vector> segments; + + static double LUT1[LUT_SIZE + 1]; + static double LUT2[LUT_SIZE + 1]; + static bool LUT_Initialized; + + void MyRGB2LabFast(); + void ComputeGradientMapByDiZenzo(); + void smoothChannel(uchar *src, uchar *smooth, double sigma); + void validateEdgeSegments(); + void testSegment(int i, int index1, int index2); + void extractNewSegments(); + double NFA(double prob, int len); + + static void fixEdgeSegments(std::vector> map, int noPixels); + + static void InitColorEDLib(); +}; + +#endif // ! _EDColor_ diff --git a/EDLib.h b/EDLib.h new file mode 100644 index 0000000..4f7b8ad --- /dev/null +++ b/EDLib.h @@ -0,0 +1,10 @@ +#ifndef _EDLib_ +#define _EDLib_ + +#include "ED.h" +#include "EDPF.h" +#include "EDLines.h" +#include "EDCircles.h" +#include "EDColor.h" + +#endif \ No newline at end of file diff --git a/EDLines.cpp b/EDLines.cpp new file mode 100644 index 0000000..0478c62 --- /dev/null +++ b/EDLines.cpp @@ -0,0 +1,1244 @@ +#include "EDLines.h" +#include "EDColor.h" +#include "NFA.h" + +using namespace cv; +using namespace std; + +EDLines::EDLines(Mat srcImage , double _line_error, int _min_line_len, double _max_distance_between_two_lines , double _max_error) + :ED(srcImage, SOBEL_OPERATOR, 36, 8) +{ + min_line_len = _min_line_len; + line_error = _line_error; + max_distance_between_two_lines = _max_distance_between_two_lines; + max_error = _max_error; + + if(min_line_len == -1) // 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! + min_line_len = 9; + + + + // Temporary buffers used during line fitting + double *x = new double[(width+height) * 8]; + double *y = new double[(width+height) * 8]; + + linesNo = 0; + + // Use the whole segment + for (int segmentNumber = 0; segmentNumber < segmentPoints.size(); segmentNumber++) { + int k = 0; + std::vector segment = segmentPoints[segmentNumber]; + for (int k = 0; k < segment.size(); k++) { + x[k] = segment[k].x; + y[k] = segment[k].y; + } + SplitSegment2Lines(x, y, segment.size(), segmentNumber); + } + + /*----------- JOIN COLLINEAR LINES ----------------*/ + JoinCollinearLines(); + + /*----------- VALIDATE LINES ----------------*/ +#define PRECISON_ANGLE 22.5 + prec = (PRECISON_ANGLE / 180)*M_PI; + double prob = 0.125; +#undef PRECISON_ANGLE + + double logNT = 2.0*(log10((double)width) + log10((double)height)); + + int lutSize = (width + height) / 8; + nfa = new NFALUT(lutSize, prob, logNT); // create look up table + + ValidateLineSegments(); + + // Delete redundant space from lines + // Pop them back + int size = lines.size(); + for (int i = 1; i <= size - linesNo; i++) + lines.pop_back(); + + for (int i = 0; i segment = segmentPoints[segmentNumber]; + for (int k = 0; k < segment.size(); k++) { + x[k] = segment[k].x; + y[k] = segment[k].y; + } + SplitSegment2Lines(x, y, segment.size(), segmentNumber); + } + + /*----------- JOIN COLLINEAR LINES ----------------*/ + JoinCollinearLines(); + + /*----------- VALIDATE LINES ----------------*/ +#define PRECISON_ANGLE 22.5 + prec = (PRECISON_ANGLE / 180)*M_PI; + double prob = 0.125; +#undef PRECISON_ANGLE + + double logNT = 2.0*(log10((double)width) + log10((double)height)); + + int lutSize = (width + height) / 8; + nfa = new NFALUT(lutSize, prob, logNT); // create look up table + + ValidateLineSegments(); + + // Delete redundant space from lines + // Pop them back + int size = lines.size(); + for (int i = 1; i <= size - linesNo; i++) + lines.pop_back(); + + + for (int i = 0; i segment = segmentPoints[segmentNumber]; + for (int k = 0; k < segment.size(); k++) { + x[k] = segment[k].x; + y[k] = segment[k].y; + } + SplitSegment2Lines(x, y, segment.size(), segmentNumber); + } + + /*----------- JOIN COLLINEAR LINES ----------------*/ + JoinCollinearLines(); + + /*----------- VALIDATE LINES ----------------*/ +#define PRECISON_ANGLE 22.5 + prec = (PRECISON_ANGLE / 180)*M_PI; + double prob = 0.125; +#undef PRECISON_ANGLE + + double logNT = 2.0*(log10((double)width) + log10((double)height)); + + int lutSize = (width + height) / 8; + nfa = new NFALUT(lutSize, prob, logNT); // create look up table + + // Since edge segments are validated in ed color, + // Validation is not performed again in line segment detection + // TODO :: further validation might be applied. + // ValidateLineSegments(); + + // Delete redundant space from lines + // Pop them back + int size = lines.size(); + for (int i = 1; i <= size - linesNo; i++) + lines.pop_back(); + + + for (int i = 0; i EDLines::getLines() +{ + return linePoints; +} + +int EDLines::getLinesNo() +{ + return linesNo; +} + +Mat EDLines::getLineImage() +{ + Mat lineImage = Mat(height, width, CV_8UC1, Scalar(255)); + for (int i = 0; i < linesNo; i++) { + line(lineImage, linePoints[i].start, linePoints[i].end, Scalar(0), 1, LINE_AA, 0); + } + + return lineImage; +} + +Mat EDLines::drawOnImage() +{ + Mat colorImage = Mat(height, width, CV_8UC1, srcImg); + cvtColor(colorImage, colorImage, COLOR_GRAY2BGR); + for (int i = 0; i < linesNo; i++) { + line(colorImage, linePoints[i].start, linePoints[i].end, Scalar(0, 255, 0), 1, LINE_AA, 0); // draw lines as green on image + } + + return colorImage; +} + +//----------------------------------------------------------------------------------------- +// Computes the minimum line length using the NFA formula given width & height values +int EDLines::ComputeMinLineLength() { + // The reason we are dividing the theoretical minimum line length by 2 is because + // we now test short line segments by a line support region rectangle having width=2. + // This means that within a line support region rectangle for a line segment of length "l" + // there are "2*l" many pixels. Thus, a line segment of length "l" has a chance of getting + // validated by NFA. + + double logNT = 2.0*(log10((double)width) + log10((double)height)); + return (int) round((-logNT / log10(0.125))*0.5); +} //end-ComputeMinLineLength + +//----------------------------------------------------------------- +// Given a full segment of pixels, splits the chain to lines +// This code is used when we use the whole segment of pixels +// +void EDLines::SplitSegment2Lines(double * x, double * y, int noPixels, int segmentNo) +{ + + // First pixel of the line segment within the segment of points + int firstPixelIndex = 0; + + while (noPixels >= min_line_len) { + // Start by fitting a line to MIN_LINE_LEN pixels + bool valid = false; + double lastA, lastB, error; + int lastInvert; + + while (noPixels >= min_line_len) { + LineFit(x, y, min_line_len, lastA, lastB, error, lastInvert); + if (error <= 0.5) { valid = true; break; } + +#if 1 + noPixels -= 1; // Go slowly + x += 1; y += 1; + firstPixelIndex += 1; +#else + noPixels -= 2; // Go faster (for speed) + x += 2; y += 2; + firstPixelIndex += 2; +#endif + } //end-while + + if (valid == false) return; + + // Now try to extend this line + int index = min_line_len; + int len = min_line_len; + + while (index < noPixels) { + int startIndex = index; + int lastGoodIndex = index - 1; + int goodPixelCount = 0; + int badPixelCount = 0; + while (index < noPixels) { + double d = ComputeMinDistance(x[index], y[index], lastA, lastB, lastInvert); + + if (d <= line_error) { + lastGoodIndex = index; + goodPixelCount++; + badPixelCount = 0; + + } + else { + badPixelCount++; + if (badPixelCount >= 5) break; + } //end-if + + index++; + } //end-while + + if (goodPixelCount >= 2) { + len += lastGoodIndex - startIndex + 1; + LineFit(x, y, len, lastA, lastB, lastInvert); // faster LineFit + index = lastGoodIndex + 1; + } // end-if + + if (goodPixelCount < 2 || index >= noPixels) { + // End of a line segment. Compute the end points + double sx, sy, ex, ey; + + int index = 0; + while (ComputeMinDistance(x[index], y[index], lastA, lastB, lastInvert) > line_error) index++; + ComputeClosestPoint(x[index], y[index], lastA, lastB, lastInvert, sx, sy); + int noSkippedPixels = index; + + index = lastGoodIndex; + while (ComputeMinDistance(x[index], y[index], lastA, lastB, lastInvert) > line_error) index--; + ComputeClosestPoint(x[index], y[index], lastA, lastB, lastInvert, ex, ey); + + // Add the line segment to lines + lines.push_back(LineSegment(lastA, lastB, lastInvert, sx, sy, ex, ey, segmentNo, firstPixelIndex + noSkippedPixels, index - noSkippedPixels + 1)); + linesNo++; + len = index + 1; + break; + } //end-else + } //end-while + + noPixels -= len; + x += len; + y += len; + firstPixelIndex += len; + } //end-while +} + +//------------------------------------------------------------------ +// Goes over the original line segments and combines collinear lines that belong to the same segment +// +void EDLines::JoinCollinearLines() +{ + int lastLineIndex = -1; //Index of the last line in the joined lines + int i = 0; + while (i < linesNo) { + int segmentNo = lines[i].segmentNo; + + lastLineIndex++; + if (lastLineIndex != i) + lines[lastLineIndex] = lines[i]; + + int firstLineIndex = lastLineIndex; // Index of the first line in this segment + + int count = 1; + for (int j = i + 1; j< linesNo; j++) { + if (lines[j].segmentNo != segmentNo) break; + + // Try to combine this line with the previous line in this segment + if (TryToJoinTwoLineSegments(&lines[lastLineIndex], &lines[j], + lastLineIndex) == false) { + lastLineIndex++; + if (lastLineIndex != j) + lines[lastLineIndex] = lines[j]; + + } //end-if + + count++; + } //end-for + + // Try to join the first & last line of this segment + if (firstLineIndex != lastLineIndex) { + if (TryToJoinTwoLineSegments(&lines[firstLineIndex], &lines[lastLineIndex], + firstLineIndex)) { + lastLineIndex--; + } //end-if + } //end-if + + i += count; + } //end-while + + linesNo = lastLineIndex + 1; +} + +void EDLines::ValidateLineSegments() +{ + + int *x = new int[(width + height) * 4]; + int *y = new int[(width + height) * 4]; + + int noValidLines = 0; + int eraseOffset = 0; + for (int i = 0; i< linesNo; i++) { + LineSegment *ls = &lines[i]; + + // Compute Line's angle + double lineAngle; + + if (ls->invert == 0) { + // y = a + bx + lineAngle = atan(ls->b); + + } + else { + // x = a + by + lineAngle = atan(1.0 / ls->b); + } //end-else + + if (lineAngle < 0) lineAngle += M_PI; + + Point *pixels = &(segmentPoints[ls->segmentNo][0]); + int noPixels = ls->len; + + bool valid = false; + + // Accept very long lines without testing. They are almost never invalidated. + if (ls->len >= 80) { + valid = true; + + // Validate short line segments by a line support region rectangle having width=2 + } + else if (ls->len <= 25) { + valid = ValidateLineSegmentRect( x, y, ls); + + } + else { + // Longer line segments are first validated by a line support region rectangle having width=1 (for speed) + // If the line segment is still invalid, then a line support region rectangle having width=2 is tried + // If the line segment fails both tests, it is discarded + int aligned = 0; + int count = 0; + for (int j = 0; j= height - 1 || c <= 0 || c >= width - 1) continue; + + count++; + + // compute gx & gy using the simple [-1 -1 -1] + // [ 1 1 1] filter in both directions + // Faster method below + // A B C + // D x E + // F G H + // gx = (C-A) + (E-D) + (H-F) + // gy = (F-A) + (G-B) + (H-C) + // + // To make this faster: + // com1 = (H-A) + // com2 = (C-F) + // Then: gx = com1 + com2 + (E-D) = (H-A) + (C-F) + (E-D) = (C-A) + (E-D) + (H-F) + // gy = com2 - com1 + (G-B) = (H-A) - (C-F) + (G-B) = (F-A) + (G-B) + (H-C) + // + int com1 = srcImg[(r + 1)*width + c + 1] - srcImg[(r - 1)*width + c - 1]; + int com2 = srcImg[(r - 1)*width + c + 1] - srcImg[(r + 1)*width + c - 1]; + + int gx = com1 + com2 + srcImg[r*width + c + 1] - srcImg[r*width + c - 1]; + int gy = com1 - com2 + srcImg[(r + 1)*width + c] - srcImg[(r - 1)*width + c]; + + double pixelAngle = nfa->myAtan2((double)gx, (double)-gy); + double diff = fabs(lineAngle - pixelAngle); + + if (diff <= prec || diff >= M_PI - prec) aligned++; + } //end-for + + // Check validation by NFA computation (fast due to LUT) + valid = nfa->checkValidationByNFA(count, aligned); + if (valid == false) valid = ValidateLineSegmentRect(x, y, ls); + } //end-else + + if (valid) { + if (i != noValidLines) lines[noValidLines] = lines[i]; + noValidLines++; + } + else { + invalidLines.push_back(lines[i]); + } //end-else + } //end-for + + linesNo = noValidLines; + + delete x; + delete y; +} + +bool EDLines::ValidateLineSegmentRect(int * x, int * y, LineSegment * ls) +{ + + // Compute Line's angle + double lineAngle; + + if (ls->invert == 0) { + // y = a + bx + lineAngle = atan(ls->b); + + } + else { + // x = a + by + lineAngle = atan(1.0 / ls->b); + } //end-else + + if (lineAngle < 0) lineAngle += M_PI; + + int noPoints = 0; + + // Enumerate all pixels that fall within the bounding rectangle + EnumerateRectPoints(ls->sx, ls->sy, ls->ex, ls->ey, x, y, &noPoints); + + int count = 0; + int aligned = 0; + + for (int i = 0; i= height - 1 || c <= 0 || c >= width - 1) continue; + + count++; + + // compute gx & gy using the simple [-1 -1 -1] + // [ 1 1 1] filter in both directions + // Faster method below + // A B C + // D x E + // F G H + // gx = (C-A) + (E-D) + (H-F) + // gy = (F-A) + (G-B) + (H-C) + // + // To make this faster: + // com1 = (H-A) + // com2 = (C-F) + // Then: gx = com1 + com2 + (E-D) = (H-A) + (C-F) + (E-D) = (C-A) + (E-D) + (H-F) + // gy = com2 - com1 + (G-B) = (H-A) - (C-F) + (G-B) = (F-A) + (G-B) + (H-C) + // + int com1 = srcImg[(r + 1)*width + c + 1] - srcImg[(r - 1)*width + c - 1]; + int com2 = srcImg[(r - 1)*width + c + 1] - srcImg[(r + 1)*width + c - 1]; + + int gx = com1 + com2 + srcImg[r*width + c + 1] - srcImg[r*width + c - 1]; + int gy = com1 - com2 + srcImg[(r + 1)*width + c] - srcImg[(r - 1)*width + c]; + double pixelAngle = nfa->myAtan2((double)gx, (double)-gy); + + double diff = fabs(lineAngle - pixelAngle); + + if (diff <= prec || diff >= M_PI - prec) aligned++; + } //end-for + + return nfa->checkValidationByNFA(count, aligned); +} + + + +double EDLines::ComputeMinDistance(double x1, double y1, double a, double b, int invert) +{ + double x2, y2; + + if (invert == 0) { + if (b == 0) { + x2 = x1; + y2 = a; + + } + else { + // Let the line passing through (x1, y1) that is perpendicular to a+bx be c+dx + double d = -1.0 / (b); + double c = y1 - d*x1; + + x2 = (a - c) / (d - b); + y2 = a + b*x2; + } //end-else + + } + else { + /// invert = 1 + if (b == 0) { + x2 = a; + y2 = y1; + + } + else { + // Let the line passing through (x1, y1) that is perpendicular to a+by be c+dy + double d = -1.0 / (b); + double c = x1 - d*y1; + + y2 = (a - c) / (d - b); + x2 = a + b*y2; + } //end-else + } //end-else + + return sqrt((x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2)); +} + +//--------------------------------------------------------------------------------- +// Given a point (x1, y1) and a line equation y=a+bx (invert=0) OR x=a+by (invert=1) +// Computes the (x2, y2) on the line that is closest to (x1, y1) +// +void EDLines::ComputeClosestPoint(double x1, double y1, double a, double b, int invert, double &xOut, double &yOut) +{ + double x2, y2; + + if (invert == 0) { + if (b == 0) { + x2 = x1; + y2 = a; + + } + else { + // Let the line passing through (x1, y1) that is perpendicular to a+bx be c+dx + double d = -1.0 / (b); + double c = y1 - d*x1; + + x2 = (a - c) / (d - b); + y2 = a + b*x2; + } //end-else + + } + else { + /// invert = 1 + if (b == 0) { + x2 = a; + y2 = y1; + + } + else { + // Let the line passing through (x1, y1) that is perpendicular to a+by be c+dy + double d = -1.0 / (b); + double c = x1 - d*y1; + + y2 = (a - c) / (d - b); + x2 = a + b*y2; + } //end-else + } //end-else + + xOut = x2; + yOut = y2; +} + +//----------------------------------------------------------------------------------- +// Fits a line of the form y=a+bx (invert == 0) OR x=a+by (invert == 1) +// Assumes that the direction of the line is known by a previous computation +// +void EDLines::LineFit(double * x, double * y, int count, double &a, double &b, int invert) +{ + if (count<2) return; + + double S = count, Sx = 0.0, Sy = 0.0, Sxx = 0.0, Sxy = 0.0; + for (int i = 0; i max_distance_between_two_lines) return false; + + // Compute line lengths. Use the longer one as the ground truth + double dx = ls1->sx - ls1->ex; + double dy = ls1->sy - ls1->ey; + double prevLen = sqrt(dx*dx + dy*dy); + + dx = ls2->sx - ls2->ex; + dy = ls2->sy - ls2->ey; + double nextLen = sqrt(dx*dx + dy*dy); + + // Use the longer line as the ground truth + LineSegment *shorter = ls1; + LineSegment *longer = ls2; + + if (prevLen > nextLen) { shorter = ls2; longer = ls1; } + +#if 0 + // Use 5 points to check for collinearity +#define POINT_COUNT 5 + double decr = 1.0 / (POINT_COUNT - 1); + double alpha = 1.0; + dist = 0.0; + + while (alpha >= 0.0) { + double px = alpha*shorter->sx + (1.0 - alpha)*shorter->ex; + double py = alpha*shorter->sy + (1.0 - alpha)*shorter->ey; + + dist += ComputeMinDistance(px, py, longer->a, longer->b, longer->invert); + + alpha -= decr; + } //end-while + + dist /= POINT_COUNT; + +#undef POINT_COUNT + +#else + // Just use 3 points to check for collinearity + dist = ComputeMinDistance(shorter->sx, shorter->sy, longer->a, longer->b, longer->invert); + dist += ComputeMinDistance((shorter->sx + shorter->ex) / 2.0, (shorter->sy + shorter->ey) / 2.0, longer->a, longer->b, longer->invert); + dist += ComputeMinDistance(shorter->ex, shorter->ey, longer->a, longer->b, longer->invert); + + dist /= 3.0; +#endif + + if (dist > max_error) return false; + +#if 0 + // Update the end points of ls1 + if (which == 0) { // SS + ls1->sx = ls2->ex; + ls1->sy = ls2->ey; + + } + else if (which == 1) { // SE + ls1->sx = ls2->sx; + ls1->sy = ls2->sy; + + } + else if (which == 2) { // ES + ls1->ex = ls2->ex; + ls1->ey = ls2->ey; + + } + else { // EE + ls1->ex = ls2->sx; + ls1->ey = ls2->sy; + } //end-else + +#else + /// 4 cases: 1:(s1, s2), 2:(s1, e2), 3:(e1, s2), 4:(e1, e2) + + /// case 1: (s1, s2) + dx = fabs(ls1->sx - ls2->sx); + dy = fabs(ls1->sy - ls2->sy); + double d = dx + dy; + double max = d; + which = 1; + + /// case 2: (s1, e2) + dx = fabs(ls1->sx - ls2->ex); + dy = fabs(ls1->sy - ls2->ey); + d = dx + dy; + if (d > max) { + max = d; + which = 2; + } //end-if + + /// case 3: (e1, s2) + dx = fabs(ls1->ex - ls2->sx); + dy = fabs(ls1->ey - ls2->sy); + d = dx + dy; + if (d > max) { + max = d; + which = 3; + } //end-if + + /// case 4: (e1, e2) + dx = fabs(ls1->ex - ls2->ex); + dy = fabs(ls1->ey - ls2->ey); + d = dx + dy; + if (d > max) { + max = d; + which = 4; + } //end-if + + if (which == 1) { + // (s1, s2) + ls1->ex = ls2->sx; + ls1->ey = ls2->sy; + + } + else if (which == 2) { + // (s1, e2) + ls1->ex = ls2->ex; + ls1->ey = ls2->ey; + + } + else if (which == 3) { + // (e1, s2) + ls1->sx = ls2->sx; + ls1->sy = ls2->sy; + + } + else { + // (e1, e2) + ls1->sx = ls1->ex; + ls1->sy = ls1->ey; + + ls1->ex = ls2->ex; + ls1->ey = ls2->ey; + } //end-else + +#endif + + // Update the first line's parameters + if (ls1->firstPixelIndex + ls1->len + 5 >= ls2->firstPixelIndex) ls1->len += ls2->len; + else if (ls2->len > ls1->len) { + ls1->firstPixelIndex = ls2->firstPixelIndex; + ls1->len = ls2->len; + } //end-if + + UpdateLineParameters(ls1); + lines[changeIndex] = *ls1; + + return true; +} + +//------------------------------------------------------------------------------- +// Computes the minimum distance between the end points of two lines +// +double EDLines::ComputeMinDistanceBetweenTwoLines(LineSegment * ls1, LineSegment * ls2, int * pwhich) +{ + double dx = ls1->sx - ls2->sx; + double dy = ls1->sy - ls2->sy; + double d = sqrt(dx*dx + dy*dy); + double min = d; + int which = SS; + + dx = ls1->sx - ls2->ex; + dy = ls1->sy - ls2->ey; + d = sqrt(dx*dx + dy*dy); + if (d < min) { min = d; which = SE; } + + dx = ls1->ex - ls2->sx; + dy = ls1->ey - ls2->sy; + d = sqrt(dx*dx + dy*dy); + if (d < min) { min = d; which = ES; } + + dx = ls1->ex - ls2->ex; + dy = ls1->ey - ls2->ey; + d = sqrt(dx*dx + dy*dy); + if (d < min) { min = d; which = EE; } + + if (pwhich) *pwhich = which; + return min; +} + +//----------------------------------------------------------------------------------- +// Uses the two end points (sx, sy)----(ex, ey) of the line segment +// and computes the line that passes through these points (a, b, invert) +// +void EDLines::UpdateLineParameters(LineSegment * ls) +{ + double dx = ls->ex - ls->sx; + double dy = ls->ey - ls->sy; + + if (fabs(dx) >= fabs(dy)) { + /// Line will be of the form y = a + bx + ls->invert = 0; + if (fabs(dy) < 1e-3) { ls->b = 0; ls->a = (ls->sy + ls->ey) / 2; } + else { + ls->b = dy / dx; + ls->a = ls->sy - (ls->b)*ls->sx; + } //end-else + + } + else { + /// Line will be of the form x = a + by + ls->invert = 1; + if (fabs(dx) < 1e-3) { ls->b = 0; ls->a = (ls->sx + ls->ex) / 2; } + else { + ls->b = dx / dy; + ls->a = ls->sx - (ls->b)*ls->sy; + } //end-else + } //end-else +} + +void EDLines::EnumerateRectPoints(double sx, double sy, double ex, double ey, int ptsx[], int ptsy[], int * pNoPoints) +{ + double vxTmp[4], vyTmp[4]; + double vx[4], vy[4]; + int n, offset; + + double x1 = sx; + double y1 = sy; + double x2 = ex; + double y2 = ey; + double width = 2; + + double dx = x2 - x1; + double dy = y2 - y1; + double vLen = sqrt(dx*dx + dy*dy); + + // make unit vector + dx = dx / vLen; + dy = dy / vLen; + + /* build list of rectangle corners ordered + in a circular way around the rectangle */ + vxTmp[0] = x1 - dy * width / 2.0; + vyTmp[0] = y1 + dx * width / 2.0; + vxTmp[1] = x2 - dy * width / 2.0; + vyTmp[1] = y2 + dx * width / 2.0; + vxTmp[2] = x2 + dy * width / 2.0; + vyTmp[2] = y2 - dx * width / 2.0; + vxTmp[3] = x1 + dy * width / 2.0; + vyTmp[3] = y1 - dx * width / 2.0; + + /* compute rotation of index of corners needed so that the first + point has the smaller x. + + if one side is vertical, thus two corners have the same smaller x + value, the one with the largest y value is selected as the first. + */ + if (x1 < x2 && y1 <= y2) offset = 0; + else if (x1 >= x2 && y1 < y2) offset = 1; + else if (x1 > x2 && y1 >= y2) offset = 2; + else offset = 3; + + /* apply rotation of index. */ + for (n = 0; n<4; n++) { + vx[n] = vxTmp[(offset + n) % 4]; + vy[n] = vyTmp[(offset + n) % 4]; + } //end-for + + /* Set a initial condition. + + The values are set to values that will cause 'ri_inc' (that will + be called immediately) to initialize correctly the first 'column' + and compute the limits 'ys' and 'ye'. + + 'y' is set to the integer value of vy[0], the starting corner. + + 'ys' and 'ye' are set to very small values, so 'ri_inc' will + notice that it needs to start a new 'column'. + + The smaller integer coordinate inside of the rectangle is + 'ceil(vx[0])'. The current 'x' value is set to that value minus + one, so 'ri_inc' (that will increase x by one) will advance to + the first 'column'. + */ + int x = (int)ceil(vx[0]) - 1; + int y = (int)ceil(vy[0]); + double ys = -DBL_MAX, ye = -DBL_MAX; + + int noPoints = 0; + while (1) { + /* if not at end of exploration, + increase y value for next pixel in the 'column' */ + y++; + + /* if the end of the current 'column' is reached, + and it is not the end of exploration, + advance to the next 'column' */ + while (y > ye && x <= vx[2]) { + /* increase x, next 'column' */ + x++; + + /* if end of exploration, return */ + if (x > vx[2]) break; + + /* update lower y limit (start) for the new 'column'. + + We need to interpolate the y value that corresponds to the + lower side of the rectangle. The first thing is to decide if + the corresponding side is + + vx[0],vy[0] to vx[3],vy[3] or + vx[3],vy[3] to vx[2],vy[2] + + Then, the side is interpolated for the x value of the + 'column'. But, if the side is vertical (as it could happen if + the rectangle is vertical and we are dealing with the first + or last 'columns') then we pick the lower value of the side + by using 'inter_low'. + */ + if ((double)x < vx[3]) { + /* interpolation */ + if (fabs(vx[0] - vx[3]) <= 0.01) { + if (vy[0]vy[3]) ys = vy[3]; + else ys = vy[0] + (x - vx[0]) * (vy[3] - vy[0]) / (vx[3] - vx[0]); + } + else + ys = vy[0] + (x - vx[0]) * (vy[3] - vy[0]) / (vx[3] - vx[0]); + + } + else { + /* interpolation */ + if (fabs(vx[3] - vx[2]) <= 0.01) { + if (vy[3]vy[2]) ys = vy[2]; + else ys = vy[3] + (x - vx[3]) * (y2 - vy[3]) / (vx[2] - vx[3]); + } + else + ys = vy[3] + (x - vx[3]) * (vy[2] - vy[3]) / (vx[2] - vx[3]); + } //end-else + + /* update upper y limit (end) for the new 'column'. + + We need to interpolate the y value that corresponds to the + upper side of the rectangle. The first thing is to decide if + the corresponding side is + + vx[0],vy[0] to vx[1],vy[1] or + vx[1],vy[1] to vx[2],vy[2] + + Then, the side is interpolated for the x value of the + 'column'. But, if the side is vertical (as it could happen if + the rectangle is vertical and we are dealing with the first + or last 'columns') then we pick the lower value of the side + by using 'inter_low'. + */ + if ((double)x < vx[1]) { + /* interpolation */ + if (fabs(vx[0] - vx[1]) <= 0.01) { + if (vy[0]vy[1]) ye = vy[0]; + else ye = vy[0] + (x - vx[0]) * (vy[1] - vy[0]) / (vx[1] - vx[0]); + } + else + ye = vy[0] + (x - vx[0]) * (vy[1] - vy[0]) / (vx[1] - vx[0]); + + } + else { + /* interpolation */ + if (fabs(vx[1] - vx[2]) <= 0.01) { + if (vy[1]vy[2]) ye = vy[1]; + else ye = vy[1] + (x - vx[1]) * (vy[2] - vy[1]) / (vx[2] - vx[1]); + } + else + ye = vy[1] + (x - vx[1]) * (vy[2] - vy[1]) / (vx[2] - vx[1]); + } //end-else + + /* new y */ + y = (int)ceil(ys); + } //end-while + + // Are we done? + if (x > vx[2]) break; + + ptsx[noPoints] = x; + ptsy[noPoints] = y; + noPoints++; + } //end-while + + *pNoPoints = noPoints; +} + +void EDLines::SplitSegment2Lines(double * x, double * y, int noPixels, int segmentNo, vector &lines, int min_line_len, double line_error) +{ + // First pixel of the line segment within the segment of points + int firstPixelIndex = 0; + + while (noPixels >= min_line_len) { + // Start by fitting a line to MIN_LINE_LEN pixels + bool valid = false; + double lastA, lastB, error; + int lastInvert; + + while (noPixels >= min_line_len) { + LineFit(x, y, min_line_len, lastA, lastB, error, lastInvert); + if (error <= 0.5) { valid = true; break; } + +#if 1 + noPixels -= 1; // Go slowly + x += 1; y += 1; + firstPixelIndex += 1; +#else + noPixels -= 2; // Go faster (for speed) + x += 2; y += 2; + firstPixelIndex += 2; +#endif + } //end-while + + if (valid == false) return; + + // Now try to extend this line + int index = min_line_len; + int len = min_line_len; + + while (index < noPixels) { + int startIndex = index; + int lastGoodIndex = index - 1; + int goodPixelCount = 0; + int badPixelCount = 0; + while (index < noPixels) { + double d = ComputeMinDistance(x[index], y[index], lastA, lastB, lastInvert); + + if (d <= line_error) { + lastGoodIndex = index; + goodPixelCount++; + badPixelCount = 0; + + } + else { + badPixelCount++; + if (badPixelCount >= 5) break; + } //end-if + + index++; + } //end-while + + if (goodPixelCount >= 2) { + len += lastGoodIndex - startIndex + 1; + LineFit(x, y, len, lastA, lastB, lastInvert); // faster LineFit + index = lastGoodIndex + 1; + } // end-if + + if (goodPixelCount < 2 || index >= noPixels) { + // End of a line segment. Compute the end points + double sx, sy, ex, ey; + + int index = 0; + while (ComputeMinDistance(x[index], y[index], lastA, lastB, lastInvert) > line_error) index++; + ComputeClosestPoint(x[index], y[index], lastA, lastB, lastInvert, sx, sy); + int noSkippedPixels = index; + + index = lastGoodIndex; + while (ComputeMinDistance(x[index], y[index], lastA, lastB, lastInvert) > line_error) index--; + ComputeClosestPoint(x[index], y[index], lastA, lastB, lastInvert, ex, ey); + + // Add the line segment to lines + lines.push_back(LineSegment(lastA, lastB, lastInvert, sx, sy, ex, ey, segmentNo, firstPixelIndex + noSkippedPixels, index - noSkippedPixels + 1)); + //linesNo++; + len = index + 1; + break; + } //end-else + } //end-while + + noPixels -= len; + x += len; + y += len; + firstPixelIndex += len; + } //end-while +} diff --git a/EDLines.h b/EDLines.h new file mode 100644 index 0000000..4640809 --- /dev/null +++ b/EDLines.h @@ -0,0 +1,128 @@ +/************************************************************************************************************** +* EDLines source codes. +* Copyright (C) 2016, Cuneyt Akinlar & Cihan Topal +* E-mails of the authors: cuneytakinlar@gmail.com, cihant@anadolu.edu.tr +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. + +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. + +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +* +* By using this library you are implicitly assumed to have accepted all of the above statements, +* and accept to cite the following papers: +* +* [1] C. Akinlar and C. Topal, “EDLines: A Real-time Line Segment Detector with a False Detection Control,” +* Pattern Recognition Letters, 32(13), 1633-1642, DOI: 10.1016/j.patrec.2011.06.001 (2011). +* +* [2] C. Akinlar and C. Topal, “EDLines: Realtime Line Segment Detection by Edge Drawing (ED),” +* IEEE Int’l Conf. on Image Processing (ICIP), Sep. 2011. +**************************************************************************************************************/ + +#ifndef _EDLines_ +#define _EDLines_ + +#include "ED.h" +#include "EDColor.h" +#include "NFA.h" + +#define SS 0 +#define SE 1 +#define ES 2 +#define EE 3 + +// light weight struct for Start & End coordinates of the line segment +struct LS { + cv::Point2d start; + cv::Point2d end; + + LS(cv::Point2d _start, cv::Point2d _end) + { + start = _start; + end = _end; + } +}; + + +struct LineSegment { + double a, b; // y = a + bx (if invert = 0) || x = a + by (if invert = 1) + int invert; + + double sx, sy; // starting x & y coordinates + double ex, ey; // ending x & y coordinates + + int segmentNo; // Edge segment that this line belongs to + int firstPixelIndex; // Index of the first pixel within the segment of pixels + int len; // No of pixels making up the line segment + + LineSegment(double _a, double _b, int _invert, double _sx, double _sy, double _ex, double _ey, int _segmentNo, int _firstPixelIndex, int _len) { + a = _a; + b = _b; + invert = _invert; + sx = _sx; + sy = _sy; + ex = _ex; + ey = _ey; + segmentNo = _segmentNo; + firstPixelIndex = _firstPixelIndex; + len = _len; + } +}; + + +class EDLines : public ED { +public: + EDLines(cv::Mat srcImage, double _line_error = 1.0, int _min_line_len = -1, double _max_distance_between_two_lines = 6.0, double _max_error = 1.3); + EDLines(ED obj, double _line_error = 1.0, int _min_line_len = -1, double _max_distance_between_two_lines = 6.0, double _max_error = 1.3); + EDLines(EDColor obj, double _line_error = 1.0, int _min_line_len = -1, double _max_distance_between_two_lines = 6.0, double _max_error = 1.3); + EDLines(); + + std::vector getLines(); + int getLinesNo(); + cv::Mat getLineImage(); + cv::Mat drawOnImage(); + + // EDCircle uses this one + static void SplitSegment2Lines(double *x, double *y, int noPixels, int segmentNo, std::vector &lines, int min_line_len = 6, double line_error = 1.0); + +private: + std::vector lines; + std::vector invalidLines; + std::vector linePoints; + int linesNo; + int min_line_len; + double line_error; + double max_distance_between_two_lines; + double max_error; + double prec; + NFALUT *nfa; + + + int ComputeMinLineLength(); + void SplitSegment2Lines(double *x, double *y, int noPixels, int segmentNo); + void JoinCollinearLines(); + + void ValidateLineSegments(); + bool ValidateLineSegmentRect(int *x, int *y, LineSegment *ls); + bool TryToJoinTwoLineSegments(LineSegment *ls1, LineSegment *ls2, int changeIndex); + + static double ComputeMinDistance(double x1, double y1, double a, double b, int invert); + static void ComputeClosestPoint(double x1, double y1, double a, double b, int invert, double &xOut, double &yOut); + static void LineFit(double *x, double *y, int count, double &a, double &b, int invert); + static void LineFit(double *x, double *y, int count, double &a, double &b, double &e, int &invert); + static double ComputeMinDistanceBetweenTwoLines(LineSegment *ls1, LineSegment *ls2, int *pwhich); + static void UpdateLineParameters(LineSegment *ls); + static void EnumerateRectPoints(double sx, double sy, double ex, double ey,int ptsx[], int ptsy[], int *pNoPoints); + + // Utility math functions + +}; + +#endif \ No newline at end of file diff --git a/EDPF.cpp b/EDPF.cpp new file mode 100644 index 0000000..705cd4d --- /dev/null +++ b/EDPF.cpp @@ -0,0 +1,244 @@ +#include "EDPF.h" + +using namespace cv; +using namespace std; + +EDPF::EDPF(Mat srcImage) + :ED(srcImage, PREWITT_OPERATOR, 11, 3) +{ + // Validate Edge Segments + sigma /= 2.5; + GaussianBlur(srcImage, smoothImage, Size(), sigma); // calculate kernel from sigma + + validateEdgeSegments(); +} + +EDPF::EDPF(ED obj) + :ED(obj) +{ + // Validate Edge Segments + sigma /= 2.5; + GaussianBlur(srcImage, smoothImage, Size(), sigma); // calculate kernel from sigma + + validateEdgeSegments(); +} + +EDPF::EDPF(EDColor obj) + :ED(obj) +{ +} + +void EDPF::validateEdgeSegments() +{ + divForTestSegment = 2.25; // Some magic number :-) + memset(edgeImg, 0, width*height); // clear edge image + + H = new double[MAX_GRAD_VALUE]; + memset(H, 0, sizeof(double)*MAX_GRAD_VALUE); + + gradImg = ComputePrewitt3x3(); + + // Compute np: # of segment pieces +#if 1 + // Does this underestimate the number of pieces of edge segments? + // What's the correct value? + np = 0; + for (int i = 0; i < segmentNos; i++) { + int len = segmentPoints[i].size(); + np += (len*(len - 1)) / 2; + } //end-for + + // np *= 32; +#elif 0 + // This definitely overestimates the number of pieces of edge segments + int np = 0; + for (int i = 0; i < segmentNos; i++) { + np += segmentPoints[i].size(); + } //end-for + np = (np*(np - 1)) / 2; +#endif + + // Validate segments + for (int i = 0; i< segmentNos; i++) { + TestSegment(i, 0, segmentPoints[i].size() - 1); + } //end-for + + ExtractNewSegments(); + + // clean space + delete[] H; + delete[] gradImg; +} + +short * EDPF::ComputePrewitt3x3() +{ + short *gradImg = new short[width*height]; + memset(gradImg, 0, sizeof(short)*width*height); + + int *grads = new int[MAX_GRAD_VALUE]; + memset(grads, 0, sizeof(int)*MAX_GRAD_VALUE); + + for (int i = 1; i0; i--) + grads[i - 1] += grads[i]; + + for (int i = 0; i < MAX_GRAD_VALUE; i++) + H[i] = (double)grads[i] / ((double)size); + + delete[] grads; + return gradImg; +} + +//---------------------------------------------------------------------------------- +// Resursive validation using half of the pixels as suggested by DMM algorithm +// We take pixels at Nyquist distance, i.e., 2 (as suggested by DMM) +// +void EDPF::TestSegment(int i, int index1, int index2) +{ + + int chainLen = index2 - index1 + 1; + if (chainLen < minPathLen) + return; + + // Test from index1 to index2. If OK, then we are done. Otherwise, split into two and + // recursively test the left & right halves + + // First find the min. gradient along the segment + int minGrad = 1 << 30; + int minGradIndex; + for (int k = index1; k <= index2; k++) { + int r = segmentPoints[i][k].y; + int c = segmentPoints[i][k].x; + if (gradImg[r*width + c] < minGrad) { minGrad = gradImg[r*width + c]; minGradIndex = k; } + } //end-for + + // Compute nfa + double nfa = NFA(H[minGrad], (int)(chainLen / divForTestSegment)); + + if (nfa <= EPSILON) { + for (int k = index1; k <= index2; k++) { + int r = segmentPoints[i][k].y; + int c = segmentPoints[i][k].x; + + edgeImg[r*width + c] = 255; + } //end-for + + return; + } //end-if + + // Split into two halves. We divide at the point where the gradient is the minimum + int end = minGradIndex - 1; + while (end > index1) { + int r = segmentPoints[i][end].y; + int c = segmentPoints[i][end].x; + + if (gradImg[r*width + c] <= minGrad) end--; + else break; + } //end-while + + int start = minGradIndex + 1; + while (start < index2) { + int r = segmentPoints[i][start].y; + int c = segmentPoints[i][start].x; + + if (gradImg[r*width + c] <= minGrad) start++; + else break; + } //end-while + + TestSegment(i, index1, end); + 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 +// +void EDPF::ExtractNewSegments() +{ + //vector *segments = &segmentPoints[segmentNos]; + vector< vector > validSegments; + int noSegments = 0; + + for (int i = 0; i < segmentNos; i++) { + int start = 0; + while (start < segmentPoints[i].size()) { + + while (start < segmentPoints[i].size()) { + int r = segmentPoints[i][start].y; + int c = segmentPoints[i][start].x; + + if (edgeImg[r*width + c]) break; + start++; + } //end-while + + int end = start + 1; + while (end < segmentPoints[i].size()) { + int r = segmentPoints[i][end].y; + int c = segmentPoints[i][end].x; + + if (edgeImg[r*width + c] == 0) break; + end++; + } //end-while + + int len = end - start; + 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()); + vector subVec(&segmentPoints[i][start], &segmentPoints[i][end - 1]); + validSegments[noSegments] = subVec; + noSegments++; + } //end-else + + start = end + 1; + } //end-while + } //end-for + + // Copy to ed + segmentPoints = validSegments; + + segmentNos = noSegments; +} + +//--------------------------------------------------------------------------- +// Number of false alarms code as suggested by Desolneux, Moisan and Morel (DMM) +// +double EDPF::NFA(double prob, int len) +{ + double nfa = np; + for (int i = 0; i EPSILON; i++) + nfa *= prob; + + return nfa; +} \ No newline at end of file diff --git a/EDPF.h b/EDPF.h new file mode 100644 index 0000000..06503a5 --- /dev/null +++ b/EDPF.h @@ -0,0 +1,55 @@ +/************************************************************************************************************** +* Edge Drawing (ED) and Edge Drawing Parameter Free (EDPF) source codes. +* Copyright (C) 2016, Cuneyt Akinlar & Cihan Topal +* E-mails of the authors: cuneytakinlar@gmail.com, cihant@anadolu.edu.tr +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. + +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. + +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . + +* By using this library you are implicitly assumed to have accepted all of the above statements, +* and accept to cite the following papers: +* +* [1] C. Topal and C. Akinlar, “Edge Drawing: A Combined Real-Time Edge and Segment Detector,” +* Journal of Visual Communication and Image Representation, 23(6), 862-872, DOI: 10.1016/j.jvcir.2012.05.004 (2012). +* +* [2] C. Akinlar and C. Topal, “EDPF: A Real-time Parameter-free Edge Segment Detector with a False Detection Control,” +* International Journal of Pattern Recognition and Artificial Intelligence, 26(1), DOI: 10.1142/S0218001412550026 (2012). +**************************************************************************************************************/ + +#ifndef _EDPF_ +#define _EDPF_ + +#include "ED.h" + +#define MAX_GRAD_VALUE 128*256 +#define EPSILON 1.0 + +class EDPF : public ED { +public: + EDPF(cv::Mat srcImage); + EDPF(ED obj); + EDPF(EDColor obj); +private: + double divForTestSegment; + double *H; + int np; + short *gradImg; + + void validateEdgeSegments(); + short *ComputePrewitt3x3(); // differs from base class's prewit function (calculates H) + void TestSegment(int i, int index1, int index2); + void ExtractNewSegments(); + double NFA(double prob, int len); +}; + +#endif // ! _EDPF_ diff --git a/NFA.cpp b/NFA.cpp new file mode 100644 index 0000000..80a67cd --- /dev/null +++ b/NFA.cpp @@ -0,0 +1,240 @@ +#include "NFA.h" +#include +#include + +NFALUT::NFALUT(int size, double _prob, double _logNT) +{ + LUTSize = size; + LUT = new int[LUTSize]; + + prob = _prob; + logNT = _logNT; + + LUT[0] = 1; + int j = 1; + for (int i = 1; i= 0) break; + } //end-while + + if (ret < 0) continue; + } //end-if + + LUT[i] = j; + } //end-for +} + + +NFALUT::~NFALUT() +{ + delete[] LUT; +} + +bool NFALUT::checkValidationByNFA(int n, int k) +{ + if (n >= LUTSize) + return nfa(n, k) >= 0.0; + else + return k >= LUT[n]; +} + +double NFALUT::myAtan2(double yy, double xx) +{ + static double LUT[MAX_LUT_SIZE + 1]; + static bool tableInited = false; + if (!tableInited) { + for (int i = 0; i <= MAX_LUT_SIZE; i++) { + LUT[i] = atan((double)i / MAX_LUT_SIZE); + } //end-for + + tableInited = true; + } //end-if + + double y = fabs(yy); + double x = fabs(xx); + + bool invert = false; + if (y > x) { + double t = x; + x = y; + y = t; + invert = true; + } //end-if + + double ratio; + if (x == 0) // avoid division error + x = 0.000001; + + ratio = y / x; + + double angle = LUT[(int)(ratio*MAX_LUT_SIZE)]; + + if (xx >= 0) { + if (yy >= 0) { + // I. quadrant + if (invert) angle = M_PI / 2 - angle; + + } + else { + // IV. quadrant + if (invert == false) angle = M_PI - angle; + else angle = M_PI / 2 + angle; + } //end-else + + } + else { + if (yy >= 0) { + /// II. quadrant + if (invert == false) angle = M_PI - angle; + else angle = M_PI / 2 + angle; + + } + else { + /// III. quadrant + if (invert) angle = M_PI / 2 - angle; + } //end-else + } //end-else + + return angle; +} + +double NFALUT::nfa(int n, int k) +{ + 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 */ + } //end-if + + /* 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) * (ii. + 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; + } //end-if + } //end-for + + return -log10(bin_tail) - logNT; +} + +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++) + { + a -= log(x + (double)n); + b += q[n] * pow(x, (double)n); + } + return a + log(b); +} + +double NFALUT::log_gamma_windschitl(double x) +{ + 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; + + /* 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); +} diff --git a/NFA.h b/NFA.h new file mode 100644 index 0000000..930ee6a --- /dev/null +++ b/NFA.h @@ -0,0 +1,51 @@ +#ifndef _NFA_ +#define _NFA_ + +#define TABSIZE 100000 + +//---------------------------------------------- +// Fast arctan2 using a lookup table +// +#define MAX_LUT_SIZE 1024 + +#ifndef TRUE +#define TRUE 1 +#endif /* !TRUE */ + +/** ln(10) */ +#ifndef M_LN10 +#define M_LN10 2.30258509299404568402 +#endif /* !M_LN10 */ + +/** PI */ +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif /* !M_PI */ + +#define RELATIVE_ERROR_FACTOR 100.0 + +// Lookup table (LUT) for NFA computation +class NFALUT { +public: + + NFALUT(int size, double _prob, double _logNT); + ~NFALUT(); + + int *LUT; // look up table + int LUTSize; + + double prob; + double logNT; + + bool checkValidationByNFA(int n, int k); + static double myAtan2(double yy, double xx); + +private: + double nfa(int n, int 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); +}; + +#endif \ No newline at end of file diff --git a/billiard.jpg b/billiard.jpg new file mode 100644 index 0000000..e2020a6 Binary files /dev/null and b/billiard.jpg differ diff --git a/test.cpp b/test.cpp new file mode 100644 index 0000000..d42d01f --- /dev/null +++ b/test.cpp @@ -0,0 +1,109 @@ +#include "EDLib.h" +#include + +using namespace cv; +using namespace std; + +int main() +{ + //***************************** ED Edge Segment Detection ***************************** + //Detection of edge segments from an input image + Mat testImg = imread("billiard.jpg", 0); + imshow("Source Image", testImg); + + //Call ED constructor + ED testED = ED(testImg, SOBEL_OPERATOR, 36, 8, 1, 10, 1.0, true); // apply ED algorithm + + //Show resulting edge image + Mat edgeImg = testED.getEdgeImage(); + imshow("Edge Image - PRESS ANY KEY TO CONTINUE", edgeImg); + waitKey(); + + //Output number of segments + int noSegments = testED.getSegmentNo(); + std::cout << "Number of edge segments: " << noSegments << std::endl; + + //Get edges in segment form (getSortedSegments() gives segments sorted w.r.t. legnths) + std::vector< std::vector > segments = testED.getSegments(); + + + //***************************** EDLINES Line Segment Detection ***************************** + //Detection of line segments from the same image + EDLines testEDLines = EDLines(testImg); + Mat lineImg = testEDLines.getLineImage(); //draws on an empty image + imshow("Line Image 1 - PRESS ANY KEY TO CONTINUE", lineImg); + + //Detection of lines segments from edge segments instead of input image + //Therefore, redundant detection of edge segmens can be avoided + testEDLines = EDLines(testED); + lineImg = testEDLines.drawOnImage(); //draws on the input image + imshow("Line Image 2 - PRESS ANY KEY TO CONTINUE", lineImg); + + //Acquiring line information, i.e. start & end points + vector lines = testEDLines.getLines(); + int noLines = testEDLines.getLinesNo(); + std::cout << "Number of line segments: " << noLines << std::endl; + + waitKey(); + + //************************** EDPF Parameter-free Edge Segment Detection ************************** + // Detection of edge segments with parameter free ED (EDPF) + + EDPF testEDPF = EDPF(testImg); + Mat edgePFImage = testEDPF.getEdgeImage(); + imshow("Edge Image Parameter Free", edgePFImage); + cout << "Number of edge segments found by EDPF: " << testEDPF.getSegmentNo() << endl; + waitKey(); + + //***************************** EDCIRCLES Circle Segment Detection ***************************** + //Detection of circles directly from the input image + + EDCircles testEDCircles = EDCircles(testImg); + Mat circleImg = testEDCircles.drawResult(false, ImageStyle::CIRCLES); + imshow("Circle Image 1", circleImg); + + //Detection of circles from already available EDPF or ED image + testEDCircles = EDCircles(testEDPF); + + //Get circle information as [cx, cy, r] + vector circles = testEDCircles.getCircles(); + + //Get ellipse information as [cx, cy, a, b, theta] + vector ellipses = testEDCircles.getEllipses(); + + //Circles and ellipses will be indicated in green and red, resp. + circleImg = testEDCircles.drawResult(true, ImageStyle::BOTH); + imshow("CIRCLES and ELLIPSES RESULT IMAGE", circleImg); + + int noCircles = testEDCircles.getCirclesNo(); + std::cout << "Number of circles: " << noCircles << std::endl; + waitKey(); + + //*********************** EDCOLOR Edge Segment Detection from Color Images ********************** + + Mat colorImg = imread("billiard.jpg"); + //Mat colorImg = imread("billiardNoise.jpg"); + EDColor testEDColor = EDColor(colorImg, 36, 4, 1.5, true); //last parameter for validation + imshow("Color Edge Image - PRESS ANY KEY TO QUIT", testEDColor.getEdgeImage()); + cout << "Number of edge segments detected by EDColor: " << testEDColor.getSegmentNo() << endl; + waitKey(); + + // get lines from color image + EDLines colorLine = EDLines(testEDColor); + imshow("Color Line", colorLine.getLineImage()); + std::cout << "Number of line segments: " << colorLine.getLinesNo() << std::endl; + waitKey(); + + // get circles from color image + EDCircles colorCircle = EDCircles(testEDColor); + // TO DO :: drawResult doesnt overlay (onImage = true) when input is from EDColor + circleImg = colorCircle.drawResult(false, ImageStyle::BOTH); + imshow("Color Circle", circleImg); + std::cout << "Number of line segments: " << colorCircle.getCirclesNo() << std::endl; + waitKey(); + + return 0; +} + + +