Subversion Repositories gelsvn

Rev

Rev 640 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 640 Rev 642
Line 6... Line 6...
6 
//  Copyright 2012 __MyCompanyName__. All rights reserved.
 6 
//  Copyright 2012 __MyCompanyName__. All rights reserved.
7 
//
 7 
//
8 
#include <queue>
 8 
#include <queue>
9 
#include <iomanip>
 9 
#include <iomanip>
10 
 
 10 
 
- 
 
 11 
 
- 
 
 12 
#include "harmonics.h"
11 
#include "polarize.h"
 13 
#include "polarize.h"
12 
#include <CGLA/Vec2d.h>
 14 
#include <CGLA/Vec2d.h>
13 
#include <LinAlg/LapackFunc.h>
 15 
#include <LinAlg/LapackFunc.h>
14 
#include <HMesh/triangulate.h>
 16 
#include <HMesh/triangulate.h>
15 
#include <HMesh/obj_save.h>
 17 
#include <HMesh/obj_save.h>
Line 22... Line 24...
22 
 
 24 
 
23 
using namespace CGLA;
 25 
using namespace CGLA;
24 
using namespace std;
 26 
using namespace std;
25 
using namespace HMesh;
 27 
using namespace HMesh;
26 
 
 28 
 
27 
 
 - 
 
28 
 
 - 
 
29 
inline bool same_level(double a, double b) {return abs(a-b) < 1e-10;}
 29 
inline bool same_level(double a, double b) {return abs(a-b) < 1e-10;}
30 
 
 30 
 
- 
 
 31 
enum VertexStatus {
- 
 
 32 
    RAW, POLE, SADDLE, COOKED
- 
 
 33 
};
31 
 
 34 
 
- 
 
 35 
void kill_wedges(Manifold& m, const VertexAttributeVector<double>& fun, const VertexAttributeVector<VertexStatus>& status)
- 
 
 36 
{
- 
 
 37 
    for(auto hid : m.halfedges())
32 
struct LevelSetInfo
 38 
        if(m.in_use(hid)){
- 
 
 39 
            Walker w = m.walker(hid);
- 
 
 40 
            if(same_level(fun[w.vertex()], fun[w.opp().vertex()]) &&
- 
 
 41 
               status[w.vertex()] == COOKED &&
- 
 
 42 
               status[w.opp().vertex()] == COOKED &&
- 
 
 43 
               no_edges(m, w.face())==3 &&
- 
 
 44 
               status[w.next().vertex()] == COOKED &&
- 
 
 45 
               precond_collapse_edge(m, hid))
- 
 
 46 
                m.collapse_edge(hid);
- 
 
 47
- 
 
 48 
    }
- 
 
 49 
}
- 
 
 50 
 
- 
 
 51 
void high_gradient_triangulate(Manifold& m, FaceID _f, const VertexAttributeVector<double>& fun)
33 
{
 52 
{
34 
    int id;
 - 
 
35 
    int no_vertices;
 - 
 
36 
    Vec3d avg_pos;
 - 
 
37 
    double rad;
 - 
 
38 
    double length;
 - 
 
39 
    int components;
 - 
 
40 
    HalfEdgeID h;
 - 
 
41 
    double fun_value;
 - 
 
42
53
- 
 
 54 
    queue<FaceID> fq;
- 
 
 55 
    fq.push(_f);
- 
 
 56 
    while(!fq.empty()) {
- 
 
 57 
        FaceID f = fq.front();
43 
    void print()
 58 
        fq.pop();
- 
 
 59 
        if(no_edges(m, f)<=3)
- 
 
 60 
            continue;
- 
 
 61 
        // Create a vector of vertices.
- 
 
 62 
        vector<VertexID> verts;
- 
 
 63 
        circulate_face_ccw(m, f, [&](VertexID v){verts.push_back(v);});
- 
 
 64
- 
 
 65 
        // Find vertex pairs that may be connected.
- 
 
 66 
        vector<pair<int,int> > vpairs;
- 
 
 67 
        const int N = verts.size();
- 
 
 68 
        for(int i = 0; i < N - 2; ++i){
- 
 
 69 
            for(int j = i + 2; j < N; ++j){
- 
 
 70 
                if(verts[i] != verts[j] && !connected(m, verts[i], verts[j]))
- 
 
 71 
                    vpairs.push_back(pair<int,int>(i, j));
- 
 
 72 
            }
- 
 
 73 
        }
- 
 
 74 
        if(vpairs.empty()){
- 
 
 75 
            cout << "Warning: could not triangulate a face."
- 
 
 76 
            << "Probably a vertex appears more than one time in other vertex's one-ring" << endl;
- 
 
 77 
            continue;
- 
 
 78 
        }
- 
 
 79
- 
 
 80 
        /* For all vertex pairs, find the edge lengths. Combine the
- 
 
 81 
         vertices forming the shortest edge. */
- 
 
 82
- 
 
 83 
        float max_grad=-FLT_MAX;
- 
 
 84 
        size_t max_k = -1;
- 
 
 85 
        for(size_t k = 0; k < vpairs.size(); ++k){
- 
 
 86 
            int i = vpairs[k].first;
- 
 
 87 
            int j = vpairs[k].second;
- 
 
 88 
            double grad = abs(fun[verts[i]]-fun[verts[j]])/length(m.pos(verts[i]) - m.pos(verts[j]));
- 
 
 89
- 
 
 90 
            if(grad>max_grad){
- 
 
 91 
                max_grad = grad;
- 
 
 92 
                max_k = k;
- 
 
 93 
            }
- 
 
 94 
        }
- 
 
 95 
        assert(max_k != -1);
- 
 
 96
- 
 
 97 
        {
- 
 
 98 
            int i = vpairs[max_k].first;
- 
 
 99 
            int j = vpairs[max_k].second;
- 
 
 100 
            if(m.in_use(verts[i]) && m.in_use(verts[j])) {
- 
 
 101 
                FaceID f2 = m.split_face_by_edge(f, verts[i], verts[j]);
- 
 
 102 
                fq.push(f);
- 
 
 103 
                fq.push(f2);
- 
 
 104 
            }
- 
 
 105 
        }
- 
 
 106 
    }
- 
 
 107 
}
- 
 
 108 
 
- 
 
 109 
FaceID merge_face_in_1_ring(Manifold& m, VertexID v)
- 
 
 110 
{
- 
 
 111 
    Walker w=m.walker(v);
- 
 
 112 
    HalfEdgeID h = w.next().halfedge();
- 
 
 113 
    if(m.remove_vertex(v))
44 
    {
 114 
    {
45 
        cout
 - 
 
46 
        << " id : " << id
 - 
 
47 
        << " no vertices : " << no_vertices
 - 
 
48 
        << " avg pos : " << avg_pos
 - 
 
49 
        << " rad : " << rad
 - 
 
50 
        << " length : " << length
 115 
        return m.close_hole(h);
51 
        << " # : " << components
 - 
 
52 
        << " F : " << fun_value << endl;
 - 
 
53 
    }
 116 
    }
- 
 
 117 
    return InvalidFaceID;
54 
};
 118 
}
55 
 
 119 
 
