FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2005 Imetric 3D GmbH                                    *
3
 *                                                                         *
4
 *   This file is part of the FreeCAD CAx development system.              *
5
 *                                                                         *
6
 *   This library is free software; you can redistribute it and/or         *
7
 *   modify it under the terms of the GNU Library General Public           *
8
 *   License as published by the Free Software Foundation; either          *
9
 *   version 2 of the License, or (at your option) any later version.      *
10
 *                                                                         *
11
 *   This library  is distributed in the hope that it will be useful,      *
12
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
14
 *   GNU Library General Public License for more details.                  *
15
 *                                                                         *
16
 *   You should have received a copy of the GNU Library General Public     *
17
 *   License along with this library; see the file COPYING.LIB. If not,    *
18
 *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
19
 *   Suite 330, Boston, MA  02111-1307, USA                                *
20
 *                                                                         *
21
 ***************************************************************************/
22

23
#include "PreCompiled.h"
24

25
#ifndef _PreComp_
26
#include <algorithm>
27
#include <vector>
28
#endif
29

30
#include <Base/Matrix.h>
31
#include <Base/Sequencer.h>
32

33
#include "Algorithm.h"
34
#include "Approximation.h"
35
#include "Evaluation.h"
36
#include "Functional.h"
37
#include "Grid.h"
38
#include "Iterator.h"
39
#include "TopoAlgorithm.h"
40

41

42
using namespace MeshCore;
43

44

45
MeshOrientationVisitor::MeshOrientationVisitor() = default;
46

47
bool MeshOrientationVisitor::Visit(const MeshFacet& rclFacet,
48
                                   const MeshFacet& rclFrom,
49
                                   FacetIndex ulFInd,
50
                                   unsigned long ulLevel)
51
{
52
    (void)ulFInd;
53
    (void)ulLevel;
54
    if (!rclFrom.HasSameOrientation(rclFacet)) {
55
        _nonuniformOrientation = true;
56
        return false;
57
    }
58

59
    return true;
60
}
61

62
bool MeshOrientationVisitor::HasNonUnifomOrientedFacets() const
63
{
64
    return _nonuniformOrientation;
65
}
66

67
MeshOrientationCollector::MeshOrientationCollector(std::vector<FacetIndex>& aulIndices,
68
                                                   std::vector<FacetIndex>& aulComplement)
69
    : _aulIndices(aulIndices)
70
    , _aulComplement(aulComplement)
71
{}
72

73
bool MeshOrientationCollector::Visit(const MeshFacet& rclFacet,
74
                                     const MeshFacet& rclFrom,
75
                                     FacetIndex ulFInd,
76
                                     unsigned long ulLevel)
77
{
78
    (void)ulLevel;
79
    // different orientation of rclFacet and rclFrom
80
    if (!rclFacet.HasSameOrientation(rclFrom)) {
81
        // is not marked as false oriented
82
        if (!rclFrom.IsFlag(MeshFacet::TMP0)) {
83
            // mark this facet as false oriented
84
            rclFacet.SetFlag(MeshFacet::TMP0);
85
            _aulIndices.push_back(ulFInd);
86
        }
87
        else {
88
            _aulComplement.push_back(ulFInd);
89
        }
90
    }
91
    else {
92
        // same orientation but if the neighbour rclFrom is false oriented
93
        // then rclFrom is also false oriented
94
        if (rclFrom.IsFlag(MeshFacet::TMP0)) {
95
            // mark this facet as false oriented
96
            rclFacet.SetFlag(MeshFacet::TMP0);
97
            _aulIndices.push_back(ulFInd);
98
        }
99
        else {
100
            _aulComplement.push_back(ulFInd);
101
        }
102
    }
103

104
    return true;
105
}
106

107
MeshSameOrientationCollector::MeshSameOrientationCollector(std::vector<FacetIndex>& aulIndices)
108
    : _aulIndices(aulIndices)
109
{}
110

111
bool MeshSameOrientationCollector::Visit(const MeshFacet& rclFacet,
112
                                         const MeshFacet& rclFrom,
113
                                         FacetIndex ulFInd,
114
                                         unsigned long ulLevel)
115
{
116
    // different orientation of rclFacet and rclFrom
117
    (void)ulLevel;
118
    if (rclFacet.HasSameOrientation(rclFrom)) {
119
        _aulIndices.push_back(ulFInd);
120
    }
121

122
    return true;
123
}
124

125
// ----------------------------------------------------
126

127
MeshEvalOrientation::MeshEvalOrientation(const MeshKernel& rclM)
128
    : MeshEvaluation(rclM)
129
{}
130

131
bool MeshEvalOrientation::Evaluate()
132
{
133
    const MeshFacetArray& rFAry = _rclMesh.GetFacets();
134
    MeshFacetArray::_TConstIterator iBeg = rFAry.begin();
135
    MeshFacetArray::_TConstIterator iEnd = rFAry.end();
136
    for (MeshFacetArray::_TConstIterator it = iBeg; it != iEnd; ++it) {
137
        for (int i = 0; i < 3; i++) {
138
            if (it->_aulNeighbours[i] != FACET_INDEX_MAX) {
139
                const MeshFacet& rclFacet = iBeg[it->_aulNeighbours[i]];
140
                for (int j = 0; j < 3; j++) {
141
                    if (it->_aulPoints[i] == rclFacet._aulPoints[j]) {
142
                        if ((it->_aulPoints[(i + 1) % 3] == rclFacet._aulPoints[(j + 1) % 3])
143
                            || (it->_aulPoints[(i + 2) % 3] == rclFacet._aulPoints[(j + 2) % 3])) {
144
                            return false;  // adjacent face with wrong orientation
145
                        }
146
                    }
147
                }
148
            }
149
        }
150
    }
151

152
    return true;
153
}
154

