Subversion Repositories gelsvn

//

//  polarize.cpp

//  GEL

//

//  Created by J. Andreas Bærentzen on 18/03/12.

//  Copyright 2012 __MyCompanyName__. All rights reserved.

//

#include <queue>

#include <iomanip>

#include "harmonics.h"

#include "polarize.h"

#include <CGLA/Vec2d.h>

#include <CGLA/Mat4x4d.h>

#include <CGLA/Vec3i.h>

#include <LinAlg/LapackFunc.h>

#include <HMesh/triangulate.h>

#include <HMesh/obj_save.h>

#include <HMesh/curvature.h>

#include <HMesh/quadric_simplify.h>

#include <HMesh/mesh_optimization.h>

#include <HMesh/smooth.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>

#include <complex>

using namespace CGLA;

using namespace Geometry;

using namespace std;

using namespace HMesh;

using namespace GLGraphics;

inline bool same_level(double a, double b) {return abs(a-b) < 1e-10;}

enum VertexStatus {

    RAW, COOKED

};

void shortest_edge_triangulate_face(Manifold& m, FaceID f0, const VertexAttributeVector<int>& level_set_id_vertex)

{

    queue<FaceID> face_queue;

    face_queue.push(f0);

    while(!face_queue.empty())

    {

        FaceID f = face_queue.front();

        face_queue.pop();

        // Create a vector of vertices.

        vector<VertexID> verts;

        for(Walker w = m.walker(f); !w.full_circle(); w = w.circulate_face_ccw())

        {

            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;

        // Find vertex pairs that may be connected.

        vector<pair<int,int> > vpairs;

        const int N = verts.size();

        for(int i = 0; i < N - 2; ++i){

            for(int j = i + 2; j < N; ++j){

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

                   !connected(m, verts[i], verts[j]) &&

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

                   )

                    vpairs.push_back(pair<int,int>(i, j));

            }

        }

        if(vpairs.empty()){

            cout << "Warning: could not triangulate a face."

            << "Probably a vertex appears more than one time in other vertex's one-ring" << endl;

            continue;

        }

        /* For all vertex pairs, find the edge lengths. Combine the

         vertices forming the shortest edge. */

        float min_len=FLT_MAX;

        int min_k = -1;

        for(size_t k = 0; k < vpairs.size(); ++k){

            int i = vpairs[k].first;

            int j = vpairs[k].second;

            float len = length(m.pos(verts[i]) - m.pos(verts[j]));

            //            cout << m.pos(verts[i]) << m.pos(verts[j]) << len << endl;

            if(len<min_len){

                min_len = len;

                min_k = k;

            }

        }

        assert(min_k != -1);

        if(min_k>=0)

        {

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

            int i = vpairs[min_k].first;

            int j = vpairs[min_k].second;

            //           cout << i << " " << j << " " << min_k << " " << vpairs.size() << endl;

            FaceID f_new = m.split_face_by_edge(f, verts[i], verts[j]);

            if(no_edges(m, f)>3)

                face_queue.push(f);

            if(no_edges(m, f_new)>3)

                face_queue.push(f_new);

        }

        else

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

    }

}

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

{

    vector<HalfEdgeID> hvec;

    for(HalfEdgeID hid : m.halfedges())

    {

        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()])

            hvec.push_back(hid);

    }

    ValencyEnergy mae;

    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 compute_edge_weights(const Manifold& m, HalfEdgeAttributeVector<double>& edge_weights, FaceAttributeVector<int>& included)

{

    edge_weights = HalfEdgeAttributeVector<double>(m.allocated_halfedges(), 0);

    for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f)

        if(included[*f])

        {

            for(Walker wv = m.walker(*f); !wv.full_circle(); wv = wv.circulate_face_ccw())

            {

                HalfEdgeID h = wv.halfedge();

                Vec3d p1(m.pos(wv.vertex()));

                Vec3d p2(m.pos(wv.next().vertex()));

                Vec3d p0(m.pos(wv.opp().vertex()));

                double ang = acos(min(1.0, max(-1.0, dot(normalize(p1-p0), normalize(p2-p0)))));

                double ang_opp = acos(min(1.0, max(-1.0, dot(normalize(p2-p1), normalize(p0-p1)))));

                double l = (p1-p0).length();

                edge_weights[h] += tan(ang/2) / l;

                edge_weights[wv.opp().halfedge()] += tan(ang_opp/2) / l;

            }

        }

}

template<typename T>

void smooth_fun(const Manifold& m,

                VertexAttributeVector<int>& nailed,

                VertexAttributeVector<T>& fun, int iter=100)

{

    HalfEdgeAttributeVector<double> edge_weights;

    FaceAttributeVector<int> included(m.allocated_faces(),1);

    compute_edge_weights(m,edge_weights, included);

    VertexAttributeVector<T> new_fun(m.no_vertices());

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

    {

        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)

            if(!nailed[*v])

            {

                double w_sum = 0;

                new_fun[*v] = T(0);

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

                {

                    double w = edge_weights[wv.halfedge()];

                    new_fun[*v] += w * fun[wv.vertex()];

                    w_sum += w;

                }

                new_fun[*v] /= w_sum;

            }

            else

                new_fun[*v] = fun[*v];

        fun = new_fun;

    }

}

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

                                const VertexAttributeVector<double>& fun,

                                const VertexAttributeVector<VertexStatus>& status,

                                VertexAttributeVector<int>& level_id)

{

    for(auto vid : vertices)

    {

        if(level_id[vid]>0)

            continue;

        queue<VertexID> vq;

        vq.push(vid);

        while(!vq.empty())

        {

            VertexID v = vq.front();

            vq.pop();

            level_id[v] = cur_id;

            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()]))

                    if(level_id[w.vertex()]==0)

                        vq.push(w.vertex());

                    else if(level_id[w.vertex()] != cur_id)

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

            }

        }

        ++cur_id;

    }

}

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

{

    for(FaceID fid : m.faces())

        for(Walker w = m.walker(fid);!w.full_circle(); w = w.next())

            if(touched[w.vertex()])

            {

                Walker w0 = w;

                w = w.next().next();

                do

                {

                    if(touched[w.vertex()] && w.next().halfedge() != w0.halfedge())

                    {

                        m.split_face_by_edge(fid, w0.vertex(), w.vertex());

                        break;

                    }

                    w = w.next();

                }

                while(!w.full_circle());

                break;

            }

}

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);

    return new_verts;

}

void smooth_cooked_vertices(Manifold& m, const vector<VertexID>& vertices,

                            VertexAttributeVector<VertexStatus>& status, int iter=7)

{

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

    {

        VertexAttributeVector<Vec3d> npos = m.positions_attribute_vector();

        VertexAttributeVector<int> cnt(m.allocated_vertices(), 1);

        for(auto vid : vertices)

            if(m.in_use(vid))

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

                    if(status[vid]==COOKED) {

                        npos[vid] += m.pos(w.vertex());

                        ++cnt[vid];

                    }

        for(auto vid : vertices)

            if(m.in_use(vid))

                m.pos(vid) = (npos[vid])/cnt[vid];

    }

}

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);

        }

}

int ssb(double a, double b)

{

    const double THRESH = 1e-6;

    double delta = a-b;

    if(delta<-THRESH) return -1;

    if(delta>THRESH) return 1;

    return 0;

}

void remove_vertices(Manifold& m, VertexAttributeVector<double>& fun,

                     const VertexAttributeVector<VertexStatus>& status,

                     const VertexAttributeVector<int>& level_id, double rmin, double rmax)

{

    vector<VertexID> vertices_to_push;

    for(auto vid: m.vertices())

        if(status[vid]==RAW &&

           fun[vid] >= rmin &&

           fun[vid] <= rmax)

            vertices_to_push.push_back(vid);

    random_shuffle(begin(vertices_to_push), end(vertices_to_push));

    for(auto vid : vertices_to_push)

    {

        FaceID fid = m.merge_one_ring(vid);

        if(fid != InvalidFaceID)

            shortest_edge_triangulate_face(m, fid, level_id);

    }

}

void polarize_mesh(Manifold& m, VertexAttributeVector<double>& fun, double vmin, double vmax, const int divisions, VertexAttributeVector<Vec2d>& parametrization)

{

    VertexAttributeVector<VertexStatus> status(m.allocated_vertices(), RAW);

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

    int cur_id = 1;

    double interval = (vmax-vmin)/(divisions+1);

    double cut_val = vmin + 0.065*(vmax-vmin);

    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);

    }

    for(int iter=0;iter<divisions+2;++iter)

    {

        cut_val += interval;

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

        int old_id = cur_id;

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

        remove_vertices(m, fun, status, level_id, cut_val-1.01*interval, cut_val);

        vector<VertexID> loop_verts;

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

        {

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

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

        }

        remove_wedges(m, level_id);

        remove_v5_vertices(m, level_id, old_id);

        smooth_cooked_vertices(m, loop_verts, status);

    }

    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)

        {

            double l_sum=0;

            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);

    for(auto vid : m_tmp.vertices())

        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) *= w;

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

        }

        float avg = 0;

        for(auto hid : m_in.halfedges())

            avg += length(m_in, hid);

        avg /= m_in.no_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)));

            if(df.in_domain(pos))

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

            else

                cout << "outside domain" << endl;

        }

    }

}

void make_height_fun(const HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,

                     double& vmin, double& vmax)

{

    VertexIDIterator vid = m.vertices_begin();

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

    smooth_fun(m, nailed, fun,0);

    vmin = vmax = m.pos(*vid)[1];

    for(; vid != m.vertices_end(); ++vid)

    {

        double v = dot(m.pos(*vid),Vec3d(0.0,1,0.00));

        fun[*vid] = v;

        vmin = min(v, vmin);

        vmax = max(v, vmax);

    }

}

void make_adf_fun(HMesh::Manifold& m, double t, HMesh::VertexAttributeVector<double>& F,

                  double& vmin, double& vmax)

{

    Harmonics harm(m);

    harm.compute_adf(F, t);

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

    vmin = vmax = F[*m.vertices_begin()];

    for(auto vid : m.vertices())

    {

        vmin = min(F[vid], vmin);

        vmax = max(F[vid], vmax);

    }

    cout << vmin << vmax << endl;

}