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