155
unsigned long MeshEvalOrientation::HasFalsePositives(const std::vector<FacetIndex>& inds) const
156
{
157
    // All faces with wrong orientation (i.e. adjacent faces with a normal flip and their
158
    // neighbours) build a segment and are marked as TMP0. Now we check all border faces of the
159
    // segments with their correct neighbours if there was really a normal flip. If there is no
160
    // normal flip we have a false positive. False-positives can occur if the mesh structure has
161
    // some defects which let the region-grow algorithm fail to detect the faces with wrong
162
    // orientation.
163
    const MeshFacetArray& rFAry = _rclMesh.GetFacets();
164
    MeshFacetArray::_TConstIterator iBeg = rFAry.begin();
165
    for (FacetIndex it : inds) {
166
        const MeshFacet& f = iBeg[it];
167
        for (FacetIndex nbIndex : f._aulNeighbours) {
168
            if (nbIndex != FACET_INDEX_MAX) {
169
                const MeshFacet& n = iBeg[nbIndex];
170
                if (f.IsFlag(MeshFacet::TMP0) && !n.IsFlag(MeshFacet::TMP0)) {
171
                    for (int j = 0; j < 3; j++) {
172
                        if (f.HasSameOrientation(n)) {
173
                            // adjacent face with same orientation => false positive
174
                            return nbIndex;
175
                        }
176
                    }
177
                }
178
            }
179
        }
180
    }
181

182
    return FACET_INDEX_MAX;
183
}
184

185
std::vector<FacetIndex> MeshEvalOrientation::GetIndices() const
186
{
187
    FacetIndex ulStartFacet {}, ulVisited {};
188

189
    if (_rclMesh.CountFacets() == 0) {
190
        return {};
191
    }
192

193
    // reset VISIT flags
194
    MeshAlgorithm cAlg(_rclMesh);
195
    cAlg.ResetFacetFlag(MeshFacet::VISIT);
196
    cAlg.ResetFacetFlag(MeshFacet::TMP0);
197

198
    const MeshFacetArray& rFAry = _rclMesh.GetFacets();
199
    MeshFacetArray::_TConstIterator iTri = rFAry.begin();
200
    MeshFacetArray::_TConstIterator iBeg = rFAry.begin();
201
    MeshFacetArray::_TConstIterator iEnd = rFAry.end();
202

203
    ulStartFacet = 0;
204

205
    std::vector<FacetIndex> uIndices, uComplement;
206
    MeshOrientationCollector clHarmonizer(uIndices, uComplement);
207

208
    while (ulStartFacet != FACET_INDEX_MAX) {
209
        unsigned long wrongFacets = uIndices.size();
210

211
        uComplement.clear();
212
        uComplement.push_back(ulStartFacet);
213
        ulVisited = _rclMesh.VisitNeighbourFacets(clHarmonizer, ulStartFacet) + 1;
214

215
        // In the currently visited component we have found less than 40% as correct
216
        // oriented and the rest as false oriented. So, we decide that it should be the other
217
        // way round and swap the indices of this component.
218
        if (uComplement.size() < static_cast<unsigned long>(0.4f * static_cast<float>(ulVisited))) {
219
            uIndices.erase(uIndices.begin() + wrongFacets, uIndices.end());
220
            uIndices.insert(uIndices.end(), uComplement.begin(), uComplement.end());
221
        }
222

223
        // if the mesh consists of several topologic independent components
224
        // We can search from position 'iTri' on because all elements _before_ are already visited
225
        // what we know from the previous iteration.
226
        MeshIsNotFlag<MeshFacet> flag;
227
        iTri = std::find_if(iTri, iEnd, [flag](const MeshFacet& f) {
228
            return flag(f, MeshFacet::VISIT);
229
        });
230

231
        if (iTri < iEnd) {
232
            ulStartFacet = iTri - iBeg;
233
        }
234
        else {
235
            ulStartFacet = FACET_INDEX_MAX;
236
        }
237
    }
238

239
    // in some very rare cases where we have some strange artifacts in the mesh structure
240
    // we get false-positives. If we find some we check all 'invalid' faces again
241
    cAlg.ResetFacetFlag(MeshFacet::TMP0);
242
    cAlg.SetFacetsFlag(uIndices, MeshFacet::TMP0);
243
    ulStartFacet = HasFalsePositives(uIndices);
244
    while (ulStartFacet != FACET_INDEX_MAX) {
245
        cAlg.ResetFacetsFlag(uIndices, MeshFacet::VISIT);
246
        std::vector<FacetIndex> falsePos;
247
        MeshSameOrientationCollector coll(falsePos);
248
        _rclMesh.VisitNeighbourFacets(coll, ulStartFacet);
249

250
        std::sort(uIndices.begin(), uIndices.end());
251
        std::sort(falsePos.begin(), falsePos.end());
252

253
        std::vector<FacetIndex> diff;
254
        std::back_insert_iterator<std::vector<FacetIndex>> biit(diff);
255
        std::set_difference(uIndices.begin(),
256
                            uIndices.end(),
257
                            falsePos.begin(),
258
                            falsePos.end(),
259
                            biit);
260
        uIndices = diff;
261

262
        cAlg.ResetFacetFlag(MeshFacet::TMP0);
263
        cAlg.SetFacetsFlag(uIndices, MeshFacet::TMP0);
264
        FacetIndex current = ulStartFacet;
265
        ulStartFacet = HasFalsePositives(uIndices);
266
        if (current == ulStartFacet) {
267
            break;  // avoid an endless loop
268
        }
269
    }
270

271
    return uIndices;
272
}
273

274
MeshFixOrientation::MeshFixOrientation(MeshKernel& rclM)
275
    : MeshValidation(rclM)
276
{}
277

278
bool MeshFixOrientation::Fixup()
279
{
280
    MeshTopoAlgorithm(_rclMesh).HarmonizeNormals();
281
    return MeshEvalOrientation(_rclMesh).Evaluate();
282
}
283

284
// ----------------------------------------------------
285

286
MeshEvalSolid::MeshEvalSolid(const MeshKernel& rclM)
287
    : MeshEvaluation(rclM)
288
{}
289

290
bool MeshEvalSolid::Evaluate()
291
{
292
    std::vector<MeshGeomEdge> edges;
293
    _rclMesh.GetEdges(edges);
294
    for (const auto& it : edges) {
295
        if (it._bBorder) {
296
            return false;
297
        }
298
    }
299

300
    return true;
301
}
302

303
// ----------------------------------------------------
304

