|
MeshKit
1.0
|
|
MeshKit
1.0
|
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 }
1.7.6.1 
