Subversion Repositories gelsvn

#include <CGLA/Mat4x4d.h>

#include <CGLA/Vec3i.h>

#include <HMesh/dual.h>

#include <Geometry/GridAlgorithm.h>

#include <Geometry/build_bbtree.h>

#include <Geometry/RGrid.h>

#include <Geometry/TrilinFilter.h>

#include <GLGraphics/ManifoldRenderer.h>

using namespace Geometry;

using namespace GLGraphics;

    queue<FaceID> face_queue;

    {

        {

            FaceID fa = w.face();

            FaceID fb = f;

            assert(fa==fb);

            verts.push_back(w.vertex());

        }

        // If there are just three we are done.

        if(verts.size() == 3) continue;

                if(verts[i] != verts[j] &&

                   (level_set_id_vertex[verts[i]] == 0 || level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[j]])

                   )

        if(min_k>=0)

            // Split faces along edge whose midpoint is closest to isovalue

                face_queue.push(f_new);

        else

            cout << "failed to triangle. Suspect NaN vertex positions!!" << endl;

        Walker w = m.walker(hid);

        if(hid< w.opp().halfedge() &&

           ls_id[w.vertex()] > 0 &&

           ls_id[w.opp().vertex()] > 0 &&

           ls_id[w.vertex()] != ls_id[w.opp().vertex()])

    bool did_work = false;

    int iter=0;

    do {

        did_work = false;

        random_shuffle(begin(hvec), end(hvec));

        for(auto hid : hvec)

        {

            Walker w =m.walker(hid);

            if(mae.delta_energy(m, hid)<=0 &&

               ls_id[w.next().vertex()] != ls_id[w.opp().next().vertex()]) {

                m.flip_edge(hid);

                did_work = true;

            }

        }

    }

    while (did_work && ++iter<1000);

    cout << "flip iter" << iter << endl;

void label_connected_components(Manifold&m, int& cur_id, const vector<VertexID>& vertices,

                                const VertexAttributeVector<double>& fun,

                                VertexAttributeVector<int>& level_id)

        vq.push(vid);

            for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_ccw())

                if(status[w.vertex()]==COOKED && same_level(fun[v], fun[w.vertex()]))

                        cout << "previously labeled vertex encountered" << endl;

            }

        ++cur_id;

}

void connect_touched(Manifold& m, const VertexAttributeVector<int>& touched)

}

void remove_shorts(Manifold& m, VertexAttributeVector<int> ls_id)

{

    double avghlen=0;

    for(auto hid : m.halfedges())

        avghlen += length(m, hid);

    avghlen /= m.no_halfedges();

    for(auto hid : m.halfedges())

        if(m.in_use(hid))

        {

            Walker w = m.walker(hid);

            if(ls_id[w.vertex()] > 0 && 	

               ls_id[w.vertex()] == ls_id[w.opp().vertex()] &&

               length(m, w.halfedge()) < 0.1 * avghlen &&

               precond_collapse_edge(m, hid))

                m.collapse_edge(hid);

        }

}

vector<VertexID> separate_edges(Manifold& m, const vector<VertexID>& vertices, int id,

                                VertexAttributeVector<double>& fun,

                                VertexAttributeVector<VertexStatus>& status,

                                VertexAttributeVector<int>& ls_id)

{

    VertexAttributeVector<int> touched(m.allocated_vertices(),0);

    vector<VertexID> new_vertices;

    // Cut all edges to vertices on other LS ID

    for(auto vid: vertices)

        if(ls_id[vid]==id)

        {

            vector<HalfEdgeID> edges_to_split;

            for(Walker w = m.walker(vid); !w.full_circle(); w = w.circulate_vertex_ccw())

            {

                int id = ls_id[w.vertex()];

                if(id > 0 & id != ls_id[vid])

                    edges_to_split.push_back(w.halfedge());

            }

            for (auto hid : edges_to_split)

            {

                VertexID vnew = m.split_edge(hid);

                fun[vnew] = fun[vid];

                m.pos(vnew) = m.pos(vid);

                status[vnew] = COOKED;

                ls_id[vnew] = ls_id[vid];

                touched[vnew] = 1;

                new_vertices.push_back(vnew);

            }

        }

    // Label old vertices as plain (uncooked)

    for(auto vid: vertices)

        if(ls_id[vid]==id)

        {

            status[vid] = RAW;

            ls_id[vid] = 0;

        }

    connect_touched(m, touched);

    return new_vertices;

}