305
namespace MeshCore
306
{
307

308
struct Edge_Index
309
{
310
    PointIndex p0, p1;
311
    FacetIndex f;
312
};
313

314
struct Edge_Less
315
{
316
    bool operator()(const Edge_Index& x, const Edge_Index& y) const
317
    {
318
        if (x.p0 < y.p0) {
319
            return true;
320
        }
321
        else if (x.p0 > y.p0) {
322
            return false;
323
        }
324
        else if (x.p1 < y.p1) {
325
            return true;
326
        }
327
        else if (x.p1 > y.p1) {
328
            return false;
329
        }
330
        return false;
331
    }
332
};
333

334
}  // namespace MeshCore
335

336
bool MeshEvalTopology::Evaluate()
337
{
338
    // Using and sorting a vector seems to be faster and more memory-efficient
339
    // than a map.
340
    const MeshFacetArray& rclFAry = _rclMesh.GetFacets();
341
    std::vector<Edge_Index> edges;
342
    edges.reserve(3 * rclFAry.size());
343

344
    // build up an array of edges
345
    MeshFacetArray::_TConstIterator pI;
346
    Base::SequencerLauncher seq("Checking topology...", rclFAry.size());
347
    for (pI = rclFAry.begin(); pI != rclFAry.end(); ++pI) {
348
        for (int i = 0; i < 3; i++) {
349
            Edge_Index item {};
350
            item.p0 = std::min<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
351
            item.p1 = std::max<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
352
            item.f = pI - rclFAry.begin();
353
            edges.push_back(item);
354
        }
355

356
        seq.next();
357
    }
358

359
    // sort the edges
360
    std::sort(edges.begin(), edges.end(), Edge_Less());
361

362
    // search for non-manifold edges
363
    PointIndex p0 = POINT_INDEX_MAX, p1 = POINT_INDEX_MAX;
364
    nonManifoldList.clear();
365
    nonManifoldFacets.clear();
366

367
    int count = 0;
368
    std::vector<FacetIndex> facets;
369
    std::vector<Edge_Index>::iterator pE;
370
    for (pE = edges.begin(); pE != edges.end(); ++pE) {
371
        if (p0 == pE->p0 && p1 == pE->p1) {
372
            count++;
373
            facets.push_back(pE->f);
374
        }
375
        else {
376
            if (count > 2) {
377
                // Edge that is shared by more than 2 facets
378
                nonManifoldList.emplace_back(p0, p1);
379
                nonManifoldFacets.push_back(facets);
380
            }
381

382
            p0 = pE->p0;
383
            p1 = pE->p1;
384
            facets.clear();
385
            facets.push_back(pE->f);
386
            count = 1;
387
        }
388
    }
389

390
    return nonManifoldList.empty();
391
}
392

393
// generate indexed edge list which tangents non-manifolds
394
void MeshEvalTopology::GetFacetManifolds(std::vector<FacetIndex>& raclFacetIndList) const
395
{
396
    raclFacetIndList.clear();
397
    const MeshFacetArray& rclFAry = _rclMesh.GetFacets();
398
    MeshFacetArray::_TConstIterator pI;
399

400
    for (pI = rclFAry.begin(); pI != rclFAry.end(); ++pI) {
401
        for (int i = 0; i < 3; i++) {
402
            PointIndex ulPt0 = std::min<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
403
            PointIndex ulPt1 = std::max<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
404
            std::pair<PointIndex, PointIndex> edge = std::make_pair(ulPt0, ulPt1);
405

406
            if (std::find(nonManifoldList.begin(), nonManifoldList.end(), edge)
407
                != nonManifoldList.end()) {
408
                raclFacetIndList.push_back(pI - rclFAry.begin());
409
            }
410
        }
411
    }
412
}
413

414
unsigned long MeshEvalTopology::CountManifolds() const
415
{
416
    return nonManifoldList.size();
417
}
418

419
bool MeshFixTopology::Fixup()
420
{
421
#if 0
422
    MeshEvalTopology eval(_rclMesh);
423
    if (!eval.Evaluate()) {
424
        eval.GetFacetManifolds(deletedFaces);
425

426
        // remove duplicates
427
        std::sort(deletedFaces.begin(), deletedFaces.end());
428
        deletedFaces.erase(std::unique(deletedFaces.begin(), deletedFaces.end()), deletedFaces.end());
429

430
        _rclMesh.DeleteFacets(deletedFaces);
431
    }
432
#else
433
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
434
    deletedFaces.reserve(3 * nonManifoldList.size());  // allocate some memory
435
    for (const auto& it : nonManifoldList) {
436
        std::vector<FacetIndex> non_mf;
437
        non_mf.reserve(it.size());
438
        for (FacetIndex jt : it) {
439
            // facet is only connected with one edge and there causes a non-manifold
440
            unsigned short numOpenEdges = rFaces[jt].CountOpenEdges();
441
            if (numOpenEdges == 2) {
442
                non_mf.push_back(jt);
443
            }
444
            else if (rFaces[jt].IsDegenerated()) {
445
                non_mf.push_back(jt);
446
            }
447
        }
448

449
        // are we able to repair the non-manifold edge by not removing all facets?
450
        if (it.size() - non_mf.size() == 2) {
451
            deletedFaces.insert(deletedFaces.end(), non_mf.begin(), non_mf.end());
452
        }
453
        else {
454
            deletedFaces.insert(deletedFaces.end(), it.begin(), it.end());
455
        }
456
    }
457

458
    if (!deletedFaces.empty()) {
459
        // remove duplicates
460
        std::sort(deletedFaces.begin(), deletedFaces.end());
461
        deletedFaces.erase(std::unique(deletedFaces.begin(), deletedFaces.end()),
462
                           deletedFaces.end());
463

464
        _rclMesh.DeleteFacets(deletedFaces);
465
        _rclMesh.RebuildNeighbours();
466
    }
467
#endif
468

469
    return true;
470
}
471

472
// ---------------------------------------------------------
473