- 
 
 120 
 
56 
void compute_edge_weights(Manifold& m, HalfEdgeAttributeVector<double>& edge_weights, FaceAttributeVector<int>& included)
 121 
void compute_edge_weights(const Manifold& m, HalfEdgeAttributeVector<double>& edge_weights, FaceAttributeVector<int>& included)
57 
{
 122 
{
58 
    edge_weights = HalfEdgeAttributeVector<double>(m.allocated_halfedges(), 0);
 123 
    edge_weights = HalfEdgeAttributeVector<double>(m.allocated_halfedges(), 0);
59 
    for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f)
 124 
    for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f)
60 
        if(included[*f])
 125 
        if(included[*f])
61 
        {
 126 
        {
Line 73... Line 138...
73 
            }
 138 
            }
74 
        }
 139 
        }
75 
}
 140 
}
76 
 
 141 
 
77 
template<typename T>
 142 
template<typename T>
78 
void smooth_fun(Manifold& m,
 143 
void smooth_fun(const Manifold& m,
79 
                VertexAttributeVector<int>& nailed,
 144 
                VertexAttributeVector<int>& nailed,
80 
                VertexAttributeVector<T>& fun, int iter=100)
 145 
                VertexAttributeVector<T>& fun, int iter=100)
81 
{
 146 
{
82 
    HalfEdgeAttributeVector<double> edge_weights;
 147 
    HalfEdgeAttributeVector<double> edge_weights;
83 
    FaceAttributeVector<int> included(m.allocated_faces(),1);
 148 
    FaceAttributeVector<int> included(m.allocated_faces(),1);
Line 102... Line 167...
102 
                new_fun[*v] = fun[*v];
 167 
                new_fun[*v] = fun[*v];
103 
        fun = new_fun;
 168 
        fun = new_fun;
104 
    }
 169 
    }
105 
}
 170 
}
106 
 
 171 
 
107 
void segment_manifold(Manifold& m, HalfEdgeAttributeVector<int>& ls_id,
 172 
vector<VertexID> cut_mesh(Manifold& m, VertexAttributeVector<double>& fun, double cut_val, VertexAttributeVector<VertexStatus>& status)
108 
                      FaceAttributeVector<int>& segmentation,
 - 
 
109 
                      vector<vector<int>>& boundaries)
 - 
 
