SingularityJudger/main.cpp

#include <iostream>
#include "SingularityJudger.h"
#include "fstream"
#include "QString"
#include "QList"
#include "QRegularExpression"
#include "gauss_map.h"
#include "bvh.h"
#include "utils.h"
#include "unordered_set"
#include "unordered_map"

array2<glm::vec3> getPtsFromStr(QString srfData);
void printCtrPtsAsQuadruples(const RationalSurface<float> &s);
void printCtrPtsAsQuadruples(const Surface<double> &s);

void sampleTimeTest(const RationalSurface<float> &s_, int sampleLevel) {
    // 由于ParaSolid那边曲面参数为double类型，这里也要保证是double类型（控制变量）
    RationalSurface<double> s;
    vector<double> knots_u(s_.knots_u.size());
    vector<double> knots_v(s_.knots_v.size());
    array2<double> weights(s_.weights.rows(), s_.weights.cols());
    array2<glm::vec<3, double>> control_points(s_.control_points.rows(), s_.control_points.cols());

    for(int i = 0; i < s_.knots_u.size(); i++) knots_u[i] = (double)s_.knots_u[i];
    for(int i = 0; i < s_.knots_v.size(); i++) knots_v[i] = (double)s_.knots_v[i];
    for(int i = 0; i < s_.weights.rows(); i++) {
        for(int j = 0; j < s_.weights.cols(); j++) {
            weights(i, j) = (double)s_.weights(i, j);
        }
    }
    for(int i = 0; i < s_.control_points.rows(); i++) {
        for(int j = 0; j < s_.control_points.cols(); j++) {
            control_points(i, j).x = (double)s_.control_points(i, j).x;
            control_points(i, j).y = (double)s_.control_points(i, j).y;
            control_points(i, j).z = (double)s_.control_points(i, j).z;
        }
    }

    s.knots_u = knots_u;
    s.knots_v = knots_v;
    s.weights = weights;
    s.control_points = control_points;
    s.degree_u = s_.degree_u;
    s.degree_v = s_.degree_v;
    auto sampleCnt = int(pow(2, sampleLevel - 1) + 1);
    auto s_first_u = *(s.knots_u.begin());
    auto s_first_v = *(s.knots_v.begin());
    auto s_step_u = (*(s.knots_u.end() - 1) - s_first_u) / double (sampleCnt - 1);
    auto s_step_v = (*(s.knots_v.end() - 1) - s_first_v) / double (sampleCnt - 1);
    // 为了分别测试赋值和求梯度的时间，这里把它们分开写了

    auto startMomEval = std::chrono::steady_clock::now();
    for (int i = 0; i < sampleCnt; i++) {
        auto u = s_first_u + s_step_u * float(i);
        for (int j = 0; j < sampleCnt; j++) {
            auto v = s_first_v + s_step_v * float(j);
            auto eval = tinynurbs::surfacePoint(s, u, v);
//            printf("(%d, %d) --> (%g, %g, %g)\n", i, j, eval.x, eval.y, eval.z);
        }
    }
    auto endMomEval = std::chrono::steady_clock::now();
    printf("time cost of evaluation: %lf ms\n", std::chrono::duration<double, std::milli>(endMomEval - startMomEval).count());

    auto startMomDer = std::chrono::steady_clock::now();
    for (int i = 0; i < sampleCnt; i++) {
        auto u = s_first_u + s_step_u * float(i);
        for (int j = 0; j < sampleCnt; j++) {
            auto v = s_first_v + s_step_v * float(j);
            auto der = tinynurbs::surfaceDerivatives(s, 1, u, v);
//            if(der(0, 0) == res.evaluation[i][j]) cout<<"amazing"<<endl;
//            else cout<<"what??? ("<<res.evaluation[i][j].x<<" "<<res.evaluation[i][j].y<<" "<<res.evaluation[i][j].z<<") | ("<<der(0, 0).x<<" "<<der(0, 0).y<<" "<<der(0, 0).z<<")"<<endl;
//            res.tangent_u[i][j] = der(1, 0);
//            res.tangent_v[i][j] = der(0, 1);
        }
    }
    auto endMomDer = std::chrono::steady_clock::now();
    printf("time cost of derivatives: %lf ms\n", std::chrono::duration<double, std::milli>(endMomDer - startMomDer).count());

   auto startMomScdDer = std::chrono::steady_clock::now();
    for (int i = 0; i < sampleCnt; i++) {
        auto u = s_first_u + s_step_u * float(i);
        for (int j = 0; j < sampleCnt; j++) {
            auto v = s_first_v + s_step_v * float(j);
            auto der = tinynurbs::surfaceDerivatives(s, 2, u, v);
        }
    }
    auto endMomScdDer = std::chrono::steady_clock::now();
    printf("time cost of second derivatives: %lf ms\n", std::chrono::duration<double, std::milli>(endMomScdDer - startMomScdDer).count());

}