474
bool MeshEvalPointManifolds::Evaluate()
475
{
476
    this->nonManifoldPoints.clear();
477
    this->facetsOfNonManifoldPoints.clear();
478

479
    MeshCore::MeshRefPointToPoints vv_it(_rclMesh);
480
    MeshCore::MeshRefPointToFacets vf_it(_rclMesh);
481

482
    unsigned long ctPoints = _rclMesh.CountPoints();
483
    for (PointIndex index = 0; index < ctPoints; index++) {
484
        // get the local neighbourhood of the point
485
        const std::set<FacetIndex>& nf = vf_it[index];
486
        const std::set<PointIndex>& np = vv_it[index];
487

488
        std::set<unsigned long>::size_type sp {}, sf {};
489
        sp = np.size();
490
        sf = nf.size();
491
        // for an inner point the number of adjacent points is equal to the number of shared faces
492
        // for a boundary point the number of adjacent points is higher by one than the number of
493
        // shared faces for a non-manifold point the number of adjacent points is higher by more
494
        // than one than the number of shared faces
495
        if (sp > sf + 1) {
496
            nonManifoldPoints.push_back(index);
497
            std::vector<FacetIndex> faces;
498
            faces.insert(faces.end(), nf.begin(), nf.end());
499
            this->facetsOfNonManifoldPoints.push_back(faces);
500
        }
501
    }
502

503
    return this->nonManifoldPoints.empty();
504
}
505

506
void MeshEvalPointManifolds::GetFacetIndices(std::vector<FacetIndex>& facets) const
507
{
508
    std::list<std::vector<FacetIndex>>::const_iterator it;
509
    for (it = facetsOfNonManifoldPoints.begin(); it != facetsOfNonManifoldPoints.end(); ++it) {
510
        facets.insert(facets.end(), it->begin(), it->end());
511
    }
512

513
    if (!facets.empty()) {
514
        // remove duplicates
515
        std::sort(facets.begin(), facets.end());
516
        facets.erase(std::unique(facets.begin(), facets.end()), facets.end());
517
    }
518
}
519

520
// ---------------------------------------------------------
521

522
bool MeshEvalSingleFacet::Evaluate()
523
{
524
    // get all non-manifolds
525
    (void)MeshEvalTopology::Evaluate();
526
    /*
527
      // for each (multiple) single linked facet there should
528
      // exist two valid facets sharing the same edge
529
      // so make facet 1 neighbour of facet 2 and vice versa
530
      const std::vector<MeshFacet>& rclFAry = _rclMesh.GetFacets();
531
      std::vector<MeshFacet>::const_iterator pI;
532

533
      std::vector<std::list<unsigned long> > aclMf = _aclManifoldList;
534
      _aclManifoldList.clear();
535

536
      std::map<std::pair<unsigned long, unsigned long>, std::list<unsigned long> > aclHits;
537
      std::map<std::pair<unsigned long, unsigned long>, std::list<unsigned long> >::iterator pEdge;
538

539
      // search for single links (a non-manifold edge and two open edges)
540
      //
541
      //
542
      // build edge <=> facet map
543
      for (pI = rclFAry.begin(); pI != rclFAry.end(); ++pI)
544
      {
545
        for (int i = 0; i < 3; i++)
546
        {
547
          unsigned long ulPt0 = std::min<unsigned long>(pI->_aulPoints[i], pI->_aulPoints[(i+1)%3]);
548
          unsigned long ulPt1 = std::max<unsigned long>(pI->_aulPoints[i], pI->_aulPoints[(i+1)%3]);
549
          aclHits[std::pair<unsigned long, unsigned long>(ulPt0, ulPt1)].push_front(pI -
550
      rclFAry.begin());
551
        }
552
      }
553

554
      // now search for single links
555
      for (std::vector<std::list<unsigned long> >::const_iterator pMF = aclMf.begin(); pMF !=
556
      aclMf.end(); ++pMF)
557
      {
558
        std::list<unsigned long> aulManifolds;
559
        for (std::list<unsigned long>::const_iterator pF = pMF->begin(); pF != pMF->end(); ++pF)
560
        {
561
          const MeshFacet& rclF = rclFAry[*pF];
562

563
          unsigned long ulCtNeighbours=0;
564
          for (int i = 0; i < 3; i++)
565
          {
566
            unsigned long ulPt0 = std::min<unsigned long>(rclF._aulPoints[i],
567
      rclF._aulPoints[(i+1)%3]); unsigned long ulPt1 = std::max<unsigned long>(rclF._aulPoints[i],
568
      rclF._aulPoints[(i+1)%3]); std::pair<unsigned long, unsigned long> clEdge(ulPt0, ulPt1);
569

570
            // number of facets sharing this edge
571
            ulCtNeighbours += aclHits[clEdge].size();
572
          }
573

574
          // single linked found
575
          if (ulCtNeighbours == pMF->size() + 2)
576
            aulManifolds.push_front(*pF);
577
        }
578

579
        if ( aulManifolds.size() > 0 )
580
          _aclManifoldList.push_back(aulManifolds);
581
      }
582
    */
583
    return (nonManifoldList.empty());
584
}
585

586
bool MeshFixSingleFacet::Fixup()
587
{
588
    std::vector<FacetIndex> aulInvalids;
589
    for (const auto& it : _raclManifoldList) {
590
        for (FacetIndex it2 : it) {
591
            aulInvalids.push_back(it2);
592
        }
593
    }
594

595
    _rclMesh.DeleteFacets(aulInvalids);
596
    return true;
597
}
598

599
// ----------------------------------------------------------------
600