vector<VertexID> cut_mesh(Manifold& m, VertexAttributeVector<double>& fun, double cut_val,

                          VertexAttributeVector<VertexStatus>& status)

{

    cout << "cutting @ " << cut_val << endl;

    vector<VertexID> new_verts;

    vector<HalfEdgeID> hidvec;

    VertexAttributeVector<int> touched(m.allocated_vertices(), 0);

    for(HalfEdgeID hid : m.halfedges())

    {

        Walker w = m.walker(hid);

        double bval = fun[w.vertex()];

        double aval = fun[w.opp().vertex()];

        if(aval<bval && aval < cut_val && cut_val <= bval)

        {

            double t = (cut_val-aval)/(bval-aval);

            Vec3d p = t*m.pos(w.vertex()) + (1.0-t)*m.pos(w.opp().vertex());

            VertexID vnew = m.split_edge(hid);

            m.pos(vnew) = p;

            status[vnew] = COOKED;

            fun[vnew] = cut_val;

            touched[vnew] = 1;

            new_verts.push_back(vnew);

        }

    }

    connect_touched(m, touched);

void remove_v5_vertices(Manifold& m, VertexAttributeVector<int>& ls_id, int id_ceiling)

{

    for(auto vid: m.vertices())

    {

        int id = ls_id[vid];

        if(id>0 && id<id_ceiling)

        {

            vector<HalfEdgeID> hit_list;

            for(Walker w = m.walker(vid); !w.full_circle(); w=w.circulate_vertex_ccw())

                if(ls_id[w.vertex()] >= id_ceiling)

                    hit_list.push_back(w.halfedge());

            if(hit_list.size()>1)

            {

                sort(begin(hit_list),end(hit_list),

                     [&m](HalfEdgeID a, HalfEdgeID b){return length(m,a)>length(m, b);});

                for(int i=1;i<hit_list.size();++i)

                {

                    HalfEdgeID hid = hit_list[i];

                    Walker w = m.walker(hid);

                    HalfEdgeID surplus_h = w.opp().prev().halfedge();

                    if(precond_collapse_edge(m,surplus_h))

                    {

                        m.collapse_edge(surplus_h,true);

                        m.merge_faces(m.walker(hid).face(), hid);

                    }

                }

            }

        }

    }

}

void remove_wedges(Manifold& m, VertexAttributeVector<int>& ls_id)

{

    for(auto hid : m.halfedges())

        if(m.in_use(hid))

        {

            Walker w = m.walker(hid);

            if(ls_id[w.vertex()]>0 &&

               ls_id[w.vertex()] == ls_id[w.opp().vertex()] &&

               ls_id[w.next().vertex()]>0 &&

               no_edges(m, w.face())==3 &&

               precond_collapse_edge(m, hid))

                m.collapse_edge(hid, true);

        }

}

    VertexAttributeVector<int> level_id(m.allocated_vertices(), 0);

    vector<VertexID> cut_verts = cut_mesh(m, fun, cut_val, status);

    label_connected_components(m, cur_id, cut_verts, fun, status, level_id);

    remove_vertices(m, fun, status, level_id, vmin, cut_val);

    flip_edges(m, level_id);

    for(int id=1;id<cur_id;++id) {

        vector<VertexID> loop_verts = separate_edges(m, cut_verts, id, fun, status, level_id);

        smooth_cooked_vertices(m, loop_verts, status, 50);

    }

        int old_id = cur_id;

        vector<VertexID> loop_verts;

        for(int id=old_id;id<cur_id;++id)

        {

            loop_verts.insert(loop_verts.end(), tmp.begin(),tmp.end());

        }

        remove_v5_vertices(m, level_id, old_id);

    dual(m);

    m.cleanup();

}

void simplify_polar_mesh(Manifold& m, double frac)