void sampleTimeTestNonRational(const RationalSurface<float> &s_, int sampleLevel) {
    // 由于ParaSolid那边曲面参数为double类型，这里也要保证是double类型（控制变量）
    Surface<double> s;
    vector<double> knots_u(s_.knots_u.size());
    vector<double> knots_v(s_.knots_v.size());
    array2<glm::vec<3, double>> control_points(s_.control_points.rows(), s_.control_points.cols());

    for(int i = 0; i < s_.knots_u.size(); i++) knots_u[i] = (double)s_.knots_u[i];
    for(int i = 0; i < s_.knots_v.size(); i++) knots_v[i] = (double)s_.knots_v[i];
    for(int i = 0; i < s_.control_points.rows(); i++) {
        for(int j = 0; j < s_.control_points.cols(); j++) {
            control_points(i, j).x = (double)s_.control_points(i, j).x;
            control_points(i, j).y = (double)s_.control_points(i, j).y;
            control_points(i, j).z = (double)s_.control_points(i, j).z;
        }
    }

    s.knots_u = knots_u;
    s.knots_v = knots_v;
    s.control_points = control_points;
    s.degree_u = s_.degree_u;
    s.degree_v = s_.degree_v;
    auto sampleCnt = int(pow(2, sampleLevel - 1) + 1);
    auto s_first_u = *(s.knots_u.begin());
    auto s_first_v = *(s.knots_v.begin());
    auto s_step_u = (*(s.knots_u.end() - 1) - s_first_u) / double (sampleCnt - 1);
    auto s_step_v = (*(s.knots_v.end() - 1) - s_first_v) / double (sampleCnt - 1);
    // 为了分别测试赋值和求梯度的时间，这里把它们分开写了

    printCtrPtsAsQuadruples(s);

    auto startMomEval = std::chrono::steady_clock::now();
    for (int i = 0; i < sampleCnt; i++) {
        auto u = s_first_u + s_step_u * float(i);
        for (int j = 0; j < sampleCnt; j++) {
            auto v = s_first_v + s_step_v * float(j);
            auto eval = tinynurbs::surfacePoint(s, u, v);
//            res.evaluation[i][j] = tinynurbs::surfacePoint(s, u, v);
        }
    }
    auto endMomEval = std::chrono::steady_clock::now();
    printf("time cost of evaluation: %lf ms\n", std::chrono::duration<double, std::milli>(endMomEval - startMomEval).count());

    auto startMomDer = std::chrono::steady_clock::now();
    for (int i = 0; i < sampleCnt; i++) {
        auto u = s_first_u + s_step_u * float(i);
        for (int j = 0; j < sampleCnt; j++) {
            auto v = s_first_v + s_step_v * float(j);
            auto der = tinynurbs::surfaceDerivatives(s, 1, u, v);
//            if(der(0, 0) == res.evaluation[i][j]) cout<<"amazing"<<endl;
//            else cout<<"what??? ("<<res.evaluation[i][j].x<<" "<<res.evaluation[i][j].y<<" "<<res.evaluation[i][j].z<<") | ("<<der(0, 0).x<<" "<<der(0, 0).y<<" "<<der(0, 0).z<<")"<<endl;
//            res.tangent_u[i][j] = der(1, 0);
//            res.tangent_v[i][j] = der(0, 1);
        }
    }
    auto endMomDer = std::chrono::steady_clock::now();
    printf("time cost of derivatives: %lf ms\n", std::chrono::duration<double, std::milli>(endMomDer - startMomDer).count());


    auto startMomScdDer = std::chrono::steady_clock::now();
    for (int i = 0; i < sampleCnt; i++) {
        auto u = s_first_u + s_step_u * float(i);
        for (int j = 0; j < sampleCnt; j++) {
            auto v = s_first_v + s_step_v * float(j);
            auto der = tinynurbs::surfaceDerivatives(s, 2, u, v);
        }
    }
    auto endMomScdDer = std::chrono::steady_clock::now();
    printf("time cost of second derivatives: %lf ms\n", std::chrono::duration<double, std::milli>(endMomScdDer - startMomScdDer).count());

}

int dirs[4][2] = {{-1, 0},
                  {0,  1},
                  {1,  0},
                  {0,  -1}};