601
bool MeshEvalSelfIntersection::Evaluate()
602
{
603
    // Contains bounding boxes for every facet
604
    std::vector<Base::BoundBox3f> boxes;
605

606
    // Splits the mesh using grid for speeding up the calculation
607
    MeshFacetGrid cMeshFacetGrid(_rclMesh);
608
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
609
    MeshGridIterator clGridIter(cMeshFacetGrid);
610
    unsigned long ulGridX {}, ulGridY {}, ulGridZ {};
611
    cMeshFacetGrid.GetCtGrids(ulGridX, ulGridY, ulGridZ);
612

613
    MeshFacetIterator cMFI(_rclMesh);
614
    for (cMFI.Begin(); cMFI.More(); cMFI.Next()) {
615
        boxes.push_back((*cMFI).GetBoundBox());
616
    }
617

618
    // Calculates the intersections
619
    Base::SequencerLauncher seq("Checking for self-intersections...", ulGridX * ulGridY * ulGridZ);
620
    for (clGridIter.Init(); clGridIter.More(); clGridIter.Next()) {
621
        // Get the facet indices, belonging to the current grid unit
622
        std::vector<FacetIndex> aulGridElements;
623
        clGridIter.GetElements(aulGridElements);
624

625
        seq.next();
626
        if (aulGridElements.empty()) {
627
            continue;
628
        }
629

630
        MeshGeomFacet facet1, facet2;
631
        Base::Vector3f pt1, pt2;
632
        for (std::vector<FacetIndex>::iterator it = aulGridElements.begin();
633
             it != aulGridElements.end();
634
             ++it) {
635
            const Base::BoundBox3f& box1 = boxes[*it];
636
            cMFI.Set(*it);
637
            facet1 = *cMFI;
638
            const MeshFacet& rface1 = rFaces[*it];
639
            for (std::vector<FacetIndex>::iterator jt = it; jt != aulGridElements.end(); ++jt) {
640
                if (jt == it) {  // the identical facet
641
                    continue;
642
                }
643
                // If the facets share a common vertex we do not check for self-intersections
644
                // because they could but usually do not intersect each other and the algorithm
645
                // below would detect false-positives, otherwise
646
                const MeshFacet& rface2 = rFaces[*jt];
647
                if (rface1._aulPoints[0] == rface2._aulPoints[0]
648
                    || rface1._aulPoints[0] == rface2._aulPoints[1]
649
                    || rface1._aulPoints[0] == rface2._aulPoints[2]) {
650
                    continue;  // ignore facets sharing a common vertex
651
                }
652
                if (rface1._aulPoints[1] == rface2._aulPoints[0]
653
                    || rface1._aulPoints[1] == rface2._aulPoints[1]
654
                    || rface1._aulPoints[1] == rface2._aulPoints[2]) {
655
                    continue;  // ignore facets sharing a common vertex
656
                }
657
                if (rface1._aulPoints[2] == rface2._aulPoints[0]
658
                    || rface1._aulPoints[2] == rface2._aulPoints[1]
659
                    || rface1._aulPoints[2] == rface2._aulPoints[2]) {
660
                    continue;  // ignore facets sharing a common vertex
661
                }
662

663
                const Base::BoundBox3f& box2 = boxes[*jt];
664
                if (box1 && box2) {
665
                    cMFI.Set(*jt);
666
                    facet2 = *cMFI;
667
                    int ret = facet1.IntersectWithFacet(facet2, pt1, pt2);
668
                    if (ret == 2) {
669
                        // abort after the first detected self-intersection
670
                        return false;
671
                    }
672
                }
673
            }
674
        }
675
    }
676

677
    return true;
678
}
679

680
void MeshEvalSelfIntersection::GetIntersections(
681
    const std::vector<std::pair<FacetIndex, FacetIndex>>& indices,
682
    std::vector<std::pair<Base::Vector3f, Base::Vector3f>>& intersection) const
683
{
684
    intersection.reserve(indices.size());
685
    MeshFacetIterator cMF1(_rclMesh);
686
    MeshFacetIterator cMF2(_rclMesh);
687

688
    Base::Vector3f pt1, pt2;
689
    std::vector<std::pair<FacetIndex, FacetIndex>>::const_iterator it;
690
    for (it = indices.begin(); it != indices.end(); ++it) {
691
        cMF1.Set(it->first);
692
        cMF2.Set(it->second);
693

694
        Base::BoundBox3f box1 = cMF1->GetBoundBox();
695
        Base::BoundBox3f box2 = cMF2->GetBoundBox();
696
        if (box1 && box2) {
697
            int ret = cMF1->IntersectWithFacet(*cMF2, pt1, pt2);
698
            if (ret == 2) {
699
                intersection.emplace_back(pt1, pt2);
700
            }
701
        }
702
    }
703
}
704

705
void MeshEvalSelfIntersection::GetIntersections(
706
    std::vector<std::pair<FacetIndex, FacetIndex>>& intersection) const
707
{
708
    // Contains bounding boxes for every facet
709
    std::vector<Base::BoundBox3f> boxes;
710
    // intersection.clear();
711

712
    // Splits the mesh using grid for speeding up the calculation
713
    MeshFacetGrid cMeshFacetGrid(_rclMesh);
714
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
715
    MeshGridIterator clGridIter(cMeshFacetGrid);
716
    unsigned long ulGridX {}, ulGridY {}, ulGridZ {};
717
    cMeshFacetGrid.GetCtGrids(ulGridX, ulGridY, ulGridZ);
718

719
    MeshFacetIterator cMFI(_rclMesh);
720
    for (cMFI.Begin(); cMFI.More(); cMFI.Next()) {
721
        boxes.push_back((*cMFI).GetBoundBox());
722
    }
723

724
    // Calculates the intersections
725
    Base::SequencerLauncher seq("Checking for self-intersections...", ulGridX * ulGridY * ulGridZ);
726
    for (clGridIter.Init(); clGridIter.More(); clGridIter.Next()) {
727
        // Get the facet indices, belonging to the current grid unit
728
        std::vector<FacetIndex> aulGridElements;
729
        clGridIter.GetElements(aulGridElements);
730

731
        seq.next(true);
732
        if (aulGridElements.empty()) {
733
            continue;
734
        }
735

736
        MeshGeomFacet facet1, facet2;
737
        Base::Vector3f pt1, pt2;
738
        for (std::vector<FacetIndex>::iterator it = aulGridElements.begin();
739
             it != aulGridElements.end();
740
             ++it) {
741
            const Base::BoundBox3f& box1 = boxes[*it];
742
            cMFI.Set(*it);
743
            facet1 = *cMFI;
744
            const MeshFacet& rface1 = rFaces[*it];
745
            for (std::vector<FacetIndex>::iterator jt = it; jt != aulGridElements.end(); ++jt) {
746
                if (jt == it) {  // the identical facet
747
                    continue;
748
                }
749
                // If the facets share a common vertex we do not check for self-intersections
750
                // because they could but usually do not intersect each other and the algorithm
751
                // below would detect false-positives, otherwise
752
                const MeshFacet& rface2 = rFaces[*jt];
753
                if (rface1._aulPoints[0] == rface2._aulPoints[0]
754
                    || rface1._aulPoints[0] == rface2._aulPoints[1]
755
                    || rface1._aulPoints[0] == rface2._aulPoints[2]) {
756
                    continue;  // ignore facets sharing a common vertex
757
                }
758
                if (rface1._aulPoints[1] == rface2._aulPoints[0]
759
                    || rface1._aulPoints[1] == rface2._aulPoints[1]
760
                    || rface1._aulPoints[1] == rface2._aulPoints[2]) {
761
                    continue;  // ignore facets sharing a common vertex
762
                }
763
                if (rface1._aulPoints[2] == rface2._aulPoints[0]
764
                    || rface1._aulPoints[2] == rface2._aulPoints[1]
765
                    || rface1._aulPoints[2] == rface2._aulPoints[2]) {
766
                    continue;  // ignore facets sharing a common vertex
767
                }
768

769
                const Base::BoundBox3f& box2 = boxes[*jt];
770
                if (box1 && box2) {
771
                    cMFI.Set(*jt);
772
                    facet2 = *cMFI;
773
                    int ret = facet1.IntersectWithFacet(facet2, pt1, pt2);
774
                    if (ret == 2) {
775
                        intersection.emplace_back(*it, *jt);
776
                    }
777
                }
778
            }
779
        }
780
    }
781
}
782