110 
{
 173 
{
111 
    segmentation = FaceAttributeVector<int>(m.no_faces(), -1);
 174 
    cout << "Cutting at " << cut_val << endl;
- 
 
 175
- 
 
 176 
    vector<VertexID> new_verts;
112 
    int SEG_NO = 0;
 177 
    vector<HalfEdgeID> hidvec;
113 
    for(auto f = m.faces_begin(); f != m.faces_end(); ++f)
 178 
    for(HalfEdgeID hid : m.halfedges())
114 
    {
 179 
    {
- 
 
 180 
        Walker w = m.walker(hid);
115 
        if (segmentation[*f] == -1)
 181 
        double bval = fun[w.vertex()];
- 
 
 182 
        double aval = fun[w.opp().vertex()];
- 
 
 183 
        if(aval<bval && aval < cut_val && cut_val <= bval)
116 
        {
 184 
        {
- 
 
 185 
            double t = (cut_val-aval)/(bval-aval);
- 
 
 186 
            Vec3d p = t*m.pos(w.vertex()) + (1.0-t)*m.pos(w.opp().vertex());
- 
 
 187 
            VertexID vnew = m.split_edge(hid);
117 
            queue<FaceID> q;
 188 
            m.pos(vnew) = p;
118 
            q.push(*f);
 189 
            status[vnew] = COOKED;
119 
            vector<int> bound(0);
 190 
            fun[vnew] = cut_val;
120 
            while (!q.empty())
 191 
            new_verts.push_back(vnew);
- 
 
 192 
        }
- 
 
 193 
    }
- 
 
 194
- 
 
 195 
    for(FaceID fid : m.faces())
- 
 
 196 
        for(Walker w = m.walker(fid);!w.full_circle(); w = w.next())
- 
 
 197 
            if(status[w.vertex()] == COOKED && same_level(fun[w.vertex()],cut_val))
121 
            {
 198 
            {
122 
                FaceID face = q.front();
 199 
                Walker w0 = w;
123 
                q.pop();
 200 
                w = w.next().next();
124 
                segmentation[face] = SEG_NO;
 201 
                do
125 
                for(Walker w = m.walker(face); !w.full_circle(); w = w.next())
 - 
 
126 
                {
 202 
                {
127 
                    if(ls_id[w.halfedge()] == -1)
 203 
                    if(status[w.vertex()] == COOKED &&
- 
 
 204 
                       w.next().halfedge() != w0.halfedge() &&
- 
 
 205 
                       same_level(fun[w.vertex()],cut_val))
128 
                    {
 206 
                    {
129 
                        FaceID fopp = w.opp().face();
 207 
                        m.split_face_by_edge(fid, w0.vertex(), w.vertex());
130 
                        if(segmentation[fopp] == -1) {
 - 
 
131 
                            q.push(fopp);
 - 
 
132 
                        }
 208 
                        break;
133 
                    }
 - 
 
134 
                    else {
 - 
 
135 
                        bound.push_back(ls_id[w.halfedge()]);
 - 
 
136 
                    }
 209 
                    }
137
210 
                    w = w.next();
138 
                }
 211 
                }
- 
 
 212 
                while(!w.full_circle());
- 
 
 213 
                break;
139 
            }
 214 
            }
140 
            sort(begin(bound), end(bound));
 - 
 
141 
            auto new_end = unique(begin(bound), end(bound));
 - 
 
142 
            boundaries.push_back(vector<int>(begin(bound),new_end));
 - 
 
143 
            cout << "Boundaries of seg# " << SEG_NO << " : ";
 - 
 
144 
            for(auto be:boundaries[SEG_NO])
 - 
 
145 
                cout << be << " ";
 - 
 
146 
            cout << endl;
 - 
 
147 
            SEG_NO += 1;
 - 
 
148 
        }
 - 
 
149
- 
 
150 
    }
 - 
 
151 
}
 - 
 
152 
 
 - 
 
153 
 
 - 
 
154 
void shortest_edge_triangulate_face(Manifold& m, FaceID f0, VertexAttributeVector<int>& level_set_id_vertex)
 - 
 
155 
{
 - 
 
156 
    queue<FaceID> face_queue;
 - 
 
157
- 
 
158 
    face_queue.push(f0);
 - 
 
159
- 
 
160 
    while(!face_queue.empty())
 - 
 
161 
    {
 - 
 
162 
        FaceID f = face_queue.front();
 - 
 
163 
        face_queue.pop();
 - 
 
164
- 
 
165 
        // Create a vector of vertices.
 - 
 
166 
        vector<VertexID> verts;
 215 
    return new_verts;
167 
        for(Walker w = m.walker(f); !w.full_circle(); w = w.circulate_face_ccw())
 - 
 
168 
        {
 - 
 
169 
            FaceID fa = w.face();
 - 
 
170 
            FaceID fb = f;
 - 
 
171 
            assert(fa==fb);
 - 
 
172 
            verts.push_back(w.vertex());
 - 
 
173 
        }
 - 
 
174 
        // If there are just three we are done.
 - 
 
175 
        if(verts.size() == 3) continue;
 - 
 
176
- 
 
177 
        // Find vertex pairs that may be connected.
 - 
 
178 
        vector<pair<int,int> > vpairs;
 - 
 
179 
        const int N = verts.size();
 - 
 
180 
        for(int i = 0; i < N - 2; ++i){
 - 
 
181 
            for(int j = i + 2; j < N; ++j){
 - 
 
182 
                if(verts[i] != verts[j] &&
 - 
 
183 
                   !connected(m, verts[i], verts[j]) &&
 - 
 
184 
                   (level_set_id_vertex[verts[i]] == 0 || level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[j]]) &&
 - 
 
185 
                   valency(m, verts[i])<5 && valency(m,verts[j])<5
 - 
 
186 
                   )
 - 
 
187 
                    vpairs.push_back(pair<int,int>(i, j));
 - 
 
188 
            }
 - 
 
189 
        }
 - 
 
190 
        if(vpairs.empty()){
 - 
 
191 
            cout << "Warning: could not triangulate a face."
 - 
 
192 
            << "Probably a vertex appears more than one time in other vertex's one-ring" << endl;
 - 
 
193 
            continue;
 - 
 
194 
        }
 - 
 
195
- 
 
196 
        /* For all vertex pairs, find the edge lengths. Combine the
 - 
 
197 
         vertices forming the shortest edge. */
 - 
 
198
- 
 
199 
        float min_len=FLT_MAX;
 - 
 
200 
        int min_k = -1;
 - 
 
201 
        for(size_t k = 0; k < vpairs.size(); ++k){
 - 
 
202 
            int i = vpairs[k].first;
 - 
 
203 
            int j = vpairs[k].second;
 - 
 
204 
            float len = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));
 - 
 
205
- 
 
206 
            if(len<min_len){
 - 
 
207 
                min_len = len;
 - 
 
208 
                min_k = k;
 - 
 
209 
            }
 - 
 
210 
        }
 - 
 
211 
        assert(min_k != -1);
 - 
 
212
- 
 
213 
        {
 - 
 
214 
            // Split faces along edge whose midpoint is closest to isovalue
 - 
 
215 
            int i = vpairs[min_k].first;
 - 
 
216 
            int j = vpairs[min_k].second;
 - 
 
217 
            FaceID f_new = m.split_face_by_edge(f, verts[i], verts[j]);
 - 
 
218
- 
 
219 
            if(no_edges(m, f)>3)
 - 
 
220 
                face_queue.push(f);
 - 
 
221 
            if(no_edges(m, f_new)>3)
 - 
 
222 
                face_queue.push(f_new);
 - 
 
223 
        }
 - 
 
224
- 
 
225 
    }
 - 
 
226 
}
 216 
}
227 
 
 217 
 
228 
 
 - 
 
229 
void shortest_edge_split_face(Manifold& m, FaceID f0, VertexAttributeVector<int>& level_set_id_vertex)
 218 
void smooth_cooked_vertices(Manifold& m, const vector<VertexID>& vertices, VertexAttributeVector<double>& fun,
- 
 
 219 
                            VertexAttributeVector<VertexStatus>& status, int iter=1)
