MeshKit  1.0
Go to the documentation of this file.
```00001 // \$Id: quadrics.cxx,v 1.5 1997/10/01 14:07:28 garland Exp \$
00002
00003 //#include "qslim.h"
00004 #include <assert.h>
00006 #include "3D.h"
00007 #include "primitives.h"
00008
00010 //
00011 // Primitive quadric construction and evaluation routines
00012 //
00013
00014 //
00015 // Construct a quadric to evaluate the squared distance of any point
00016 // to the given point v.  Naturally, the iso-surfaces are just spheres
00017 // centered at v.
00018 //
00020 {
00021     Mat4 L(Mat4::identity);
00022
00023     L(0,3) = -v[0];
00024     L(1,3) = -v[1];
00025     L(2,3) = -v[2];
00026     L(3,3) = v*v;
00027
00028     L(3,0) = L(0,3);
00029     L(3,1) = L(1,3);
00030     L(3,2) = L(2,3);
00031
00032     return L;
00033 }
00035 {
00036     Mat4 L(Mat4::identity);
00037     double v[3];
00038     mb->get_coords(&vert, 1 , v);
00039     L(0,3) = -v[0];
00040     L(1,3) = -v[1];
00041     L(2,3) = -v[2];
00042     L(3,3) = v[0]*v[0] + v[1]*v[1] + v[2]*v[2] ; //  v*v;
00043
00044     L(3,0) = L(0,3);
00045     L(3,1) = L(1,3);
00046     L(3,2) = L(2,3);
00047
00048     return L;
00049 }
00050 //
00051 // Construct a quadric to evaluate the squared distance of any point
00052 // to the given plane [ax+by+cz+d = 0].  This is the "fundamental error
00053 // quadric" discussed in the paper.
00054 //
00055 Mat4 quadrix_plane_constraint(double a, double b, double c, double d)
00056 {
00057     Mat4 K(Mat4::zero);
00058
00059     K(0,0) = a*a;   K(0,1) = a*b;   K(0,2) = a*c;  K(0,3) = a*d;
00060     K(1,0) =K(0,1); K(1,1) = b*b;   K(1,2) = b*c;  K(1,3) = b*d;
00061     K(2,0) =K(0,2); K(2,1) =K(1,2); K(2,2) = c*c;  K(2,3) = c*d;
00062     K(3,0) =K(0,3); K(3,1) =K(1,3); K(3,2) =K(2,3);K(3,3) = d*d;
00063
00064     return K;
00065 }
00066
00067 //
00068 // Define some other convenient ways for constructing these plane quadrics.
00069 //
00070 Mat4 quadrix_plane_constraint(const Vec3& n, double d)
00071 {
00072     return quadrix_plane_constraint(n[X], n[Y], n[Z], d);
00073 }
00074
00077 {
00078     // const Plane& p = T.plane();
00079     Plane p=trianglePlane(mb, triangle);
00080     double a,b,c,d;
00081     p.coeffs(&a, &b, &c, &d);
00082
00083     return quadrix_plane_constraint(a, b, c, d);
00084 }
00085
00086 Mat4 quadrix_plane_constraint(const Vec3& v1, const Vec3& v2, const Vec3& v3)
00087 {
00088     Plane P(v1,v2,v3);
00089     double a,b,c,d;
00090     P.coeffs(&a, &b, &c, &d);
00091
00092     return quadrix_plane_constraint(a, b, c, d);
00093 }
00094
00095 double quadrix_evaluate_vertex(const Vec3& v, const Mat4& K)
00096 {
00097     double x=v[X], y=v[Y], z=v[Z];
00098
00099 #ifndef VECTOR_COST_EVALUATION
00100     //
00101     // This is the fast way of computing (v^T Q v).
00102     //
00103     return x*x*K(0,0) + 2*x*y*K(0,1) + 2*x*z*K(0,2) + 2*x*K(0,3)
00104                       + y*y*K(1,1)   + 2*y*z*K(1,2) + 2*y*K(1,3)
00105                                      + z*z*K(2,2)   + 2*z*K(2,3)
00106                                                     + K(3,3);
00107 #else
00108     //
00109     // The equivalent thing using matrix/vector operations.
00110     // It's a lot clearer, but it's also slower.
00111     //
00112     Vec4 v2(x,y,z,1);
00113     return v2*(K*v2);
00114 #endif
00115 }
00116
00117
00118
00120 //
00121 // Routines for computing discontinuity constraints
00122 //
00123
00124 //bool is_border(Edge *e )
00125 bool is_border(moab::EntityHandle eh )
00126 {
00127    //  use nb of triangles connected to check if number of tri is 1
00130                                       moab::Interface::UNION);
00132          return true;
00133     return false;
00134 }
00135
00136 bool check_for_discontinuity(moab::EntityHandle eh) //Edge *e)
00137 {
00138     return is_border(eh);
00139 }
00140
00142     /*Edge *edge*/, const Vec3& n)
00143 {
00144     //Vec3& org = *edge->org();
00145     //Vec3& dest = *edge->dest();
00146     // to do: Vec3 origin of edge
00147         // is the orientation important for an edge?
00148         const moab::EntityHandle * conn;
00149         int num_nodes;
00150         mb->get_connectivity(mbe, conn, num_nodes);
00151     Vec3 dest = getVec3FromMBVertex(mb, conn[1]); //
00152     Vec3 org = getVec3FromMBVertex(mb, conn[0]);
00153     Vec3 e = dest - org;
00154
00155     Vec3 n2 = e ^ n;
00156     unitize(n2);
00157
00158     double d = -n2 * org;
00160 }
00161
00162
00163 Mat4 quadrix_discontinuity_constraint(moab::EntityHandle mbedge)// Edge *edge)
00164 {
00165     Mat4 D(Mat4::zero);
00166
00167 /*
00168     face_buffer& faces = edge->faceUses();
00169
00170     for(int i=0; i<faces.length(); i++)
00172 */
00173     // to do : to get the connected faces to the edge
00174
00177                                                                   moab::Interface::UNION);
00178         assert(rval == moab::MB_SUCCESS);
00179
00180         for (unsigned int i=0; i<adjFaces.size(); i++)
00181         {
00183                 Plane p = trianglePlane( mb, F);
00185         }
00186     return D;
00187 }
00188
00189
00190
00192 //
00193 // Routines for computing contraction target
00194 //
00195
00197                             const Vec3& v1, const Vec3& v2,
00198                             Vec3& candidate)
00199 {
00200
00201     Vec3 v3 = (v1 + v2) / 2;
00202
00203     bool try_midpoint = opts.placement_policy > PLACE_ENDPOINTS;
00204
00205     double c1 = quadrix_evaluate_vertex(v1, K);
00206     double c2 = quadrix_evaluate_vertex(v2, K);
00207     double c3;
00208     if( try_midpoint ) c3 = quadrix_evaluate_vertex(v3, K);
00209
00210     if( c1<c2 )
00211     {
00212         if( try_midpoint && c3<c1 )
00213             candidate=v3;
00214         else
00215             candidate=v1;
00216     }
00217     else
00218     {
00219         if( try_midpoint && c3<c2 )
00220             candidate=v3;
00221         else
00222             candidate=v2;
00223     }
00224
00225     return true;
00226 }
00227
00229                            const Vec3& v1, const Vec3& v2,
00230                            Vec3& candidate)
00231 {
00232     Vec3 d = v1-v2;
00233
00234     Vec3 Qv2 = Q*v2;
00235     Vec3 Qd  = Q*d;
00236
00237     double denom = 2*d*Qd;
00238
00239     if( denom == 0.0 )
00240         return false;
00241
00242     double a = (d*Qv2 + v2*Qd) / denom;
00243
00244     if( a<0.0 ) a=0.0;
00245     if( a>1.0 ) a=1.0;
00246
00247
00248     candidate = a*d + v2;
00249     return true;
00250 }
00251
00252 bool quadrix_find_best_fit(const Mat4& Q, Vec3& candidate)
00253 {
00254     Mat4 K = Q;
00255     K(3,0) = K(3,1) = K(3,2) = 0.0;  K(3,3) = 1;
00256
00257
00258     Mat4 M;
00259     double det = K.inverse(M);
00260     if( FEQ(det, 0.0, 1e-12) )
00261         return false;
00262
00263
00264 #ifdef SAFETY
00265     //
00266     // The homogeneous division SHOULDN'T be necessary.
00267     // But, when we're being SAFE, we do it anyway just in case.
00268     //
00269     candidate[X] = M(0,3)/M(3,3);
00270     candidate[Y] = M(1,3)/M(3,3);
00271     candidate[Z] = M(2,3)/M(3,3);
00272 #else
00273     candidate[X] = M(0,3);
00274     candidate[Y] = M(1,3);
00275     candidate[Z] = M(2,3);
00276 #endif
00277
00278     return true;
00279 }
00280
00281
00283                          moab::EntityHandle v1, //Vertex *v1,
00284                          moab::EntityHandle v2, // Vertex *v2,
00285                          Vec3& candidate)
00286 {
00287     int policy = opts.placement_policy;
00288
00289     //
00290     // This analytic boundary preservation isn't really necessary.  The
00291     // boundary constraint quadrics are quite effective.  But, I've left it
00292     // in anyway.
00293     //
00294     Vec3 vec1 = getVec3FromMBVertex(mb, v1);
00295     Vec3 vec2 = getVec3FromMBVertex(mb, v2);
00296     if( opts.will_preserve_boundaries )
00297     {
00298         int c1 = classifyVertex(mb, v1);
00299         int c2 = classifyVertex(mb, v2);
00300
00301         // if both are on boundary, put a high penalty cost
00302         if (c1>0 && c2>0)
00303                 return 1.e11;// greater than quality error
00304         if( c1 > c2 )
00305         {
00306                 candidate = vec1;
00308         }
00309         else if( c2 > c1 )
00310                 {
00311                         //candidate = *v2;
00312                         candidate = vec2;
00314                 }
00315                 else if( c1>0 && policy>PLACE_LINE )
00316                         policy = PLACE_LINE;
00317
00318         if( policy == PLACE_OPTIMAL ) assert(c1==0 && c2==0);
00319     }
00320
00321     switch( policy )
00322     {
00323     case PLACE_OPTIMAL:
00325                 break;
00326
00327     case PLACE_LINE:
00328         if( quadrix_find_line_fit(Q, vec1, vec2, candidate) )
00329                 break;
00330
00331     default: