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 "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;}

struct LevelSetInfo

{

    int id;

    int no_vertices;

    Vec3d avg_pos;

    double rad;

    double length;

    int components;

    HalfEdgeID h;

    double fun_value;

    void print()

    {

        cout

        << " id : " << id

        << " no vertices : " << no_vertices

        << " avg pos : " << avg_pos

        << " rad : " << rad

        << " length : " << length

        << " # : " << components

        << " F : " << fun_value << endl;

    }

};

void compute_edge_weights(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(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 segment_manifold(Manifold& m, HalfEdgeAttributeVector<int>& ls_id,

                      FaceAttributeVector<int>& segmentation,

                      vector<vector<int>>& boundaries)

{

    segmentation = FaceAttributeVector<int>(m.no_faces(), -1);

    int SEG_NO = 0;

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

    {

        if (segmentation[*f] == -1)

        {

            queue<FaceID> q;

            q.push(*f);

            vector<int> bound(0);

            while (!q.empty())

            {

                FaceID face = q.front();

                q.pop();

                segmentation[face] = SEG_NO;

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

                {

                    if(ls_id[w.halfedge()] == -1)

                    {

                        FaceID fopp = w.opp().face();

                        if(segmentation[fopp] == -1) {

                            q.push(fopp);

                        }

                    }

                    else {

                        bound.push_back(ls_id[w.halfedge()]);

                    }

                }

            }

            sort(begin(bound), end(bound));

            auto new_end = unique(begin(bound), end(bound));

            boundaries.push_back(vector<int>(begin(bound),new_end));

            cout << "Boundaries of seg# " << SEG_NO << " : ";

            for(auto be:boundaries[SEG_NO])

                cout << be << " ";

            cout << endl;

            SEG_NO += 1;

        }

    }

}

void shortest_edge_triangulate_face(Manifold& m, FaceID f0, 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]]) &&

                   valency(m, verts[i])<5 && valency(m,verts[j])<5

                   )

                    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 = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));

            if(len<min_len){

                min_len = len;

                min_k = k;

            }

        }

        assert(min_k != -1);

        {

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

            int i = vpairs[min_k].first;

            int j = vpairs[min_k].second;

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

        }

    }

}

void shortest_edge_split_face(Manifold& m, FaceID f0, 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())

        {

            verts.push_back(w.vertex());

        }

        // Find vertex pairs that may be connected.

        vector<pair<int,int> > vpairs;

        const int N = verts.size();

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

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

                int jj = (j+i)%N;

                if(!connected(m, verts[i], verts[jj]) &&

                   (level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[jj]]) &&

                   (level_set_id_vertex[verts[(i+1)%N]] == level_set_id_vertex[verts[i]]) &&

                   (level_set_id_vertex[verts[(i+N-1)%N]] == level_set_id_vertex[verts[i]]) &&

                   (level_set_id_vertex[verts[(jj+1)%N]] == level_set_id_vertex[verts[jj]]) &&

                   (level_set_id_vertex[verts[(jj+N-1)%N]] == level_set_id_vertex[verts[jj]]))

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

            }

        }

        if(vpairs.empty()) 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 = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));

            if(len<min_len){

                min_len = len;

                min_k = k;

            }

        }

        assert(min_k != -1);

        {

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

            int i = vpairs[min_k].first;

            int j = vpairs[min_k].second;

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

            if(no_edges(m, f)>5)

                face_queue.push(f);

            if(no_edges(m, f_new)>5)

                face_queue.push(f_new);

        }

    }

}

struct EdgeQElem {

    float priority;

    HalfEdgeID he;

    EdgeQElem(float _priority, HalfEdgeID _he): priority(_priority), he(_he) {}

};

bool operator<(const EdgeQElem& l, const EdgeQElem& r)

{

    return l.priority < r.priority;

}

class FunctionalDifference: public EnergyFun

{

    VertexAttributeVector<float>& fun;

    VertexAttributeVector<int>& status;

public:

    FunctionalDifference(VertexAttributeVector<float>& _fun, VertexAttributeVector<int>& _status): fun(_fun), status(_status) {}

    virtual double delta_energy(const Manifold& m, HalfEdgeID h) const

    {

        Walker w = m.walker(h);

        if(status[w.vertex()] == 1 && status[w.opp().vertex()]==1)

            return DBL_MAX;

        return sqr_length(m.pos(w.next().vertex())-m.pos(w.opp().next().vertex()))/(1e-6+abs(fun[w.next().vertex()]-fun[w.opp().next().vertex()])) - sqr_length(m.pos(w.vertex())-m.pos(w.opp().vertex()))/(1e-6+abs(fun[w.vertex()]-fun[w.opp().vertex()]));

    }

};

class TriangleQualityValence: public EnergyFun

{

    VertexAttributeVector<int>& idv;

    MinAngleEnergy mae;

    ValencyEnergy vae;

public:

    TriangleQualityValence(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}

    virtual double delta_energy(const Manifold& m, HalfEdgeID h) const

    {

        Walker w = m.walker(h);

        if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] ||

           idv[w.vertex()] == idv[w.opp().vertex()])

            return DBL_MAX;

        VertexID v1 = w.opp().vertex();

        VertexID v2 = w.vertex();

        VertexID vo1 = w.next().vertex();

        VertexID vo2 = w.opp().next().vertex();

        int val1  = valency(m, v1);

        int val2  = valency(m, v2);

        int valo1 = valency(m, vo1);

        int valo2 = valency(m, vo2);

        // The optimal valency is four for a boundary vertex

        // and six elsewhere.

        int t1 = boundary(m, v1) ? 4 : 6;

        int t2 = boundary(m, v2) ? 4 : 6;

        int to1 = boundary(m, vo1) ? 4 : 6;

        int to2 = boundary(m, vo2) ? 4 : 6;

        int before =

        max(max(sqr(val1-t1),sqr(val2-t2)), max(sqr(valo1-to1), sqr(valo2-to2)));

        int after =

        max(max(sqr(valo1+1-to1),sqr(val1-1-t1)), max(sqr(val2-1-t2),sqr(valo2+1-to2)));

        return static_cast<double>(after-before);

        //        return vae.delta_energy(m,h);

        float la = length(m.pos(w.next().vertex())-m.pos(w.opp().next().vertex()));

        float lb = length(m.pos(w.vertex())-m.pos(w.opp().vertex()));

        return la-lb;

        return mae.delta_energy(m,h);

    }

};

class TriangleQuality: public EnergyFun

{

    VertexAttributeVector<int>& idv;

    MinAngleEnergy mae;

    ValencyEnergy vae;

public:

    TriangleQuality(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}

    virtual double delta_energy(const Manifold& m, HalfEdgeID h) const

    {

        Walker w = m.walker(h);

        if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] ||

           idv[w.vertex()] == idv[w.opp().vertex()])

            return DBL_MAX;

        return mae.delta_energy(m,h);

    }

};

Vec3d grad(HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun, HMesh::FaceID f)

{

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

        return Vec3d(0.0);

    Vec3d n = normal(m, f);

    Vec3d gsum(0.0);

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

    {

        VertexID avid = w.opp().vertex();

        VertexID bvid = w.vertex();

        VertexID cvid = w.next().vertex();

        Vec3d gdir = normalize(cross(n, m.pos(bvid) - m.pos(avid)));

        double l = dot(m.pos(cvid)-m.pos(avid), gdir);

        gdir *= fun[cvid]/l;

        gsum += gdir;

    }

    return gsum;

}

Vec2d extend_fun2(HMesh::Manifold& m, HMesh::HalfEdgeID h,

                  HMesh::VertexAttributeVector<double>& fun,

                  HMesh::VertexAttributeVector<Vec2d>& fun2)

{

    Walker w = m.walker(h);

    FaceID f = w.face();

    Vec3d a = m.pos(w.opp().vertex());

    Vec3d b = m.pos(w.vertex());

    Vec3d c = m.pos(w.next().vertex());

    Vec3d g = grad(m, fun, f);

    Vec3d n = normal(m, f);

    Vec3d N = normalize(cross(g, n));

    float dot_nba = dot(N,b-a);

    Vec2d cval;

    if(dot_nba > 1e-5)

    {

        float t = dot(N, c-a)/dot_nba;

        if(t>0 && t < 1) {

            Vec2d aval = fun2[w.opp().vertex()];

            Vec2d bval = fun2[w.vertex()];

            return (1-t)*aval + t*bval;

        }

    }

    return Vec2d(0);

}

double solve_for_orthogonal_gradients(HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,

                                      HalfEdgeID h, double beta, double gamma)

{

    Vec3d n = normal(m, m.walker(h).face());

    Vec3d uvw[3];

    Vec3d abc;

    Walker w = m.walker(h);

    VertexID vid[] = {w.next().vertex(), w.opp().vertex(),w.vertex()};

    for(int i=0; !w.full_circle(); w = w.next(),++  i)

    {

        uvw[i] = normalize(cross(n, m.pos(vid[(i+1)%3]) - m.pos(vid[(i+2)%3])));

        uvw[i] /= dot(m.pos(vid[i])-m.pos(vid[(i+1)%3]), uvw[i]);

        abc[i] = fun[vid[i]];

    }

    Vec3d A(dot(uvw[0],uvw[0]),dot(uvw[0],uvw[1]),dot(uvw[0],uvw[2]));

    Vec3d B(dot(uvw[1],uvw[0]),dot(uvw[1],uvw[1]),dot(uvw[1],uvw[2]));

    Vec3d C(dot(uvw[2],uvw[0]),dot(uvw[2],uvw[1]),dot(uvw[2],uvw[2]));

    return - (beta * dot(B, abc) + gamma * dot(C, abc))/dot(A, abc);

}

void recompute_circle_param(Manifold& m, VertexID v, VertexAttributeVector<double>& fun, VertexAttributeVector<Vec2d>& circle_pos, double f_origin)

{

    Vec2d c_new(0);

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

    {

        Walker w_base = w.next();

        Vec2d X = circle_pos[w_base.opp().vertex()];

        Vec2d Y(-X[1], X[0]);

        Vec2d vec = circle_pos[w_base.vertex()];

        double angle = atan2(max(-1.0,min(1.0,dot(vec,Y))),

                             max(-1.0,min(1.0,dot(vec,X))));

        double new_angle = solve_for_orthogonal_gradients(m, fun, w_base.halfedge(), 0, angle);

        new_angle = max(-M_PI/4, min(M_PI/4,new_angle));

        c_new += cos(new_angle)*X + sin(new_angle)*Y;

    }

    circle_pos[v] = normalize(c_new);

}

template<typename T>

void orthogonal_trajectories(Manifold& m, VertexID v, VertexAttributeVector<double>& f, double min_length_fun,

                             VertexAttributeVector<T>& h_in, VertexAttributeVector<T>& h_out)

{

    h_out[v] = T(0);

    T num(0);

    double denom = 0.0001;

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

    {

        double weight = area(m, w.face());

        Vec3d a = m.pos(v);

        Vec3d b = m.pos(w.vertex());

        Vec3d c = m.pos(w.next().vertex());

        Vec3d fvec(f[v],f[w.vertex()], f[w.next().vertex()]);

        if(sqr_length(fvec)>0)

            fvec.normalize();

        double f_a = fvec[0];

        double f_b = fvec[1];

        double f_c = fvec[2];

        T h_b = h_in[w.vertex()];

        T h_c = h_in[w.next().vertex()];

        Vec3d e_a = c-b;

        Vec3d e_b = a-c;

        Vec3d e_c = b-a;

        Vec3d N = normalize(cross(e_b,e_c)+cross(e_c,e_a)+cross(e_a,e_b));

        Vec3d g_a = cross(N, normalize(e_a));

        g_a /= dot(g_a, e_b);

        Vec3d g_b = cross(N, normalize(e_b));

        g_b /= dot(g_b, e_c);

        Vec3d g_c = cross(N, normalize(e_c));

        g_c /= dot(g_c, e_a);

        num -= weight * (h_b * (f_a*dot(g_a, g_b) + f_b*dot(g_b,g_b) + f_c*dot(g_b,g_c)) +

                         h_c * (f_a*dot(g_a, g_c) + f_b*dot(g_b,g_c) + f_c*dot(g_c, g_c)));

        denom += weight * (f_a * dot(g_a,g_a)+ f_b*dot(g_a,g_b) +f_c*dot(g_a, g_c));

    }

    h_out[v] = num / denom;

}

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

{

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

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

    // ----------------------------------------

    cout << "Tracing level set curves" << endl;

    for(int cut_no=1; cut_no<divisions+1;++cut_no)

    {

        double cut_val = vmin + cut_no * interval;

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

                fun[vnew] = cut_val;

            }

        }

        for(FaceID fid : m.faces())

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

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

                    {

                        Walker w0 = w;

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

                        do

                        {

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

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

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

                            {

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

                                break;

                            }

                            w = w.next();

                        }

                        while(!w.full_circle());

                        break;

                    }

    }

    shortest_edge_triangulate(m);

    // ----------------------------

    cout << "Numbering the level sets" << endl;

    vector<LevelSetInfo> level_set_info;

    HalfEdgeAttributeVector<int> level_set_id(m.allocated_halfedges(), -1);

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

    int no_id=0;

    for(HalfEdgeID hid : m.halfedges())

    {

        Walker w = m.walker(hid);

        if(status[w.vertex()] == 1 && status[w.opp().vertex()] == 1 &&

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

           level_set_id[w.halfedge()] == -1)

        {

            LevelSetInfo lsi;

            lsi.id = no_id;

            lsi.fun_value = fun[w.vertex()];

            lsi.components = 1;

            lsi.no_vertices = 0;

            double level_set_length = 0;

            while(level_set_id[w.halfedge()] != no_id)

            {

                level_set_length += length(m,w.halfedge());

                level_set_id[w.halfedge()] = no_id;

                level_set_id[w.opp().halfedge()] = no_id;

                level_set_id_vertex[w.vertex()] = no_id;

                w = w.next();

                while(!(status[w.vertex()] == 1 && same_level(fun[w.vertex()], fun[w.opp().vertex()])))

                    w = w.circulate_vertex_cw();

                lsi.no_vertices += 1;

            }

            lsi.h = w.halfedge();

            lsi.length = level_set_length;

            level_set_info.push_back(lsi);

            level_set_info.back().print();

            ++no_id;

            assert(level_set_info.size()==no_id);

        }

    }

    cout << "ids" << no_id << " " << level_set_info.size() << endl;

    for(int i=0;i<level_set_info.size(); ++i)

        for(int j=0;j<level_set_info.size(); ++j)

            if(i != j && same_level(level_set_info[i].fun_value, level_set_info[j].fun_value))

                level_set_info[i].components += 1;

    HalfEdgeID min_length_h;

    double min_length_fun;

    double min_length=FLT_MAX;

    int min_length_id = -1;

    for(int i=0;i<level_set_info.size(); ++i)

    {

        LevelSetInfo& lsi = level_set_info[i];

        if(lsi.components==1 && lsi.length<min_length)

        {

            min_length_h = lsi.h;

            min_length_fun = lsi.fun_value;

            min_length = lsi.length;

            min_length_id = lsi.id;

        }

    }

    // ----------------------------

    cout << "Remove vertices not on level set curves" << endl;

}

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

                     double& vmin, double& vmax)

{

    VertexIDIterator vid = m.vertices_begin();

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

    }

}