MeshKit  1.0
SurfProHarmonicMap.cpp
Go to the documentation of this file.
00001 #include "SurfProHarmonicMap.hpp"
00002 #include <iostream>
00003 #include <math.h>
00004 #include <map>
00005 
00006 namespace MeshKit {
00007 
00008 SurfProHarmonicMap::SurfProHarmonicMap(MKCore* core, iBase_EntityHandle s, iBase_EntityHandle t, iBase_EntityHandle v)
00009 {
00010   mk_core = core;
00011   source.gFaceHandle = s;
00012   target.gFaceHandle = t;
00013   volume = v;
00014   
00015   irel_pair = mk_core->irel_pair();
00016   irel_instance = mk_core->irel_instance();
00017   igeom_instance = mk_core->igeom_instance();
00018   imesh_instance = mk_core->imesh_instance();
00019 
00020   g_err = igeom_instance->getTagHandle("GLOBAL_ID", global_geom_tag);
00021   IBERRCHK(g_err, "Trouble get the geom_id_tag for 'GLOBAL_ID'.");
00022   m_err = imesh_instance->getTagHandle("GLOBAL_ID", global_mesh_tag);
00023   IBERRCHK(m_err, "Trouble get the mesh_id_tag for 'GLOBAL_ID'.");
00024   g_err = igeom_instance->createTag("HARMONIC_SURF_PROJECTION", 1, iBase_INTEGER, harmonic_surf_pro);
00025   IBERRCHK(g_err, "Trouble create the tag handle.");
00026   m_err = imesh_instance->createTag("HARMONIC_FACET_MESH", 1, iBase_INTEGER, facet_mesh_tag);
00027   IBERRCHK(m_err, "Trouble create the tag handle.");
00028   preprocessing();
00029   getFacets();
00030 }
00031 
00032 void SurfProHarmonicMap::projection()
00033 {
00034         //compute the normal vector
00035         vector<Vector3D> normal_src, normal_tgt;
00036         normal_src.resize(src_facet_tri.size());
00037         normal_tgt.resize(tgt_facet_tri.size());
00038         for (unsigned int i = 0; i < src_facet_tri.size(); i++){
00039                 computeNormalTri(src_facet_tri[i], normal_src[i], source);
00040                 //std::cout << "source index = " << i << "\tnrml = {" << normal_src[i][0] << "," << normal_src[i][1] << "," << normal_src[i][2] << "}\n";
00041         }
00042         for (unsigned int i = 0; i < tgt_facet_tri.size(); i++){
00043                 computeNormalTri(tgt_facet_tri[i], normal_tgt[i], target);
00044                 //std::cout << "target index = " << i << "\tnrml = {" << normal_tgt[i][0] << "," << normal_tgt[i][1] << "," << normal_tgt[i][2] << "}\n";
00045         }
00046 
00047         //loop over the interior nodes of quad mesh on the source surface and compute barycentric coordinates in the physical source surface
00048         for (vector<Vertex>::iterator it = quad_mesh_src.begin(); it != quad_mesh_src.end(); it++){
00049                 if (it->onBoundary)
00050                         continue;
00051         //std::cout << "index=" << it->index  << "\tsrc quad node xyz = {" << it->xyz[0] << "," << it->xyz[1] << "," << it->xyz[2] << "}\t";
00052         Vector3D bary_coord_src, bary_coord_tgt;
00053                 int tri_index = findFacetTri(src_facet_tri, normal_src, it->xyz, bary_coord_src);
00054                 //compute the positions of all the quad mesh nodes on the unit disk
00055                 it->uv[0] = 0.0;
00056                 it->uv[1] = 0.0;
00057                 for (int i = 0; i < 3; i++){
00058                         it->uv[0] += bary_coord_src[i]*src_facet_tri[tri_index].connect[i]->uv[0];
00059                         it->uv[1] += bary_coord_src[i]*src_facet_tri[tri_index].connect[i]->uv[1];
00060                 }
00061         //std::cout << "\t\tsrc uv={" << it->uv[0] << "," << it->uv[1] << "}  facet index=" << tri_index << "\tv1={" << src_facet_tri[tri_index].connect[0]->xyz[0] << "," << src_facet_tri[tri_index].connect[0]->xyz[1] << "," << src_facet_tri[tri_index].connect[0]->xyz[2] << "} v2={" << src_facet_tri[tri_index].connect[1]->xyz[0] << "," << src_facet_tri[tri_index].connect[1]->xyz[1] << "," << src_facet_tri[tri_index].connect[1]->xyz[2] << "} v3={" << src_facet_tri[tri_index].connect[2]->xyz[0] << "," << src_facet_tri[tri_index].connect[2]->xyz[1] << "," << src_facet_tri[tri_index].connect[2]->xyz[2] << "}\n"; 
00062                 //compute the barycentric coordinates of those quad mesh nodes on the unit disk of target surface
00063                 tri_index = findFacetTri(tgt_facet_tri, it->uv, bary_coord_tgt);
00064                 Vector3D xyz(0.0);      
00065                 for (int i = 0; i < 3; i++)
00066                         for (int j = 0; j < 3; j++)
00067                                 xyz[j] += bary_coord_tgt[i]*tgt_facet_tri[tri_index].connect[i]->xyz[j];
00068                 
00069 
00070                 g_err = igeom_instance->getEntClosestPtTrimmed(target.gFaceHandle, xyz[0], xyz[1], xyz[2], quad_mesh_tgt[it->index].xyz[0], quad_mesh_tgt[it->index].xyz[1], quad_mesh_tgt[it->index].xyz[2]);
00071                 IBERRCHK(g_err, "Trouble get the trimmed closed point positions on the target surface!");
00072 
00073          //std::cout << "\t\ttgt uvw={" << bary_coord_tgt[0] << "," << bary_coord_tgt[1] << ","  << bary_coord_tgt[2] << "}  facet index=" << tri_index << "\tv1={" << tgt_facet_tri[tri_index].connect[0]->xyz[0] << "," << tgt_facet_tri[tri_index].connect[0]->xyz[1] << "," << tgt_facet_tri[tri_index].connect[0]->xyz[2] << "} v2={" << tgt_facet_tri[tri_index].connect[1]->xyz[0] << "," << tgt_facet_tri[tri_index].connect[1]->xyz[1] << "," << tgt_facet_tri[tri_index].connect[1]->xyz[2] << "} v3={" << tgt_facet_tri[tri_index].connect[2]->xyz[0] << "," << tgt_facet_tri[tri_index].connect[2]->xyz[1] << "," << tgt_facet_tri[tri_index].connect[2]->xyz[2] << "}\n";
00074          //std::cout << "\t\tprojected pts xyz = {" << quad_mesh_tgt[it->index].xyz[0] << "," << quad_mesh_tgt[it->index].xyz[1] << "," << quad_mesh_tgt[it->index].xyz[2] << "}\ttmp = {" << xyz[0] << "," << xyz[1] << "," << xyz[2] << "}\n";
00075         }
00076         /*
00077         vector<iBase_EntityHandle> nodes(quad_mesh_tgt.size());
00078         for (unsigned int i = 0; i < quad_mesh_src.size(); i++){
00079                 if (quad_mesh_src[i].onBoundary)
00080                         continue;
00081                 //create vtx
00082                 m_err = imesh_instance->createVtx(quad_mesh_tgt[i].xyz[0], quad_mesh_tgt[i].xyz[1], quad_mesh_tgt[i].xyz[2], nodes[i]);
00083                 IBERRCHK(m_err, "Trouble create vtx ent!");
00084         }
00085 
00086         
00087         for (unsigned int i = 0; i < facelist.size(); i++){
00088                 bool is_boundary = true;
00089                 for (int j = 0; j < 4; j++)
00090                         is_boundary = is_boundary && !facelist[i].connect[j]->onBoundary;
00091                 if (is_boundary){
00092                         iBase_EntityHandle quads;
00093                         vector<iBase_EntityHandle> tmp;
00094                         tmp.push_back(nodes[facelist[i].connect[0]->index]);
00095                         tmp.push_back(nodes[facelist[i].connect[1]->index]);
00096                         tmp.push_back(nodes[facelist[i].connect[2]->index]);
00097                         tmp.push_back(nodes[facelist[i].connect[3]->index]);
00098                         m_err = imesh_instance->createEnt(iMesh_QUADRILATERAL, &tmp[0], 4, quads);
00099                         IBERRCHK(m_err, "Trouble create an face entity!");                                              
00100                 }
00101         }
00102         */
00103         cleanup();
00104 }
00105 
00106 bool SurfProHarmonicMap::ComputeBarycentric(Vector3D a, Vector3D b, Vector3D c, Vector3D xyz, Vector3D &uvw)
00107 {
00108         //Compute vectors;
00109         Vector3D v0 = b - a;
00110         Vector3D v1 = c - a;
00111         Vector3D v2 = xyz - a;
00112         //compute dot products
00113         double dot00 = v0[0]*v0[0] + v0[1]*v0[1] + v0[2]*v0[2];
00114         double dot01 = v0[0]*v1[0] + v0[1]*v1[1] + v0[2]*v1[2];
00115         double dot02 = v0[0]*v2[0] + v0[1]*v2[1] + v0[2]*v2[2];
00116         double dot11 = v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2];
00117         double dot12 = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
00118         //compute barycentric coordinates
00119         double invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01);
00120         uvw[1] = (dot11 * dot02 - dot01 * dot12) * invDenom;
00121         uvw[2] = (dot00 * dot12 - dot01 * dot02) * invDenom;
00122         uvw[0] = 1.0 - uvw[1] - uvw[2];
00123 
00124         //check if the point is in the triangle
00125         return (uvw[0]>=-eps)&&(uvw[1]>=-eps)&&(uvw[0]+uvw[1]<=1.0+eps);
00126 }
00127 
00128 void SurfProHarmonicMap::test()
00129 {
00130         vector<iBase_EntityHandle> test_nodes(src_facet_v.size()), test_edges(src_facet_e.size()), test_tri(src_facet_tri.size());
00131         for (vector<Vertex>::iterator it = src_facet_v.begin(); it != src_facet_v.end(); it++){
00132                 m_err = imesh_instance->createVtx(it->uv[0], it->uv[1], 200.0, test_nodes[it->index]);
00133                 IBERRCHK(m_err, "Trouble create a test vertex!");
00134         }
00135         for (vector<Face>::iterator it = src_facet_tri.begin(); it != src_facet_tri.end(); it++){
00136                 vector<iBase_EntityHandle> tmp;
00137                 tmp.push_back(test_nodes[it->connect[0]->index]);
00138                 tmp.push_back(test_nodes[it->connect[1]->index]);
00139                 tmp.push_back(test_nodes[it->connect[2]->index]);
00140                 m_err = imesh_instance->createEnt(iMesh_TRIANGLE, &tmp[0], tmp.size(), test_tri[it->index]);
00141                 IBERRCHK(m_err, "Trouble create a triangle mesh entity!");
00142         }
00143         test_nodes.clear();
00144         test_edges.clear();
00145         test_tri.clear();
00146         test_nodes.resize(tgt_facet_v.size());
00147         test_tri.resize(tgt_facet_tri.size());
00148 
00149         for (vector<Vertex>::iterator it = tgt_facet_v.begin(); it != tgt_facet_v.end(); it++){
00150                 m_err = imesh_instance->createVtx(it->uv[0], it->uv[1], -200.0, test_nodes[it->index]);
00151                 IBERRCHK(m_err, "Trouble create a test vertex!");
00152                 //std::cout << "unit disk vertex index = " << it->index << "\tuv = {" << it->uv[0] << ", " << it->uv[1] << "}\tis_boundary = " << it->onBoundary << std::endl;
00153         }
00154         for (vector<Face>::iterator it = tgt_facet_tri.begin(); it != tgt_facet_tri.end(); it++){
00155                 vector<iBase_EntityHandle> tmp;
00156                 tmp.push_back(test_nodes[it->connect[0]->index]);
00157                 tmp.push_back(test_nodes[it->connect[1]->index]);
00158                 tmp.push_back(test_nodes[it->connect[2]->index]);
00159                 m_err = imesh_instance->createEnt(iMesh_TRIANGLE, &tmp[0], tmp.size(), test_tri[it->index]);
00160                 IBERRCHK(m_err, "Trouble create a triangle mesh entity!");
00161         }       
00162 }
00163 
00164 void SurfProHarmonicMap::cleanup(){
00165         vector<iBase_EntityHandle> ents_to_del;
00166         for (vector<Edge>::iterator it = src_facet_e.begin(); it != src_facet_e.end(); it++)
00167                 ents_to_del.push_back(it->gEdgeHandle);
00168         for (vector<Edge>::iterator it = tgt_facet_e.begin(); it != tgt_facet_e.end(); it++)
00169                 ents_to_del.push_back(it->gEdgeHandle);
00170         for (vector<Face>::iterator it = src_facet_tri.begin(); it != src_facet_tri.end(); it++)
00171                 ents_to_del.push_back(it->gFaceHandle);
00172         for (vector<Face>::iterator it = tgt_facet_tri.begin(); it != tgt_facet_tri.end(); it++)
00173                 ents_to_del.push_back(it->gFaceHandle);
00174         for (vector<Vertex>::iterator it = src_facet_v.begin(); it != src_facet_v.end(); it++)
00175                 ents_to_del.push_back(it->gVertexHandle);
00176         for (vector<Vertex>::iterator it = tgt_facet_v.begin(); it != tgt_facet_v.end(); it++)
00177                 ents_to_del.push_back(it->gVertexHandle);
00178         m_err = imesh_instance->rmvArrTag(&ents_to_del[0], ents_to_del.size(), facet_mesh_tag);
00179         IBERRCHK(m_err, "Trouble remove the tag data!");
00180         m_err = imesh_instance->deleteEntArr(&ents_to_del[0], ents_to_del.size());
00181         IBERRCHK(m_err, "Trouble delete entities!");
00182         g_err = igeom_instance->destroyTag(harmonic_surf_pro, true);
00183         IBERRCHK(g_err, "Trouble delete a tag!");
00184         m_err = imesh_instance->destroyTag(facet_mesh_tag, true);
00185         IBERRCHK(g_err, "Trouble delete a tag!");
00186                 
00187 }
00188 
00189 int SurfProHarmonicMap::findFacetTri(vector<Face> &facet_tri, Vector2D uv, Vector3D &uvw){
00190         for (unsigned int i = 0; i < facet_tri.size(); i++){
00191                 if (ComputeBarycentric(facet_tri[i].connect[0]->uv, facet_tri[i].connect[1]->uv, facet_tri[i].connect[2]->uv, uv, uvw))
00192                         return i;
00193         }
00194 
00195         return -1;
00196 }
00197 
00198 bool SurfProHarmonicMap::ComputeBarycentric(Vector2D a, Vector2D b, Vector2D c, Vector2D xy, Vector3D &uvw)
00199 {
00200         //Compute vectors;
00201         Vector2D v0 = b - a;
00202         Vector2D v1 = c - a;
00203         Vector2D v2 = xy - a;
00204         //compute dot products
00205         double dot00 = v0[0]*v0[0] + v0[1]*v0[1];
00206         double dot01 = v0[0]*v1[0] + v0[1]*v1[1];
00207         double dot02 = v0[0]*v2[0] + v0[1]*v2[1];
00208         double dot11 = v1[0]*v1[0] + v1[1]*v1[1];
00209         double dot12 = v1[0]*v2[0] + v1[1]*v2[1];
00210         //compute barycentric coordinates
00211         double invDenom = 1.0 / (dot00 * dot11 - dot01 * dot01);
00212         uvw[1] = (dot11 * dot02 - dot01 * dot12) * invDenom;
00213         uvw[2] = (dot00 * dot12 - dot01 * dot02) * invDenom;
00214         uvw[0] = 1.0 - uvw[1] - uvw[2];
00215 
00216         //check if the point is in the triangle
00217         return (uvw[0]>=-eps)&&(uvw[1]>=-eps)&&(uvw[0]+uvw[1]<=1+eps);
00218 }
00219 //compute the normal of a triangle
00220 void SurfProHarmonicMap::computeNormalTri(Face &tri, Vector3D& nrml, Face surf)
00221 {
00222         Vector3D a = tri.connect[0]->xyz, b = tri.connect[1]->xyz, c = tri.connect[2]->xyz;
00223         Vector3D va = b - a, vb = c - b;
00224         nrml[0] = va[1]*vb[2] - va[2]*vb[1];
00225         nrml[1] = va[2]*vb[0] - va[0]*vb[2];
00226         nrml[2] = va[0]*vb[1] - va[1]*vb[0];
00227         double len = sqrt(nrml[0]*nrml[0] + nrml[1]*nrml[1] + nrml[2]*nrml[2]);
00228         for (int i = 0; i < 3; i++)
00229                 nrml[i] /= len;
00230 
00231         Vector3D surf_nrml(0.0);
00232         g_err = igeom_instance->getEntNrmlXYZ(surf.gFaceHandle, tri.connect[0]->xyz[0], tri.connect[0]->xyz[1], tri.connect[0]->xyz[2], surf_nrml[0], surf_nrml[1], surf_nrml[2]);
00233         IBERRCHK(g_err, "Trouble get the surface normal at a given position!");
00234         double dotproduct = surf_nrml[0]*nrml[0] + surf_nrml[1]*nrml[1] + surf_nrml[2]*nrml[2];
00235         if (dotproduct < 0){
00236                 for (int j = 0; j < 3; j++)
00237                         nrml[j] = -1.0*nrml[j];
00238                 Vertex *tmp_vtx = tri.connect[1];
00239                 tri.connect[1] = tri.connect[2];
00240                 tri.connect[2] = tmp_vtx;
00241         }
00242 }
00243 //find a specific triangle where the point is located
00244 int SurfProHarmonicMap::findFacetTri(vector<Face> &facet_tri, vector<Vector3D> nrml, Vector3D xyz, Vector3D &uvw)
00245 {
00246         for (unsigned int i = 0; i < facet_tri.size(); i++){
00247                 Vector3D prj_pts;
00248                 prjPtsToTri(facet_tri[i], xyz, nrml[i], prj_pts);
00249                 
00250                 //bool test_bool = ComputeBarycentric(facet_tri[i].connect[0]->xyz, facet_tri[i].connect[1]->xyz, facet_tri[i].connect[2]->xyz, prj_pts, uvw);
00251                 //std::cout << "\t\t\tindex = " << i << "\tbool = "<< test_bool << "\tbary uvw = {" << uvw[0] << "," << uvw[1] << "," << uvw[2] << "}\n";
00252 
00253                 if (ComputeBarycentric(facet_tri[i].connect[0]->xyz, facet_tri[i].connect[1]->xyz, facet_tri[i].connect[2]->xyz, prj_pts, uvw)){
00254                         //std::cout << "projected pts = {" << prj_pts[0] << "," << prj_pts[1] << "," << prj_pts[2] << "}\n";
00255                         //double x = uvw[0]*facet_tri[i].connect[0]->xyz[0]+uvw[1]*facet_tri[i].connect[1]->xyz[0]+uvw[2]*facet_tri[i].connect[2]->xyz[0];
00256                         //double y = uvw[0]*facet_tri[i].connect[0]->xyz[1]+uvw[1]*facet_tri[i].connect[1]->xyz[1]+uvw[2]*facet_tri[i].connect[2]->xyz[1];
00257                         //double z = uvw[0]*facet_tri[i].connect[0]->xyz[2]+uvw[1]*facet_tri[i].connect[1]->xyz[2]+uvw[2]*facet_tri[i].connect[2]->xyz[2];
00258                         //std::cout << "\t\t\tx = " << x << "  y = " << y << "  z = " << z << std::endl;
00259                         return i;
00260 
00261                 }
00262         }
00263 
00264         return -1;
00265 }
00266 //project the point onto a plane determined by triangle
00267 void SurfProHarmonicMap::prjPtsToTri(Face tri, Vector3D pts, Vector3D nrml, Vector3D &xyz)
00268 {
00269         Vector3D origin = tri.connect[0]->xyz;
00270         Vector3D v = pts - origin;
00271         double dist = v[0]*nrml[0] + v[1]*nrml[1] + v[2]*nrml[2];
00272         xyz = origin + v - dist*nrml;
00273 }
00274 
00275 void SurfProHarmonicMap::setMeshData(vector<Vertex> &s, vector<Vertex> &t, vector<Face> &f){
00276         quad_mesh_src.insert(quad_mesh_src.begin(), s.begin(), s.end());
00277         quad_mesh_tgt.insert(quad_mesh_tgt.begin(), t.begin(), t.end());
00278         facelist.insert(facelist.begin(), f.begin(), f.end());  
00279 }
00280 
00281 void SurfProHarmonicMap::getMeshData(vector<Vertex> &v){
00282         int index = -1;
00283         for (vector<Vertex>::iterator it = quad_mesh_tgt.begin(); it != quad_mesh_tgt.end(); it++){
00284                 index++;
00285                 if (it->onBoundary)
00286                         continue;
00287                 v[index].xyz[0] = it->xyz[0];
00288                 v[index].xyz[1] = it->xyz[1];
00289                 v[index].xyz[2] = it->xyz[2];
00290         }
00291 
00292 }
00293 
00294 void SurfProHarmonicMap::match()
00295 {
00296         //outmost boundary will be distributed on the unit disk
00297         boundaryDistribution();
00298         ComputeWeight();
00299         //extra processing of interior loops
00300 
00301     HarmonicMapper hm_src(mk_core, src_facet_v, src_facet_tri, src_facet_e, adj_src);
00302         hm_src.execute();
00303         hm_src.getUV(src_facet_v);
00304         
00305         LocateBoundaryNodesTarget();
00306 
00307         HarmonicMapper hm_tgt(mk_core, tgt_facet_v, tgt_facet_tri, tgt_facet_e, adj_tgt);
00308         hm_tgt.execute();
00309         hm_tgt.getUV(tgt_facet_v);
00310 
00311         //test();
00312 }
00313 
00314 void SurfProHarmonicMap::LocateBoundaryNodesTarget()
00315 {
00316         std::set<int> set_edges, set_vertices;
00317         for (unsigned int i = 0; i < target.connEdges.size(); i++)
00318                 set_edges.insert(target.connEdges[i]->index);
00319         for (unsigned int i = 0; i < source.connect.size(); i++)
00320                 set_vertices.insert(source.connect[i]->index);
00321 
00322     //match the vertices between the source and target surface.
00323     //std::map<int,int> tgt_src_v, tgt_src_e;
00324         std::vector<int> tgt_src_v(vertices.size(), -1), tgt_src_e(edges.size(), -1);   
00325     int count = -1;
00326         std::set<int> bak_edges;
00327         bak_edges.insert(set_edges.begin(), set_edges.end());
00328     for (unsigned int mm = 0; mm < target.vertexloops.size(); mm++){
00329         for (unsigned int i = 0; i < target.vertexloops[mm].size(); i++){
00330                     count++;
00331                         if (target.vertexloops[mm].size()==1){
00332                                 tgt_src_v[target.vertexloops[mm][0]] = source.vertexloops[mm][0];
00333                                 continue;
00334                         }
00335         
00336                     bool is_found = false;
00337                     int edge_index = -1;
00338                     int pre_vertex = target.vertexloops[mm][i];
00339                     int next_vertex = -1;
00340                     while(!is_found){
00341                             vector<iBase_EntityHandle> adj;
00342                             //get the adjacent edges of a vertex
00343                             g_err = igeom_instance->getEntAdj(vertices[pre_vertex].gVertexHandle, iBase_EDGE, adj);
00344                             IBERRCHK(g_err, "Trouble get the adjacent edges around a vertex");
00345                             //find the edge on the linking surface                      
00346                             if (edge_index != -1){
00347                                     //get edges on the linking surfaces
00348                                     set_edges.clear();
00349                                     std::vector<iBase_EntityHandle> adj_faces;
00350                                     g_err = igeom_instance->getEntAdj(edges[edge_index].gEdgeHandle, iBase_FACE, adj_faces);
00351                                     IBERRCHK(g_err, "Trouble get the adjacent faces around an edge");
00352                                     for (std::vector<iBase_EntityHandle>::iterator it = adj_faces.begin(); it != adj_faces.end(); it++){
00353                                             vector<iBase_EntityHandle> adj_edges;
00354                                             g_err = igeom_instance->getEntAdj(*it, iBase_EDGE, adj_edges);
00355                                             IBERRCHK(g_err, "Trouble get the adjacent edges around a face!");
00356                                             for (vector<iBase_EntityHandle>::iterator it_e = adj_edges.begin(); it_e != adj_edges.end(); it_e++){
00357                                                     int intdata = -1;
00358                                                     g_err = igeom_instance->getIntData(*it_e, harmonic_surf_pro, intdata);
00359                                                     IBERRCHK(g_err, "Trouble get the int data for an edge");
00360                                                         if (set_edges.find(intdata)==set_edges.end())
00361                                                                 set_edges.insert(intdata);
00362                                             }
00363                                     }
00364                             }
00365                             //proceed to the next linking edge
00366                             is_found = true;
00367                             for (vector<iBase_EntityHandle>::iterator it = adj.begin(); it != adj.end(); it++){
00368                                     g_err = igeom_instance->getIntData(*it, harmonic_surf_pro, edge_index);
00369                                     IBERRCHK(g_err, "Trouble get the adjacent edges of a vertex!");
00370                                     if (set_edges.find(edge_index)==set_edges.end()){//this is an edge on the linking surface.
00371                                             is_found = is_found && false;
00372                                             break;
00373                                     }
00374                                     else{
00375                                                 edge_index = -1;
00376                                             is_found = is_found && true;
00377                                             continue;
00378                                     }
00379                             } 
00380                                 //check whether next vertex is on the source surface or not
00381                             if (is_found){
00382                                                                  
00383                                    break;
00384                                 }
00385                             //find the next vertex connected by the linking edge
00386                             if (edges[edge_index].connect[0]->index == pre_vertex)
00387                                     next_vertex = edges[edge_index].connect[1]->index;
00388                             else
00389                                     next_vertex = edges[edge_index].connect[0]->index;
00390                             pre_vertex = next_vertex;
00391                     }
00392                         tgt_src_v[target.vertexloops[mm][i]] = next_vertex;
00393                         set_edges.clear();
00394                         set_edges.insert(bak_edges.begin(), bak_edges.end());
00395         }
00396     }
00397         //done with the mapping of vertices between the source and target surfaces;
00398     //test the vertex mapping between source and target
00399     count = -1;
00400     for (unsigned int mm = 0; mm < target.vertexloops.size(); mm++){
00401         for (unsigned int i = 0; i < target.vertexloops[mm].size(); i++){
00402                         int tgt_edge_index = target.edgeloops[mm][i];
00403                         if (target.vertexloops[mm].size()==1){
00404                                 tgt_src_e[tgt_edge_index] = source.edgeloops[mm][0];
00405                                 continue;
00406                         }
00407                         int vtx_a = tgt_src_v[edges[tgt_edge_index].connect[0]->index], vtx_b = tgt_src_v[edges[tgt_edge_index].connect[1]->index];
00408                         for (unsigned int j = 0; j < source.connEdges.size(); j++){
00409                                 int vtx_c = (source.connEdges[j]->connect[0])->index, vtx_d = (source.connEdges[j]->connect[1])->index;
00410                                 if (((vtx_c == vtx_a)&&(vtx_d == vtx_b))||((vtx_c == vtx_b)&&(vtx_d == vtx_a))){
00411                                         tgt_src_e[tgt_edge_index] = source.connEdges[j]->index;
00412                                         break;
00413                                 }
00414                         } 
00415         }
00416     }
00417     
00418     count = -1;
00419     map<int, int> src_vtx_facet;
00420     vector<vector<int> > src_nodes_list;
00421     vector<vector<double> > src_nodes_u;
00422     src_nodes_list.resize(edges.size());
00423         src_nodes_u.resize(edges.size());
00424     for (unsigned int i = 0; i < source.edgeloops.size(); i++){
00425 
00426         for (unsigned int j = 0; j < source.edgeloops[i].size(); j++){
00427            count++;
00428                    double u;
00429            for (list<int>::iterator it = src_geom_facet[count+source.connect.size()].begin(); it != src_geom_facet[count+source.connect.size()].end(); it++){
00430                 src_nodes_list[source.edgeloops[i][j]].push_back(*it);
00431                 g_err = igeom_instance->getEntXYZtoU(edges[source.edgeloops[i][j]].gEdgeHandle, src_facet_v[*it].xyz[0], src_facet_v[*it].xyz[1], src_facet_v[*it].xyz[2], u);
00432                 IBERRCHK(g_err, "Trouble get the u coordinates!");
00433                 src_nodes_u[source.edgeloops[i][j]].push_back(u);
00434            }
00435         }
00436     }
00437         count = -1;
00438         for (unsigned int i = 0; i < source.vertexloops.size(); i++){
00439                 for (unsigned int j = 0; j < source.vertexloops[i].size(); j++){
00440                         count++;
00441                         src_vtx_facet[source.vertexloops[i][j]] = *(src_geom_facet[count].begin());
00442 
00443                 }
00444         }
00445     //match the boundary facet nodes between source and target
00446         count = -1;
00447     for (unsigned int i = 0; i < target.vertexloops.size(); i++){
00448         //loop i
00449                 vector<iBase_EntityHandle> tmp_src, tmp_tgt;
00450                 vector<double> dist_source, dist_target;
00451                 for (unsigned int j = 0; j < target.edgeloops[i].size(); j++){
00452                         tmp_tgt.push_back(edges[target.edgeloops[i][j]].gEdgeHandle);
00453                         tmp_src.push_back(edges[tgt_src_e[target.edgeloops[i][j]]].gEdgeHandle);
00454                 }
00455                 dist_source.resize(tmp_src.size());
00456                 g_err = igeom_instance->measure(&tmp_src[0], tmp_src.size(), &dist_source[0]);
00457                 IBERRCHK(g_err, "Trouble get the measure of geometric edges!");
00458                 
00459                 dist_target.resize(tmp_tgt.size());
00460                 g_err = igeom_instance->measure(&tmp_tgt[0], tmp_tgt.size(), &dist_target[0]);
00461                 IBERRCHK(g_err, "Trouble get the measure of geometric edges!");     
00462                 for (unsigned int j = 0; j < target.vertexloops[i].size(); j++){
00463                         count++;            
00464                         int current_tgt_vtx = target.vertexloops[i][j], current_src_vtx = tgt_src_v[current_tgt_vtx];
00465             int next_tgt_vtx = target.vertexloops[i][(j+1+target.vertexloops[i].size())%target.vertexloops[i].size()], next_src_vtx = tgt_src_v[next_tgt_vtx];
00466             int tgt_edge_index = target.edgeloops[i][j];
00467             int src_edge_index = tgt_src_e[tgt_edge_index];
00468                         std::cout << "current target vtx = " << current_tgt_vtx << "\tcurrent source vtx = " << current_src_vtx << "\tnext_tgt_vtx = " << next_tgt_vtx << "\tnext_src_vtx = " << next_src_vtx << "\ttarget edge index = " << tgt_edge_index << "\tsource edge index = " << std::endl;
00469             double u1, u2, u3, u4, u;
00470             g_err = igeom_instance->getEntXYZtoU(edges[tgt_edge_index].gEdgeHandle, vertices[current_tgt_vtx].xyz[0], vertices[current_tgt_vtx].xyz[1], vertices[current_tgt_vtx].xyz[2], u1);
00471             IBERRCHK(g_err, "Trouble get the parametric coordinate of a vertex on an edge!");
00472             g_err = igeom_instance->getEntXYZtoU(edges[tgt_edge_index].gEdgeHandle, vertices[next_tgt_vtx].xyz[0], vertices[next_tgt_vtx].xyz[1], vertices[next_tgt_vtx].xyz[2], u2);
00473             IBERRCHK(g_err, "Trouble get the parametric coordinate of a vertex on an edge!");
00474             g_err = igeom_instance->getEntXYZtoU(edges[src_edge_index].gEdgeHandle, vertices[current_src_vtx].xyz[0], vertices[current_src_vtx].xyz[1], vertices[current_src_vtx].xyz[2], u3);
00475             IBERRCHK(g_err, "Trouble get the parametric coordinate of a vertex on an edge!");
00476             g_err = igeom_instance->getEntXYZtoU(edges[src_edge_index].gEdgeHandle, vertices[next_src_vtx].xyz[0], vertices[next_src_vtx].xyz[1], vertices[next_src_vtx].xyz[2], u4);
00477             IBERRCHK(g_err, "Trouble get the parametric coordinate of a vertex on an edge!");
00478                         int count_pts = 0;
00479                         int tgt_vtx_index = *(tgt_geom_facet[count].begin());
00480                         tgt_facet_v[tgt_vtx_index].onBoundary = true;
00481                         tgt_facet_v[tgt_vtx_index].uv[0] = src_facet_v[src_vtx_facet[current_src_vtx]].uv[0];
00482                         tgt_facet_v[tgt_vtx_index].uv[1] = src_facet_v[src_vtx_facet[current_src_vtx]].uv[1];
00483                         std::cout << "src index = " << src_vtx_facet[current_src_vtx] << std::endl;
00484 
00485                         std::vector<double> lists_u;
00486                         std::vector<int> index_list;
00487                         int start_index = -1;
00488                         if (target.vertexloops[i].size() == 1){
00489                                 //get the source vertex facet
00490                                 std::vector<int> pointlist;
00491                                 pointlist.push_back(*(src_geom_facet[count].begin()));
00492                                 for (list<int>::iterator it = src_geom_facet[count+source.connect.size()].begin(); it != src_geom_facet[count+source.connect.size()].end(); it++)
00493                                         pointlist.push_back(*it);
00494                                 std::sort(pointlist.begin(), pointlist.end());
00495                                 start_index = find(pointlist.begin(), pointlist.end(), *(src_geom_facet[count].begin()))-pointlist.begin();
00496                                 for (std::vector<int>::iterator it = pointlist.begin(); it != pointlist.end(); it++){
00497                                         index_list.push_back(*it);
00498                                         g_err = igeom_instance->getEntXYZtoU(edges[src_edge_index].gEdgeHandle, src_facet_v[*it].xyz[0], src_facet_v[*it].xyz[1], src_facet_v[*it].xyz[2], u);
00499                             IBERRCHK(g_err, "Trouble get the parametric coordinate of a vertex on an edge!");
00500                                         lists_u.push_back(u);
00501                                 }
00502                         }                                               
00503                         for (list<int>::iterator it = tgt_geom_facet[count+target.connect.size()].begin(); it != tgt_geom_facet[count+target.connect.size()].end(); it++){
00504                                 count_pts++;
00505                 g_err = igeom_instance->getEntXYZtoU(edges[tgt_edge_index].gEdgeHandle, tgt_facet_v[*it].xyz[0], tgt_facet_v[*it].xyz[1], tgt_facet_v[*it].xyz[2], u);
00506                 IBERRCHK(g_err, "Trouble get the parametric coordinate of a vertex on an edge!");
00507                 
00508         
00509                 tgt_facet_v[*it].onBoundary = true;
00510                 //how to compute the corresponding uv coordinates of each nodes on the source H_1
00511                 //loop over facet nodes
00512                 //cases to be discussed: (1) no facet nodes (2) last facet nodes (3) first facet nodes
00513                                 int pre_facet_v = -1, next_facet_v = -1;
00514                                 double pre_facet_u = 0.0, next_facet_u = 0.0, u_src;
00515                                 if (src_nodes_u[src_edge_index].size() == 0){//no facet nodes
00516                                         //it will be a straight line between current_src_vtx and next_src_vtx
00517                                         pre_facet_v = src_vtx_facet[current_src_vtx];
00518                                         next_facet_v = src_vtx_facet[next_src_vtx];
00519                                         pre_facet_u = u3;
00520                                         next_facet_u = u4;
00521                                 }
00522                                 else if (target.vertexloops[i].size() == 1){//circle case
00523                                         double tgt_min, tgt_max, src_min, src_max;
00524                                         g_err = igeom_instance->getEntURange(edges[target.edgeloops[i][0]].gEdgeHandle, tgt_min, tgt_max);
00525                                         IBERRCHK(g_err, "Trouble get the parametric range of U on an edge!");
00526                                         g_err = igeom_instance->getEntURange(edges[source.edgeloops[i][0]].gEdgeHandle, src_min, src_max);
00527                                         IBERRCHK(g_err, "Trouble get the parametric range of U on an edge!");
00528                                         
00529                                         //get edge sense w.r.t. face
00530                                         int dotproduct = 0, sense_tgt, sense_src;
00531                                         g_err = igeom_instance->getEgFcSense(edges[target.edgeloops[i][0]].gEdgeHandle, target.gFaceHandle, sense_tgt);
00532                                         IBERRCHK(g_err, "Trouble get the edge sense w.r.t. the target surface!");
00533                                         g_err = igeom_instance->getEgFcSense(edges[source.edgeloops[i][0]].gEdgeHandle, source.gFaceHandle, sense_src);
00534                                         IBERRCHK(g_err, "Trouble get the edge sense w.r.t. the source surface!");
00535                                         dotproduct = sense_tgt*sense_src;
00536                                         if (dotproduct > 0)
00537                                                 u_src = (src_min - src_max)*(u - tgt_min)/(tgt_max - tgt_min) + src_max;
00538                                         else
00539                                                 u_src = (src_max - src_min)*(u - tgt_min)/(tgt_max - tgt_min) + src_min;
00540                                                 
00541                                         //note(very important):usually edge facets on the target surface are in the reverse order as those on the source surface. 
00542                                         int index_count = 0;
00543                                         //determine whether lists_u[start_index] has max or min value
00544                                         bool is_max = (lists_u[start_index] >= lists_u[(start_index+1)%lists_u.size()])&&(lists_u[start_index] >= lists_u[(start_index-1)%lists_u.size()]);
00545                                         if (lists_u[(start_index+1)%lists_u.size()] <= lists_u[(start_index+2)%lists_u.size()]){                                                
00546                                                 if ((!is_max)&&(u_src <= lists_u[start_index])){
00547                                                         pre_facet_u = src_min - src_max + lists_u[(start_index + lists_u.size()-1)%lists_u.size()];
00548                                                         next_facet_u = lists_u[start_index];
00549                                                         pre_facet_v = index_list[(start_index + lists_u.size()-1)%lists_u.size()];
00550                                                         next_facet_v = index_list[start_index];
00551                                                 }
00552                                                 else if ((!is_max)&&(u_src >= lists_u[(start_index + lists_u.size()-1)%lists_u.size()])){
00553                                                         pre_facet_u = lists_u[(start_index + lists_u.size()-1)%lists_u.size()];
00554                                                         next_facet_u = lists_u[start_index]-src_min+src_max;
00555                                                         pre_facet_v = index_list[(start_index + lists_u.size()-1)%lists_u.size()];
00556                                                         next_facet_v = index_list[start_index];
00557                                                 }
00558                                                 else{
00559                                                         int it_index = start_index;
00560                                                         if (is_max) it_index += 1;                                                                                      
00561                                                         for (; index_count < lists_u.size(); it_index++, index_count++){
00562                                                                 if (u_src <= lists_u[it_index%lists_u.size()]){
00563                                                                         pre_facet_u = lists_u[(it_index-1+lists_u.size())%lists_u.size()];
00564                                                                         next_facet_u = lists_u[it_index%lists_u.size()];
00565                                                                         pre_facet_v = index_list[(it_index-1+lists_u.size())%lists_u.size()];
00566                                                                         next_facet_v = index_list[it_index%lists_u.size()];
00567                                                                         break;
00568                                                                 }
00569                                                         }
00570                                                 }
00571                                         }
00572                                         else{//lists_u is decreasing
00573                                                 //determine if lists_u has max or min values
00574                                                 if ((!is_max)&&(u_src <= lists_u[start_index])){
00575                                                         pre_facet_u = src_min - src_max + lists_u[(start_index +1)%lists_u.size()];
00576                                                         next_facet_u = lists_u[start_index];
00577                                                         pre_facet_v = index_list[(start_index+1)%lists_u.size()];
00578                                                         next_facet_v = index_list[start_index];
00579                                                 }
00580                                                 else if ((!is_max)&&(u_src >= lists_u[(start_index +1)%lists_u.size()])){
00581                                                         pre_facet_u = lists_u[(start_index + 1)%lists_u.size()];
00582                                                         next_facet_u = lists_u[start_index]-src_min+src_max;
00583                                                         pre_facet_v = index_list[(start_index +1)%lists_u.size()];
00584                                                         next_facet_v = index_list[start_index];
00585                                                 }
00586                                                 else{
00587                                                         int it_index = start_index;
00588                                                         if (is_max) it_index -= 1;                                                                                      
00589                                                         for (; index_count < lists_u.size(); it_index--, index_count++){
00590                                                                 if (u_src <= lists_u[(it_index+lists_u.size())%lists_u.size()]){
00591                                                                         pre_facet_u = lists_u[(it_index+1+lists_u.size())%lists_u.size()];
00592                                                                         next_facet_u = lists_u[(it_index+lists_u.size())%lists_u.size()];
00593                                                                         pre_facet_v = index_list[(it_index+1+lists_u.size())%lists_u.size()];
00594                                                                         next_facet_v = index_list[(it_index+lists_u.size())%lists_u.size()];
00595                                                                         break;
00596                                                                 }
00597                                                         }
00598                                                 }
00599 
00600                                         }                                       
00601                                 }
00602                                 else{
00603                                         u_src = (u4 - u3)*(u - u1)/(u2 - u1) + u3;
00604                                         for (int m = src_nodes_u[src_edge_index].size()-1; m >= 0; m--){
00605                                                 if (u4 > u3){
00606                                                         if (m == 0)//it is located before the first item
00607                                                                 if (u_src >= src_nodes_u[src_edge_index][0]){
00608                                                                         next_facet_v = src_vtx_facet[next_src_vtx];
00609                                                                         next_facet_u = u4;
00610                                                                         pre_facet_v = src_nodes_list[src_edge_index][0];
00611                                                                         pre_facet_u = src_nodes_u[src_edge_index][0];
00612                                                                         break;          
00613                                                                 }
00614                                                         if (m == (src_nodes_u[src_edge_index].size()-1)){//it is located after the last item
00615                                                                 int tmp_index = src_nodes_list[src_edge_index].size() -1;
00616                                                                 if (u_src <= src_nodes_u[src_edge_index][tmp_index]){
00617                                                                         next_facet_v = src_nodes_list[src_edge_index][tmp_index];
00618                                                                         next_facet_u = src_nodes_u[src_edge_index][tmp_index];
00619                                                                         pre_facet_v = src_vtx_facet[current_src_vtx];
00620                                                                         pre_facet_u = u3;
00621                                                                         break;
00622                                                                 }
00623                                                         }                               
00624                                                         if (src_nodes_u[src_edge_index][m] >= u_src){//it is located between items
00625                                                                 if ((m+1)==src_nodes_list[src_edge_index].size()){
00626                                                                         pre_facet_v = src_vtx_facet[current_src_vtx];
00627                                                                         pre_facet_u = u3;
00628                                                                 }
00629                                                                 else{
00630                                                                         pre_facet_v = src_nodes_list[src_edge_index][m+1];
00631                                                                         pre_facet_u = src_nodes_u[src_edge_index][m+1];
00632                                                                 }
00633                                                                 
00634                                                                 next_facet_v = src_nodes_list[src_edge_index][m];
00635                                                                 next_facet_u = src_nodes_u[src_edge_index][m];
00636                                                                 break;
00637                                                         }
00638                                                         
00639                                                 }
00640                                                 else{
00641                                                         if (m == 0)//it is located before the first item
00642                                                                 if (u_src <= src_nodes_u[src_edge_index][0]){
00643                                                                         next_facet_v = src_vtx_facet[next_src_vtx];
00644                                                                         next_facet_u = u4;
00645                                                                         pre_facet_v = src_nodes_list[src_edge_index][0];
00646                                                                         pre_facet_u = src_nodes_u[src_edge_index][0];
00647                                                                         break;          
00648                                                                 }
00649                                                         if (m == (src_nodes_u[src_edge_index].size()-1)){//it is located after the last item
00650                                                                 int tmp_index = src_nodes_list[src_edge_index].size() -1;
00651                                                                 if (u_src >= src_nodes_u[src_edge_index][tmp_index]){
00652                                                                         next_facet_v = src_nodes_list[src_edge_index][tmp_index];
00653                                                                         next_facet_u = src_nodes_u[src_edge_index][tmp_index];
00654                                                                         pre_facet_v = src_vtx_facet[current_src_vtx];
00655                                                                         pre_facet_u = u3;
00656                                                                         break;
00657                                                                 }
00658                                                         }
00659                                                         if (src_nodes_u[src_edge_index][m] <= u_src){
00660                                                                 next_facet_v = src_nodes_list[src_edge_index][m];
00661                                                                 next_facet_u = src_nodes_u[src_edge_index][m];
00662                                                                 if ((m+1)==src_nodes_list[src_edge_index].size()){
00663                                                                         pre_facet_v = src_vtx_facet[current_src_vtx];
00664                                                                         pre_facet_u = u3;
00665                                                                 }
00666                                                                 else{
00667                                                                         pre_facet_v = src_nodes_list[src_edge_index][m+1];
00668                                                                         pre_facet_u = src_nodes_u[src_edge_index][m+1];
00669                                                                 }
00670                                                                 break;
00671                                                         }
00672                                                 }
00673                                         }
00674                                 }
00675                                 double alpha = (u_src - pre_facet_u)/(next_facet_u - pre_facet_u);                              
00676                                 tgt_facet_v[*it].uv[0] = (1.0 - alpha)*src_facet_v[pre_facet_v].uv[0]+alpha*src_facet_v[next_facet_v].uv[0];
00677                                 tgt_facet_v[*it].uv[1] = (1.0 - alpha)*src_facet_v[pre_facet_v].uv[1]+alpha*src_facet_v[next_facet_v].uv[1];
00678                                 std::cout << std::endl;
00679                         }
00680         }
00681     }
00682 }
00683 
00684 void SurfProHarmonicMap::ComputeWeight()
00685 {
00686         //create the facet mesh on the source surface
00687         int count = -1;
00688         for (vector<Vertex>::iterator it = src_facet_v.begin(); it != src_facet_v.end(); it++){
00689                 m_err = imesh_instance->createVtx(it->xyz[0], it->xyz[1], it->xyz[2], it->gVertexHandle);
00690                 IBERRCHK(m_err, "Trouble create a mesh vertex!");
00691                 count++;
00692                 m_err = imesh_instance->setIntData(it->gVertexHandle, facet_mesh_tag, count);
00693                 IBERRCHK(m_err, "Trouble set the int data for a facet vertex");
00694         }
00695         //create the facet edge entities on the source surface
00696         vector<set<int> > edge_connect;
00697         edge_connect.resize(src_facet_v.size());
00698         count = -1;
00699         for (vector<Face>::iterator it = src_facet_tri.begin(); it != src_facet_tri.end(); it++){
00700                 int a = it->connect[0]->index, b = it->connect[1]->index, c = it->connect[2]->index;
00701                 if (edge_connect[a].find(b)==edge_connect[a].end())
00702                         addEdgeToList(a, b, count, edge_connect, src_facet_e, src_facet_v);
00703                 if (edge_connect[a].find(c)==edge_connect[a].end())
00704                         addEdgeToList(a, c, count, edge_connect, src_facet_e, src_facet_v);
00705                 if (edge_connect[b].find(c)==edge_connect[b].end())
00706                         addEdgeToList(b, c, count, edge_connect, src_facet_e, src_facet_v);
00707         }
00708         size_src_e = src_facet_e.size();
00709         //create the facet face entities on the source surface
00710         count = -1;
00711         for (vector<Face>::iterator it = src_facet_tri.begin(); it != src_facet_tri.end(); it++){
00712                 vector<iBase_EntityHandle> nodes;
00713                 for (unsigned int i = 0; i < it->connect.size(); i++)
00714                         nodes.push_back((it->connect[i])->gVertexHandle);
00715                 m_err = imesh_instance->createEnt(iMesh_TRIANGLE, &nodes[0], nodes.size(), it->gFaceHandle);
00716                 IBERRCHK(m_err, "Trouble create a triangle mesh entity on the source surface!");
00717                 count++;
00718                 m_err = imesh_instance->setIntData(it->gFaceHandle, facet_mesh_tag, count);
00719                 IBERRCHK(m_err, "Trouble set the int data for a facet vertex");
00720         }
00721         //add extra stuff in order to seal the interior loops
00722         addExtra(source, src_facet_v, src_facet_e, src_facet_tri, src_geom_facet, size_src_v);
00723 
00724         //get the adjacent edges of every vertex on the source surface
00725         adj_src.resize(src_facet_v.size());
00726         for (vector<Vertex>::iterator it = src_facet_v.begin(); it != src_facet_v.end(); it++){
00727                 vector<iBase_EntityHandle> adj_edges;
00728                 m_err = imesh_instance->getEntAdj(it->gVertexHandle, iBase_EDGE, adj_edges);
00729                 IBERRCHK(m_err, "Trouble get the adjacent edges around a vertex!");
00730                 int intdata = -1;
00731                 for (vector<iBase_EntityHandle>::iterator it_e = adj_edges.begin(); it_e != adj_edges.end(); it_e++){
00732                         m_err = imesh_instance->getIntData(*it_e, facet_mesh_tag, intdata);
00733                         IBERRCHK(m_err, "Trouble get the int data for the adjacent edges!");
00734                         adj_src[it->index].insert(intdata);
00735                 }
00736         }
00737         //compute the weight for edges on the source surface
00738         computeEdgeWeight(src_facet_e, src_facet_tri, src_facet_v);
00739 
00740         //create the facet mesh on the target surface
00741         count = -1;
00742         for (vector<Vertex>::iterator it = tgt_facet_v.begin(); it != tgt_facet_v.end(); it++){
00743                 m_err = imesh_instance->createVtx(it->xyz[0], it->xyz[1], it->xyz[2], it->gVertexHandle);
00744                 IBERRCHK(m_err, "Trouble create a mesh vertex");
00745                 count++;
00746                 m_err = imesh_instance->setIntData(it->gVertexHandle, facet_mesh_tag, count);
00747                 IBERRCHK(m_err, "Trouble set the int data for a facet vertex");
00748         }
00749         //get the facet edges on the target surface
00750         edge_connect.clear();
00751         edge_connect.resize(tgt_facet_v.size());
00752         count = -1;
00753         for (vector<Face>::iterator it = tgt_facet_tri.begin(); it != tgt_facet_tri.end(); it++){
00754                 int a = it->connect[0]->index, b = it->connect[1]->index, c = it->connect[2]->index;
00755                 if (edge_connect[a].find(b)==edge_connect[a].end())
00756                         addEdgeToList(a, b, count, edge_connect, tgt_facet_e, tgt_facet_v);
00757                 if (edge_connect[a].find(c)==edge_connect[a].end())
00758                         addEdgeToList(a, c, count, edge_connect, tgt_facet_e, tgt_facet_v);
00759                 if (edge_connect[b].find(c)==edge_connect[b].end())
00760                         addEdgeToList(b, c, count, edge_connect, tgt_facet_e, tgt_facet_v);
00761         }
00762         size_tgt_e = tgt_facet_e.size();
00763         //get the facet faces on the target surface
00764         count = -1;
00765         for (vector<Face>::iterator it = tgt_facet_tri.begin(); it != tgt_facet_tri.end(); it++){
00766                 vector<iBase_EntityHandle> nodes;
00767                 for (unsigned int i = 0; i < it->connect.size(); i++)
00768                         nodes.push_back((it->connect[i])->gVertexHandle);
00769                 m_err = imesh_instance->createEnt(iMesh_TRIANGLE, &nodes[0], nodes.size(), it->gFaceHandle);
00770                 IBERRCHK(m_err, "Trouble create a triangle mesh entity on the target surface!");
00771                 count++;
00772                 m_err = imesh_instance->setIntData(it->gFaceHandle, facet_mesh_tag, count);
00773                 IBERRCHK(m_err, "Trouble set the int data for a facet vertex!");
00774         }
00775         //add extra stuff in order to seal the interior loops
00776         //addExtra(target, tgt_facet_v, tgt_facet_e, tgt_facet_tri, tgt_geom_facet, size_tgt_v);
00777         //get the adjacent edges of every vertex on the target surface
00778         adj_tgt.resize(tgt_facet_v.size());
00779         for (vector<Vertex>::iterator it = tgt_facet_v.begin(); it != tgt_facet_v.end(); it++){
00780                 vector<iBase_EntityHandle> adj_edges;
00781                 m_err = imesh_instance->getEntAdj(it->gVertexHandle, iBase_EDGE, adj_edges);
00782                 IBERRCHK(m_err, "Trouble get the adjacent edges around a vertex!");
00783                 int intdata = -1;
00784                 for (vector<iBase_EntityHandle>::iterator it_e = adj_edges.begin(); it_e != adj_edges.end(); it_e++){
00785                         m_err = imesh_instance->getIntData(*it_e, facet_mesh_tag, intdata);
00786                         IBERRCHK(m_err, "Trouble get the int data for the adjacent edges!");
00787                         adj_tgt[it->index].insert(intdata);
00788                 }
00789         }
00790         //compute the edge weight on the target surface
00791         computeEdgeWeight(tgt_facet_e, tgt_facet_tri, tgt_facet_v);
00792 
00793 }
00794 
00795 void SurfProHarmonicMap::addExtra(Face f, vector<Vertex> &facet_v, vector<Edge> &facet_e, vector<Face> &facet_tri, vector<list<int> > geom_facet, int size_facet_v)
00796 {
00797         //add extra nodes to the list
00798         if (f.vertexloops.size() > 1){
00799                 int count = 0;
00800                 for (unsigned int i = 1; i < f.vertexloops.size(); i++){
00801                         double xyz[3] = {0.0, 0.0, 0.0};
00802                         count += f.vertexloops[i-1].size();
00803                         int vtx_center = size_facet_v+i-1;
00804                         vector<int> pointlist;
00805                         for (unsigned int j = 0; j < f.vertexloops[i].size(); j++){
00806                                 pointlist.push_back(*(geom_facet[count+j].begin()));
00807                                 for (int m = 0; m < 3; m++)
00808                                         xyz[m] += facet_v[*(geom_facet[count+j].begin())].xyz[m];
00809                                 for (list<int>::iterator it = geom_facet[count+j+f.connect.size()].begin(); it != geom_facet[count+j+f.connect.size()].end(); it++){
00810                                         pointlist.push_back(*it);
00811                                         for (int m = 0; m < 3; m++)
00812                                                 xyz[m] += facet_v[*it].xyz[m];
00813                                 }
00814                         }
00815                         //done
00816                         std::sort(pointlist.begin(), pointlist.end());
00817                         //add extra vertex
00818                         for (int m = 0; m < 3; m++)
00819                                  facet_v[vtx_center].xyz[m] = xyz[m]/double(pointlist.size());
00820                         m_err = imesh_instance->createVtx(facet_v[vtx_center].xyz[0], facet_v[vtx_center].xyz[1], facet_v[vtx_center].xyz[2], facet_v[vtx_center].gVertexHandle);
00821                         IBERRCHK(m_err, "Trouble create a facet vertex on the interior loops of source surface!");
00822                         m_err = imesh_instance->setIntData(facet_v[vtx_center].gVertexHandle, facet_mesh_tag, facet_v[vtx_center].index);
00823                         IBERRCHK(m_err, "Trouble set the int data for a facet vertex");
00824                         facet_v[vtx_center].index = vtx_center;
00825                         facet_v[vtx_center].onBoundary = false;
00826                         //add extra facet edges and adjacent edges
00827                         int ent_index = src_facet_e.size();
00828                         vector<iBase_EntityHandle> tmp;                 
00829                         for (vector<int>::iterator it = pointlist.begin(); it != pointlist.end(); it++, ent_index++){
00830                                 //std::cout << "test ent_index = " << ent_index << std::endl;
00831                                 Edge tmp_edge;          
00832                                 tmp_edge.index = ent_index;
00833                                 tmp_edge.connect.resize(2);
00834                                 tmp_edge.connect[0] = &facet_v[*it];
00835                                 tmp_edge.connect[1] = &facet_v[vtx_center];
00836                                 tmp.clear();
00837                                 tmp.push_back(facet_v[*it].gVertexHandle);
00838                                 tmp.push_back(facet_v[vtx_center].gVertexHandle);
00839                                 m_err = imesh_instance->createEnt(iMesh_LINE_SEGMENT, &tmp[0], 2, tmp_edge.gEdgeHandle);
00840                                 IBERRCHK(m_err, "Trouble create a extra facet edge on the interior loops of source surface!");
00841                                 m_err = imesh_instance->setIntData(tmp_edge.gEdgeHandle, facet_mesh_tag, tmp_edge.index);
00842                                 IBERRCHK(m_err, "Trouble set the int data for a facet vertex");
00843                                 
00844                                 facet_e.push_back(tmp_edge);
00845                         }
00846                         //add extra facet faces
00847                         ent_index = facet_tri.size();
00848                         //src_facet_tri.resize(ent_index+pointlist.size());
00849                         for (unsigned int j = 0; j < pointlist.size(); j++, ent_index++){
00850                                 Face tmp_f;
00851                                 int vtx_index[3];
00852                                 if (j == (pointlist.size()-1)){
00853                                         vtx_index[0] = pointlist[j];
00854                                         vtx_index[1] = pointlist[0];
00855                                 }
00856                                 else{
00857                                         vtx_index[0] = pointlist[j];
00858                                         vtx_index[1] = pointlist[j+1];
00859                                 }
00860                                 vtx_index[2] = vtx_center;
00861                                 tmp_f.connect.resize(3);
00862                                 tmp.clear();
00863                                 for (int m = 0; m < 3; m++){
00864                                         tmp_f.connect[m] = &facet_v[vtx_index[m]];
00865                                         tmp.push_back(facet_v[vtx_index[m]].gVertexHandle);
00866                                 }
00867                                 tmp_f.index = ent_index;
00868                                 m_err = imesh_instance->createEnt(iMesh_TRIANGLE, &tmp[0], tmp.size(), tmp_f.gFaceHandle);
00869                                 IBERRCHK(m_err, "Trouble create a extra triangle facet face entity!");
00870                                 m_err = imesh_instance->setIntData(tmp_f.gFaceHandle, facet_mesh_tag, tmp_f.index);
00871                                 IBERRCHK(m_err, "Trouble set the int data for a facet vertex");
00872                                 facet_tri.push_back(tmp_f);
00873                         }                       
00874                 }
00875         }
00876 
00877 
00878 
00879 }
00880 
00881 void SurfProHarmonicMap::computeEdgeWeight(vector<Edge> &f_edges, vector<Face> &f, vector<Vertex> f_v)
00882 {
00883         //compute the weight 
00884         for (vector<Edge>::iterator it = f_edges.begin(); it != f_edges.end(); it++){
00885                 vector<iBase_EntityHandle> adj_faces;
00886                 m_err = imesh_instance->getEntAdj(it->gEdgeHandle, iBase_FACE, adj_faces);
00887                 IBERRCHK(m_err, "Trouble get the adjacent faces for an facet edge!");           
00888                 assert(adj_faces.size()<=2);
00889                 double weight = 0.0;
00890                 int a = (it->connect[0])->index, b = (it->connect[1])->index;
00891                 for (vector<iBase_EntityHandle>::iterator it_ent = adj_faces.begin(); it_ent != adj_faces.end(); it_ent++){
00892                         int intdata = -1, c = -1;
00893                         m_err = imesh_instance->getIntData(*it_ent, facet_mesh_tag, intdata);
00894                         IBERRCHK(m_err, "Trouble get the int data for the adjacent face!");
00895                         if (((f[intdata].connect[0]->index == a)&&(f[intdata].connect[1]->index==b))||((f[intdata].connect[1]->index == a)&&(f[intdata].connect[0]->index==b)))
00896                                 c = f[intdata].connect[2]->index;
00897                         else if (((f[intdata].connect[0]->index == a)&&(f[intdata].connect[2]->index==b))||((f[intdata].connect[2]->index == a)&&(f[intdata].connect[0]->index==b)))
00898                                 c = f[intdata].connect[1]->index;
00899                         else if (((f[intdata].connect[1]->index == a)&&(f[intdata].connect[2]->index==b))||((f[intdata].connect[2]->index == a)&&(f[intdata].connect[1]->index==b)))
00900                                 c = f[intdata].connect[0]->index;
00901                         else
00902                                 continue;
00903                         double vec1[3], vec2[3];
00904                         for (int k = 0; k < 3; k++){
00905                                 vec1[k] = f_v[a].xyz[k] - f_v[c].xyz[k];
00906                                 vec2[k] = f_v[b].xyz[k] - f_v[c].xyz[k];
00907                         }
00908                         double cos_value = (vec1[0]*vec2[0] + vec1[1]*vec2[1] + vec1[2]*vec2[2])/sqrt((vec1[0]*vec1[0]+vec1[1]*vec1[1]+vec1[2]*vec1[2])*(vec2[0]*vec2[0]+vec2[1]*vec2[1]+vec2[2]*vec2[2]));
00909                         //std::cout << "ctan(angle) = " << cos_value/sqrt(1-pow(cos_value,2)) << std::endl;
00910                         weight += cos_value/sqrt(1-pow(cos_value,2));
00911                 }
00912                 it->e = 0.5*weight;
00913         }
00914 
00915 }
00916 
00917 void SurfProHarmonicMap::addEdgeToList(int a, int b, int &count, vector<set<int> > &edge_connect, vector<Edge> &f_edges, vector<Vertex> &facet_v)
00918 {
00919         count++;
00920         f_edges.resize(count+1);
00921         edge_connect[a].insert(b);
00922         edge_connect[b].insert(a);
00923         f_edges[count].connect.resize(2);
00924         f_edges[count].connect[0] = &facet_v[a];
00925         f_edges[count].connect[1] = &facet_v[b];
00926         f_edges[count].index = count;
00927         vector<iBase_EntityHandle> nodes;
00928         nodes.push_back(facet_v[a].gVertexHandle);
00929         nodes.push_back(facet_v[b].gVertexHandle);
00930         m_err = imesh_instance->createEnt(iMesh_LINE_SEGMENT, &nodes[0], nodes.size(), f_edges[count].gEdgeHandle);
00931         IBERRCHK(m_err, "Trouble create a facet line segment!");
00932         m_err = imesh_instance->setIntData(f_edges[count].gEdgeHandle, facet_mesh_tag, count);
00933         IBERRCHK(m_err, "Trouble set the int data for a facet edge!");
00934         
00935 }
00936 
00937 void SurfProHarmonicMap::boundaryDistribution()
00938 {
00939         //loop over the geometry entities on the outmost boundary
00940         //1. get the total distance of the outmost boundary loop
00941         vector<iBase_EntityHandle> tmp;
00942         vector<double> dist_source, dist_target;
00943         double total_dist_source = 0.0, total_dist_target = 0.0;
00944         for (unsigned int i = 0; i < source.edgeloops[0].size(); i++)
00945                 tmp.push_back(edges[source.edgeloops[0][i]].gEdgeHandle);
00946         dist_source.resize(tmp.size());
00947         g_err = igeom_instance->measure(&tmp[0], tmp.size(), &dist_source[0]);
00948         IBERRCHK(g_err, "Trouble get the measure of geometric edges!");
00949         tmp.clear();
00950         for (unsigned int i = 0; i < target.edgeloops[0].size(); i++)
00951                 tmp.push_back(edges[target.edgeloops[0][i]].gEdgeHandle);
00952         dist_target.resize(tmp.size());
00953         g_err = igeom_instance->measure(&tmp[0], tmp.size(), &dist_target[0]);
00954         IBERRCHK(g_err, "Trouble get the measure of geometric edges!");
00955         
00956         //Initialization of u,v coordinates
00957         for (unsigned int i = 0; i < src_facet_v.size(); i++){
00958                 src_facet_v[i].uv[0] = 0.0;
00959                 src_facet_v[i].uv[1] = 0.0;
00960         }
00961         for (unsigned int i = 0; i < tgt_facet_v.size(); i++){
00962                 tgt_facet_v[i].uv[0] = 0.0;
00963                 tgt_facet_v[i].uv[1] = 0.0;
00964         }               
00965 
00966         //assign u,v coordinates for facet vertices on the outmost boundary onto the unit disk, both source and target surface
00967         for (vector<Vertex>::iterator it = src_facet_v.begin(); it != src_facet_v.end(); it++)
00968                 it->onBoundary = false;
00969         for (vector<Vertex>::iterator it = tgt_facet_v.begin(); it != tgt_facet_v.end(); it++)
00970                 it->onBoundary = false;
00971         boundaryUnitDisk(source, dist_source, src_facet_v, src_geom_facet);
00972 }
00973 
00974 
00975 void SurfProHarmonicMap::boundaryUnitDisk(Face f, vector<double> dist, vector<Vertex> &facet_v, std::vector<std::list<int> > geom_facet)
00976 {
00977         double local_dist = 0.0, dist_sum = 0.0, len_curve = 0.0;
00978         int start_vertex = f.vertexloops[0][0];
00979         for (std::vector<double>::iterator it = dist.begin(); it != dist.end(); it++)
00980                 dist_sum += *it;
00981         vector<int> pointlist;
00982         for (unsigned int i = 0; i < f.vertexloops[0].size(); i++){
00983                 pointlist.push_back(*(geom_facet[i].begin()));
00984                 for (std::list<int>::iterator it = geom_facet[i+f.connect.size()].begin(); it != geom_facet[i+f.connect.size()].end(); it++)
00985                         pointlist.push_back(*it);
00986         }
00987         std::sort(pointlist.begin(), pointlist.end());
00988         for (unsigned int i = 0; i < pointlist.size(); i++){
00989                 int now = pointlist[i];
00990                 if (i > 0){
00991                         double tmp_dist = sqrt(pow(facet_v[now].xyz[0]-facet_v[pointlist[i-1]].xyz[0],2)+pow(facet_v[now].xyz[1]-facet_v[pointlist[i-1]].xyz[1],2)+pow(facet_v[now].xyz[2]-facet_v[pointlist[i-1]].xyz[2],2));          
00992                         local_dist += tmp_dist;
00993                 }
00994                 //get the facet on a geometry vertex of a loop  
00995                 facet_v[now].uv[0] = 100.0*cos(2*PI*local_dist/dist_sum);
00996                 facet_v[now].uv[1] = 100.0*sin(2*PI*local_dist/dist_sum);
00997                 facet_v[now].onBoundary = true;
00998         }
00999 }
01000 
01001 void SurfProHarmonicMap::getFacets()
01002 {
01003         SimpleArray<double> src_coords, tgt_coords;
01004         SimpleArray<int> src_facets, tgt_facets;
01005         double dist_tolerance = 1.0e-1;
01006         int err;
01007         //facets on the source surface
01008         iGeom_getFacets(igeom_instance->instance(), source.gFaceHandle, dist_tolerance, ARRAY_INOUT(src_coords), ARRAY_INOUT(src_facets), &err);
01009         assert(!err);
01010         src_facet_v.resize(src_coords.size()/3+source.vertexloops.size()-1);
01011         src_facet_tri.resize(src_facets.size()/3);
01012         size_src_v = src_coords.size()/3;
01013         size_src_f = src_facets.size()/3;
01014         
01015         for (unsigned int i = 0; i < src_coords.size(); i += 3){
01016                 src_facet_v[i/3].index = i/3;
01017                 for (int k = 0; k < 3; k++)
01018                         src_facet_v[i/3].xyz[k] = src_coords[i+k];
01019         }
01020         for (unsigned int i = 0; i < src_facets.size(); i += 3){
01021                 src_facet_tri[i/3].index = i/3;
01022                 src_facet_tri[i/3].connect.resize(3);
01023                 for (int k = 0; k < 3; k++)
01024                         src_facet_tri[i/3].connect[k] = &src_facet_v[src_facets[i+k]];
01025         }
01026         //facets on the target surface
01027         iGeom_getFacets(igeom_instance->instance(), target.gFaceHandle, dist_tolerance, ARRAY_INOUT(tgt_coords), ARRAY_INOUT(tgt_facets), &err);
01028         assert(!err);
01029         //tgt_facet_v.resize(tgt_coords.size()/3+target.vertexloops.size()-1);
01030         tgt_facet_v.resize(tgt_coords.size()/3);
01031         tgt_facet_tri.resize(tgt_facets.size()/3);
01032         size_tgt_v = tgt_coords.size()/3;
01033         size_tgt_f = tgt_facet_tri.size();
01034         for (unsigned int i = 0; i < tgt_coords.size(); i+=3){
01035                 tgt_facet_v[i/3].index = i/3;
01036                 for (int k = 0; k < 3; k++)
01037                         tgt_facet_v[i/3].xyz[k] = tgt_coords[i+k];
01038                 tgt_facet_v[i/3].uv[0] = 0.0;
01039                 tgt_facet_v[i/3].uv[1] = 0.0;
01040         }
01041         for (unsigned int i = 0; i < tgt_facets.size(); i += 3){
01042                 tgt_facet_tri[i/3].index = i/3;
01043                 tgt_facet_tri[i/3].connect.resize(3);
01044                 for (int k = 0; k < 3; k++)
01045                         tgt_facet_tri[i/3].connect[k] = &tgt_facet_v[tgt_facets[i+k]];
01046 
01047         }
01048 
01049         //match the facet mesh with geometry(source surface and target surface)
01050         //map-->vertices, edges, faces on the source surface
01051         
01052         std::cout << "starting the source surface mapping\n";
01053         MapFacetGeom(source, src_facet_v, src_facet_geom, src_geom_facet, size_src_v);
01054         std::cout << "starting the target surface mapping\n";
01055         MapFacetGeom(target, tgt_facet_v, tgt_facet_geom, tgt_geom_facet, size_tgt_v);  
01056 }
01057 
01058 void SurfProHarmonicMap::MapFacetGeom(Face f_surf, vector<Vertex> facet_node, std::map<int, int> &map_data, vector<list<int> > &geom_facet, int size_facet_v)
01059 {
01060         geom_facet.resize(f_surf.connect.size()+f_surf.connEdges.size()+1);
01061         
01062         for (unsigned int i = 0; i < size_facet_v; i++){
01063                 facet_node[i].index = i;
01064                 int is_on = false;
01065 
01066                 map_data[i] = -1;               
01067 
01068                 iGeom_isPositionOn(igeom_instance->instance(), f_surf.gFaceHandle, facet_node[i].xyz[0], facet_node[i].xyz[1], facet_node[i].xyz[2], &is_on);
01069                 if (!is_on){            
01070                         double proj_xyz[3];
01071                         g_err = igeom_instance->getEntClosestPtTrimmed(f_surf.gFaceHandle, facet_node[i].xyz[0], facet_node[i].xyz[1], facet_node[i].xyz[2], proj_xyz[0], proj_xyz[1], proj_xyz[2]);
01072                         IBERRCHK(g_err, "Trouble get the closest point on a geometry entity!");
01073                 }
01074                 is_on = false;
01075 
01076                 //check whether the facet mesh is on vertices or not.
01077                 for (unsigned int j = 0; j < f_surf.connect.size(); j++){                       
01078                         iGeom_isPositionOn(igeom_instance->instance(), f_surf.connect[j]->gVertexHandle, facet_node[i].xyz[0], facet_node[i].xyz[1], facet_node[i].xyz[2], &is_on);                     
01079                         if (is_on){
01080                                 map_data[i] = j;
01081                                 geom_facet[j].push_back(i);             
01082                                 break;
01083                         }
01084                 }
01085                 if (is_on)
01086                         continue;
01087                 //check whether the facet mesh is on edges or not
01088                 for (unsigned int j = 0; j < f_surf.connEdges.size(); j++){
01089                         Vector3D xyz_f(0.0), xyz_e(0.0);
01090                         g_err = igeom_instance->getEntClosestPtTrimmed((f_surf.connEdges[j])->gEdgeHandle, facet_node[i].xyz[0], facet_node[i].xyz[1], facet_node[i].xyz[2], xyz_e[0], xyz_e[1], xyz_e[2]);
01091                         IBERRCHK(g_err, "Trouble get the trimmed closest point on an edge!");
01092                         g_err = igeom_instance->getEntClosestPtTrimmed(f_surf.gFaceHandle, facet_node[i].xyz[0], facet_node[i].xyz[1], facet_node[i].xyz[2], xyz_f[0], xyz_f[1], xyz_f[2]);
01093                         IBERRCHK(g_err, "Trouble get the trimmed closest point on a face!");
01094                         Vector3D diff = xyz_f - xyz_e;
01095                         double dist_tol = sqrt(diff[0]*diff[0]+diff[1]*diff[1]+diff[2]*diff[2]);
01096                         if ( dist_tol < dist_tolerance){
01097                                 iGeom_isPositionOn(igeom_instance->instance(), f_surf.connEdges[j]->gEdgeHandle, xyz_e[0], xyz_e[1], xyz_e[2], &is_on);
01098                         }
01099                         if (is_on){
01100                                 map_data[i] = f_surf.connect.size()+j;
01101                                 geom_facet[f_surf.connect.size()+j].push_back(i);                               
01102                                 break;
01103                         }
01104                 }
01105                 if (is_on)
01106                         continue;
01107                 //check whether the facet mesh is on the surface or not
01108                 iGeom_isPositionOn(igeom_instance->instance(), f_surf.gFaceHandle, facet_node[i].xyz[0], facet_node[i].xyz[1], facet_node[i].xyz[2], &is_on);           
01109                 if (is_on){
01110                         map_data[i] = f_surf.connect.size()+f_surf.connEdges.size();
01111                         geom_facet[f_surf.connect.size()+f_surf.connEdges.size()].push_back(i);
01112                 }
01113                 else{
01114                         std::cout << "Something is wrong for the facet mesh!\txyz = {" << facet_node[i].xyz[0] << "," << facet_node[i].xyz[1] << "," << facet_node[i].xyz[2] << "}\n";
01115                         
01116                 }
01117         }
01118 }
01119 
01120 void SurfProHarmonicMap::preprocessing()
01121 {
01122         vector<iBase_EntityHandle> v;
01123         vector<iBase_EntityHandle> e;
01124         vector<iBase_EntityHandle> f;
01125         g_err = igeom_instance->getEntAdj(volume, iBase_VERTEX, v);
01126         IBERRCHK(g_err, "Trouble get the adjacent vertices on a volume!");
01127         g_err = igeom_instance->getEntAdj(volume, iBase_EDGE, e);
01128         IBERRCHK(g_err, "Trouble get the adjacent edges on a volume!");
01129         g_err = igeom_instance->getEntAdj(volume, iBase_FACE, f);
01130         IBERRCHK(g_err, "Trouble get the adjacent faces on a volume!");
01131         vertices.resize(v.size());
01132         for (unsigned int i = 0; i < v.size(); i++){
01133                 vertices[i].gVertexHandle = v[i];
01134                 g_err = igeom_instance->getVtxCoord(v[i], vertices[i].xyz[0], vertices[i].xyz[1], vertices[i].xyz[2]);
01135                 IBERRCHK(g_err, "Trouble get the xyz coordinates!");
01136                 vertices[i].index = i;
01137                 g_err = igeom_instance->getIntData(vertices[i].gVertexHandle, global_geom_tag, vertices[i].id);
01138                 IBERRCHK(g_err, "Trouble get the int data for vertex entity");
01139                 g_err = igeom_instance->setIntData(vertices[i].gVertexHandle, harmonic_surf_pro, i);
01140                 IBERRCHK(g_err, "Trouble set the int data for vertex entity");  
01141         }
01142         edges.resize(e.size());
01143         for (unsigned int i = 0; i < e.size(); i++){
01144                 edges[i].index = i;
01145                 edges[i].gEdgeHandle = e[i];
01146                 g_err = igeom_instance->getIntData(edges[i].gEdgeHandle, global_geom_tag, edges[i].id);
01147                 IBERRCHK(g_err, "Trouble get the int data for geometrical edge entity!");
01148                 vector<iBase_EntityHandle> adjs;
01149                 g_err = igeom_instance->getEntAdj(e[i], iBase_VERTEX, adjs);
01150                 IBERRCHK(g_err, "Trouble get the adjacent vertices of a geometrical edge!");
01151                 for (unsigned int j = 0; j < adjs.size(); j++){
01152                         int tmp;
01153                         g_err = igeom_instance->getIntData(adjs[j], harmonic_surf_pro, tmp);
01154                         IBERRCHK(g_err, "Trouble get the int data for a geometrical vertex!");
01155                         edges[i].connect.push_back(&vertices[tmp]);
01156                 }
01157                 g_err = igeom_instance->setIntData(edges[i].gEdgeHandle, harmonic_surf_pro, i);
01158                 IBERRCHK(g_err, "Trouble set the int data for a geometrical edge!");
01159         }
01160         link.resize(f.size());
01161         for (unsigned int i = 0; i < f.size(); i++){
01162                 if (f[i] == source.gFaceHandle){
01163                         addFaceToList(f[i], source, i, false);
01164                         GetGeomLoops(source, source.vertexloops, source.edgeloops);
01165                         postProcessGeomLoops(source);
01166                         adjustVtxEdges(source);
01167                         source.src_tgt_link = 0;
01168                         link[i].index = -1;
01169                         continue;
01170                 }                       
01171                 if (f[i] == target.gFaceHandle){
01172                         addFaceToList(f[i], target, i, false);
01173                         GetGeomLoops(target, target.vertexloops, target.edgeloops);
01174                         postProcessGeomLoops(target);
01175                         adjustVtxEdges(target);
01176                         target.src_tgt_link = 1;
01177                         link[i].index = -2;
01178                         continue;
01179                 }
01180                 link[i].src_tgt_link = 2;               
01181                 addFaceToList(f[i], link[i], i, true);
01182 
01183         }
01184         
01185                 
01186 }
01187 
01188 void SurfProHarmonicMap::adjustVtxEdges(Face &f)
01189 {
01190         vector<int> tmp;
01191         for (unsigned int mm = 0; mm < f.vertexloops.size(); mm++)
01192                 for (unsigned int nn = 0; nn < f.vertexloops[mm].size(); nn++)
01193                         tmp.push_back(f.vertexloops[mm][nn]);
01194         f.connect.clear();
01195         for (unsigned mm = 0; mm < tmp.size(); mm++)
01196                 f.connect.push_back(&vertices[tmp[mm]]);
01197         tmp.clear();
01198         for (unsigned int mm = 0; mm < f.edgeloops.size(); mm++)
01199                 for (unsigned int nn = 0; nn < f.edgeloops[mm].size(); nn++)
01200                         tmp.push_back(f.edgeloops[mm][nn]);
01201         f.connEdges.clear();
01202         for (unsigned mm = 0; mm < tmp.size(); mm++)
01203                 f.connEdges.push_back(&edges[tmp[mm]]);
01204 }
01205 
01206 void SurfProHarmonicMap::addFaceToList(iBase_EntityHandle entity, Face& f, int index, bool is_set_int)
01207 {
01208         f.gFaceHandle = entity;
01209         if (is_set_int){
01210                 g_err = igeom_instance->setIntData(f.gFaceHandle, harmonic_surf_pro, index);
01211                 IBERRCHK(g_err, "Trouble set the int data for geometrical face entity!");
01212                 f.index = index;
01213         }
01214         
01215         //get the adjacent vertices on a face
01216         vector<iBase_EntityHandle> tmp;
01217         g_err = igeom_instance->getEntAdj(f.gFaceHandle, iBase_VERTEX, tmp);
01218         IBERRCHK(g_err, "Trouble get the adjacent vertices on a geometrical face!");
01219         for (unsigned int i = 0; i < tmp.size(); i++){
01220                 int tmpint = -1;
01221                 g_err = igeom_instance->getIntData(tmp[i], harmonic_surf_pro, tmpint);
01222                 IBERRCHK(g_err, "Trouble get the int data for a vertex!");
01223                 f.connect.push_back(&vertices[tmpint]);
01224         }
01225         //get the adjacent edges on a face
01226         tmp.clear();
01227         g_err = igeom_instance->getEntAdj(f.gFaceHandle, iBase_EDGE, tmp);
01228         IBERRCHK(g_err, "Trouble get the adjacent edges on a geometrical face!");
01229         for (unsigned int i = 0; i < tmp.size(); i++){
01230                 int tmpint = -1;
01231                 g_err = igeom_instance->getIntData(tmp[i], harmonic_surf_pro, tmpint);
01232                 IBERRCHK(g_err, "Trouble get the int data for a vertex!");
01233                 f.connEdges.push_back(&edges[tmpint]);
01234         }
01235         
01236 }
01237 
01238 //post process the geometric loops on the source/target surface
01239 //problem not to be solved: matching the loops between source and target surfaces should be done!!!!! 
01240 void SurfProHarmonicMap::postProcessGeomLoops(Face& surf)
01241 {
01242         //make sure the outmost boundary loop is at [0];
01243         int outmost_index = -1;
01244         double volume = -1.0;
01245         if (surf.edgeloops.size()<=1)
01246                 return;
01247         for (unsigned int j = 0; j < surf.edgeloops.size(); j++){
01248                 std::vector<iBase_EntityHandle> entities;
01249                 double mincorner[3] = {1.0e10,1.0e10,1.0e10}, maxcorner[3] = {-1.0e10,-1.0e10,-1.0e10};
01250                 for (unsigned int i = 0; i < surf.edgeloops[j].size(); i++){
01251                         double t_min[3], t_max[3];
01252                         entities.clear();
01253                         entities.push_back(edges[surf.edgeloops[j][i]].gEdgeHandle);
01254                         g_err = igeom_instance->getArrBoundBox(&entities[0], entities.size(), iBase_StorageOrder_MIN, &t_min[0], &t_max[0]);
01255                         IBERRCHK(g_err, "Trouble get the bound box for an array of entities");
01256                         if (fabs(mincorner[0])==fabs(mincorner[1])&&fabs(mincorner[1])==fabs(mincorner[2])&&fabs(maxcorner[0])==fabs(maxcorner[1])&&fabs(maxcorner[1])==fabs(mincorner[2])){
01257                                 for (int m = 0; m < 3; m++){
01258                                         mincorner[m] = t_min[m];
01259                                         maxcorner[m] = t_max[m];
01260                                 }
01261                         }
01262                         else{
01263                                 if (t_min[0] < mincorner[0])
01264                                         mincorner[0] = t_min[0];
01265                                 if (t_min[1] < mincorner[1])                            
01266                                         mincorner[1] = t_min[1];
01267                                 if (t_min[2] < mincorner[2])                                    
01268                                         mincorner[2] = t_min[2];
01269                                 if (t_max[0] > maxcorner[0])
01270                                         maxcorner[0] = t_max[0];
01271                                 if (t_max[1] > maxcorner[1])
01272                                         maxcorner[1] = t_max[1];
01273                                 if (t_max[2] > maxcorner[2])
01274                                         maxcorner[2] = t_max[2];                                
01275                         }                       
01276                 }
01277                 //g_err = igeom_instance->getArrBoundBox(&entities[0], entities.size(), iBase_StorageOrder_MIN, &mincorner[0], &maxcorner[0]);
01278                 //IBERRCHK(g_err, "Trouble get the bound box for an array of entities");
01279                 double len_x = fabs(maxcorner[0]-mincorner[0]), len_y = fabs(maxcorner[1]-mincorner[1]), len_z = fabs(maxcorner[2]-mincorner[2]);
01280                 if (fabs(maxcorner[0]-mincorner[0]) < 1.0e-5)
01281                         len_x = 1.0;
01282                 if (fabs(maxcorner[1]-mincorner[1]) < 1.0e-5)
01283                         len_y = 1.0;
01284                 if (fabs(maxcorner[2]-mincorner[2]) < 1.0e-5)
01285                         len_z = 1.0;
01286                 double tmp_volume = len_x*len_y*len_z;
01287                 if ( tmp_volume > volume){
01288                         volume = tmp_volume;
01289                         outmost_index = j;
01290                 }
01291         }
01292     if (outmost_index == 0)
01293                 return;
01294         //exchange loop[0] and loop[outmost_index]
01295         std::vector<int> tmp_loop1, tmp_loop2;
01296         for (unsigned int j = 0; j < surf.vertexloops[0].size(); j++){
01297                 tmp_loop1.push_back(surf.vertexloops[0][j]);
01298                 tmp_loop2.push_back(surf.edgeloops[0][j]);
01299         }
01300         surf.vertexloops[0].clear();
01301         surf.edgeloops[0].clear();
01302         for (unsigned int j = 0; j < surf.vertexloops[outmost_index].size(); j++){
01303                 surf.vertexloops[0].push_back(surf.vertexloops[outmost_index][j]);
01304                 surf.edgeloops[0].push_back(surf.edgeloops[outmost_index][j]);
01305         }
01306         surf.vertexloops[outmost_index].clear();
01307         surf.edgeloops[outmost_index].clear();
01308 
01309         surf.vertexloops[outmost_index].insert(surf.vertexloops[outmost_index].begin(), tmp_loop1.begin(), tmp_loop1.end());
01310         surf.edgeloops[outmost_index].insert(surf.edgeloops[outmost_index].begin(), tmp_loop2.begin(), tmp_loop2.end());
01311         
01312 }
01313 
01314 //get the geometric loops on target surfaces
01315 void SurfProHarmonicMap::GetGeomLoops(Face surf, vector<vector<int> > &loops_vertex, vector<vector<int> > &loops_edge)
01316 {
01317         //collect edges' information
01318         set<int> edge_index;
01319         for (unsigned int i = 0; i < surf.connEdges.size(); i++)
01320                 edge_index.insert(surf.connEdges[i]->index);
01321 
01322         int     CurrentNum = 0;
01323         int TotalNum = surf.connect.size();
01324         
01325         while(CurrentNum < TotalNum){
01326                 int start_vertex = 0, next_vertex = -1, begin_vertex = 0;
01327                 int start_edge = 0, next_edge = -1;
01328                 
01329                 //get the starting edge
01330                 start_vertex = edges[*(edge_index.begin())].connect[0]->index;
01331                 start_edge = *(edge_index.begin());
01332 
01333                 //get the edge sense and make sure that the nodes are oriented correctly
01334                 int edge_sense = 0;
01335                 int first_vertex = start_vertex, second_vertex;
01336                 if (edges[*(edge_index.begin())].connect.size()==2){
01337                         second_vertex = edges[*(edge_index.begin())].connect[1]->index;
01338                         g_err = igeom_instance->getEgVtxSense(edges[start_edge].gEdgeHandle, vertices[first_vertex].gVertexHandle, vertices[second_vertex].gVertexHandle, edge_sense);
01339                         IBERRCHK(g_err, "Trouble get the edge sense.");
01340                         int face_sense = 0;
01341                         g_err = igeom_instance->getEgFcSense(edges[start_edge].gEdgeHandle,  surf.gFaceHandle, face_sense);
01342                         IBERRCHK(g_err, "Trouble get the face sense.");
01343                         if (face_sense*edge_sense < 0){
01344                                 start_vertex = edges[*(edge_index.begin())].connect[1]->index;
01345                         }
01346                 }
01347                 //initialization
01348                 loops_vertex.resize(loops_vertex.size()+1);
01349                 loops_edge.resize(loops_edge.size()+1);
01350                 begin_vertex = start_vertex;
01351                 //find the geometric loop               
01352                 while(next_vertex != begin_vertex){
01353                         loops_vertex[loops_vertex.size()-1].push_back(start_vertex);
01354                         loops_edge[loops_edge.size()-1].push_back(start_edge);
01355                         edge_index.erase(start_edge);
01356                         if (edges[start_edge].connect.size()==1)//circle case
01357                                 break;
01358                         if (start_vertex == edges[start_edge].connect[0]->index)
01359                                 next_vertex = edges[start_edge].connect[1]->index;
01360                         else
01361                                 next_vertex = edges[start_edge].connect[0]->index;
01362                         //find the adjacent geometric edge
01363                         vector<iBase_EntityHandle> adj_edges;
01364                         g_err = igeom_instance->getEntAdj(vertices[next_vertex].gVertexHandle, iBase_EDGE, adj_edges);
01365                         IBERRCHK(g_err, "Trouble get the adjacent edges around a geometric vertex.");
01366                         //remove unnecessary geometric edges
01367                         assert(adj_edges.size()>=1);
01368                         int intdata = -1;
01369                         int next_edge_index = -1;
01370                         for (int i = adj_edges.size()-1; i >= 0; i--){
01371                                 g_err = igeom_instance->getIntData(adj_edges[i], harmonic_surf_pro, intdata);
01372                                 if (edge_index.find(intdata)==edge_index.end())//remove that geometric edge
01373                                         continue;
01374                                 else{
01375                                         next_edge_index = intdata;
01376                                         break;
01377                                 }
01378                         }
01379                         if (next_edge_index == -1)
01380                                 break;
01381                         next_edge = intdata;
01382                         start_vertex = next_vertex;
01383                         start_edge = next_edge;
01384                 }
01385                 CurrentNum += loops_vertex[loops_vertex.size()-1].size();
01386         }
01387 }
01388 
01389 
01390 SurfProHarmonicMap::~SurfProHarmonicMap()
01391 {
01392 }
01393 
01394 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines