Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
254 views
in Technique[技术] by (71.8m points)

opencv - eigen3 select is too slow. am i make something wrong?

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.


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)
等待大神答复

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...