1/* 2Copyright (c) 2007, Michael Kazhdan 3All rights reserved. 4 5Redistribution and use in source and binary forms, with or without modification, 6are permitted provided that the following conditions are met: 7 8Redistributions of source code must retain the above copyright notice, this list of 9conditions and the following disclaimer. Redistributions in binary form must reproduce 10the above copyright notice, this list of conditions and the following disclaimer 11in the documentation and/or other materials provided with the distribution. 12 13Neither the name of the Johns Hopkins University nor the names of its contributors 14may be used to endorse or promote products derived from this software without specific 15prior written permission. 16 17THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY 18EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO THE IMPLIED WARRANTIES 19OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT 20SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 21INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 22TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 23BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 24CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 25ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH 26DAMAGE. 27*/ 28////////////////////////////// 29// MinimalAreaTriangulation // 30////////////////////////////// 31template <class Real> 32MinimalAreaTriangulation<Real>::MinimalAreaTriangulation(void) 33{ 34 bestTriangulation=NULL; 35 midPoint=NULL; 36} 37template <class Real> 38MinimalAreaTriangulation<Real>::~MinimalAreaTriangulation(void) 39{ 40 if(bestTriangulation) 41 delete[] bestTriangulation; 42 bestTriangulation=NULL; 43 if(midPoint) 44 delete[] midPoint; 45 midPoint=NULL; 46} 47template <class Real> 48void MinimalAreaTriangulation<Real>::GetTriangulation(const std::vector<Point3D<Real> >& vertices,std::vector<TriangleIndex>& triangles) 49{ 50 if(vertices.size()==3) 51 { 52 triangles.resize(1); 53 triangles[0].idx[0]=0; 54 triangles[0].idx[1]=1; 55 triangles[0].idx[2]=2; 56 return; 57 } 58 else if(vertices.size()==4) 59 { 60 TriangleIndex tIndex[2][2]; 61 Real area[2]; 62 63 area[0]=area[1]=0; 64 triangles.resize(2); 65 66 tIndex[0][0].idx[0]=0; 67 tIndex[0][0].idx[1]=1; 68 tIndex[0][0].idx[2]=2; 69 tIndex[0][1].idx[0]=2; 70 tIndex[0][1].idx[1]=3; 71 tIndex[0][1].idx[2]=0; 72 73 tIndex[1][0].idx[0]=0; 74 tIndex[1][0].idx[1]=1; 75 tIndex[1][0].idx[2]=3; 76 tIndex[1][1].idx[0]=3; 77 tIndex[1][1].idx[1]=1; 78 tIndex[1][1].idx[2]=2; 79 80 Point3D<Real> n,p1,p2; 81 for(int i=0;i<2;i++) 82 for(int j=0;j<2;j++) 83 { 84 p1=vertices[tIndex[i][j].idx[1]]-vertices[tIndex[i][j].idx[0]]; 85 p2=vertices[tIndex[i][j].idx[2]]-vertices[tIndex[i][j].idx[0]]; 86 CrossProduct(p1,p2,n); 87 area[i] += Real( Length(n) ); 88 } 89 if(area[0]>area[1]) 90 { 91 triangles[0]=tIndex[1][0]; 92 triangles[1]=tIndex[1][1]; 93 } 94 else 95 { 96 triangles[0]=tIndex[0][0]; 97 triangles[1]=tIndex[0][1]; 98 } 99 return; 100 } 101 if(bestTriangulation) 102 delete[] bestTriangulation; 103 if(midPoint) 104 delete[] midPoint; 105 bestTriangulation=NULL; 106 midPoint=NULL; 107 size_t eCount=vertices.size(); 108 bestTriangulation=new Real[eCount*eCount]; 109 midPoint=new int[eCount*eCount]; 110 for(size_t i=0;i<eCount*eCount;i++) 111 bestTriangulation[i]=-1; 112 memset(midPoint,-1,sizeof(int)*eCount*eCount); 113 GetArea(0,1,vertices); 114 triangles.clear(); 115 GetTriangulation(0,1,vertices,triangles); 116} 117template <class Real> 118Real MinimalAreaTriangulation<Real>::GetArea(const std::vector<Point3D<Real> >& vertices) 119{ 120 if(bestTriangulation) 121 delete[] bestTriangulation; 122 if(midPoint) 123 delete[] midPoint; 124 bestTriangulation=NULL; 125 midPoint=NULL; 126 int eCount=vertices.size(); 127 bestTriangulation=new double[eCount*eCount]; 128 midPoint=new int[eCount*eCount]; 129 for(int i=0;i<eCount*eCount;i++) 130 bestTriangulation[i]=-1; 131 memset(midPoint,-1,sizeof(int)*eCount*eCount); 132 return GetArea(0,1,vertices); 133} 134 135template<class Real> 136void MinimalAreaTriangulation<Real>::GetTriangulation( 137 const size_t& i, 138 const size_t& j, 139 const std::vector<Point3D<Real> >& vertices, 140 std::vector<TriangleIndex>& triangles 141) { 142 TriangleIndex tIndex; 143 size_t eCount=vertices.size(); 144 int ii=int(i); // [Bruno Levy] ii is an int, note a size_t, 145 // it is set 5 lines later to midPoint[] that is an int 146 // (else test if(ii>=0) is always true !!) 147 if(i<j) 148 ii+=int(eCount); 149 if(j+1>=size_t(ii)) 150 return; 151 ii=midPoint[i*eCount+j]; 152 if(ii>=0) 153 { 154 tIndex.idx[0] = int( i ); 155 tIndex.idx[1] = int( j ); 156 tIndex.idx[2] = int( ii ); 157 triangles.push_back(tIndex); 158 GetTriangulation(i,ii,vertices,triangles); 159 GetTriangulation(ii,j,vertices,triangles); 160 } 161} 162 163template<class Real> 164Real MinimalAreaTriangulation<Real>::GetArea(const size_t& i,const size_t& j,const std::vector<Point3D<Real> >& vertices) 165{ 166 Real a=FLT_MAX,temp; 167 size_t eCount=vertices.size(); 168 size_t idx=i*eCount+j; 169 size_t ii=i; 170 if(i<j) 171 ii+=eCount; 172 if(j+1>=ii) 173 { 174 bestTriangulation[idx]=0; 175 return 0; 176 } 177 if(midPoint[idx]!=-1) 178 return bestTriangulation[idx]; 179 int mid=-1; 180 for(size_t r=j+1;r<ii;r++) 181 { 182 size_t rr=r%eCount; 183 size_t idx1=i*eCount+rr,idx2=rr*eCount+j; 184 Point3D<Real> p,p1,p2; 185 p1=vertices[i]-vertices[rr]; 186 p2=vertices[j]-vertices[rr]; 187 CrossProduct(p1,p2,p); 188 temp = Real( Length(p) ); 189 if(bestTriangulation[idx1]>=0) 190 { 191 temp+=bestTriangulation[idx1]; 192 if(temp>a) 193 continue; 194 if(bestTriangulation[idx2]>0) 195 temp+=bestTriangulation[idx2]; 196 else 197 temp+=GetArea(rr,j,vertices); 198 } 199 else 200 { 201 if(bestTriangulation[idx2]>=0) 202 temp+=bestTriangulation[idx2]; 203 else 204 temp+=GetArea(rr,j,vertices); 205 if(temp>a) 206 continue; 207 temp+=GetArea(i,rr,vertices); 208 } 209 210 if(temp<a) 211 { 212 a=temp; 213 mid=int(rr); 214 } 215 } 216 bestTriangulation[idx]=a; 217 midPoint[idx]=mid; 218 219 return a; 220} 221