783
std::vector<FacetIndex> MeshFixSelfIntersection::GetFacets() const
784
{
785
    std::vector<FacetIndex> indices;
786
    const MeshFacetArray& rFaces = _rclMesh.GetFacets();
787
    for (const auto& it : selfIntersectons) {
788
        unsigned short numOpenEdges1 = rFaces[it.first].CountOpenEdges();
789
        unsigned short numOpenEdges2 = rFaces[it.second].CountOpenEdges();
790

791
        // often we have only single or border facets that intersect other facets
792
        // in this case remove only these facets and keep the other one
793
        if (numOpenEdges1 == 0 && numOpenEdges2 > 0) {
794
            indices.push_back(it.second);
795
        }
796
        else if (numOpenEdges1 > 0 && numOpenEdges2 == 0) {
797
            indices.push_back(it.first);
798
        }
799
        else {
800
            indices.push_back(it.first);
801
            indices.push_back(it.second);
802
        }
803
    }
804

805
    // remove duplicates
806
    std::sort(indices.begin(), indices.end());
807
    indices.erase(std::unique(indices.begin(), indices.end()), indices.end());
808

809
    return indices;
810
}
811

812
bool MeshFixSelfIntersection::Fixup()
813
{
814
    _rclMesh.DeleteFacets(GetFacets());
815
    return true;
816
}
817

818
// ----------------------------------------------------------------
819

820
bool MeshEvalNeighbourhood::Evaluate()
821
{
822
    // Note: If more than two facets are attached to the edge then we have a
823
    // non-manifold edge here.
824
    // This means that the neighbourhood cannot be valid, for sure. But we just
825
    // want to check whether the neighbourhood is valid for topologic correctly
826
    // edges and thus we ignore this case.
827
    // Non-manifolds are an own category of errors and are handled by the class
828
    // MeshEvalTopology.
829
    //
830
    // Using and sorting a vector seems to be faster and more memory-efficient
831
    // than a map.
832
    const MeshFacetArray& rclFAry = _rclMesh.GetFacets();
833
    std::vector<Edge_Index> edges;
834
    edges.reserve(3 * rclFAry.size());
835

836
    // build up an array of edges
837
    MeshFacetArray::_TConstIterator pI;
838
    Base::SequencerLauncher seq("Checking indices...", rclFAry.size());
839
    for (pI = rclFAry.begin(); pI != rclFAry.end(); ++pI) {
840
        for (int i = 0; i < 3; i++) {
841
            Edge_Index item {};
842
            item.p0 = std::min<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
843
            item.p1 = std::max<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
844
            item.f = pI - rclFAry.begin();
845
            edges.push_back(item);
846
        }
847

848
        seq.next();
849
    }
850

851
    // sort the edges
852
    std::sort(edges.begin(), edges.end(), Edge_Less());
853

854
    PointIndex p0 = POINT_INDEX_MAX, p1 = POINT_INDEX_MAX;
855
    PointIndex f0 = FACET_INDEX_MAX, f1 = FACET_INDEX_MAX;
856
    int count = 0;
857
    std::vector<Edge_Index>::iterator pE;
858
    for (pE = edges.begin(); pE != edges.end(); ++pE) {
859
        if (p0 == pE->p0 && p1 == pE->p1) {
860
            f1 = pE->f;
861
            count++;
862
        }
863
        else {
864
            // we handle only the cases for 1 and 2, for all higher
865
            // values we have a non-manifold that is ignored here
866
            if (count == 2) {
867
                const MeshFacet& rFace0 = rclFAry[f0];
868
                const MeshFacet& rFace1 = rclFAry[f1];
869
                unsigned short side0 = rFace0.Side(p0, p1);
870
                unsigned short side1 = rFace1.Side(p0, p1);
871
                // Check whether rFace0 and rFace1 reference each other as
872
                // neighbours
873
                if (rFace0._aulNeighbours[side0] != f1 || rFace1._aulNeighbours[side1] != f0) {
874
                    return false;
875
                }
876
            }
877
            else if (count == 1) {
878
                const MeshFacet& rFace = rclFAry[f0];
879
                unsigned short side = rFace.Side(p0, p1);
880
                // should be "open edge" but isn't marked as such
881
                if (rFace._aulNeighbours[side] != FACET_INDEX_MAX) {
882
                    return false;
883
                }
884
            }
885

886
            p0 = pE->p0;
887
            p1 = pE->p1;
888
            f0 = pE->f;
889
            count = 1;
890
        }
891
    }
892

893
    return true;
894
}
895

