source: src/tesselation.cpp@ c38826

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 c38826 was af2c424, checked in by Frederik Heber <heber@…>, 14 years ago

LinkedCell constructor rewritten.

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