Commit 0d05dc44 authored by Teseo Schneider's avatar Teseo Schneider
Browse files

Merge branch 'master' of github.com:wildmeshing/fTetWild

parents cc6df16d 7668be9d
Loading
Loading
Loading
Loading
+19 −0
Original line number Diff line number Diff line
@@ -155,6 +155,25 @@ namespace floatTetWild {
            return false;
        }

        inline bool is_out_sf_envelope(const Vector3& p, const Scalar eps_2,
                                       GEO::index_t& prev_facet, double& sq_dist, GEO::vec3& nearest_p) const {
            GEO::vec3 geo_p(p[0], p[1], p[2]);
            return is_out_sf_envelope(geo_p, eps_2, prev_facet, sq_dist, nearest_p);
        }
        inline bool is_out_sf_envelope(const GEO::vec3& geo_p, const Scalar eps_2,
                                       GEO::index_t& prev_facet, double& sq_dist, GEO::vec3& nearest_p) const {
            if (prev_facet != GEO::NO_FACET) {
                get_point_facet_nearest_point(sf_mesh, geo_p, prev_facet, nearest_p, sq_dist);
            }
            if (Scalar(sq_dist) > eps_2) {
                sf_tree.facet_in_envelope_with_hint(geo_p, eps_2, prev_facet, nearest_p, sq_dist);
            }

            if (Scalar(sq_dist) > eps_2)
                return true;
            return false;
        }

        inline bool is_out_b_envelope(const Vector3 &p, const Scalar eps_2, GEO::index_t &prev_facet) const {
            GEO::vec3 nearest_p;
            double sq_dist;
+113 −0
Original line number Diff line number Diff line
@@ -618,9 +618,13 @@ bool floatTetWild::is_out_envelope(const Mesh& mesh, int v_id, const Vector3& ne
                        vs[k] = mesh.tet_vertices[mesh.tets[t_id][mod4(j + 1 + k)]].pos;
                }

#ifdef STORE_SAMPLE_POINTS
                ps.clear();
                sample_triangle(vs, ps, mesh.params.dd);
                bool is_out = tree.is_out_sf_envelope(ps, mesh.params.eps_2, prev_facet);
#else
                bool is_out = sample_triangle_and_check_is_out(vs, mesh.params.dd, mesh.params.eps_2, tree, prev_facet);
#endif
                if(!mesh.params.envelope_log.empty() && envelope_log_csv_cnt < 1e5){
                    std::ostringstream ss;
                    ss<<std::setprecision(17);
@@ -788,6 +792,115 @@ void floatTetWild::sample_triangle(const std::array<Vector3, 3>& vs, std::vector
    }
}

bool floatTetWild::sample_triangle_and_check_is_out(const std::array<Vector3, 3>& vs, Scalar sampling_dist,
        Scalar eps_2, const AABBWrapper& tree, GEO::index_t& prev_facet){
    GEO::vec3 nearest_point;
    double sq_dist = std::numeric_limits<double>::max();

    Scalar sqrt3_2 = std::sqrt(3) / 2;

    std::array<Scalar, 3> ls;
    for (int i = 0; i < 3; i++) {
        ls[i] = (vs[i] - vs[mod3(i + 1)]).squaredNorm();
    }
    auto min_max = std::minmax_element(ls.begin(), ls.end());
    int min_i = min_max.first - ls.begin();
    int max_i = min_max.second - ls.begin();
    Scalar N = sqrt(ls[max_i]) / sampling_dist;
    if (N <= 1) {
        for (int i = 0; i < 3; i++) {
//            ps.push_back(GEO::vec3(vs[i][0], vs[i][1], vs[i][2]));
            if (tree.is_out_sf_envelope(vs[i], eps_2, prev_facet, sq_dist, nearest_point))
                return true;
        }
//        return;
        return false;
    }
    if (N == int(N))
        N -= 1;

    GEO::vec3 v0(vs[max_i][0], vs[max_i][1], vs[max_i][2]);
    GEO::vec3 v1(vs[mod3(max_i + 1)][0], vs[mod3(max_i + 1)][1], vs[mod3(max_i + 1)][2]);
    GEO::vec3 v2(vs[mod3(max_i + 2)][0], vs[mod3(max_i + 2)][1], vs[mod3(max_i + 2)][2]);

    GEO::vec3 n_v0v1 = GEO::normalize(v1 - v0);
    for (int n = 0; n <= N; n++) {
//        ps.push_back(v0 + n_v0v1 * sampling_dist * n);
        if (tree.is_out_sf_envelope(v0 + n_v0v1 * sampling_dist * n, eps_2, prev_facet, sq_dist, nearest_point))
            return true;
    }
//    ps.push_back(v1);
    if (tree.is_out_sf_envelope(v1, eps_2, prev_facet, sq_dist, nearest_point))
        return true;

    Scalar h = GEO::distance(GEO::dot((v2 - v0), (v1 - v0)) * (v1 - v0) / ls[max_i] + v0, v2);
    int M = h / (sqrt3_2 * sampling_dist);
    if (M < 1) {
//        ps.push_back(v2);
//        return;
        return tree.is_out_sf_envelope(v2, eps_2, prev_facet, sq_dist, nearest_point);
    }

    GEO::vec3 n_v0v2 = GEO::normalize(v2 - v0);
    GEO::vec3 n_v1v2 = GEO::normalize(v2 - v1);
    Scalar tan_v0, tan_v1, sin_v0, sin_v1;
    sin_v0 = GEO::length(GEO::cross((v2 - v0), (v1 - v0))) / (GEO::distance(v0, v2) * GEO::distance(v0, v1));
    tan_v0 = GEO::length(GEO::cross((v2 - v0), (v1 - v0))) / GEO::dot((v2 - v0), (v1 - v0));
    tan_v1 = GEO::length(GEO::cross((v2 - v1), (v0 - v1))) / GEO::dot((v2 - v1), (v0 - v1));
    sin_v1 = GEO::length(GEO::cross((v2 - v1), (v0 - v1))) / (GEO::distance(v1, v2) * GEO::distance(v0, v1));

    for (int m = 1; m <= M; m++) {
        int n = sqrt3_2 / tan_v0 * m + 0.5;
        int n1 = sqrt3_2 / tan_v0 * m;
        if (m % 2 == 0 && n == n1) {
            n += 1;
        }
        GEO::vec3 v0_m = v0 + m * sqrt3_2 * sampling_dist / sin_v0 * n_v0v2;
        GEO::vec3 v1_m = v1 + m * sqrt3_2 * sampling_dist / sin_v1 * n_v1v2;
        if (GEO::distance(v0_m, v1_m) <= sampling_dist)
            break;

        Scalar delta_d = ((n + (m % 2) / 2.0) - m * sqrt3_2 / tan_v0) * sampling_dist;
        GEO::vec3 v = v0_m + delta_d * n_v0v1;
        int N1 = GEO::distance(v, v1_m) / sampling_dist;
        for (int i = 0; i <= N1; i++) {
//            ps.push_back(v + i * n_v0v1 * sampling_dist);
            if (tree.is_out_sf_envelope(v + i * n_v0v1 * sampling_dist, eps_2, prev_facet, sq_dist, nearest_point))
                return true;
        }
    }
//    ps.push_back(v2);
    if (tree.is_out_sf_envelope(v2, eps_2, prev_facet, sq_dist, nearest_point))
        return true;

    //sample edges
    N = sqrt(ls[mod3(max_i + 1)]) / sampling_dist;
    if (N > 1) {
        if (N == int(N))
            N -= 1;
        GEO::vec3 n_v1v2 = GEO::normalize(v2 - v1);
        for (int n = 1; n <= N; n++) {
//            ps.push_back(v1 + n_v1v2 * sampling_dist * n);
            if (tree.is_out_sf_envelope(v1 + n_v1v2 * sampling_dist * n, eps_2, prev_facet, sq_dist, nearest_point))
                return true;
        }
    }

    N = sqrt(ls[mod3(max_i + 2)]) / sampling_dist;
    if (N > 1) {
        if (N == int(N))
            N -= 1;
        GEO::vec3 n_v2v0 = GEO::normalize(v0 - v2);
        for (int n = 1; n <= N; n++) {
//            ps.push_back(v2 + n_v2v0 * sampling_dist * n);
            if (tree.is_out_sf_envelope(v2 + n_v2v0 * sampling_dist * n, eps_2, prev_facet, sq_dist, nearest_point))
                return true;
        }
    }

    return false;
}

void floatTetWild::get_new_tet_slots(Mesh& mesh, int n, std::vector<int>& new_conn_tets) {
    int cnt = 0;
    for (int i = mesh.t_empty_start; i < mesh.tets.size(); i++) {
+2 −0
Original line number Diff line number Diff line
@@ -56,6 +56,8 @@ namespace floatTetWild {
    bool is_out_envelope(const Mesh& mesh, int v_id, const Vector3& new_pos, const AABBWrapper& tree);
    bool is_out_boundary_envelope(const Mesh& mesh, int v_id, const Vector3& new_pos, const AABBWrapper& tree);
    void sample_triangle(const std::array<Vector3, 3>& vs, std::vector<GEO::vec3>& ps, Scalar sampling_dist);
    bool sample_triangle_and_check_is_out(const std::array<Vector3, 3>& vs, Scalar sampling_dist,
            Scalar eps_2, const AABBWrapper& tree, GEO::index_t& prev_facet);

    bool is_bbox_edge(const Mesh& mesh, int v1_id, int v2_id, const std::vector<int>& n12_t_ids);
    bool is_surface_edge(const Mesh& mesh, int v1_id, int v2_id, const std::vector<int>& n12_t_ids);
+4 −0
Original line number Diff line number Diff line
@@ -882,9 +882,13 @@ floatTetWild::Scalar floatTetWild::get_angle_cos(const Vector3& p, const Vector3
}

bool floatTetWild::is_out_envelope(const std::array<Vector3, 3>& vs, const AABBWrapper& tree, const Parameters& params) {
#ifdef STORE_SAMPLE_POINTS
    std::vector<GEO::vec3> ps;
    sample_triangle(vs, ps, params.dd_simplification);
    return tree.is_out_sf_envelope(ps, params.eps_2_simplification);
#else
    return sample_triangle_and_check_is_out(vs, params.dd_simplification, params.eps_2_simplification, tree, GEO::NO_FACET);
#endif

    // GEO::vec3 init_point(vs[0][0], vs[0][1], vs[0][2]);
    // GEO::vec3 nearest_point;
+12 −3
Original line number Diff line number Diff line
@@ -2001,10 +2001,15 @@ bool floatTetWild::insert_boundary_edges(const std::vector<Vector3> &input_verti
                            new_tets, new_track_surface_fs, modified_t_ids)) {
            bool is_inside_envelope = true;
            for (auto &f: cut_fs) {
#ifdef STORE_SAMPLE_POINTS
                std::vector<GEO::vec3> ps;
                sample_triangle({{mesh.tet_vertices[f[0]].pos, mesh.tet_vertices[f[1]].pos,
                                         mesh.tet_vertices[f[2]].pos}}, ps, mesh.params.dd);
                if (tree.is_out_sf_envelope(ps, mesh.params.eps_2)) {
#else
                if(sample_triangle_and_check_is_out({{mesh.tet_vertices[f[0]].pos, mesh.tet_vertices[f[1]].pos,
                                         mesh.tet_vertices[f[2]].pos}}, mesh.params.dd, mesh.params.eps_2, tree, GEO::NO_FACET)){
#endif
                    is_inside_envelope = false;
                    break;
                }
@@ -2583,12 +2588,16 @@ void floatTetWild::mark_surface_fs(const std::vector<Vector3> &input_vertices, c
//                if(is_all_out)
//                    continue;
//                //
                double eps = (mesh.params.eps + mesh.params.eps_simplification) / 2;
                double eps_2 = (mesh.params.eps + mesh.params.eps_simplification) / 2;
                double dd = (mesh.params.dd + mesh.params.dd_simplification) / 2;
                eps *= eps;
                eps_2 *= eps_2;
#ifdef STORE_SAMPLE_POINTS
                std::vector<GEO::vec3> ps;
                sample_triangle({{tp1_3d, tp2_3d, tp3_3d}}, ps, dd);
                if (tree.is_out_sf_envelope(ps, eps))
                if (tree.is_out_sf_envelope(ps, eps_2))
#else
                if(sample_triangle_and_check_is_out({{tp1_3d, tp2_3d, tp3_3d}}, dd, eps_2, tree, GEO::NO_FACET))
#endif
                    continue;
                else
                    ff_id = track_surface_fs[t_id][j].front();
Loading