896
std::vector<FacetIndex> MeshEvalNeighbourhood::GetIndices() const
897
{
898
    std::vector<FacetIndex> inds;
899
    const MeshFacetArray& rclFAry = _rclMesh.GetFacets();
900
    std::vector<Edge_Index> edges;
901
    edges.reserve(3 * rclFAry.size());
902

903
    // build up an array of edges
904
    MeshFacetArray::_TConstIterator pI;
905
    Base::SequencerLauncher seq("Checking indices...", rclFAry.size());
906
    for (pI = rclFAry.begin(); pI != rclFAry.end(); ++pI) {
907
        for (int i = 0; i < 3; i++) {
908
            Edge_Index item {};
909
            item.p0 = std::min<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
910
            item.p1 = std::max<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
911
            item.f = pI - rclFAry.begin();
912
            edges.push_back(item);
913
        }
914

915
        seq.next();
916
    }
917

918
    // sort the edges
919
    std::sort(edges.begin(), edges.end(), Edge_Less());
920

921
    PointIndex p0 = POINT_INDEX_MAX, p1 = POINT_INDEX_MAX;
922
    PointIndex f0 = FACET_INDEX_MAX, f1 = FACET_INDEX_MAX;
923
    int count = 0;
924
    std::vector<Edge_Index>::iterator pE;
925
    for (pE = edges.begin(); pE != edges.end(); ++pE) {
926
        if (p0 == pE->p0 && p1 == pE->p1) {
927
            f1 = pE->f;
928
            count++;
929
        }
930
        else {
931
            // we handle only the cases for 1 and 2, for all higher
932
            // values we have a non-manifold that is ignored here
933
            if (count == 2) {
934
                const MeshFacet& rFace0 = rclFAry[f0];
935
                const MeshFacet& rFace1 = rclFAry[f1];
936
                unsigned short side0 = rFace0.Side(p0, p1);
937
                unsigned short side1 = rFace1.Side(p0, p1);
938
                // Check whether rFace0 and rFace1 reference each other as
939
                // neighbours
940
                if (rFace0._aulNeighbours[side0] != f1 || rFace1._aulNeighbours[side1] != f0) {
941
                    inds.push_back(f0);
942
                    inds.push_back(f1);
943
                }
944
            }
945
            else if (count == 1) {
946
                const MeshFacet& rFace = rclFAry[f0];
947
                unsigned short side = rFace.Side(p0, p1);
948
                // should be "open edge" but isn't marked as such
949
                if (rFace._aulNeighbours[side] != FACET_INDEX_MAX) {
950
                    inds.push_back(f0);
951
                }
952
            }
953

954
            p0 = pE->p0;
955
            p1 = pE->p1;
956
            f0 = pE->f;
957
            count = 1;
958
        }
959
    }
960

961
    // remove duplicates
962
    std::sort(inds.begin(), inds.end());
963
    inds.erase(std::unique(inds.begin(), inds.end()), inds.end());
964

965
    return inds;
966
}
967

968
bool MeshFixNeighbourhood::Fixup()
969
{
970
    _rclMesh.RebuildNeighbours();
971
    return true;
972
}
973

974
void MeshKernel::RebuildNeighbours(FacetIndex index)
975
{
976
    std::vector<Edge_Index> edges;
977
    edges.reserve(3 * (this->_aclFacetArray.size() - index));
978

979
    // build up an array of edges
980
    MeshFacetArray::_TConstIterator pI;
981
    MeshFacetArray::_TConstIterator pB = this->_aclFacetArray.begin();
982
    for (pI = pB + index; pI != this->_aclFacetArray.end(); ++pI) {
983
        for (int i = 0; i < 3; i++) {
984
            Edge_Index item {};
985
            item.p0 = std::min<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
986
            item.p1 = std::max<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
987
            item.f = pI - pB;
988
            edges.push_back(item);
989
        }
990
    }
991

992
    // sort the edges
993
    // std::sort(edges.begin(), edges.end(), Edge_Less());
994
    int threads = int(std::thread::hardware_concurrency());
995
    MeshCore::parallel_sort(edges.begin(), edges.end(), Edge_Less(), threads);
996

997
    PointIndex p0 = POINT_INDEX_MAX, p1 = POINT_INDEX_MAX;
998
    PointIndex f0 = FACET_INDEX_MAX, f1 = FACET_INDEX_MAX;
999
    int count = 0;
1000
    std::vector<Edge_Index>::iterator pE;
1001
    for (pE = edges.begin(); pE != edges.end(); ++pE) {
1002
        if (p0 == pE->p0 && p1 == pE->p1) {
1003
            f1 = pE->f;
1004
            count++;
1005
        }
1006
        else {
1007
            // we handle only the cases for 1 and 2, for all higher
1008
            // values we have a non-manifold that is ignored here
1009
            if (count == 2) {
1010
                MeshFacet& rFace0 = this->_aclFacetArray[f0];
1011
                MeshFacet& rFace1 = this->_aclFacetArray[f1];
1012
                unsigned short side0 = rFace0.Side(p0, p1);
1013
                unsigned short side1 = rFace1.Side(p0, p1);
1014
                rFace0._aulNeighbours[side0] = f1;
1015
                rFace1._aulNeighbours[side1] = f0;
1016
            }
1017
            else if (count == 1) {
1018
                MeshFacet& rFace = this->_aclFacetArray[f0];
1019
                unsigned short side = rFace.Side(p0, p1);
1020
                rFace._aulNeighbours[side] = FACET_INDEX_MAX;
1021
            }
1022

1023
            p0 = pE->p0;
1024
            p1 = pE->p1;
1025
            f0 = pE->f;
1026
            count = 1;
1027
        }
1028
    }
1029

1030
    // we handle only the cases for 1 and 2, for all higher
1031
    // values we have a non-manifold that is ignored here
1032
    if (count == 2) {
1033
        MeshFacet& rFace0 = this->_aclFacetArray[f0];
1034
        MeshFacet& rFace1 = this->_aclFacetArray[f1];
1035
        unsigned short side0 = rFace0.Side(p0, p1);
1036
        unsigned short side1 = rFace1.Side(p0, p1);
1037
        rFace0._aulNeighbours[side0] = f1;
1038
        rFace1._aulNeighbours[side1] = f0;
1039
    }
1040
    else if (count == 1) {
1041
        MeshFacet& rFace = this->_aclFacetArray[f0];
1042
        unsigned short side = rFace.Side(p0, p1);
1043
        rFace._aulNeighbours[side] = FACET_INDEX_MAX;
1044
    }
1045
}
1046

1047
void MeshKernel::RebuildNeighbours()
1048
{
1049
    // complete rebuild
1050
    RebuildNeighbours(0);
1051
}
1052

1053
// ----------------------------------------------------------------
1054

1055
MeshEigensystem::MeshEigensystem(const MeshKernel& rclB)
1056
    : MeshEvaluation(rclB)
1057
    , _cU(1.0f, 0.0f, 0.0f)
1058
    , _cV(0.0f, 1.0f, 0.0f)
1059
    , _cW(0.0f, 0.0f, 1.0f)