230 
{
 220 
{
231 
    queue<FaceID> face_queue;
 - 
 
232
- 
 
233 
    face_queue.push(f0);
 - 
 
234
- 
 
235 
    while(!face_queue.empty())
 221 
    for(int i=0;i<iter;++i)
236 
    {
 222 
    {
237 
        FaceID f = face_queue.front();
 223 
        VertexAttributeVector<Vec3d> npos(m.allocated_vertices(), Vec3d(0));
238 
        face_queue.pop();
 224 
        VertexAttributeVector<int> cnt(m.allocated_vertices(), 0);
239
- 
 
240 
        // Create a vector of vertices.
 225 
        for(auto vid : vertices)
241 
        vector<VertexID> verts;
 226 
           if(m.in_use(vid))
242 
        for(Walker w = m.walker(f); !w.full_circle(); w = w.circulate_face_ccw())
 227 
            for(Walker w = m.walker(vid); !w.full_circle(); w = w.circulate_vertex_ccw())
243 
        {
 - 
 
244 
            verts.push_back(w.vertex());
 - 
 
245 
        }
 - 
 
246
- 
 
247 
        // Find vertex pairs that may be connected.
 - 
 
248 
        vector<pair<int,int> > vpairs;
 - 
 
249 
        const int N = verts.size();
 - 
 
250 
        for(int i = 0; i < N ; ++i){
 - 
 
251 
            for(int j = 3; j < N-2; ++j){
 - 
 
252 
                int jj = (j+i)%N;
 228 
                if(status[vid]==COOKED) {
253 
                if(!connected(m, verts[i], verts[jj]) &&
 - 
 
254 
                   (level_set_id_vertex[verts[i]] != level_set_id_vertex[verts[jj]]) &&
 - 
 
255 
                   (level_set_id_vertex[verts[(i+1)%N]] == level_set_id_vertex[verts[i]]) &&
 - 
 
256 
                   (level_set_id_vertex[verts[(i+N-1)%N]] == level_set_id_vertex[verts[i]]) &&
 - 
 
257 
                   (level_set_id_vertex[verts[(jj+1)%N]] == level_set_id_vertex[verts[jj]]) &&
 - 
 
258 
                   (level_set_id_vertex[verts[(jj+N-1)%N]] == level_set_id_vertex[verts[jj]]))
 - 
 
259 
                    vpairs.push_back(pair<int,int>(i, jj));
 229 
                    npos[vid] += m.pos(w.vertex());
260 
            }
 - 
 
261 
        }
 - 
 
262 
        if(vpairs.empty()) continue;
 - 
 
263
- 
 
264
- 
 
265 
        /* For all vertex pairs, find the edge lengths. Combine the
 - 
 
266 
         vertices forming the shortest edge. */
 - 
 
267
- 
 
268 
        float min_len=FLT_MAX;
 - 
 
269 
        int min_k = -1;
 - 
 
270 
        for(size_t k = 0; k < vpairs.size(); ++k){
 - 
 
271 
            int i = vpairs[k].first;
 - 
 
272 
            int j = vpairs[k].second;
 - 
 
273 
            float len = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));
 - 
 
274
- 
 
275 
            if(len<min_len){
 - 
 
276 
                min_len = len;
 230 
                    ++cnt[vid];
277 
                min_k = k;
 - 
 
278 
            }
 231 
                }
279 
        }
 - 
 
280 
        assert(min_k != -1);
 232 
        for(auto vid : vertices)
281
- 
 
282 
        {
 - 
 
283 
            // Split faces along edge whose midpoint is closest to isovalue
 - 
 
284 
            int i = vpairs[min_k].first;
 - 
 
285 
            int j = vpairs[min_k].second;
 - 
 
286 
            FaceID f_new = m.split_face_by_edge(f, verts[i], verts[j]);
 - 
 
287 
            if(no_edges(m, f)>5)
 233 
            if(m.in_use(vid))
288 
                face_queue.push(f);
 - 
 
289 
            if(no_edges(m, f_new)>5)
 - 
 
290 
                face_queue.push(f_new);
 234 
            m.pos(vid) = (npos[vid])/cnt[vid];
291 
        }
 - 
 
292
- 
 
293 
    }
 235 
    }
294 
}
 236 
}
295 
 
 237 
 
296 
 
 - 
 
297 
 
 - 
 
298 
struct EdgeQElem {
 238 
int ssb(double a, double b)
299 
    float priority;
 - 
 
300 
    HalfEdgeID he;
 - 
 
301 
    EdgeQElem(float _priority, HalfEdgeID _he): priority(_priority), he(_he) {}
 - 
 
302 
};
 - 
 
303 
 
 - 
 
304 
bool operator<(const EdgeQElem& l, const EdgeQElem& r)
 - 
 
305 
{
 239 
{
- 
 
 240 
    const double THRESH = 1e-6;
- 
 
 241 
    double delta = a-b;
- 
 
 242 
    if(delta<-THRESH) return -1;
306 
    return l.priority < r.priority;
 243 
    if(delta>THRESH) return 1;
- 
 
 244 
    return 0;
307 
}
 245 
}
308 
 
 246 
 
309 
class FunctionalDifference: public EnergyFun
 - 
 
310 
{
 - 
 
311 
    VertexAttributeVector<float>& fun;
 - 
 
312 
    VertexAttributeVector<int>& status;
 - 
 
313 
public:
 - 
 
314 
    FunctionalDifference(VertexAttributeVector<float>& _fun, VertexAttributeVector<int>& _status): fun(_fun), status(_status) {}
 247 
void classify_vertices(Manifold& m, VertexAttributeVector<double>& fun, VertexAttributeVector<VertexStatus>& status)
315 
    virtual double delta_energy(const Manifold& m, HalfEdgeID h) const
 - 
 
316 
    {
 - 
 
317 
        Walker w = m.walker(h);
 - 
 
318 
        if(status[w.vertex()] == 1 && status[w.opp().vertex()]==1)
 - 
 
319 
            return DBL_MAX;
 - 
 
320 
        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()]));
 - 
 
321 
    }
 - 
 
322 
};
 - 
 
323 
 
 - 
 
324 
class TriangleQualityValence: public EnergyFun
 - 
 
