Point Cloud Library (PCL) 1.13.1
Loading...
Searching...
No Matches
mat.hpp
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//////////////////////////////
31
32namespace pcl
33{
34 namespace poisson
35 {
36
37 template <class Real>
39 {
40 bestTriangulation=NULL;
41 midPoint=NULL;
42 }
43 template <class Real>
45 {
46 delete[] bestTriangulation;
47 bestTriangulation=NULL;
48 delete[] midPoint;
49 midPoint=NULL;
50 }
51 template <class Real>
52 void MinimalAreaTriangulation<Real>::GetTriangulation(const std::vector<Point3D<Real> >& vertices,std::vector<TriangleIndex>& triangles)
53 {
54 if(vertices.size()==3)
55 {
56 triangles.resize(1);
57 triangles[0].idx[0]=0;
58 triangles[0].idx[1]=1;
59 triangles[0].idx[2]=2;
60 return;
61 }
62 else if(vertices.size()==4)
63 {
64 TriangleIndex tIndex[2][2];
65 Real area[2];
66
67 area[0]=area[1]=0;
68 triangles.resize(2);
69
70 tIndex[0][0].idx[0]=0;
71 tIndex[0][0].idx[1]=1;
72 tIndex[0][0].idx[2]=2;
73 tIndex[0][1].idx[0]=2;
74 tIndex[0][1].idx[1]=3;
75 tIndex[0][1].idx[2]=0;
76
77 tIndex[1][0].idx[0]=0;
78 tIndex[1][0].idx[1]=1;
79 tIndex[1][0].idx[2]=3;
80 tIndex[1][1].idx[0]=3;
81 tIndex[1][1].idx[1]=1;
82 tIndex[1][1].idx[2]=2;
83
84 Point3D<Real> n,p1,p2;
85 for(int i=0;i<2;i++)
86 for(int j=0;j<2;j++)
87 {
88 p1=vertices[tIndex[i][j].idx[1]]-vertices[tIndex[i][j].idx[0]];
89 p2=vertices[tIndex[i][j].idx[2]]-vertices[tIndex[i][j].idx[0]];
90 CrossProduct(p1,p2,n);
91 area[i] += Real( Length(n) );
92 }
93 if(area[0]>area[1])
94 {
95 triangles[0]=tIndex[1][0];
96 triangles[1]=tIndex[1][1];
97 }
98 else
99 {
100 triangles[0]=tIndex[0][0];
101 triangles[1]=tIndex[0][1];
102 }
103 return;
104 }
105 delete[] bestTriangulation;
106 delete[] midPoint;
107 bestTriangulation=NULL;
108 midPoint=NULL;
109 std::size_t eCount=vertices.size();
110 bestTriangulation=new Real[eCount*eCount];
111 midPoint=new int[eCount*eCount];
112 for(std::size_t i=0;i<eCount*eCount;i++)
113 bestTriangulation[i]=-1;
114 memset(midPoint,-1,sizeof(int)*eCount*eCount);
115 GetArea(0,1,vertices);
116 triangles.clear();
117 GetTriangulation(0,1,vertices,triangles);
118 }
119 template <class Real>
121 {
122 delete[] bestTriangulation;
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 template<class Real>
135 void MinimalAreaTriangulation<Real>::GetTriangulation(const std::size_t& i,const std::size_t& j,const std::vector<Point3D<Real> >& vertices,std::vector<TriangleIndex>& triangles)
136 {
137 TriangleIndex tIndex;
138 std::size_t eCount=vertices.size();
139 std::size_t ii=i;
140 if(i<j)
141 ii+=eCount;
142 if(j+1>=ii)
143 return;
144 ii=midPoint[i*eCount+j];
145 if(ii>=0)
146 {
147 tIndex.idx[0] = int( i );
148 tIndex.idx[1] = int( j );
149 tIndex.idx[2] = int( ii );
150 triangles.push_back(tIndex);
151 GetTriangulation(i,ii,vertices,triangles);
152 GetTriangulation(ii,j,vertices,triangles);
153 }
154 }
155
156 template<class Real>
157 Real MinimalAreaTriangulation<Real>::GetArea(const std::size_t& i,const std::size_t& j,const std::vector<Point3D<Real> >& vertices)
158 {
159 Real a=FLT_MAX,temp;
160 std::size_t eCount=vertices.size();
161 std::size_t idx=i*eCount+j;
162 std::size_t ii=i;
163 if(i<j)
164 ii+=eCount;
165 if(j+1>=ii)
166 {
167 bestTriangulation[idx]=0;
168 return 0;
169 }
170 if(midPoint[idx]!=-1)
171 return bestTriangulation[idx];
172 int mid=-1;
173 for(std::size_t r=j+1;r<ii;r++)
174 {
175 std::size_t rr=r%eCount;
176 std::size_t idx1=i*eCount+rr,idx2=rr*eCount+j;
177 Point3D<Real> p,p1,p2;
178 p1=vertices[i]-vertices[rr];
179 p2=vertices[j]-vertices[rr];
180 CrossProduct(p1,p2,p);
181 temp = Real( Length(p) );
182 if(bestTriangulation[idx1]>=0)
183 {
184 temp+=bestTriangulation[idx1];
185 if(temp>a)
186 continue;
187 if(bestTriangulation[idx2]>0)
188 temp+=bestTriangulation[idx2];
189 else
190 temp+=GetArea(rr,j,vertices);
191 }
192 else
193 {
194 if(bestTriangulation[idx2]>=0)
195 temp+=bestTriangulation[idx2];
196 else
197 temp+=GetArea(rr,j,vertices);
198 if(temp>a)
199 continue;
200 temp+=GetArea(i,rr,vertices);
201 }
202
203 if(temp<a)
204 {
205 a=temp;
206 mid=int(rr);
207 }
208 }
209 bestTriangulation[idx]=a;
210 midPoint[idx]=mid;
211
212 return a;
213 }
214
215
216 }
217}
void CrossProduct(const Point3D< Real > &p1, const Point3D< Real > &p2, Point3D< Real > &p)
Definition geometry.hpp:74
double Length(const Point3D< Real > &p)
Definition geometry.hpp:63