1060
{
1061
    // use the values of world coordinates as default
1062
    Base::BoundBox3f box = _rclMesh.GetBoundBox();
1063
    _fU = box.LengthX();
1064
    _fV = box.LengthY();
1065
    _fW = box.LengthZ();
1066
}
1067

1068
Base::Matrix4D MeshEigensystem::Transform() const
1069
{
1070
    // x,y,c ... vectors
1071
    // R,Q   ... matrices (R is orthonormal so its transposed(=inverse) is equal to Q)
1072
    //
1073
    // from local (x) to world (y,c) coordinates we have the equation
1074
    // y = R * x  + c
1075
    //     <==>
1076
    // x = Q * y - Q * c
1077
    Base::Matrix4D clTMat;
1078
    // rotation part
1079
    clTMat[0][0] = double(_cU.x);
1080
    clTMat[0][1] = double(_cU.y);
1081
    clTMat[0][2] = double(_cU.z);
1082
    clTMat[0][3] = 0.0;
1083
    clTMat[1][0] = double(_cV.x);
1084
    clTMat[1][1] = double(_cV.y);
1085
    clTMat[1][2] = double(_cV.z);
1086
    clTMat[1][3] = 0.0;
1087
    clTMat[2][0] = double(_cW.x);
1088
    clTMat[2][1] = double(_cW.y);
1089
    clTMat[2][2] = double(_cW.z);
1090
    clTMat[2][3] = 0.0;
1091
    clTMat[3][0] = 0.0;
1092
    clTMat[3][1] = 0.0;
1093
    clTMat[3][2] = 0.0;
1094
    clTMat[3][3] = 1.0;
1095

1096
    Base::Vector3f c(_cC);
1097
    c = clTMat * c;
1098

1099
    // translation part
1100
    clTMat[0][3] = double(-c.x);
1101
    clTMat[1][3] = double(-c.y);
1102
    clTMat[2][3] = double(-c.z);
1103

1104
    return clTMat;
1105
}
1106

1107
bool MeshEigensystem::Evaluate()
1108
{
1109
    CalculateLocalSystem();
1110

1111
    float xmin = 0.0f, xmax = 0.0f, ymin = 0.0f, ymax = 0.0f, zmin = 0.0f, zmax = 0.0f;
1112

1113
    Base::Vector3f clVect, clProj;
1114
    float fH {};
1115

1116
    const MeshPointArray& aclPoints = _rclMesh.GetPoints();
1117
    for (const auto& it : aclPoints) {
1118
        // u-direction
1119
        clVect = it - _cC;
1120
        clProj.ProjectToLine(clVect, _cU);
1121
        clVect = clVect + clProj;
1122
        fH = clVect.Length();
1123

1124
        // point vectors in the same direction ?
1125
        if ((clVect * _cU) < 0.0f) {
1126
            fH = -fH;
1127
        }
1128

1129
        xmax = std::max<float>(xmax, fH);
1130
        xmin = std::min<float>(xmin, fH);
1131

1132
        // v-direction
1133
        clVect = it - _cC;
1134
        clProj.ProjectToLine(clVect, _cV);
1135
        clVect = clVect + clProj;
1136
        fH = clVect.Length();
1137

1138
        // point vectors in the same direction ?
1139
        if ((clVect * _cV) < 0.0f) {
1140
            fH = -fH;
1141
        }
1142

1143
        ymax = std::max<float>(ymax, fH);
1144
        ymin = std::min<float>(ymin, fH);
1145

1146
        // w-direction
1147
        clVect = it - _cC;
1148
        clProj.ProjectToLine(clVect, _cW);
1149
        clVect = clVect + clProj;
1150
        fH = clVect.Length();
1151

1152
        // point vectors in the same direction ?
1153
        if ((clVect * _cW) < 0.0f) {
1154
            fH = -fH;
1155
        }
1156

1157
        zmax = std::max<float>(zmax, fH);
1158
        zmin = std::min<float>(zmin, fH);
1159
    }
1160

1161
    _fU = xmax - xmin;
1162
    _fV = ymax - ymin;
1163
    _fW = zmax - zmin;
1164

1165
    return false;  // to call Fixup() if needed
1166
}
1167

1168
Base::Vector3f MeshEigensystem::GetBoundings() const
1169
{
1170
    return Base::Vector3f(_fU, _fV, _fW);
1171
}
1172

1173
void MeshEigensystem::CalculateLocalSystem()
1174
{
1175
    // at least one facet is needed
1176
    if (_rclMesh.CountFacets() < 1) {
1177
        return;  // cannot continue calculation
1178
    }
1179

1180
    const MeshPointArray& aclPoints = _rclMesh.GetPoints();
1181
    MeshPointArray::_TConstIterator it;
1182

1183
    PlaneFit planeFit;
1184
    for (it = aclPoints.begin(); it != aclPoints.end(); ++it) {
1185
        planeFit.AddPoint(*it);
1186
    }
1187

1188
    planeFit.Fit();
1189
    _cC = planeFit.GetBase();
1190
    _cU = planeFit.GetDirU();
1191
    _cV = planeFit.GetDirV();
1192
    _cW = planeFit.GetNormal();
1193

1194
    // set the sign for the vectors
1195
    float fSumU {0.0F}, fSumV {0.0F}, fSumW {0.0F};
1196
    for (it = aclPoints.begin(); it != aclPoints.end(); ++it) {
1197
        float fU = _cU * (*it - _cC);
1198
        float fV = _cV * (*it - _cC);
1199
        float fW = _cW * (*it - _cC);
1200
        fSumU += (fU > 0 ? fU * fU : -fU * fU);
1201
        fSumV += (fV > 0 ? fV * fV : -fV * fV);
1202
        fSumW += (fW > 0 ? fW * fW : -fW * fW);
1203
    }
1204

1205
    // avoid ambiguities concerning directions
1206
    if (fSumU < 0.0f) {
1207
        _cU *= -1.0f;
1208
    }
1209
    if (fSumV < 0.0f) {
1210
        _cV *= -1.0f;
1211
    }
1212
    if (fSumW < 0.0f) {
1213
        _cW *= -1.0f;
1214
    }
1215

1216
    if ((_cU % _cV) * _cW < 0.0f) {
1217
        _cW = -_cW;  // make a right-handed system
1218
    }
1219
}
1220