Branch data Line data Source code
1 : : // $Id: quadrics.cxx,v 1.5 1997/10/01 14:07:28 garland Exp $
2 : :
3 : : //#include "qslim.h"
4 : : #include <assert.h>
5 : : #include "quadrics.h"
6 : : #include "3D.h"
7 : : #include "primitives.h"
8 : :
9 : : ////////////////////////////////////////////////////////////////////////
10 : : //
11 : : // Primitive quadric construction and evaluation routines
12 : : //
13 : :
14 : : //
15 : : // Construct a quadric to evaluate the squared distance of any point
16 : : // to the given point v. Naturally, the iso-surfaces are just spheres
17 : : // centered at v.
18 : : //
19 : 0 : Mat4 quadrix_vertex_constraint(const Vec3& v)
20 : : {
21 : 0 : Mat4 L(Mat4::identity);
22 : :
23 : 0 : L(0,3) = -v[0];
24 : 0 : L(1,3) = -v[1];
25 : 0 : L(2,3) = -v[2];
26 : 0 : L(3,3) = v*v;
27 : :
28 : 0 : L(3,0) = L(0,3);
29 : 0 : L(3,1) = L(1,3);
30 : 0 : L(3,2) = L(2,3);
31 : :
32 : 0 : return L;
33 : : }
34 : 0 : Mat4 quadrix_vertex_constraint(moab::EntityHandle vert)
35 : : {
36 [ # # ]: 0 : Mat4 L(Mat4::identity);
37 : : double v[3];
38 [ # # ]: 0 : mb->get_coords(&vert, 1 , v);
39 [ # # ]: 0 : L(0,3) = -v[0];
40 [ # # ]: 0 : L(1,3) = -v[1];
41 [ # # ]: 0 : L(2,3) = -v[2];
42 [ # # ]: 0 : L(3,3) = v[0]*v[0] + v[1]*v[1] + v[2]*v[2] ; // v*v;
43 : :
44 [ # # ][ # # ]: 0 : L(3,0) = L(0,3);
45 [ # # ][ # # ]: 0 : L(3,1) = L(1,3);
46 [ # # ][ # # ]: 0 : L(3,2) = L(2,3);
47 : :
48 : 0 : return L;
49 : : }
50 : : //
51 : : // Construct a quadric to evaluate the squared distance of any point
52 : : // to the given plane [ax+by+cz+d = 0]. This is the "fundamental error
53 : : // quadric" discussed in the paper.
54 : : //
55 : 10400 : Mat4 quadrix_plane_constraint(double a, double b, double c, double d)
56 : : {
57 : 10400 : Mat4 K(Mat4::zero);
58 : :
59 : 10400 : K(0,0) = a*a; K(0,1) = a*b; K(0,2) = a*c; K(0,3) = a*d;
60 : 10400 : K(1,0) =K(0,1); K(1,1) = b*b; K(1,2) = b*c; K(1,3) = b*d;
61 : 10400 : K(2,0) =K(0,2); K(2,1) =K(1,2); K(2,2) = c*c; K(2,3) = c*d;
62 : 10400 : K(3,0) =K(0,3); K(3,1) =K(1,3); K(3,2) =K(2,3);K(3,3) = d*d;
63 : :
64 : 10400 : return K;
65 : : }
66 : :
67 : : //
68 : : // Define some other convenient ways for constructing these plane quadrics.
69 : : //
70 : 400 : Mat4 quadrix_plane_constraint(const Vec3& n, double d)
71 : : {
72 : 400 : return quadrix_plane_constraint(n[X], n[Y], n[Z], d);
73 : : }
74 : :
75 : : //Mat4 quadrix_plane_constraint(Face& T)
76 : 10000 : Mat4 quadrix_plane_constraint(moab::EntityHandle triangle)
77 : : {
78 : : // const Plane& p = T.plane();
79 [ + - ]: 10000 : Plane p=trianglePlane(mb, triangle);
80 : : double a,b,c,d;
81 [ + - ]: 10000 : p.coeffs(&a, &b, &c, &d);
82 : :
83 [ + - ]: 10000 : return quadrix_plane_constraint(a, b, c, d);
84 : : }
85 : :
86 : 0 : Mat4 quadrix_plane_constraint(const Vec3& v1, const Vec3& v2, const Vec3& v3)
87 : : {
88 [ # # ]: 0 : Plane P(v1,v2,v3);
89 : : double a,b,c,d;
90 [ # # ]: 0 : P.coeffs(&a, &b, &c, &d);
91 : :
92 [ # # ]: 0 : return quadrix_plane_constraint(a, b, c, d);
93 : : }
94 : :
95 : 60616 : double quadrix_evaluate_vertex(const Vec3& v, const Mat4& K)
96 : : {
97 : 60616 : double x=v[X], y=v[Y], z=v[Z];
98 : :
99 : : #ifndef VECTOR_COST_EVALUATION
100 : : //
101 : : // This is the fast way of computing (v^T Q v).
102 : : //
103 : 60616 : return x*x*K(0,0) + 2*x*y*K(0,1) + 2*x*z*K(0,2) + 2*x*K(0,3)
104 : 60616 : + y*y*K(1,1) + 2*y*z*K(1,2) + 2*y*K(1,3)
105 : 60616 : + z*z*K(2,2) + 2*z*K(2,3)
106 : 60616 : + K(3,3);
107 : : #else
108 : : //
109 : : // The equivalent thing using matrix/vector operations.
110 : : // It's a lot clearer, but it's also slower.
111 : : //
112 : : Vec4 v2(x,y,z,1);
113 : : return v2*(K*v2);
114 : : #endif
115 : : }
116 : :
117 : :
118 : :
119 : : ////////////////////////////////////////////////////////////////////////
120 : : //
121 : : // Routines for computing discontinuity constraints
122 : : //
123 : :
124 : : //bool is_border(Edge *e )
125 : 15200 : bool is_border(moab::EntityHandle eh )
126 : : {
127 : : // use nb of triangles connected to check if number of tri is 1
128 [ + - ]: 15200 : std::vector<moab::EntityHandle> adjTri;
129 : : moab::ErrorCode rval = mb->get_adjacencies(&eh, 1, 2, false, adjTri,
130 [ + - ]: 15200 : moab::Interface::UNION);
131 [ + - ][ + + ]: 15200 : if (moab::MB_SUCCESS==rval && adjTri.size()==1)
[ + + ]
132 : 400 : return true;
133 : 15200 : return false;
134 : : }
135 : :
136 : 0 : bool check_for_discontinuity(moab::EntityHandle eh) //Edge *e)
137 : : {
138 : 0 : return is_border(eh);
139 : : }
140 : :
141 : 400 : Mat4 quadrix_discontinuity_constraint( moab::EntityHandle mbe
142 : : /*Edge *edge*/, const Vec3& n)
143 : : {
144 : : //Vec3& org = *edge->org();
145 : : //Vec3& dest = *edge->dest();
146 : : // to do: Vec3 origin of edge
147 : : // is the orientation important for an edge?
148 : : const moab::EntityHandle * conn;
149 : : int num_nodes;
150 [ + - ]: 400 : mb->get_connectivity(mbe, conn, num_nodes);
151 [ + - ]: 400 : Vec3 dest = getVec3FromMBVertex(mb, conn[1]); //
152 [ + - ]: 400 : Vec3 org = getVec3FromMBVertex(mb, conn[0]);
153 [ + - ]: 400 : Vec3 e = dest - org;
154 : :
155 [ + - ]: 400 : Vec3 n2 = e ^ n;
156 [ + - ]: 400 : unitize(n2);
157 : :
158 [ + - ][ + - ]: 400 : double d = -n2 * org;
159 [ + - ]: 400 : return quadrix_plane_constraint(n2, d);
160 : : }
161 : :
162 : :
163 : 400 : Mat4 quadrix_discontinuity_constraint(moab::EntityHandle mbedge)// Edge *edge)
164 : : {
165 [ + - ]: 400 : Mat4 D(Mat4::zero);
166 : :
167 : : /*
168 : : face_buffer& faces = edge->faceUses();
169 : :
170 : : for(int i=0; i<faces.length(); i++)
171 : : D += quadrix_discontinuity_constraint(edge,faces(i)->plane().normal());
172 : : */
173 : : // to do : to get the connected faces to the edge
174 : :
175 [ + - ]: 400 : std::vector<moab::EntityHandle> adjFaces;
176 : : moab::ErrorCode rval = mb->get_adjacencies(&mbedge, 1, 2, false, adjFaces,
177 [ + - ]: 400 : moab::Interface::UNION);
178 [ - + ][ # # ]: 400 : assert(rval == moab::MB_SUCCESS);
179 : :
180 [ + + ]: 800 : for (unsigned int i=0; i<adjFaces.size(); i++)
181 : : {
182 [ + - ]: 400 : moab::EntityHandle F = adjFaces[i];
183 [ + - ]: 400 : Plane p = trianglePlane( mb, F);
184 [ + - ][ + - ]: 400 : D += quadrix_discontinuity_constraint(mbedge, p.normal());
[ + - ]
185 : : }
186 : 400 : return D;
187 : : }
188 : :
189 : :
190 : :
191 : : ////////////////////////////////////////////////////////////////////////
192 : : //
193 : : // Routines for computing contraction target
194 : : //
195 : :
196 : 0 : bool quadrix_find_local_fit(const Mat4& K,
197 : : const Vec3& v1, const Vec3& v2,
198 : : Vec3& candidate)
199 : : {
200 : :
201 [ # # ][ # # ]: 0 : Vec3 v3 = (v1 + v2) / 2;
202 : :
203 : 0 : bool try_midpoint = opts.placement_policy > PLACE_ENDPOINTS;
204 : :
205 [ # # ]: 0 : double c1 = quadrix_evaluate_vertex(v1, K);
206 [ # # ]: 0 : double c2 = quadrix_evaluate_vertex(v2, K);
207 : : double c3;
208 [ # # ][ # # ]: 0 : if( try_midpoint ) c3 = quadrix_evaluate_vertex(v3, K);
209 : :
210 [ # # ]: 0 : if( c1<c2 )
211 : : {
212 [ # # ][ # # ]: 0 : if( try_midpoint && c3<c1 )
213 [ # # ]: 0 : candidate=v3;
214 : : else
215 [ # # ]: 0 : candidate=v1;
216 : : }
217 : : else
218 : : {
219 [ # # ][ # # ]: 0 : if( try_midpoint && c3<c2 )
220 [ # # ]: 0 : candidate=v3;
221 : : else
222 [ # # ]: 0 : candidate=v2;
223 : : }
224 : :
225 : 0 : return true;
226 : : }
227 : :
228 : 2 : bool quadrix_find_line_fit(const Mat4& Q,
229 : : const Vec3& v1, const Vec3& v2,
230 : : Vec3& candidate)
231 : : {
232 [ + - ]: 2 : Vec3 d = v1-v2;
233 : :
234 [ + - ]: 2 : Vec3 Qv2 = Q*v2;
235 [ + - ]: 2 : Vec3 Qd = Q*d;
236 : :
237 [ + - ][ + - ]: 2 : double denom = 2*d*Qd;
238 : :
239 [ - + ]: 2 : if( denom == 0.0 )
240 : 0 : return false;
241 : :
242 [ + - ][ + - ]: 2 : double a = (d*Qv2 + v2*Qd) / denom;
243 : :
244 [ + - ]: 2 : if( a<0.0 ) a=0.0;
245 [ - + ]: 2 : if( a>1.0 ) a=1.0;
246 : :
247 : :
248 [ + - ][ + - ]: 2 : candidate = a*d + v2;
[ + - ]
249 : 2 : return true;
250 : : }
251 : :
252 : 60616 : bool quadrix_find_best_fit(const Mat4& Q, Vec3& candidate)
253 : : {
254 [ + - ]: 60616 : Mat4 K = Q;
255 [ + - ][ + - ]: 60616 : K(3,0) = K(3,1) = K(3,2) = 0.0; K(3,3) = 1;
[ + - ][ + - ]
256 : :
257 : :
258 [ + - ]: 60616 : Mat4 M;
259 [ + - ]: 60616 : double det = K.inverse(M);
260 [ + - ][ + + ]: 60616 : if( FEQ(det, 0.0, 1e-12) )
261 : 2 : return false;
262 : :
263 : :
264 : : #ifdef SAFETY
265 : : //
266 : : // The homogeneous division SHOULDN'T be necessary.
267 : : // But, when we're being SAFE, we do it anyway just in case.
268 : : //
269 : : candidate[X] = M(0,3)/M(3,3);
270 : : candidate[Y] = M(1,3)/M(3,3);
271 : : candidate[Z] = M(2,3)/M(3,3);
272 : : #else
273 [ + - ][ + - ]: 60614 : candidate[X] = M(0,3);
274 [ + - ][ + - ]: 60614 : candidate[Y] = M(1,3);
275 [ + - ][ + - ]: 60614 : candidate[Z] = M(2,3);
276 : : #endif
277 : :
278 : 60616 : return true;
279 : : }
280 : :
281 : :
282 : 60616 : double quadrix_pair_target(const Mat4& Q,
283 : : moab::EntityHandle v1, //Vertex *v1,
284 : : moab::EntityHandle v2, // Vertex *v2,
285 : : Vec3& candidate)
286 : : {
287 : 60616 : int policy = opts.placement_policy;
288 : :
289 : : //
290 : : // This analytic boundary preservation isn't really necessary. The
291 : : // boundary constraint quadrics are quite effective. But, I've left it
292 : : // in anyway.
293 : : //
294 [ + - ]: 60616 : Vec3 vec1 = getVec3FromMBVertex(mb, v1);
295 [ + - ]: 60616 : Vec3 vec2 = getVec3FromMBVertex(mb, v2);
296 [ - + ]: 60616 : if( opts.will_preserve_boundaries )
297 : : {
298 [ # # ]: 0 : int c1 = classifyVertex(mb, v1);
299 [ # # ]: 0 : int c2 = classifyVertex(mb, v2);
300 : :
301 : : // if both are on boundary, put a high penalty cost
302 [ # # ][ # # ]: 0 : if (c1>0 && c2>0)
303 : 0 : return 1.e11;// greater than quality error
304 [ # # ]: 0 : if( c1 > c2 )
305 : : {
306 [ # # ]: 0 : candidate = vec1;
307 [ # # ]: 0 : return quadrix_evaluate_vertex(candidate, Q);
308 : : }
309 [ # # ]: 0 : else if( c2 > c1 )
310 : : {
311 : : //candidate = *v2;
312 [ # # ]: 0 : candidate = vec2;
313 [ # # ]: 0 : return quadrix_evaluate_vertex(candidate, Q);
314 : : }
315 [ # # ][ # # ]: 0 : else if( c1>0 && policy>PLACE_LINE )
316 : 0 : policy = PLACE_LINE;
317 : :
318 [ # # ][ # # ]: 0 : if( policy == PLACE_OPTIMAL ) assert(c1==0 && c2==0);
[ # # ][ # # ]
319 : : }
320 : :
321 [ + - - ]: 60616 : switch( policy )
322 : : {
323 : : case PLACE_OPTIMAL:
324 [ + - ][ + + ]: 60616 : if( quadrix_find_best_fit(Q, candidate) )
325 : 60614 : break;
326 : :
327 : : case PLACE_LINE:
328 [ + - ][ + - ]: 2 : if( quadrix_find_line_fit(Q, vec1, vec2, candidate) )
329 : 2 : break;
330 : :
331 : : default:
332 [ # # ]: 0 : quadrix_find_local_fit(Q, vec1, vec2, candidate);
333 : 0 : break;
334 : : }
335 : :
336 [ + - ]: 60616 : return quadrix_evaluate_vertex(candidate, Q);
337 [ + - ][ + - ]: 312 : }
|
