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