1 | /*
|
---|
2 | * Project: MoleCuilder
|
---|
3 | * Description: creates and alters molecular systems
|
---|
4 | * Copyright (C) 2010 University of Bonn. All rights reserved.
|
---|
5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
---|
6 | */
|
---|
7 |
|
---|
8 | /*
|
---|
9 | * TesselationHelpers.cpp
|
---|
10 | *
|
---|
11 | * Created on: Aug 3, 2009
|
---|
12 | * Author: heber
|
---|
13 | */
|
---|
14 |
|
---|
15 | // include config.h
|
---|
16 | #ifdef HAVE_CONFIG_H
|
---|
17 | #include <config.h>
|
---|
18 | #endif
|
---|
19 |
|
---|
20 | #include "Helpers/MemDebug.hpp"
|
---|
21 |
|
---|
22 | #include <fstream>
|
---|
23 |
|
---|
24 | #include "BoundaryLineSet.hpp"
|
---|
25 | #include "BoundaryPointSet.hpp"
|
---|
26 | #include "BoundaryPolygonSet.hpp"
|
---|
27 | #include "BoundaryTriangleSet.hpp"
|
---|
28 | #include "CandidateForTesselation.hpp"
|
---|
29 | #include "Helpers/Info.hpp"
|
---|
30 | #include "linkedcell.hpp"
|
---|
31 | #include "LinearAlgebra/linearsystemofequations.hpp"
|
---|
32 | #include "Helpers/Log.hpp"
|
---|
33 | #include "tesselation.hpp"
|
---|
34 | #include "tesselationhelpers.hpp"
|
---|
35 | #include "LinearAlgebra/Vector.hpp"
|
---|
36 | #include "LinearAlgebra/Line.hpp"
|
---|
37 | #include "LinearAlgebra/vector_ops.hpp"
|
---|
38 | #include "Helpers/Verbose.hpp"
|
---|
39 | #include "LinearAlgebra/Plane.hpp"
|
---|
40 | #include "LinearAlgebra/Matrix.hpp"
|
---|
41 |
|
---|
42 | void GetSphere(Vector * const center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS)
|
---|
43 | {
|
---|
44 | Info FunctionInfo(__func__);
|
---|
45 | Matrix mat;
|
---|
46 | double m11, m12, m13, m14;
|
---|
47 |
|
---|
48 | for(int i=0;i<3;i++) {
|
---|
49 | mat.set(i, 0, a[i]);
|
---|
50 | mat.set(i, 1, b[i]);
|
---|
51 | mat.set(i, 2, c[i]);
|
---|
52 | }
|
---|
53 | m11 = mat.determinant();
|
---|
54 |
|
---|
55 | for(int i=0;i<3;i++) {
|
---|
56 | mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
|
---|
57 | mat.set(i, 1, b[i]);
|
---|
58 | mat.set(i, 2, c[i]);
|
---|
59 | }
|
---|
60 | m12 = mat.determinant();
|
---|
61 |
|
---|
62 | for(int i=0;i<3;i++) {
|
---|
63 | mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
|
---|
64 | mat.set(i, 1, a[i]);
|
---|
65 | mat.set(i, 2, c[i]);
|
---|
66 | }
|
---|
67 | m13 = mat.determinant();
|
---|
68 |
|
---|
69 | for(int i=0;i<3;i++) {
|
---|
70 | mat.set(i, 0, a[i]*a[i] + b[i]*b[i] + c[i]*c[i]);
|
---|
71 | mat.set(i, 1, a[i]);
|
---|
72 | mat.set(i, 2, b[i]);
|
---|
73 | }
|
---|
74 | m14 = mat.determinant();
|
---|
75 |
|
---|
76 | if (fabs(m11) < MYEPSILON)
|
---|
77 | DoeLog(1) && (eLog()<< Verbose(1) << "three points are colinear." << endl);
|
---|
78 |
|
---|
79 | center->at(0) = 0.5 * m12/ m11;
|
---|
80 | center->at(1) = -0.5 * m13/ m11;
|
---|
81 | center->at(2) = 0.5 * m14/ m11;
|
---|
82 |
|
---|
83 | if (fabs(a.distance(*center) - RADIUS) > MYEPSILON)
|
---|
84 | DoeLog(1) && (eLog()<< Verbose(1) << "The given center is further way by " << fabs(a.distance(*center) - RADIUS) << " from a than RADIUS." << endl);
|
---|
85 | };
|
---|
86 |
|
---|
87 |
|
---|
88 |
|
---|
89 | /**
|
---|
90 | * Function returns center of sphere with RADIUS, which rests on points a, b, c
|
---|
91 | * @param Center this vector will be used for return
|
---|
92 | * @param a vector first point of triangle
|
---|
93 | * @param b vector second point of triangle
|
---|
94 | * @param c vector third point of triangle
|
---|
95 | * @param *Umkreismittelpunkt new center point of circumference
|
---|
96 | * @param Direction vector indicates up/down
|
---|
97 | * @param AlternativeDirection Vector, needed in case the triangles have 90 deg angle
|
---|
98 | * @param Halfplaneindicator double indicates whether Direction is up or down
|
---|
99 | * @param AlternativeIndicator double indicates in case of orthogonal triangles which direction of AlternativeDirection is suitable
|
---|
100 | * @param alpha double angle at a
|
---|
101 | * @param beta double, angle at b
|
---|
102 | * @param gamma, double, angle at c
|
---|
103 | * @param Radius, double
|
---|
104 | * @param Umkreisradius double radius of circumscribing circle
|
---|
105 | */
|
---|
106 | void GetCenterOfSphere(Vector* const & Center, const Vector &a, const Vector &b, const Vector &c, Vector * const NewUmkreismittelpunkt, const Vector* const Direction, const Vector* const AlternativeDirection,
|
---|
107 | const double HalfplaneIndicator, const double AlternativeIndicator, const double alpha, const double beta, const double gamma, const double RADIUS, const double Umkreisradius)
|
---|
108 | {
|
---|
109 | Info FunctionInfo(__func__);
|
---|
110 | Vector TempNormal, helper;
|
---|
111 | double Restradius;
|
---|
112 | Vector OtherCenter;
|
---|
113 | Center->Zero();
|
---|
114 | helper = sin(2.*alpha) * a;
|
---|
115 | (*Center) += helper;
|
---|
116 | helper = sin(2.*beta) * b;
|
---|
117 | (*Center) += helper;
|
---|
118 | helper = sin(2.*gamma) * c;
|
---|
119 | (*Center) += helper;
|
---|
120 | //*Center = a * sin(2.*alpha) + b * sin(2.*beta) + c * sin(2.*gamma) ;
|
---|
121 | Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
|
---|
122 | (*NewUmkreismittelpunkt) = (*Center);
|
---|
123 | DoLog(1) && (Log() << Verbose(1) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n");
|
---|
124 | // Here we calculated center of circumscribing circle, using barycentric coordinates
|
---|
125 | DoLog(1) && (Log() << Verbose(1) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n");
|
---|
126 |
|
---|
127 | TempNormal = a - b;
|
---|
128 | helper = a - c;
|
---|
129 | TempNormal.VectorProduct(helper);
|
---|
130 | if (fabs(HalfplaneIndicator) < MYEPSILON)
|
---|
131 | {
|
---|
132 | if ((TempNormal.ScalarProduct(*AlternativeDirection) <0 && AlternativeIndicator >0) || (TempNormal.ScalarProduct(*AlternativeDirection) >0 && AlternativeIndicator <0))
|
---|
133 | {
|
---|
134 | TempNormal *= -1;
|
---|
135 | }
|
---|
136 | }
|
---|
137 | else
|
---|
138 | {
|
---|
139 | if (((TempNormal.ScalarProduct(*Direction)<0) && (HalfplaneIndicator >0)) || ((TempNormal.ScalarProduct(*Direction)>0) && (HalfplaneIndicator<0)))
|
---|
140 | {
|
---|
141 | TempNormal *= -1;
|
---|
142 | }
|
---|
143 | }
|
---|
144 |
|
---|
145 | TempNormal.Normalize();
|
---|
146 | Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius);
|
---|
147 | DoLog(1) && (Log() << Verbose(1) << "Height of center of circumference to center of sphere is " << Restradius << ".\n");
|
---|
148 | TempNormal.Scale(Restradius);
|
---|
149 | DoLog(1) && (Log() << Verbose(1) << "Shift vector to sphere of circumference is " << TempNormal << ".\n");
|
---|
150 | (*Center) += TempNormal;
|
---|
151 | DoLog(1) && (Log() << Verbose(1) << "Center of sphere of circumference is " << *Center << ".\n");
|
---|
152 | GetSphere(&OtherCenter, a, b, c, RADIUS);
|
---|
153 | DoLog(1) && (Log() << Verbose(1) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n");
|
---|
154 | };
|
---|
155 |
|
---|
156 |
|
---|
157 | /** Constructs the center of the circumcircle defined by three points \a *a, \a *b and \a *c.
|
---|
158 | * \param *Center new center on return
|
---|
159 | * \param *a first point
|
---|
160 | * \param *b second point
|
---|
161 | * \param *c third point
|
---|
162 | */
|
---|
163 | void GetCenterofCircumcircle(Vector &Center, const Vector &a, const Vector &b, const Vector &c)
|
---|
164 | {
|
---|
165 | Info FunctionInfo(__func__);
|
---|
166 | Vector helper;
|
---|
167 | Vector SideA = b - c;
|
---|
168 | Vector SideB = c - a;
|
---|
169 | Vector SideC = a - b;
|
---|
170 |
|
---|
171 | helper[0] = SideA.NormSquared()*(SideB.NormSquared()+SideC.NormSquared() - SideA.NormSquared());
|
---|
172 | helper[1] = SideB.NormSquared()*(SideC.NormSquared()+SideA.NormSquared() - SideB.NormSquared());
|
---|
173 | helper[2] = SideC.NormSquared()*(SideA.NormSquared()+SideB.NormSquared() - SideC.NormSquared());
|
---|
174 |
|
---|
175 | Center.Zero();
|
---|
176 | Center += helper[0] * a;
|
---|
177 | Center += helper[1] * b;
|
---|
178 | Center += helper[2] * c;
|
---|
179 | if (fabs(helper[0]+helper[1]+helper[2]) > MYEPSILON)
|
---|
180 | Center.Scale(1./(helper[0]+helper[1]+helper[2]));
|
---|
181 | Log() << Verbose(1) << "INFO: Center (2nd algo) is at " << Center << "." << endl;
|
---|
182 | };
|
---|
183 |
|
---|
184 | /** Returns the parameter "path length" for a given \a NewSphereCenter relative to \a OldSphereCenter on a circle on the plane \a CirclePlaneNormal with center \a CircleCenter and radius \a CircleRadius.
|
---|
185 | * Test whether the \a NewSphereCenter is really on the given plane and in distance \a CircleRadius from \a CircleCenter.
|
---|
186 | * It calculates the angle, making it unique on [0,2.*M_PI) by comparing to SearchDirection.
|
---|
187 | * Also the new center is invalid if it the same as the old one and does not lie right above (\a NormalVector) the base line (\a CircleCenter).
|
---|
188 | * \param CircleCenter Center of the parameter circle
|
---|
189 | * \param CirclePlaneNormal normal vector to plane of the parameter circle
|
---|
190 | * \param CircleRadius radius of the parameter circle
|
---|
191 | * \param NewSphereCenter new center of a circumcircle
|
---|
192 | * \param OldSphereCenter old center of a circumcircle, defining the zero "path length" on the parameter circle
|
---|
193 | * \param NormalVector normal vector
|
---|
194 | * \param SearchDirection search direction to make angle unique on return.
|
---|
195 | * \param HULLEPSILON machine precision for tesselation points
|
---|
196 | * \return Angle between \a NewSphereCenter and \a OldSphereCenter relative to \a CircleCenter, 2.*M_PI if one test fails
|
---|
197 | */
|
---|
198 | double GetPathLengthonCircumCircle(const Vector &CircleCenter, const Vector &CirclePlaneNormal, const double CircleRadius, const Vector &NewSphereCenter, const Vector &OldSphereCenter, const Vector &NormalVector, const Vector &SearchDirection, const double HULLEPSILON)
|
---|
199 | {
|
---|
200 | Info FunctionInfo(__func__);
|
---|
201 | Vector helper;
|
---|
202 | double radius, alpha;
|
---|
203 |
|
---|
204 | Vector RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
|
---|
205 | Vector RelativeNewSphereCenter = NewSphereCenter - CircleCenter;
|
---|
206 | helper = RelativeNewSphereCenter;
|
---|
207 | // test whether new center is on the parameter circle's plane
|
---|
208 | if (fabs(helper.ScalarProduct(CirclePlaneNormal)) > HULLEPSILON) {
|
---|
209 | DoeLog(1) && (eLog()<< Verbose(1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(CirclePlaneNormal)) << "!" << endl);
|
---|
210 | helper.ProjectOntoPlane(CirclePlaneNormal);
|
---|
211 | }
|
---|
212 | radius = helper.NormSquared();
|
---|
213 | // test whether the new center vector has length of CircleRadius
|
---|
214 | if (fabs(radius - CircleRadius) > HULLEPSILON)
|
---|
215 | DoeLog(1) && (eLog()<< Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl);
|
---|
216 | alpha = helper.Angle(RelativeOldSphereCenter);
|
---|
217 | // make the angle unique by checking the halfplanes/search direction
|
---|
218 | if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON) // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
|
---|
219 | alpha = 2.*M_PI - alpha;
|
---|
220 | DoLog(1) && (Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << "." << endl);
|
---|
221 | radius = helper.distance(RelativeOldSphereCenter);
|
---|
222 | helper.ProjectOntoPlane(NormalVector);
|
---|
223 | // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
|
---|
224 | if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
|
---|
225 | DoLog(1) && (Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl);
|
---|
226 | return alpha;
|
---|
227 | } else {
|
---|
228 | DoLog(1) && (Log() << Verbose(1) << "INFO: NewSphereCenter " << RelativeNewSphereCenter << " is too close to RelativeOldSphereCenter" << RelativeOldSphereCenter << "." << endl);
|
---|
229 | return 2.*M_PI;
|
---|
230 | }
|
---|
231 | };
|
---|
232 |
|
---|
233 | struct Intersection {
|
---|
234 | Vector x1;
|
---|
235 | Vector x2;
|
---|
236 | Vector x3;
|
---|
237 | Vector x4;
|
---|
238 | };
|
---|
239 |
|
---|
240 | /** Gets the angle between a point and a reference relative to the provided center.
|
---|
241 | * We have two shanks point and reference between which the angle is calculated
|
---|
242 | * and by scalar product with OrthogonalVector we decide the interval.
|
---|
243 | * @param point to calculate the angle for
|
---|
244 | * @param reference to which to calculate the angle
|
---|
245 | * @param OrthogonalVector points in direction of [pi,2pi] interval
|
---|
246 | *
|
---|
247 | * @return angle between point and reference
|
---|
248 | */
|
---|
249 | double GetAngle(const Vector &point, const Vector &reference, const Vector &OrthogonalVector)
|
---|
250 | {
|
---|
251 | Info FunctionInfo(__func__);
|
---|
252 | if (reference.IsZero())
|
---|
253 | return M_PI;
|
---|
254 |
|
---|
255 | // calculate both angles and correct with in-plane vector
|
---|
256 | if (point.IsZero())
|
---|
257 | return M_PI;
|
---|
258 | double phi = point.Angle(reference);
|
---|
259 | if (OrthogonalVector.ScalarProduct(point) > 0) {
|
---|
260 | phi = 2.*M_PI - phi;
|
---|
261 | }
|
---|
262 |
|
---|
263 | DoLog(1) && (Log() << Verbose(1) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl);
|
---|
264 |
|
---|
265 | return phi;
|
---|
266 | }
|
---|
267 |
|
---|
268 |
|
---|
269 | /** Calculates the volume of a general tetraeder.
|
---|
270 | * \param *a first vector
|
---|
271 | * \param *b second vector
|
---|
272 | * \param *c third vector
|
---|
273 | * \param *d fourth vector
|
---|
274 | * \return \f$ \frac{1}{6} \cdot ((a-d) \times (a-c) \cdot (a-b)) \f$
|
---|
275 | */
|
---|
276 | double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d)
|
---|
277 | {
|
---|
278 | Info FunctionInfo(__func__);
|
---|
279 | Vector Point, TetraederVector[3];
|
---|
280 | double volume;
|
---|
281 |
|
---|
282 | TetraederVector[0] = a;
|
---|
283 | TetraederVector[1] = b;
|
---|
284 | TetraederVector[2] = c;
|
---|
285 | for (int j=0;j<3;j++)
|
---|
286 | TetraederVector[j].SubtractVector(d);
|
---|
287 | Point = TetraederVector[0];
|
---|
288 | Point.VectorProduct(TetraederVector[1]);
|
---|
289 | volume = 1./6. * fabs(Point.ScalarProduct(TetraederVector[2]));
|
---|
290 | return volume;
|
---|
291 | };
|
---|
292 |
|
---|
293 | /** Calculates the area of a general triangle.
|
---|
294 | * We use the Heron's formula of area, [Bronstein, S. 138]
|
---|
295 | * \param &A first vector
|
---|
296 | * \param &B second vector
|
---|
297 | * \param &C third vector
|
---|
298 | * \return \f$ \frac{1}{6} \cdot ((a-d) \times (a-c) \cdot (a-b)) \f$
|
---|
299 | */
|
---|
300 | double CalculateAreaofGeneralTriangle(const Vector &A, const Vector &B, const Vector &C)
|
---|
301 | {
|
---|
302 | Info FunctionInfo(__func__);
|
---|
303 |
|
---|
304 | const double sidea = B.distance(C);
|
---|
305 | const double sideb = A.distance(C);
|
---|
306 | const double sidec = A.distance(B);
|
---|
307 | const double s = (sidea+sideb+sidec)/2.;
|
---|
308 |
|
---|
309 | const double area = sqrt(s*(s-sidea)*(s-sideb)*(s-sidec));
|
---|
310 | return area;
|
---|
311 | };
|
---|
312 |
|
---|
313 |
|
---|
314 | /** Checks for a new special triangle whether one of its edges is already present with one one triangle connected.
|
---|
315 | * This enforces that special triangles (i.e. degenerated ones) should at last close the open-edge frontier and not
|
---|
316 | * make it bigger (i.e. closing one (the baseline) and opening two new ones).
|
---|
317 | * \param TPS[3] nodes of the triangle
|
---|
318 | * \return true - there is such a line (i.e. creation of degenerated triangle is valid), false - no such line (don't create)
|
---|
319 | */
|
---|
320 | bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet * const nodes[3])
|
---|
321 | {
|
---|
322 | Info FunctionInfo(__func__);
|
---|
323 | bool result = false;
|
---|
324 | int counter = 0;
|
---|
325 |
|
---|
326 | // check all three points
|
---|
327 | for (int i=0;i<3;i++)
|
---|
328 | for (int j=i+1; j<3; j++) {
|
---|
329 | if (nodes[i] == NULL) {
|
---|
330 | DoLog(1) && (Log() << Verbose(1) << "Node nr. " << i << " is not yet present." << endl);
|
---|
331 | result = true;
|
---|
332 | } else if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) { // there already is a line
|
---|
333 | LineMap::const_iterator FindLine;
|
---|
334 | pair<LineMap::const_iterator,LineMap::const_iterator> FindPair;
|
---|
335 | FindPair = nodes[i]->lines.equal_range(nodes[j]->node->nr);
|
---|
336 | for (FindLine = FindPair.first; FindLine != FindPair.second; ++FindLine) {
|
---|
337 | // If there is a line with less than two attached triangles, we don't need a new line.
|
---|
338 | if (FindLine->second->triangles.size() < 2) {
|
---|
339 | counter++;
|
---|
340 | break; // increase counter only once per edge
|
---|
341 | }
|
---|
342 | }
|
---|
343 | } else { // no line
|
---|
344 | DoLog(1) && (Log() << Verbose(1) << "The line between " << *nodes[i] << " and " << *nodes[j] << " is not yet present, hence no need for a degenerate triangle." << endl);
|
---|
345 | result = true;
|
---|
346 | }
|
---|
347 | }
|
---|
348 | if ((!result) && (counter > 1)) {
|
---|
349 | DoLog(1) && (Log() << Verbose(1) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl);
|
---|
350 | result = true;
|
---|
351 | }
|
---|
352 | return result;
|
---|
353 | };
|
---|
354 |
|
---|
355 |
|
---|
356 | ///** Sort function for the candidate list.
|
---|
357 | // */
|
---|
358 | //bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2)
|
---|
359 | //{
|
---|
360 | // Info FunctionInfo(__func__);
|
---|
361 | // Vector BaseLineVector, OrthogonalVector, helper;
|
---|
362 | // if (candidate1->BaseLine != candidate2->BaseLine) { // sanity check
|
---|
363 | // DoeLog(1) && (eLog()<< Verbose(1) << "sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl);
|
---|
364 | // //return false;
|
---|
365 | // exit(1);
|
---|
366 | // }
|
---|
367 | // // create baseline vector
|
---|
368 | // BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
|
---|
369 | // BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
|
---|
370 | // BaseLineVector.Normalize();
|
---|
371 | //
|
---|
372 | // // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!)
|
---|
373 | // helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
|
---|
374 | // helper.SubtractVector(candidate1->point->node);
|
---|
375 | // OrthogonalVector.CopyVector(&helper);
|
---|
376 | // helper.VectorProduct(&BaseLineVector);
|
---|
377 | // OrthogonalVector.SubtractVector(&helper);
|
---|
378 | // OrthogonalVector.Normalize();
|
---|
379 | //
|
---|
380 | // // calculate both angles and correct with in-plane vector
|
---|
381 | // helper.CopyVector(candidate1->point->node);
|
---|
382 | // helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
|
---|
383 | // double phi = BaseLineVector.Angle(&helper);
|
---|
384 | // if (OrthogonalVector.ScalarProduct(&helper) > 0) {
|
---|
385 | // phi = 2.*M_PI - phi;
|
---|
386 | // }
|
---|
387 | // helper.CopyVector(candidate2->point->node);
|
---|
388 | // helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
|
---|
389 | // double psi = BaseLineVector.Angle(&helper);
|
---|
390 | // if (OrthogonalVector.ScalarProduct(&helper) > 0) {
|
---|
391 | // psi = 2.*M_PI - psi;
|
---|
392 | // }
|
---|
393 | //
|
---|
394 | // Log() << Verbose(1) << *candidate1->point << " has angle " << phi << endl;
|
---|
395 | // Log() << Verbose(1) << *candidate2->point << " has angle " << psi << endl;
|
---|
396 | //
|
---|
397 | // // return comparison
|
---|
398 | // return phi < psi;
|
---|
399 | //};
|
---|
400 |
|
---|
401 | /**
|
---|
402 | * Finds the point which is second closest to the provided one.
|
---|
403 | *
|
---|
404 | * @param Point to which to find the second closest other point
|
---|
405 | * @param linked cell structure
|
---|
406 | *
|
---|
407 | * @return point which is second closest to the provided one
|
---|
408 | */
|
---|
409 | TesselPoint* FindSecondClosestTesselPoint(const Vector& Point, const LinkedCell* const LC)
|
---|
410 | {
|
---|
411 | Info FunctionInfo(__func__);
|
---|
412 | TesselPoint* closestPoint = NULL;
|
---|
413 | TesselPoint* secondClosestPoint = NULL;
|
---|
414 | double distance = 1e16;
|
---|
415 | double secondDistance = 1e16;
|
---|
416 | Vector helper;
|
---|
417 | int N[NDIM], Nlower[NDIM], Nupper[NDIM];
|
---|
418 |
|
---|
419 | LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
|
---|
420 | for(int i=0;i<NDIM;i++) // store indices of this cell
|
---|
421 | N[i] = LC->n[i];
|
---|
422 | DoLog(1) && (Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl);
|
---|
423 |
|
---|
424 | LC->GetNeighbourBounds(Nlower, Nupper);
|
---|
425 | //Log() << Verbose(1) << endl;
|
---|
426 | for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
|
---|
427 | for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
|
---|
428 | for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
|
---|
429 | const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
|
---|
430 | //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
|
---|
431 | if (List != NULL) {
|
---|
432 | for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
|
---|
433 | helper = (Point) - ((*Runner)->getPosition());
|
---|
434 | double currentNorm = helper. Norm();
|
---|
435 | if (currentNorm < distance) {
|
---|
436 | // remember second point
|
---|
437 | secondDistance = distance;
|
---|
438 | secondClosestPoint = closestPoint;
|
---|
439 | // mark down new closest point
|
---|
440 | distance = currentNorm;
|
---|
441 | closestPoint = (*Runner);
|
---|
442 | //Log() << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *secondClosestPoint << "." << endl;
|
---|
443 | }
|
---|
444 | }
|
---|
445 | } else {
|
---|
446 | DoeLog(1) && (eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
|
---|
447 | }
|
---|
448 | }
|
---|
449 |
|
---|
450 | return secondClosestPoint;
|
---|
451 | };
|
---|
452 |
|
---|
453 | /**
|
---|
454 | * Finds the point which is closest to the provided one.
|
---|
455 | *
|
---|
456 | * @param Point to which to find the closest other point
|
---|
457 | * @param SecondPoint the second closest other point on return, NULL if none found
|
---|
458 | * @param linked cell structure
|
---|
459 | *
|
---|
460 | * @return point which is closest to the provided one, NULL if none found
|
---|
461 | */
|
---|
462 | TesselPoint* FindClosestTesselPoint(const Vector& Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
|
---|
463 | {
|
---|
464 | Info FunctionInfo(__func__);
|
---|
465 | TesselPoint* closestPoint = NULL;
|
---|
466 | SecondPoint = NULL;
|
---|
467 | double distance = 1e16;
|
---|
468 | double secondDistance = 1e16;
|
---|
469 | Vector helper;
|
---|
470 | int N[NDIM], Nlower[NDIM], Nupper[NDIM];
|
---|
471 |
|
---|
472 | LC->SetIndexToVector(Point); // ignore status as we calculate bounds below sensibly
|
---|
473 | for(int i=0;i<NDIM;i++) // store indices of this cell
|
---|
474 | N[i] = LC->n[i];
|
---|
475 | DoLog(1) && (Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl);
|
---|
476 |
|
---|
477 | LC->GetNeighbourBounds(Nlower, Nupper);
|
---|
478 | //Log() << Verbose(1) << endl;
|
---|
479 | for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
|
---|
480 | for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
|
---|
481 | for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
|
---|
482 | const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
|
---|
483 | //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
|
---|
484 | if (List != NULL) {
|
---|
485 | for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
|
---|
486 | helper = (Point) - ((*Runner)->getPosition());
|
---|
487 | double currentNorm = helper.NormSquared();
|
---|
488 | if (currentNorm < distance) {
|
---|
489 | secondDistance = distance;
|
---|
490 | SecondPoint = closestPoint;
|
---|
491 | distance = currentNorm;
|
---|
492 | closestPoint = (*Runner);
|
---|
493 | //Log() << Verbose(1) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
|
---|
494 | } else if (currentNorm < secondDistance) {
|
---|
495 | secondDistance = currentNorm;
|
---|
496 | SecondPoint = (*Runner);
|
---|
497 | //Log() << Verbose(1) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
|
---|
498 | }
|
---|
499 | }
|
---|
500 | } else {
|
---|
501 | DoeLog(1) && (eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
|
---|
502 | }
|
---|
503 | }
|
---|
504 | // output
|
---|
505 | if (closestPoint != NULL) {
|
---|
506 | DoLog(1) && (Log() << Verbose(1) << "Closest point is " << *closestPoint);
|
---|
507 | if (SecondPoint != NULL)
|
---|
508 | DoLog(0) && (Log() << Verbose(0) << " and second closest is " << *SecondPoint);
|
---|
509 | DoLog(0) && (Log() << Verbose(0) << "." << endl);
|
---|
510 | }
|
---|
511 | return closestPoint;
|
---|
512 | };
|
---|
513 |
|
---|
514 | /** Returns the closest point on \a *Base with respect to \a *OtherBase.
|
---|
515 | * \param *out output stream for debugging
|
---|
516 | * \param *Base reference line
|
---|
517 | * \param *OtherBase other base line
|
---|
518 | * \return Vector on reference line that has closest distance
|
---|
519 | */
|
---|
520 | Vector * GetClosestPointBetweenLine(const BoundaryLineSet * const Base, const BoundaryLineSet * const OtherBase)
|
---|
521 | {
|
---|
522 | Info FunctionInfo(__func__);
|
---|
523 | // construct the plane of the two baselines (i.e. take both their directional vectors)
|
---|
524 | Vector Baseline = (Base->endpoints[1]->node->getPosition()) - (Base->endpoints[0]->node->getPosition());
|
---|
525 | Vector OtherBaseline = (OtherBase->endpoints[1]->node->getPosition()) - (OtherBase->endpoints[0]->node->getPosition());
|
---|
526 | Vector Normal = Baseline;
|
---|
527 | Normal.VectorProduct(OtherBaseline);
|
---|
528 | Normal.Normalize();
|
---|
529 | DoLog(1) && (Log() << Verbose(1) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl);
|
---|
530 |
|
---|
531 | // project one offset point of OtherBase onto this plane (and add plane offset vector)
|
---|
532 | Vector NewOffset = (OtherBase->endpoints[0]->node->getPosition()) - (Base->endpoints[0]->node->getPosition());
|
---|
533 | NewOffset.ProjectOntoPlane(Normal);
|
---|
534 | NewOffset += (Base->endpoints[0]->node->getPosition());
|
---|
535 | Vector NewDirection = NewOffset + OtherBaseline;
|
---|
536 |
|
---|
537 | // calculate the intersection between this projected baseline and Base
|
---|
538 | Vector *Intersection = new Vector;
|
---|
539 | Line line1 = makeLineThrough((Base->endpoints[0]->node->getPosition()),(Base->endpoints[1]->node->getPosition()));
|
---|
540 | Line line2 = makeLineThrough(NewOffset, NewDirection);
|
---|
541 | *Intersection = line1.getIntersection(line2);
|
---|
542 | Normal = (*Intersection) - (Base->endpoints[0]->node->getPosition());
|
---|
543 | DoLog(1) && (Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << "." << endl);
|
---|
544 |
|
---|
545 | return Intersection;
|
---|
546 | };
|
---|
547 |
|
---|
548 | /** Returns the distance to the plane defined by \a *triangle
|
---|
549 | * \param *out output stream for debugging
|
---|
550 | * \param *x Vector to calculate distance to
|
---|
551 | * \param *triangle triangle defining plane
|
---|
552 | * \return distance between \a *x and plane defined by \a *triangle, -1 - if something went wrong
|
---|
553 | */
|
---|
554 | double DistanceToTrianglePlane(const Vector *x, const BoundaryTriangleSet * const triangle)
|
---|
555 | {
|
---|
556 | Info FunctionInfo(__func__);
|
---|
557 | double distance = 0.;
|
---|
558 | if (x == NULL) {
|
---|
559 | return -1;
|
---|
560 | }
|
---|
561 | distance = x->DistanceToSpace(triangle->getPlane());
|
---|
562 | return distance;
|
---|
563 | };
|
---|
564 |
|
---|
565 | /** Creates the objects in a VRML file.
|
---|
566 | * \param *out output stream for debugging
|
---|
567 | * \param *vrmlfile output stream for tecplot data
|
---|
568 | * \param *Tess Tesselation structure with constructed triangles
|
---|
569 | * \param *mol molecule structure with atom positions
|
---|
570 | */
|
---|
571 | void WriteVrmlFile(ofstream * const vrmlfile, const Tesselation * const Tess, const PointCloud * const cloud)
|
---|
572 | {
|
---|
573 | Info FunctionInfo(__func__);
|
---|
574 | TesselPoint *Walker = NULL;
|
---|
575 | int i;
|
---|
576 | Vector *center = cloud->GetCenter();
|
---|
577 | if (vrmlfile != NULL) {
|
---|
578 | //Log() << Verbose(1) << "Writing Raster3D file ... ";
|
---|
579 | *vrmlfile << "#VRML V2.0 utf8" << endl;
|
---|
580 | *vrmlfile << "#Created by molecuilder" << endl;
|
---|
581 | *vrmlfile << "#All atoms as spheres" << endl;
|
---|
582 | cloud->GoToFirst();
|
---|
583 | while (!cloud->IsEnd()) {
|
---|
584 | Walker = cloud->GetPoint();
|
---|
585 | *vrmlfile << "Sphere {" << endl << " "; // 2 is sphere type
|
---|
586 | for (i=0;i<NDIM;i++)
|
---|
587 | *vrmlfile << Walker->at(i)-center->at(i) << " ";
|
---|
588 | *vrmlfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
|
---|
589 | cloud->GoToNext();
|
---|
590 | }
|
---|
591 |
|
---|
592 | *vrmlfile << "# All tesselation triangles" << endl;
|
---|
593 | for (TriangleMap::const_iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
|
---|
594 | *vrmlfile << "1" << endl << " "; // 1 is triangle type
|
---|
595 | for (i=0;i<3;i++) { // print each node
|
---|
596 | for (int j=0;j<NDIM;j++) // and for each node all NDIM coordinates
|
---|
597 | *vrmlfile << TriangleRunner->second->endpoints[i]->node->at(j)-center->at(j) << " ";
|
---|
598 | *vrmlfile << "\t";
|
---|
599 | }
|
---|
600 | *vrmlfile << "1. 0. 0." << endl; // red as colour
|
---|
601 | *vrmlfile << "18" << endl << " 0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
|
---|
602 | }
|
---|
603 | } else {
|
---|
604 | DoeLog(1) && (eLog()<< Verbose(1) << "Given vrmlfile is " << vrmlfile << "." << endl);
|
---|
605 | }
|
---|
606 | delete(center);
|
---|
607 | };
|
---|
608 |
|
---|
609 | /** Writes additionally the current sphere (i.e. the last triangle to file).
|
---|
610 | * \param *out output stream for debugging
|
---|
611 | * \param *rasterfile output stream for tecplot data
|
---|
612 | * \param *Tess Tesselation structure with constructed triangles
|
---|
613 | * \param *mol molecule structure with atom positions
|
---|
614 | */
|
---|
615 | void IncludeSphereinRaster3D(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
|
---|
616 | {
|
---|
617 | Info FunctionInfo(__func__);
|
---|
618 | Vector helper;
|
---|
619 |
|
---|
620 | if (Tess->LastTriangle != NULL) {
|
---|
621 | // include the current position of the virtual sphere in the temporary raster3d file
|
---|
622 | Vector *center = cloud->GetCenter();
|
---|
623 | // make the circumsphere's center absolute again
|
---|
624 | Vector helper = (1./3.) * ((Tess->LastTriangle->endpoints[0]->node->getPosition()) +
|
---|
625 | (Tess->LastTriangle->endpoints[1]->node->getPosition()) +
|
---|
626 | (Tess->LastTriangle->endpoints[2]->node->getPosition()));
|
---|
627 | helper -= (*center);
|
---|
628 | // and add to file plus translucency object
|
---|
629 | *rasterfile << "# current virtual sphere\n";
|
---|
630 | *rasterfile << "8\n 25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0\n";
|
---|
631 | *rasterfile << "2\n " << helper[0] << " " << helper[1] << " " << helper[2] << "\t" << 5. << "\t1 0 0\n";
|
---|
632 | *rasterfile << "9\n terminating special property\n";
|
---|
633 | delete(center);
|
---|
634 | }
|
---|
635 | };
|
---|
636 |
|
---|
637 | /** Creates the objects in a raster3d file (renderable with a header.r3d).
|
---|
638 | * \param *out output stream for debugging
|
---|
639 | * \param *rasterfile output stream for tecplot data
|
---|
640 | * \param *Tess Tesselation structure with constructed triangles
|
---|
641 | * \param *mol molecule structure with atom positions
|
---|
642 | */
|
---|
643 | void WriteRaster3dFile(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
|
---|
644 | {
|
---|
645 | Info FunctionInfo(__func__);
|
---|
646 | TesselPoint *Walker = NULL;
|
---|
647 | int i;
|
---|
648 | Vector *center = cloud->GetCenter();
|
---|
649 | if (rasterfile != NULL) {
|
---|
650 | //Log() << Verbose(1) << "Writing Raster3D file ... ";
|
---|
651 | *rasterfile << "# Raster3D object description, created by MoleCuilder" << endl;
|
---|
652 | *rasterfile << "@header.r3d" << endl;
|
---|
653 | *rasterfile << "# All atoms as spheres" << endl;
|
---|
654 | cloud->GoToFirst();
|
---|
655 | while (!cloud->IsEnd()) {
|
---|
656 | Walker = cloud->GetPoint();
|
---|
657 | *rasterfile << "2" << endl << " "; // 2 is sphere type
|
---|
658 | for (int j=0;j<NDIM;j++) { // and for each node all NDIM coordinates
|
---|
659 | const double tmp = Walker->at(j)-center->at(j);
|
---|
660 | *rasterfile << ((fabs(tmp) < MYEPSILON) ? 0 : tmp) << " ";
|
---|
661 | }
|
---|
662 | *rasterfile << "\t0.1\t1. 1. 1." << endl; // radius 0.05 and white as colour
|
---|
663 | cloud->GoToNext();
|
---|
664 | }
|
---|
665 |
|
---|
666 | *rasterfile << "# All tesselation triangles" << endl;
|
---|
667 | *rasterfile << "8\n 25. -1. 1. 1. 1. 0.0 0 0 0 2\n SOLID 1.0 0.0 0.0\n BACKFACE 0.3 0.3 1.0 0 0\n";
|
---|
668 | for (TriangleMap::const_iterator TriangleRunner = Tess->TrianglesOnBoundary.begin(); TriangleRunner != Tess->TrianglesOnBoundary.end(); TriangleRunner++) {
|
---|
669 | *rasterfile << "1" << endl << " "; // 1 is triangle type
|
---|
670 | for (i=0;i<3;i++) { // print each node
|
---|
671 | for (int j=0;j<NDIM;j++) { // and for each node all NDIM coordinates
|
---|
672 | const double tmp = TriangleRunner->second->endpoints[i]->node->at(j)-center->at(j);
|
---|
673 | *rasterfile << ((fabs(tmp) < MYEPSILON) ? 0 : tmp) << " ";
|
---|
674 | }
|
---|
675 | *rasterfile << "\t";
|
---|
676 | }
|
---|
677 | *rasterfile << "1. 0. 0." << endl; // red as colour
|
---|
678 | //*rasterfile << "18" << endl << " 0.5 0.5 0.5" << endl; // 18 is transparency type for previous object
|
---|
679 | }
|
---|
680 | *rasterfile << "9\n# terminating special property\n";
|
---|
681 | } else {
|
---|
682 | DoeLog(1) && (eLog()<< Verbose(1) << "Given rasterfile is " << rasterfile << "." << endl);
|
---|
683 | }
|
---|
684 | IncludeSphereinRaster3D(rasterfile, Tess, cloud);
|
---|
685 | delete(center);
|
---|
686 | };
|
---|
687 |
|
---|
688 | /** This function creates the tecplot file, displaying the tesselation of the hull.
|
---|
689 | * \param *out output stream for debugging
|
---|
690 | * \param *tecplot output stream for tecplot data
|
---|
691 | * \param N arbitrary number to differentiate various zones in the tecplot format
|
---|
692 | */
|
---|
693 | void WriteTecplotFile(ofstream * const tecplot, const Tesselation * const TesselStruct, const PointCloud * const cloud, const int N)
|
---|
694 | {
|
---|
695 | Info FunctionInfo(__func__);
|
---|
696 | if ((tecplot != NULL) && (TesselStruct != NULL)) {
|
---|
697 | // write header
|
---|
698 | *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
|
---|
699 | *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl;
|
---|
700 | *tecplot << "ZONE T=\"";
|
---|
701 | if (N < 0) {
|
---|
702 | *tecplot << cloud->GetName();
|
---|
703 | } else {
|
---|
704 | *tecplot << N << "-";
|
---|
705 | if (TesselStruct->LastTriangle != NULL) {
|
---|
706 | for (int i=0;i<3;i++)
|
---|
707 | *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->getName();
|
---|
708 | } else {
|
---|
709 | *tecplot << "none";
|
---|
710 | }
|
---|
711 | }
|
---|
712 | *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
|
---|
713 | const int MaxId=cloud->GetMaxId();
|
---|
714 | int *LookupList = new int[MaxId];
|
---|
715 | for (int i=0; i< MaxId ; i++){
|
---|
716 | LookupList[i] = -1;
|
---|
717 | }
|
---|
718 |
|
---|
719 | // print atom coordinates
|
---|
720 | int Counter = 1;
|
---|
721 | TesselPoint *Walker = NULL;
|
---|
722 | for (PointMap::const_iterator target = TesselStruct->PointsOnBoundary.begin(); target != TesselStruct->PointsOnBoundary.end(); ++target) {
|
---|
723 | Walker = target->second->node;
|
---|
724 | LookupList[Walker->nr] = Counter++;
|
---|
725 | for (int i=0;i<NDIM;i++) {
|
---|
726 | const double tmp = Walker->at(i);
|
---|
727 | *tecplot << ((fabs(tmp) < MYEPSILON) ? 0 : tmp) << " ";
|
---|
728 | }
|
---|
729 | *tecplot << target->second->value << endl;
|
---|
730 | }
|
---|
731 | *tecplot << endl;
|
---|
732 | // print connectivity
|
---|
733 | DoLog(1) && (Log() << Verbose(1) << "The following triangles were created:" << endl);
|
---|
734 | for (TriangleMap::const_iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
|
---|
735 | DoLog(1) && (Log() << Verbose(1) << " " << runner->second->endpoints[0]->node->getName() << "<->" << runner->second->endpoints[1]->node->getName() << "<->" << runner->second->endpoints[2]->node->getName() << endl);
|
---|
736 | *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
|
---|
737 | }
|
---|
738 | delete[] (LookupList);
|
---|
739 | }
|
---|
740 | };
|
---|
741 |
|
---|
742 | /** Calculates the concavity for each of the BoundaryPointSet's in a Tesselation.
|
---|
743 | * Sets BoundaryPointSet::value equal to the number of connected lines that are not convex.
|
---|
744 | * \param *out output stream for debugging
|
---|
745 | * \param *TesselStruct pointer to Tesselation structure
|
---|
746 | */
|
---|
747 | void CalculateConcavityPerBoundaryPoint(const Tesselation * const TesselStruct)
|
---|
748 | {
|
---|
749 | Info FunctionInfo(__func__);
|
---|
750 | class BoundaryPointSet *point = NULL;
|
---|
751 | class BoundaryLineSet *line = NULL;
|
---|
752 | class BoundaryTriangleSet *triangle = NULL;
|
---|
753 | double ConcavityPerLine = 0.;
|
---|
754 | double ConcavityPerTriangle = 0.;
|
---|
755 | double area = 0.;
|
---|
756 | double totalarea = 0.;
|
---|
757 |
|
---|
758 | for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
|
---|
759 | point = PointRunner->second;
|
---|
760 | DoLog(1) && (Log() << Verbose(1) << "INFO: Current point is " << *point << "." << endl);
|
---|
761 |
|
---|
762 | // calculate mean concavity over all connected line
|
---|
763 | ConcavityPerLine = 0.;
|
---|
764 | for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
|
---|
765 | line = LineRunner->second;
|
---|
766 | //Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
|
---|
767 | ConcavityPerLine -= line->CalculateConvexity();
|
---|
768 | }
|
---|
769 | ConcavityPerLine /= point->lines.size();
|
---|
770 |
|
---|
771 | // weigh with total area of the surrounding triangles
|
---|
772 | totalarea = 0.;
|
---|
773 | TriangleSet *triangles = TesselStruct->GetAllTriangles(PointRunner->second);
|
---|
774 | for (TriangleSet::iterator TriangleRunner = triangles->begin(); TriangleRunner != triangles->end(); ++TriangleRunner) {
|
---|
775 | totalarea += CalculateAreaofGeneralTriangle((*TriangleRunner)->endpoints[0]->node->getPosition() , (*TriangleRunner)->endpoints[1]->node->getPosition() , (*TriangleRunner)->endpoints[2]->node->getPosition());
|
---|
776 | }
|
---|
777 | ConcavityPerLine *= totalarea;
|
---|
778 |
|
---|
779 | // calculate mean concavity over all attached triangles
|
---|
780 | ConcavityPerTriangle = 0.;
|
---|
781 | for (TriangleSet::const_iterator TriangleRunner = triangles->begin(); TriangleRunner != triangles->end(); ++TriangleRunner) {
|
---|
782 | line = (*TriangleRunner)->GetThirdLine(PointRunner->second);
|
---|
783 | triangle = line->GetOtherTriangle(*TriangleRunner);
|
---|
784 | area = CalculateAreaofGeneralTriangle(triangle->endpoints[0]->node->getPosition() , triangle->endpoints[1]->node->getPosition() , triangle->endpoints[2]->node->getPosition());
|
---|
785 | area += CalculateAreaofGeneralTriangle((*TriangleRunner)->endpoints[0]->node->getPosition() , (*TriangleRunner)->endpoints[1]->node->getPosition() , (*TriangleRunner)->endpoints[2]->node->getPosition());
|
---|
786 | area *= -line->CalculateConvexity();
|
---|
787 | if (area > 0)
|
---|
788 | ConcavityPerTriangle += area;
|
---|
789 | // else
|
---|
790 | // ConcavityPerTriangle -= area;
|
---|
791 | }
|
---|
792 | ConcavityPerTriangle /= triangles->size()/totalarea;
|
---|
793 | delete(triangles);
|
---|
794 |
|
---|
795 | // add up
|
---|
796 | point->value = ConcavityPerLine + ConcavityPerTriangle;
|
---|
797 | }
|
---|
798 | };
|
---|
799 |
|
---|
800 |
|
---|
801 |
|
---|
802 | /** Calculates the concavity for each of the BoundaryPointSet's in a Tesselation.
|
---|
803 | * Sets BoundaryPointSet::value equal to the nearest distance to convex envelope.
|
---|
804 | * \param *out output stream for debugging
|
---|
805 | * \param *TesselStruct pointer to Tesselation structure
|
---|
806 | * \param *Convex pointer to convex Tesselation structure as reference
|
---|
807 | */
|
---|
808 | void CalculateConstrictionPerBoundaryPoint(const Tesselation * const TesselStruct, const Tesselation * const Convex)
|
---|
809 | {
|
---|
810 | Info FunctionInfo(__func__);
|
---|
811 | double distance = 0.;
|
---|
812 |
|
---|
813 | for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
|
---|
814 | DoeLog(1) && (eLog() << Verbose(1) << "INFO: Current point is " << * PointRunner->second << "." << endl);
|
---|
815 |
|
---|
816 | distance = 0.;
|
---|
817 | for (TriangleMap::const_iterator TriangleRunner = Convex->TrianglesOnBoundary.begin(); TriangleRunner != Convex->TrianglesOnBoundary.end(); TriangleRunner++) {
|
---|
818 | const double CurrentDistance = Convex->GetDistanceSquaredToTriangle(PointRunner->second->node->getPosition() , TriangleRunner->second);
|
---|
819 | if (CurrentDistance < distance)
|
---|
820 | distance = CurrentDistance;
|
---|
821 | }
|
---|
822 |
|
---|
823 | PointRunner->second->value = distance;
|
---|
824 | }
|
---|
825 | };
|
---|
826 |
|
---|
827 | /** Checks whether each BoundaryLineSet in the Tesselation has two triangles.
|
---|
828 | * \param *out output stream for debugging
|
---|
829 | * \param *TesselStruct
|
---|
830 | * \return true - all have exactly two triangles, false - some not, list is printed to screen
|
---|
831 | */
|
---|
832 | bool CheckListOfBaselines(const Tesselation * const TesselStruct)
|
---|
833 | {
|
---|
834 | Info FunctionInfo(__func__);
|
---|
835 | LineMap::const_iterator testline;
|
---|
836 | bool result = false;
|
---|
837 | int counter = 0;
|
---|
838 |
|
---|
839 | DoLog(1) && (Log() << Verbose(1) << "Check: List of Baselines with not two connected triangles:" << endl);
|
---|
840 | for (testline = TesselStruct->LinesOnBoundary.begin(); testline != TesselStruct->LinesOnBoundary.end(); testline++) {
|
---|
841 | if (testline->second->triangles.size() != 2) {
|
---|
842 | DoLog(2) && (Log() << Verbose(2) << *testline->second << "\t" << testline->second->triangles.size() << endl);
|
---|
843 | counter++;
|
---|
844 | }
|
---|
845 | }
|
---|
846 | if (counter == 0) {
|
---|
847 | DoLog(1) && (Log() << Verbose(1) << "None." << endl);
|
---|
848 | result = true;
|
---|
849 | }
|
---|
850 | return result;
|
---|
851 | }
|
---|
852 |
|
---|
853 | /** Counts the number of triangle pairs that contain the given polygon.
|
---|
854 | * \param *P polygon with endpoints to look for
|
---|
855 | * \param *T set of triangles to create pairs from containing \a *P
|
---|
856 | */
|
---|
857 | int CountTrianglePairContainingPolygon(const BoundaryPolygonSet * const P, const TriangleSet * const T)
|
---|
858 | {
|
---|
859 | Info FunctionInfo(__func__);
|
---|
860 | // check number of endpoints in *P
|
---|
861 | if (P->endpoints.size() != 4) {
|
---|
862 | DoeLog(1) && (eLog()<< Verbose(1) << "CountTrianglePairContainingPolygon works only on polygons with 4 nodes!" << endl);
|
---|
863 | return 0;
|
---|
864 | }
|
---|
865 |
|
---|
866 | // check number of triangles in *T
|
---|
867 | if (T->size() < 2) {
|
---|
868 | DoeLog(1) && (eLog()<< Verbose(1) << "Not enough triangles to have pairs!" << endl);
|
---|
869 | return 0;
|
---|
870 | }
|
---|
871 |
|
---|
872 | DoLog(0) && (Log() << Verbose(0) << "Polygon is " << *P << endl);
|
---|
873 | // create each pair, get the endpoints and check whether *P is contained.
|
---|
874 | int counter = 0;
|
---|
875 | PointSet Trianglenodes;
|
---|
876 | class BoundaryPolygonSet PairTrianglenodes;
|
---|
877 | for(TriangleSet::iterator Walker = T->begin(); Walker != T->end(); Walker++) {
|
---|
878 | for (int i=0;i<3;i++)
|
---|
879 | Trianglenodes.insert((*Walker)->endpoints[i]);
|
---|
880 |
|
---|
881 | for(TriangleSet::iterator PairWalker = Walker; PairWalker != T->end(); PairWalker++) {
|
---|
882 | if (Walker != PairWalker) { // skip first
|
---|
883 | PairTrianglenodes.endpoints = Trianglenodes;
|
---|
884 | for (int i=0;i<3;i++)
|
---|
885 | PairTrianglenodes.endpoints.insert((*PairWalker)->endpoints[i]);
|
---|
886 | const int size = PairTrianglenodes.endpoints.size();
|
---|
887 | if (size == 4) {
|
---|
888 | DoLog(0) && (Log() << Verbose(0) << " Current pair of triangles: " << **Walker << "," << **PairWalker << " with " << size << " distinct endpoints:" << PairTrianglenodes << endl);
|
---|
889 | // now check
|
---|
890 | if (PairTrianglenodes.ContainsPresentTupel(P)) {
|
---|
891 | counter++;
|
---|
892 | DoLog(0) && (Log() << Verbose(0) << " ACCEPT: Matches with " << *P << endl);
|
---|
893 | } else {
|
---|
894 | DoLog(0) && (Log() << Verbose(0) << " REJECT: No match with " << *P << endl);
|
---|
895 | }
|
---|
896 | } else {
|
---|
897 | DoLog(0) && (Log() << Verbose(0) << " REJECT: Less than four endpoints." << endl);
|
---|
898 | }
|
---|
899 | }
|
---|
900 | }
|
---|
901 | Trianglenodes.clear();
|
---|
902 | }
|
---|
903 | return counter;
|
---|
904 | };
|
---|
905 |
|
---|
906 | /** Checks whether two give polygons have two or more points in common.
|
---|
907 | * \param *P1 first polygon
|
---|
908 | * \param *P2 second polygon
|
---|
909 | * \return true - are connected, false = are note
|
---|
910 | */
|
---|
911 | bool ArePolygonsEdgeConnected(const BoundaryPolygonSet * const P1, const BoundaryPolygonSet * const P2)
|
---|
912 | {
|
---|
913 | Info FunctionInfo(__func__);
|
---|
914 | int counter = 0;
|
---|
915 | for(PointSet::const_iterator Runner = P1->endpoints.begin(); Runner != P1->endpoints.end(); Runner++) {
|
---|
916 | if (P2->ContainsBoundaryPoint((*Runner))) {
|
---|
917 | counter++;
|
---|
918 | DoLog(1) && (Log() << Verbose(1) << *(*Runner) << " of second polygon is found in the first one." << endl);
|
---|
919 | return true;
|
---|
920 | }
|
---|
921 | }
|
---|
922 | return false;
|
---|
923 | };
|
---|
924 |
|
---|
925 | /** Combines second into the first and deletes the second.
|
---|
926 | * \param *P1 first polygon, contains all nodes on return
|
---|
927 | * \param *&P2 second polygon, is deleted.
|
---|
928 | */
|
---|
929 | void CombinePolygons(BoundaryPolygonSet * const P1, BoundaryPolygonSet * &P2)
|
---|
930 | {
|
---|
931 | Info FunctionInfo(__func__);
|
---|
932 | pair <PointSet::iterator, bool> Tester;
|
---|
933 | for(PointSet::iterator Runner = P2->endpoints.begin(); Runner != P2->endpoints.end(); Runner++) {
|
---|
934 | Tester = P1->endpoints.insert((*Runner));
|
---|
935 | if (Tester.second)
|
---|
936 | DoLog(0) && (Log() << Verbose(0) << "Inserting endpoint " << *(*Runner) << " into first polygon." << endl);
|
---|
937 | }
|
---|
938 | P2->endpoints.clear();
|
---|
939 | delete(P2);
|
---|
940 | };
|
---|
941 |
|
---|