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 <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 <complex>

using namespace CGLA;

using namespace std;

using namespace HMesh;

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

enum VertexStatus {

    RAW, POLE, SADDLE, COOKED

};

void kill_wedges(Manifold& m, const VertexAttributeVector<double>& fun, const VertexAttributeVector<VertexStatus>& status)

{

    for(auto hid : m.halfedges())

        if(m.in_use(hid)){

            Walker w = m.walker(hid);

            if(same_level(fun[w.vertex()], fun[w.opp().vertex()]) &&

               status[w.vertex()] == COOKED &&

               status[w.opp().vertex()] == COOKED &&

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

               status[w.next().vertex()] == COOKED &&

               precond_collapse_edge(m, hid))

                m.collapse_edge(hid);

    }

}

void high_gradient_triangulate(Manifold& m, FaceID _f, const VertexAttributeVector<double>& fun)

{

    queue<FaceID> fq;

    fq.push(_f);

    while(!fq.empty()) {

        FaceID f = fq.front();

        fq.pop();

        if(no_edges(m, f)<=3)

            continue;

        // Create a vector of vertices.

        vector<VertexID> verts;

        circulate_face_ccw(m, f, [&](VertexID v){verts.push_back(v);});

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

                    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 max_grad=-FLT_MAX;

        size_t max_k = -1;

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

            int i = vpairs[k].first;

            int j = vpairs[k].second;

            double grad = abs(fun[verts[i]]-fun[verts[j]])/length(m.pos(verts[i]) - m.pos(verts[j]));

            if(grad>max_grad){

                max_grad = grad;

                max_k = k;

            }

        }

        assert(max_k != -1);

        {

            int i = vpairs[max_k].first;

            int j = vpairs[max_k].second;

            if(m.in_use(verts[i]) && m.in_use(verts[j])) {

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

                fq.push(f);

                fq.push(f2);

            }

        }

    }

}

FaceID merge_face_in_1_ring(Manifold& m, VertexID v)

{

    Walker w=m.walker(v);

    HalfEdgeID h = w.next().halfedge();

    if(m.remove_vertex(v))

    {

        return m.close_hole(h);

    }

    return InvalidFaceID;

}

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;

    }

}

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

{

    cout << "Cutting at " << cut_val << endl;

    vector<VertexID> new_verts;

    vector<HalfEdgeID> hidvec;

    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;

            new_verts.push_back(vnew);

        }

    }

    for(FaceID fid : m.faces())

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

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

            {

                Walker w0 = w;

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

                do

                {

                    if(status[w.vertex()] == COOKED &&

                       w.next().halfedge() != w0.halfedge() &&

                       same_level(fun[w.vertex()],cut_val))

                    {

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

                        break;

                    }

                    w = w.next();

                }

                while(!w.full_circle());

                break;

            }

    return new_verts;

}

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

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

{

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

    {

        VertexAttributeVector<Vec3d> npos(m.allocated_vertices(), Vec3d(0));

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

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

    }

}

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 classify_vertices(Manifold& m, VertexAttributeVector<double>& fun, VertexAttributeVector<VertexStatus>& status)

{

    for(auto vid : m.vertices())

    {

        double val = fun[vid];

        auto w = m.walker(vid);

        int above_below = ssb(fun[w.vertex()],val);

        int changes = 0;

        w = w.circulate_vertex_ccw();

        for(; !w.full_circle();w = w.circulate_vertex_ccw())

        {

            int above_below_new = ssb(fun[w.vertex()],val);

            if(above_below != above_below_new)

                ++changes;

            above_below=above_below_new;

        }

        if(changes==0)

            status[vid]=POLE;

        else if(changes>=4)

            status[vid]=SADDLE;

        else if(changes==2)

            status[vid]=RAW;

        else assert(0);

    }

}

void remove_vertices(Manifold& m, VertexAttributeVector<double>& fun, VertexAttributeVector<VertexStatus>& status, 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 = merge_face_in_1_ring(m,vid);

        if(fid != InvalidFaceID)

            high_gradient_triangulate(m, fid, fun);

    }

}

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

    classify_vertices(m, fun, status);

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

    double base_cut_val = vmin + 0.02*(vmax-vmin);

    vector<VertexID> new_verts = cut_mesh(m, fun, base_cut_val, status);

    smooth_cooked_vertices(m, new_verts, fun, status, 1);

    cout << "Base " << base_cut_val << endl;

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

    {

        int _iter=iter;//%2==0?iter/2:-1-iter/2;

        double cut_val = base_cut_val + interval * _iter;

        cout << "i,c " << _iter << "," << cut_val << endl;

        remove_vertices(m, fun, status, cut_val, cut_val+interval);

        vector<VertexID> new_verts = cut_mesh(m, fun,_iter<0?cut_val:cut_val+interval, status);

        kill_wedges(m, fun, status);

        smooth_cooked_vertices(m, new_verts, fun, status, 1);

    }

    int C=0, R=0, P=0,S=0;

    for(auto vid : m.vertices())

    {

        switch(status[vid]) {

            case RAW: ++R; break;

            case COOKED: ++C; break;

            case POLE: ++P; break;

            case SADDLE: ++S; break;

        }

    }

    cout << "R=" << R << " C=" << C << " P=" << P << " S=" << S << 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;

}