source: src/tesselation.cpp@ a0064e

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since a0064e was a0064e, checked in by Frederik Heber <heber@…>, 14 years ago

Moved stuff in src/Helpers and src/Patterns out to stand-alone project CodePatterns.

Makefile.am's:

  • CodePatterns added to AM_LDFLAGS and AM_CFLAGS
  • libMoleCuilderHelpers removed

Helpers/...

  • defs.hpp include prefixed with Helpers/
  • helpers.?pp removed lots of old, unused functions: bound, ask_value, ...
  • fast_functions.hpp has lots of functions removed as well.
  • all other files and unit tests moved to project CodePatterns.

Patterns/...

  • added ax_codepatterns.m4 containing AM_PATH_CODEPATTERNS to configure.ac
  • Property mode set to 100644
File size: 189.3 KB
Line 
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 * tesselation.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#include <iomanip>
24
25#include "Helpers/helpers.hpp"
26#include "Helpers/Info.hpp"
27#include "BoundaryPointSet.hpp"
28#include "BoundaryLineSet.hpp"
29#include "BoundaryTriangleSet.hpp"
30#include "BoundaryPolygonSet.hpp"
31#include "TesselPoint.hpp"
32#include "CandidateForTesselation.hpp"
33#include "PointCloud.hpp"
34#include "linkedcell.hpp"
35#include "Helpers/Log.hpp"
36#include "tesselation.hpp"
37#include "tesselationhelpers.hpp"
38#include "triangleintersectionlist.hpp"
39#include "LinearAlgebra/Vector.hpp"
40#include "LinearAlgebra/Line.hpp"
41#include "LinearAlgebra/vector_ops.hpp"
42#include "Helpers/Verbose.hpp"
43#include "LinearAlgebra/Plane.hpp"
44#include "Exceptions/LinearDependenceException.hpp"
45#include "Helpers/Assert.hpp"
46
47#include "Helpers/fast_functions.hpp"
48
49class molecule;
50
51const char *TecplotSuffix=".dat";
52const char *Raster3DSuffix=".r3d";
53const char *VRMLSUffix=".wrl";
54
55const double ParallelEpsilon=1e-3;
56const double Tesselation::HULLEPSILON = 1e-9;
57
58/** Constructor of class Tesselation.
59 */
60Tesselation::Tesselation() :
61 PointsOnBoundaryCount(0),
62 LinesOnBoundaryCount(0),
63 TrianglesOnBoundaryCount(0),
64 LastTriangle(NULL),
65 TriangleFilesWritten(0),
66 InternalPointer(PointsOnBoundary.begin())
67{
68 Info FunctionInfo(__func__);
69}
70;
71
72/** Destructor of class Tesselation.
73 * We have to free all points, lines and triangles.
74 */
75Tesselation::~Tesselation()
76{
77 Info FunctionInfo(__func__);
78 DoLog(0) && (Log() << Verbose(0) << "Free'ing TesselStruct ... " << endl);
79 for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) {
80 if (runner->second != NULL) {
81 delete (runner->second);
82 runner->second = NULL;
83 } else
84 DoeLog(1) && (eLog() << Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl);
85 }
86 DoLog(0) && (Log() << Verbose(0) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl);
87}
88;
89
90TesselPoint * Tesselation::getValue(const_iterator &rhs) const
91{
92 return (*rhs).second->node;
93}
94
95TesselPoint * Tesselation::getValue(iterator &rhs) const
96{
97 return (*rhs).second->node;
98}
99
100/** PointCloud implementation of GetCenter
101 * Uses PointsOnBoundary and STL stuff.
102 */
103Vector * Tesselation::GetCenter(ofstream *out) const
104{
105 Info FunctionInfo(__func__);
106 Vector *Center = new Vector(0., 0., 0.);
107 int num = 0;
108 for (GoToFirst(); (!IsEnd()); GoToNext()) {
109 (*Center) += (GetPoint()->getPosition());
110 num++;
111 }
112 Center->Scale(1. / num);
113 return Center;
114}
115;
116
117/** PointCloud implementation of GoPoint
118 * Uses PointsOnBoundary and STL stuff.
119 */
120TesselPoint * Tesselation::GetPoint() const
121{
122 Info FunctionInfo(__func__);
123 return (InternalPointer->second->node);
124}
125;
126
127/** PointCloud implementation of GoToNext.
128 * Uses PointsOnBoundary and STL stuff.
129 */
130void Tesselation::GoToNext() const
131{
132 Info FunctionInfo(__func__);
133 if (InternalPointer != PointsOnBoundary.end())
134 InternalPointer++;
135}
136;
137
138/** PointCloud implementation of GoToFirst.
139 * Uses PointsOnBoundary and STL stuff.
140 */
141void Tesselation::GoToFirst() const
142{
143 Info FunctionInfo(__func__);
144 InternalPointer = PointsOnBoundary.begin();
145}
146;
147
148/** PointCloud implementation of IsEmpty.
149 * Uses PointsOnBoundary and STL stuff.
150 */
151bool Tesselation::IsEmpty() const
152{
153 Info FunctionInfo(__func__);
154 return (PointsOnBoundary.empty());
155}
156;
157
158/** PointCloud implementation of IsLast.
159 * Uses PointsOnBoundary and STL stuff.
160 */
161bool Tesselation::IsEnd() const
162{
163 Info FunctionInfo(__func__);
164 return (InternalPointer == PointsOnBoundary.end());
165}
166;
167
168/** Gueses first starting triangle of the convex envelope.
169 * We guess the starting triangle by taking the smallest distance between two points and looking for a fitting third.
170 * \param *out output stream for debugging
171 * \param PointsOnBoundary set of boundary points defining the convex envelope of the cluster
172 */
173void Tesselation::GuessStartingTriangle()
174{
175 Info FunctionInfo(__func__);
176 // 4b. create a starting triangle
177 // 4b1. create all distances
178 DistanceMultiMap DistanceMMap;
179 double distance, tmp;
180 Vector PlaneVector, TrialVector;
181 PointMap::iterator A, B, C; // three nodes of the first triangle
182 A = PointsOnBoundary.begin(); // the first may be chosen arbitrarily
183
184 // with A chosen, take each pair B,C and sort
185 if (A != PointsOnBoundary.end()) {
186 B = A;
187 B++;
188 for (; B != PointsOnBoundary.end(); B++) {
189 C = B;
190 C++;
191 for (; C != PointsOnBoundary.end(); C++) {
192 tmp = A->second->node->DistanceSquared(B->second->node->getPosition());
193 distance = tmp * tmp;
194 tmp = A->second->node->DistanceSquared(C->second->node->getPosition());
195 distance += tmp * tmp;
196 tmp = B->second->node->DistanceSquared(C->second->node->getPosition());
197 distance += tmp * tmp;
198 DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C)));
199 }
200 }
201 }
202 // // listing distances
203 // Log() << Verbose(1) << "Listing DistanceMMap:";
204 // for(DistanceMultiMap::iterator runner = DistanceMMap.begin(); runner != DistanceMMap.end(); runner++) {
205 // Log() << Verbose(0) << " " << runner->first << "(" << *runner->second.first->second << ", " << *runner->second.second->second << ")";
206 // }
207 // Log() << Verbose(0) << endl;
208 // 4b2. pick three baselines forming a triangle
209 // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
210 DistanceMultiMap::iterator baseline = DistanceMMap.begin();
211 for (; baseline != DistanceMMap.end(); baseline++) {
212 // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate
213 // 2. next, we have to check whether all points reside on only one side of the triangle
214 // 3. construct plane vector
215 PlaneVector = Plane(A->second->node->getPosition(),
216 baseline->second.first->second->node->getPosition(),
217 baseline->second.second->second->node->getPosition()).getNormal();
218 DoLog(2) && (Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl);
219 // 4. loop over all points
220 double sign = 0.;
221 PointMap::iterator checker = PointsOnBoundary.begin();
222 for (; checker != PointsOnBoundary.end(); checker++) {
223 // (neglecting A,B,C)
224 if ((checker == A) || (checker == baseline->second.first) || (checker == baseline->second.second))
225 continue;
226 // 4a. project onto plane vector
227 TrialVector = (checker->second->node->getPosition() - A->second->node->getPosition());
228 distance = TrialVector.ScalarProduct(PlaneVector);
229 if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok
230 continue;
231 DoLog(2) && (Log() << Verbose(2) << "Projection of " << checker->second->node->getName() << " yields distance of " << distance << "." << endl);
232 tmp = distance / fabs(distance);
233 // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle)
234 if ((sign != 0) && (tmp != sign)) {
235 // 4c. If so, break 4. loop and continue with next candidate in 1. loop
236 DoLog(2) && (Log() << Verbose(2) << "Current candidates: " << A->second->node->getName() << "," << baseline->second.first->second->node->getName() << "," << baseline->second.second->second->node->getName() << " leaves " << checker->second->node->getName() << " outside the convex hull." << endl);
237 break;
238 } else { // note the sign for later
239 DoLog(2) && (Log() << Verbose(2) << "Current candidates: " << A->second->node->getName() << "," << baseline->second.first->second->node->getName() << "," << baseline->second.second->second->node->getName() << " leave " << checker->second->node->getName() << " inside the convex hull." << endl);
240 sign = tmp;
241 }
242 // 4d. Check whether the point is inside the triangle (check distance to each node
243 tmp = checker->second->node->DistanceSquared(A->second->node->getPosition());
244 int innerpoint = 0;
245 if ((tmp < A->second->node->DistanceSquared(baseline->second.first->second->node->getPosition())) && (tmp < A->second->node->DistanceSquared(baseline->second.second->second->node->getPosition())))
246 innerpoint++;
247 tmp = checker->second->node->DistanceSquared(baseline->second.first->second->node->getPosition());
248 if ((tmp < baseline->second.first->second->node->DistanceSquared(A->second->node->getPosition())) && (tmp < baseline->second.first->second->node->DistanceSquared(baseline->second.second->second->node->getPosition())))
249 innerpoint++;
250 tmp = checker->second->node->DistanceSquared(baseline->second.second->second->node->getPosition());
251 if ((tmp < baseline->second.second->second->node->DistanceSquared(baseline->second.first->second->node->getPosition())) && (tmp < baseline->second.second->second->node->DistanceSquared(A->second->node->getPosition())))
252 innerpoint++;
253 // 4e. If so, break 4. loop and continue with next candidate in 1. loop
254 if (innerpoint == 3)
255 break;
256 }
257 // 5. come this far, all on same side? Then break 1. loop and construct triangle
258 if (checker == PointsOnBoundary.end()) {
259 DoLog(2) && (Log() << Verbose(2) << "Looks like we have a candidate!" << endl);
260 break;
261 }
262 }
263 if (baseline != DistanceMMap.end()) {
264 BPS[0] = baseline->second.first->second;
265 BPS[1] = baseline->second.second->second;
266 BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
267 BPS[0] = A->second;
268 BPS[1] = baseline->second.second->second;
269 BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
270 BPS[0] = baseline->second.first->second;
271 BPS[1] = A->second;
272 BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
273
274 // 4b3. insert created triangle
275 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
276 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
277 TrianglesOnBoundaryCount++;
278 for (int i = 0; i < NDIM; i++) {
279 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i]));
280 LinesOnBoundaryCount++;
281 }
282
283 DoLog(1) && (Log() << Verbose(1) << "Starting triangle is " << *BTS << "." << endl);
284 } else {
285 DoeLog(0) && (eLog() << Verbose(0) << "No starting triangle found." << endl);
286 }
287}
288;
289
290/** Tesselates the convex envelope of a cluster from a single starting triangle.
291 * The starting triangle is made out of three baselines. Each line in the final tesselated cluster may belong to at most
292 * 2 triangles. Hence, we go through all current lines:
293 * -# if the lines contains to only one triangle
294 * -# We search all points in the boundary
295 * -# if the triangle is in forward direction of the baseline (at most 90 degrees angle between vector orthogonal to
296 * baseline in triangle plane pointing out of the triangle and normal vector of new triangle)
297 * -# if the triangle with the baseline and the current point has the smallest of angles (comparison between normal vectors)
298 * -# then we have a new triangle, whose baselines we again add (or increase their TriangleCount)
299 * \param *out output stream for debugging
300 * \param *configuration for IsAngstroem
301 * \param *cloud cluster of points
302 */
303void Tesselation::TesselateOnBoundary(const PointCloud * const cloud)
304{
305 Info FunctionInfo(__func__);
306 bool flag;
307 PointMap::iterator winner;
308 class BoundaryPointSet *peak = NULL;
309 double SmallestAngle, TempAngle;
310 Vector NormalVector, VirtualNormalVector, CenterVector, TempVector, helper, PropagationVector, *Center = NULL;
311 LineMap::iterator LineChecker[2];
312
313 Center = cloud->GetCenter();
314 // create a first tesselation with the given BoundaryPoints
315 do {
316 flag = false;
317 for (LineMap::iterator baseline = LinesOnBoundary.begin(); baseline != LinesOnBoundary.end(); baseline++)
318 if (baseline->second->triangles.size() == 1) {
319 // 5a. go through each boundary point if not _both_ edges between either endpoint of the current line and this point exist (and belong to 2 triangles)
320 SmallestAngle = M_PI;
321
322 // get peak point with respect to this base line's only triangle
323 BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far
324 DoLog(0) && (Log() << Verbose(0) << "Current baseline is between " << *(baseline->second) << "." << endl);
325 for (int i = 0; i < 3; i++)
326 if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1]))
327 peak = BTS->endpoints[i];
328 DoLog(1) && (Log() << Verbose(1) << " and has peak " << *peak << "." << endl);
329
330 // prepare some auxiliary vectors
331 Vector BaseLineCenter, BaseLine;
332 BaseLineCenter = 0.5 * ((baseline->second->endpoints[0]->node->getPosition()) +
333 (baseline->second->endpoints[1]->node->getPosition()));
334 BaseLine = (baseline->second->endpoints[0]->node->getPosition()) - (baseline->second->endpoints[1]->node->getPosition());
335
336 // offset to center of triangle
337 CenterVector.Zero();
338 for (int i = 0; i < 3; i++)
339 CenterVector += BTS->getEndpoint(i);
340 CenterVector.Scale(1. / 3.);
341 DoLog(2) && (Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl);
342
343 // normal vector of triangle
344 NormalVector = (*Center) - CenterVector;
345 BTS->GetNormalVector(NormalVector);
346 NormalVector = BTS->NormalVector;
347 DoLog(2) && (Log() << Verbose(2) << "NormalVector of base triangle is " << NormalVector << endl);
348
349 // vector in propagation direction (out of triangle)
350 // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection)
351 PropagationVector = Plane(BaseLine, NormalVector,0).getNormal();
352 TempVector = CenterVector - (baseline->second->endpoints[0]->node->getPosition()); // TempVector is vector on triangle plane pointing from one baseline egde towards center!
353 //Log() << Verbose(0) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;
354 if (PropagationVector.ScalarProduct(TempVector) > 0) // make sure normal propagation vector points outward from baseline
355 PropagationVector.Scale(-1.);
356 DoLog(2) && (Log() << Verbose(2) << "PropagationVector of base triangle is " << PropagationVector << endl);
357 winner = PointsOnBoundary.end();
358
359 // loop over all points and calculate angle between normal vector of new and present triangle
360 for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) {
361 if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints
362 DoLog(1) && (Log() << Verbose(1) << "Target point is " << *(target->second) << ":" << endl);
363
364 // first check direction, so that triangles don't intersect
365 VirtualNormalVector = (target->second->node->getPosition()) - BaseLineCenter;
366 VirtualNormalVector.ProjectOntoPlane(NormalVector);
367 TempAngle = VirtualNormalVector.Angle(PropagationVector);
368 DoLog(2) && (Log() << Verbose(2) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl);
369 if (TempAngle > (M_PI / 2.)) { // no bends bigger than Pi/2 (90 degrees)
370 DoLog(2) && (Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl);
371 continue;
372 } else
373 DoLog(2) && (Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl);
374
375 // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle)
376 LineChecker[0] = baseline->second->endpoints[0]->lines.find(target->first);
377 LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first);
378 if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->triangles.size() == 2))) {
379 DoLog(2) && (Log() << Verbose(2) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->triangles.size() << " triangles." << endl);
380 continue;
381 }
382 if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->triangles.size() == 2))) {
383 DoLog(2) && (Log() << Verbose(2) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->triangles.size() << " triangles." << endl);
384 continue;
385 }
386
387 // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint)
388 if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) {
389 DoLog(4) && (Log() << Verbose(4) << "Current target is peak!" << endl);
390 continue;
391 }
392
393 // check for linear dependence
394 TempVector = (baseline->second->endpoints[0]->node->getPosition()) - (target->second->node->getPosition());
395 helper = (baseline->second->endpoints[1]->node->getPosition()) - (target->second->node->getPosition());
396 helper.ProjectOntoPlane(TempVector);
397 if (fabs(helper.NormSquared()) < MYEPSILON) {
398 DoLog(2) && (Log() << Verbose(2) << "Chosen set of vectors is linear dependent." << endl);
399 continue;
400 }
401
402 // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle
403 flag = true;
404 VirtualNormalVector = Plane((baseline->second->endpoints[0]->node->getPosition()),
405 (baseline->second->endpoints[1]->node->getPosition()),
406 (target->second->node->getPosition())).getNormal();
407 TempVector = (1./3.) * ((baseline->second->endpoints[0]->node->getPosition()) +
408 (baseline->second->endpoints[1]->node->getPosition()) +
409 (target->second->node->getPosition()));
410 TempVector -= (*Center);
411 // make it always point outward
412 if (VirtualNormalVector.ScalarProduct(TempVector) < 0)
413 VirtualNormalVector.Scale(-1.);
414 // calculate angle
415 TempAngle = NormalVector.Angle(VirtualNormalVector);
416 DoLog(2) && (Log() << Verbose(2) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl);
417 if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner
418 SmallestAngle = TempAngle;
419 winner = target;
420 DoLog(2) && (Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl);
421 } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle)
422 // hence, check the angles to some normal direction from our base line but in this common plane of both targets...
423 helper = (target->second->node->getPosition()) - BaseLineCenter;
424 helper.ProjectOntoPlane(BaseLine);
425 // ...the one with the smaller angle is the better candidate
426 TempVector = (target->second->node->getPosition()) - BaseLineCenter;
427 TempVector.ProjectOntoPlane(VirtualNormalVector);
428 TempAngle = TempVector.Angle(helper);
429 TempVector = (winner->second->node->getPosition()) - BaseLineCenter;
430 TempVector.ProjectOntoPlane(VirtualNormalVector);
431 if (TempAngle < TempVector.Angle(helper)) {
432 TempAngle = NormalVector.Angle(VirtualNormalVector);
433 SmallestAngle = TempAngle;
434 winner = target;
435 DoLog(2) && (Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl);
436 } else
437 DoLog(2) && (Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl);
438 } else
439 DoLog(2) && (Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl);
440 }
441 } // end of loop over all boundary points
442
443 // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle
444 if (winner != PointsOnBoundary.end()) {
445 DoLog(0) && (Log() << Verbose(0) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl);
446 // create the lins of not yet present
447 BLS[0] = baseline->second;
448 // 5c. add lines to the line set if those were new (not yet part of a triangle), delete lines that belong to two triangles)
449 LineChecker[0] = baseline->second->endpoints[0]->lines.find(winner->first);
450 LineChecker[1] = baseline->second->endpoints[1]->lines.find(winner->first);
451 if (LineChecker[0] == baseline->second->endpoints[0]->lines.end()) { // create
452 BPS[0] = baseline->second->endpoints[0];
453 BPS[1] = winner->second;
454 BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
455 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[1]));
456 LinesOnBoundaryCount++;
457 } else
458 BLS[1] = LineChecker[0]->second;
459 if (LineChecker[1] == baseline->second->endpoints[1]->lines.end()) { // create
460 BPS[0] = baseline->second->endpoints[1];
461 BPS[1] = winner->second;
462 BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
463 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[2]));
464 LinesOnBoundaryCount++;
465 } else
466 BLS[2] = LineChecker[1]->second;
467 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
468 BTS->GetCenter(helper);
469 helper -= (*Center);
470 helper *= -1;
471 BTS->GetNormalVector(helper);
472 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
473 TrianglesOnBoundaryCount++;
474 } else {
475 DoeLog(2) && (eLog() << Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl);
476 }
477
478 // 5d. If the set of lines is not yet empty, go to 5. and continue
479 } else
480 DoLog(0) && (Log() << Verbose(0) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl);
481 } while (flag);
482
483 // exit
484 delete (Center);
485}
486;
487
488/** Inserts all points outside of the tesselated surface into it by adding new triangles.
489 * \param *out output stream for debugging
490 * \param *cloud cluster of points
491 * \param *LC LinkedCell structure to find nearest point quickly
492 * \return true - all straddling points insert, false - something went wrong
493 */
494bool Tesselation::InsertStraddlingPoints(const PointCloud *cloud, const LinkedCell *LC)
495{
496 Info FunctionInfo(__func__);
497 Vector Intersection, Normal;
498 TesselPoint *Walker = NULL;
499 Vector *Center = cloud->GetCenter();
500 TriangleList *triangles = NULL;
501 bool AddFlag = false;
502 LinkedCell *BoundaryPoints = NULL;
503 bool SuccessFlag = true;
504
505 cloud->GoToFirst();
506 BoundaryPoints = new LinkedCell(*this, 5.);
507 while (!cloud->IsEnd()) { // we only have to go once through all points, as boundary can become only bigger
508 if (AddFlag) {
509 delete (BoundaryPoints);
510 BoundaryPoints = new LinkedCell(*this, 5.);
511 AddFlag = false;
512 }
513 Walker = cloud->GetPoint();
514 DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Walker << "." << endl);
515 // get the next triangle
516 triangles = FindClosestTrianglesToVector(Walker->getPosition(), BoundaryPoints);
517 if (triangles != NULL)
518 BTS = triangles->front();
519 else
520 BTS = NULL;
521 delete triangles;
522 if ((BTS == NULL) || (BTS->ContainsBoundaryPoint(Walker))) {
523 DoLog(0) && (Log() << Verbose(0) << "No triangles found, probably a tesselation point itself." << endl);
524 cloud->GoToNext();
525 continue;
526 } else {
527 }
528 DoLog(0) && (Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl);
529 // get the intersection point
530 if (BTS->GetIntersectionInsideTriangle(*Center, Walker->getPosition(), Intersection)) {
531 DoLog(0) && (Log() << Verbose(0) << "We have an intersection at " << Intersection << "." << endl);
532 // we have the intersection, check whether in- or outside of boundary
533 if ((Center->DistanceSquared(Walker->getPosition()) - Center->DistanceSquared(Intersection)) < -MYEPSILON) {
534 // inside, next!
535 DoLog(0) && (Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl);
536 } else {
537 // outside!
538 DoLog(0) && (Log() << Verbose(0) << *Walker << " is outside wrt triangle " << *BTS << "." << endl);
539 class BoundaryLineSet *OldLines[3], *NewLines[3];
540 class BoundaryPointSet *OldPoints[3], *NewPoint;
541 // store the three old lines and old points
542 for (int i = 0; i < 3; i++) {
543 OldLines[i] = BTS->lines[i];
544 OldPoints[i] = BTS->endpoints[i];
545 }
546 Normal = BTS->NormalVector;
547 // add Walker to boundary points
548 DoLog(0) && (Log() << Verbose(0) << "Adding " << *Walker << " to BoundaryPoints." << endl);
549 AddFlag = true;
550 if (AddBoundaryPoint(Walker, 0))
551 NewPoint = BPS[0];
552 else
553 continue;
554 // remove triangle
555 DoLog(0) && (Log() << Verbose(0) << "Erasing triangle " << *BTS << "." << endl);
556 TrianglesOnBoundary.erase(BTS->Nr);
557 delete (BTS);
558 // create three new boundary lines
559 for (int i = 0; i < 3; i++) {
560 BPS[0] = NewPoint;
561 BPS[1] = OldPoints[i];
562 NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount);
563 DoLog(1) && (Log() << Verbose(1) << "Creating new line " << *NewLines[i] << "." << endl);
564 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one
565 LinesOnBoundaryCount++;
566 }
567 // create three new triangle with new point
568 for (int i = 0; i < 3; i++) { // find all baselines
569 BLS[0] = OldLines[i];
570 int n = 1;
571 for (int j = 0; j < 3; j++) {
572 if (NewLines[j]->IsConnectedTo(BLS[0])) {
573 if (n > 2) {
574 DoeLog(2) && (eLog() << Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl);
575 return false;
576 } else
577 BLS[n++] = NewLines[j];
578 }
579 }
580 // create the triangle
581 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
582 Normal.Scale(-1.);
583 BTS->GetNormalVector(Normal);
584 Normal.Scale(-1.);
585 DoLog(0) && (Log() << Verbose(0) << "Created new triangle " << *BTS << "." << endl);
586 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
587 TrianglesOnBoundaryCount++;
588 }
589 }
590 } else { // something is wrong with FindClosestTriangleToPoint!
591 DoeLog(1) && (eLog() << Verbose(1) << "The closest triangle did not produce an intersection!" << endl);
592 SuccessFlag = false;
593 break;
594 }
595 cloud->GoToNext();
596 }
597
598 // exit
599 delete (Center);
600 delete (BoundaryPoints);
601 return SuccessFlag;
602}
603;
604
605/** Adds a point to the tesselation::PointsOnBoundary list.
606 * \param *Walker point to add
607 * \param n TesselStruct::BPS index to put pointer into
608 * \return true - new point was added, false - point already present
609 */
610bool Tesselation::AddBoundaryPoint(TesselPoint * Walker, const int n)
611{
612 Info FunctionInfo(__func__);
613 PointTestPair InsertUnique;
614 BPS[n] = new class BoundaryPointSet(Walker);
615 InsertUnique = PointsOnBoundary.insert(PointPair(Walker->nr, BPS[n]));
616 if (InsertUnique.second) { // if new point was not present before, increase counter
617 PointsOnBoundaryCount++;
618 return true;
619 } else {
620 delete (BPS[n]);
621 BPS[n] = InsertUnique.first->second;
622 return false;
623 }
624}
625;
626
627/** Adds point to Tesselation::PointsOnBoundary if not yet present.
628 * Tesselation::TPS is set to either this new BoundaryPointSet or to the existing one of not unique.
629 * @param Candidate point to add
630 * @param n index for this point in Tesselation::TPS array
631 */
632void Tesselation::AddTesselationPoint(TesselPoint* Candidate, const int n)
633{
634 Info FunctionInfo(__func__);
635 PointTestPair InsertUnique;
636 TPS[n] = new class BoundaryPointSet(Candidate);
637 InsertUnique = PointsOnBoundary.insert(PointPair(Candidate->nr, TPS[n]));
638 if (InsertUnique.second) { // if new point was not present before, increase counter
639 PointsOnBoundaryCount++;
640 } else {
641 delete TPS[n];
642 DoLog(0) && (Log() << Verbose(0) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl);
643 TPS[n] = (InsertUnique.first)->second;
644 }
645}
646;
647
648/** Sets point to a present Tesselation::PointsOnBoundary.
649 * Tesselation::TPS is set to the existing one or NULL if not found.
650 * @param Candidate point to set to
651 * @param n index for this point in Tesselation::TPS array
652 */
653void Tesselation::SetTesselationPoint(TesselPoint* Candidate, const int n) const
654{
655 Info FunctionInfo(__func__);
656 PointMap::const_iterator FindPoint = PointsOnBoundary.find(Candidate->nr);
657 if (FindPoint != PointsOnBoundary.end())
658 TPS[n] = FindPoint->second;
659 else
660 TPS[n] = NULL;
661}
662;
663
664/** Function tries to add line from current Points in BPS to BoundaryLineSet.
665 * If successful it raises the line count and inserts the new line into the BLS,
666 * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one.
667 * @param *OptCenter desired OptCenter if there are more than one candidate line
668 * @param *candidate third point of the triangle to be, for checking between multiple open line candidates
669 * @param *a first endpoint
670 * @param *b second endpoint
671 * @param n index of Tesselation::BLS giving the line with both endpoints
672 */
673void Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
674{
675 bool insertNewLine = true;
676 LineMap::iterator FindLine = a->lines.find(b->node->nr);
677 BoundaryLineSet *WinningLine = NULL;
678 if (FindLine != a->lines.end()) {
679 DoLog(1) && (Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl);
680
681 pair<LineMap::iterator, LineMap::iterator> FindPair;
682 FindPair = a->lines.equal_range(b->node->nr);
683
684 for (FindLine = FindPair.first; (FindLine != FindPair.second) && (insertNewLine); FindLine++) {
685 DoLog(1) && (Log() << Verbose(1) << "INFO: Checking line " << *(FindLine->second) << " ..." << endl);
686 // If there is a line with less than two attached triangles, we don't need a new line.
687 if (FindLine->second->triangles.size() == 1) {
688 CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
689 if (!Finder->second->pointlist.empty())
690 DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << "." << endl);
691 else
692 DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate." << endl);
693 // get open line
694 for (TesselPointList::const_iterator CandidateChecker = Finder->second->pointlist.begin(); CandidateChecker != Finder->second->pointlist.end(); ++CandidateChecker) {
695 if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(Finder->second->OptCenter) < MYEPSILON )) { // stop searching if candidate matches
696 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Candidate " << *(*CandidateChecker) << " has the right center " << Finder->second->OptCenter << "." << endl);
697 insertNewLine = false;
698 WinningLine = FindLine->second;
699 break;
700 } else {
701 DoLog(1) && (Log() << Verbose(1) << "REJECT: Candidate " << *(*CandidateChecker) << "'s center " << Finder->second->OptCenter << " does not match desired on " << *OptCenter << "." << endl);
702 }
703 }
704 }
705 }
706 }
707
708 if (insertNewLine) {
709 AddNewTesselationTriangleLine(a, b, n);
710 } else {
711 AddExistingTesselationTriangleLine(WinningLine, n);
712 }
713}
714;
715
716/**
717 * Adds lines from each of the current points in the BPS to BoundaryLineSet.
718 * Raises the line count and inserts the new line into the BLS.
719 *
720 * @param *a first endpoint
721 * @param *b second endpoint
722 * @param n index of Tesselation::BLS giving the line with both endpoints
723 */
724void Tesselation::AddNewTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
725{
726 Info FunctionInfo(__func__);
727 DoLog(0) && (Log() << Verbose(0) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl);
728 BPS[0] = a;
729 BPS[1] = b;
730 BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); // this also adds the line to the local maps
731 // add line to global map
732 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n]));
733 // increase counter
734 LinesOnBoundaryCount++;
735 // also add to open lines
736 CandidateForTesselation *CFT = new CandidateForTesselation(BLS[n]);
737 OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT));
738}
739;
740
741/** Uses an existing line for a new triangle.
742 * Sets Tesselation::BLS[\a n] and removes the lines from Tesselation::OpenLines.
743 * \param *FindLine the line to add
744 * \param n index of the line to set in Tesselation::BLS
745 */
746void Tesselation::AddExistingTesselationTriangleLine(class BoundaryLineSet *Line, int n)
747{
748 Info FunctionInfo(__func__);
749 DoLog(0) && (Log() << Verbose(0) << "Using existing line " << *Line << endl);
750
751 // set endpoints and line
752 BPS[0] = Line->endpoints[0];
753 BPS[1] = Line->endpoints[1];
754 BLS[n] = Line;
755 // remove existing line from OpenLines
756 CandidateMap::iterator CandidateLine = OpenLines.find(BLS[n]);
757 if (CandidateLine != OpenLines.end()) {
758 DoLog(1) && (Log() << Verbose(1) << " Removing line from OpenLines." << endl);
759 delete (CandidateLine->second);
760 OpenLines.erase(CandidateLine);
761 } else {
762 DoeLog(1) && (eLog() << Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl);
763 }
764}
765;
766
767/** Function adds triangle to global list.
768 * Furthermore, the triangle receives the next free id and id counter \a TrianglesOnBoundaryCount is increased.
769 */
770void Tesselation::AddTesselationTriangle()
771{
772 Info FunctionInfo(__func__);
773 DoLog(1) && (Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl);
774
775 // add triangle to global map
776 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS));
777 TrianglesOnBoundaryCount++;
778
779 // set as last new triangle
780 LastTriangle = BTS;
781
782 // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
783}
784;
785
786/** Function adds triangle to global list.
787 * Furthermore, the triangle number is set to \a nr.
788 * \param nr triangle number
789 */
790void Tesselation::AddTesselationTriangle(const int nr)
791{
792 Info FunctionInfo(__func__);
793 DoLog(0) && (Log() << Verbose(0) << "Adding triangle to global TrianglesOnBoundary map." << endl);
794
795 // add triangle to global map
796 TrianglesOnBoundary.insert(TrianglePair(nr, BTS));
797
798 // set as last new triangle
799 LastTriangle = BTS;
800
801 // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet
802}
803;
804
805/** Removes a triangle from the tesselation.
806 * Removes itself from the TriangleMap's of its lines, calls for them RemoveTriangleLine() if they are no more connected.
807 * Removes itself from memory.
808 * \param *triangle to remove
809 */
810void Tesselation::RemoveTesselationTriangle(class BoundaryTriangleSet *triangle)
811{
812 Info FunctionInfo(__func__);
813 if (triangle == NULL)
814 return;
815 for (int i = 0; i < 3; i++) {
816 if (triangle->lines[i] != NULL) {
817 DoLog(0) && (Log() << Verbose(0) << "Removing triangle Nr." << triangle->Nr << " in line " << *triangle->lines[i] << "." << endl);
818 triangle->lines[i]->triangles.erase(triangle->Nr);
819 if (triangle->lines[i]->triangles.empty()) {
820 DoLog(0) && (Log() << Verbose(0) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl);
821 RemoveTesselationLine(triangle->lines[i]);
822 } else {
823 DoLog(0) && (Log() << Verbose(0) << *triangle->lines[i] << " is still attached to another triangle: " << endl);
824 OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> (triangle->lines[i], NULL));
825 for (TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++)
826 DoLog(0) && (Log() << Verbose(0) << "\t[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t");
827 DoLog(0) && (Log() << Verbose(0) << endl);
828 // for (int j=0;j<2;j++) {
829 // Log() << Verbose(0) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": ";
830 // for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++)
831 // Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t";
832 // Log() << Verbose(0) << endl;
833 // }
834 }
835 triangle->lines[i] = NULL; // free'd or not: disconnect
836 } else
837 DoeLog(1) && (eLog() << Verbose(1) << "This line " << i << " has already been free'd." << endl);
838 }
839
840 if (TrianglesOnBoundary.erase(triangle->Nr))
841 DoLog(0) && (Log() << Verbose(0) << "Removing triangle Nr. " << triangle->Nr << "." << endl);
842 delete (triangle);
843}
844;
845
846/** Removes a line from the tesselation.
847 * Removes itself from each endpoints' LineMap, then removes itself from global LinesOnBoundary list and free's the line.
848 * \param *line line to remove
849 */
850void Tesselation::RemoveTesselationLine(class BoundaryLineSet *line)
851{
852 Info FunctionInfo(__func__);
853 int Numbers[2];
854
855 if (line == NULL)
856 return;
857 // get other endpoint number for finding copies of same line
858 if (line->endpoints[1] != NULL)
859 Numbers[0] = line->endpoints[1]->Nr;
860 else
861 Numbers[0] = -1;
862 if (line->endpoints[0] != NULL)
863 Numbers[1] = line->endpoints[0]->Nr;
864 else
865 Numbers[1] = -1;
866
867 for (int i = 0; i < 2; i++) {
868 if (line->endpoints[i] != NULL) {
869 if (Numbers[i] != -1) { // as there may be multiple lines with same endpoints, we have to go through each and find in the endpoint's line list this line set
870 pair<LineMap::iterator, LineMap::iterator> erasor = line->endpoints[i]->lines.equal_range(Numbers[i]);
871 for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++)
872 if ((*Runner).second == line) {
873 DoLog(0) && (Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl);
874 line->endpoints[i]->lines.erase(Runner);
875 break;
876 }
877 } else { // there's just a single line left
878 if (line->endpoints[i]->lines.erase(line->Nr))
879 DoLog(0) && (Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl);
880 }
881 if (line->endpoints[i]->lines.empty()) {
882 DoLog(0) && (Log() << Verbose(0) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl);
883 RemoveTesselationPoint(line->endpoints[i]);
884 } else {
885 DoLog(0) && (Log() << Verbose(0) << *line->endpoints[i] << " has still lines it's attached to: ");
886 for (LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++)
887 DoLog(0) && (Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t");
888 DoLog(0) && (Log() << Verbose(0) << endl);
889 }
890 line->endpoints[i] = NULL; // free'd or not: disconnect
891 } else
892 DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << i << " has already been free'd." << endl);
893 }
894 if (!line->triangles.empty())
895 DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *line << " am still connected to some triangles." << endl);
896
897 if (LinesOnBoundary.erase(line->Nr))
898 DoLog(0) && (Log() << Verbose(0) << "Removing line Nr. " << line->Nr << "." << endl);
899 delete (line);
900}
901;
902
903/** Removes a point from the tesselation.
904 * Checks whether there are still lines connected, removes from global PointsOnBoundary list, then free's the point.
905 * \note If a point should be removed, while keep the tesselated surface intact (i.e. closed), use RemovePointFromTesselatedSurface()
906 * \param *point point to remove
907 */
908void Tesselation::RemoveTesselationPoint(class BoundaryPointSet *point)
909{
910 Info FunctionInfo(__func__);
911 if (point == NULL)
912 return;
913 if (PointsOnBoundary.erase(point->Nr))
914 DoLog(0) && (Log() << Verbose(0) << "Removing point Nr. " << point->Nr << "." << endl);
915 delete (point);
916}
917;
918
919/** Checks validity of a given sphere of a candidate line.
920 * \sa CandidateForTesselation::CheckValidity(), which is more evolved.
921 * We check CandidateForTesselation::OtherOptCenter
922 * \param &CandidateLine contains other degenerated candidates which we have to subtract as well
923 * \param RADIUS radius of sphere
924 * \param *LC LinkedCell structure with other atoms
925 * \return true - candidate triangle is degenerated, false - candidate triangle is not degenerated
926 */
927bool Tesselation::CheckDegeneracy(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC) const
928{
929 Info FunctionInfo(__func__);
930
931 DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl);
932 bool flag = true;
933
934 DoLog(1) && (Log() << Verbose(1) << "Check by: draw sphere {" << CandidateLine.OtherOptCenter[0] << " " << CandidateLine.OtherOptCenter[1] << " " << CandidateLine.OtherOptCenter[2] << "} radius " << RADIUS << " resolution 30" << endl);
935 // get all points inside the sphere
936 TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, &CandidateLine.OtherOptCenter);
937
938 DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
939 for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
940 DoLog(1) && (Log() << Verbose(1) << " " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << "." << endl);
941
942 // remove triangles's endpoints
943 for (int i = 0; i < 2; i++)
944 ListofPoints->remove(CandidateLine.BaseLine->endpoints[i]->node);
945
946 // remove other candidates
947 for (TesselPointList::const_iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); ++Runner)
948 ListofPoints->remove(*Runner);
949
950 // check for other points
951 if (!ListofPoints->empty()) {
952 DoLog(1) && (Log() << Verbose(1) << "CheckDegeneracy: There are still " << ListofPoints->size() << " points inside the sphere." << endl);
953 flag = false;
954 DoLog(1) && (Log() << Verbose(1) << "External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":" << endl);
955 for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner)
956 DoLog(1) && (Log() << Verbose(1) << " " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << "." << endl);
957 }
958 delete (ListofPoints);
959
960 return flag;
961}
962;
963
964/** Checks whether the triangle consisting of the three points is already present.
965 * Searches for the points in Tesselation::PointsOnBoundary and checks their
966 * lines. If any of the three edges already has two triangles attached, false is
967 * returned.
968 * \param *out output stream for debugging
969 * \param *Candidates endpoints of the triangle candidate
970 * \return integer 0 if no triangle exists, 1 if one triangle exists, 2 if two
971 * triangles exist which is the maximum for three points
972 */
973int Tesselation::CheckPresenceOfTriangle(TesselPoint *Candidates[3]) const
974{
975 Info FunctionInfo(__func__);
976 int adjacentTriangleCount = 0;
977 class BoundaryPointSet *Points[3];
978
979 // builds a triangle point set (Points) of the end points
980 for (int i = 0; i < 3; i++) {
981 PointMap::const_iterator FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
982 if (FindPoint != PointsOnBoundary.end()) {
983 Points[i] = FindPoint->second;
984 } else {
985 Points[i] = NULL;
986 }
987 }
988
989 // checks lines between the points in the Points for their adjacent triangles
990 for (int i = 0; i < 3; i++) {
991 if (Points[i] != NULL) {
992 for (int j = i; j < 3; j++) {
993 if (Points[j] != NULL) {
994 LineMap::const_iterator FindLine = Points[i]->lines.find(Points[j]->node->nr);
995 for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) {
996 TriangleMap *triangles = &FindLine->second->triangles;
997 DoLog(1) && (Log() << Verbose(1) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl);
998 for (TriangleMap::const_iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) {
999 if (FindTriangle->second->IsPresentTupel(Points)) {
1000 adjacentTriangleCount++;
1001 }
1002 }
1003 DoLog(1) && (Log() << Verbose(1) << "end." << endl);
1004 }
1005 // Only one of the triangle lines must be considered for the triangle count.
1006 //Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
1007 //return adjacentTriangleCount;
1008 }
1009 }
1010 }
1011 }
1012
1013 DoLog(0) && (Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl);
1014 return adjacentTriangleCount;
1015}
1016;
1017
1018/** Checks whether the triangle consisting of the three points is already present.
1019 * Searches for the points in Tesselation::PointsOnBoundary and checks their
1020 * lines. If any of the three edges already has two triangles attached, false is
1021 * returned.
1022 * \param *out output stream for debugging
1023 * \param *Candidates endpoints of the triangle candidate
1024 * \return NULL - none found or pointer to triangle
1025 */
1026class BoundaryTriangleSet * Tesselation::GetPresentTriangle(TesselPoint *Candidates[3])
1027{
1028 Info FunctionInfo(__func__);
1029 class BoundaryTriangleSet *triangle = NULL;
1030 class BoundaryPointSet *Points[3];
1031
1032 // builds a triangle point set (Points) of the end points
1033 for (int i = 0; i < 3; i++) {
1034 PointMap::iterator FindPoint = PointsOnBoundary.find(Candidates[i]->nr);
1035 if (FindPoint != PointsOnBoundary.end()) {
1036 Points[i] = FindPoint->second;
1037 } else {
1038 Points[i] = NULL;
1039 }
1040 }
1041
1042 // checks lines between the points in the Points for their adjacent triangles
1043 for (int i = 0; i < 3; i++) {
1044 if (Points[i] != NULL) {
1045 for (int j = i; j < 3; j++) {
1046 if (Points[j] != NULL) {
1047 LineMap::iterator FindLine = Points[i]->lines.find(Points[j]->node->nr);
1048 for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) {
1049 TriangleMap *triangles = &FindLine->second->triangles;
1050 for (TriangleMap::iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) {
1051 if (FindTriangle->second->IsPresentTupel(Points)) {
1052 if ((triangle == NULL) || (triangle->Nr > FindTriangle->second->Nr))
1053 triangle = FindTriangle->second;
1054 }
1055 }
1056 }
1057 // Only one of the triangle lines must be considered for the triangle count.
1058 //Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;
1059 //return adjacentTriangleCount;
1060 }
1061 }
1062 }
1063 }
1064
1065 return triangle;
1066}
1067;
1068
1069/** Finds the starting triangle for FindNonConvexBorder().
1070 * Looks at the outermost point per axis, then FindSecondPointForTesselation()
1071 * for the second and FindNextSuitablePointViaAngleOfSphere() for the third
1072 * point are called.
1073 * \param *out output stream for debugging
1074 * \param RADIUS radius of virtual rolling sphere
1075 * \param *LC LinkedCell structure with neighbouring TesselPoint's
1076 * \return true - a starting triangle has been created, false - no valid triple of points found
1077 */
1078bool Tesselation::FindStartingTriangle(const double RADIUS, const LinkedCell *LC)
1079{
1080 Info FunctionInfo(__func__);
1081 int i = 0;
1082 TesselPoint* MaxPoint[NDIM];
1083 TesselPoint* Temporary;
1084 double maxCoordinate[NDIM];
1085 BoundaryLineSet *BaseLine = NULL;
1086 Vector helper;
1087 Vector Chord;
1088 Vector SearchDirection;
1089 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers
1090 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
1091 Vector SphereCenter;
1092 Vector NormalVector;
1093
1094 NormalVector.Zero();
1095
1096 for (i = 0; i < 3; i++) {
1097 MaxPoint[i] = NULL;
1098 maxCoordinate[i] = -1;
1099 }
1100
1101 // 1. searching topmost point with respect to each axis
1102 for (int i = 0; i < NDIM; i++) { // each axis
1103 LC->n[i] = LC->N[i] - 1; // current axis is topmost cell
1104 const int map[NDIM] = {i, (i + 1) % NDIM, (i + 2) % NDIM};
1105 for (LC->n[map[1]] = 0; LC->n[map[1]] < LC->N[map[1]]; LC->n[map[1]]++)
1106 for (LC->n[map[2]] = 0; LC->n[map[2]] < LC->N[map[2]]; LC->n[map[2]]++) {
1107 const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
1108 //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
1109 if (List != NULL) {
1110 for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
1111 if ((*Runner)->at(map[0]) > maxCoordinate[map[0]]) {
1112 DoLog(1) && (Log() << Verbose(1) << "New maximal for axis " << map[0] << " node is " << *(*Runner) << " at " << (*Runner)->getPosition() << "." << endl);
1113 maxCoordinate[map[0]] = (*Runner)->at(map[0]);
1114 MaxPoint[map[0]] = (*Runner);
1115 }
1116 }
1117 } else {
1118 DoeLog(1) && (eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
1119 }
1120 }
1121 }
1122
1123 DoLog(1) && (Log() << Verbose(1) << "Found maximum coordinates: ");
1124 for (int i = 0; i < NDIM; i++)
1125 DoLog(0) && (Log() << Verbose(0) << i << ": " << *MaxPoint[i] << "\t");
1126 DoLog(0) && (Log() << Verbose(0) << endl);
1127
1128 BTS = NULL;
1129 for (int k = 0; k < NDIM; k++) {
1130 NormalVector.Zero();
1131 NormalVector[k] = 1.;
1132 BaseLine = new BoundaryLineSet();
1133 BaseLine->endpoints[0] = new BoundaryPointSet(MaxPoint[k]);
1134 DoLog(0) && (Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl);
1135
1136 double ShortestAngle;
1137 ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant.
1138
1139 Temporary = NULL;
1140 FindSecondPointForTesselation(BaseLine->endpoints[0]->node, NormalVector, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_...
1141 if (Temporary == NULL) {
1142 // have we found a second point?
1143 delete BaseLine;
1144 continue;
1145 }
1146 BaseLine->endpoints[1] = new BoundaryPointSet(Temporary);
1147
1148 // construct center of circle
1149 CircleCenter = 0.5 * ((BaseLine->endpoints[0]->node->getPosition()) + (BaseLine->endpoints[1]->node->getPosition()));
1150
1151 // construct normal vector of circle
1152 CirclePlaneNormal = (BaseLine->endpoints[0]->node->getPosition()) - (BaseLine->endpoints[1]->node->getPosition());
1153
1154 double radius = CirclePlaneNormal.NormSquared();
1155 double CircleRadius = sqrt(RADIUS * RADIUS - radius / 4.);
1156
1157 NormalVector.ProjectOntoPlane(CirclePlaneNormal);
1158 NormalVector.Normalize();
1159 ShortestAngle = 2. * M_PI; // This will indicate the quadrant.
1160
1161 SphereCenter = (CircleRadius * NormalVector) + CircleCenter;
1162 // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized)
1163
1164 // look in one direction of baseline for initial candidate
1165 SearchDirection = Plane(CirclePlaneNormal, NormalVector,0).getNormal(); // whether we look "left" first or "right" first is not important ...
1166
1167 // adding point 1 and point 2 and add the line between them
1168 DoLog(0) && (Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl);
1169 DoLog(0) && (Log() << Verbose(0) << "Found second point is at " << *BaseLine->endpoints[1]->node << ".\n");
1170
1171 //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n";
1172 CandidateForTesselation OptCandidates(BaseLine);
1173 FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC);
1174 DoLog(0) && (Log() << Verbose(0) << "List of third Points is:" << endl);
1175 for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) {
1176 DoLog(0) && (Log() << Verbose(0) << " " << *(*it) << endl);
1177 }
1178 if (!OptCandidates.pointlist.empty()) {
1179 BTS = NULL;
1180 AddCandidatePolygon(OptCandidates, RADIUS, LC);
1181 } else {
1182 delete BaseLine;
1183 continue;
1184 }
1185
1186 if (BTS != NULL) { // we have created one starting triangle
1187 delete BaseLine;
1188 break;
1189 } else {
1190 // remove all candidates from the list and then the list itself
1191 OptCandidates.pointlist.clear();
1192 }
1193 delete BaseLine;
1194 }
1195
1196 return (BTS != NULL);
1197}
1198;
1199
1200/** Checks for a given baseline and a third point candidate whether baselines of the found triangle don't have even better candidates.
1201 * This is supposed to prevent early closing of the tesselation.
1202 * \param CandidateLine CandidateForTesselation with baseline and shortestangle , i.e. not \a *OptCandidate
1203 * \param *ThirdNode third point in triangle, not in BoundaryLineSet::endpoints
1204 * \param RADIUS radius of sphere
1205 * \param *LC LinkedCell structure
1206 * \return true - there is a better candidate (smaller angle than \a ShortestAngle), false - no better TesselPoint candidate found
1207 */
1208//bool Tesselation::HasOtherBaselineBetterCandidate(CandidateForTesselation &CandidateLine, const TesselPoint * const ThirdNode, double RADIUS, const LinkedCell * const LC) const
1209//{
1210// Info FunctionInfo(__func__);
1211// bool result = false;
1212// Vector CircleCenter;
1213// Vector CirclePlaneNormal;
1214// Vector OldSphereCenter;
1215// Vector SearchDirection;
1216// Vector helper;
1217// TesselPoint *OtherOptCandidate = NULL;
1218// double OtherShortestAngle = 2.*M_PI; // This will indicate the quadrant.
1219// double radius, CircleRadius;
1220// BoundaryLineSet *Line = NULL;
1221// BoundaryTriangleSet *T = NULL;
1222//
1223// // check both other lines
1224// PointMap::const_iterator FindPoint = PointsOnBoundary.find(ThirdNode->nr);
1225// if (FindPoint != PointsOnBoundary.end()) {
1226// for (int i=0;i<2;i++) {
1227// LineMap::const_iterator FindLine = (FindPoint->second)->lines.find(BaseRay->endpoints[0]->node->nr);
1228// if (FindLine != (FindPoint->second)->lines.end()) {
1229// Line = FindLine->second;
1230// Log() << Verbose(0) << "Found line " << *Line << "." << endl;
1231// if (Line->triangles.size() == 1) {
1232// T = Line->triangles.begin()->second;
1233// // construct center of circle
1234// CircleCenter.CopyVector(Line->endpoints[0]->node->node);
1235// CircleCenter.AddVector(Line->endpoints[1]->node->node);
1236// CircleCenter.Scale(0.5);
1237//
1238// // construct normal vector of circle
1239// CirclePlaneNormal.CopyVector(Line->endpoints[0]->node->node);
1240// CirclePlaneNormal.SubtractVector(Line->endpoints[1]->node->node);
1241//
1242// // calculate squared radius of circle
1243// radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal);
1244// if (radius/4. < RADIUS*RADIUS) {
1245// CircleRadius = RADIUS*RADIUS - radius/4.;
1246// CirclePlaneNormal.Normalize();
1247// //Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;
1248//
1249// // construct old center
1250// GetCenterofCircumcircle(&OldSphereCenter, *T->endpoints[0]->node->node, *T->endpoints[1]->node->node, *T->endpoints[2]->node->node);
1251// helper.CopyVector(&T->NormalVector); // normal vector ensures that this is correct center of the two possible ones
1252// radius = Line->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter);
1253// helper.Scale(sqrt(RADIUS*RADIUS - radius));
1254// OldSphereCenter.AddVector(&helper);
1255// OldSphereCenter.SubtractVector(&CircleCenter);
1256// //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;
1257//
1258// // construct SearchDirection
1259// SearchDirection.MakeNormalVector(&T->NormalVector, &CirclePlaneNormal);
1260// helper.CopyVector(Line->endpoints[0]->node->node);
1261// helper.SubtractVector(ThirdNode->node);
1262// if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
1263// SearchDirection.Scale(-1.);
1264// SearchDirection.ProjectOntoPlane(&OldSphereCenter);
1265// SearchDirection.Normalize();
1266// Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;
1267// if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) {
1268// // rotated the wrong way!
1269// DoeLog(1) && (eLog()<< Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);
1270// }
1271//
1272// // add third point
1273// FindThirdPointForTesselation(T->NormalVector, SearchDirection, OldSphereCenter, OptCandidates, ThirdNode, RADIUS, LC);
1274// for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); ++it) {
1275// if (((*it) == BaseRay->endpoints[0]->node) || ((*it) == BaseRay->endpoints[1]->node)) // skip if it's the same triangle than suggested
1276// continue;
1277// Log() << Verbose(0) << " Third point candidate is " << (*it)
1278// << " with circumsphere's center at " << (*it)->OptCenter << "." << endl;
1279// Log() << Verbose(0) << " Baseline is " << *BaseRay << endl;
1280//
1281// // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2)
1282// TesselPoint *PointCandidates[3];
1283// PointCandidates[0] = (*it);
1284// PointCandidates[1] = BaseRay->endpoints[0]->node;
1285// PointCandidates[2] = BaseRay->endpoints[1]->node;
1286// bool check=false;
1287// int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates);
1288// // If there is no triangle, add it regularly.
1289// if (existentTrianglesCount == 0) {
1290// SetTesselationPoint((*it), 0);
1291// SetTesselationPoint(BaseRay->endpoints[0]->node, 1);
1292// SetTesselationPoint(BaseRay->endpoints[1]->node, 2);
1293//
1294// if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) {
1295// OtherOptCandidate = (*it);
1296// check = true;
1297// }
1298// } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time.
1299// SetTesselationPoint((*it), 0);
1300// SetTesselationPoint(BaseRay->endpoints[0]->node, 1);
1301// SetTesselationPoint(BaseRay->endpoints[1]->node, 2);
1302//
1303// // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1)
1304// // i.e. at least one of the three lines must be present with TriangleCount <= 1
1305// if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS)) {
1306// OtherOptCandidate = (*it);
1307// check = true;
1308// }
1309// }
1310//
1311// if (check) {
1312// if (ShortestAngle > OtherShortestAngle) {
1313// Log() << Verbose(0) << "There is a better candidate than " << *ThirdNode << " with " << ShortestAngle << " from baseline " << *Line << ": " << *OtherOptCandidate << " with " << OtherShortestAngle << "." << endl;
1314// result = true;
1315// break;
1316// }
1317// }
1318// }
1319// delete(OptCandidates);
1320// if (result)
1321// break;
1322// } else {
1323// Log() << Verbose(0) << "Circumcircle for base line " << *Line << " and base triangle " << T << " is too big!" << endl;
1324// }
1325// } else {
1326// DoeLog(2) && (eLog()<< Verbose(2) << "Baseline is connected to two triangles already?" << endl);
1327// }
1328// } else {
1329// Log() << Verbose(1) << "No present baseline between " << BaseRay->endpoints[0] << " and candidate " << *ThirdNode << "." << endl;
1330// }
1331// }
1332// } else {
1333// DoeLog(1) && (eLog()<< Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl);
1334// }
1335//
1336// return result;
1337//};
1338
1339/** This function finds a triangle to a line, adjacent to an existing one.
1340 * @param out output stream for debugging
1341 * @param CandidateLine current cadndiate baseline to search from
1342 * @param T current triangle which \a Line is edge of
1343 * @param RADIUS radius of the rolling ball
1344 * @param N number of found triangles
1345 * @param *LC LinkedCell structure with neighbouring points
1346 */
1347bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, const BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)
1348{
1349 Info FunctionInfo(__func__);
1350 Vector CircleCenter;
1351 Vector CirclePlaneNormal;
1352 Vector RelativeSphereCenter;
1353 Vector SearchDirection;
1354 Vector helper;
1355 BoundaryPointSet *ThirdPoint = NULL;
1356 LineMap::iterator testline;
1357 double radius, CircleRadius;
1358
1359 for (int i = 0; i < 3; i++)
1360 if ((T.endpoints[i] != CandidateLine.BaseLine->endpoints[0]) && (T.endpoints[i] != CandidateLine.BaseLine->endpoints[1])) {
1361 ThirdPoint = T.endpoints[i];
1362 break;
1363 }
1364 DoLog(0) && (Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " with ThirdPoint " << *ThirdPoint << " of triangle " << T << "." << endl);
1365
1366 CandidateLine.T = &T;
1367
1368 // construct center of circle
1369 CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
1370 (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
1371
1372 // construct normal vector of circle
1373 CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
1374 (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
1375
1376 // calculate squared radius of circle
1377 radius = CirclePlaneNormal.ScalarProduct(CirclePlaneNormal);
1378 if (radius / 4. < RADIUS * RADIUS) {
1379 // construct relative sphere center with now known CircleCenter
1380 RelativeSphereCenter = T.SphereCenter - CircleCenter;
1381
1382 CircleRadius = RADIUS * RADIUS - radius / 4.;
1383 CirclePlaneNormal.Normalize();
1384 DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);
1385
1386 DoLog(1) && (Log() << Verbose(1) << "INFO: OldSphereCenter is at " << T.SphereCenter << "." << endl);
1387
1388 // construct SearchDirection and an "outward pointer"
1389 SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal,0).getNormal();
1390 helper = CircleCenter - (ThirdPoint->node->getPosition());
1391 if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards!
1392 SearchDirection.Scale(-1.);
1393 DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
1394 if (fabs(RelativeSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) {
1395 // rotated the wrong way!
1396 DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);
1397 }
1398
1399 // add third point
1400 FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, ThirdPoint, RADIUS, LC);
1401
1402 } else {
1403 DoLog(0) && (Log() << Verbose(0) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!" << endl);
1404 }
1405
1406 if (CandidateLine.pointlist.empty()) {
1407 DoeLog(2) && (eLog() << Verbose(2) << "Could not find a suitable candidate." << endl);
1408 return false;
1409 }
1410 DoLog(0) && (Log() << Verbose(0) << "Third Points are: " << endl);
1411 for (TesselPointList::iterator it = CandidateLine.pointlist.begin(); it != CandidateLine.pointlist.end(); ++it) {
1412 DoLog(0) && (Log() << Verbose(0) << " " << *(*it) << endl);
1413 }
1414
1415 return true;
1416}
1417;
1418
1419/** Walks through Tesselation::OpenLines() and finds candidates for newly created ones.
1420 * \param *&LCList atoms in LinkedCell list
1421 * \param RADIUS radius of the virtual sphere
1422 * \return true - for all open lines without candidates so far, a candidate has been found,
1423 * false - at least one open line without candidate still
1424 */
1425bool Tesselation::FindCandidatesforOpenLines(const double RADIUS, const LinkedCell *&LCList)
1426{
1427 bool TesselationFailFlag = true;
1428 CandidateForTesselation *baseline = NULL;
1429 BoundaryTriangleSet *T = NULL;
1430
1431 for (CandidateMap::iterator Runner = OpenLines.begin(); Runner != OpenLines.end(); Runner++) {
1432 baseline = Runner->second;
1433 if (baseline->pointlist.empty()) {
1434 ASSERT((baseline->BaseLine->triangles.size() == 1),"Open line without exactly one attached triangle");
1435 T = (((baseline->BaseLine->triangles.begin()))->second);
1436 DoLog(1) && (Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl);
1437 TesselationFailFlag = TesselationFailFlag && FindNextSuitableTriangle(*baseline, *T, RADIUS, LCList); //the line is there, so there is a triangle, but only one.
1438 }
1439 }
1440 return TesselationFailFlag;
1441}
1442;
1443
1444/** Adds the present line and candidate point from \a &CandidateLine to the Tesselation.
1445 * \param CandidateLine triangle to add
1446 * \param RADIUS Radius of sphere
1447 * \param *LC LinkedCell structure
1448 * \NOTE we need the copy operator here as the original CandidateForTesselation is removed in
1449 * AddTesselationLine() in AddCandidateTriangle()
1450 */
1451void Tesselation::AddCandidatePolygon(CandidateForTesselation CandidateLine, const double RADIUS, const LinkedCell *LC)
1452{
1453 Info FunctionInfo(__func__);
1454 Vector Center;
1455 TesselPoint * const TurningPoint = CandidateLine.BaseLine->endpoints[0]->node;
1456 TesselPointList::iterator Runner;
1457 TesselPointList::iterator Sprinter;
1458
1459 // fill the set of neighbours
1460 TesselPointSet SetOfNeighbours;
1461
1462 SetOfNeighbours.insert(CandidateLine.BaseLine->endpoints[1]->node);
1463 for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)
1464 SetOfNeighbours.insert(*Runner);
1465 TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->getPosition());
1466
1467 DoLog(0) && (Log() << Verbose(0) << "List of Candidates for Turning Point " << *TurningPoint << ":" << endl);
1468 for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner)
1469 DoLog(0) && (Log() << Verbose(0) << " " << **TesselRunner << endl);
1470
1471 // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles)
1472 Runner = connectedClosestPoints->begin();
1473 Sprinter = Runner;
1474 Sprinter++;
1475 while (Sprinter != connectedClosestPoints->end()) {
1476 DoLog(0) && (Log() << Verbose(0) << "Current Runner is " << *(*Runner) << " and sprinter is " << *(*Sprinter) << "." << endl);
1477
1478 AddTesselationPoint(TurningPoint, 0);
1479 AddTesselationPoint(*Runner, 1);
1480 AddTesselationPoint(*Sprinter, 2);
1481
1482 AddCandidateTriangle(CandidateLine, Opt);
1483
1484 Runner = Sprinter;
1485 Sprinter++;
1486 if (Sprinter != connectedClosestPoints->end()) {
1487 // fill the internal open lines with its respective candidate (otherwise lines in degenerate case are not picked)
1488 FindDegeneratedCandidatesforOpenLines(*Sprinter, &CandidateLine.OptCenter); // Assume BTS contains last triangle
1489 DoLog(0) && (Log() << Verbose(0) << " There are still more triangles to add." << endl);
1490 }
1491 // pick candidates for other open lines as well
1492 FindCandidatesforOpenLines(RADIUS, LC);
1493
1494 // check whether we add a degenerate or a normal triangle
1495 if (CheckDegeneracy(CandidateLine, RADIUS, LC)) {
1496 // add normal and degenerate triangles
1497 DoLog(1) && (Log() << Verbose(1) << "Triangle of endpoints " << *TPS[0] << "," << *TPS[1] << " and " << *TPS[2] << " is degenerated, adding both sides." << endl);
1498 AddCandidateTriangle(CandidateLine, OtherOpt);
1499
1500 if (Sprinter != connectedClosestPoints->end()) {
1501 // fill the internal open lines with its respective candidate (otherwise lines in degenerate case are not picked)
1502 FindDegeneratedCandidatesforOpenLines(*Sprinter, &CandidateLine.OtherOptCenter);
1503 }
1504 // pick candidates for other open lines as well
1505 FindCandidatesforOpenLines(RADIUS, LC);
1506 }
1507 }
1508 delete (connectedClosestPoints);
1509};
1510
1511/** for polygons (multiple candidates for a baseline) sets internal edges to the correct next candidate.
1512 * \param *Sprinter next candidate to which internal open lines are set
1513 * \param *OptCenter OptCenter for this candidate
1514 */
1515void Tesselation::FindDegeneratedCandidatesforOpenLines(TesselPoint * const Sprinter, const Vector * const OptCenter)
1516{
1517 Info FunctionInfo(__func__);
1518
1519 pair<LineMap::iterator, LineMap::iterator> FindPair = TPS[0]->lines.equal_range(TPS[2]->node->nr);
1520 for (LineMap::const_iterator FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) {
1521 DoLog(1) && (Log() << Verbose(1) << "INFO: Checking line " << *(FindLine->second) << " ..." << endl);
1522 // If there is a line with less than two attached triangles, we don't need a new line.
1523 if (FindLine->second->triangles.size() == 1) {
1524 CandidateMap::iterator Finder = OpenLines.find(FindLine->second);
1525 if (!Finder->second->pointlist.empty())
1526 DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << "." << endl);
1527 else {
1528 DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate, setting to next Sprinter" << (*Sprinter) << endl);
1529 Finder->second->T = BTS; // is last triangle
1530 Finder->second->pointlist.push_back(Sprinter);
1531 Finder->second->ShortestAngle = 0.;
1532 Finder->second->OptCenter = *OptCenter;
1533 }
1534 }
1535 }
1536};
1537
1538/** If a given \a *triangle is degenerated, this adds both sides.
1539 * i.e. the triangle with same BoundaryPointSet's but NormalVector in opposite direction.
1540 * Note that endpoints are stored in Tesselation::TPS
1541 * \param CandidateLine CanddiateForTesselation structure for the desired BoundaryLine
1542 * \param RADIUS radius of sphere
1543 * \param *LC pointer to LinkedCell structure
1544 */
1545void Tesselation::AddDegeneratedTriangle(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC)
1546{
1547 Info FunctionInfo(__func__);
1548 Vector Center;
1549 CandidateMap::const_iterator CandidateCheck = OpenLines.end();
1550 BoundaryTriangleSet *triangle = NULL;
1551
1552 /// 1. Create or pick the lines for the first triangle
1553 DoLog(0) && (Log() << Verbose(0) << "INFO: Creating/Picking lines for first triangle ..." << endl);
1554 for (int i = 0; i < 3; i++) {
1555 BLS[i] = NULL;
1556 DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
1557 AddTesselationLine(&CandidateLine.OptCenter, TPS[(i + 2) % 3], TPS[(i + 0) % 3], TPS[(i + 1) % 3], i);
1558 }
1559
1560 /// 2. create the first triangle and NormalVector and so on
1561 DoLog(0) && (Log() << Verbose(0) << "INFO: Adding first triangle with center at " << CandidateLine.OptCenter << " ..." << endl);
1562 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
1563 AddTesselationTriangle();
1564
1565 // create normal vector
1566 BTS->GetCenter(Center);
1567 Center -= CandidateLine.OptCenter;
1568 BTS->SphereCenter = CandidateLine.OptCenter;
1569 BTS->GetNormalVector(Center);
1570 // give some verbose output about the whole procedure
1571 if (CandidateLine.T != NULL)
1572 DoLog(0) && (Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl);
1573 else
1574 DoLog(0) && (Log() << Verbose(0) << "--> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl);
1575 triangle = BTS;
1576
1577 /// 3. Gather candidates for each new line
1578 DoLog(0) && (Log() << Verbose(0) << "INFO: Adding candidates to new lines ..." << endl);
1579 for (int i = 0; i < 3; i++) {
1580 DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
1581 CandidateCheck = OpenLines.find(BLS[i]);
1582 if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) {
1583 if (CandidateCheck->second->T == NULL)
1584 CandidateCheck->second->T = triangle;
1585 FindNextSuitableTriangle(*(CandidateCheck->second), *CandidateCheck->second->T, RADIUS, LC);
1586 }
1587 }
1588
1589 /// 4. Create or pick the lines for the second triangle
1590 DoLog(0) && (Log() << Verbose(0) << "INFO: Creating/Picking lines for second triangle ..." << endl);
1591 for (int i = 0; i < 3; i++) {
1592 DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
1593 AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i + 2) % 3], TPS[(i + 0) % 3], TPS[(i + 1) % 3], i);
1594 }
1595
1596 /// 5. create the second triangle and NormalVector and so on
1597 DoLog(0) && (Log() << Verbose(0) << "INFO: Adding second triangle with center at " << CandidateLine.OtherOptCenter << " ..." << endl);
1598 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
1599 AddTesselationTriangle();
1600
1601 BTS->SphereCenter = CandidateLine.OtherOptCenter;
1602 // create normal vector in other direction
1603 BTS->GetNormalVector(triangle->NormalVector);
1604 BTS->NormalVector.Scale(-1.);
1605 // give some verbose output about the whole procedure
1606 if (CandidateLine.T != NULL)
1607 DoLog(0) && (Log() << Verbose(0) << "--> New degenerate triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl);
1608 else
1609 DoLog(0) && (Log() << Verbose(0) << "--> New degenerate starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl);
1610
1611 /// 6. Adding triangle to new lines
1612 DoLog(0) && (Log() << Verbose(0) << "INFO: Adding second triangles to new lines ..." << endl);
1613 for (int i = 0; i < 3; i++) {
1614 DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl);
1615 CandidateCheck = OpenLines.find(BLS[i]);
1616 if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) {
1617 if (CandidateCheck->second->T == NULL)
1618 CandidateCheck->second->T = BTS;
1619 }
1620 }
1621}
1622;
1623
1624/** Adds a triangle to the Tesselation structure from three given TesselPoint's.
1625 * Note that endpoints are in Tesselation::TPS.
1626 * \param CandidateLine CandidateForTesselation structure contains other information
1627 * \param type which opt center to add (i.e. which side) and thus which NormalVector to take
1628 */
1629void Tesselation::AddCandidateTriangle(CandidateForTesselation &CandidateLine, enum centers type)
1630{
1631 Info FunctionInfo(__func__);
1632 Vector Center;
1633 Vector *OptCenter = (type == Opt) ? &CandidateLine.OptCenter : &CandidateLine.OtherOptCenter;
1634
1635 // add the lines
1636 AddTesselationLine(OptCenter, TPS[2], TPS[0], TPS[1], 0);
1637 AddTesselationLine(OptCenter, TPS[1], TPS[0], TPS[2], 1);
1638 AddTesselationLine(OptCenter, TPS[0], TPS[1], TPS[2], 2);
1639
1640 // add the triangles
1641 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
1642 AddTesselationTriangle();
1643
1644 // create normal vector
1645 BTS->GetCenter(Center);
1646 Center.SubtractVector(*OptCenter);
1647 BTS->SphereCenter = *OptCenter;
1648 BTS->GetNormalVector(Center);
1649
1650 // give some verbose output about the whole procedure
1651 if (CandidateLine.T != NULL)
1652 DoLog(0) && (Log() << Verbose(0) << "--> New" << ((type == OtherOpt) ? " degenerate " : " ") << "triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl);
1653 else
1654 DoLog(0) && (Log() << Verbose(0) << "--> New" << ((type == OtherOpt) ? " degenerate " : " ") << "starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl);
1655}
1656;
1657
1658/** Checks whether the quadragon of the two triangles connect to \a *Base is convex.
1659 * We look whether the closest point on \a *Base with respect to the other baseline is outside
1660 * of the segment formed by both endpoints (concave) or not (convex).
1661 * \param *out output stream for debugging
1662 * \param *Base line to be flipped
1663 * \return NULL - convex, otherwise endpoint that makes it concave
1664 */
1665class BoundaryPointSet *Tesselation::IsConvexRectangle(class BoundaryLineSet *Base)
1666{
1667 Info FunctionInfo(__func__);
1668 class BoundaryPointSet *Spot = NULL;
1669 class BoundaryLineSet *OtherBase;
1670 Vector *ClosestPoint;
1671
1672 int m = 0;
1673 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
1674 for (int j = 0; j < 3; j++) // all of their endpoints and baselines
1675 if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
1676 BPS[m++] = runner->second->endpoints[j];
1677 OtherBase = new class BoundaryLineSet(BPS, -1);
1678
1679 DoLog(1) && (Log() << Verbose(1) << "INFO: Current base line is " << *Base << "." << endl);
1680 DoLog(1) && (Log() << Verbose(1) << "INFO: Other base line is " << *OtherBase << "." << endl);
1681
1682 // get the closest point on each line to the other line
1683 ClosestPoint = GetClosestPointBetweenLine(Base, OtherBase);
1684
1685 // delete the temporary other base line
1686 delete (OtherBase);
1687
1688 // get the distance vector from Base line to OtherBase line
1689 Vector DistanceToIntersection[2], BaseLine;
1690 double distance[2];
1691 BaseLine = (Base->endpoints[1]->node->getPosition()) - (Base->endpoints[0]->node->getPosition());
1692 for (int i = 0; i < 2; i++) {
1693 DistanceToIntersection[i] = (*ClosestPoint) - (Base->endpoints[i]->node->getPosition());
1694 distance[i] = BaseLine.ScalarProduct(DistanceToIntersection[i]);
1695 }
1696 delete (ClosestPoint);
1697 if ((distance[0] * distance[1]) > 0) { // have same sign?
1698 DoLog(1) && (Log() << Verbose(1) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1] << ". " << *Base << "' rectangle is concave." << endl);
1699 if (distance[0] < distance[1]) {
1700 Spot = Base->endpoints[0];
1701 } else {
1702 Spot = Base->endpoints[1];
1703 }
1704 return Spot;
1705 } else { // different sign, i.e. we are in between
1706 DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl);
1707 return NULL;
1708 }
1709
1710}
1711;
1712
1713void Tesselation::PrintAllBoundaryPoints(ofstream *out) const
1714{
1715 Info FunctionInfo(__func__);
1716 // print all lines
1717 DoLog(0) && (Log() << Verbose(0) << "Printing all boundary points for debugging:" << endl);
1718 for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin(); PointRunner != PointsOnBoundary.end(); PointRunner++)
1719 DoLog(0) && (Log() << Verbose(0) << *(PointRunner->second) << endl);
1720}
1721;
1722
1723void Tesselation::PrintAllBoundaryLines(ofstream *out) const
1724{
1725 Info FunctionInfo(__func__);
1726 // print all lines
1727 DoLog(0) && (Log() << Verbose(0) << "Printing all boundary lines for debugging:" << endl);
1728 for (LineMap::const_iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++)
1729 DoLog(0) && (Log() << Verbose(0) << *(LineRunner->second) << endl);
1730}
1731;
1732
1733void Tesselation::PrintAllBoundaryTriangles(ofstream *out) const
1734{
1735 Info FunctionInfo(__func__);
1736 // print all triangles
1737 DoLog(0) && (Log() << Verbose(0) << "Printing all boundary triangles for debugging:" << endl);
1738 for (TriangleMap::const_iterator TriangleRunner = TrianglesOnBoundary.begin(); TriangleRunner != TrianglesOnBoundary.end(); TriangleRunner++)
1739 DoLog(0) && (Log() << Verbose(0) << *(TriangleRunner->second) << endl);
1740}
1741;
1742
1743/** For a given boundary line \a *Base and its two triangles, picks the central baseline that is "higher".
1744 * \param *out output stream for debugging
1745 * \param *Base line to be flipped
1746 * \return volume change due to flipping (0 - then no flipped occured)
1747 */
1748double Tesselation::PickFarthestofTwoBaselines(class BoundaryLineSet *Base)
1749{
1750 Info FunctionInfo(__func__);
1751 class BoundaryLineSet *OtherBase;
1752 Vector *ClosestPoint[2];
1753 double volume;
1754
1755 int m = 0;
1756 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
1757 for (int j = 0; j < 3; j++) // all of their endpoints and baselines
1758 if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
1759 BPS[m++] = runner->second->endpoints[j];
1760 OtherBase = new class BoundaryLineSet(BPS, -1);
1761
1762 DoLog(0) && (Log() << Verbose(0) << "INFO: Current base line is " << *Base << "." << endl);
1763 DoLog(0) && (Log() << Verbose(0) << "INFO: Other base line is " << *OtherBase << "." << endl);
1764
1765 // get the closest point on each line to the other line
1766 ClosestPoint[0] = GetClosestPointBetweenLine(Base, OtherBase);
1767 ClosestPoint[1] = GetClosestPointBetweenLine(OtherBase, Base);
1768
1769 // get the distance vector from Base line to OtherBase line
1770 Vector Distance = (*ClosestPoint[1]) - (*ClosestPoint[0]);
1771
1772 // calculate volume
1773 volume = CalculateVolumeofGeneralTetraeder(Base->endpoints[1]->node->getPosition(), OtherBase->endpoints[0]->node->getPosition(), OtherBase->endpoints[1]->node->getPosition(), Base->endpoints[0]->node->getPosition());
1774
1775 // delete the temporary other base line and the closest points
1776 delete (ClosestPoint[0]);
1777 delete (ClosestPoint[1]);
1778 delete (OtherBase);
1779
1780 if (Distance.NormSquared() < MYEPSILON) { // check for intersection
1781 DoLog(0) && (Log() << Verbose(0) << "REJECT: Both lines have an intersection: Nothing to do." << endl);
1782 return false;
1783 } else { // check for sign against BaseLineNormal
1784 Vector BaseLineNormal;
1785 BaseLineNormal.Zero();
1786 if (Base->triangles.size() < 2) {
1787 DoeLog(1) && (eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
1788 return 0.;
1789 }
1790 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
1791 DoLog(1) && (Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl);
1792 BaseLineNormal += (runner->second->NormalVector);
1793 }
1794 BaseLineNormal.Scale(1. / 2.);
1795
1796 if (Distance.ScalarProduct(BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip
1797 DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl);
1798 // calculate volume summand as a general tetraeder
1799 return volume;
1800 } else { // Base higher than OtherBase -> do nothing
1801 DoLog(0) && (Log() << Verbose(0) << "REJECT: Base line is higher: Nothing to do." << endl);
1802 return 0.;
1803 }
1804 }
1805}
1806;
1807
1808/** For a given baseline and its two connected triangles, flips the baseline.
1809 * I.e. we create the new baseline between the other two endpoints of these four
1810 * endpoints and reconstruct the two triangles accordingly.
1811 * \param *out output stream for debugging
1812 * \param *Base line to be flipped
1813 * \return pointer to allocated new baseline - flipping successful, NULL - something went awry
1814 */
1815class BoundaryLineSet * Tesselation::FlipBaseline(class BoundaryLineSet *Base)
1816{
1817 Info FunctionInfo(__func__);
1818 class BoundaryLineSet *OldLines[4], *NewLine;
1819 class BoundaryPointSet *OldPoints[2];
1820 Vector BaseLineNormal;
1821 int OldTriangleNrs[2], OldBaseLineNr;
1822 int i, m;
1823
1824 // calculate NormalVector for later use
1825 BaseLineNormal.Zero();
1826 if (Base->triangles.size() < 2) {
1827 DoeLog(1) && (eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);
1828 return NULL;
1829 }
1830 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {
1831 DoLog(1) && (Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl);
1832 BaseLineNormal += (runner->second->NormalVector);
1833 }
1834 BaseLineNormal.Scale(-1. / 2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()
1835
1836 // get the two triangles
1837 // gather four endpoints and four lines
1838 for (int j = 0; j < 4; j++)
1839 OldLines[j] = NULL;
1840 for (int j = 0; j < 2; j++)
1841 OldPoints[j] = NULL;
1842 i = 0;
1843 m = 0;
1844 DoLog(0) && (Log() << Verbose(0) << "The four old lines are: ");
1845 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
1846 for (int j = 0; j < 3; j++) // all of their endpoints and baselines
1847 if (runner->second->lines[j] != Base) { // pick not the central baseline
1848 OldLines[i++] = runner->second->lines[j];
1849 DoLog(0) && (Log() << Verbose(0) << *runner->second->lines[j] << "\t");
1850 }
1851 DoLog(0) && (Log() << Verbose(0) << endl);
1852 DoLog(0) && (Log() << Verbose(0) << "The two old points are: ");
1853 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)
1854 for (int j = 0; j < 3; j++) // all of their endpoints and baselines
1855 if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) { // and neither of its endpoints
1856 OldPoints[m++] = runner->second->endpoints[j];
1857 DoLog(0) && (Log() << Verbose(0) << *runner->second->endpoints[j] << "\t");
1858 }
1859 DoLog(0) && (Log() << Verbose(0) << endl);
1860
1861 // check whether everything is in place to create new lines and triangles
1862 if (i < 4) {
1863 DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl);
1864 return NULL;
1865 }
1866 for (int j = 0; j < 4; j++)
1867 if (OldLines[j] == NULL) {
1868 DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl);
1869 return NULL;
1870 }
1871 for (int j = 0; j < 2; j++)
1872 if (OldPoints[j] == NULL) {
1873 DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough endpoints!" << endl);
1874 return NULL;
1875 }
1876
1877 // remove triangles and baseline removes itself
1878 DoLog(0) && (Log() << Verbose(0) << "INFO: Deleting baseline " << *Base << " from global list." << endl);
1879 OldBaseLineNr = Base->Nr;
1880 m = 0;
1881 // first obtain all triangle to delete ... (otherwise we pull the carpet (Base) from under the for-loop's feet)
1882 list <BoundaryTriangleSet *> TrianglesOfBase;
1883 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); ++runner)
1884 TrianglesOfBase.push_back(runner->second);
1885 // .. then delete each triangle (which deletes the line as well)
1886 for (list <BoundaryTriangleSet *>::iterator runner = TrianglesOfBase.begin(); !TrianglesOfBase.empty(); runner = TrianglesOfBase.begin()) {
1887 DoLog(0) && (Log() << Verbose(0) << "INFO: Deleting triangle " << *(*runner) << "." << endl);
1888 OldTriangleNrs[m++] = (*runner)->Nr;
1889 RemoveTesselationTriangle((*runner));
1890 TrianglesOfBase.erase(runner);
1891 }
1892
1893 // construct new baseline (with same number as old one)
1894 BPS[0] = OldPoints[0];
1895 BPS[1] = OldPoints[1];
1896 NewLine = new class BoundaryLineSet(BPS, OldBaseLineNr);
1897 LinesOnBoundary.insert(LinePair(OldBaseLineNr, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one
1898 DoLog(0) && (Log() << Verbose(0) << "INFO: Created new baseline " << *NewLine << "." << endl);
1899
1900 // construct new triangles with flipped baseline
1901 i = -1;
1902 if (OldLines[0]->IsConnectedTo(OldLines[2]))
1903 i = 2;
1904 if (OldLines[0]->IsConnectedTo(OldLines[3]))
1905 i = 3;
1906 if (i != -1) {
1907 BLS[0] = OldLines[0];
1908 BLS[1] = OldLines[i];
1909 BLS[2] = NewLine;
1910 BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[0]);
1911 BTS->GetNormalVector(BaseLineNormal);
1912 AddTesselationTriangle(OldTriangleNrs[0]);
1913 DoLog(0) && (Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl);
1914
1915 BLS[0] = (i == 2 ? OldLines[3] : OldLines[2]);
1916 BLS[1] = OldLines[1];
1917 BLS[2] = NewLine;
1918 BTS = new class BoundaryTriangleSet(BLS, OldTriangleNrs[1]);
1919 BTS->GetNormalVector(BaseLineNormal);
1920 AddTesselationTriangle(OldTriangleNrs[1]);
1921 DoLog(0) && (Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl);
1922 } else {
1923 DoeLog(0) && (eLog() << Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl);
1924 return NULL;
1925 }
1926
1927 return NewLine;
1928}
1929;
1930
1931/** Finds the second point of starting triangle.
1932 * \param *a first node
1933 * \param Oben vector indicating the outside
1934 * \param OptCandidate reference to recommended candidate on return
1935 * \param Storage[3] array storing angles and other candidate information
1936 * \param RADIUS radius of virtual sphere
1937 * \param *LC LinkedCell structure with neighbouring points
1938 */
1939void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC)
1940{
1941 Info FunctionInfo(__func__);
1942 Vector AngleCheck;
1943 class TesselPoint* Candidate = NULL;
1944 double norm = -1.;
1945 double angle = 0.;
1946 int N[NDIM];
1947 int Nlower[NDIM];
1948 int Nupper[NDIM];
1949
1950 if (LC->SetIndexToNode(a)) { // get cell for the starting point
1951 for (int i = 0; i < NDIM; i++) // store indices of this cell
1952 N[i] = LC->n[i];
1953 } else {
1954 DoeLog(1) && (eLog() << Verbose(1) << "Point " << *a << " is not found in cell " << LC->index << "." << endl);
1955 return;
1956 }
1957 // then go through the current and all neighbouring cells and check the contained points for possible candidates
1958 for (int i = 0; i < NDIM; i++) {
1959 Nlower[i] = ((N[i] - 1) >= 0) ? N[i] - 1 : 0;
1960 Nupper[i] = ((N[i] + 1) < LC->N[i]) ? N[i] + 1 : LC->N[i] - 1;
1961 }
1962 DoLog(0) && (Log() << Verbose(0) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :" << " [" << Nlower[0] << "," << Nupper[0] << "], " << " [" << Nlower[1] << "," << Nupper[1] << "], " << " [" << Nlower[2] << "," << Nupper[2] << "], " << endl);
1963
1964 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
1965 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
1966 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
1967 const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
1968 //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
1969 if (List != NULL) {
1970 for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
1971 Candidate = (*Runner);
1972 // check if we only have one unique point yet ...
1973 if (a != Candidate) {
1974 // Calculate center of the circle with radius RADIUS through points a and Candidate
1975 Vector OrthogonalizedOben, aCandidate, Center;
1976 double distance, scaleFactor;
1977
1978 OrthogonalizedOben = Oben;
1979 aCandidate = (a->getPosition()) - (Candidate->getPosition());
1980 OrthogonalizedOben.ProjectOntoPlane(aCandidate);
1981 OrthogonalizedOben.Normalize();
1982 distance = 0.5 * aCandidate.Norm();
1983 scaleFactor = sqrt(((RADIUS * RADIUS) - (distance * distance)));
1984 OrthogonalizedOben.Scale(scaleFactor);
1985
1986 Center = 0.5 * ((Candidate->getPosition()) + (a->getPosition()));
1987 Center += OrthogonalizedOben;
1988
1989 AngleCheck = Center - (a->getPosition());
1990 norm = aCandidate.Norm();
1991 // second point shall have smallest angle with respect to Oben vector
1992 if (norm < RADIUS * 2.) {
1993 angle = AngleCheck.Angle(Oben);
1994 if (angle < Storage[0]) {
1995 //Log() << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);
1996 DoLog(1) && (Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n");
1997 OptCandidate = Candidate;
1998 Storage[0] = angle;
1999 //Log() << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);
2000 } else {
2001 //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *OptCandidate << endl;
2002 }
2003 } else {
2004 //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;
2005 }
2006 } else {
2007 //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;
2008 }
2009 }
2010 } else {
2011 DoLog(0) && (Log() << Verbose(0) << "Linked cell list is empty." << endl);
2012 }
2013 }
2014}
2015;
2016
2017/** This recursive function finds a third point, to form a triangle with two given ones.
2018 * Note that this function is for the starting triangle.
2019 * The idea is as follows: A sphere with fixed radius is (almost) uniquely defined in space by three points
2020 * that sit on its boundary. Hence, when two points are given and we look for the (next) third point, then
2021 * the center of the sphere is still fixed up to a single parameter. The band of possible values
2022 * describes a circle in 3D-space. The old center of the sphere for the current base triangle gives
2023 * us the "null" on this circle, the new center of the candidate point will be some way along this
2024 * circle. The shorter the way the better is the candidate. Note that the direction is clearly given
2025 * by the normal vector of the base triangle that always points outwards by construction.
2026 * Hence, we construct a Center of this circle which sits right in the middle of the current base line.
2027 * We construct the normal vector that defines the plane this circle lies in, it is just in the
2028 * direction of the baseline. And finally, we need the radius of the circle, which is given by the rest
2029 * with respect to the length of the baseline and the sphere's fixed \a RADIUS.
2030 * Note that there is one difficulty: The circumcircle is uniquely defined, but for the circumsphere's center
2031 * there are two possibilities which becomes clear from the construction as seen below. Hence, we must check
2032 * both.
2033 * Note also that the acos() function is not unique on [0, 2.*M_PI). Hence, we need an additional check
2034 * to decide for one of the two possible angles. Therefore we need a SearchDirection and to make this check
2035 * sensible we need OldSphereCenter to be orthogonal to it. Either we construct SearchDirection orthogonal
2036 * right away, or -- what we do here -- we rotate the relative sphere centers such that this orthogonality
2037 * holds. Then, the normalized projection onto the SearchDirection is either +1 or -1 and thus states whether
2038 * the angle is uniquely in either (0,M_PI] or [M_PI, 2.*M_PI).
2039 * @param NormalVector normal direction of the base triangle (here the unit axis vector, \sa FindStartingTriangle())
2040 * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine
2041 * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle
2042 * @param CandidateLine CandidateForTesselation with the current base line and list of candidates and ShortestAngle
2043 * @param ThirdPoint third point to avoid in search
2044 * @param RADIUS radius of sphere
2045 * @param *LC LinkedCell structure with neighbouring points
2046 */
2047void Tesselation::FindThirdPointForTesselation(const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet * const ThirdPoint, const double RADIUS, const LinkedCell *LC) const
2048{
2049 Info FunctionInfo(__func__);
2050 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers
2051 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in
2052 Vector SphereCenter;
2053 Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility
2054 Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility
2055 Vector NewNormalVector; // normal vector of the Candidate's triangle
2056 Vector helper, OptCandidateCenter, OtherOptCandidateCenter;
2057 Vector RelativeOldSphereCenter;
2058 Vector NewPlaneCenter;
2059 double CircleRadius; // radius of this circle
2060 double radius;
2061 double otherradius;
2062 double alpha, Otheralpha; // angles (i.e. parameter for the circle).
2063 int N[NDIM], Nlower[NDIM], Nupper[NDIM];
2064 TesselPoint *Candidate = NULL;
2065
2066 DoLog(1) && (Log() << Verbose(1) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl);
2067
2068 // copy old center
2069 CandidateLine.OldCenter = OldSphereCenter;
2070 CandidateLine.ThirdPoint = ThirdPoint;
2071 CandidateLine.pointlist.clear();
2072
2073 // construct center of circle
2074 CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) +
2075 (CandidateLine.BaseLine->endpoints[1]->node->getPosition()));
2076
2077 // construct normal vector of circle
2078 CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) -
2079 (CandidateLine.BaseLine->endpoints[1]->node->getPosition());
2080
2081 RelativeOldSphereCenter = OldSphereCenter - CircleCenter;
2082
2083 // calculate squared radius TesselPoint *ThirdPoint,f circle
2084 radius = CirclePlaneNormal.NormSquared() / 4.;
2085 if (radius < RADIUS * RADIUS) {
2086 CircleRadius = RADIUS * RADIUS - radius;
2087 CirclePlaneNormal.Normalize();
2088 DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);
2089
2090 // test whether old center is on the band's plane
2091 if (fabs(RelativeOldSphereCenter.ScalarProduct(CirclePlaneNormal)) > HULLEPSILON) {
2092 DoeLog(1) && (eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(CirclePlaneNormal)) << "!" << endl);
2093 RelativeOldSphereCenter.ProjectOntoPlane(CirclePlaneNormal);
2094 }
2095 radius = RelativeOldSphereCenter.NormSquared();
2096 if (fabs(radius - CircleRadius) < HULLEPSILON) {
2097 DoLog(1) && (Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl);
2098
2099 // check SearchDirection
2100 DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
2101 if (fabs(RelativeOldSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) { // rotated the wrong way!
2102 DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl);
2103 }
2104
2105 // get cell for the starting point
2106 if (LC->SetIndexToVector(CircleCenter)) {
2107 for (int i = 0; i < NDIM; i++) // store indices of this cell
2108 N[i] = LC->n[i];
2109 //Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
2110 } else {
2111 DoeLog(1) && (eLog() << Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl);
2112 return;
2113 }
2114 // then go through the current and all neighbouring cells and check the contained points for possible candidates
2115 //Log() << Verbose(1) << "LC Intervals:";
2116 for (int i = 0; i < NDIM; i++) {
2117 Nlower[i] = ((N[i] - 1) >= 0) ? N[i] - 1 : 0;
2118 Nupper[i] = ((N[i] + 1) < LC->N[i]) ? N[i] + 1 : LC->N[i] - 1;
2119 //Log() << Verbose(0) << " [" << Nlower[i] << "," << Nupper[i] << "] ";
2120 }
2121 //Log() << Verbose(0) << endl;
2122 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
2123 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
2124 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
2125 const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
2126 //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;
2127 if (List != NULL) {
2128 for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
2129 Candidate = (*Runner);
2130
2131 // check for three unique points
2132 DoLog(2) && (Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " for BaseLine " << *CandidateLine.BaseLine << " with OldSphereCenter " << OldSphereCenter << "." << endl);
2133 if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node)) {
2134
2135 // find center on the plane
2136 GetCenterofCircumcircle(NewPlaneCenter, CandidateLine.BaseLine->endpoints[0]->node->getPosition(), CandidateLine.BaseLine->endpoints[1]->node->getPosition(), Candidate->getPosition());
2137 DoLog(1) && (Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl);
2138
2139 try {
2140 NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()),
2141 (CandidateLine.BaseLine->endpoints[1]->node->getPosition()),
2142 (Candidate->getPosition())).getNormal();
2143 DoLog(1) && (Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl);
2144 radius = CandidateLine.BaseLine->endpoints[0]->node->DistanceSquared(NewPlaneCenter);
2145 DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);
2146 DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);
2147 DoLog(1) && (Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl);
2148 if (radius < RADIUS * RADIUS) {
2149 otherradius = CandidateLine.BaseLine->endpoints[1]->node->DistanceSquared(NewPlaneCenter);
2150 if (fabs(radius - otherradius) < HULLEPSILON) {
2151 // construct both new centers
2152 NewSphereCenter = NewPlaneCenter;
2153 OtherNewSphereCenter= NewPlaneCenter;
2154 helper = NewNormalVector;
2155 helper.Scale(sqrt(RADIUS * RADIUS - radius));
2156 DoLog(2) && (Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl);
2157 NewSphereCenter += helper;
2158 DoLog(2) && (Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl);
2159 // OtherNewSphereCenter is created by the same vector just in the other direction
2160 helper.Scale(-1.);
2161 OtherNewSphereCenter += helper;
2162 DoLog(2) && (Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl);
2163 alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection, HULLEPSILON);
2164 Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection, HULLEPSILON);
2165 if ((ThirdPoint != NULL) && (Candidate == ThirdPoint->node)) { // in that case only the other circlecenter is valid
2166 if (OldSphereCenter.DistanceSquared(NewSphereCenter) < OldSphereCenter.DistanceSquared(OtherNewSphereCenter))
2167 alpha = Otheralpha;
2168 } else
2169 alpha = min(alpha, Otheralpha);
2170 // if there is a better candidate, drop the current list and add the new candidate
2171 // otherwise ignore the new candidate and keep the list
2172 if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) {
2173 if (fabs(alpha - Otheralpha) > MYEPSILON) {
2174 CandidateLine.OptCenter = NewSphereCenter;
2175 CandidateLine.OtherOptCenter = OtherNewSphereCenter;
2176 } else {
2177 CandidateLine.OptCenter = OtherNewSphereCenter;
2178 CandidateLine.OtherOptCenter = NewSphereCenter;
2179 }
2180 // if there is an equal candidate, add it to the list without clearing the list
2181 if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) {
2182 CandidateLine.pointlist.push_back(Candidate);
2183 DoLog(0) && (Log() << Verbose(0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl);
2184 } else {
2185 // remove all candidates from the list and then the list itself
2186 CandidateLine.pointlist.clear();
2187 CandidateLine.pointlist.push_back(Candidate);
2188 DoLog(0) && (Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl);
2189 }
2190 CandidateLine.ShortestAngle = alpha;
2191 DoLog(0) && (Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << endl);
2192 } else {
2193 if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) {
2194 DoLog(1) && (Log() << Verbose(1) << "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl);
2195 } else {
2196 DoLog(1) && (Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl);
2197 }
2198 }
2199 } else {
2200 DoeLog(0) && (eLog() << Verbose(1) << "REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius - otherradius) << endl);
2201 }
2202 } else {
2203 DoLog(1) && (Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl);
2204 }
2205 }
2206 catch (LinearDependenceException &excp){
2207 Log() << Verbose(1) << excp;
2208 Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;
2209 }
2210 } else {
2211 if (ThirdPoint != NULL) {
2212 DoLog(1) && (Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << "." << endl);
2213 } else {
2214 DoLog(1) && (Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << "." << endl);
2215 }
2216 }
2217 }
2218 }
2219 }
2220 } else {
2221 DoeLog(1) && (eLog() << Verbose(1) << "The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl);
2222 }
2223 } else {
2224 if (ThirdPoint != NULL)
2225 DoLog(1) && (Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!" << endl);
2226 else
2227 DoLog(1) && (Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!" << endl);
2228 }
2229
2230 DoLog(1) && (Log() << Verbose(1) << "INFO: Sorting candidate list ..." << endl);
2231 if (CandidateLine.pointlist.size() > 1) {
2232 CandidateLine.pointlist.unique();
2233 CandidateLine.pointlist.sort(); //SortCandidates);
2234 }
2235
2236 if ((!CandidateLine.pointlist.empty()) && (!CandidateLine.CheckValidity(RADIUS, LC))) {
2237 DoeLog(0) && (eLog() << Verbose(0) << "There were other points contained in the rolling sphere as well!" << endl);
2238 performCriticalExit();
2239 }
2240}
2241;
2242
2243/** Finds the endpoint two lines are sharing.
2244 * \param *line1 first line
2245 * \param *line2 second line
2246 * \return point which is shared or NULL if none
2247 */
2248class BoundaryPointSet *Tesselation::GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const
2249{
2250 Info FunctionInfo(__func__);
2251 const BoundaryLineSet * lines[2] = { line1, line2 };
2252 class BoundaryPointSet *node = NULL;
2253 PointMap OrderMap;
2254 PointTestPair OrderTest;
2255 for (int i = 0; i < 2; i++)
2256 // for both lines
2257 for (int j = 0; j < 2; j++) { // for both endpoints
2258 OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j]));
2259 if (!OrderTest.second) { // if insertion fails, we have common endpoint
2260 node = OrderTest.first->second;
2261 DoLog(1) && (Log() << Verbose(1) << "Common endpoint of lines " << *line1 << " and " << *line2 << " is: " << *node << "." << endl);
2262 j = 2;
2263 i = 2;
2264 break;
2265 }
2266 }
2267 return node;
2268}
2269;
2270
2271/** Finds the boundary points that are closest to a given Vector \a *x.
2272 * \param *out output stream for debugging
2273 * \param *x Vector to look from
2274 * \return map of BoundaryPointSet of closest points sorted by squared distance or NULL.
2275 */
2276DistanceToPointMap * Tesselation::FindClosestBoundaryPointsToVector(const Vector &x, const LinkedCell* LC) const
2277{
2278 Info FunctionInfo(__func__);
2279 PointMap::const_iterator FindPoint;
2280 int N[NDIM], Nlower[NDIM], Nupper[NDIM];
2281
2282 if (LinesOnBoundary.empty()) {
2283 DoeLog(1) && (eLog() << Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl);
2284 return NULL;
2285 }
2286
2287 // gather all points close to the desired one
2288 LC->SetIndexToVector(x); // ignore status as we calculate bounds below sensibly
2289 for (int i = 0; i < NDIM; i++) // store indices of this cell
2290 N[i] = LC->n[i];
2291 DoLog(1) && (Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl);
2292 DistanceToPointMap * points = new DistanceToPointMap;
2293 LC->GetNeighbourBounds(Nlower, Nupper);
2294 //Log() << Verbose(1) << endl;
2295 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
2296 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
2297 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
2298 const LinkedCell::LinkedNodes *List = LC->GetCurrentCell();
2299 //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
2300 if (List != NULL) {
2301 for (LinkedCell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
2302 FindPoint = PointsOnBoundary.find((*Runner)->nr);
2303 if (FindPoint != PointsOnBoundary.end()) {
2304 points->insert(DistanceToPointPair(FindPoint->second->node->DistanceSquared(x), FindPoint->second));
2305 DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl);
2306 }
2307 }
2308 } else {
2309 DoeLog(1) && (eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);
2310 }
2311 }
2312
2313 // check whether we found some points
2314 if (points->empty()) {
2315 DoeLog(1) && (eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
2316 delete (points);
2317 return NULL;
2318 }
2319 return points;
2320}
2321;
2322
2323/** Finds the boundary line that is closest to a given Vector \a *x.
2324 * \param *out output stream for debugging
2325 * \param *x Vector to look from
2326 * \return closest BoundaryLineSet or NULL in degenerate case.
2327 */
2328BoundaryLineSet * Tesselation::FindClosestBoundaryLineToVector(const Vector &x, const LinkedCell* LC) const
2329{
2330 Info FunctionInfo(__func__);
2331 // get closest points
2332 DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x, LC);
2333 if (points == NULL) {
2334 DoeLog(1) && (eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
2335 return NULL;
2336 }
2337
2338 // for each point, check its lines, remember closest
2339 DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryLine to " << x << " ... " << endl);
2340 BoundaryLineSet *ClosestLine = NULL;
2341 double MinDistance = -1.;
2342 Vector helper;
2343 Vector Center;
2344 Vector BaseLine;
2345 for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
2346 for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
2347 // calculate closest point on line to desired point
2348 helper = 0.5 * (((LineRunner->second)->endpoints[0]->node->getPosition()) +
2349 ((LineRunner->second)->endpoints[1]->node->getPosition()));
2350 Center = (x) - helper;
2351 BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
2352 ((LineRunner->second)->endpoints[1]->node->getPosition());
2353 Center.ProjectOntoPlane(BaseLine);
2354 const double distance = Center.NormSquared();
2355 if ((ClosestLine == NULL) || (distance < MinDistance)) {
2356 // additionally calculate intersection on line (whether it's on the line section or not)
2357 helper = (x) - ((LineRunner->second)->endpoints[0]->node->getPosition()) - Center;
2358 const double lengthA = helper.ScalarProduct(BaseLine);
2359 helper = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition()) - Center;
2360 const double lengthB = helper.ScalarProduct(BaseLine);
2361 if (lengthB * lengthA < 0) { // if have different sign
2362 ClosestLine = LineRunner->second;
2363 MinDistance = distance;
2364 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl);
2365 } else {
2366 DoLog(1) && (Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl);
2367 }
2368 } else {
2369 DoLog(1) && (Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl);
2370 }
2371 }
2372 }
2373 delete (points);
2374 // check whether closest line is "too close" :), then it's inside
2375 if (ClosestLine == NULL) {
2376 DoLog(0) && (Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl);
2377 return NULL;
2378 }
2379 return ClosestLine;
2380}
2381;
2382
2383/** Finds the triangle that is closest to a given Vector \a *x.
2384 * \param *out output stream for debugging
2385 * \param *x Vector to look from
2386 * \return BoundaryTriangleSet of nearest triangle or NULL.
2387 */
2388TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector &x, const LinkedCell* LC) const
2389{
2390 Info FunctionInfo(__func__);
2391 // get closest points
2392 DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x, LC);
2393 if (points == NULL) {
2394 DoeLog(1) && (eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl);
2395 return NULL;
2396 }
2397
2398 // for each point, check its lines, remember closest
2399 DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << x << " ... " << endl);
2400 LineSet ClosestLines;
2401 double MinDistance = 1e+16;
2402 Vector BaseLineIntersection;
2403 Vector Center;
2404 Vector BaseLine;
2405 Vector BaseLineCenter;
2406 for (DistanceToPointMap::iterator Runner = points->begin(); Runner != points->end(); Runner++) {
2407 for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) {
2408
2409 BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) -
2410 ((LineRunner->second)->endpoints[1]->node->getPosition());
2411 const double lengthBase = BaseLine.NormSquared();
2412
2413 BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[0]->node->getPosition());
2414 const double lengthEndA = BaseLineIntersection.NormSquared();
2415
2416 BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition());
2417 const double lengthEndB = BaseLineIntersection.NormSquared();
2418
2419 if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) { // intersection would be outside, take closer endpoint
2420 const double lengthEnd = std::min(lengthEndA, lengthEndB);
2421 if (lengthEnd - MinDistance < -MYEPSILON) { // new best line
2422 ClosestLines.clear();
2423 ClosestLines.insert(LineRunner->second);
2424 MinDistance = lengthEnd;
2425 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << "." << endl);
2426 } else if (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate
2427 ClosestLines.insert(LineRunner->second);
2428 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl);
2429 } else { // line is worse
2430 DoLog(1) && (Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl);
2431 }
2432 } else { // intersection is closer, calculate
2433 // calculate closest point on line to desired point
2434 BaseLineIntersection = (x) - ((LineRunner->second)->endpoints[1]->node->getPosition());
2435 Center = BaseLineIntersection;
2436 Center.ProjectOntoPlane(BaseLine);
2437 BaseLineIntersection -= Center;
2438 const double distance = BaseLineIntersection.NormSquared();
2439 if (Center.NormSquared() > BaseLine.NormSquared()) {
2440 DoeLog(0) && (eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl);
2441 }
2442 if ((ClosestLines.empty()) || (distance < MinDistance)) {
2443 ClosestLines.insert(LineRunner->second);
2444 MinDistance = distance;
2445 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " with projected distance " << MinDistance << "." << endl);
2446 } else {
2447 DoLog(2) && (Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl);
2448 }
2449 }
2450 }
2451 }
2452 delete (points);
2453
2454 // check whether closest line is "too close" :), then it's inside
2455 if (ClosestLines.empty()) {
2456 DoLog(0) && (Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl);
2457 return NULL;
2458 }
2459 TriangleList * candidates = new TriangleList;
2460 for (LineSet::iterator LineRunner = ClosestLines.begin(); LineRunner != ClosestLines.end(); LineRunner++)
2461 for (TriangleMap::iterator Runner = (*LineRunner)->triangles.begin(); Runner != (*LineRunner)->triangles.end(); Runner++) {
2462 candidates->push_back(Runner->second);
2463 }
2464 return candidates;
2465}
2466;
2467
2468/** Finds closest triangle to a point.
2469 * This basically just takes care of the degenerate case, which is not handled in FindClosestTrianglesToPoint().
2470 * \param *out output stream for debugging
2471 * \param *x Vector to look from
2472 * \param &distance contains found distance on return
2473 * \return list of BoundaryTriangleSet of nearest triangles or NULL.
2474 */
2475class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector &x, const LinkedCell* LC) const
2476{
2477 Info FunctionInfo(__func__);
2478 class BoundaryTriangleSet *result = NULL;
2479 TriangleList *triangles = FindClosestTrianglesToVector(x, LC);
2480 TriangleList candidates;
2481 Vector Center;
2482 Vector helper;
2483
2484 if ((triangles == NULL) || (triangles->empty()))
2485 return NULL;
2486
2487 // go through all and pick the one with the best alignment to x
2488 double MinAlignment = 2. * M_PI;
2489 for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) {
2490 (*Runner)->GetCenter(Center);
2491 helper = (x) - Center;
2492 const double Alignment = helper.Angle((*Runner)->NormalVector);
2493 if (Alignment < MinAlignment) {
2494 result = *Runner;
2495 MinAlignment = Alignment;
2496 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl);
2497 } else {
2498 DoLog(1) && (Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl);
2499 }
2500 }
2501 delete (triangles);
2502
2503 return result;
2504}
2505;
2506
2507/** Checks whether the provided Vector is within the Tesselation structure.
2508 * Basically calls Tesselation::GetDistanceToSurface() and checks the sign of the return value.
2509 * @param point of which to check the position
2510 * @param *LC LinkedCell structure
2511 *
2512 * @return true if the point is inside the Tesselation structure, false otherwise
2513 */
2514bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const
2515{
2516 Info FunctionInfo(__func__);
2517 TriangleIntersectionList Intersections(Point, this, LC);
2518
2519 return Intersections.IsInside();
2520}
2521;
2522
2523/** Returns the distance to the surface given by the tesselation.
2524 * Calls FindClosestTriangleToVector() and checks whether the resulting triangle's BoundaryTriangleSet#NormalVector points
2525 * towards or away from the given \a &Point. Additionally, we check whether it's normal to the normal vector, i.e. on the
2526 * closest triangle's plane. Then, we have to check whether \a Point is inside the triangle or not to determine whether it's
2527 * an inside or outside point. This is done by calling BoundaryTriangleSet::GetIntersectionInsideTriangle().
2528 * In the end we additionally find the point on the triangle who was smallest distance to \a Point:
2529 * -# Separate distance from point to center in vector in NormalDirection and on the triangle plane.
2530 * -# Check whether vector on triangle plane points inside the triangle or crosses triangle bounds.
2531 * -# If inside, take it to calculate closest distance
2532 * -# If not, take intersection with BoundaryLine as distance
2533 *
2534 * @note distance is squared despite it still contains a sign to determine in-/outside!
2535 *
2536 * @param point of which to check the position
2537 * @param *LC LinkedCell structure
2538 *
2539 * @return >0 if outside, ==0 if on surface, <0 if inside
2540 */
2541double Tesselation::GetDistanceSquaredToTriangle(const Vector &Point, const BoundaryTriangleSet* const triangle) const
2542{
2543 Info FunctionInfo(__func__);
2544 Vector Center;
2545 Vector helper;
2546 Vector DistanceToCenter;
2547 Vector Intersection;
2548 double distance = 0.;
2549
2550 if (triangle == NULL) {// is boundary point or only point in point cloud?
2551 DoLog(1) && (Log() << Verbose(1) << "No triangle given!" << endl);
2552 return -1.;
2553 } else {
2554 DoLog(1) && (Log() << Verbose(1) << "INFO: Closest triangle found is " << *triangle << " with normal vector " << triangle->NormalVector << "." << endl);
2555 }
2556
2557 triangle->GetCenter(Center);
2558 DoLog(2) && (Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl);
2559 DistanceToCenter = Center - Point;
2560 DoLog(2) && (Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl);
2561
2562 // check whether we are on boundary
2563 if (fabs(DistanceToCenter.ScalarProduct(triangle->NormalVector)) < MYEPSILON) {
2564 // calculate whether inside of triangle
2565 DistanceToCenter = Point + triangle->NormalVector; // points outside
2566 Center = Point - triangle->NormalVector; // points towards MolCenter
2567 DoLog(1) && (Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl);
2568 if (triangle->GetIntersectionInsideTriangle(Center, DistanceToCenter, Intersection)) {
2569 DoLog(1) && (Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl);
2570 return 0.;
2571 } else {
2572 DoLog(1) && (Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl);
2573 return false;
2574 }
2575 } else {
2576 // calculate smallest distance
2577 distance = triangle->GetClosestPointInsideTriangle(Point, Intersection);
2578 DoLog(1) && (Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl);
2579
2580 // then check direction to boundary
2581 if (DistanceToCenter.ScalarProduct(triangle->NormalVector) > MYEPSILON) {
2582 DoLog(1) && (Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl);
2583 return -distance;
2584 } else {
2585 DoLog(1) && (Log() << Verbose(1) << Point << " is NOT an inner point, " << distance << " above surface." << endl);
2586 return +distance;
2587 }
2588 }
2589}
2590;
2591
2592/** Calculates minimum distance from \a&Point to a tesselated surface.
2593 * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle().
2594 * \param &Point point to calculate distance from
2595 * \param *LC needed for finding closest points fast
2596 * \return distance squared to closest point on surface
2597 */
2598double Tesselation::GetDistanceToSurface(const Vector &Point, const LinkedCell* const LC) const
2599{
2600 Info FunctionInfo(__func__);
2601 TriangleIntersectionList Intersections(Point, this, LC);
2602
2603 return Intersections.GetSmallestDistance();
2604}
2605;
2606
2607/** Calculates minimum distance from \a&Point to a tesselated surface.
2608 * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle().
2609 * \param &Point point to calculate distance from
2610 * \param *LC needed for finding closest points fast
2611 * \return distance squared to closest point on surface
2612 */
2613BoundaryTriangleSet * Tesselation::GetClosestTriangleOnSurface(const Vector &Point, const LinkedCell* const LC) const
2614{
2615 Info FunctionInfo(__func__);
2616 TriangleIntersectionList Intersections(Point, this, LC);
2617
2618 return Intersections.GetClosestTriangle();
2619}
2620;
2621
2622/** Gets all points connected to the provided point by triangulation lines.
2623 *
2624 * @param *Point of which get all connected points
2625 *
2626 * @return set of the all points linked to the provided one
2627 */
2628TesselPointSet * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const
2629{
2630 Info FunctionInfo(__func__);
2631 TesselPointSet *connectedPoints = new TesselPointSet;
2632 class BoundaryPointSet *ReferencePoint = NULL;
2633 TesselPoint* current;
2634 bool takePoint = false;
2635 // find the respective boundary point
2636 PointMap::const_iterator PointRunner = PointsOnBoundary.find(Point->nr);
2637 if (PointRunner != PointsOnBoundary.end()) {
2638 ReferencePoint = PointRunner->second;
2639 } else {
2640 DoeLog(2) && (eLog() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);
2641 ReferencePoint = NULL;
2642 }
2643
2644 // little trick so that we look just through lines connect to the BoundaryPoint
2645 // OR fall-back to look through all lines if there is no such BoundaryPoint
2646 const LineMap *Lines;
2647 ;
2648 if (ReferencePoint != NULL)
2649 Lines = &(ReferencePoint->lines);
2650 else
2651 Lines = &LinesOnBoundary;
2652 LineMap::const_iterator findLines = Lines->begin();
2653 while (findLines != Lines->end()) {
2654 takePoint = false;
2655
2656 if (findLines->second->endpoints[0]->Nr == Point->nr) {
2657 takePoint = true;
2658 current = findLines->second->endpoints[1]->node;
2659 } else if (findLines->second->endpoints[1]->Nr == Point->nr) {
2660 takePoint = true;
2661 current = findLines->second->endpoints[0]->node;
2662 }
2663
2664 if (takePoint) {
2665 DoLog(1) && (Log() << Verbose(1) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl);
2666 connectedPoints->insert(current);
2667 }
2668
2669 findLines++;
2670 }
2671
2672 if (connectedPoints->empty()) { // if have not found any points
2673 DoeLog(1) && (eLog() << Verbose(1) << "We have not found any connected points to " << *Point << "." << endl);
2674 return NULL;
2675 }
2676
2677 return connectedPoints;
2678}
2679;
2680
2681/** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
2682 * Maps them down onto the plane designated by the axis \a *Point and \a *Reference. The center of all points
2683 * connected in the tesselation to \a *Point is mapped to spherical coordinates with the zero angle being given
2684 * by the mapped down \a *Reference. Hence, the biggest and the smallest angles are those of the two shanks of the
2685 * triangle we are looking for.
2686 *
2687 * @param *out output stream for debugging
2688 * @param *SetOfNeighbours all points for which the angle should be calculated
2689 * @param *Point of which get all connected points
2690 * @param *Reference Reference vector for zero angle or NULL for no preference
2691 * @return list of the all points linked to the provided one
2692 */
2693TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector &Reference) const
2694{
2695 Info FunctionInfo(__func__);
2696 map<double, TesselPoint*> anglesOfPoints;
2697 TesselPointList *connectedCircle = new TesselPointList;
2698 Vector PlaneNormal;
2699 Vector AngleZero;
2700 Vector OrthogonalVector;
2701 Vector helper;
2702 const TesselPoint * const TrianglePoints[3] = { Point, NULL, NULL };
2703 TriangleList *triangles = NULL;
2704
2705 if (SetOfNeighbours == NULL) {
2706 DoeLog(2) && (eLog() << Verbose(2) << "Could not find any connected points!" << endl);
2707 delete (connectedCircle);
2708 return NULL;
2709 }
2710
2711 // calculate central point
2712 triangles = FindTriangles(TrianglePoints);
2713 if ((triangles != NULL) && (!triangles->empty())) {
2714 for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++)
2715 PlaneNormal += (*Runner)->NormalVector;
2716 } else {
2717 DoeLog(0) && (eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl);
2718 performCriticalExit();
2719 }
2720 PlaneNormal.Scale(1.0 / triangles->size());
2721 DoLog(1) && (Log() << Verbose(1) << "INFO: Calculated PlaneNormal of all circle points is " << PlaneNormal << "." << endl);
2722 PlaneNormal.Normalize();
2723
2724 // construct one orthogonal vector
2725 AngleZero = (Reference) - (Point->getPosition());
2726 AngleZero.ProjectOntoPlane(PlaneNormal);
2727 if ((AngleZero.NormSquared() < MYEPSILON)) {
2728 DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << (*SetOfNeighbours->begin())->getPosition() << " as angle 0 referencer." << endl);
2729 AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
2730 AngleZero.ProjectOntoPlane(PlaneNormal);
2731 if (AngleZero.NormSquared() < MYEPSILON) {
2732 DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
2733 performCriticalExit();
2734 }
2735 }
2736 DoLog(1) && (Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl);
2737 if (AngleZero.NormSquared() > MYEPSILON)
2738 OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
2739 else
2740 OrthogonalVector.MakeNormalTo(PlaneNormal);
2741 DoLog(1) && (Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl);
2742
2743 // go through all connected points and calculate angle
2744 for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
2745 helper = ((*listRunner)->getPosition()) - (Point->getPosition());
2746 helper.ProjectOntoPlane(PlaneNormal);
2747 double angle = GetAngle(helper, AngleZero, OrthogonalVector);
2748 DoLog(0) && (Log() << Verbose(0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl);
2749 anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner)));
2750 }
2751
2752 for (map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
2753 connectedCircle->push_back(AngleRunner->second);
2754 }
2755
2756 return connectedCircle;
2757}
2758
2759/** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point.
2760 * Maps them down onto the plane designated by the axis \a *Point and \a *Reference. The center of all points
2761 * connected in the tesselation to \a *Point is mapped to spherical coordinates with the zero angle being given
2762 * by the mapped down \a *Reference. Hence, the biggest and the smallest angles are those of the two shanks of the
2763 * triangle we are looking for.
2764 *
2765 * @param *SetOfNeighbours all points for which the angle should be calculated
2766 * @param *Point of which get all connected points
2767 * @param *Reference Reference vector for zero angle or (0,0,0) for no preference
2768 * @return list of the all points linked to the provided one
2769 */
2770TesselPointList * Tesselation::GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector &Reference) const
2771{
2772 Info FunctionInfo(__func__);
2773 map<double, TesselPoint*> anglesOfPoints;
2774 TesselPointList *connectedCircle = new TesselPointList;
2775 Vector center;
2776 Vector PlaneNormal;
2777 Vector AngleZero;
2778 Vector OrthogonalVector;
2779 Vector helper;
2780
2781 if (SetOfNeighbours == NULL) {
2782 DoeLog(2) && (eLog() << Verbose(2) << "Could not find any connected points!" << endl);
2783 delete (connectedCircle);
2784 return NULL;
2785 }
2786
2787 // check whether there's something to do
2788 if (SetOfNeighbours->size() < 3) {
2789 for (TesselPointSet::iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)
2790 connectedCircle->push_back(*TesselRunner);
2791 return connectedCircle;
2792 }
2793
2794 DoLog(1) && (Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << Reference << "." << endl);
2795 // calculate central point
2796 TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin();
2797 TesselPointSet::const_iterator TesselB = SetOfNeighbours->begin();
2798 TesselPointSet::const_iterator TesselC = SetOfNeighbours->begin();
2799 TesselB++;
2800 TesselC++;
2801 TesselC++;
2802 int counter = 0;
2803 while (TesselC != SetOfNeighbours->end()) {
2804 helper = Plane(((*TesselA)->getPosition()),
2805 ((*TesselB)->getPosition()),
2806 ((*TesselC)->getPosition())).getNormal();
2807 DoLog(0) && (Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl);
2808 counter++;
2809 TesselA++;
2810 TesselB++;
2811 TesselC++;
2812 PlaneNormal += helper;
2813 }
2814 //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size()
2815 // << "; scale factor " << counter;
2816 PlaneNormal.Scale(1.0 / (double) counter);
2817 // Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;
2818 //
2819 // // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points
2820 // PlaneNormal.CopyVector(Point->node);
2821 // PlaneNormal.SubtractVector(&center);
2822 // PlaneNormal.Normalize();
2823 DoLog(1) && (Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl);
2824
2825 // construct one orthogonal vector
2826 if (!Reference.IsZero()) {
2827 AngleZero = (Reference) - (Point->getPosition());
2828 AngleZero.ProjectOntoPlane(PlaneNormal);
2829 }
2830 if ((Reference.IsZero()) || (AngleZero.NormSquared() < MYEPSILON )) {
2831 DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << (*SetOfNeighbours->begin())->getPosition() << " as angle 0 referencer." << endl);
2832 AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition());
2833 AngleZero.ProjectOntoPlane(PlaneNormal);
2834 if (AngleZero.NormSquared() < MYEPSILON) {
2835 DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);
2836 performCriticalExit();
2837 }
2838 }
2839 DoLog(1) && (Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl);
2840 if (AngleZero.NormSquared() > MYEPSILON)
2841 OrthogonalVector = Plane(PlaneNormal, AngleZero,0).getNormal();
2842 else
2843 OrthogonalVector.MakeNormalTo(PlaneNormal);
2844 DoLog(1) && (Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl);
2845
2846 // go through all connected points and calculate angle
2847 pair<map<double, TesselPoint*>::iterator, bool> InserterTest;
2848 for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {
2849 helper = ((*listRunner)->getPosition()) - (Point->getPosition());
2850 helper.ProjectOntoPlane(PlaneNormal);
2851 double angle = GetAngle(helper, AngleZero, OrthogonalVector);
2852 if (angle > M_PI) // the correction is of no use here (and not desired)
2853 angle = 2. * M_PI - angle;
2854 DoLog(0) && (Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl);
2855 InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner)));
2856 if (!InserterTest.second) {
2857 DoeLog(0) && (eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl);
2858 performCriticalExit();
2859 }
2860 }
2861
2862 for (map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {
2863 connectedCircle->push_back(AngleRunner->second);
2864 }
2865
2866 return connectedCircle;
2867}
2868
2869/** Gets all points connected to the provided point by triangulation lines, ordered such that we walk along a closed path.
2870 *
2871 * @param *out output stream for debugging
2872 * @param *Point of which get all connected points
2873 * @return list of the all points linked to the provided one
2874 */
2875ListOfTesselPointList * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const
2876{
2877 Info FunctionInfo(__func__);
2878 map<double, TesselPoint*> anglesOfPoints;
2879 list<TesselPointList *> *ListOfPaths = new list<TesselPointList *> ;
2880 TesselPointList *connectedPath = NULL;
2881 Vector center;
2882 Vector PlaneNormal;
2883 Vector AngleZero;
2884 Vector OrthogonalVector;
2885 Vector helper;
2886 class BoundaryPointSet *ReferencePoint = NULL;
2887 class BoundaryPointSet *CurrentPoint = NULL;
2888 class BoundaryTriangleSet *triangle = NULL;
2889 class BoundaryLineSet *CurrentLine = NULL;
2890 class BoundaryLineSet *StartLine = NULL;
2891 // find the respective boundary point
2892 PointMap::const_iterator PointRunner = PointsOnBoundary.find(Point->nr);
2893 if (PointRunner != PointsOnBoundary.end()) {
2894 ReferencePoint = PointRunner->second;
2895 } else {
2896 DoeLog(1) && (eLog() << Verbose(1) << "GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);
2897 return NULL;
2898 }
2899
2900 map<class BoundaryLineSet *, bool> TouchedLine;
2901 map<class BoundaryTriangleSet *, bool> TouchedTriangle;
2902 map<class BoundaryLineSet *, bool>::iterator LineRunner;
2903 map<class BoundaryTriangleSet *, bool>::iterator TriangleRunner;
2904 for (LineMap::iterator Runner = ReferencePoint->lines.begin(); Runner != ReferencePoint->lines.end(); Runner++) {
2905 TouchedLine.insert(pair<class BoundaryLineSet *, bool> (Runner->second, false));
2906 for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++)
2907 TouchedTriangle.insert(pair<class BoundaryTriangleSet *, bool> (Sprinter->second, false));
2908 }
2909 if (!ReferencePoint->lines.empty()) {
2910 for (LineMap::iterator runner = ReferencePoint->lines.begin(); runner != ReferencePoint->lines.end(); runner++) {
2911 LineRunner = TouchedLine.find(runner->second);
2912 if (LineRunner == TouchedLine.end()) {
2913 DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *runner->second << " in the touched list." << endl);
2914 } else if (!LineRunner->second) {
2915 LineRunner->second = true;
2916 connectedPath = new TesselPointList;
2917 triangle = NULL;
2918 CurrentLine = runner->second;
2919 StartLine = CurrentLine;
2920 CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
2921 DoLog(1) && (Log() << Verbose(1) << "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl);
2922 do {
2923 // push current one
2924 DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl);
2925 connectedPath->push_back(CurrentPoint->node);
2926
2927 // find next triangle
2928 for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) {
2929 DoLog(1) && (Log() << Verbose(1) << "INFO: Inspecting triangle " << *Runner->second << "." << endl);
2930 if ((Runner->second != triangle)) { // look for first triangle not equal to old one
2931 triangle = Runner->second;
2932 TriangleRunner = TouchedTriangle.find(triangle);
2933 if (TriangleRunner != TouchedTriangle.end()) {
2934 if (!TriangleRunner->second) {
2935 TriangleRunner->second = true;
2936 DoLog(1) && (Log() << Verbose(1) << "INFO: Connecting triangle is " << *triangle << "." << endl);
2937 break;
2938 } else {
2939 DoLog(1) && (Log() << Verbose(1) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl);
2940 triangle = NULL;
2941 }
2942 } else {
2943 DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *triangle << " in the touched list." << endl);
2944 triangle = NULL;
2945 }
2946 }
2947 }
2948 if (triangle == NULL)
2949 break;
2950 // find next line
2951 for (int i = 0; i < 3; i++) {
2952 if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point
2953 CurrentLine = triangle->lines[i];
2954 DoLog(1) && (Log() << Verbose(1) << "INFO: Connecting line is " << *CurrentLine << "." << endl);
2955 break;
2956 }
2957 }
2958 LineRunner = TouchedLine.find(CurrentLine);
2959 if (LineRunner == TouchedLine.end())
2960 DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *CurrentLine << " in the touched list." << endl);
2961 else
2962 LineRunner->second = true;
2963 // find next point
2964 CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint);
2965
2966 } while (CurrentLine != StartLine);
2967 // last point is missing, as it's on start line
2968 DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl);
2969 if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back())
2970 connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node);
2971
2972 ListOfPaths->push_back(connectedPath);
2973 } else {
2974 DoLog(1) && (Log() << Verbose(1) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl);
2975 }
2976 }
2977 } else {
2978 DoeLog(1) && (eLog() << Verbose(1) << "There are no lines attached to " << *ReferencePoint << "." << endl);
2979 }
2980
2981 return ListOfPaths;
2982}
2983
2984/** Gets all closed paths on the circle of points connected to the provided point by triangulation lines, if this very point is removed.
2985 * From GetPathsOfConnectedPoints() extracts all single loops of intracrossing paths in the list of closed paths.
2986 * @param *out output stream for debugging
2987 * @param *Point of which get all connected points
2988 * @return list of the closed paths
2989 */
2990ListOfTesselPointList * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const
2991{
2992 Info FunctionInfo(__func__);
2993 list<TesselPointList *> *ListofPaths = GetPathsOfConnectedPoints(Point);
2994 list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *> ;
2995 TesselPointList *connectedPath = NULL;
2996 TesselPointList *newPath = NULL;
2997 int count = 0;
2998 TesselPointList::iterator CircleRunner;
2999 TesselPointList::iterator CircleStart;
3000
3001 for (list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {
3002 connectedPath = *ListRunner;
3003
3004 DoLog(1) && (Log() << Verbose(1) << "INFO: Current path is " << connectedPath << "." << endl);
3005
3006 // go through list, look for reappearance of starting Point and count
3007 CircleStart = connectedPath->begin();
3008 // go through list, look for reappearance of starting Point and create list
3009 TesselPointList::iterator Marker = CircleStart;
3010 for (CircleRunner = CircleStart; CircleRunner != connectedPath->end(); CircleRunner++) {
3011 if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point
3012 // we have a closed circle from Marker to new Marker
3013 DoLog(1) && (Log() << Verbose(1) << count + 1 << ". closed path consists of: ");
3014 newPath = new TesselPointList;
3015 TesselPointList::iterator CircleSprinter = Marker;
3016 for (; CircleSprinter != CircleRunner; CircleSprinter++) {
3017 newPath->push_back(*CircleSprinter);
3018 DoLog(0) && (Log() << Verbose(0) << (**CircleSprinter) << " <-> ");
3019 }
3020 DoLog(0) && (Log() << Verbose(0) << ".." << endl);
3021 count++;
3022 Marker = CircleRunner;
3023
3024 // add to list
3025 ListofClosedPaths->push_back(newPath);
3026 }
3027 }
3028 }
3029 DoLog(1) && (Log() << Verbose(1) << "INFO: " << count << " closed additional path(s) have been created." << endl);
3030
3031 // delete list of paths
3032 while (!ListofPaths->empty()) {
3033 connectedPath = *(ListofPaths->begin());
3034 ListofPaths->remove(connectedPath);
3035 delete (connectedPath);
3036 }
3037 delete (ListofPaths);
3038
3039 // exit
3040 return ListofClosedPaths;
3041}
3042;
3043
3044/** Gets all belonging triangles for a given BoundaryPointSet.
3045 * \param *out output stream for debugging
3046 * \param *Point BoundaryPoint
3047 * \return pointer to allocated list of triangles
3048 */
3049TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const
3050{
3051 Info FunctionInfo(__func__);
3052 TriangleSet *connectedTriangles = new TriangleSet;
3053
3054 if (Point == NULL) {
3055 DoeLog(1) && (eLog() << Verbose(1) << "Point given is NULL." << endl);
3056 } else {
3057 // go through its lines and insert all triangles
3058 for (LineMap::const_iterator LineRunner = Point->lines.begin(); LineRunner != Point->lines.end(); LineRunner++)
3059 for (TriangleMap::iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
3060 connectedTriangles->insert(TriangleRunner->second);
3061 }
3062 }
3063
3064 return connectedTriangles;
3065}
3066;
3067
3068/** Removes a boundary point from the envelope while keeping it closed.
3069 * We remove the old triangles connected to the point and re-create new triangles to close the surface following this ansatz:
3070 * -# a closed path(s) of boundary points surrounding the point to be removed is constructed
3071 * -# on each closed path, we pick three adjacent points, create a triangle with them and subtract the middle point from the path
3072 * -# we advance two points (i.e. the next triangle will start at the ending point of the last triangle) and continue as before
3073 * -# the surface is closed, when the path is empty
3074 * Thereby, we (hopefully) make sure that the removed points remains beneath the surface (this is checked via IsInnerPoint eventually).
3075 * \param *out output stream for debugging
3076 * \param *point point to be removed
3077 * \return volume added to the volume inside the tesselated surface by the removal
3078 */
3079double Tesselation::RemovePointFromTesselatedSurface(class BoundaryPointSet *point)
3080{
3081 class BoundaryLineSet *line = NULL;
3082 class BoundaryTriangleSet *triangle = NULL;
3083 Vector OldPoint, NormalVector;
3084 double volume = 0;
3085 int count = 0;
3086
3087 if (point == NULL) {
3088 DoeLog(1) && (eLog() << Verbose(1) << "Cannot remove the point " << point << ", it's NULL!" << endl);
3089 return 0.;
3090 } else
3091 DoLog(0) && (Log() << Verbose(0) << "Removing point " << *point << " from tesselated boundary ..." << endl);
3092
3093 // copy old location for the volume
3094 OldPoint = (point->node->getPosition());
3095
3096 // get list of connected points
3097 if (point->lines.empty()) {
3098 DoeLog(1) && (eLog() << Verbose(1) << "Cannot remove the point " << *point << ", it's connected to no lines!" << endl);
3099 return 0.;
3100 }
3101
3102 list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
3103 TesselPointList *connectedPath = NULL;
3104
3105 // gather all triangles
3106 for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
3107 count += LineRunner->second->triangles.size();
3108 TriangleMap Candidates;
3109 for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
3110 line = LineRunner->second;
3111 for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
3112 triangle = TriangleRunner->second;
3113 Candidates.insert(TrianglePair(triangle->Nr, triangle));
3114 }
3115 }
3116
3117 // remove all triangles
3118 count = 0;
3119 NormalVector.Zero();
3120 for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
3121 DoLog(1) && (Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl);
3122 NormalVector -= Runner->second->NormalVector; // has to point inward
3123 RemoveTesselationTriangle(Runner->second);
3124 count++;
3125 }
3126 DoLog(1) && (Log() << Verbose(1) << count << " triangles were removed." << endl);
3127
3128 list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
3129 list<TesselPointList *>::iterator ListRunner = ListAdvance;
3130 TriangleMap::iterator NumberRunner = Candidates.begin();
3131 TesselPointList::iterator StartNode, MiddleNode, EndNode;
3132 double angle;
3133 double smallestangle;
3134 Vector Point, Reference, OrthogonalVector;
3135 if (count > 2) { // less than three triangles, then nothing will be created
3136 class TesselPoint *TriangleCandidates[3];
3137 count = 0;
3138 for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
3139 if (ListAdvance != ListOfClosedPaths->end())
3140 ListAdvance++;
3141
3142 connectedPath = *ListRunner;
3143 // re-create all triangles by going through connected points list
3144 LineList NewLines;
3145 for (; !connectedPath->empty();) {
3146 // search middle node with widest angle to next neighbours
3147 EndNode = connectedPath->end();
3148 smallestangle = 0.;
3149 for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
3150 DoLog(1) && (Log() << Verbose(1) << "INFO: MiddleNode is " << **MiddleNode << "." << endl);
3151 // construct vectors to next and previous neighbour
3152 StartNode = MiddleNode;
3153 if (StartNode == connectedPath->begin())
3154 StartNode = connectedPath->end();
3155 StartNode--;
3156 //Log() << Verbose(3) << "INFO: StartNode is " << **StartNode << "." << endl;
3157 Point = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
3158 StartNode = MiddleNode;
3159 StartNode++;
3160 if (StartNode == connectedPath->end())
3161 StartNode = connectedPath->begin();
3162 //Log() << Verbose(3) << "INFO: EndNode is " << **StartNode << "." << endl;
3163 Reference = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
3164 OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
3165 OrthogonalVector.MakeNormalTo(Reference);
3166 angle = GetAngle(Point, Reference, OrthogonalVector);
3167 //if (angle < M_PI) // no wrong-sided triangles, please?
3168 if (fabs(angle - M_PI) < fabs(smallestangle - M_PI)) { // get straightest angle (i.e. construct those triangles with smallest area first)
3169 smallestangle = angle;
3170 EndNode = MiddleNode;
3171 }
3172 }
3173 MiddleNode = EndNode;
3174 if (MiddleNode == connectedPath->end()) {
3175 DoeLog(0) && (eLog() << Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl);
3176 performCriticalExit();
3177 }
3178 StartNode = MiddleNode;
3179 if (StartNode == connectedPath->begin())
3180 StartNode = connectedPath->end();
3181 StartNode--;
3182 EndNode++;
3183 if (EndNode == connectedPath->end())
3184 EndNode = connectedPath->begin();
3185 DoLog(2) && (Log() << Verbose(2) << "INFO: StartNode is " << **StartNode << "." << endl);
3186 DoLog(2) && (Log() << Verbose(2) << "INFO: MiddleNode is " << **MiddleNode << "." << endl);
3187 DoLog(2) && (Log() << Verbose(2) << "INFO: EndNode is " << **EndNode << "." << endl);
3188 DoLog(1) && (Log() << Verbose(1) << "INFO: Attempting to create triangle " << (*StartNode)->getName() << ", " << (*MiddleNode)->getName() << " and " << (*EndNode)->getName() << "." << endl);
3189 TriangleCandidates[0] = *StartNode;
3190 TriangleCandidates[1] = *MiddleNode;
3191 TriangleCandidates[2] = *EndNode;
3192 triangle = GetPresentTriangle(TriangleCandidates);
3193 if (triangle != NULL) {
3194 DoeLog(0) && (eLog() << Verbose(0) << "New triangle already present, skipping!" << endl);
3195 StartNode++;
3196 MiddleNode++;
3197 EndNode++;
3198 if (StartNode == connectedPath->end())
3199 StartNode = connectedPath->begin();
3200 if (MiddleNode == connectedPath->end())
3201 MiddleNode = connectedPath->begin();
3202 if (EndNode == connectedPath->end())
3203 EndNode = connectedPath->begin();
3204 continue;
3205 }
3206 DoLog(3) && (Log() << Verbose(3) << "Adding new triangle points." << endl);
3207 AddTesselationPoint(*StartNode, 0);
3208 AddTesselationPoint(*MiddleNode, 1);
3209 AddTesselationPoint(*EndNode, 2);
3210 DoLog(3) && (Log() << Verbose(3) << "Adding new triangle lines." << endl);
3211 AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
3212 AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
3213 NewLines.push_back(BLS[1]);
3214 AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
3215 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
3216 BTS->GetNormalVector(NormalVector);
3217 AddTesselationTriangle();
3218 // calculate volume summand as a general tetraeder
3219 volume += CalculateVolumeofGeneralTetraeder(TPS[0]->node->getPosition(), TPS[1]->node->getPosition(), TPS[2]->node->getPosition(), OldPoint);
3220 // advance number
3221 count++;
3222
3223 // prepare nodes for next triangle
3224 StartNode = EndNode;
3225 DoLog(2) && (Log() << Verbose(2) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl);
3226 connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles)
3227 if (connectedPath->size() == 2) { // we are done
3228 connectedPath->remove(*StartNode); // remove the start node
3229 connectedPath->remove(*EndNode); // remove the end node
3230 break;
3231 } else if (connectedPath->size() < 2) { // something's gone wrong!
3232 DoeLog(0) && (eLog() << Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl);
3233 performCriticalExit();
3234 } else {
3235 MiddleNode = StartNode;
3236 MiddleNode++;
3237 if (MiddleNode == connectedPath->end())
3238 MiddleNode = connectedPath->begin();
3239 EndNode = MiddleNode;
3240 EndNode++;
3241 if (EndNode == connectedPath->end())
3242 EndNode = connectedPath->begin();
3243 }
3244 }
3245 // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing)
3246 if (NewLines.size() > 1) {
3247 LineList::iterator Candidate;
3248 class BoundaryLineSet *OtherBase = NULL;
3249 double tmp, maxgain;
3250 do {
3251 maxgain = 0;
3252 for (LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {
3253 tmp = PickFarthestofTwoBaselines(*Runner);
3254 if (maxgain < tmp) {
3255 maxgain = tmp;
3256 Candidate = Runner;
3257 }
3258 }
3259 if (maxgain != 0) {
3260 volume += maxgain;
3261 DoLog(1) && (Log() << Verbose(1) << "Flipping baseline with highest volume" << **Candidate << "." << endl);
3262 OtherBase = FlipBaseline(*Candidate);
3263 NewLines.erase(Candidate);
3264 NewLines.push_back(OtherBase);
3265 }
3266 } while (maxgain != 0.);
3267 }
3268
3269 ListOfClosedPaths->remove(connectedPath);
3270 delete (connectedPath);
3271 }
3272 DoLog(0) && (Log() << Verbose(0) << count << " triangles were created." << endl);
3273 } else {
3274 while (!ListOfClosedPaths->empty()) {
3275 ListRunner = ListOfClosedPaths->begin();
3276 connectedPath = *ListRunner;
3277 ListOfClosedPaths->remove(connectedPath);
3278 delete (connectedPath);
3279 }
3280 DoLog(0) && (Log() << Verbose(0) << "No need to create any triangles." << endl);
3281 }
3282 delete (ListOfClosedPaths);
3283
3284 DoLog(0) && (Log() << Verbose(0) << "Removed volume is " << volume << "." << endl);
3285
3286 return volume;
3287}
3288;
3289
3290/**
3291 * Finds triangles belonging to the three provided points.
3292 *
3293 * @param *Points[3] list, is expected to contain three points (NULL means wildcard)
3294 *
3295 * @return triangles which belong to the provided points, will be empty if there are none,
3296 * will usually be one, in case of degeneration, there will be two
3297 */
3298TriangleList *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const
3299{
3300 Info FunctionInfo(__func__);
3301 TriangleList *result = new TriangleList;
3302 LineMap::const_iterator FindLine;
3303 TriangleMap::const_iterator FindTriangle;
3304 class BoundaryPointSet *TrianglePoints[3];
3305 size_t NoOfWildcards = 0;
3306
3307 for (int i = 0; i < 3; i++) {
3308 if (Points[i] == NULL) {
3309 NoOfWildcards++;
3310 TrianglePoints[i] = NULL;
3311 } else {
3312 PointMap::const_iterator FindPoint = PointsOnBoundary.find(Points[i]->nr);
3313 if (FindPoint != PointsOnBoundary.end()) {
3314 TrianglePoints[i] = FindPoint->second;
3315 } else {
3316 TrianglePoints[i] = NULL;
3317 }
3318 }
3319 }
3320
3321 switch (NoOfWildcards) {
3322 case 0: // checks lines between the points in the Points for their adjacent triangles
3323 for (int i = 0; i < 3; i++) {
3324 if (TrianglePoints[i] != NULL) {
3325 for (int j = i + 1; j < 3; j++) {
3326 if (TrianglePoints[j] != NULL) {
3327 for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap!
3328 (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr); FindLine++) {
3329 for (FindTriangle = FindLine->second->triangles.begin(); FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
3330 if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
3331 result->push_back(FindTriangle->second);
3332 }
3333 }
3334 }
3335 // Is it sufficient to consider one of the triangle lines for this.
3336 return result;
3337 }
3338 }
3339 }
3340 }
3341 break;
3342 case 1: // copy all triangles of the respective line
3343 {
3344 int i = 0;
3345 for (; i < 3; i++)
3346 if (TrianglePoints[i] == NULL)
3347 break;
3348 for (FindLine = TrianglePoints[(i + 1) % 3]->lines.find(TrianglePoints[(i + 2) % 3]->node->nr); // is a multimap!
3349 (FindLine != TrianglePoints[(i + 1) % 3]->lines.end()) && (FindLine->first == TrianglePoints[(i + 2) % 3]->node->nr); FindLine++) {
3350 for (FindTriangle = FindLine->second->triangles.begin(); FindTriangle != FindLine->second->triangles.end(); FindTriangle++) {
3351 if (FindTriangle->second->IsPresentTupel(TrianglePoints)) {
3352 result->push_back(FindTriangle->second);
3353 }
3354 }
3355 }
3356 break;
3357 }
3358 case 2: // copy all triangles of the respective point
3359 {
3360 int i = 0;
3361 for (; i < 3; i++)
3362 if (TrianglePoints[i] != NULL)
3363 break;
3364 for (LineMap::const_iterator line = TrianglePoints[i]->lines.begin(); line != TrianglePoints[i]->lines.end(); line++)
3365 for (TriangleMap::const_iterator triangle = line->second->triangles.begin(); triangle != line->second->triangles.end(); triangle++)
3366 result->push_back(triangle->second);
3367 result->sort();
3368 result->unique();
3369 break;
3370 }
3371 case 3: // copy all triangles
3372 {
3373 for (TriangleMap::const_iterator triangle = TrianglesOnBoundary.begin(); triangle != TrianglesOnBoundary.end(); triangle++)
3374 result->push_back(triangle->second);
3375 break;
3376 }
3377 default:
3378 DoeLog(0) && (eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl);
3379 performCriticalExit();
3380 break;
3381 }
3382
3383 return result;
3384}
3385
3386struct BoundaryLineSetCompare
3387{
3388 bool operator()(const BoundaryLineSet * const a, const BoundaryLineSet * const b)
3389 {
3390 int lowerNra = -1;
3391 int lowerNrb = -1;
3392
3393 if (a->endpoints[0] < a->endpoints[1])
3394 lowerNra = 0;
3395 else
3396 lowerNra = 1;
3397
3398 if (b->endpoints[0] < b->endpoints[1])
3399 lowerNrb = 0;
3400 else
3401 lowerNrb = 1;
3402
3403 if (a->endpoints[lowerNra] < b->endpoints[lowerNrb])
3404 return true;
3405 else if (a->endpoints[lowerNra] > b->endpoints[lowerNrb])
3406 return false;
3407 else { // both lower-numbered endpoints are the same ...
3408 if (a->endpoints[(lowerNra + 1) % 2] < b->endpoints[(lowerNrb + 1) % 2])
3409 return true;
3410 else if (a->endpoints[(lowerNra + 1) % 2] > b->endpoints[(lowerNrb + 1) % 2])
3411 return false;
3412 }
3413 return false;
3414 }
3415 ;
3416};
3417
3418#define UniqueLines set < class BoundaryLineSet *, BoundaryLineSetCompare>
3419
3420/**
3421 * Finds all degenerated lines within the tesselation structure.
3422 *
3423 * @return map of keys of degenerated line pairs, each line occurs twice
3424 * in the list, once as key and once as value
3425 */
3426IndexToIndex * Tesselation::FindAllDegeneratedLines()
3427{
3428 Info FunctionInfo(__func__);
3429 UniqueLines AllLines;
3430 IndexToIndex * DegeneratedLines = new IndexToIndex;
3431
3432 // sanity check
3433 if (LinesOnBoundary.empty()) {
3434 DoeLog(2) && (eLog() << Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure.");
3435 return DegeneratedLines;
3436 }
3437 LineMap::iterator LineRunner1;
3438 pair<UniqueLines::iterator, bool> tester;
3439 for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) {
3440 tester = AllLines.insert(LineRunner1->second);
3441 if (!tester.second) { // found degenerated line
3442 DegeneratedLines->insert(pair<int, int> (LineRunner1->second->Nr, (*tester.first)->Nr));
3443 DegeneratedLines->insert(pair<int, int> ((*tester.first)->Nr, LineRunner1->second->Nr));
3444 }
3445 }
3446
3447 AllLines.clear();
3448
3449 DoLog(0) && (Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl);
3450 IndexToIndex::iterator it;
3451 for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) {
3452 const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first);
3453 const LineMap::const_iterator Line2 = LinesOnBoundary.find((*it).second);
3454 if (Line1 != LinesOnBoundary.end() && Line2 != LinesOnBoundary.end())
3455 DoLog(0) && (Log() << Verbose(0) << *Line1->second << " => " << *Line2->second << endl);
3456 else
3457 DoeLog(1) && (eLog() << Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl);
3458 }
3459
3460 return DegeneratedLines;
3461}
3462
3463/**
3464 * Finds all degenerated triangles within the tesselation structure.
3465 *
3466 * @return map of keys of degenerated triangle pairs, each triangle occurs twice
3467 * in the list, once as key and once as value
3468 */
3469IndexToIndex * Tesselation::FindAllDegeneratedTriangles()
3470{
3471 Info FunctionInfo(__func__);
3472 IndexToIndex * DegeneratedLines = FindAllDegeneratedLines();
3473 IndexToIndex * DegeneratedTriangles = new IndexToIndex;
3474 TriangleMap::iterator TriangleRunner1, TriangleRunner2;
3475 LineMap::iterator Liner;
3476 class BoundaryLineSet *line1 = NULL, *line2 = NULL;
3477
3478 for (IndexToIndex::iterator LineRunner = DegeneratedLines->begin(); LineRunner != DegeneratedLines->end(); ++LineRunner) {
3479 // run over both lines' triangles
3480 Liner = LinesOnBoundary.find(LineRunner->first);
3481 if (Liner != LinesOnBoundary.end())
3482 line1 = Liner->second;
3483 Liner = LinesOnBoundary.find(LineRunner->second);
3484 if (Liner != LinesOnBoundary.end())
3485 line2 = Liner->second;
3486 for (TriangleRunner1 = line1->triangles.begin(); TriangleRunner1 != line1->triangles.end(); ++TriangleRunner1) {
3487 for (TriangleRunner2 = line2->triangles.begin(); TriangleRunner2 != line2->triangles.end(); ++TriangleRunner2) {
3488 if ((TriangleRunner1->second != TriangleRunner2->second) && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) {
3489 DegeneratedTriangles->insert(pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr));
3490 DegeneratedTriangles->insert(pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr));
3491 }
3492 }
3493 }
3494 }
3495 delete (DegeneratedLines);
3496
3497 DoLog(0) && (Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl);
3498 for (IndexToIndex::iterator it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++)
3499 DoLog(0) && (Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl);
3500
3501 return DegeneratedTriangles;
3502}
3503
3504/**
3505 * Purges degenerated triangles from the tesselation structure if they are not
3506 * necessary to keep a single point within the structure.
3507 */
3508void Tesselation::RemoveDegeneratedTriangles()
3509{
3510 Info FunctionInfo(__func__);
3511 IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles();
3512 TriangleMap::iterator finder;
3513 BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL;
3514 int count = 0;
3515
3516 // iterate over all degenerated triangles
3517 for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); !DegeneratedTriangles->empty(); TriangleKeyRunner = DegeneratedTriangles->begin()) {
3518 DoLog(0) && (Log() << Verbose(0) << "Checking presence of triangles " << TriangleKeyRunner->first << " and " << TriangleKeyRunner->second << "." << endl);
3519 // both ways are stored in the map, only use one
3520 if (TriangleKeyRunner->first > TriangleKeyRunner->second)
3521 continue;
3522
3523 // determine from the keys in the map the two _present_ triangles
3524 finder = TrianglesOnBoundary.find(TriangleKeyRunner->first);
3525 if (finder != TrianglesOnBoundary.end())
3526 triangle = finder->second;
3527 else
3528 continue;
3529 finder = TrianglesOnBoundary.find(TriangleKeyRunner->second);
3530 if (finder != TrianglesOnBoundary.end())
3531 partnerTriangle = finder->second;
3532 else
3533 continue;
3534
3535 // determine which lines are shared by the two triangles
3536 bool trianglesShareLine = false;
3537 for (int i = 0; i < 3; ++i)
3538 for (int j = 0; j < 3; ++j)
3539 trianglesShareLine = trianglesShareLine || triangle->lines[i] == partnerTriangle->lines[j];
3540
3541 if (trianglesShareLine && (triangle->endpoints[1]->LinesCount > 2) && (triangle->endpoints[2]->LinesCount > 2) && (triangle->endpoints[0]->LinesCount > 2)) {
3542 // check whether we have to fix lines
3543 BoundaryTriangleSet *Othertriangle = NULL;
3544 BoundaryTriangleSet *OtherpartnerTriangle = NULL;
3545 TriangleMap::iterator TriangleRunner;
3546 for (int i = 0; i < 3; ++i)
3547 for (int j = 0; j < 3; ++j)
3548 if (triangle->lines[i] != partnerTriangle->lines[j]) {
3549 // get the other two triangles
3550 for (TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); ++TriangleRunner)
3551 if (TriangleRunner->second != triangle) {
3552 Othertriangle = TriangleRunner->second;
3553 }
3554 for (TriangleRunner = partnerTriangle->lines[i]->triangles.begin(); TriangleRunner != partnerTriangle->lines[i]->triangles.end(); ++TriangleRunner)
3555 if (TriangleRunner->second != partnerTriangle) {
3556 OtherpartnerTriangle = TriangleRunner->second;
3557 }
3558 /// interchanges their lines so that triangle->lines[i] == partnerTriangle->lines[j]
3559 // the line of triangle receives the degenerated ones
3560 triangle->lines[i]->triangles.erase(Othertriangle->Nr);
3561 triangle->lines[i]->triangles.insert(TrianglePair(partnerTriangle->Nr, partnerTriangle));
3562 for (int k = 0; k < 3; k++)
3563 if (triangle->lines[i] == Othertriangle->lines[k]) {
3564 Othertriangle->lines[k] = partnerTriangle->lines[j];
3565 break;
3566 }
3567 // the line of partnerTriangle receives the non-degenerated ones
3568 partnerTriangle->lines[j]->triangles.erase(partnerTriangle->Nr);
3569 partnerTriangle->lines[j]->triangles.insert(TrianglePair(Othertriangle->Nr, Othertriangle));
3570 partnerTriangle->lines[j] = triangle->lines[i];
3571 }
3572
3573 // erase the pair
3574 count += (int) DegeneratedTriangles->erase(triangle->Nr);
3575 DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl);
3576 RemoveTesselationTriangle(triangle);
3577 count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr);
3578 DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl);
3579 RemoveTesselationTriangle(partnerTriangle);
3580 } else {
3581 DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle << " and its partner " << *partnerTriangle << " because it is essential for at" << " least one of the endpoints to be kept in the tesselation structure." << endl);
3582 }
3583 }
3584 delete (DegeneratedTriangles);
3585 if (count > 0)
3586 LastTriangle = NULL;
3587
3588 DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl);
3589}
3590
3591/** Adds an outside Tesselpoint to the envelope via (two) degenerated triangles.
3592 * We look for the closest point on the boundary, we look through its connected boundary lines and
3593 * seek the one with the minimum angle between its center point and the new point and this base line.
3594 * We open up the line by adding a degenerated triangle, whose other side closes the base line again.
3595 * \param *out output stream for debugging
3596 * \param *point point to add
3597 * \param *LC Linked Cell structure to find nearest point
3598 */
3599void Tesselation::AddBoundaryPointByDegeneratedTriangle(class TesselPoint *point, LinkedCell *LC)
3600{
3601 Info FunctionInfo(__func__);
3602 // find nearest boundary point
3603 class TesselPoint *BackupPoint = NULL;
3604 class TesselPoint *NearestPoint = FindClosestTesselPoint(point->getPosition(), BackupPoint, LC);
3605 class BoundaryPointSet *NearestBoundaryPoint = NULL;
3606 PointMap::iterator PointRunner;
3607
3608 if (NearestPoint == point)
3609 NearestPoint = BackupPoint;
3610 PointRunner = PointsOnBoundary.find(NearestPoint->nr);
3611 if (PointRunner != PointsOnBoundary.end()) {
3612 NearestBoundaryPoint = PointRunner->second;
3613 } else {
3614 DoeLog(1) && (eLog() << Verbose(1) << "I cannot find the boundary point." << endl);
3615 return;
3616 }
3617 DoLog(0) && (Log() << Verbose(0) << "Nearest point on boundary is " << NearestPoint->getName() << "." << endl);
3618
3619 // go through its lines and find the best one to split
3620 Vector CenterToPoint;
3621 Vector BaseLine;
3622 double angle, BestAngle = 0.;
3623 class BoundaryLineSet *BestLine = NULL;
3624 for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) {
3625 BaseLine = (Runner->second->endpoints[0]->node->getPosition()) -
3626 (Runner->second->endpoints[1]->node->getPosition());
3627 CenterToPoint = 0.5 * ((Runner->second->endpoints[0]->node->getPosition()) +
3628 (Runner->second->endpoints[1]->node->getPosition()));
3629 CenterToPoint -= (point->getPosition());
3630 angle = CenterToPoint.Angle(BaseLine);
3631 if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) {
3632 BestAngle = angle;
3633 BestLine = Runner->second;
3634 }
3635 }
3636
3637 // remove one triangle from the chosen line
3638 class BoundaryTriangleSet *TempTriangle = (BestLine->triangles.begin())->second;
3639 BestLine->triangles.erase(TempTriangle->Nr);
3640 int nr = -1;
3641 for (int i = 0; i < 3; i++) {
3642 if (TempTriangle->lines[i] == BestLine) {
3643 nr = i;
3644 break;
3645 }
3646 }
3647
3648 // create new triangle to connect point (connects automatically with the missing spot of the chosen line)
3649 DoLog(2) && (Log() << Verbose(2) << "Adding new triangle points." << endl);
3650 AddTesselationPoint((BestLine->endpoints[0]->node), 0);
3651 AddTesselationPoint((BestLine->endpoints[1]->node), 1);
3652 AddTesselationPoint(point, 2);
3653 DoLog(2) && (Log() << Verbose(2) << "Adding new triangle lines." << endl);
3654 AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
3655 AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
3656 AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
3657 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
3658 BTS->GetNormalVector(TempTriangle->NormalVector);
3659 BTS->NormalVector.Scale(-1.);
3660 DoLog(1) && (Log() << Verbose(1) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl);
3661 AddTesselationTriangle();
3662
3663 // create other side of this triangle and close both new sides of the first created triangle
3664 DoLog(2) && (Log() << Verbose(2) << "Adding new triangle points." << endl);
3665 AddTesselationPoint((BestLine->endpoints[0]->node), 0);
3666 AddTesselationPoint((BestLine->endpoints[1]->node), 1);
3667 AddTesselationPoint(point, 2);
3668 DoLog(2) && (Log() << Verbose(2) << "Adding new triangle lines." << endl);
3669 AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
3670 AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
3671 AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
3672 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
3673 BTS->GetNormalVector(TempTriangle->NormalVector);
3674 DoLog(1) && (Log() << Verbose(1) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl);
3675 AddTesselationTriangle();
3676
3677 // add removed triangle to the last open line of the second triangle
3678 for (int i = 0; i < 3; i++) { // look for the same line as BestLine (only it's its degenerated companion)
3679 if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) {
3680 if (BestLine == BTS->lines[i]) {
3681 DoeLog(0) && (eLog() << Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl);
3682 performCriticalExit();
3683 }
3684 BTS->lines[i]->triangles.insert(pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle));
3685 TempTriangle->lines[nr] = BTS->lines[i];
3686 break;
3687 }
3688 }
3689}
3690;
3691
3692/** Writes the envelope to file.
3693 * \param *out otuput stream for debugging
3694 * \param *filename basename of output file
3695 * \param *cloud PointCloud structure with all nodes
3696 */
3697void Tesselation::Output(const char *filename, const PointCloud * const cloud)
3698{
3699 Info FunctionInfo(__func__);
3700 ofstream *tempstream = NULL;
3701 string NameofTempFile;
3702 string NumberName;
3703
3704 if (LastTriangle != NULL) {
3705 stringstream sstr;
3706 sstr << "-"<< TrianglesOnBoundary.size() << "-" << LastTriangle->getEndpointName(0) << "_" << LastTriangle->getEndpointName(1) << "_" << LastTriangle->getEndpointName(2);
3707 NumberName = sstr.str();
3708 if (DoTecplotOutput) {
3709 string NameofTempFile(filename);
3710 NameofTempFile.append(NumberName);
3711 for (size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
3712 NameofTempFile.erase(npos, 1);
3713 NameofTempFile.append(TecplotSuffix);
3714 DoLog(0) && (Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n");
3715 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
3716 WriteTecplotFile(tempstream, this, cloud, TriangleFilesWritten);
3717 tempstream->close();
3718 tempstream->flush();
3719 delete (tempstream);
3720 }
3721
3722 if (DoRaster3DOutput) {
3723 string NameofTempFile(filename);
3724 NameofTempFile.append(NumberName);
3725 for (size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))
3726 NameofTempFile.erase(npos, 1);
3727 NameofTempFile.append(Raster3DSuffix);
3728 DoLog(0) && (Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n");
3729 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
3730 WriteRaster3dFile(tempstream, this, cloud);
3731 IncludeSphereinRaster3D(tempstream, this, cloud);
3732 tempstream->close();
3733 tempstream->flush();
3734 delete (tempstream);
3735 }
3736 }
3737 if (DoTecplotOutput || DoRaster3DOutput)
3738 TriangleFilesWritten++;
3739}
3740;
3741
3742struct BoundaryPolygonSetCompare
3743{
3744 bool operator()(const BoundaryPolygonSet * s1, const BoundaryPolygonSet * s2) const
3745 {
3746 if (s1->endpoints.size() < s2->endpoints.size())
3747 return true;
3748 else if (s1->endpoints.size() > s2->endpoints.size())
3749 return false;
3750 else { // equality of number of endpoints
3751 PointSet::const_iterator Walker1 = s1->endpoints.begin();
3752 PointSet::const_iterator Walker2 = s2->endpoints.begin();
3753 while ((Walker1 != s1->endpoints.end()) || (Walker2 != s2->endpoints.end())) {
3754 if ((*Walker1)->Nr < (*Walker2)->Nr)
3755 return true;
3756 else if ((*Walker1)->Nr > (*Walker2)->Nr)
3757 return false;
3758 Walker1++;
3759 Walker2++;
3760 }
3761 return false;
3762 }
3763 }
3764};
3765
3766#define UniquePolygonSet set < BoundaryPolygonSet *, BoundaryPolygonSetCompare>
3767
3768/** Finds all degenerated polygons and calls ReTesselateDegeneratedPolygon()/
3769 * \return number of polygons found
3770 */
3771int Tesselation::CorrectAllDegeneratedPolygons()
3772{
3773 Info FunctionInfo(__func__);
3774 /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector
3775 IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles();
3776 set<BoundaryPointSet *> EndpointCandidateList;
3777 pair<set<BoundaryPointSet *>::iterator, bool> InsertionTester;
3778 pair<map<int, Vector *>::iterator, bool> TriangleInsertionTester;
3779 for (PointMap::const_iterator Runner = PointsOnBoundary.begin(); Runner != PointsOnBoundary.end(); Runner++) {
3780 DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Runner->second << "." << endl);
3781 map<int, Vector *> TriangleVectors;
3782 // gather all NormalVectors
3783 DoLog(1) && (Log() << Verbose(1) << "Gathering triangles ..." << endl);
3784 for (LineMap::const_iterator LineRunner = (Runner->second)->lines.begin(); LineRunner != (Runner->second)->lines.end(); LineRunner++)
3785 for (TriangleMap::const_iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {
3786 if (DegeneratedTriangles->find(TriangleRunner->second->Nr) == DegeneratedTriangles->end()) {
3787 TriangleInsertionTester = TriangleVectors.insert(pair<int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)));
3788 if (TriangleInsertionTester.second)
3789 DoLog(1) && (Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl);
3790 } else {
3791 DoLog(1) && (Log() << Verbose(1) << " NOT adding triangle " << *(TriangleRunner->second) << " as it's a simply degenerated one." << endl);
3792 }
3793 }
3794 // check whether there are two that are parallel
3795 DoLog(1) && (Log() << Verbose(1) << "Finding two parallel triangles ..." << endl);
3796 for (map<int, Vector *>::iterator VectorWalker = TriangleVectors.begin(); VectorWalker != TriangleVectors.end(); VectorWalker++)
3797 for (map<int, Vector *>::iterator VectorRunner = VectorWalker; VectorRunner != TriangleVectors.end(); VectorRunner++)
3798 if (VectorWalker != VectorRunner) { // skip equals
3799 const double SCP = VectorWalker->second->ScalarProduct(*VectorRunner->second); // ScalarProduct should result in -1. for degenerated triangles
3800 DoLog(1) && (Log() << Verbose(1) << "Checking " << *VectorWalker->second << " against " << *VectorRunner->second << ": " << SCP << endl);
3801 if (fabs(SCP + 1.) < ParallelEpsilon) {
3802 InsertionTester = EndpointCandidateList.insert((Runner->second));
3803 if (InsertionTester.second)
3804 DoLog(0) && (Log() << Verbose(0) << " Adding " << *Runner->second << " to endpoint candidate list." << endl);
3805 // and break out of both loops
3806 VectorWalker = TriangleVectors.end();
3807 VectorRunner = TriangleVectors.end();
3808 break;
3809 }
3810 }
3811 }
3812 delete DegeneratedTriangles;
3813
3814 /// 3. Find connected endpoint candidates and put them into a polygon
3815 UniquePolygonSet ListofDegeneratedPolygons;
3816 BoundaryPointSet *Walker = NULL;
3817 BoundaryPointSet *OtherWalker = NULL;
3818 BoundaryPolygonSet *Current = NULL;
3819 stack<BoundaryPointSet*> ToCheckConnecteds;
3820 while (!EndpointCandidateList.empty()) {
3821 Walker = *(EndpointCandidateList.begin());
3822 if (Current == NULL) { // create a new polygon with current candidate
3823 DoLog(0) && (Log() << Verbose(0) << "Starting new polygon set at point " << *Walker << endl);
3824 Current = new BoundaryPolygonSet;
3825 Current->endpoints.insert(Walker);
3826 EndpointCandidateList.erase(Walker);
3827 ToCheckConnecteds.push(Walker);
3828 }
3829
3830 // go through to-check stack
3831 while (!ToCheckConnecteds.empty()) {
3832 Walker = ToCheckConnecteds.top(); // fetch ...
3833 ToCheckConnecteds.pop(); // ... and remove
3834 for (LineMap::const_iterator LineWalker = Walker->lines.begin(); LineWalker != Walker->lines.end(); LineWalker++) {
3835 OtherWalker = (LineWalker->second)->GetOtherEndpoint(Walker);
3836 DoLog(1) && (Log() << Verbose(1) << "Checking " << *OtherWalker << endl);
3837 set<BoundaryPointSet *>::iterator Finder = EndpointCandidateList.find(OtherWalker);
3838 if (Finder != EndpointCandidateList.end()) { // found a connected partner
3839 DoLog(1) && (Log() << Verbose(1) << " Adding to polygon." << endl);
3840 Current->endpoints.insert(OtherWalker);
3841 EndpointCandidateList.erase(Finder); // remove from candidates
3842 ToCheckConnecteds.push(OtherWalker); // but check its partners too
3843 } else {
3844 DoLog(1) && (Log() << Verbose(1) << " is not connected to " << *Walker << endl);
3845 }
3846 }
3847 }
3848
3849 DoLog(0) && (Log() << Verbose(0) << "Final polygon is " << *Current << endl);
3850 ListofDegeneratedPolygons.insert(Current);
3851 Current = NULL;
3852 }
3853
3854 const int counter = ListofDegeneratedPolygons.size();
3855
3856 DoLog(0) && (Log() << Verbose(0) << "The following " << counter << " degenerated polygons have been found: " << endl);
3857 for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++)
3858 DoLog(0) && (Log() << Verbose(0) << " " << **PolygonRunner << endl);
3859
3860 /// 4. Go through all these degenerated polygons
3861 for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++) {
3862 stack<int> TriangleNrs;
3863 Vector NormalVector;
3864 /// 4a. Gather all triangles of this polygon
3865 TriangleSet *T = (*PolygonRunner)->GetAllContainedTrianglesFromEndpoints();
3866
3867 // check whether number is bigger than 2, otherwise it's just a simply degenerated one and nothing to do.
3868 if (T->size() == 2) {
3869 DoLog(1) && (Log() << Verbose(1) << " Skipping degenerated polygon, is just a (already simply degenerated) triangle." << endl);
3870 delete (T);
3871 continue;
3872 }
3873
3874 // check whether number is even
3875 // If this case occurs, we have to think about it!
3876 // The Problem is probably due to two degenerated polygons being connected by a bridging, non-degenerated polygon, as somehow one node has
3877 // connections to either polygon ...
3878 if (T->size() % 2 != 0) {
3879 DoeLog(0) && (eLog() << Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl);
3880 performCriticalExit();
3881 }
3882 TriangleSet::iterator TriangleWalker = T->begin(); // is the inner iterator
3883 /// 4a. Get NormalVector for one side (this is "front")
3884 NormalVector = (*TriangleWalker)->NormalVector;
3885 DoLog(1) && (Log() << Verbose(1) << "\"front\" defining triangle is " << **TriangleWalker << " and Normal vector of \"front\" side is " << NormalVector << endl);
3886 TriangleWalker++;
3887 TriangleSet::iterator TriangleSprinter = TriangleWalker; // is the inner advanced iterator
3888 /// 4b. Remove all triangles whose NormalVector is in opposite direction (i.e. "back")
3889 BoundaryTriangleSet *triangle = NULL;
3890 while (TriangleSprinter != T->end()) {
3891 TriangleWalker = TriangleSprinter;
3892 triangle = *TriangleWalker;
3893 TriangleSprinter++;
3894 DoLog(1) && (Log() << Verbose(1) << "Current triangle to test for removal: " << *triangle << endl);
3895 if (triangle->NormalVector.ScalarProduct(NormalVector) < 0) { // if from other side, then delete and remove from list
3896 DoLog(1) && (Log() << Verbose(1) << " Removing ... " << endl);
3897 TriangleNrs.push(triangle->Nr);
3898 T->erase(TriangleWalker);
3899 RemoveTesselationTriangle(triangle);
3900 } else
3901 DoLog(1) && (Log() << Verbose(1) << " Keeping ... " << endl);
3902 }
3903 /// 4c. Copy all "front" triangles but with inverse NormalVector
3904 TriangleWalker = T->begin();
3905 while (TriangleWalker != T->end()) { // go through all front triangles
3906 DoLog(1) && (Log() << Verbose(1) << " Re-creating triangle " << **TriangleWalker << " with NormalVector " << (*TriangleWalker)->NormalVector << endl);
3907 for (int i = 0; i < 3; i++)
3908 AddTesselationPoint((*TriangleWalker)->endpoints[i]->node, i);
3909 AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
3910 AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
3911 AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
3912 if (TriangleNrs.empty())
3913 DoeLog(0) && (eLog() << Verbose(0) << "No more free triangle numbers!" << endl);
3914 BTS = new BoundaryTriangleSet(BLS, TriangleNrs.top()); // copy triangle ...
3915 AddTesselationTriangle(); // ... and add
3916 TriangleNrs.pop();
3917 BTS->NormalVector = -1 * (*TriangleWalker)->NormalVector;
3918 TriangleWalker++;
3919 }
3920 if (!TriangleNrs.empty()) {
3921 DoeLog(0) && (eLog() << Verbose(0) << "There have been less triangles created than removed!" << endl);
3922 }
3923 delete (T); // remove the triangleset
3924 }
3925 IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();
3926 DoLog(0) && (Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl);
3927 IndexToIndex::iterator it;
3928 for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)
3929 DoLog(0) && (Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl);
3930 delete (SimplyDegeneratedTriangles);
3931 /// 5. exit
3932 UniquePolygonSet::iterator PolygonRunner;
3933 while (!ListofDegeneratedPolygons.empty()) {
3934 PolygonRunner = ListofDegeneratedPolygons.begin();
3935 delete (*PolygonRunner);
3936 ListofDegeneratedPolygons.erase(PolygonRunner);
3937 }
3938
3939 return counter;
3940}
3941;
Note: See TracBrowser for help on using the repository browser.