parent
d9f3ef31c8
commit
4b345d96ef
Binary file not shown.
@ -0,0 +1,275 @@
|
|||||||
|
#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.87 // 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<mCircle> getCircles();
|
||||||
|
std::vector<mEllipse> getEllipses();
|
||||||
|
int getCirclesNo();
|
||||||
|
int getEllipsesNo();
|
||||||
|
|
||||||
|
private:
|
||||||
|
int noEllipses;
|
||||||
|
int noCircles;
|
||||||
|
std::vector<mCircle> circles;
|
||||||
|
std::vector<mEllipse> 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<LineSegment> 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_
|
||||||
|
|
||||||
@ -0,0 +1,77 @@
|
|||||||
|
#ifndef _EDColor_
|
||||||
|
#define _EDColor_
|
||||||
|
|
||||||
|
#include <opencv2/opencv.hpp>
|
||||||
|
|
||||||
|
// 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<std::vector<cv::Point>> 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<std::vector<cv::Point>> segments;
|
||||||
|
|
||||||
|
static double LUT1[LUT_SIZE + 1];
|
||||||
|
static double LUT2[LUT_SIZE + 1];
|
||||||
|
static bool LUT_Initialized;
|
||||||
|
|
||||||
|
void MyRGB2LabFast();
|
||||||
|
void ComputeGradientMapByDiZenzo();
|
||||||
|
void smoothChannel(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<std::vector<cv::Point>> map, int noPixels);
|
||||||
|
|
||||||
|
static void InitColorEDLib();
|
||||||
|
};
|
||||||
|
|
||||||
|
#endif // ! _EDColor_
|
||||||
@ -0,0 +1,27 @@
|
|||||||
|
#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_
|
||||||
@ -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
|
||||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@ -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<std::vector<cv::Point>> 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; i<width*height; i++) {
|
||||||
|
red = redImg[i] / 255.0;
|
||||||
|
green = greenImg[i] / 255.0;
|
||||||
|
blue = blueImg[i] / 255.0;
|
||||||
|
|
||||||
|
red = LUT1[(int)(red*LUT_SIZE + 0.5)];
|
||||||
|
green = LUT1[(int)(green*LUT_SIZE + 0.5)];
|
||||||
|
blue = LUT1[(int)(blue*LUT_SIZE + 0.5)];
|
||||||
|
|
||||||
|
red = red * 100;
|
||||||
|
green = green * 100;
|
||||||
|
blue = blue * 100;
|
||||||
|
|
||||||
|
//Observer. = 2? Illuminant = D65
|
||||||
|
x = red*0.4124564 + green*0.3575761 + blue*0.1804375;
|
||||||
|
y = red*0.2126729 + green*0.7151522 + blue*0.0721750;
|
||||||
|
z = red*0.0193339 + green*0.1191920 + blue*0.9503041;
|
||||||
|
|
||||||
|
// Now xyz 2 Lab
|
||||||
|
double refX = 95.047;
|
||||||
|
double refY = 100.000;
|
||||||
|
double refZ = 108.883;
|
||||||
|
|
||||||
|
x = x / refX; //ref_X = 95.047 Observer= 2? Illuminant= D65
|
||||||
|
y = y / refY; //ref_Y = 100.000
|
||||||
|
z = z / refZ; //ref_Z = 108.883
|
||||||
|
|
||||||
|
x = LUT2[(int)(x*LUT_SIZE + 0.5)];
|
||||||
|
y = LUT2[(int)(y*LUT_SIZE + 0.5)];
|
||||||
|
z = LUT2[(int)(z*LUT_SIZE + 0.5)];
|
||||||
|
|
||||||
|
L[i] = (116.0*y) - 16;
|
||||||
|
a[i] = 500 * (x / y);
|
||||||
|
b[i] = 200 * (y - z);
|
||||||
|
} //end-for
|
||||||
|
|
||||||
|
// Scale L to [0-255]
|
||||||
|
double min = 1e10;
|
||||||
|
double max = -1e10;
|
||||||
|
for (int i = 0; i<width*height; i++) {
|
||||||
|
if (L[i]<min) min = L[i];
|
||||||
|
else if (L[i]>max) max = L[i];
|
||||||
|
} //end-for
|
||||||
|
|
||||||
|
double scale = 255.0 / (max - min);
|
||||||
|
for (int i = 0; i<width*height; i++) { L_Img[i] = (unsigned char)((L[i] - min)*scale); }
|
||||||
|
|
||||||
|
// Scale a to [0-255]
|
||||||
|
min = 1e10;
|
||||||
|
max = -1e10;
|
||||||
|
for (int i = 0; i<width*height; i++) {
|
||||||
|
if (a[i]<min) min = a[i];
|
||||||
|
else if (a[i]>max) max = a[i];
|
||||||
|
} //end-for
|
||||||
|
|
||||||
|
scale = 255.0 / (max - min);
|
||||||
|
for (int i = 0; i<width*height; i++) { a_Img[i] = (unsigned char)((a[i] - min)*scale); }
|
||||||
|
|
||||||
|
// Scale b to [0-255]
|
||||||
|
min = 1e10;
|
||||||
|
max = -1e10;
|
||||||
|
for (int i = 0; i<width*height; i++) {
|
||||||
|
if (b[i]<min) min = b[i];
|
||||||
|
else if (b[i]>max) max = b[i];
|
||||||
|
} //end-for
|
||||||
|
|
||||||
|
scale = 255.0 / (max - min);
|
||||||
|
for (int i = 0; i<width*height; i++) { b_Img[i] = (unsigned char)((b[i] - min)*scale); }
|
||||||
|
|
||||||
|
|
||||||
|
// clean space
|
||||||
|
delete[] L;
|
||||||
|
delete[] a;
|
||||||
|
delete[] b;
|
||||||
|
}
|
||||||
|
|
||||||
|
void EDColor::ComputeGradientMapByDiZenzo()
|
||||||
|
{
|
||||||
|
memset(gradImg, 0, sizeof(short)*width*height);
|
||||||
|
|
||||||
|
int max = 0;
|
||||||
|
|
||||||
|
for (int i = 1; i < height - 1; i++) {
|
||||||
|
for (int j = 1; j < width - 1; j++) {
|
||||||
|
#if 1
|
||||||
|
// Prewitt for channel1
|
||||||
|
int com1 = smooth_L[(i + 1)*width + j + 1] - smooth_L[(i - 1)*width + j - 1];
|
||||||
|
int com2 = smooth_L[(i - 1)*width + j + 1] - smooth_L[(i + 1)*width + j - 1];
|
||||||
|
|
||||||
|
int gxCh1 = com1 + com2 + (smooth_L[i*width + j + 1] - smooth_L[i*width + j - 1]);
|
||||||
|
int gyCh1 = com1 - com2 + (smooth_L[(i + 1)*width + j] - smooth_L[(i - 1)*width + j]);
|
||||||
|
|
||||||
|
// Prewitt for channel2
|
||||||
|
com1 = smooth_a[(i + 1)*width + j + 1] - smooth_a[(i - 1)*width + j - 1];
|
||||||
|
com2 = smooth_a[(i - 1)*width + j + 1] - smooth_a[(i + 1)*width + j - 1];
|
||||||
|
|
||||||
|
int gxCh2 = com1 + com2 + (smooth_a[i*width + j + 1] - smooth_a[i*width + j - 1]);
|
||||||
|
int gyCh2 = com1 - com2 + (smooth_a[(i + 1)*width + j] - smooth_a[(i - 1)*width + j]);
|
||||||
|
|
||||||
|
// Prewitt for channel3
|
||||||
|
com1 = smooth_b[(i + 1)*width + j + 1] - smooth_b[(i - 1)*width + j - 1];
|
||||||
|
com2 = smooth_b[(i - 1)*width + j + 1] - smooth_b[(i + 1)*width + j - 1];
|
||||||
|
|
||||||
|
int gxCh3 = com1 + com2 + (smooth_b[i*width + j + 1] - smooth_b[i*width + j - 1]);
|
||||||
|
int gyCh3 = com1 - com2 + (smooth_b[(i + 1)*width + j] - smooth_b[(i - 1)*width + j]);
|
||||||
|
#else
|
||||||
|
// Sobel for channel1
|
||||||
|
int com1 = smooth_L[(i + 1)*width + j + 1] - smooth_L[(i - 1)*width + j - 1];
|
||||||
|
int com2 = smooth_L[(i - 1)*width + j + 1] - smooth_L[(i + 1)*width + j - 1];
|
||||||
|
|
||||||
|
int gxCh1 = com1 + com2 + 2 * (smooth_L[i*width + j + 1] - smooth_L[i*width + j - 1]);
|
||||||
|
int gyCh1 = com1 - com2 + 2 * (smooth_L[(i + 1)*width + j] - smooth_L[(i - 1)*width + j]);
|
||||||
|
|
||||||
|
// Sobel for channel2
|
||||||
|
com1 = smooth_a[(i + 1)*width + j + 1] - smooth_a[(i - 1)*width + j - 1];
|
||||||
|
com2 = smooth_a[(i - 1)*width + j + 1] - smooth_a[(i + 1)*width + j - 1];
|
||||||
|
|
||||||
|
int gxCh2 = com1 + com2 + 2 * (smooth_a[i*width + j + 1] - smooth_a[i*width + j - 1]);
|
||||||
|
int gyCh2 = com1 - com2 + 2 * (smooth_a[(i + 1)*width + j] - smooth_a[(i - 1)*width + j]);
|
||||||
|
|
||||||
|
// Sobel for channel3
|
||||||
|
com1 = smooth_b[(i + 1)*width + j + 1] - smooth_b[(i - 1)*width + j - 1];
|
||||||
|
com2 = smooth_b[(i - 1)*width + j + 1] - smooth_b[(i + 1)*width + j - 1];
|
||||||
|
|
||||||
|
int gxCh3 = com1 + com2 + 2 * (smooth_b[i*width + j + 1] - smooth_b[i*width + j - 1]);
|
||||||
|
int gyCh3 = com1 - com2 + 2 * (smooth_b[(i + 1)*width + j] - smooth_b[(i - 1)*width + j]);
|
||||||
|
#endif
|
||||||
|
int gxx = gxCh1*gxCh1 + gxCh2*gxCh2 + gxCh3*gxCh3;
|
||||||
|
int gyy = gyCh1*gyCh1 + gyCh2*gyCh2 + gyCh3*gyCh3;
|
||||||
|
int gxy = gxCh1*gyCh1 + gxCh2*gyCh2 + gxCh3*gyCh3;
|
||||||
|
|
||||||
|
#if 1
|
||||||
|
// Di Zenzo's formulas from Gonzales & Woods - Page 337
|
||||||
|
double theta = atan2(2.0*gxy, (double)(gxx - gyy)) / 2; // Gradient Direction
|
||||||
|
int grad = (int)(sqrt(((gxx + gyy) + (gxx - gyy)*cos(2 * theta) + 2 * gxy*sin(2 * theta)) / 2.0) + 0.5); // Gradient Magnitude
|
||||||
|
#else
|
||||||
|
// Koschan & Abidi - 2005 - Signal Processing Magazine
|
||||||
|
double theta = atan2(2.0*gxy, (double)(gxx - gyy)) / 2; // Gradient Direction
|
||||||
|
|
||||||
|
double cosTheta = cos(theta);
|
||||||
|
double sinTheta = sin(theta);
|
||||||
|
int grad = (int)(sqrt(gxx*cosTheta*cosTheta + 2 * gxy*sinTheta*cosTheta + gyy*sinTheta*sinTheta) + 0.5); // Gradient Magnitude
|
||||||
|
#endif
|
||||||
|
|
||||||
|
// Gradient is perpendicular to the edge passing through the pixel
|
||||||
|
if (theta >= -3.14159 / 4 && theta <= 3.14159 / 4)
|
||||||
|
dirImg[i*width + j] = EDGE_VERTICAL;
|
||||||
|
else
|
||||||
|
dirImg[i*width + j] = EDGE_HORIZONTAL;
|
||||||
|
|
||||||
|
gradImg[i*width + j] = 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; i<width*height; i++)
|
||||||
|
gradImg[i] = (short)(gradImg[i] * scale);
|
||||||
|
}
|
||||||
|
|
||||||
|
void EDColor::smoothChannel(uchar *src, uchar *smooth, double sigma)
|
||||||
|
{
|
||||||
|
Mat srcImage = Mat(height, width, CV_8UC1, src);
|
||||||
|
Mat smoothImage = Mat(height, width, CV_8UC1, smooth);
|
||||||
|
|
||||||
|
if (sigma == 1.0)
|
||||||
|
GaussianBlur(srcImage, smoothImage, Size(5, 5), 1);
|
||||||
|
else if (sigma == 1.5)
|
||||||
|
GaussianBlur(srcImage, smoothImage, Size(7, 7), 1.5); // seems to be better?
|
||||||
|
else
|
||||||
|
GaussianBlur(srcImage, smoothImage, Size(), sigma);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
//--------------------------------------------------------------------------------------------------------------------
|
||||||
|
// Validate the edge segments using the Helmholtz principle (for color images) channel1, channel2 and channel3 images
|
||||||
|
//
|
||||||
|
void EDColor::validateEdgeSegments()
|
||||||
|
{
|
||||||
|
int maxGradValue = MAX_GRAD_VALUE;
|
||||||
|
H = new double[maxGradValue];
|
||||||
|
memset(H, 0, sizeof(double)*maxGradValue);
|
||||||
|
|
||||||
|
memset(edgeImg, 0, width*height); // clear edge image
|
||||||
|
|
||||||
|
// Compute the gradient
|
||||||
|
memset(gradImg, 0, sizeof(short)*width*height); // reset gradient Image pixels to zero
|
||||||
|
|
||||||
|
int *grads = new int[maxGradValue];
|
||||||
|
memset(grads, 0, sizeof(int)*maxGradValue);
|
||||||
|
|
||||||
|
for (int i = 1; i<height - 1; i++) {
|
||||||
|
for (int j = 1; j<width - 1; j++) {
|
||||||
|
// Gradient for channel1
|
||||||
|
int com1 = smooth_L[(i + 1)*width + j + 1] - smooth_L[(i - 1)*width + j - 1];
|
||||||
|
int com2 = smooth_L[(i - 1)*width + j + 1] - smooth_L[(i + 1)*width + j - 1];
|
||||||
|
|
||||||
|
int gxCh1 = abs(com1 + com2 + (smooth_L[i*width + j + 1] - smooth_L[i*width + j - 1]));
|
||||||
|
int gyCh1 = abs(com1 - com2 + (smooth_L[(i + 1)*width + j] - smooth_L[(i - 1)*width + j]));
|
||||||
|
int ch1Grad = gxCh1 + gyCh1;
|
||||||
|
|
||||||
|
// Gradient for channel2
|
||||||
|
com1 = smooth_a[(i + 1)*width + j + 1] - smooth_a[(i - 1)*width + j - 1];
|
||||||
|
com2 = smooth_a[(i - 1)*width + j + 1] - smooth_a[(i + 1)*width + j - 1];
|
||||||
|
|
||||||
|
int gxCh2 = abs(com1 + com2 + (smooth_a[i*width + j + 1] - smooth_a[i*width + j - 1]));
|
||||||
|
int gyCh2 = abs(com1 - com2 + (smooth_a[(i + 1)*width + j] - smooth_a[(i - 1)*width + j]));
|
||||||
|
int ch2Grad = gxCh2 + gyCh2;
|
||||||
|
|
||||||
|
// Gradient for channel3
|
||||||
|
com1 = smooth_b[(i + 1)*width + j + 1] - smooth_b[(i - 1)*width + j - 1];
|
||||||
|
com2 = smooth_b[(i - 1)*width + j + 1] - smooth_b[(i + 1)*width + j - 1];
|
||||||
|
|
||||||
|
int gxCh3 = abs(com1 + com2 + (smooth_b[i*width + j + 1] - smooth_b[i*width + j - 1]));
|
||||||
|
int gyCh3 = abs(com1 - com2 + (smooth_b[(i + 1)*width + j] - smooth_b[(i - 1)*width + j]));
|
||||||
|
int ch3Grad = gxCh3 + gyCh3;
|
||||||
|
|
||||||
|
// Take average
|
||||||
|
int grad = (ch1Grad + ch2Grad + ch3Grad + 2) / 3;
|
||||||
|
|
||||||
|
gradImg[i*width + j] = grad;
|
||||||
|
grads[grad]++;
|
||||||
|
} //end-for
|
||||||
|
} //end-for
|
||||||
|
|
||||||
|
Mat gradImage = Mat(height, width, CV_16SC1, gradImg);
|
||||||
|
imwrite("newGrad.pgm", gradImage);
|
||||||
|
|
||||||
|
// Compute probability function H
|
||||||
|
int size = (width - 2)*(height - 2);
|
||||||
|
// size -= grads[0];
|
||||||
|
|
||||||
|
for (int i = maxGradValue - 1; i>0; i--)
|
||||||
|
grads[i - 1] += grads[i];
|
||||||
|
|
||||||
|
for (int i = 0; i<maxGradValue; i++)
|
||||||
|
H[i] = (double)grads[i] / ((double)size);
|
||||||
|
|
||||||
|
// Compute np: # of segment pieces
|
||||||
|
np = 0;
|
||||||
|
for (int i = 0; i<segments.size(); i++) {
|
||||||
|
int len = segments[i].size();
|
||||||
|
np += (len*(len - 1)) / 2;
|
||||||
|
} //end-for
|
||||||
|
|
||||||
|
// Validate segments
|
||||||
|
for (int i = 0; i< segments.size(); i++) {
|
||||||
|
testSegment(i, 0, segments[i].size() - 1);
|
||||||
|
} //end-for
|
||||||
|
|
||||||
|
// clear space
|
||||||
|
delete[] H;
|
||||||
|
delete[] grads;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
//----------------------------------------------------------------------------------
|
||||||
|
// 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 EDColor::testSegment(int i, int index1, int index2)
|
||||||
|
{
|
||||||
|
|
||||||
|
int chainLen = index2 - index1 + 1;
|
||||||
|
if (chainLen < MIN_PATH_LEN)
|
||||||
|
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 = segments[i][k].y;
|
||||||
|
int c = segments[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 = segments[i][k].y;
|
||||||
|
int c = segments[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 = 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<Point> > 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<Point>());
|
||||||
|
vector<Point> 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<len && nfa > 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<std::vector<cv::Point>> 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};
|
||||||
File diff suppressed because it is too large
Load Diff
@ -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; i<height - 1; i++) {
|
||||||
|
for (int j = 1; j<width - 1; j++) {
|
||||||
|
// Prewitt Operator in horizontal and vertical direction
|
||||||
|
// 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 = com1 - com2 + (G-B) = (H-A) - (C-F) + (G-B) = (F-A) + (G-B) + (H-C)
|
||||||
|
//
|
||||||
|
int com1 = srcImg[(i + 1)*width + j + 1] - srcImg[(i - 1)*width + j - 1];
|
||||||
|
int com2 = srcImg[(i - 1)*width + j + 1] - srcImg[(i + 1)*width + j - 1];
|
||||||
|
|
||||||
|
int gx = abs(com1 + com2 + (srcImg[i*width + j + 1] - srcImg[i*width + j - 1]));
|
||||||
|
int gy = abs(com1 - com2 + (srcImg[(i + 1)*width + j] - srcImg[(i - 1)*width + j]));
|
||||||
|
|
||||||
|
int g = gx + gy;
|
||||||
|
|
||||||
|
gradImg[i*width + j] = g;
|
||||||
|
grads[g]++;
|
||||||
|
} // end-for
|
||||||
|
} //end-for
|
||||||
|
|
||||||
|
// Compute probability function H
|
||||||
|
int size = (width - 2)*(height - 2);
|
||||||
|
|
||||||
|
for (int i = MAX_GRAD_VALUE - 1; i>0; 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<Point> *segments = &segmentPoints[segmentNos];
|
||||||
|
vector< vector<Point> > 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<Point>());
|
||||||
|
vector<Point> 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<len && nfa > EPSILON; i++)
|
||||||
|
nfa *= prob;
|
||||||
|
|
||||||
|
return nfa;
|
||||||
|
}
|
||||||
@ -0,0 +1,240 @@
|
|||||||
|
#include "NFA.h"
|
||||||
|
#include <math.h>
|
||||||
|
#include <float.h>
|
||||||
|
//Number of False Alarms
|
||||||
|
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<LUTSize; i++) {
|
||||||
|
LUT[i] = LUTSize + 1;
|
||||||
|
double ret = nfa(i, j);
|
||||||
|
if (ret < 0) {
|
||||||
|
while (j < i) {
|
||||||
|
j++;
|
||||||
|
ret = nfa(i, j);
|
||||||
|
if (ret >= 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) * (i<TABSIZE ?
|
||||||
|
(inv[i] != 0.0 ? inv[i] : (inv[i] = 1.0 / (double)i)) :
|
||||||
|
1.0 / (double)i);
|
||||||
|
|
||||||
|
mult_term = bin_term * p_term;
|
||||||
|
term *= mult_term;
|
||||||
|
bin_tail += term;
|
||||||
|
|
||||||
|
if (bin_term<1.0) {
|
||||||
|
/* When bin_term<1 then mult_term_j<mult_term_i for j>i.
|
||||||
|
Then, the error on the binomial tail when truncated at
|
||||||
|
the i term can be bounded by a geometric series of form
|
||||||
|
term_i * sum mult_term_i^j. */
|
||||||
|
err = term * ((1.0 - pow(mult_term, (double)(n - i + 1))) /
|
||||||
|
(1.0 - mult_term) - 1.0);
|
||||||
|
|
||||||
|
/* One wants an error at most of tolerance*final_result, or:
|
||||||
|
tolerance * abs(-log10(bin_tail)-logNT).
|
||||||
|
Now, the error that can be accepted on bin_tail is
|
||||||
|
given by tolerance*final_result divided by the derivative
|
||||||
|
of -log10(x) when x=bin_tail. that is:
|
||||||
|
tolerance * abs(-log10(bin_tail)-logNT) / (1/bin_tail)
|
||||||
|
Finally, we truncate the tail if the error is less than:
|
||||||
|
tolerance * abs(-log10(bin_tail)-logNT) * bin_tail */
|
||||||
|
if (err < tolerance * fabs(-log10(bin_tail) - logNT) * bin_tail) break;
|
||||||
|
} //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);
|
||||||
|
}
|
||||||
@ -0,0 +1,32 @@
|
|||||||
|
#ifndef _H_QFUNCTIONTRANSFER_H_
|
||||||
|
#define _H_QFUNCTIONTRANSFER_H_
|
||||||
|
|
||||||
|
#include <QObject>
|
||||||
|
#include <functional>
|
||||||
|
|
||||||
|
class QFunctionTransfer : public QObject
|
||||||
|
{
|
||||||
|
Q_OBJECT
|
||||||
|
|
||||||
|
public:
|
||||||
|
explicit QFunctionTransfer(QObject *parent = 0) {
|
||||||
|
qRegisterMetaType<std::tr1::function<void()>>("std::tr1::function<void()>");
|
||||||
|
connect(this, &QFunctionTransfer::comming, this, &QFunctionTransfer::exec, Qt::QueuedConnection);
|
||||||
|
};
|
||||||
|
void execInMain(std::tr1::function<void()> f) {
|
||||||
|
emit this->comming(f);
|
||||||
|
};
|
||||||
|
static QFunctionTransfer* Instance() {
|
||||||
|
static QFunctionTransfer ins;
|
||||||
|
return &ins;
|
||||||
|
}
|
||||||
|
signals:
|
||||||
|
void comming(std::tr1::function<void()> f);
|
||||||
|
public slots:
|
||||||
|
void exec(std::tr1::function<void()> f)
|
||||||
|
{
|
||||||
|
f();
|
||||||
|
};
|
||||||
|
};
|
||||||
|
|
||||||
|
#endif
|
||||||
@ -0,0 +1,68 @@
|
|||||||
|
<?xml version="1.0" encoding="utf-8"?>
|
||||||
|
<Project ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
|
||||||
|
<ItemGroup>
|
||||||
|
<Filter Include="Source Files">
|
||||||
|
<UniqueIdentifier>{4FC737F1-C7A5-4376-A066-2A32D752A2FF}</UniqueIdentifier>
|
||||||
|
<Extensions>cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx</Extensions>
|
||||||
|
</Filter>
|
||||||
|
<Filter Include="Header Files">
|
||||||
|
<UniqueIdentifier>{93995380-89BD-4b04-88EB-625FBE52EBFB}</UniqueIdentifier>
|
||||||
|
<Extensions>h;hh;hpp;hxx;hm;inl;inc;xsd</Extensions>
|
||||||
|
</Filter>
|
||||||
|
<Filter Include="Resource Files">
|
||||||
|
<UniqueIdentifier>{D9D6E242-F8AF-46E4-B9FD-80ECBC20BA3E}</UniqueIdentifier>
|
||||||
|
<Extensions>qrc;*</Extensions>
|
||||||
|
<ParseFiles>false</ParseFiles>
|
||||||
|
</Filter>
|
||||||
|
<Filter Include="Resource Files">
|
||||||
|
<UniqueIdentifier>{D9D6E242-F8AF-46E4-B9FD-80ECBC20BA3E}</UniqueIdentifier>
|
||||||
|
<Extensions>qrc;*</Extensions>
|
||||||
|
<ParseFiles>false</ParseFiles>
|
||||||
|
</Filter>
|
||||||
|
<Filter Include="Generated Files">
|
||||||
|
<UniqueIdentifier>{71ED8ED8-ACB9-4CE9-BBE1-E00B30144E11}</UniqueIdentifier>
|
||||||
|
<Extensions>moc;h;cpp</Extensions>
|
||||||
|
<SourceControlFiles>False</SourceControlFiles>
|
||||||
|
</Filter>
|
||||||
|
</ItemGroup>
|
||||||
|
<ItemGroup>
|
||||||
|
<ClInclude Include="..\..\src\edcircle\include\ED.h">
|
||||||
|
<Filter>Header Files</Filter>
|
||||||
|
</ClInclude>
|
||||||
|
<ClInclude Include="..\..\src\edcircle\include\EDCircles.h">
|
||||||
|
<Filter>Header Files</Filter>
|
||||||
|
</ClInclude>
|
||||||
|
<ClInclude Include="..\..\src\edcircle\include\EDColor.h">
|
||||||
|
<Filter>Header Files</Filter>
|
||||||
|
</ClInclude>
|
||||||
|
<ClInclude Include="..\..\src\edcircle\include\EDLines.h">
|
||||||
|
<Filter>Header Files</Filter>
|
||||||
|
</ClInclude>
|
||||||
|
<ClInclude Include="..\..\src\edcircle\include\EDPF.h">
|
||||||
|
<Filter>Header Files</Filter>
|
||||||
|
</ClInclude>
|
||||||
|
<ClInclude Include="..\..\src\edcircle\include\NFA.h">
|
||||||
|
<Filter>Header Files</Filter>
|
||||||
|
</ClInclude>
|
||||||
|
</ItemGroup>
|
||||||
|
<ItemGroup>
|
||||||
|
<ClCompile Include="..\..\src\edcircle\src\ED.cpp">
|
||||||
|
<Filter>Source Files</Filter>
|
||||||
|
</ClCompile>
|
||||||
|
<ClCompile Include="..\..\src\edcircle\src\EDCircles.cpp">
|
||||||
|
<Filter>Source Files</Filter>
|
||||||
|
</ClCompile>
|
||||||
|
<ClCompile Include="..\..\src\edcircle\src\EDColor.cpp">
|
||||||
|
<Filter>Source Files</Filter>
|
||||||
|
</ClCompile>
|
||||||
|
<ClCompile Include="..\..\src\edcircle\src\EDLines.cpp">
|
||||||
|
<Filter>Source Files</Filter>
|
||||||
|
</ClCompile>
|
||||||
|
<ClCompile Include="..\..\src\edcircle\src\EDPF.cpp">
|
||||||
|
<Filter>Source Files</Filter>
|
||||||
|
</ClCompile>
|
||||||
|
<ClCompile Include="..\..\src\edcircle\src\NFA.cpp">
|
||||||
|
<Filter>Source Files</Filter>
|
||||||
|
</ClCompile>
|
||||||
|
</ItemGroup>
|
||||||
|
</Project>
|
||||||
@ -0,0 +1,131 @@
|
|||||||
|
<?xml version="1.0" encoding="utf-8"?>
|
||||||
|
<Project DefaultTargets="Build" ToolsVersion="15.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
|
||||||
|
<ItemGroup Label="ProjectConfigurations">
|
||||||
|
<ProjectConfiguration Include="Debug|Win32">
|
||||||
|
<Configuration>Debug</Configuration>
|
||||||
|
<Platform>Win32</Platform>
|
||||||
|
</ProjectConfiguration>
|
||||||
|
<ProjectConfiguration Include="Release|Win32">
|
||||||
|
<Configuration>Release</Configuration>
|
||||||
|
<Platform>Win32</Platform>
|
||||||
|
</ProjectConfiguration>
|
||||||
|
<ProjectConfiguration Include="Debug|x64">
|
||||||
|
<Configuration>Debug</Configuration>
|
||||||
|
<Platform>x64</Platform>
|
||||||
|
</ProjectConfiguration>
|
||||||
|
<ProjectConfiguration Include="Release|x64">
|
||||||
|
<Configuration>Release</Configuration>
|
||||||
|
<Platform>x64</Platform>
|
||||||
|
</ProjectConfiguration>
|
||||||
|
</ItemGroup>
|
||||||
|
<PropertyGroup Label="Globals">
|
||||||
|
<VCProjectVersion>15.0</VCProjectVersion>
|
||||||
|
<ProjectGuid>{4FDC5305-11FB-496A-BEBE-828569051007}</ProjectGuid>
|
||||||
|
<RootNamespace>My1EDCircle</RootNamespace>
|
||||||
|
<WindowsTargetPlatformVersion>10.0.17763.0</WindowsTargetPlatformVersion>
|
||||||
|
</PropertyGroup>
|
||||||
|
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.Default.props" />
|
||||||
|
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'" Label="Configuration">
|
||||||
|
<ConfigurationType>Application</ConfigurationType>
|
||||||
|
<UseDebugLibraries>true</UseDebugLibraries>
|
||||||
|
<PlatformToolset>v141</PlatformToolset>
|
||||||
|
<CharacterSet>MultiByte</CharacterSet>
|
||||||
|
</PropertyGroup>
|
||||||
|
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'" Label="Configuration">
|
||||||
|
<ConfigurationType>Application</ConfigurationType>
|
||||||
|
<UseDebugLibraries>false</UseDebugLibraries>
|
||||||
|
<PlatformToolset>v141</PlatformToolset>
|
||||||
|
<WholeProgramOptimization>true</WholeProgramOptimization>
|
||||||
|
<CharacterSet>MultiByte</CharacterSet>
|
||||||
|
</PropertyGroup>
|
||||||
|
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
|
||||||
|
<ConfigurationType>Application</ConfigurationType>
|
||||||
|
<UseDebugLibraries>true</UseDebugLibraries>
|
||||||
|
<PlatformToolset>v141</PlatformToolset>
|
||||||
|
<CharacterSet>MultiByte</CharacterSet>
|
||||||
|
</PropertyGroup>
|
||||||
|
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
|
||||||
|
<ConfigurationType>Application</ConfigurationType>
|
||||||
|
<UseDebugLibraries>false</UseDebugLibraries>
|
||||||
|
<PlatformToolset>v141</PlatformToolset>
|
||||||
|
<WholeProgramOptimization>true</WholeProgramOptimization>
|
||||||
|
<CharacterSet>MultiByte</CharacterSet>
|
||||||
|
</PropertyGroup>
|
||||||
|
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.props" />
|
||||||
|
<ImportGroup Label="ExtensionSettings">
|
||||||
|
</ImportGroup>
|
||||||
|
<ImportGroup Label="Shared">
|
||||||
|
</ImportGroup>
|
||||||
|
<ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
|
||||||
|
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
|
||||||
|
</ImportGroup>
|
||||||
|
<ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
|
||||||
|
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
|
||||||
|
</ImportGroup>
|
||||||
|
<ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
|
||||||
|
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
|
||||||
|
</ImportGroup>
|
||||||
|
<ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
|
||||||
|
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
|
||||||
|
</ImportGroup>
|
||||||
|
<PropertyGroup Label="UserMacros" />
|
||||||
|
<PropertyGroup />
|
||||||
|
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
|
||||||
|
<ClCompile>
|
||||||
|
<WarningLevel>Level3</WarningLevel>
|
||||||
|
<Optimization>Disabled</Optimization>
|
||||||
|
<SDLCheck>true</SDLCheck>
|
||||||
|
<ConformanceMode>true</ConformanceMode>
|
||||||
|
</ClCompile>
|
||||||
|
</ItemDefinitionGroup>
|
||||||
|
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
|
||||||
|
<ClCompile>
|
||||||
|
<WarningLevel>Level3</WarningLevel>
|
||||||
|
<Optimization>Disabled</Optimization>
|
||||||
|
<SDLCheck>true</SDLCheck>
|
||||||
|
<ConformanceMode>true</ConformanceMode>
|
||||||
|
<AdditionalIncludeDirectories>.\..\..\..\lpopencv\build\include;.\..\..\..\lpopencv\build\include\opencv;.\..\..\..\lpopencv\build\include\opencv2;F:\LeaperProject\1EDCircle\srccode\include;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
|
||||||
|
</ClCompile>
|
||||||
|
<Link>
|
||||||
|
<AdditionalLibraryDirectories>.\..\..\..\lpopencv\build\x64\vc15\lib;F:\LeaperProject\1EDCircle\tp_2017\x64\Debug;%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
|
||||||
|
<AdditionalDependencies>opencv_world341d.lib;EDCircled.lib;%(AdditionalDependencies)</AdditionalDependencies>
|
||||||
|
</Link>
|
||||||
|
</ItemDefinitionGroup>
|
||||||
|
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
|
||||||
|
<ClCompile>
|
||||||
|
<WarningLevel>Level3</WarningLevel>
|
||||||
|
<Optimization>MaxSpeed</Optimization>
|
||||||
|
<FunctionLevelLinking>true</FunctionLevelLinking>
|
||||||
|
<IntrinsicFunctions>true</IntrinsicFunctions>
|
||||||
|
<SDLCheck>true</SDLCheck>
|
||||||
|
<ConformanceMode>true</ConformanceMode>
|
||||||
|
</ClCompile>
|
||||||
|
<Link>
|
||||||
|
<EnableCOMDATFolding>true</EnableCOMDATFolding>
|
||||||
|
<OptimizeReferences>true</OptimizeReferences>
|
||||||
|
</Link>
|
||||||
|
</ItemDefinitionGroup>
|
||||||
|
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
|
||||||
|
<ClCompile>
|
||||||
|
<WarningLevel>Level3</WarningLevel>
|
||||||
|
<Optimization>MaxSpeed</Optimization>
|
||||||
|
<FunctionLevelLinking>true</FunctionLevelLinking>
|
||||||
|
<IntrinsicFunctions>true</IntrinsicFunctions>
|
||||||
|
<SDLCheck>true</SDLCheck>
|
||||||
|
<ConformanceMode>true</ConformanceMode>
|
||||||
|
<AdditionalIncludeDirectories>.\..\..\lpopencv\build\include;.\..\..\lpopencv\build\include\opencv;.\..\..\lpopencv\build\include\opencv2;.\..\EDCircle\include;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
|
||||||
|
</ClCompile>
|
||||||
|
<Link>
|
||||||
|
<EnableCOMDATFolding>true</EnableCOMDATFolding>
|
||||||
|
<OptimizeReferences>true</OptimizeReferences>
|
||||||
|
<AdditionalDependencies>opencv_world341.lib;EDCircle.lib;%(AdditionalDependencies)</AdditionalDependencies>
|
||||||
|
<AdditionalLibraryDirectories>.\..\..\lpopencv\build\x64\vc15\lib;.\..\x64\Debug;%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
|
||||||
|
</Link>
|
||||||
|
</ItemDefinitionGroup>
|
||||||
|
<ItemGroup>
|
||||||
|
<ClCompile Include="1_EDCircle.cpp" />
|
||||||
|
</ItemGroup>
|
||||||
|
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
|
||||||
|
<ImportGroup Label="ExtensionTargets">
|
||||||
|
</ImportGroup>
|
||||||
|
</Project>
|
||||||
@ -0,0 +1,22 @@
|
|||||||
|
<?xml version="1.0" encoding="utf-8"?>
|
||||||
|
<Project ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
|
||||||
|
<ItemGroup>
|
||||||
|
<Filter Include="源文件">
|
||||||
|
<UniqueIdentifier>{4FC737F1-C7A5-4376-A066-2A32D752A2FF}</UniqueIdentifier>
|
||||||
|
<Extensions>cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx</Extensions>
|
||||||
|
</Filter>
|
||||||
|
<Filter Include="头文件">
|
||||||
|
<UniqueIdentifier>{93995380-89BD-4b04-88EB-625FBE52EBFB}</UniqueIdentifier>
|
||||||
|
<Extensions>h;hh;hpp;hxx;hm;inl;inc;ipp;xsd</Extensions>
|
||||||
|
</Filter>
|
||||||
|
<Filter Include="资源文件">
|
||||||
|
<UniqueIdentifier>{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}</UniqueIdentifier>
|
||||||
|
<Extensions>rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms</Extensions>
|
||||||
|
</Filter>
|
||||||
|
</ItemGroup>
|
||||||
|
<ItemGroup>
|
||||||
|
<ClCompile Include="1_EDCircle.cpp">
|
||||||
|
<Filter>源文件</Filter>
|
||||||
|
</ClCompile>
|
||||||
|
</ItemGroup>
|
||||||
|
</Project>
|
||||||
Loading…
Reference in New Issue