Blender V2.61 - r43446

mvmcoords.cpp

Go to the documentation of this file.
00001 
00004 /******************************************************************************
00005  *
00006 // El'Beem - the visual lattice boltzmann freesurface simulator
00007 // All code distributed as part of El'Beem is covered by the version 2 of the 
00008 // GNU General Public License. See the file COPYING for details.  
00009 //
00010 // Copyright 2008 Nils Thuerey , Richard Keiser, Mark Pauly, Ulrich Ruede
00011 //
00012  *
00013  * Mean Value Mesh Coords class
00014  *
00015  *****************************************************************************/
00016 
00017 #include "mvmcoords.h"
00018 #include <algorithm>
00019 using std::vector;
00020 
00021 void MeanValueMeshCoords::clear() 
00022 {
00023     mVertices.resize(0);
00024     mNumVerts = 0;
00025 }
00026 
00027 void MeanValueMeshCoords::calculateMVMCs(vector<ntlVec3Gfx> &reference_vertices, vector<ntlTriangle> &tris, 
00028         vector<ntlVec3Gfx> &points, gfxReal numweights)
00029 {
00030     clear();
00031     mvmTransferPoint tds;
00032     int mem = 0;
00033     int i = 0;
00034 
00035     mNumVerts = (int)reference_vertices.size();
00036 
00037     for (vector<ntlVec3Gfx>::iterator iter = points.begin(); iter != points.end(); ++iter, ++i) {
00038         /*
00039         if(i%(points.size()/10)==1) debMsgStd("MeanValueMeshCoords::calculateMVMCs",DM_MSG,"Computing weights, points: "<<i<<"/"<<points.size(),5 );
00040         */
00041         tds.lastpos = *iter;
00042         tds.weights.resize(0);  // clear
00043         computeWeights(reference_vertices, tris, tds, numweights);
00044         mem += (int)tds.weights.size();
00045         mVertices.push_back(tds);
00046     }    
00047     int mbmem = mem * sizeof(mvmFloat) / (1024*1024);
00048     debMsgStd("MeanValueMeshCoords::calculateMVMCs",DM_MSG,"vertices:"<<mNumVerts<<" points:"<<points.size()<<" weights:"<<mem<<", wmem:"<<mbmem<<"MB ",7 );
00049 }
00050 
00051 // from: mean value coordinates for closed triangular meshes
00052 // attention: fails if a point is exactly (or very close) to a vertex
00053 void MeanValueMeshCoords::computeWeights(vector<ntlVec3Gfx> &reference_vertices, vector<ntlTriangle>& tris,
00054         mvmTransferPoint& tds, gfxReal numweights)
00055 {
00056     const bool mvmFullDebug=false;
00057     //const ntlVec3Gfx cEPS = 1.0e-6;
00058     const mvmFloat cEPS = 1.0e-14;
00059 
00060     //mvmFloat d[3], s[3], phi[3],c[3];
00061     ntlVec3d u[3],c,d,s,phi;
00062     int indices[3];
00063 
00064     for (int i = 0; i < (int)reference_vertices.size(); ++i) {
00065         tds.weights.push_back(mvmIndexWeight(i, 0.0));
00066     }
00067 
00068     // for each triangle
00069     //for (vector<ntlTriangle>::iterator iter = tris.begin(); iter != tris.end();) {
00070     for(int t=0; t<(int)tris.size(); t++) {
00071 
00072         for (int i = 0; i < 3; ++i) { //, ++iter) {
00073             indices[i] = tris[t].getPoints()[i];
00074             u[i] = vec2D(reference_vertices[ indices[i] ]-tds.lastpos);
00075             d[i] = normalize(u[i]); //.normalize();        
00076             //assert(d[i] != 0.);
00077             if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","t"<<t<<" i"<<indices[i] //<<" lp"<<tds.lastpos
00078                     <<" v"<<reference_vertices[indices[i]]<<" u"<<u[i]<<" ");
00079             // on vertex!
00080             //? if(d[i]<=0.) continue;
00081         }
00082         //for (int i = 0; i < 3; ++i) { errMsg("III"," "<<i     <<" i"<<indices[i]<<reference_vertices[ indices[i] ] ); }
00083 
00084         // arcsin is not needed, see paper
00085         phi[0] = 2.*asin( (mvmFloat)(0.5* norm(u[1]-u[2]) )  );
00086         phi[1] = 2.*asin( (mvmFloat)(0.5* norm(u[0]-u[2]) )  );
00087         phi[2] = 2.*asin( (mvmFloat)(0.5* norm(u[0]-u[1]) )  );
00088         mvmFloat h = (phi[0] + phi[1] + phi[2])*0.5;
00089         if (M_PI-h < cEPS) {
00090             if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","point on triangle");
00091             tds.weights.resize(0);
00092             tds.weights.push_back( mvmIndexWeight(indices[0], sin(phi[0])*d[1]*d[2]));
00093             tds.weights.push_back( mvmIndexWeight(indices[1], sin(phi[1])*d[0]*d[2]));
00094             tds.weights.push_back( mvmIndexWeight(indices[2], sin(phi[2])*d[1]*d[0]));
00095             break;
00096         }
00097         mvmFloat sinh = 2.*sin(h);
00098         c[0] = (sinh*sin(h-phi[0]))/(sin(phi[1])*sin(phi[2]))-1.;
00099         c[1] = (sinh*sin(h-phi[1]))/(sin(phi[0])*sin(phi[2]))-1.;    
00100         c[2] = (sinh*sin(h-phi[2]))/(sin(phi[0])*sin(phi[1]))-1.;   
00101         if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","c="<<c<<" phi="<<phi<<" d="<<d);
00102         //if (c[0] > 1. || c[0] < 0. || c[1] > 1. || c[1] < 0. || c[2] > 1. || c[2] < 0.) continue;
00103 
00104         s[0] = sqrt((float)(1.-c[0]*c[0]));
00105         s[1] = sqrt((float)(1.-c[1]*c[1]));
00106         s[2] = sqrt((float)(1.-c[2]*c[2]));
00107 
00108         if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","s");
00109         if (s[0] <= cEPS || s[1] <= cEPS || s[2] <= cEPS) {
00110             //MSG("position lies outside the triangle on the same plane -> ignore it");
00111             continue;
00112         }
00113         const mvmFloat u0x = u[0][0];
00114         const mvmFloat u0y = u[0][1];
00115         const mvmFloat u0z = u[0][2];
00116         const mvmFloat u1x = u[1][0];
00117         const mvmFloat u1y = u[1][1];
00118         const mvmFloat u1z = u[1][2];
00119         const mvmFloat u2x = u[2][0];
00120         const mvmFloat u2y = u[2][1];
00121         const mvmFloat u2z = u[2][2];
00122         mvmFloat det = u0x*u1y*u2z - u0x*u1z*u2y + u0y*u1z*u2x - u0y*u1x*u2z + u0z*u1x*u2y - u0z*u1y*u2x;
00123         //assert(det != 0.);
00124         if (det < 0.) {
00125             s[0] = -s[0];
00126             s[1] = -s[1];
00127             s[2] = -s[2];
00128         }
00129 
00130         tds.weights[indices[0]].weight += (phi[0]-c[1]*phi[2]-c[2]*phi[1])/(d[0]*sin(phi[1])*s[2]);
00131         tds.weights[indices[1]].weight += (phi[1]-c[2]*phi[0]-c[0]*phi[2])/(d[1]*sin(phi[2])*s[0]);
00132         tds.weights[indices[2]].weight += (phi[2]-c[0]*phi[1]-c[1]*phi[0])/(d[2]*sin(phi[0])*s[1]);
00133         if(mvmFullDebug) { errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[0]<<" o"<<tds.weights[indices[0]].weight);
00134         errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[1]<<" o"<<tds.weights[indices[1]].weight);
00135         errMsg("MeanValueMeshCoords::computeWeights","i"<<indices[2]<<" o"<<tds.weights[indices[2]].weight);
00136         errMsg("MeanValueMeshCoords::computeWeights","\n\n\n"); }
00137     }
00138 
00139     //sort weights
00140     if((numweights>0.)&& (numweights<1.) ) {
00141     //if( ((int)tds.weights.size() > maxNumWeights) && (maxNumWeights > 0) ) {
00142       int maxNumWeights = (int)(tds.weights.size()*numweights);
00143         if(maxNumWeights<=0) maxNumWeights = 1;
00144         std::sort(tds.weights.begin(), tds.weights.end(), std::greater<mvmIndexWeight>());
00145         // only use maxNumWeights-th largest weights
00146         tds.weights.resize(maxNumWeights);
00147     }
00148 
00149     // normalize weights
00150     mvmFloat totalWeight = 0.;
00151     for (vector<mvmIndexWeight>::const_iterator witer = tds.weights.begin();
00152             witer != tds.weights.end(); ++witer) {
00153         totalWeight += witer->weight;
00154     }
00155     mvmFloat invTotalWeight;
00156     if (totalWeight == 0.) {
00157         if(mvmFullDebug) errMsg("MeanValueMeshCoords::computeWeights","totalWeight == 0");
00158         invTotalWeight = 0.0;
00159     } else {
00160         invTotalWeight = 1.0/totalWeight;
00161     }
00162 
00163     for (vector<mvmIndexWeight>::iterator viter = tds.weights.begin();
00164             viter != tds.weights.end(); ++viter) {
00165         viter->weight *= invTotalWeight;  
00166         //assert(finite(viter->weight) != 0);
00167         if(!finite(viter->weight)) viter->weight=0.;
00168     }
00169 }
00170 
00171 void MeanValueMeshCoords::transfer(vector<ntlVec3Gfx> &vertices, vector<ntlVec3Gfx>& displacements) 
00172 {
00173     displacements.resize(0);
00174 
00175     //debMsgStd("MeanValueMeshCoords::transfer",DM_MSG,"vertices:"<<mNumVerts<<" curr_verts:"<<vertices.size()<<" ",7 );
00176     if((int)vertices.size() != mNumVerts) {
00177         errMsg("MeanValueMeshCoords::transfer","Different no of verts: "<<vertices.size()<<" vs "<<mNumVerts);
00178         return;
00179     }
00180 
00181     for (vector<mvmTransferPoint>::iterator titer = mVertices.begin(); titer != mVertices.end(); ++titer) {
00182         mvmTransferPoint &tds = *titer;
00183         ntlVec3Gfx newpos(0.0);
00184 
00185         for (vector<mvmIndexWeight>::iterator witer = tds.weights.begin();
00186                 witer != tds.weights.end(); ++witer) {
00187             newpos += vertices[witer->index] * witer->weight;
00188             //errMsg("transfer","np"<<newpos<<" v"<<vertices[witer->index]<<" w"<< witer->weight);
00189         }
00190 
00191         displacements.push_back(newpos);
00192         //displacements.push_back(newpos - tds.lastpos);
00193         //tds.lastpos = newpos;
00194     }
00195 }
00196