325 
{
 248 
{
326 
    VertexAttributeVector<int>& idv;
 249 
    for(auto vid : m.vertices())
327 
    MinAngleEnergy mae;
 - 
 
328 
    ValencyEnergy vae;
 - 
 
329 
public:
 - 
 
330 
    TriangleQualityValence(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}
 - 
 
331 
    virtual double delta_energy(const Manifold& m, HalfEdgeID h) const
 - 
 
332 
    {
 250 
    {
333 
        Walker w = m.walker(h);
 251 
        double val = fun[vid];
334 
        if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] ||
 - 
 
335 
           idv[w.vertex()] == idv[w.opp().vertex()])
 - 
 
336 
            return DBL_MAX;
 - 
 
337
- 
 
338 
        VertexID v1 = w.opp().vertex();
 - 
 
339 
        VertexID v2 = w.vertex();
 252 
        auto w = m.walker(vid);
340 
        VertexID vo1 = w.next().vertex();
 253 
        int above_below = ssb(fun[w.vertex()],val);
341 
        VertexID vo2 = w.opp().next().vertex();
 - 
 
342
- 
 
343 
        int val1  = valency(m, v1);
 - 
 
344 
        int val2  = valency(m, v2);
 254 
        int changes = 0;
345 
        int valo1 = valency(m, vo1);
 - 
 
346 
        int valo2 = valency(m, vo2);
 255 
        w = w.circulate_vertex_ccw();
347
- 
 
348 
        // The optimal valency is four for a boundary vertex
 256 
        for(; !w.full_circle();w = w.circulate_vertex_ccw())
349 
        // and six elsewhere.
 - 
 
350 
        int t1 = boundary(m, v1) ? 4 : 6;
 - 
 
351 
        int t2 = boundary(m, v2) ? 4 : 6;
 - 
 
352 
        int to1 = boundary(m, vo1) ? 4 : 6;
 - 
 
353 
        int to2 = boundary(m, vo2) ? 4 : 6;
 - 
 
354
- 
 
355 
        int before =
 - 
 
356 
        max(max(sqr(val1-t1),sqr(val2-t2)), max(sqr(valo1-to1), sqr(valo2-to2)));
 - 
 
357 
        int after =
 - 
 
358 
        max(max(sqr(valo1+1-to1),sqr(val1-1-t1)), max(sqr(val2-1-t2),sqr(valo2+1-to2)));
 - 
 
359
- 
 
360 
        return static_cast<double>(after-before);
 - 
 
361
257 
        {
362 
        //        return vae.delta_energy(m,h);
 - 
 
363 
        float la = length(m.pos(w.next().vertex())-m.pos(w.opp().next().vertex()));
 - 
 
364 
        float lb = length(m.pos(w.vertex())-m.pos(w.opp().vertex()));
 258 
            int above_below_new = ssb(fun[w.vertex()],val);
365 
        return la-lb;
 - 
 
366 
        return mae.delta_energy(m,h);
 - 
 
367 
    }
 - 
 
368 
};
 - 
 
369 
 
 - 
 
370 
class TriangleQuality: public EnergyFun
 - 
 
371 
{
 - 
 
372 
    VertexAttributeVector<int>& idv;
 - 
 
373 
    MinAngleEnergy mae;
 - 
 
374 
    ValencyEnergy vae;
 - 
 
375 
public:
 - 
 
376 
    TriangleQuality(VertexAttributeVector<int>& _idv): idv(_idv), mae(-1) {}
 - 
 
377 
    virtual double delta_energy(const Manifold& m, HalfEdgeID h) const
 - 
 
378 
    {
 - 
 
379 
        Walker w = m.walker(h);
 - 
 
380 
        if(idv[w.next().vertex()] == idv[w.opp().next().vertex()] ||
 - 
 
381 
           idv[w.vertex()] == idv[w.opp().vertex()])
 259 
            if(above_below != above_below_new)
382 
            return DBL_MAX;
 260 
                ++changes;
383 
        return mae.delta_energy(m,h);
 - 
 
384 
    }
 - 
 
385 
};
 - 
 
386 
 
 - 
 
387 
Vec3d grad(HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun, HMesh::FaceID f)
 - 
 
388 
{
 - 
 
389 
    if(no_edges(m,f) != 3)
 - 
 
390 
        return Vec3d(0.0);
 - 
 
391
- 
 
392
- 
 
393 
    Vec3d n = normal(m, f);
 - 
 
394
- 
 
395 
    Vec3d gsum(0.0);
 - 
 
396 
    for(Walker w = m.walker(f); !w.full_circle(); w = w.next())
 - 
 
397 
    {
 - 
 
398 
        VertexID avid = w.opp().vertex();
 - 
 
399 
        VertexID bvid = w.vertex();
 - 
 
400 
        VertexID cvid = w.next().vertex();
 - 
 
401 
        Vec3d gdir = normalize(cross(n, m.pos(bvid) - m.pos(avid)));
 - 
 
402 
        double l = dot(m.pos(cvid)-m.pos(avid), gdir);
 - 
 
403 
        gdir *= fun[cvid]/l;
 - 
 
404 
        gsum += gdir;
 - 
 
405 
    }
 - 
 
406 
    return gsum;
 - 
 
407 
}
 - 
 
408 
 
 - 
 
409 
Vec2d extend_fun2(HMesh::Manifold& m, HMesh::HalfEdgeID h,
 - 
 
410 
                  HMesh::VertexAttributeVector<double>& fun,
 - 
 
411 
                  HMesh::VertexAttributeVector<Vec2d>& fun2)
 - 
 
412 
{
 - 
 
413 
    Walker w = m.walker(h);
 - 
 
414 
    FaceID f = w.face();
 - 
 
415 
    Vec3d a = m.pos(w.opp().vertex());
 - 
 
416 
    Vec3d b = m.pos(w.vertex());
 - 
 
417 
    Vec3d c = m.pos(w.next().vertex());
 - 
 
418 
    Vec3d g = grad(m, fun, f);
 - 
 
419 
    Vec3d n = normal(m, f);
 - 
 
420 
    Vec3d N = normalize(cross(g, n));
 - 
 
421 
    float dot_nba = dot(N,b-a);
 - 
 
422 
    Vec2d cval;
 - 
 
423 
    if(dot_nba > 1e-5)
 - 
 
424 
    {
 - 
 
425 
        float t = dot(N, c-a)/dot_nba;
 - 
 
426 
        if(t>0 && t < 1) {
 - 
 
427 
            Vec2d aval = fun2[w.opp().vertex()];
 - 
 
428 
            Vec2d bval = fun2[w.vertex()];
 261 
            above_below=above_below_new;
429 
            return (1-t)*aval + t*bval;
 - 
 
430 
        }
 262 
        }
- 
 
 263 
        if(changes==0)
- 
 
 264 
            status[vid]=POLE;
- 
 
 265 
        else if(changes>=4)
- 
 
 266 
            status[vid]=SADDLE;
- 
 
 267 
        else if(changes==2)
- 
 
 268 
            status[vid]=RAW;
- 
 
 269 
        else assert(0);
431 
    }
 270 
    }
432 
    return Vec2d(0);
 - 
 
433 
}
 271 
}
434 
 
 272 
 
435 
 
 - 
 
436 
double solve_for_orthogonal_gradients(HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,
 273 
void remove_vertices(Manifold& m, VertexAttributeVector<double>& fun, VertexAttributeVector<VertexStatus>& status, double rmin, double rmax)
437 
                                      HalfEdgeID h, double beta, double gamma)
 - 
 
438 
{
 274 
{
439 
    Vec3d n = normal(m, m.walker(h).face());
 275 
    vector<VertexID> vertices_to_push;
440
276 
    for(auto vid: m.vertices())
441 
    Vec3d uvw[3];
 277 
        if(status[vid]==RAW &&
442 
    Vec3d abc;
 278 
           fun[vid] >= rmin &&
443 
    Walker w = m.walker(h);
 279 
           fun[vid] <= rmax)
444 
    VertexID vid[] = {w.next().vertex(), w.opp().vertex(),w.vertex()};
 280 
            vertices_to_push.push_back(vid);
445
281
446 
    for(int i=0; !w.full_circle(); w = w.next(),++  i)
 - 
 
447 
    {
 - 
 
448 
        uvw[i] = normalize(cross(n, m.pos(vid[(i+1)%3]) - m.pos(vid[(i+2)%3])));
 - 
 
449 
        uvw[i] /= dot(m.pos(vid[i])-m.pos(vid[(i+1)%3]), uvw[i]);
 282 
    random_shuffle(begin(vertices_to_push), end(vertices_to_push));
450 
        abc[i] = fun[vid[i]];
 - 
 
451 
    }
 - 
 
452 
    Vec3d A(dot(uvw[0],uvw[0]),dot(uvw[0],uvw[1]),dot(uvw[0],uvw[2]));
 - 
 
453 
    Vec3d B(dot(uvw[1],uvw[0]),dot(uvw[1],uvw[1]),dot(uvw[1],uvw[2]));
 - 
 
454 
    Vec3d C(dot(uvw[2],uvw[0]),dot(uvw[2],uvw[1]),dot(uvw[2],uvw[2]));
 - 
 
455
283
456 
    return - (beta * dot(B, abc) + gamma * dot(C, abc))/dot(A, abc);
 - 
 
457 
}
 - 
 
458 
 
 - 
 
459 
 
 - 
 
460 
void recompute_circle_param(Manifold& m, VertexID v, VertexAttributeVector<double>& fun, VertexAttributeVector<Vec2d>& circle_pos, double f_origin)
 - 
 
461 
{
 - 
 
462 
    Vec2d c_new(0);
 - 
 
463 
    for(auto w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_ccw())
 284 
    for(auto vid : vertices_to_push)
464 
    {
 285 
    {
465 
        Walker w_base = w.next();
 - 
 
466 
        Vec2d X = circle_pos[w_base.opp().vertex()];
 286 
        FaceID fid = merge_face_in_1_ring(m,vid);
467 
        Vec2d Y(-X[1], X[0]);
 287 
        if(fid != InvalidFaceID)
468 
        Vec2d vec = circle_pos[w_base.vertex()];
 - 
 
469 
        double angle = atan2(max(-1.0,min(1.0,dot(vec,Y))),
 - 
 
470 
                             max(-1.0,min(1.0,dot(vec,X))));
 - 
 
471 
        double new_angle = solve_for_orthogonal_gradients(m, fun, w_base.halfedge(), 0, angle);
 - 
 
472 
        new_angle = max(-M_PI/4, min(M_PI/4,new_angle));
 - 
 
473 
        c_new += cos(new_angle)*X + sin(new_angle)*Y;
 288 
            high_gradient_triangulate(m, fid, fun);
474 
    }
 289 
    }
475 
    circle_pos[v] = normalize(c_new);
 - 
 
476 
}
 290 
}
477 
 
 291 
 
478 
template<typename T>
 - 
 
479 
void orthogonal_trajectories(Manifold& m, VertexID v, VertexAttributeVector<double>& f, double min_length_fun,
 - 
 
480 
                             VertexAttributeVector<T>& h_in, VertexAttributeVector<T>& h_out)
 - 
 
481 
{
 - 
 
482 
    h_out[v] = T(0);
 - 
 
483 
 
 - 
 
484 
    T num(0);
 - 
 
485 
    double denom = 0.0001;
 - 
 
486 
    for(Walker w = m.walker(v); !w.full_circle(); w = w.circulate_vertex_ccw())
 - 
 
487 
    {
 - 
 
488 
        double weight = area(m, w.face());
 - 
 
489 
 
 - 
 
490 
        Vec3d a = m.pos(v);
 - 
 
491 
        Vec3d b = m.pos(w.vertex());
 - 
 
492 
        Vec3d c = m.pos(w.next().vertex());
 - 
 
493
- 
 
494
- 
 
495 
        Vec3d fvec(f[v],f[w.vertex()], f[w.next().vertex()]);
 - 
 
496 
        if(sqr_length(fvec)>0)
 - 
 
497 
            fvec.normalize();
 - 
 
498 
        double f_a = fvec[0];
 - 
 
499 
        double f_b = fvec[1];
 - 
 
500 
        double f_c = fvec[2];
 - 
 
501
- 
 
502 
        T h_b = h_in[w.vertex()];
 - 
 
503 
        T h_c = h_in[w.next().vertex()];
 - 
 
504
- 
 
505 
        Vec3d e_a = c-b;
 - 
 
506 
        Vec3d e_b = a-c;
 - 
 
507 
        Vec3d e_c = b-a;
 - 
 
508
- 
 
509 
        Vec3d N = normalize(cross(e_b,e_c)+cross(e_c,e_a)+cross(e_a,e_b));
 - 
 
510 
        Vec3d g_a = cross(N, normalize(e_a));
 - 
 
511 
        g_a /= dot(g_a, e_b);
 - 
 
512 
 
 - 
 
513 
        Vec3d g_b = cross(N, normalize(e_b));
 - 
 
514 
        g_b /= dot(g_b, e_c);
 - 
 
515 
 
 - 
 
516 
        Vec3d g_c = cross(N, normalize(e_c));
 - 
 
517 
        g_c /= dot(g_c, e_a);
 - 
 
518
- 
 
519 
        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)) +
 - 
 
520 
                         h_c * (f_a*dot(g_a, g_c) + f_b*dot(g_b,g_c) + f_c*dot(g_c, g_c)));
 - 
 
521 
        denom += weight * (f_a * dot(g_a,g_a)+ f_b*dot(g_a,g_b) +f_c*dot(g_a, g_c));
 - 
 
522 
 
 - 
 
523 
    }
 - 
 
524 
    h_out[v] = num / denom;
 - 
 
525 
}
 - 
 
526 
 
 292 
 
527 
void polarize_mesh(Manifold& m, VertexAttributeVector<double>& fun, double vmin, double vmax, const int divisions, VertexAttributeVector<Vec2d>& parametrization)
 293 
void polarize_mesh(Manifold& m, VertexAttributeVector<double>& fun, double vmin, double vmax, const int divisions, VertexAttributeVector<Vec2d>& parametrization)
528 
{
 294 
{
529 
    double interval = (vmax-vmin)/(divisions+1);
 - 
 
530
- 
 
531 
    VertexAttributeVector<int> status(m.allocated_vertices(), 0);
 295 
    VertexAttributeVector<VertexStatus> status(m.allocated_vertices(), RAW);
532
- 
 
533
- 
 
534 
    // ----------------------------------------
 - 
 
535 
    cout << "Tracing level set curves" << endl;
 - 
 
536
296
- 
 
 297 
    classify_vertices(m, fun, status);
537
298
538 
    for(int cut_no=1; cut_no<divisions+1;++cut_no)
 299 
    double interval = (vmax-vmin)/(divisions+1);
539 
    {
 - 
 
540 
        double cut_val = vmin + cut_no * interval;
 300 
    double base_cut_val = vmin + 0.02*(vmax-vmin);
541
- 
 
542 
        vector<HalfEdgeID> hidvec;
 - 
 
543 
        for(HalfEdgeID hid : m.halfedges())
 - 
 
544 
        {
 - 
 
545 
            Walker w = m.walker(hid);
 - 
 
546 
            double bval = fun[w.vertex()];
 - 
 
547 
            double aval = fun[w.opp().vertex()];
 - 
 
548 
            if(aval<bval && aval <= cut_val && cut_val < bval)
 - 
 
549 
            {
 - 
 
550 
                double t = (cut_val-aval)/(bval-aval);
 - 
 
551 
                Vec3d p = t*m.pos(w.vertex()) + (1.0-t)*m.pos(w.opp().vertex());
 - 
 
552 
                VertexID vnew = m.split_edge(hid);
 - 
 
553 
                m.pos(vnew) = p;
 - 
 
554 
                status[vnew] = 1;
 - 
 
555 
                fun[vnew] = cut_val;
 - 
 
556 
            }
 - 
 
557 
        }
 - 
 
558
- 
 
559 
        for(FaceID fid : m.faces())
 - 
 
560 
                for(Walker w = m.walker(fid);!w.full_circle(); w = w.next())
 - 
 
561 
                    if(status[w.vertex()] == 1 && same_level(fun[w.vertex()],cut_val))
 - 
 
562 
                    {
 - 
 
563 
                        Walker w0 = w;
 - 
 
564 
                        w = w.next().next();
 - 
 
565 
                        do
 - 
 
566 
                        {
 - 
 
567 
                            if(status[w.vertex()] == 1 &&
 - 
 
568 
                               w.next().halfedge() != w0.halfedge() &&
 - 
 
569 
                               same_level(fun[w0.vertex()],cut_val))
 - 
 
570 
                            {
 - 
 
571 
                                m.split_face_by_edge(fid, w0.vertex(), w.vertex());
 - 
 
572 
                                break;
 - 
 
573 
                            }
 - 
 
574 
                            w = w.next();
 - 
 
575 
                        }
 - 
 
576 
                        while(!w.full_circle());
 - 
 
577 
                        break;
 - 
 
578 
                    }
 - 
 
579 
    }
 - 
 
580 
 
 - 
 
581 
    shortest_edge_triangulate(m);
 - 
 
582
- 
 
583 
 
 - 
 
584
301
585 
    // ----------------------------
 - 
 
586 
    cout << "Numbering the level sets" << endl;
 302 
    vector<VertexID> new_verts = cut_mesh(m, fun, base_cut_val, status);
587 
    vector<LevelSetInfo> level_set_info;
 303 
    smooth_cooked_vertices(m, new_verts, fun, status, 1);
588 
    HalfEdgeAttributeVector<int> level_set_id(m.allocated_halfedges(), -1);
 304 
    cout << "Base " << base_cut_val << endl;
589 
    VertexAttributeVector<int> level_set_id_vertex(m.allocated_vertices(), -1);
 - 
 
590 
    int no_id=0;
 - 
 
591 
    for(HalfEdgeID hid : m.halfedges())
 305 
    for(int iter=0;iter<divisions+2;++iter)
592 
    {
 306 
    {
593 
        Walker w = m.walker(hid);
 - 
 
594 
        if(status[w.vertex()] == 1 && status[w.opp().vertex()] == 1 &&
 - 
 
595 
           same_level(fun[w.vertex()], fun[w.opp().vertex()]) &&
 - 
 
596 
           level_set_id[w.halfedge()] == -1)
 307 
        int _iter=iter;//%2==0?iter/2:-1-iter/2;
597 
        {
 - 
 
598 
            LevelSetInfo lsi;
 - 
 
599 
            lsi.id = no_id;
 - 
 
600 
            lsi.fun_value = fun[w.vertex()];
 - 
 
601 
            lsi.components = 1;
 - 
 
602 
            lsi.no_vertices = 0;
 - 
 
603 
            double level_set_length = 0;
 308 
        double cut_val = base_cut_val + interval * _iter;
604 
            while(level_set_id[w.halfedge()] != no_id)
 309 
        cout << "i,c " << _iter << "," << cut_val << endl;
605 
            {
 - 
 
606 
                level_set_length += length(m,w.halfedge());
 310 
        remove_vertices(m, fun, status, cut_val, cut_val+interval);
607 
                level_set_id[w.halfedge()] = no_id;
 - 
 
608 
                level_set_id[w.opp().halfedge()] = no_id;
 - 
 
609 
                level_set_id_vertex[w.vertex()] = no_id;
 - 
 
610 
                w = w.next();
 - 
 
611 
                while(!(status[w.vertex()] == 1 && same_level(fun[w.vertex()], fun[w.opp().vertex()])))
 311 
        vector<VertexID> new_verts = cut_mesh(m, fun,_iter<0?cut_val:cut_val+interval, status);
612 
                    w = w.circulate_vertex_cw();
 - 
 
613 
                lsi.no_vertices += 1;
 - 
 
614 
            }
 - 
 
615 
            lsi.h = w.halfedge();
 - 
 
616 
            lsi.length = level_set_length;
 - 
 
617 
            level_set_info.push_back(lsi);
 312 
        kill_wedges(m, fun, status);
618 
            level_set_info.back().print();
 - 
 
619 
            ++no_id;
 - 
 
620 
            assert(level_set_info.size()==no_id);
 313 
        smooth_cooked_vertices(m, new_verts, fun, status, 1);
621 
        }
 - 
 
622 
    }
 314 
    }
623
315
624 
    cout << "ids" << no_id << " " << level_set_info.size() << endl;
 - 
 
625 
    for(int i=0;i<level_set_info.size(); ++i)
 - 
 
626 
        for(int j=0;j<level_set_info.size(); ++j)
 - 
 
627 
            if(i != j && same_level(level_set_info[i].fun_value, level_set_info[j].fun_value))
 - 
 
628 
                level_set_info[i].components += 1;
 - 
 
629
- 
 
630 
    HalfEdgeID min_length_h;
 - 
 
631 
    double min_length_fun;
 - 
 
632 
    double min_length=FLT_MAX;
 - 
 
633 
    int min_length_id = -1;
 316 
    int C=0, R=0, P=0,S=0;
634 
    for(int i=0;i<level_set_info.size(); ++i)
 317 
    for(auto vid : m.vertices())
635 
    {
 318 
    {
636 
        LevelSetInfo& lsi = level_set_info[i];
 319 
        switch(status[vid]) {
637 
        if(lsi.components==1 && lsi.length<min_length)
 - 
 
638 
        {
 - 
 
639 
            min_length_h = lsi.h;
 320 
            case RAW: ++R; break;
640 
            min_length_fun = lsi.fun_value;
 321 
            case COOKED: ++C; break;
641 
            min_length = lsi.length;
 322 
            case POLE: ++P; break;
642 
            min_length_id = lsi.id;
 323 
            case SADDLE: ++S; break;
643
- 
 
644 
        }
 324 
        }
645 
    }
 325 
    }
646
- 
 
647 
    // ----------------------------
 - 
 
648 
    cout << "Remove vertices not on level set curves" << endl;
 326 
    cout << "R=" << R << " C=" << C << " P=" << P << " S=" << S << endl;
649
- 
 
650
- 
 
651 
}
 327 
}
652 
 
 328 
 
