32 Real Random(
void){
return Real(rand())/RAND_MAX;}
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);
48 Real l=Real(Length(p));
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];}
59 double Length(
const Point3D<Real>& p){
return sqrt(SquareLength(p));}
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]);
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];
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];
80 double Ratio=12.0/sqrt(3.0);
82 remapTable=
new int[positions.size()];
83 pointCount=
new int[positions.size()];
84 for(i=0;i<int(positions.size());i++){
88 for(i=
int(triangles.size()-1);i>=0;i--){
90 idx[j]=triangles[i].idx[j];
91 while(remapTable[idx[j]]<idx[j]){idx[j]=remapTable[idx[j]];}
93 if(idx[0]==idx[1] || idx[0]==idx[2] || idx[1]==idx[2]){
94 triangles[i]=triangles[triangles.size()-1];
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]];
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]);
108 CrossProduct(q[0],q[1],c);
111 if((d[0]+d[1]+d[2])*edgeRatio > a*Ratio){
118 if(idx[j]<idx[(j+1)%3]){
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];
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];
134 pointCount[idx1]+=pointCount[idx2];
135 remapTable[idx2]=idx1;
136 triangles[i]=triangles[triangles.size()-1];
137 triangles.pop_back();
141 for(i=0;i<int(positions.size());i++){
142 for(j=0;j<3;j++){positions[i].coords[j]/=pointCount[i];}
144 Real l=Real(Length((*normals)[i]));
145 for(j=0;j<3;j++){(*normals)[i].coords[j]/=l;}
147 if(remapTable[i]==i){
148 positions[pCount]=positions[i];
149 if(normals){(*normals)[pCount]=(*normals)[i];}
150 pointCount[i]=pCount;
154 positions.resize(pCount);
155 for(i=
int(triangles.size()-1);i>=0;i--){
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]];
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();
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];
175 double Ratio=12.0/sqrt(3.0);
177 remapTable=
new int[positions.size()];
178 pointCount=
new int[positions.size()];
179 for(i=0;i<int(positions.size());i++){
183 for(i=
int(triangles.size()-1);i>=0;i--){
185 idx[j]=triangles[i].idx[j];
186 while(remapTable[idx[j]]<idx[j]){idx[j]=remapTable[idx[j]];}
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();
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]];
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]);
203 CrossProduct(q[0],q[1],c);
206 if((d[0]+d[1]+d[2])*edgeRatio > a*Ratio){
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];
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];
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();
253 for(i=0;i<int(positions.size());i++){
254 for(j=0;j<3;j++){positions[i].coords[j]/=pointCount[i];}
256 Real l=Real(Length((*normals)[i]));
257 for(j=0;j<3;j++){(*normals)[i].coords[j]/=l;}
259 if(remapTable[i]==i){
260 positions[pCount]=positions[i];
261 if(normals){(*normals)[pCount]=(*normals)[i];}
262 pointCount[i]=pCount;
266 positions.resize(pCount);
267 for(i=
int(triangles.size()-1);i>=0;i--){
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]];
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();
288 if(p1>p2) {
return ((
long long)(p1)<<32) | ((
long long)(p2));}
289 else {
return ((
long long)(p2)<<32) | ((
long long)(p1));}
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];}
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];
310 CrossProduct(q1,q2,q);
316 factor(tIndex,p1,p2,p3);
317 return area(p1,p2,p3);
322 for(
int i=0;i<int(triangles.size());i++){a+=area(i);}
327 hash_map<long long,int>::iterator iter;
333 tIdx=int(triangles.size())-1;
337 long long e =
EdgeIndex(p[i],p[(i+1)%3]);
338 iter=edgeMap.find(e);
339 if(iter==edgeMap.end())
343 edge.pIndex[1]=p[(i+1)%3];
344 edges.push_back(edge);
345 eIdx=int(edges.size())-1;
347 edges[eIdx].tIndex[0]=tIdx;
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;}
356 if(edges[eIdx].tIndex[1]<0){edges[eIdx].tIndex[1]=tIdx;}
357 else{printf(
"Edge Triangle in use 2\n");
return 0;}
361 triangles[tIdx].eIndex[i]=eIdx;
367 double oldArea,newArea;
368 int oldP[3],oldQ[3],newP[3],newQ[3];
371 if(edges[eIndex].tIndex[0]<0 || edges[eIndex].tIndex[1]<0){
return 0;}
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;}
376 oldArea=area(oldP[0],oldP[1],oldP[2])+area(oldQ[0],oldQ[1],oldQ[2]);
378 for(idxP=0;idxP<3;idxP++){
380 for(i=0;i<3;i++){
if(oldP[idxP]==oldQ[i]){
break;}}
383 for(idxQ=0;idxQ<3;idxQ++){
385 for(i=0;i<3;i++){
if(oldP[i]==oldQ[idxQ]){
break;}}
388 if(idxP==3 || idxQ==3){
return 0;}
390 newP[1]=oldP[(idxP+1)%3];
393 newQ[1]=oldP[(idxP+2)%3];
396 newArea=area(newP[0],newP[1],newP[2])+area(newQ[0],newQ[1],newQ[2]);
397 if(oldArea<=newArea){
return 0;}
400 edgeMap.erase(
EdgeIndex(edges[eIndex].pIndex[0],edges[eIndex].pIndex[1]));
402 edges[eIndex].pIndex[0]=newP[0];
403 edges[eIndex].pIndex[1]=newQ[0];
405 edgeMap[
EdgeIndex(newP[0],newQ[0])]=eIndex;
407 for(
int i=0;i<3;i++){
409 idx=edgeMap[
EdgeIndex(newQ[i],newQ[(i+1)%3])];
410 triangles[edges[eIndex].tIndex[0]].eIndex[i]=idx;
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];}
416 idx=edgeMap[
EdgeIndex(newP[i],newP[(i+1)%3])];
417 triangles[edges[eIndex].tIndex[1]].eIndex[i]=idx;
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];}