Developer Documentation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
Geometry.inl
1 /*
2 Copyright (c) 2006, Michael Kazhdan and Matthew Bolitho
3 All rights reserved.
4 
5 Redistribution and use in source and binary forms, with or without modification,
6 are permitted provided that the following conditions are met:
7 
8 Redistributions of source code must retain the above copyright notice, this list of
9 conditions and the following disclaimer. Redistributions in binary form must reproduce
10 the above copyright notice, this list of conditions and the following disclaimer
11 in the documentation and/or other materials provided with the distribution.
12 
13 Neither the name of the Johns Hopkins University nor the names of its contributors
14 may be used to endorse or promote products derived from this software without specific
15 prior written permission.
16 
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
18 EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES
19 OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
20 SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
21 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
22 TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
23 BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
24 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
25 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
26 DAMAGE.
27 */
28 
29 #include <stdio.h>
30 
31 template<class Real>
32 Real Random(void){return Real(rand())/RAND_MAX;}
33 
34 template<class Real>
35 Point3D<Real> RandomBallPoint(void){
36  Point3D<Real> p;
37  while(1){
38  p.coords[0]=Real(1.0-2.0*Random<Real>());
39  p.coords[1]=Real(1.0-2.0*Random<Real>());
40  p.coords[2]=Real(1.0-2.0*Random<Real>());
41  double l=SquareLength(p);
42  if(l<=1){return p;}
43  }
44 }
45 template<class Real>
46 Point3D<Real> RandomSpherePoint(void){
47  Point3D<Real> p=RandomBallPoint<Real>();
48  Real l=Real(Length(p));
49  p.coords[0]/=l;
50  p.coords[1]/=l;
51  p.coords[2]/=l;
52  return p;
53 }
54 
55 template<class Real>
56 double SquareLength(const Point3D<Real>& p){return p.coords[0]*p.coords[0]+p.coords[1]*p.coords[1]+p.coords[2]*p.coords[2];}
57 
58 template<class Real>
59 double Length(const Point3D<Real>& p){return sqrt(SquareLength(p));}
60 
61 template<class Real>
62 double SquareDistance(const Point3D<Real>& p1,const Point3D<Real>& p2){
63  return (p1.coords[0]-p2.coords[0])*(p1.coords[0]-p2.coords[0])+(p1.coords[1]-p2.coords[1])*(p1.coords[1]-p2.coords[1])+(p1.coords[2]-p2.coords[2])*(p1.coords[2]-p2.coords[2]);
64 }
65 
66 template<class Real>
67 double Distance(const Point3D<Real>& p1,const Point3D<Real>& p2){return sqrt(SquareDistance(p1,p2));}
68 
69 template <class Real>
70 void CrossProduct(const Point3D<Real>& p1,const Point3D<Real>& p2,Point3D<Real>& p){
71  p.coords[0]= p1.coords[1]*p2.coords[2]-p1.coords[2]*p2.coords[1];
72  p.coords[1]=-p1.coords[0]*p2.coords[2]+p1.coords[2]*p2.coords[0];
73  p.coords[2]= p1.coords[0]*p2.coords[1]-p1.coords[1]*p2.coords[0];
74 }
75 template<class Real>
76 void EdgeCollapse(const Real& edgeRatio,std::vector<TriangleIndex>& triangles,std::vector< Point3D<Real> >& positions,std::vector< Point3D<Real> >* normals){
77  int i,j,*remapTable,*pointCount,idx[3];
78  Point3D<Real> p[3],q[2],c;
79  double d[3],a;
80  double Ratio=12.0/sqrt(3.0); // (Sum of Squares Length / Area) for and equilateral triangle
81 
82  remapTable=new int[positions.size()];
83  pointCount=new int[positions.size()];
84  for(i=0;i<int(positions.size());i++){
85  remapTable[i]=i;
86  pointCount[i]=1;
87  }
88  for(i=int(triangles.size()-1);i>=0;i--){
89  for(j=0;j<3;j++){
90  idx[j]=triangles[i].idx[j];
91  while(remapTable[idx[j]]<idx[j]){idx[j]=remapTable[idx[j]];}
92  }
93  if(idx[0]==idx[1] || idx[0]==idx[2] || idx[1]==idx[2]){
94  triangles[i]=triangles[triangles.size()-1];
95  triangles.pop_back();
96  continue;
97  }
98  for(j=0;j<3;j++){
99  p[j].coords[0]=positions[idx[j]].coords[0]/pointCount[idx[j]];
100  p[j].coords[1]=positions[idx[j]].coords[1]/pointCount[idx[j]];
101  p[j].coords[2]=positions[idx[j]].coords[2]/pointCount[idx[j]];
102  }
103  for(j=0;j<3;j++){
104  q[0].coords[j]=p[1].coords[j]-p[0].coords[j];
105  q[1].coords[j]=p[2].coords[j]-p[0].coords[j];
106  d[j]=SquareDistance(p[j],p[(j+1)%3]);
107  }
108  CrossProduct(q[0],q[1],c);
109  a=Length(c)/2;
110 
111  if((d[0]+d[1]+d[2])*edgeRatio > a*Ratio){
112  // Find the smallest edge
113  j=0;
114  if(d[1]<d[j]){j=1;}
115  if(d[2]<d[j]){j=2;}
116 
117  int idx1,idx2;
118  if(idx[j]<idx[(j+1)%3]){
119  idx1=idx[j];
120  idx2=idx[(j+1)%3];
121  }
122  else{
123  idx2=idx[j];
124  idx1=idx[(j+1)%3];
125  }
126  positions[idx1].coords[0]+=positions[idx2].coords[0];
127  positions[idx1].coords[1]+=positions[idx2].coords[1];
128  positions[idx1].coords[2]+=positions[idx2].coords[2];
129  if(normals){
130  (*normals)[idx1].coords[0]+=(*normals)[idx2].coords[0];
131  (*normals)[idx1].coords[1]+=(*normals)[idx2].coords[1];
132  (*normals)[idx1].coords[2]+=(*normals)[idx2].coords[2];
133  }
134  pointCount[idx1]+=pointCount[idx2];
135  remapTable[idx2]=idx1;
136  triangles[i]=triangles[triangles.size()-1];
137  triangles.pop_back();
138  }
139  }
140  int pCount=0;
141  for(i=0;i<int(positions.size());i++){
142  for(j=0;j<3;j++){positions[i].coords[j]/=pointCount[i];}
143  if(normals){
144  Real l=Real(Length((*normals)[i]));
145  for(j=0;j<3;j++){(*normals)[i].coords[j]/=l;}
146  }
147  if(remapTable[i]==i){ // If vertex i is being used
148  positions[pCount]=positions[i];
149  if(normals){(*normals)[pCount]=(*normals)[i];}
150  pointCount[i]=pCount;
151  pCount++;
152  }
153  }
154  positions.resize(pCount);
155  for(i=int(triangles.size()-1);i>=0;i--){
156  for(j=0;j<3;j++){
157  idx[j]=triangles[i].idx[j];
158  while(remapTable[idx[j]]<idx[j]){idx[j]=remapTable[idx[j]];}
159  triangles[i].idx[j]=pointCount[idx[j]];
160  }
161  if(idx[0]==idx[1] || idx[0]==idx[2] || idx[1]==idx[2]){
162  triangles[i]=triangles[triangles.size()-1];
163  triangles.pop_back();
164  }
165  }
166 
167  delete[] pointCount;
168  delete[] remapTable;
169 }
170 template<class Real>
171 void TriangleCollapse(const Real& edgeRatio,std::vector<TriangleIndex>& triangles,std::vector< Point3D<Real> >& positions,std::vector< Point3D<Real> >* normals){
172  int i,j,*remapTable,*pointCount,idx[3];
173  Point3D<Real> p[3],q[2],c;
174  double d[3],a;
175  double Ratio=12.0/sqrt(3.0); // (Sum of Squares Length / Area) for and equilateral triangle
176 
177  remapTable=new int[positions.size()];
178  pointCount=new int[positions.size()];
179  for(i=0;i<int(positions.size());i++){
180  remapTable[i]=i;
181  pointCount[i]=1;
182  }
183  for(i=int(triangles.size()-1);i>=0;i--){
184  for(j=0;j<3;j++){
185  idx[j]=triangles[i].idx[j];
186  while(remapTable[idx[j]]<idx[j]){idx[j]=remapTable[idx[j]];}
187  }
188  if(idx[0]==idx[1] || idx[0]==idx[2] || idx[1]==idx[2]){
189  triangles[i]=triangles[triangles.size()-1];
190  triangles.pop_back();
191  continue;
192  }
193  for(j=0;j<3;j++){
194  p[j].coords[0]=positions[idx[j]].coords[0]/pointCount[idx[j]];
195  p[j].coords[1]=positions[idx[j]].coords[1]/pointCount[idx[j]];
196  p[j].coords[2]=positions[idx[j]].coords[2]/pointCount[idx[j]];
197  }
198  for(j=0;j<3;j++){
199  q[0].coords[j]=p[1].coords[j]-p[0].coords[j];
200  q[1].coords[j]=p[2].coords[j]-p[0].coords[j];
201  d[j]=SquareDistance(p[j],p[(j+1)%3]);
202  }
203  CrossProduct(q[0],q[1],c);
204  a=Length(c)/2;
205 
206  if((d[0]+d[1]+d[2])*edgeRatio > a*Ratio){
207  // Find the smallest edge
208  j=0;
209  if(d[1]<d[j]){j=1;}
210  if(d[2]<d[j]){j=2;}
211 
212  int idx1,idx2,idx3;
213  if(idx[0]<idx[1]){
214  if(idx[0]<idx[2]){
215  idx1=idx[0];
216  idx2=idx[2];
217  idx3=idx[1];
218  }
219  else{
220  idx1=idx[2];
221  idx2=idx[0];
222  idx3=idx[1];
223  }
224  }
225  else{
226  if(idx[1]<idx[2]){
227  idx1=idx[1];
228  idx2=idx[2];
229  idx3=idx[0];
230  }
231  else{
232  idx1=idx[2];
233  idx2=idx[1];
234  idx3=idx[0];
235  }
236  }
237  positions[idx1].coords[0]+=positions[idx2].coords[0]+positions[idx3].coords[0];
238  positions[idx1].coords[1]+=positions[idx2].coords[1]+positions[idx3].coords[1];
239  positions[idx1].coords[2]+=positions[idx2].coords[2]+positions[idx3].coords[2];
240  if(normals){
241  (*normals)[idx1].coords[0]+=(*normals)[idx2].coords[0]+(*normals)[idx3].coords[0];
242  (*normals)[idx1].coords[1]+=(*normals)[idx2].coords[1]+(*normals)[idx3].coords[1];
243  (*normals)[idx1].coords[2]+=(*normals)[idx2].coords[2]+(*normals)[idx3].coords[2];
244  }
245  pointCount[idx1]+=pointCount[idx2]+pointCount[idx3];
246  remapTable[idx2]=idx1;
247  remapTable[idx3]=idx1;
248  triangles[i]=triangles[triangles.size()-1];
249  triangles.pop_back();
250  }
251  }
252  int pCount=0;
253  for(i=0;i<int(positions.size());i++){
254  for(j=0;j<3;j++){positions[i].coords[j]/=pointCount[i];}
255  if(normals){
256  Real l=Real(Length((*normals)[i]));
257  for(j=0;j<3;j++){(*normals)[i].coords[j]/=l;}
258  }
259  if(remapTable[i]==i){ // If vertex i is being used
260  positions[pCount]=positions[i];
261  if(normals){(*normals)[pCount]=(*normals)[i];}
262  pointCount[i]=pCount;
263  pCount++;
264  }
265  }
266  positions.resize(pCount);
267  for(i=int(triangles.size()-1);i>=0;i--){
268  for(j=0;j<3;j++){
269  idx[j]=triangles[i].idx[j];
270  while(remapTable[idx[j]]<idx[j]){idx[j]=remapTable[idx[j]];}
271  triangles[i].idx[j]=pointCount[idx[j]];
272  }
273  if(idx[0]==idx[1] || idx[0]==idx[2] || idx[1]==idx[2]){
274  triangles[i]=triangles[triangles.size()-1];
275  triangles.pop_back();
276  }
277  }
278  delete[] pointCount;
279  delete[] remapTable;
280 }
281 
283 // Triangulation //
285 template<class Real>
286 long long Triangulation<Real>::EdgeIndex( int p1 , int p2 )
287 {
288  if(p1>p2) {return ((long long)(p1)<<32) | ((long long)(p2));}
289  else {return ((long long)(p2)<<32) | ((long long)(p1));}
290 }
291 
292 template<class Real>
293 int Triangulation<Real>::factor(int tIndex,int& p1,int& p2,int & p3){
294  if(triangles[tIndex].eIndex[0]<0 || triangles[tIndex].eIndex[1]<0 || triangles[tIndex].eIndex[2]<0){return 0;}
295  if(edges[triangles[tIndex].eIndex[0]].tIndex[0]==tIndex){p1=edges[triangles[tIndex].eIndex[0]].pIndex[0];}
296  else {p1=edges[triangles[tIndex].eIndex[0]].pIndex[1];}
297  if(edges[triangles[tIndex].eIndex[1]].tIndex[0]==tIndex){p2=edges[triangles[tIndex].eIndex[1]].pIndex[0];}
298  else {p2=edges[triangles[tIndex].eIndex[1]].pIndex[1];}
299  if(edges[triangles[tIndex].eIndex[2]].tIndex[0]==tIndex){p3=edges[triangles[tIndex].eIndex[2]].pIndex[0];}
300  else {p3=edges[triangles[tIndex].eIndex[2]].pIndex[1];}
301  return 1;
302 }
303 template<class Real>
304 double Triangulation<Real>::area(int p1,int p2,int p3){
305  Point3D<Real> q1,q2,q;
306  for(int i=0;i<3;i++){
307  q1.coords[i]=points[p2].coords[i]-points[p1].coords[i];
308  q2.coords[i]=points[p3].coords[i]-points[p1].coords[i];
309  }
310  CrossProduct(q1,q2,q);
311  return Length(q);
312 }
313 template<class Real>
314 double Triangulation<Real>::area(int tIndex){
315  int p1,p2,p3;
316  factor(tIndex,p1,p2,p3);
317  return area(p1,p2,p3);
318 }
319 template<class Real>
320 double Triangulation<Real>::area(void){
321  double a=0;
322  for(int i=0;i<int(triangles.size());i++){a+=area(i);}
323  return a;
324 }
325 template<class Real>
326 int Triangulation<Real>::addTriangle(int p1,int p2,int p3){
327  hash_map<long long,int>::iterator iter;
328  int tIdx,eIdx,p[3];
329  p[0]=p1;
330  p[1]=p2;
331  p[2]=p3;
332  triangles.push_back(TriangulationTriangle());
333  tIdx=int(triangles.size())-1;
334 
335  for(int i=0;i<3;i++)
336  {
337  long long e = EdgeIndex(p[i],p[(i+1)%3]);
338  iter=edgeMap.find(e);
339  if(iter==edgeMap.end())
340  {
341  TriangulationEdge edge;
342  edge.pIndex[0]=p[i];
343  edge.pIndex[1]=p[(i+1)%3];
344  edges.push_back(edge);
345  eIdx=int(edges.size())-1;
346  edgeMap[e]=eIdx;
347  edges[eIdx].tIndex[0]=tIdx;
348  }
349  else{
350  eIdx=edgeMap[e];
351  if(edges[eIdx].pIndex[0]==p[i]){
352  if(edges[eIdx].tIndex[0]<0){edges[eIdx].tIndex[0]=tIdx;}
353  else{printf("Edge Triangle in use 1\n");return 0;}
354  }
355  else{
356  if(edges[eIdx].tIndex[1]<0){edges[eIdx].tIndex[1]=tIdx;}
357  else{printf("Edge Triangle in use 2\n");return 0;}
358  }
359 
360  }
361  triangles[tIdx].eIndex[i]=eIdx;
362  }
363  return tIdx;
364 }
365 template<class Real>
366 int Triangulation<Real>::flipMinimize(int eIndex){
367  double oldArea,newArea;
368  int oldP[3],oldQ[3],newP[3],newQ[3];
369  TriangulationEdge newEdge;
370 
371  if(edges[eIndex].tIndex[0]<0 || edges[eIndex].tIndex[1]<0){return 0;}
372 
373  if(!factor(edges[eIndex].tIndex[0],oldP[0],oldP[1],oldP[2])){return 0;}
374  if(!factor(edges[eIndex].tIndex[1],oldQ[0],oldQ[1],oldQ[2])){return 0;}
375 
376  oldArea=area(oldP[0],oldP[1],oldP[2])+area(oldQ[0],oldQ[1],oldQ[2]);
377  int idxP,idxQ;
378  for(idxP=0;idxP<3;idxP++){
379  int i;
380  for(i=0;i<3;i++){if(oldP[idxP]==oldQ[i]){break;}}
381  if(i==3){break;}
382  }
383  for(idxQ=0;idxQ<3;idxQ++){
384  int i;
385  for(i=0;i<3;i++){if(oldP[i]==oldQ[idxQ]){break;}}
386  if(i==3){break;}
387  }
388  if(idxP==3 || idxQ==3){return 0;}
389  newP[0]=oldP[idxP];
390  newP[1]=oldP[(idxP+1)%3];
391  newP[2]=oldQ[idxQ];
392  newQ[0]=oldQ[idxQ];
393  newQ[1]=oldP[(idxP+2)%3];
394  newQ[2]=oldP[idxP];
395 
396  newArea=area(newP[0],newP[1],newP[2])+area(newQ[0],newQ[1],newQ[2]);
397  if(oldArea<=newArea){return 0;}
398 
399  // Remove the entry in the hash_table for the old edge
400  edgeMap.erase(EdgeIndex(edges[eIndex].pIndex[0],edges[eIndex].pIndex[1]));
401  // Set the new edge so that the zero-side is newQ
402  edges[eIndex].pIndex[0]=newP[0];
403  edges[eIndex].pIndex[1]=newQ[0];
404  // Insert the entry into the hash_table for the new edge
405  edgeMap[EdgeIndex(newP[0],newQ[0])]=eIndex;
406  // Update the triangle information
407  for(int i=0;i<3;i++){
408  int idx;
409  idx=edgeMap[EdgeIndex(newQ[i],newQ[(i+1)%3])];
410  triangles[edges[eIndex].tIndex[0]].eIndex[i]=idx;
411  if(idx!=eIndex){
412  if(edges[idx].tIndex[0]==edges[eIndex].tIndex[1]){edges[idx].tIndex[0]=edges[eIndex].tIndex[0];}
413  if(edges[idx].tIndex[1]==edges[eIndex].tIndex[1]){edges[idx].tIndex[1]=edges[eIndex].tIndex[0];}
414  }
415 
416  idx=edgeMap[EdgeIndex(newP[i],newP[(i+1)%3])];
417  triangles[edges[eIndex].tIndex[1]].eIndex[i]=idx;
418  if(idx!=eIndex){
419  if(edges[idx].tIndex[0]==edges[eIndex].tIndex[0]){edges[idx].tIndex[0]=edges[eIndex].tIndex[1];}
420  if(edges[idx].tIndex[1]==edges[eIndex].tIndex[0]){edges[idx].tIndex[1]=edges[eIndex].tIndex[1];}
421  }
422  }
423  return 1;
424 }