653 
void make_height_fun(const HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,
 329 
void make_height_fun(const HMesh::Manifold& m, HMesh::VertexAttributeVector<double>& fun,
654 
                     double& vmin, double& vmax)
 330 
                     double& vmin, double& vmax)
655 
{
 331 
{
656 
    VertexIDIterator vid = m.vertices_begin();
 332 
    VertexIDIterator vid = m.vertices_begin();
- 
 
 333 
    VertexAttributeVector<int> nailed(m.allocated_vertices(), 0);
- 
 
 334 
    smooth_fun(m, nailed, fun,0);
657 
    vmin = vmax = m.pos(*vid)[1];
 335 
    vmin = vmax = m.pos(*vid)[1];
658 
    for(; vid != m.vertices_end(); ++vid)
 336 
    for(; vid != m.vertices_end(); ++vid)
659 
    {
 337 
    {
660 
        double v = dot(m.pos(*vid),Vec3d(0.0,1,0.00));
 338 
        double v = dot(m.pos(*vid),Vec3d(0.0,1,0.00));
661 
        fun[*vid] = v;
 339 
        fun[*vid] = v;
662 
        vmin = min(v, vmin);
 340 
        vmin = min(v, vmin);
663 
        vmax = max(v, vmax);
 341 
        vmax = max(v, vmax);
664 
    }
 342 
    }
- 
 
 343
- 
 
 344 
}
- 
 
 345 
 
- 
 
 346 
void make_adf_fun(HMesh::Manifold& m, double t, HMesh::VertexAttributeVector<double>& F,
- 
 
 347 
                     double& vmin, double& vmax)
- 
 
 348 
{
- 
 
 349 
    Harmonics harm(m);
- 
 
 350 
    harm.compute_adf(F, t);
- 
 
 351 
    VertexAttributeVector<int> nailed(m.allocated_vertices(), 0);
- 
 
 352 
    vmin = vmax = F[*m.vertices_begin()];
- 
 
 353 
    for(auto vid : m.vertices())
- 
 
 354 
    {
- 
 
 355 
        vmin = min(F[vid], vmin);
- 
 
 356 
        vmax = max(F[vid], vmax);
- 
 
 357 
    }
- 
 
 358 
    cout << vmin << vmax << endl;
- 
 
 359
665 
}
 360 
}
666 
 
 361