MatrixXd PointsToLineDist(Matrix<double, Dynamic, 2> mpoints, Matrix<double, Dynamic, 4>lines)
{
int lenpoints = mpoints.rows();
int lenlines = lines.rows();
MatrixXd innerLines = lines.transpose();
MatrixXd x1 = innerLines.row(0);
MatrixXd y1 = innerLines.row(1);
MatrixXd x2 = innerLines.row(2);
MatrixXd y2 = innerLines.row(3);
MatrixXd A = y2 - y1;
MatrixXd B = x1 - x2;
MatrixXd C = y1.array()*x2.array() - x1.array()*y2.array();
MatrixXd x0 = mpoints.col(0);
MatrixXd y0 = mpoints.col(1);
MatrixXd den = (A.array().pow(2) + B.array().pow(2)).sqrt();
//MatrixXd denm = den.replicate(lenpoints,1);
MatrixXd denm = MatrixXd::Ones(lenpoints, 1)*den;
//A = A.replicate(lenpoints,1);
//B = B.replicate(lenpoints,1);
//C = C.replicate(lenpoints,1);
//x0 = x0.replicate(1, lenlines);
//y0 = y0.replicate(1, lenlines);
//MatrixXd numerator = x0.replicate(1, lenlines).array()*A.replicate(lenpoints, 1).array() + y0.replicate(1, lenlines).array()*B.replicate(lenpoints, 1).array() + C.replicate(lenpoints, 1).array();
MatrixXd numerator = (x0*MatrixXd::Ones(1, lenlines)).array()*(MatrixXd::Ones(lenpoints, 1)*A).array() + (y0*MatrixXd::Ones(1, lenlines)).array()*(MatrixXd::Ones(lenpoints, 1)*B).array() + (MatrixXd::Ones(lenpoints, 1)*C).array();
MatrixXd dist = (numerator.array() / denm.array()).abs();
return dist;
}
MatrixXd PointsToSegDist(Matrix<double, Dynamic, 2> mpoints, Matrix<double, Dynamic, 4>lines)
{
struct timeb t1, t2;
long t;
ftime(&t1);
MatrixXd dp1 = PointsToLineDist(mpoints, lines);
ftime(&t2);
t = (t2.time - t1.time) * 1000 + (t2.millitm - t1.millitm);
cout << "PointsToLineDist eigen cost: " << t << endl;
int lenpoints = mpoints.rows();
int lenlines = lines.rows();
ftime(&t1);
//lines = lines';
MatrixXd innerlines = lines.transpose();
MatrixXd x0 = mpoints.col(0)*MatrixXd::Ones(1, lenlines);// .replicate(1, lenlines);
MatrixXd y0 = mpoints.col(1)*MatrixXd::Ones(1, lenlines);// .replicate(1, lenlines);
MatrixXd v1x = MatrixXd::Ones(lenpoints, 1)*innerlines.row(0);// .replicate(lenpoints, 1);
MatrixXd v1y = MatrixXd::Ones(lenpoints, 1)*innerlines.row(1);// innerlines.row(1).replicate(lenpoints, 1);
MatrixXd v2x = MatrixXd::Ones(lenpoints, 1)*innerlines.row(2);// innerlines.row(2).replicate(lenpoints, 1);
MatrixXd v2y = MatrixXd::Ones(lenpoints, 1)*innerlines.row(3);// innerlines.row(3).replicate(lenpoints, 1);
MatrixXd dsp = ((x0 - v1x).array().pow(2) + (y0 - v1y).array().pow(2)).sqrt();
MatrixXd dep = ((x0 - v2x).array().pow(2) + (y0 - v2y).array().pow(2)).sqrt();
MatrixXd dse = ((v1x - v2x).array().pow(2) + (v1y - v2y).array().pow(2)).sqrt();
ftime(&t2);
t = (t2.time - t1.time) * 1000 + (t2.millitm - t1.millitm);
cout << "dep dsp cost: " << t << endl;
auto type1 = (dse.array().pow(2) + dep.array().pow(2)).sqrt() <= dsp.array();
auto type2 = (dse.array().pow(2) + dsp.array().pow(2)).sqrt() <= dep.array();
auto type3 = type1!=true && type2!=true;
ftime(&t1);
MatrixXd distance = MatrixXd::Zero(lenpoints, lenlines);
distance = type1.select(dep, distance);
distance = type2.select(dsp, distance);
distance = type3.select(dp1, distance);
ftime(&t2);
t = (t2.time - t1.time) * 1000 + (t2.millitm - t1.millitm);
cout << "select cost: " << t <<endl;
return distance;
}
cv::Mat PointsToLineDist(cv::Mat mpoints, cv::Mat lines)
{
struct timeb t1, t2;
long t;
int lenpoints = mpoints.rows;
int lenlines = lines.rows;
cv::Mat innerLines; cv::transpose(lines, innerLines);
cv::Mat x1 = innerLines.row(0);
cv::Mat y1 = innerLines.row(1);
cv::Mat x2 = innerLines.row(2);
cv::Mat y2 = innerLines.row(3);
cv::Mat A = y2 - y1;
cv::Mat B = x1 - x2;
cv::Mat C = y1.mul(x2) - x1.mul(y2);
cv::Mat x0 = mpoints.col(0);
cv::Mat y0 = mpoints.col(1);
cv::Mat den;
cv::sqrt(A.mul(A) + B.mul(B), den);
//MatrixXd denm = den.replicate(lenpoints,1);
cv::Mat denm = cv::repeat(den, lenpoints, 1);// MatrixXd::Ones(lenpoints, 1)*den;
cv::Mat numerator = cv::repeat(x0, 1, lenlines).mul(cv::repeat(A, lenpoints, 1)) + cv::repeat(y0, 1, lenlines).mul(cv::repeat(B, lenpoints, 1)) + cv::repeat(C, lenpoints, 1);
//MatrixXd numerator = (x0*MatrixXd::Ones(1, lenlines)).array()*(MatrixXd::Ones(lenpoints, 1)*A).array() + (y0*MatrixXd::Ones(1, lenlines)).array()*(MatrixXd::Ones(lenpoints, 1)*B).array() + (MatrixXd::Ones(lenpoints, 1)*C).array();
cv::Mat dist = abs(numerator / denm);
return dist;
}
cv::Mat PointsToSegDist(cv::Mat mpoints, cv::Mat lines)
{
struct timeb t1, t2;
long t;
ftime(&t1);
cv::Mat dp1 = PointsToLineDist(mpoints, lines);
ftime(&t2);
t = (t2.time - t1.time) * 1000 + (t2.millitm - t1.millitm);
cout << "PointsToLineDist mat cost: " << t << endl;
int lenpoints = mpoints.rows;
int lenlines = lines.rows;
ftime(&t1);
//lines = lines';
cv::Mat innerLines; cv::transpose(lines, innerLines);
cv::Mat x0 = cv::repeat(mpoints.col(0), 1, lenlines);
cv::Mat y0 = cv::repeat(mpoints.col(1), 1, lenlines);// .replicate(1, lenlines);
cv::Mat v1x = cv::repeat(innerLines.row(0), lenpoints, 1);
cv::Mat v1y = cv::repeat(innerLines.row(1), lenpoints, 1);
cv::Mat v2x = cv::repeat(innerLines.row(2), lenpoints, 1);
cv::Mat v2y = cv::repeat(innerLines.row(3), lenpoints, 1);
cv::Mat dsp; cv::sqrt((x0 - v1x).mul(x0 - v1x) + (y0 - v1y).mul(y0 - v1y), dsp);
cv::Mat dep; cv::sqrt((x0 - v2x).mul(x0 - v2x) + (y0 - v2y).mul(y0 - v2y), dep);
cv::Mat dse; cv::sqrt((v1x - v2x).mul(v1x - v2x) + (v1y - v2y).mul(v1y - v2y), dse);
ftime(&t2);
t = (t2.time - t1.time) * 1000 + (t2.millitm - t1.millitm);
cout << "dep dsp cost: " << t << endl;
cv::Mat tmp1,tmp2;
cv::sqrt((dse.mul(dse) + dep.mul(dep)), tmp1);
const auto type1 = tmp1 <= dsp;
cv::sqrt((dse.mul(dse) + dsp.mul(dsp)), tmp2);
const auto type2 = tmp2 <= dep;
const auto type3 = type1 != 255 & type2 != 255;
ftime(&t1);
cv::Mat distance = cv::Mat::zeros(lenpoints, lenlines, CV_64FC1);
dep.copyTo(distance, type1);
dsp.copyTo(distance, type2);
dp1.copyTo(distance, type3);
ftime(&t2);
t = (t2.time - t1.time) * 1000 + (t2.millitm - t1.millitm);
cout << "select cost: " << t << endl;
return distance;
}
int main()
{
cv::Mat points, ob,lines;
struct timeb t1, t2;
long t;
//ob = cv::Mat::zeros(5000, 5000, CV_8UC1);
points.create(5000, 2, CV_64FC1);
lines.create(5000, 4, CV_64FC1);
cv::randu(points, 0, 20);
cv::randu(lines, 0, 20);
ftime(&t1);
cv::Mat ret = PointsToSegDist(points,lines);
ftime(&t2);
t = (t2.time - t1.time) * 1000 + (t2.millitm - t1.millitm);
cout << "PointsToSegDist mat cost: " << t << endl;
MatrixXd mxpoints, mxlines;
cv2eigen(points, mxpoints);
cv2eigen(lines, mxlines);
ftime(&t1);
MatrixXd mxret = PointsToSegDist(mxpoints, mxlines);
ftime(&t2);
t = (t2.time - t1.time) * 1000 + (t2.millitm - t1.millitm);
cout << "PointsToSegDist eigen cost: " << t << endl;
}
my computer return:
- PointsToLineDist mat cost: 1113
- dep dsp cost: 2766
- select cost: 361
- PointsToSegDist mat cost: 5222
- PointsToLineDist eigen cost: 822
- dep dsp cost: 6556
- select cost: 8529
- PointsToSegDist eigen cost: 16062
opencv condition cost 361 and eigen cost a lot.
i allso find that eigen replicate is relly slow
so i use multiply instead.