Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <queue>
34 : #include <algorithm>
35 : #include "triangle_fill.h"
36 : #include "common/math/ugmath.h"
37 : #include "lib_grid/algorithms/geom_obj_util/geom_obj_util.h"
38 : #include "lib_grid/algorithms/polychain_util.h"
39 :
40 : using namespace std;
41 :
42 : namespace ug
43 : {
44 :
45 : /** Make sure that edgeNormalsOut has the same size as poly-chain.
46 : * The poly-chain is treated as a closed poly-chain where an edge
47 : * between the last and the first entry exists.*/
48 0 : void CalculatePolychainEdgeNormals(vector2* edgeNormalsOut, vector2* polyChain,
49 : size_t polyChainSize, bool bOuterNormals)
50 : {
51 0 : if(polyChainSize < 2)
52 0 : return;
53 :
54 : // perform a first guess
55 : // in order to check whether those normals are inner or outer normals,
56 : // we're building the sum of the normal of an edge with the following edge.
57 : number dot = 0;
58 :
59 : // the edge to which we will calculate the normal
60 : vector2 e;
61 : VecSubtract(e, polyChain[1], polyChain[0]);
62 0 : VecNormalize(e, e);
63 :
64 0 : for(size_t i = 0; i < polyChainSize; ++i){
65 0 : edgeNormalsOut[i].x() = e.y();
66 0 : edgeNormalsOut[i].y() = -e.x();
67 :
68 : // calculate the next edge (the normal in the next iteration
69 : // will be calculated for this edge)
70 0 : VecSubtract(e, polyChain[(i+2)%polyChainSize],
71 0 : polyChain[(i+1)%polyChainSize]);
72 0 : VecNormalize(e, e);
73 :
74 : // Build the dot-product with the last normal
75 0 : dot += VecDot(edgeNormalsOut[i], e);
76 : }
77 :
78 : // if the dot-product is smaller than 0, the normal is an outer normal
79 : // if it is bigger than 0, the normal is an inner normal.
80 0 : if((dot < 0 && (!bOuterNormals)) ||
81 0 : (dot > 0 && bOuterNormals))
82 : {
83 : // we have to invert the normals
84 0 : for(size_t i = 0; i < polyChainSize; ++i){
85 0 : edgeNormalsOut[i].x() = -edgeNormalsOut[i].x();
86 0 : edgeNormalsOut[i].y() = -edgeNormalsOut[i].y();
87 : }
88 : }
89 : }
90 :
91 : ////////////////////////////////////////////////////////////////////////
92 : struct ChainInfo{
93 0 : ChainInfo() {}
94 : ChainInfo(int vrtInd, int vrtIndIn, int vrtIndOut) :
95 0 : myVrt(vrtInd), inVrt(vrtIndIn), outVrt(vrtIndOut),
96 : isCandidate(false), associatedDistanceSq(0) {}
97 :
98 : bool operator == (const ChainInfo& ci){
99 0 : return myVrt == ci.myVrt && inVrt == ci.inVrt && outVrt == ci.outVrt
100 0 : && isCandidate == ci.isCandidate
101 0 : && associatedDistanceSq == ci.associatedDistanceSq;
102 : }
103 :
104 : bool operator < (const ChainInfo& ci) const {return associatedDistanceSq >
105 : ci.associatedDistanceSq;}
106 :
107 : size_t myVrt;
108 : size_t inVrt;
109 : size_t outVrt;
110 : bool isCandidate;
111 : number associatedDistanceSq;
112 : };
113 :
114 0 : void UpdateChainInfo(ChainInfo& ci, vector2* polyChain, vector2* edgeNormals)
115 : {
116 : vector2 outEdge;
117 0 : VecSubtract(outEdge, polyChain[ci.outVrt], polyChain[ci.myVrt]);
118 0 : ci.isCandidate = VecDot(edgeNormals[ci.inVrt], outEdge) < -SMALL; //ci.inVrt == inEdgeIndex.
119 :
120 : // calculate relative length
121 0 : number l = VecDistance(polyChain[ci.inVrt], polyChain[ci.outVrt]);
122 0 : number t1 = VecDistance(polyChain[ci.inVrt], polyChain[ci.myVrt]);
123 0 : number t2 = VecDistance(polyChain[ci.outVrt], polyChain[ci.myVrt]);
124 : number relLen = 1;
125 0 : if(t1 + t2 > SMALL)
126 0 : relLen = l / (t1 + t2);
127 :
128 : // if the relative length is small we'll calculate the squared length as the indiceator
129 0 : if(relLen < 0.9)
130 : {
131 0 : ci.associatedDistanceSq = VecDistanceSq(polyChain[ci.inVrt],
132 : polyChain[ci.outVrt]);
133 : }
134 : else{
135 : // set it to a very big value
136 0 : ci.associatedDistanceSq = 1e+12;
137 : }
138 0 : }
139 :
140 : ////////////////////////////////////////////////////////////////////////
141 0 : bool TriangleFill(std::vector<int>& vTriIndsOut, vector2* polyChain,
142 : size_t polyChainSize, bool bTriangulateInside)
143 : {
144 : vTriIndsOut.clear();
145 :
146 0 : if(polyChainSize < 2)
147 : return false;
148 :
149 : // calculate the normals
150 0 : vector<vector2> vNormals(polyChainSize);
151 0 : CalculatePolychainEdgeNormals(&vNormals.front(), polyChain,
152 : polyChainSize, bTriangulateInside);
153 :
154 : //tmp. create edges from the edges mid-point along their normals
155 :
156 : // set up a chain that holds the in- and out- edges of each point
157 : // together with a value that determines the priority.
158 : // The closer to 0, the higher the priority
159 :
160 : // chain infos are stored in a vector
161 : // this chain always holds valid and up to date elements
162 0 : vector<ChainInfo> chainInfos(polyChainSize);
163 :
164 : // create a priority queue that holds pointers to the chain-infos.
165 : // warning: this queue will contain invalid entries later on.
166 : // this is handled by a special check before an element is used.
167 : //priority_queue<ChainInfo> qChainInfos;
168 : queue<ChainInfo> qChainInfos;
169 0 : for(size_t i = 0; i < polyChainSize; ++i){
170 0 : int inVrtInd = (polyChainSize + i - 1) % polyChainSize;//avoid negativity in %
171 0 : int outVrtInd = (i + 1) % polyChainSize;
172 0 : chainInfos[i] = ChainInfo(i, inVrtInd, outVrtInd);
173 0 : UpdateChainInfo(chainInfos[i], polyChain, &vNormals.front());
174 0 : if(chainInfos[i].isCandidate)
175 : qChainInfos.push(chainInfos[i]);
176 : }
177 :
178 : // iterate until all triangles have been created.
179 0 : int counter = polyChainSize - 2;
180 0 : while(counter > 0 && (!qChainInfos.empty()))
181 : {
182 : // get the element with top priority
183 : //ChainInfo ci = qChainInfos.top();
184 0 : ChainInfo ci = qChainInfos.front();
185 : qChainInfos.pop();
186 0 : if(!ci.isCandidate){
187 : continue;
188 : }
189 :
190 : // we have to make sure that the chain info in the
191 : // priority queue matches the actual chain-entry.
192 : if(ci == chainInfos[ci.myVrt])
193 : {
194 : // check whether the triangle is a valid triangle or not.
195 : // do this by checking for all points whether they lie on the triangle.
196 0 : vector2& p0 = polyChain[ci.inVrt];
197 0 : vector2& p1 = polyChain[ci.myVrt];
198 0 : vector2& p2 = polyChain[ci.outVrt];
199 :
200 : // the bounding rect of the points
201 : vector2 vMin;
202 : vector2 vMax;
203 0 : vMin.x() = min(p0.x(), min(p1.x(), p2.x()));
204 0 : vMin.y() = min(p0.y(), min(p1.y(), p2.y()));
205 0 : vMax.x() = max(p0.x(), max(p1.x(), p2.x()));
206 0 : vMax.y() = max(p0.y(), max(p1.y(), p2.y()));
207 :
208 : // only consider points that lie in the box
209 : //TODO: replace this with a tree lookup.
210 : bool badTri = false;
211 0 : for(size_t i = 0; i < polyChainSize; ++i){
212 0 : if(i != ci.inVrt && i != ci.outVrt && i != ci.myVrt){
213 0 : vector2& p = polyChain[i];
214 0 : if(BoxBoundProbe(p, vMin, vMax))
215 : {
216 0 : if(PointIsInsideTriangle(p, p0, p1, p2)){
217 : //UG_LOG(" bad tri in " << ci.myVrt << " - ");
218 : badTri = true;
219 : break;
220 : }
221 : }
222 : }
223 : }
224 :
225 : // if no points lie in the triangle we can create it.
226 0 : if(!badTri){
227 0 : --counter;
228 0 : if(bTriangulateInside){
229 0 : vTriIndsOut.push_back(ci.inVrt);
230 0 : vTriIndsOut.push_back(ci.myVrt);
231 0 : vTriIndsOut.push_back(ci.outVrt);
232 : }
233 : else{
234 0 : vTriIndsOut.push_back(ci.outVrt);
235 0 : vTriIndsOut.push_back(ci.myVrt);
236 0 : vTriIndsOut.push_back(ci.inVrt);
237 : }
238 :
239 : // calculate the new normal of the edge
240 : vector2 eNew;
241 : VecSubtract(eNew, polyChain[ci.outVrt], polyChain[ci.inVrt]);
242 : vector2 nNew;
243 0 : nNew.x() = eNew.y();
244 0 : nNew.y() = -eNew.x();
245 : // make sure that the normal points into the right direction.
246 : vector2 vTestNorm;
247 : VecAdd(vTestNorm, vNormals[ci.inVrt], vNormals[ci.myVrt]);
248 0 : if(VecDot(nNew, vTestNorm) < 0)
249 : VecScale(nNew, nNew, -1.);
250 0 : VecNormalize(vNormals[ci.inVrt], nNew);
251 :
252 :
253 : // we have to update the chain-infos and the queue
254 0 : chainInfos[ci.myVrt].isCandidate = false;
255 : ChainInfo& ciIn = chainInfos[ci.inVrt];
256 : ChainInfo& ciOut = chainInfos[ci.outVrt];
257 0 : ciIn.outVrt = ci.outVrt;
258 0 : ciOut.inVrt = ci.inVrt;
259 0 : UpdateChainInfo(ciIn, polyChain, &vNormals.front());
260 0 : UpdateChainInfo(ciOut, polyChain, &vNormals.front());
261 0 : if(ciIn.isCandidate)
262 : qChainInfos.push(ciIn);
263 0 : if(ciOut.isCandidate)
264 : qChainInfos.push(ciOut);
265 : }
266 : }
267 : }
268 : /*
269 : if(counter == 0){
270 : UG_LOG(" counter is empty. ");
271 : }
272 : if(qChainInfos.empty()){
273 : UG_LOG(" queue is empty. ");
274 : }
275 : */
276 : // done. return true
277 : return true;
278 0 : }
279 :
280 0 : bool TriangleFill(Grid& grid, EdgeIterator edgesBegin,
281 : EdgeIterator edgesEnd, bool bTriangulateInside)
282 : {
283 0 : if(edgesBegin == edgesEnd)
284 : return false;
285 :
286 0 : if(!grid.has_vertex_attachment(aPosition))
287 : return false;
288 :
289 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPosition);
290 :
291 : // get the vertices of the poly-chain
292 : std::vector<Vertex*> vrtPolyChain;
293 0 : CreatePolyChain(vrtPolyChain, grid, edgesBegin, edgesEnd);
294 :
295 : // a poly-chain that stores positions
296 0 : std::vector<vector3> polyChain3D(vrtPolyChain.size());
297 0 : for(size_t i = 0; i < vrtPolyChain.size(); ++i)
298 0 : polyChain3D[i] = aaPos[vrtPolyChain[i]];
299 :
300 : // perform transform to 2d
301 0 : std::vector<vector2> polyChain2D(vrtPolyChain.size());
302 0 : TransformPointSetTo2D(&polyChain2D.front(), &polyChain3D.front(),
303 : polyChain3D.size());
304 : /*
305 : //TODO: perform a more elaborated projection!
306 : for(size_t i = 0; i < vrtPolyChain.size(); ++i){
307 : vector3& v = aaPos[vrtPolyChain[i]];
308 : polyChain2D[i] = vector2(v.x(), v.y());
309 : }
310 : */
311 :
312 : // perform triangulation
313 : vector<int> triIndices;
314 0 : if(TriangleFill(triIndices, &polyChain2D.front(), polyChain2D.size(),
315 : bTriangulateInside))
316 : {
317 : // create the triangles
318 0 : for(size_t i = 0; i < triIndices.size(); i+=3){
319 0 : grid.create<Triangle>(TriangleDescriptor(vrtPolyChain[triIndices[i]],
320 0 : vrtPolyChain[triIndices[i + 1]],
321 0 : vrtPolyChain[triIndices[i + 2]]));
322 : }
323 :
324 : return true;
325 : }
326 :
327 : /*
328 : //tmp. create edges from the edges mid-point along their normals
329 : // calculate the normals
330 : vector<vector2> vNormals(polyChain2D.size());
331 : CalculatePolychainEdgeNormals(&vNormals.front(), &polyChain2D.front(),
332 : polyChain2D.size(), bTriangulateInside);
333 :
334 : for(size_t i = 0; i < polyChain2D.size(); ++i){
335 : Vertex* vrt1 = *grid.create<RegularVertex>();
336 : Vertex* vrt2 = *grid.create<RegularVertex>();
337 : Edge* e = *grid.create<RegularEdge>(EdgeDescriptor(vrt1, vrt2));
338 :
339 : VecAdd(aaPos[vrt1], aaPos[vrtPolyChain[i]], aaPos[vrtPolyChain[(i+1)%vrtPolyChain.size()]]);
340 : VecScale(aaPos[vrt1], aaPos[vrt1], 0.5);
341 :
342 : aaPos[vrt2].x() = aaPos[vrt1].x() + vNormals[i].x() * 0.1;
343 : aaPos[vrt2].y() = aaPos[vrt1].y() + vNormals[i].y() * 0.1;
344 : aaPos[vrt2].z() = aaPos[vrt1].z();
345 : }
346 : return true;
347 : */
348 :
349 : return false;
350 0 : }
351 :
352 : }// end of namespace
|
