Commit 465544ba authored by YIxin-Hu's avatar YIxin-Hu
Browse files

fix insert_boundary_edges_get_intersecting_edges_and_points()

compute intersection per edge rather than per face
parent 2070a0bb
Loading
Loading
Loading
Loading
+3 −1
Original line number Diff line number Diff line
@@ -831,7 +831,9 @@ void floatTetWild::set_intersection_sorted(const std::vector<int>& s1, const std
    v.resize(it - v.begin());
}

void floatTetWild::pausee() {
void floatTetWild::pausee(std::string msg) {
    if (!msg.empty())
        cout << msg << endl;
    cout << "Is pausing... (Enter '0' to exit and other characters to continue.)" << endl;
    char c;
    std::cin >> c;
+1 −1
Original line number Diff line number Diff line
@@ -107,7 +107,7 @@ namespace floatTetWild {
        return j%2;
    }

    void pausee();
    void pausee(std::string msg = "");

    ///////////////
    class ElementInQueue{
+232 −132
Original line number Diff line number Diff line
@@ -1430,6 +1430,8 @@ bool floatTetWild::insert_boundary_edges(const std::vector<Vector3> &input_verti

    std::vector<std::vector<std::pair<int, int>>> covered_fs_infos(input_faces.size());
    for (int i = 0; i < track_surface_fs.size(); i++) {
        if (mesh.tets[i].is_removed)
            continue;
        for (int j = 0; j < 4; j++) {
            for (int f_id: track_surface_fs[i][j])
                covered_fs_infos[f_id].push_back(std::make_pair(i, j));
@@ -1538,23 +1540,23 @@ bool floatTetWild::insert_boundary_edges(const std::vector<Vector3> &input_verti

    logger().info("uninsert boundary #e = {}/{}", b_edge_infos.size() - cnt, b_edge_infos.size());

    //fortest
    std::ofstream fout_e("b_edges1.obj");
//    std::ofstream fout_v("b_edges1_points.xyz");
    for (auto &e: b_edges) {
        fout_e << "v " << input_vertices[e[0]][0] << " "
               << input_vertices[e[0]][1] << " "
               << input_vertices[e[0]][2] << endl;
        fout_e << "v " << input_vertices[e[1]][0] << " "
               << input_vertices[e[1]][1] << " "
               << input_vertices[e[1]][2] << endl;
    }
    for (int i = 0; i < b_edges.size(); i++) {
        fout_e << "l " << i * 2 + 1 << " " << i * 2 + 2 << endl;
    }
    fout_e.close();
//    fout_v.close();
    //fortest
//    //fortest
//    std::ofstream fout_e("b_edges1.obj");
////    std::ofstream fout_v("b_edges1_points.xyz");
//    for (auto &e: b_edges) {
//        fout_e << "v " << input_vertices[e[0]][0] << " "
//               << input_vertices[e[0]][1] << " "
//               << input_vertices[e[0]][2] << endl;
//        fout_e << "v " << input_vertices[e[1]][0] << " "
//               << input_vertices[e[1]][1] << " "
//               << input_vertices[e[1]][2] << endl;
//    }
//    for (int i = 0; i < b_edges.size(); i++) {
//        fout_e << "l " << i * 2 + 1 << " " << i * 2 + 2 << endl;
//    }
//    fout_e.close();
////    fout_v.close();
//    //fortest

    cout << "b_edges.size = " << b_edges.size() << endl;

@@ -1581,7 +1583,9 @@ bool floatTetWild::insert_boundary_edges_get_intersecting_edges_and_points(
                  input_vertices[input_faces[n_f_ids.front()][1]],
                  input_vertices[input_faces[n_f_ids.front()][2]]);
    std::array<Vector2, 2> evs_2d = {{to_2d(input_vertices[e[0]], t), to_2d(input_vertices[e[1]], t)}};
    Vector3 n = (input_vertices[input_faces[n_f_ids.front()][1]] - input_vertices[input_faces[n_f_ids.front()][0]]).cross(input_vertices[input_faces[n_f_ids.front()][2]] - input_vertices[input_faces[n_f_ids.front()][0]]);
    Vector3 n = (input_vertices[input_faces[n_f_ids.front()][1]] -
                 input_vertices[input_faces[n_f_ids.front()][0]]).cross(
            input_vertices[input_faces[n_f_ids.front()][2]] - input_vertices[input_faces[n_f_ids.front()][0]]);
    n.normalize();
    const Vector3 &pp = input_vertices[input_faces[n_f_ids.front()][0]];

@@ -1631,7 +1635,8 @@ bool floatTetWild::insert_boundary_edges_get_intersecting_edges_and_points(
        }
    }

    std::vector<std::array<int, 3>> f_oris;
//    std::vector<std::array<int, 3>> f_oris;
    std::vector<int> v_oris(mesh.tet_vertices.size(), Predicates::ORI_UNKNOWN);
    while (!t_ids_queue.empty()) {
        int t_id = t_ids_queue.front();
        t_ids_queue.pop();
@@ -1661,20 +1666,40 @@ bool floatTetWild::insert_boundary_edges_get_intersecting_edges_and_points(
            int cnt_pos = 0;
            int cnt_neg = 0;
            int cnt_on = 0;
            std::array<int, 3> oris;
//            std::array<int, 3> oris;
//            for (int k = 0; k < 3; k++) {
//                oris[k] = Predicates::orient_2d(evs_2d[0], evs_2d[1], fvs_2d[k]);
//                if (oris[k] == Predicates::ORI_ZERO) {
//                    cnt_on++;
//                } else {
//                    Scalar dis_2 = p_seg_squared_dist_3d(mesh.tet_vertices[f_v_ids[k]].pos, input_vertices[e[0]],
//                                                         input_vertices[e[1]]);
//                    if (dis_2 < mesh.params.eps_2_coplanar) {
//                        oris[k] = Predicates::ORI_ZERO;
//                        cnt_on++;
//                        continue;
//                    }
//                    if (oris[k] == Predicates::ORI_POSITIVE)
//                        cnt_pos++;
//                    else
//                        cnt_neg++;
//                }
//            }
            for (int k = 0; k < 3; k++) {
                oris[k] = Predicates::orient_2d(evs_2d[0], evs_2d[1], fvs_2d[k]);
                if (oris[k] == Predicates::ORI_ZERO) {
                int &ori = v_oris[f_v_ids[k]];
                if (ori == Predicates::ORI_UNKNOWN)
                    ori = Predicates::orient_2d(evs_2d[0], evs_2d[1], fvs_2d[k]);
                if (ori == Predicates::ORI_ZERO) {
                    cnt_on++;
                } else {
                    Scalar dis_2 = p_seg_squared_dist_3d(mesh.tet_vertices[f_v_ids[k]].pos, input_vertices[e[0]],
                                                         input_vertices[e[1]]);
                    if (dis_2 < mesh.params.eps_2_coplanar) {
                        oris[k] = Predicates::ORI_ZERO;
                        ori = Predicates::ORI_ZERO;
                        cnt_on++;
                        continue;
                    }
                    if (oris[k] == Predicates::ORI_POSITIVE)
                    if (ori == Predicates::ORI_POSITIVE)
                        cnt_pos++;
                    else
                        cnt_neg++;
@@ -1683,7 +1708,7 @@ bool floatTetWild::insert_boundary_edges_get_intersecting_edges_and_points(
            if (cnt_on >= 2) {
                cut_fs.push_back(f_v_ids);
                std::sort(cut_fs.back().begin(), cut_fs.back().end());
                f_oris.push_back(oris);
//                f_oris.push_back(oris);
                continue;
            }
            if (cnt_neg == 0 || cnt_pos == 0)
@@ -1700,7 +1725,10 @@ bool floatTetWild::insert_boundary_edges_get_intersecting_edges_and_points(
            if (!is_intersected) { ///then check if there's intersection
                for (int k = 0; k < 3; k++) {
                    //if cross
                    if (!is_cross(oris[k], oris[(k + 1) % 3]))
//                    if (!is_cross(oris[k], oris[(k + 1) % 3]))
                    if(v_oris[f_v_ids[k]] == Predicates::ORI_UNKNOWN || v_oris[f_v_ids[(k + 1) % 3]] == Predicates::ORI_UNKNOWN)
                        pausee("hehe");
                    if (!is_cross(v_oris[f_v_ids[k]], v_oris[f_v_ids[(k + 1) % 3]]))
                        continue;
                    //if already know intersect
                    std::array<int, 2> tri_e = {{f_v_ids[k], f_v_ids[(k + 1) % 3]}};
@@ -1719,13 +1747,13 @@ bool floatTetWild::insert_boundary_edges_get_intersecting_edges_and_points(
                        double dis2 = (p - mesh.tet_vertices[f_v_ids[(k + 1) % 3]].pos).squaredNorm();
//                        if (dis1 < SCALAR_ZERO_2) {
                        if (dis1 < mesh.params.eps_2_coplanar) {
                            oris[k] = Predicates::ORI_ZERO;
                            v_oris[f_v_ids[k]] = Predicates::ORI_ZERO;
                            is_intersected = true;
                            break;
                        }
//                        if (dis2 < SCALAR_ZERO_2) {
                        if (dis2 < mesh.params.eps_2_coplanar) {
                            oris[k] = Predicates::ORI_ZERO;
                            v_oris[f_v_ids[k]] = Predicates::ORI_ZERO;
                            is_intersected = true;
                            break;
                        }
@@ -1739,9 +1767,9 @@ bool floatTetWild::insert_boundary_edges_get_intersecting_edges_and_points(
                    continue;
            }

            std::sort(f_v_ids.begin(), f_v_ids.end());
            cut_fs.push_back(f_v_ids);
            std::sort(cut_fs.back().begin(), cut_fs.back().end());
            f_oris.push_back(oris);
//            f_oris.push_back(oris);
            is_cut_vs[(j + 1) % 4] = true;
            is_cut_vs[(j + 2) % 4] = true;
            is_cut_vs[(j + 3) % 4] = true;
@@ -1758,90 +1786,82 @@ bool floatTetWild::insert_boundary_edges_get_intersecting_edges_and_points(
        }
    }
    //remove duplicated elements
    {
        std::vector<int> indices(cut_fs.size());
        for (int i = 0; i < cut_fs.size(); i++)
            indices[i] = i;
        std::sort(indices.begin(), indices.end(), [&cut_fs](int i1, int i2) {
            return cut_fs[i1] < cut_fs[i2];
        });
        indices.erase(std::unique(indices.begin(), indices.end(), [&cut_fs](int i1, int i2) {
            return cut_fs[i1] == cut_fs[i2];
        }), indices.end());
        std::vector<std::array<int, 3>> new_cut_fs(indices.size());
        std::vector<std::array<int, 3>> new_f_oris(indices.size());
        for (int i = 0; i < indices.size(); i++) {
            new_cut_fs[i] = cut_fs[indices[i]];
            new_f_oris[i] = f_oris[indices[i]];
        }
        cut_fs = new_cut_fs;
        f_oris = new_f_oris;
    }
//    {
//        std::vector<int> indices(cut_fs.size());
//        for (int i = 0; i < cut_fs.size(); i++)
//            indices[i] = i;
//        std::sort(indices.begin(), indices.end(), [&cut_fs](int i1, int i2) {
//            return cut_fs[i1] < cut_fs[i2];
//        });
//        indices.erase(std::unique(indices.begin(), indices.end(), [&cut_fs](int i1, int i2) {
//            return cut_fs[i1] == cut_fs[i2];
//        }), indices.end());
//        std::vector<std::array<int, 3>> new_cut_fs(indices.size());
//        std::vector<std::array<int, 3>> new_f_oris(indices.size());
//        for (int i = 0; i < indices.size(); i++) {
//            new_cut_fs[i] = cut_fs[indices[i]];
//            new_f_oris[i] = f_oris[indices[i]];
//        }
//        cut_fs = new_cut_fs;
//        f_oris = new_f_oris;
//    }
    vector_unique(cut_fs);

    for (int i = 0; i < cut_fs.size(); i++) {
        std::array<bool, 3> is_e_intersected = {{false, false, false}};
        int cnt = 0;
    std::vector<std::array<int, 2>> tet_edges;
    for (const auto &f:cut_fs) {
        for (int j = 0; j < 3; j++) {
            if (f_oris[i][j] == Predicates::ORI_ZERO) {
                cnt++;
                is_e_intersected[j] = true;
                snapped_v_ids.push_back(cut_fs[i][j]);
                continue;
            if (f[j] < f[(j + 1) % 3])
                tet_edges.push_back({{f[j], f[(j + 1) % 3]}});
            else
                tet_edges.push_back({{f[(j + 1) % 3], f[j]}});
        }
            std::array<int, 2> tri_e = {{cut_fs[i][j], cut_fs[i][(j + 1) % 3]}};
            if (tri_e[0] > tri_e[1])
                std::swap(tri_e[0], tri_e[1]);
            if (map_edge_to_intersecting_point.find(tri_e) != map_edge_to_intersecting_point.end()) {
                cnt++;
                is_e_intersected[j] = true;
                continue;
    }

            if (f_oris[i][(j + 1) % 3] == Predicates::ORI_ZERO)
                is_e_intersected[j] = true;
    vector_unique(tet_edges);
    //
    for (const auto &tet_e:tet_edges) {
        bool is_snapped = false;
        for (int j = 0; j < 2; j++) {
            if (v_oris[tet_e[j]] == Predicates::ORI_ZERO) {
                snapped_v_ids.push_back(tet_e[j]);
                is_snapped = true;
            }
//        if (cnt == 2)
        if (cnt >= 2)
            continue;

        //line - tri edges intersection
        for (int j = 0; j < 3; j++) {
            if (is_e_intersected[j])
        }
        if (is_snapped)
            continue;
            if (!is_cross(f_oris[i][j], f_oris[i][(j + 1) % 3]))
        if (map_edge_to_intersecting_point.find(tet_e) != map_edge_to_intersecting_point.end())
            continue;

//            std::array<Vector2, 2> tri_evs_2d = {{to_2d(mesh.tet_vertices[cut_fs[i][j]].pos, t),
//                                                         to_2d(mesh.tet_vertices[cut_fs[i][(j + 1) % 3]].pos, t)}};
            std::array<Vector2, 2> tri_evs_2d = {{to_2d(mesh.tet_vertices[cut_fs[i][j]].pos, n, pp, t),
                                                         to_2d(mesh.tet_vertices[cut_fs[i][(j + 1) % 3]].pos, n, pp,
                                                               t)}};

        std::array<Vector2, 2> tri_evs_2d = {{to_2d(mesh.tet_vertices[tet_e[0]].pos, n, pp, t),
                                                     to_2d(mesh.tet_vertices[tet_e[1]].pos, n, pp, t)}};
        Scalar t_seg = -1;
        if (seg_line_intersection_2d(tri_evs_2d, evs_2d, t_seg)) {
                std::array<int, 2> tri_e = {{cut_fs[i][j], cut_fs[i][(j + 1) % 3]}};
                if (tri_e[0] > tri_e[1])
                    std::swap(tri_e[0], tri_e[1]);
                points.push_back((1 - t_seg) * mesh.tet_vertices[cut_fs[i][j]].pos
                                 + t_seg * mesh.tet_vertices[cut_fs[i][(j + 1) % 3]].pos);
                map_edge_to_intersecting_point[tri_e] = points.size() - 1;
                cnt++;
            points.push_back((1 - t_seg) * mesh.tet_vertices[tet_e[0]].pos
                             + t_seg * mesh.tet_vertices[tet_e[1]].pos);
            map_edge_to_intersecting_point[tet_e] = points.size() - 1;
        }
    }
//        if (cnt != 2) {
        if (cnt < 2) {
//            cout<<"return 2, cnt = "<<cnt<<endl;
    vector_unique(snapped_v_ids);
//    cout << "map_edge_to_intersecting_point.size = " << map_edge_to_intersecting_point.size() << endl;
//    cout << "points.size = " << points.size() << endl;

//    //fortest
//    Eigen::MatrixXd V(cut_fs.size() * 3, 3), C(cut_fs.size() * 3, 3);
//    Eigen::MatrixXi F(cut_fs.size(), 3);
//    for (int i = 0; i < cut_fs.size(); i++) {
//        for (int j = 0; j < 3; j++) {
//            V.row(i * 3 + j) = mesh.tet_vertices[cut_fs[i][j]].pos;
//                    C.row(i * 3 + j) << 1, 1, 1;
//            C.row(i * 3 + j) << 0, 0, 1;
//        }
//        F.row(i) << i * 3, i * 3 + 1, i * 3 + 2;
//    }
//    igl::writeOFF("test_f.off", V, F, C);
////            Eigen::MatrixXd V(3, 3), C(3, 3);
////            Eigen::MatrixXi F(1, 3);
////            for (int j = 0; j < 3; j++) {
////                V.row(j) = mesh.tet_vertices[cut_fs[i][j]].pos;
////                C.row(j) << 0, 0, 1;
////            }
////            F.row(i) << 0, 1, 2;
////            igl::writeOFF("test_f.off", V, F, C);
//    //
//    std::ofstream fout("test_e.obj");
//    fout << "v " << input_vertices[e[0]].transpose() << endl;
@@ -1849,21 +1869,101 @@ bool floatTetWild::insert_boundary_edges_get_intersecting_edges_and_points(
//    fout << "l 1 2" << endl;
//    fout.close();
//    //
////            cout << "cnt_tmp = " << cnt_tmp << endl;
////            cout << "cnt = " << cnt << endl;
//    pausee();
//    //fortest


//    for (int i = 0; i < cut_fs.size(); i++) {
//        std::array<bool, 3> is_e_intersected = {{false, false, false}};
//        int cnt = 0;
//        for (int j = 0; j < 3; j++) {
//                cout << cut_fs[i][j] << " " << cut_fs[i][(j + 1) % 3] << ": "
//                     << seg_seg_squared_dist_3d(
//                             {{mesh.tet_vertices[cut_fs[i][j]].pos, mesh.tet_vertices[cut_fs[i][(j + 1) % 3]].pos}},
//                             {{input_vertices[e[0]], input_vertices[e[1]]}}) << endl;
//            if (f_oris[i][j] == Predicates::ORI_ZERO) {
//                cnt++;
//                is_e_intersected[j] = true;
//                snapped_v_ids.push_back(cut_fs[i][j]);
//                continue;
//            }
//            std::array<int, 2> tri_e = {{cut_fs[i][j], cut_fs[i][(j + 1) % 3]}};
//            if (tri_e[0] > tri_e[1])
//                std::swap(tri_e[0], tri_e[1]);
//            if (map_edge_to_intersecting_point.find(tri_e) != map_edge_to_intersecting_point.end()) {
//                cnt++;
//                is_e_intersected[j] = true;
//                continue;
//            }
//
//            if (f_oris[i][(j + 1) % 3] == Predicates::ORI_ZERO)
//                is_e_intersected[j] = true;
//        }
////        if (cnt == 2)
//        if (cnt >= 2)
//            continue;
//
//        //line - tri edges intersection
//        int cnt_tmp = 0;
//        for (int j = 0; j < 3; j++) {
//            if (is_e_intersected[j])
//                continue;
//            if (!is_cross(f_oris[i][j], f_oris[i][(j + 1) % 3]))
//                continue;
//
////            std::array<Vector2, 2> tri_evs_2d = {{to_2d(mesh.tet_vertices[cut_fs[i][j]].pos, t),
////                                                         to_2d(mesh.tet_vertices[cut_fs[i][(j + 1) % 3]].pos, t)}};
//            std::array<Vector2, 2> tri_evs_2d = {{to_2d(mesh.tet_vertices[cut_fs[i][j]].pos, n, pp, t),
//                                                         to_2d(mesh.tet_vertices[cut_fs[i][(j + 1) % 3]].pos, n, pp,
//                                                               t)}};
//
//            Scalar t_seg = -1;
//            if (seg_line_intersection_2d(tri_evs_2d, evs_2d, t_seg)) {
//                std::array<int, 2> tri_e = {{cut_fs[i][j], cut_fs[i][(j + 1) % 3]}};
//                if (tri_e[0] > tri_e[1])
//                    std::swap(tri_e[0], tri_e[1]);
//                points.push_back((1 - t_seg) * mesh.tet_vertices[cut_fs[i][j]].pos
//                                 + t_seg * mesh.tet_vertices[cut_fs[i][(j + 1) % 3]].pos);
//                map_edge_to_intersecting_point[tri_e] = points.size() - 1;
//                cnt++;
//                cnt_tmp++;
//            }
//        }
//        if (cnt < 2) {
//            cut_fs.erase(cut_fs.begin()+i);
//            i--;
////            //fortest
//////            Eigen::MatrixXd V(cut_fs.size() * 3, 3), C(cut_fs.size() * 3, 3);
//////            Eigen::MatrixXi F(cut_fs.size(), 3);
//////            for (int i = 0; i < cut_fs.size(); i++) {
//////                for (int j = 0; j < 3; j++) {
//////                    V.row(i * 3 + j) = mesh.tet_vertices[cut_fs[i][j]].pos;
//////                    C.row(i * 3 + j) << 0, 0, 1;
//////                }
//////                F.row(i) << i * 3, i * 3 + 1, i * 3 + 2;
//////            }
//////            igl::writeOFF("test_f.off", V, F, C);
////            Eigen::MatrixXd V(3, 3), C(3, 3);
////            Eigen::MatrixXi F(1, 3);
////            for (int j = 0; j < 3; j++) {
////                V.row(j) = mesh.tet_vertices[cut_fs[i][j]].pos;
////                C.row(j) << 0, 0, 1;
////            }
////            F.row(i) << 0, 1, 2;
////            igl::writeOFF("test_f.off", V, F, C);
////            //
////            std::ofstream fout("test_e.obj");
////            fout << "v " << input_vertices[e[0]].transpose() << endl;
////            fout << "v " << input_vertices[e[1]].transpose() << endl;
////            fout << "l 1 2" << endl;
////            fout.close();
////            //
////            cout << "cnt_tmp = " << cnt_tmp << endl;
////            cout << "cnt = " << cnt << endl;
////            pausee();
////            //fortest
////
////            return false;///has to return false here
//        }
//            cout<<"e "<<e[0]<<" "<<e[1]<<endl;
//            for (int i = 0; i < cut_fs.size(); i++) {
//                cout<<"f "<<cut_fs[i][0]<<" "<<cut_fs[i][1]<<" "<<cut_fs[i][2]<<endl;
//    }
//            pausee();

            return false;///has to return false here
        }
    }

    return true;
}