void dfs(const map<pair<int, int>, set<pair<int, int>>>& pairMap,
//         unordered_map<pair<int, int>, char> &book,
         set<pair<int, int>>& book,
//         vector<vector<pair<pair<int, int>, pair<int, int>>>> &boxGroups,
         vector<vector<pair<int, int>>> &boxGroups,
         int x, int y, int iOfGroup) {
//    book[{x, y}] = 1;
//    book.insert(pair<pair<int, int>, char>(pair<int, int>(x, y), 1));
    book.insert({x,y});
    boxGroups[iOfGroup].emplace_back(pair<int, int>(x, y));
    for (auto dir: dirs) {
        auto nx = x + dir[0], ny = y + dir[1];
        if (book.find({nx, ny}) != book.end() || pairMap.find({nx, ny}) == pairMap.end())continue;
        dfs(pairMap, book, boxGroups, nx, ny, iOfGroup);
    }
}

int main() {
    RationalSurface<float> s;
    RationalSurface<float> f;

    int level = 9;
    printf("level: %d, sample cnt: %d * %d\n", level, int(pow(2, level - 1)), int(pow(2, level - 1)));

    ifstream fin;
    fin.open(R"(intersectTest\casea1\surfaces.txt)");

    string str;
    string tmp;
    while (fin.good()) {
        getline(fin, tmp);
        str += tmp + "\n";
    }
//    cout << str << endl;

    auto infos = QString(str.c_str()).split(";");
    auto wData = infos[0];
    auto srf1Data = infos[1];
    auto srf2Data = infos[2];

    auto points1 = getPtsFromStr(srf1Data);
    auto points2 = getPtsFromStr(srf2Data);

    wData = wData.remove(0, wData.indexOf("Matrix") + 6).trimmed();
    wData.remove(0, 3);
    wData.remove(wData.size() - 3, 3);
    auto wInfos = wData.split(QRegularExpression("\\]( )*,( )*\\["));

    array2<float> weights;
    vector<float> wArray(points1.cols() * points1.rows());
    for (int i = 0; i < points1.cols(); i++) {
        auto wsInV = wInfos[i].split(QRegularExpression("( )*,( )*"));
        for (int j = 0; j < points1.rows(); j++) {
            wArray[i * points1.rows() + j] = wsInV[j].toFloat();
        }
    }
    weights = {points1.cols(), points2.rows(), wArray};

    // knot
    std::vector<float> knot_u(2 * points1.cols(), 1.);
    for (int i = 0; i < knot_u.size() / 2; i++)knot_u[i] = 0;
    std::vector<float> knot_v(2 * points1.rows(), 1.);
    for (int i = 0; i < knot_v.size() / 2; i++)knot_v[i] = 0;

    s.degree_u = knot_u.size() / 2 - 1;
    s.degree_v = knot_v.size() / 2 - 1;
    s.knots_u = knot_u;
    s.knots_v = knot_v;
    s.control_points = points1;
    s.weights = weights;

    f.degree_u = knot_u.size() / 2 - 1;
    f.degree_v = knot_v.size() / 2 - 1;
    f.knots_u = knot_u;
    f.knots_v = knot_v;
    f.control_points = points2;
    f.weights = weights;

    fin.close();

    printCtrPtsAsQuadruples(f);

    // ====================== 测试 =======================
    vector<vector<glm::vec3>> s_evaluation;
    vector<vector<glm::vec3>> f_evaluation;
    // 曲面s和f的切向量。zd*-sznmj
    vector<vector<glm::vec3>> s_tangent_v;
    vector<vector<glm::vec3>> f_tangent_u;
    const vector<vector<glm::vec3>> f_tangent_v;
    // 曲面s和f的法向量
    const vector<vector<glm::vec3>> s_normal;
    const vector<vector<glm::vec3>> f_normal;

//    sampleTimeTestNonRational(s, level);
    sampleTimeTest(s, level);

    auto mesh1 = SrfMesh(s, level);
    auto mesh2 = SrfMesh(f, level);

    BVH bvh1(mesh1.evaluation);
    BVH bvh2(mesh2.evaluation);
    bvh1.build();
    bvh2.build();
    auto intersectBoxPairs = getOverlapLeafNodes(bvh1, bvh2); // [{{u1, v1}, {u2, v2}}]
    printf("box pairs size: %lld\n", intersectBoxPairs.size());

    /**
     * 测试对bvh结果用dfs
     */
//    double timeClassify = utils::get_time();
//    map<pair<int, int>, set<pair<int, int>>> pairMap;
//    for (const auto &boxPair: intersectBoxPairs) {
//        pairMap[boxPair.first].insert(boxPair.second);
//    }
////    vector<vector<pair<pair<int, int>, pair<int, int>>>> groups;
//    vector<vector<pair<int, int>>> groups;
////    unordered_map<pair<int, int>, char> book;
//    set<pair<int, int>>book;
//    int groupNum = 0;
//    for(auto boxPair: intersectBoxPairs) {
//        if(book.find(boxPair.first) != book.end())continue;
//        groups.emplace_back();
//        dfs(pairMap, book, groups, boxPair.first.first, boxPair.first.second, groupNum++);
//    }
//    timeClassify = utils::get_time() - timeClassify;
//    printf("time cost to group boxes from BVH: %lf\n", timeClassify);
    /**
     * end of test
     */

    double time_cost = utils::get_time();
    SingularityJudger singularityJudger(s, f, mesh1, mesh2);

    pair<int, int> cellIdxFullRange = {0, pow(2, level - 1) - 1};
    auto cellsWithCriticalPts = singularityJudger.judge(intersectBoxPairs, cellIdxFullRange, cellIdxFullRange, cellIdxFullRange, cellIdxFullRange);

    time_cost = utils::get_time() - time_cost;
    printf("time cost of singularityJudger: %lf\n", time_cost);

    // ================== 测试对整个曲面，gauss能排除多少（或保留多少）==================
//    GaussMap gaussMap1(mesh1.normal);
//    GaussMap gaussMap2(mesh2.normal);
//    gaussMap1.build();
//    gaussMap2.build();
//    auto pairs = getOverlapLeafNodes(gaussMap1, gaussMap2);
//    printf("Gauss Map: keep %lld samples in totally %lld boxes\n", pairs.size(),
//           mesh1.normal.size() * mesh1.normal[0].size() * mesh2.normal.size() * mesh2.normal[0].size());

    return 0;
}