{

    int N_orig = m.no_vertices();

    HalfEdgeAttributeVector<int> touched(m.allocated_halfedges(),0);

    vector<pair<double,vector<HalfEdgeID>>> work_vector;

    for(auto hid : m.halfedges())

    {

        Walker w = m.walker(hid);

        if(!touched[hid] &&

           no_edges(m, w.face()) == 4 &&

           no_edges(m, w.opp().face()) == 4)

        {

            int n=0;

            vector<HalfEdgeID> hvec;

            bool abort=false;

            for(int i=0;i<2;++i)

            {

                Walker w0 = (i==0?w:w.opp().next().next().opp());

                while(!touched[w0.halfedge()])

                {

                    touched[w0.halfedge()] = 1;

                    touched[w0.opp().halfedge()] = 1;

                    hvec.push_back(w0.halfedge());

                    ++n;

                    l_sum += sqr(area(m,w0.face()));

                    if(valency(m,w0.vertex()) != 4 || valency(m, w0.opp().vertex()) != 4)

                        abort = true;

                    if(no_edges(m, w0.face())==4)

                        w0 = w0.next().next().opp();

                    else {

                        if (valency(m,w0.next().vertex())<5)

                            abort = true;

                        break;

                    }

                }

            }

            if(abort)

                continue;

            else

                work_vector.push_back(make_pair(l_sum/sqr(n),hvec));

        }

    }

    sort(begin(work_vector),end(work_vector));

    for(auto work : work_vector)

    {

        for(HalfEdgeID h :work.second)

            if(m.in_use(h) && precond_collapse_edge(m, h))

                m.collapse_edge(h,true);

            else cout << "Collapse failed ..." << endl;

        if(N_orig * frac > m.no_vertices())

            break;

    }

    m.cleanup();

}

using namespace Geometry;

const Mat4x4d fit_bounding_volume(const Vec3d& p0,

                                  const Vec3d& p7,

                                  float buf_reg,

                                  int vol_dim)

{

    Vec3d sz = p7 - p0;

    Vec3i dims(vol_dim);

    Vec3d scal_vec = (Vec3d(dims)-Vec3d(2*buf_reg+2))/sz;

    float scal = min(scal_vec[0], min(scal_vec[1],scal_vec[2]));

    Mat4x4d m = translation_Mat4x4d(Vec3d(0)+Vec3d(buf_reg+1));

    m *= scaling_Mat4x4d(Vec3d(scal));

    m *= translation_Mat4x4d(-p0);

    return m;

}

void smooth_and_refit(HMesh::Manifold& m_in,  HMesh::Manifold& m_ref, int iter, int vol_dim)

{

    AABBTree tree;

    Manifold m_tmp = m_ref;

    Vec3d p0,p7;

    bbox(m_tmp, p0, p7);

    Mat4x4d mat = fit_bounding_volume(p0, p7, 3, vol_dim);

    Mat4x4d imat = invert(mat);

        m_tmp.pos(vid) = mat.mul_3D_point(m_tmp.pos(vid));

    build_AABBTree(m_tmp, tree);

    cout << "Done building AABB tree" << endl;

    RGrid<float> grid(Vec3i(vol_dim),FLT_MAX);

    auto F = [&](const Vec3i& p,float& v)

    {

        v =  tree.compute_signed_distance(Vec3f(p));

    };

    for_each_voxel(grid, F);

    cout << "Made distance field" << endl;

    TrilinFilter<RGrid<float>> df(&grid);

    for (int i=0;i<iter; ++i)

        auto vsafe = m_in.positions_attribute_vector();

        laplacian_smooth(m_in,.25,10);

        for(auto vid : m_in.vertices())

        {

            double w = exp(4-valency(m_in, vid));

            m_in.pos(vid) += (1-w) * vsafe[vid];

        }

        for(auto hid : m_in.halfedges())

        auto npos = m_in.positions_attribute_vector();

        for(auto vid : m_in.vertices())

        {

            npos[vid] += 0.5 * avg * normal(m_in, vid);

        }

        m_in.positions_attribute_vector() = npos;

        for(auto vid : m_in.vertices())

        {

            Vec3f pos = Vec3f(mat.mul_3D_point(m_in.pos(vid)));

                m_in.pos(vid) = imat.mul_3D_point(Vec3d(pos - df(pos) * df.grad(pos)));

            else

                cout << "outside domain" << endl;