array2<glm::vec3> getPtsFromStr(QString srfData) {
    srfData = srfData.remove(0, srfData.indexOf("Matrix") + 6).trimmed();
    // 去掉首尾的括号
    srfData.remove(0, 1);
    srfData.remove(srfData.size() - 1, 1);
    auto srfInfos = srfData.split("{");
    auto srfDimInfos = srfInfos[0].split(QRegularExpression(",( )*"));
    auto dimU = srfDimInfos[0].toInt();
    auto dimV = srfDimInfos[1].toInt();
    auto srfPtInfo = srfInfos[1].remove(srfInfos[1].indexOf("}"), 1);
    auto srfPtsInfos = srfPtInfo.split(QRegularExpression(",( )*\\("));
    std::vector<glm::vec3> points(dimU * dimV, glm::vec3());
    for (int i = 0; i < srfPtsInfos.size(); i++) {
        auto pt = srfPtsInfos[i].split("=")[1].trimmed();
        // 去掉首尾的方括号
        pt.remove(0, 1);
        pt.remove(pt.size() - 1, 1);
//        int iu = i / dimU;
//        int iv = i % dimU;
        auto coords = pt.split(QRegularExpression(",( )*"));
        for (int k = 0; k < 3; k++) {
            float num;
            if (coords[k].contains("/")) {
                auto nums = coords[k].split("/");
                num = nums[0].trimmed().toFloat() / nums[1].trimmed().toFloat();
            } else {
                num = coords[k].toFloat();
            }
            if (k == 0) points[i].x = num;
            else if (k == 1) points[i].y = num;
            else points[i].z = num;
        }
    }
    array2<glm::vec3> res = {(size_t) dimU, (size_t) dimV, points};
    return res;
}

void printCtrPtsAsQuadruples(const RationalSurface<float> &s) {
    cout<<endl;
    for(int i = 0; i < s.control_points.rows(); i++) {
        for(int j = 0; j < s.control_points.cols(); j++) {
            auto pt = s.control_points(i, j);
            auto w = s.weights(i, j);
            cout<<pt.x * w<<", "<<pt.y * w<<", "<<pt.z * w<<", "<<s.weights(i, j)<<", "<<endl;
        }
    }
    cout<<s.control_points.rows() * s.control_points.cols() * 4 << endl;
}

void printCtrPtsAsQuadruples(const Surface<double> &s) {
    cout<<endl;
    for(int i = 0; i < s.control_points.rows(); i++) {
        for(int j = 0; j < s.control_points.cols(); j++) {
            auto pt = s.control_points(i, j);
            cout<<pt.x<<", "<<pt.y<<", "<<pt.z<<", "<<endl;
        }
    }
    cout<<s.control_points.rows() * s.control_points.cols() * 4 << endl;
}