~ubuntu-branches/ubuntu/maverick/freecad/maverick

« back to all changes in this revision

Viewing changes to src/Mod/Mesh/App/WildMagic4/Wm4TriangulateEC.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Teemu Ikonen
  • Date: 2009-07-16 18:37:41 UTC
  • Revision ID: james.westby@ubuntu.com-20090716183741-oww9kcxqrk991i1n
Tags: upstream-0.8.2237
ImportĀ upstreamĀ versionĀ 0.8.2237

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
// Wild Magic Source Code
 
2
// David Eberly
 
3
// http://www.geometrictools.com
 
4
// Copyright (c) 1998-2007
 
5
//
 
6
// This library is free software; you can redistribute it and/or modify it
 
7
// under the terms of the GNU Lesser General Public License as published by
 
8
// the Free Software Foundation; either version 2.1 of the License, or (at
 
9
// your option) any later version.  The license is available for reading at
 
10
// either of the locations:
 
11
//     http://www.gnu.org/copyleft/lgpl.html
 
12
//     http://www.geometrictools.com/License/WildMagicLicense.pdf
 
13
// The license applies to versions 0 through 4 of Wild Magic.
 
14
//
 
15
// Version: 4.0.5 (2006/10/23)
 
16
 
 
17
#include "Wm4FoundationPCH.h"
 
18
#include "Wm4TriangulateEC.h"
 
19
#include "Wm4Query2Filtered.h"
 
20
#include "Wm4Query2Int64.h"
 
21
#include "Wm4Query2TInteger.h"
 
22
#include "Wm4Query2TRational.h"
 
23
 
 
24
namespace Wm4
 
25
{
 
26
//----------------------------------------------------------------------------
 
27
template <class Real>
 
28
TriangulateEC<Real>::TriangulateEC (const Positions& rkPositions,
 
29
    Query::Type eQueryType, Real fEpsilon, Indices& rkTriangles)
 
30
{
 
31
    // No extra elements are needed for triangulating a simple polygon.
 
32
    InitializePositions(rkPositions,eQueryType,fEpsilon,0);
 
33
 
 
34
    // Triangulate the unindexed polygon.
 
35
    int iVQuantity = (int)rkPositions.size();
 
36
    const int* aiIndex = 0;
 
37
    InitializeVertices(iVQuantity,aiIndex,rkTriangles);
 
38
    DoEarClipping(iVQuantity,aiIndex,rkTriangles);
 
39
}
 
40
//----------------------------------------------------------------------------
 
41
template <class Real>
 
42
TriangulateEC<Real>::TriangulateEC (const Positions& rkPositions,
 
43
    Query::Type eQueryType, Real fEpsilon, const Indices& rkPolygon,
 
44
    Indices& rkTriangles)
 
45
{
 
46
    // No extra elements are needed for triangulating a simple polygon.
 
47
    InitializePositions(rkPositions,eQueryType,fEpsilon,0);
 
48
 
 
49
    // Triangulate the indexed polygon.
 
50
    int iVQuantity = (int)rkPolygon.size();
 
51
    const int* aiIndex = &rkPolygon[0];
 
52
    InitializeVertices(iVQuantity,aiIndex,rkTriangles);
 
53
    DoEarClipping(iVQuantity,aiIndex,rkTriangles);
 
54
}
 
55
//----------------------------------------------------------------------------
 
56
template <class Real>
 
57
TriangulateEC<Real>::TriangulateEC (const Positions& rkPositions,
 
58
    Query::Type eQueryType, Real fEpsilon, const Indices& rkOuter,
 
59
    const Indices& rkInner, Indices& rkTriangles)
 
60
{
 
61
    // Two extra elements are needed to duplicate the endpoints of the edge
 
62
    // introduced to combine outer and inner polygons.
 
63
    InitializePositions(rkPositions,eQueryType,fEpsilon,2);
 
64
 
 
65
    // Combine the outer polygon and the inner polygon into a simple polygon
 
66
    // by inserting two edges connecting mutually visible vertices, one from
 
67
    // the outer polygon and one from the inner polygon.
 
68
    int iNextElement = (int)rkPositions.size();  // next available element
 
69
    IndexMap kMap;
 
70
    Indices kCombined;
 
71
    CombinePolygons(eQueryType,fEpsilon,iNextElement,rkOuter,rkInner,kMap,
 
72
        kCombined);
 
73
 
 
74
    // The combined polygon is now in the format of a simple polygon, albeit
 
75
    // one with coincident edges.
 
76
    int iVQuantity = (int)kCombined.size();
 
77
    const int* aiIndex = &kCombined[0];
 
78
    InitializeVertices(iVQuantity,aiIndex,rkTriangles);
 
79
    DoEarClipping(iVQuantity,aiIndex,rkTriangles);
 
80
 
 
81
    // Map the duplicate indices back to the original indices.
 
82
    RemapIndices(kMap,rkTriangles);
 
83
}
 
84
//----------------------------------------------------------------------------
 
85
template <class Real>
 
86
TriangulateEC<Real>::TriangulateEC (const Positions& rkPositions,
 
87
    Query::Type eQueryType, Real fEpsilon, const Indices& rkOuter,
 
88
    const IndicesArray& rkInners, Indices& rkTriangles)
 
89
{
 
90
    // Two extra elements per inner polygon are needed to duplicate the
 
91
    // endpoints of the edges introduced to combine outer and inner polygons.
 
92
    int iNumInners = (int)rkInners.size();
 
93
    int iExtraElements = 2*iNumInners;
 
94
    InitializePositions(rkPositions,eQueryType,fEpsilon,iExtraElements);
 
95
 
 
96
    // Combine the outer polygon and the inner polygons into a simple polygon
 
97
    // by inserting two edges per inner polygon connecting mutually visible
 
98
    // vertices.
 
99
    int iNextElement = (int)rkPositions.size();
 
100
    Indices kCombined;
 
101
    IndexMap kMap;
 
102
    ProcessOuterAndInners(eQueryType,fEpsilon,rkOuter,rkInners,iNextElement,
 
103
        kMap,kCombined);
 
104
 
 
105
    // The combined polygon is now in the format of a simple polygon, albeit
 
106
    // with coincident edges.
 
107
    int iVQuantity = (int)kCombined.size();
 
108
    const int* aiIndex = &kCombined[0];
 
109
    InitializeVertices(iVQuantity,aiIndex,rkTriangles);
 
110
    DoEarClipping(iVQuantity,aiIndex,rkTriangles);
 
111
 
 
112
    // Map the duplicate indices back to the original indices.
 
113
    RemapIndices(kMap,rkTriangles);
 
114
}
 
115
//----------------------------------------------------------------------------
 
116
template <class Real>
 
117
TriangulateEC<Real>::TriangulateEC (const Positions& rkPositions,
 
118
    Query::Type eQueryType, Real fEpsilon, const Tree* pkTree,
 
119
    Indices& rkTriangles)
 
120
{
 
121
    // Two extra elements per inner polygon are needed to duplicate the
 
122
    // endpoints of the edges introduced to combine outer and inner polygons.
 
123
    int iExtraElements = GetExtraElements(pkTree);
 
124
    InitializePositions(rkPositions,eQueryType,fEpsilon,iExtraElements);
 
125
 
 
126
    int iNextElement = (int)rkPositions.size();
 
127
    IndexMap kMap;
 
128
 
 
129
    std::queue<const Tree*> kQueue;
 
130
    kQueue.push(pkTree);
 
131
    while (kQueue.size() > 0)
 
132
    {
 
133
        const Tree* pkOuterNode = kQueue.front();
 
134
        kQueue.pop();
 
135
 
 
136
        int iNumChildren = (int)pkOuterNode->Child.size();
 
137
        int iVQuantity;
 
138
        const int* aiIndex;
 
139
 
 
140
        if (iNumChildren == 0)
 
141
        {
 
142
            // The outer polygon is a simple polygon (no nested inner
 
143
            // polygons).  Triangulate the simple polygon.
 
144
            iVQuantity = (int)pkOuterNode->Polygon.size();
 
145
            aiIndex = &pkOuterNode->Polygon[0];
 
146
            InitializeVertices(iVQuantity,aiIndex,rkTriangles);
 
147
            DoEarClipping(iVQuantity,aiIndex,rkTriangles);
 
148
        }
 
149
        else
 
150
        {
 
151
            // Place the next level of outer polygon nodes on the queue for
 
152
            // triangulation.
 
153
            std::vector<std::vector<int>*> kInners(iNumChildren);
 
154
            for (int i = 0; i < iNumChildren; i++)
 
155
            {
 
156
                const Tree* pkInnerNode = pkOuterNode->Child[i];
 
157
                kInners[i] = (std::vector<int>*)&pkInnerNode->Polygon;
 
158
                int iNumGrandChildren = (int)pkInnerNode->Child.size();
 
159
                for (int j = 0; j < iNumGrandChildren; j++)
 
160
                {
 
161
                    kQueue.push(pkInnerNode->Child[j]);
 
162
                }
 
163
            }
 
164
 
 
165
            // Combine the outer polygon and the inner polygons into a
 
166
            // simple polygon by inserting two edges per inner polygon
 
167
            // connecting mutually visible vertices.
 
168
            std::vector<int> kCombined;
 
169
            ProcessOuterAndInners(eQueryType,fEpsilon,pkOuterNode->Polygon,
 
170
                kInners,iNextElement,kMap,kCombined);
 
171
 
 
172
            // The combined polygon is now in the format of a simple polygon,
 
173
            // albeit with coincident edges.
 
174
            iVQuantity = (int)kCombined.size();
 
175
            aiIndex = &kCombined[0];
 
176
            InitializeVertices(iVQuantity,aiIndex,rkTriangles);
 
177
            DoEarClipping(iVQuantity,aiIndex,rkTriangles);
 
178
        }
 
179
    }
 
180
 
 
181
    // Map the duplicate indices back to the original indices.
 
182
    RemapIndices(kMap,rkTriangles);
 
183
}
 
184
//----------------------------------------------------------------------------
 
185
template <class Real>
 
186
TriangulateEC<Real>::~TriangulateEC ()
 
187
{
 
188
    WM4_DELETE m_pkQuery;
 
189
}
 
190
//----------------------------------------------------------------------------
 
191
template <class Real>
 
192
void TriangulateEC<Real>::InitializePositions (const Positions& rkPositions,
 
193
    Query::Type eQueryType, Real fEpsilon, int iExtraElements)
 
194
{
 
195
    int iPQuantity = (int)rkPositions.size();
 
196
    assert(iPQuantity >= 3);
 
197
    int iPEQuantity = iPQuantity + iExtraElements;
 
198
    m_kSPositions.resize(iPEQuantity);
 
199
 
 
200
    if (eQueryType == Query::QT_FILTERED)
 
201
    {
 
202
        assert((Real)0.0 <= fEpsilon && fEpsilon <= (Real)1.0);
 
203
    }
 
204
 
 
205
    Vector2<Real> kMin, kMax, kRange;
 
206
    Real fScale, fRMax;
 
207
    int i;
 
208
 
 
209
    switch (eQueryType)
 
210
    {
 
211
    case Query::QT_INT64:
 
212
        // Transform the vertices to the square [0,2^{20}]^2.
 
213
        Vector2<Real>::ComputeExtremes(iPQuantity,&rkPositions[0],kMin,kMax);
 
214
        kRange = kMax - kMin;
 
215
        fRMax = (kRange[0] >= kRange[1] ? kRange[0] : kRange[1]);
 
216
        fScale = ((Real)(1 << 20))/fRMax;
 
217
        for (i = 0; i < iPQuantity; i++)
 
218
        {
 
219
            m_kSPositions[i] = (rkPositions[i] - kMin)*fScale;
 
220
        }
 
221
        m_pkQuery = WM4_NEW Query2Int64<Real>(iPEQuantity,&m_kSPositions[0]);
 
222
        return;
 
223
 
 
224
    case Query::QT_INTEGER:
 
225
        // Transform the vertices to the square [0,2^{24}]^2.
 
226
        Vector2<Real>::ComputeExtremes(iPQuantity,&rkPositions[0],kMin,kMax);
 
227
        kRange = kMax - kMin;
 
228
        fRMax = (kRange[0] >= kRange[1] ? kRange[0] : kRange[1]);
 
229
        fScale = ((Real)(1 << 24))/fRMax;
 
230
        for (i = 0; i < iPQuantity; i++)
 
231
        {
 
232
            m_kSPositions[i] = (rkPositions[i] - kMin)*fScale;
 
233
        }
 
234
        m_pkQuery = WM4_NEW Query2TInteger<Real>(iPEQuantity,
 
235
            &m_kSPositions[0]);
 
236
        return;
 
237
 
 
238
    case Query::QT_REAL:
 
239
        // Transform the vertices to the square [0,1]^2.
 
240
        Vector2<Real>::ComputeExtremes(iPQuantity,&rkPositions[0],kMin,kMax);
 
241
        kRange = kMax - kMin;
 
242
        fRMax = (kRange[0] >= kRange[1] ? kRange[0] : kRange[1]);
 
243
        fScale = ((Real)1.0)/fRMax;
 
244
        for (i = 0; i < iPQuantity; i++)
 
245
        {
 
246
            m_kSPositions[i] = (rkPositions[i] - kMin)*fScale;
 
247
        }
 
248
        m_pkQuery = WM4_NEW Query2<Real>(iPEQuantity,&m_kSPositions[0]);
 
249
        return;
 
250
 
 
251
    case Query::QT_RATIONAL:
 
252
        // No transformation of the input data.  Make a copy that can be
 
253
        // expanded when triangulating polygons with holes.
 
254
        for (i = 0; i < iPQuantity; i++)
 
255
        {
 
256
            m_kSPositions[i] = rkPositions[i];
 
257
        }
 
258
        m_pkQuery = WM4_NEW Query2TRational<Real>(iPEQuantity,
 
259
            &m_kSPositions[0]);
 
260
        return;
 
261
 
 
262
    case Query::QT_FILTERED:
 
263
        // No transformation of the input data.  Make a copy that can be
 
264
        // expanded when triangulating polygons with holes.
 
265
        for (i = 0; i < iPQuantity; i++)
 
266
        {
 
267
            m_kSPositions[i] = rkPositions[i];
 
268
        }
 
269
        m_pkQuery = WM4_NEW Query2Filtered<Real>(iPEQuantity,
 
270
            &m_kSPositions[0],fEpsilon);
 
271
        return;
 
272
    }
 
273
 
 
274
    assert(false);
 
275
}
 
276
//----------------------------------------------------------------------------
 
277
template <class Real>
 
278
void TriangulateEC<Real>::InitializeVertices (int iVQuantity,
 
279
    const int* aiIndex, std::vector<int>& rkTriangle)
 
280
{
 
281
    m_kVertex.clear();
 
282
    m_kVertex.resize(iVQuantity);
 
283
    m_iCFirst = -1;
 
284
    m_iCLast = -1;
 
285
    m_iRFirst = -1;
 
286
    m_iRLast = -1;
 
287
    m_iEFirst = -1;
 
288
    m_iELast = -1;
 
289
 
 
290
    // Create a circular list of the polygon vertices for dynamic removal of
 
291
    // vertices.
 
292
    int iVQm1 = iVQuantity - 1;
 
293
    int i;
 
294
    for (i = 0; i <= iVQm1; i++)
 
295
    {
 
296
        Vertex& rkV = V(i);
 
297
        rkV.Index = (aiIndex ? aiIndex[i] : i);
 
298
        rkV.VPrev = (i > 0 ? i-1 : iVQm1);
 
299
        rkV.VNext = (i < iVQm1 ? i+1 : 0);
 
300
    }
 
301
 
 
302
    // Create a circular list of the polygon vertices for dynamic removal of
 
303
    // vertices.  Keep track of two linear sublists, one for the convex
 
304
    // vertices and one for the reflex vertices.  This is an O(N) process
 
305
    // where N is the number of polygon vertices.
 
306
    for (i = 0; i <= iVQm1; i++)
 
307
    {
 
308
        if (IsConvex(i))
 
309
        {
 
310
            InsertAfterC(i);
 
311
        }
 
312
        else
 
313
        {
 
314
            InsertAfterR(i);
 
315
        }
 
316
    }
 
317
}
 
318
//----------------------------------------------------------------------------
 
319
template <class Real>
 
320
void TriangulateEC<Real>::DoEarClipping (int iVQuantity, const int* aiIndex,
 
321
    std::vector<int>& rkTriangle)
 
322
{
 
323
    // If the polygon is convex, just create a triangle fan.
 
324
    int i;
 
325
    if (m_iRFirst == -1)
 
326
    {
 
327
        int iVQm1 = iVQuantity - 1;
 
328
        if (aiIndex)
 
329
        {
 
330
            for (i = 1; i < iVQm1; i++)
 
331
            {
 
332
                rkTriangle.push_back(aiIndex[0]);
 
333
                rkTriangle.push_back(aiIndex[i]);
 
334
                rkTriangle.push_back(aiIndex[i+1]);
 
335
            }
 
336
        }
 
337
        else
 
338
        {
 
339
            for (i = 1; i < iVQm1; i++)
 
340
            {
 
341
                rkTriangle.push_back(0);
 
342
                rkTriangle.push_back(i);
 
343
                rkTriangle.push_back(i+1);
 
344
            }
 
345
        }
 
346
        return;
 
347
    }
 
348
 
 
349
    // Identify the ears and build a circular list of them.  Let V0, V1, and
 
350
    // V2 be consecutive vertices forming a triangle T.  The vertex V1 is an
 
351
    // ear if no other vertices of the polygon lie inside T.  Although it is
 
352
    // enough to show that V1 is not an ear by finding at least one other
 
353
    // vertex inside T, it is sufficient to search only the reflex vertices.
 
354
    // This is an O(C*R) process, where C is the number of convex vertices and
 
355
    // R is the number of reflex vertices with N = C+R.  The order is O(N^2),
 
356
    // for example when C = R = N/2.
 
357
    for (i = m_iCFirst; i != -1; i = V(i).SNext)
 
358
    {
 
359
        if (IsEar(i))
 
360
        {
 
361
            InsertEndE(i);
 
362
        }
 
363
    }
 
364
    V(m_iEFirst).EPrev = m_iELast;
 
365
    V(m_iELast).ENext = m_iEFirst;
 
366
 
 
367
    // Remove the ears, one at a time.
 
368
    while (true)
 
369
    {
 
370
        // Add the triangle with the ear to the output list of triangles.
 
371
        int iVPrev = V(m_iEFirst).VPrev;
 
372
        int iVNext = V(m_iEFirst).VNext;
 
373
        rkTriangle.push_back(V(iVPrev).Index);
 
374
        rkTriangle.push_back(V(m_iEFirst).Index);
 
375
        rkTriangle.push_back(V(iVNext).Index);
 
376
 
 
377
        // Remove the vertex corresponding to the ear.
 
378
        RemoveV(m_iEFirst);
 
379
        if (--iVQuantity == 3)
 
380
        {
 
381
            // Only one triangle remains, just remove the ear and copy it.
 
382
            m_iEFirst = RemoveE(m_iEFirst);
 
383
            iVPrev = V(m_iEFirst).VPrev;
 
384
            iVNext = V(m_iEFirst).VNext;
 
385
            rkTriangle.push_back(V(iVPrev).Index);
 
386
            rkTriangle.push_back(V(m_iEFirst).Index);
 
387
            rkTriangle.push_back(V(iVNext).Index);
 
388
            break;
 
389
        }
 
390
 
 
391
        // Removal of the ear can cause an adjacent vertex to become an ear
 
392
        // or to stop being an ear.
 
393
        Vertex& rkVPrev = V(iVPrev);
 
394
        if (rkVPrev.IsEar)
 
395
        {
 
396
            if (!IsEar(iVPrev))
 
397
            {
 
398
                RemoveE(iVPrev);
 
399
            }
 
400
        }
 
401
        else
 
402
        {
 
403
            bool bWasReflex = !rkVPrev.IsConvex;
 
404
            if (IsConvex(iVPrev))
 
405
            {
 
406
                if (bWasReflex)
 
407
                {
 
408
                    RemoveR(iVPrev);
 
409
                }
 
410
 
 
411
                if (IsEar(iVPrev))
 
412
                {
 
413
                    InsertBeforeE(iVPrev);
 
414
                }
 
415
            }
 
416
        }
 
417
 
 
418
        Vertex& rkVNext = V(iVNext);
 
419
        if (rkVNext.IsEar)
 
420
        {
 
421
            if (!IsEar(iVNext))
 
422
            {
 
423
                RemoveE(iVNext);
 
424
            }
 
425
        }
 
426
        else
 
427
        {
 
428
            bool bWasReflex = !rkVNext.IsConvex;
 
429
            if (IsConvex(iVNext))
 
430
            {
 
431
                if (bWasReflex)
 
432
                {
 
433
                    RemoveR(iVNext);
 
434
                }
 
435
 
 
436
                if (IsEar(iVNext))
 
437
                {
 
438
                    InsertAfterE(iVNext);
 
439
                }
 
440
            }
 
441
        }
 
442
 
 
443
        // Remove the ear.
 
444
        m_iEFirst = RemoveE(m_iEFirst);
 
445
    }
 
446
}
 
447
//----------------------------------------------------------------------------
 
448
template <class Real>
 
449
int TriangulateEC<Real>::TriangleQuery (const Vector2<Real>& rkPoint,
 
450
    Query::Type eQueryType, Real fEpsilon, const Vector2<Real> akSTriangle[3])
 
451
    const
 
452
{
 
453
    switch (eQueryType)
 
454
    {
 
455
    case Query::QT_INT64:
 
456
        return Query2Int64<Real>(3,akSTriangle).ToTriangle(rkPoint,0,1,2);
 
457
 
 
458
    case Query::QT_INTEGER:
 
459
        return Query2TInteger<Real>(3,akSTriangle).ToTriangle(rkPoint,0,1,2);
 
460
 
 
461
    case Query::QT_REAL:
 
462
        return Query2<Real>(3,akSTriangle).ToTriangle(rkPoint,0,1,2);
 
463
 
 
464
    case Query::QT_RATIONAL:
 
465
        return Query2TRational<Real>(3,akSTriangle).ToTriangle(rkPoint,0,1,2);
 
466
 
 
467
    case Query::QT_FILTERED:
 
468
        return Query2Filtered<Real>(3,akSTriangle,fEpsilon).ToTriangle(
 
469
            rkPoint,0,1,2);
 
470
    }
 
471
 
 
472
    assert(false);
 
473
    return 1;
 
474
}
 
475
//----------------------------------------------------------------------------
 
476
template <class Real>
 
477
void TriangulateEC<Real>::CombinePolygons (Query::Type eQueryType,
 
478
    Real fEpsilon, int iNextElement, const Indices& rkOuter,
 
479
    const Indices& rkInner, IndexMap& rkMap, Indices& rkCombined)
 
480
{
 
481
    int iOQuantity = (int)rkOuter.size();
 
482
    int iIQuantity = (int)rkInner.size();
 
483
 
 
484
    // Locate the inner-polygon vertex of maximum x-value, call this vertex M.
 
485
    Real fXMax = m_kSPositions[rkInner[0]][0];
 
486
    int iXMaxIndex = 0;
 
487
    int i;
 
488
    for (i = 1; i < iIQuantity; i++)
 
489
    {
 
490
        Real fX = m_kSPositions[rkInner[i]][0];
 
491
        if (fX > fXMax)
 
492
        {
 
493
            fXMax = fX;
 
494
            iXMaxIndex = i;
 
495
        }
 
496
    }
 
497
    Vector2<Real> kM = m_kSPositions[rkInner[iXMaxIndex]];
 
498
 
 
499
    // Find the edge whose intersection Intr with the ray M+t*(1,0) minimizes
 
500
    // the ray parameter t >= 0.
 
501
    Vector2<Real> kIntr(Math<Real>::MAX_REAL,kM[1]);
 
502
    int iV0Min = -1, iV1Min = -1, iEndMin = -1;
 
503
    int i0, i1;
 
504
    for (i0 = iOQuantity-1, i1 = 0; i1 < iOQuantity; i0 = i1++)
 
505
    {
 
506
        // Only consider edges for which the first vertex is below (or on)
 
507
        // the ray and the second vertex is above (or on) the ray.
 
508
        Vector2<Real> kDiff0 = m_kSPositions[rkOuter[i0]] - kM;
 
509
        if (kDiff0[1] > (Real)0.0)
 
510
        {
 
511
            continue;
 
512
        }
 
513
        Vector2<Real> kDiff1 = m_kSPositions[rkOuter[i1]] - kM;
 
514
        if (kDiff1[1] < (Real)0.0)
 
515
        {
 
516
            continue;
 
517
        }
 
518
 
 
519
        // At this time, diff0.y <= 0 and diff1.y >= 0.
 
520
        Real fS, fT;
 
521
        int iCurrentEndMin = -1;
 
522
        if (kDiff0[1] < (Real)0.0)
 
523
        {
 
524
            if (kDiff1[1] > (Real)0.0)
 
525
            {
 
526
                // The intersection of the edge and ray occurs at an interior
 
527
                // edge point.
 
528
                fS = kDiff0[1]/(kDiff0[1] - kDiff1[1]);
 
529
                fT = kDiff0[0] + fS*(kDiff1[0] - kDiff0[0]);
 
530
            }
 
531
            else  // diff1.y == 0
 
532
            {
 
533
                // The vertex Outer[i1] is the intersection of the edge and
 
534
                // the ray.
 
535
                fT = kDiff1[0];
 
536
                iCurrentEndMin = i1;
 
537
            }
 
538
        }
 
539
        else  // diff0.y == 0
 
540
        {
 
541
            if (kDiff1[1] > (Real)0.0)
 
542
            {
 
543
                // The vertex Outer[i0] is the intersection of the edge and
 
544
                // the ray;
 
545
                fT = kDiff0[0];
 
546
                iCurrentEndMin = i0;
 
547
            }
 
548
            else  // diff1.y == 0
 
549
            {
 
550
                if (kDiff0[0] < kDiff1[0])
 
551
                {
 
552
                    fT = kDiff0[0];
 
553
                    iCurrentEndMin = i0;
 
554
                }
 
555
                else
 
556
                {
 
557
                    fT = kDiff1[0];
 
558
                    iCurrentEndMin = i1;
 
559
                }
 
560
            }
 
561
        }
 
562
 
 
563
        if ((Real)0.0 <= fT && fT < kIntr[0])
 
564
        {
 
565
            kIntr[0] = fT;
 
566
            iV0Min = i0;
 
567
            iV1Min = i1;
 
568
            if (iCurrentEndMin == -1)
 
569
            {
 
570
                // The current closest point is an edge-interior point.
 
571
                iEndMin = -1;
 
572
            }
 
573
            else
 
574
            {
 
575
                // The current closest point is a vertex.
 
576
                iEndMin = iCurrentEndMin;
 
577
            }
 
578
        }
 
579
    }
 
580
 
 
581
    int iMaxCosIndex;
 
582
    if (iEndMin == -1)
 
583
    {
 
584
        // Select one of Outer[v0min] and Outer[v1min] that has an x-value
 
585
        // larger than M.x, call this vertex P.  The triangle <M,I,P> must
 
586
        // contain an outer-polygon vertex that is visible to M, which is
 
587
        // possibly P itself.
 
588
        Vector2<Real> akSTriangle[3];  // <P,M,I> or <P,I,M>
 
589
        int iPIndex;
 
590
        if (m_kSPositions[rkOuter[iV0Min]][0] >
 
591
            m_kSPositions[rkOuter[iV1Min]][0])
 
592
        {
 
593
            akSTriangle[0] = m_kSPositions[rkOuter[iV0Min]];
 
594
            akSTriangle[1] = kIntr;
 
595
            akSTriangle[2] = kM;
 
596
            iPIndex = iV0Min;
 
597
        }
 
598
        else
 
599
        {
 
600
            akSTriangle[0] = m_kSPositions[rkOuter[iV1Min]];
 
601
            akSTriangle[1] = kM;
 
602
            akSTriangle[2] = kIntr;
 
603
            iPIndex = iV1Min;
 
604
        }
 
605
 
 
606
        // If any outer-polygon vertices other than P are inside the triangle
 
607
        // <M,I,P>, then at least one of these vertices must be a reflex
 
608
        // vertex.  It is sufficient to locate the reflex vertex R (if any)
 
609
        // in <M,I,P> that minimizes the angle between R-M and (1,0).  The
 
610
        // data member m_pkQuery is used for the reflex query.
 
611
        Vector2<Real> kDiff = akSTriangle[0] - kM;
 
612
        Real fMaxSqrLen = kDiff.SquaredLength();
 
613
        Real fMaxCos = kDiff[0]*kDiff[0]/fMaxSqrLen;
 
614
        iMaxCosIndex = iPIndex;
 
615
        for (i = 0; i < iOQuantity; i++)
 
616
        {
 
617
            if (i == iPIndex)
 
618
            {
 
619
                continue;
 
620
            }
 
621
 
 
622
            int iCurr = rkOuter[i];
 
623
            int iPrev = rkOuter[(i+iOQuantity-1) % iOQuantity];
 
624
            int iNext = rkOuter[(i+1) % iOQuantity];
 
625
            if (m_pkQuery->ToLine(iCurr,iPrev,iNext) <= 0
 
626
            &&  TriangleQuery(m_kSPositions[iCurr],eQueryType,fEpsilon,
 
627
                    akSTriangle) <= 0)
 
628
            {
 
629
                // The vertex is reflex and inside the <M,I,P> triangle.
 
630
                kDiff = m_kSPositions[iCurr] - kM;
 
631
                Real fSqrLen = kDiff.SquaredLength();
 
632
                Real fCos = kDiff[0]*kDiff[0]/fSqrLen;
 
633
                if (fCos > fMaxCos)
 
634
                {
 
635
                    // The reflex vertex forms a smaller angle with the
 
636
                    // positive x-axis, so it becomes the new visible
 
637
                    // candidate.
 
638
                    fMaxSqrLen = fSqrLen;
 
639
                    fMaxCos = fCos;
 
640
                    iMaxCosIndex = i;
 
641
                }
 
642
                else if (fCos == fMaxCos && fSqrLen < fMaxSqrLen)
 
643
                {
 
644
                    // The reflex vertex has angle equal to the current
 
645
                    // minimum but the length is smaller, so it becomes the
 
646
                    // new visible candidate.
 
647
                    fMaxSqrLen = fSqrLen;
 
648
                    iMaxCosIndex = i;
 
649
                }
 
650
            }
 
651
        }
 
652
    }
 
653
    else
 
654
    {
 
655
        iMaxCosIndex = iEndMin;
 
656
    }
 
657
 
 
658
    // The visible vertices are Position[Inner[iXMaxIndex]] and
 
659
    // Position[Outer[iMaxCosIndex]].  Two coincident edges with these
 
660
    // endpoints are inserted to connect the outer and inner polygons into a
 
661
    // simple polygon.  Each of the two Position[] values must be duplicated,
 
662
    // because the original might be convex (or reflex) and the duplicate is
 
663
    // reflex (or convex).  The ear-clipping algorithm needs to distinguish
 
664
    // between them.
 
665
    rkCombined.resize(iOQuantity+iIQuantity+2);
 
666
    int iCIndex = 0;
 
667
    for (i = 0; i <= iMaxCosIndex; i++, iCIndex++)
 
668
    {
 
669
        rkCombined[iCIndex] = rkOuter[i];
 
670
    }
 
671
 
 
672
    for (i = 0; i < iIQuantity; i++, iCIndex++)
 
673
    {
 
674
        int j = (iXMaxIndex + i) % iIQuantity;
 
675
        rkCombined[iCIndex] = rkInner[j];
 
676
    }
 
677
 
 
678
    int iInnerIndex = rkInner[iXMaxIndex];
 
679
    m_kSPositions[iNextElement] = m_kSPositions[iInnerIndex];
 
680
    rkCombined[iCIndex] = iNextElement;
 
681
    IndexMap::iterator pkIter = rkMap.find(iInnerIndex);
 
682
    if (pkIter != rkMap.end())
 
683
    {
 
684
        iInnerIndex = pkIter->second;
 
685
    }
 
686
    rkMap[iNextElement] = iInnerIndex;
 
687
    iCIndex++;
 
688
    iNextElement++;
 
689
 
 
690
    int iOuterIndex = rkOuter[iMaxCosIndex];
 
691
    m_kSPositions[iNextElement] = m_kSPositions[iOuterIndex];
 
692
    rkCombined[iCIndex] = iNextElement;
 
693
    pkIter = rkMap.find(iOuterIndex);
 
694
    if (pkIter != rkMap.end())
 
695
    {
 
696
        iOuterIndex = pkIter->second;
 
697
    }
 
698
    rkMap[iNextElement] = iOuterIndex;
 
699
    iCIndex++;
 
700
    iNextElement++;
 
701
 
 
702
    for (i = iMaxCosIndex+1; i < iOQuantity; i++, iCIndex++)
 
703
    {
 
704
        rkCombined[iCIndex] = rkOuter[i];
 
705
    }
 
706
}
 
707
//----------------------------------------------------------------------------
 
708
template <class Real>
 
709
void TriangulateEC<Real>::ProcessOuterAndInners (Query::Type eQueryType,
 
710
    Real fEpsilon, const Indices& rkOuter, const IndicesArray& rkInners,
 
711
    int& riNextElement, IndexMap& rkMap, Indices& rkCombined)
 
712
{
 
713
    // Sort the inner polygons based on maximum x-values.
 
714
    int iNumInners = (int)rkInners.size();
 
715
    std::vector<std::pair<Real,int> > kPairs(iNumInners);
 
716
    int i;
 
717
    for (i = 0; i < iNumInners; i++)
 
718
    {
 
719
        const Indices& rkInner = *rkInners[i];
 
720
        int iVQuantity = (int)rkInner.size();
 
721
        Real fXMax = m_kSPositions[rkInner[0]][0];
 
722
        for (int j = 1; j < iVQuantity; j++)
 
723
        {
 
724
            Real fX = m_kSPositions[rkInner[j]][0];
 
725
            if (fX > fXMax)
 
726
            {
 
727
                fXMax = fX;
 
728
            }
 
729
        }
 
730
        kPairs[i].first = fXMax;
 
731
        kPairs[i].second = i;
 
732
    }
 
733
    std::sort(kPairs.begin(),kPairs.end());
 
734
 
 
735
    // Merge the inner polygons with the outer polygon.
 
736
    Indices kCurrentOuter = rkOuter;
 
737
    for (i = iNumInners-1; i >= 0; i--)
 
738
    {
 
739
        const Indices& rkInner = *rkInners[kPairs[i].second];
 
740
        Indices kCurrentCombined;
 
741
        CombinePolygons(eQueryType,fEpsilon,riNextElement,kCurrentOuter,
 
742
            rkInner,rkMap,kCurrentCombined);
 
743
        kCurrentOuter = kCurrentCombined;
 
744
        riNextElement += 2;
 
745
    }
 
746
 
 
747
    for (i = 0; i < (int)kCurrentOuter.size(); i++)
 
748
    {
 
749
        rkCombined.push_back(kCurrentOuter[i]);
 
750
    }
 
751
}
 
752
//----------------------------------------------------------------------------
 
753
template <class Real>
 
754
void TriangulateEC<Real>::RemapIndices (const IndexMap& rkMap,
 
755
    Indices& rkTriangles) const
 
756
{
 
757
    // The triangulation includes indices to the duplicated outer and inner
 
758
    // vertices.  These indices must be mapped back to the original ones.
 
759
    for (int i = 0; i < (int)rkTriangles.size(); i++)
 
760
    {
 
761
        IndexMap::const_iterator pkIter = rkMap.find(rkTriangles[i]);
 
762
        if (pkIter != rkMap.end())
 
763
        {
 
764
            rkTriangles[i] = pkIter->second;
 
765
        }
 
766
    }
 
767
}
 
768
//----------------------------------------------------------------------------
 
769
 
 
770
//----------------------------------------------------------------------------
 
771
// Vertex list handling
 
772
//----------------------------------------------------------------------------
 
773
template <class Real>
 
774
typename TriangulateEC<Real>::Vertex& TriangulateEC<Real>::V (int i)
 
775
{
 
776
    return m_kVertex[i];
 
777
}
 
778
//----------------------------------------------------------------------------
 
779
template <class Real>
 
780
bool TriangulateEC<Real>::IsConvex (int i)
 
781
{
 
782
    Vertex& rkV = V(i);
 
783
    int iCurr = rkV.Index;
 
784
    int iPrev = V(rkV.VPrev).Index;
 
785
    int iNext = V(rkV.VNext).Index;
 
786
    rkV.IsConvex = (m_pkQuery->ToLine(iCurr,iPrev,iNext) > 0);
 
787
    return rkV.IsConvex;
 
788
}
 
789
//----------------------------------------------------------------------------
 
790
template <class Real>
 
791
bool TriangulateEC<Real>::IsEar (int i)
 
792
{
 
793
    Vertex& rkV = V(i);
 
794
 
 
795
    if (m_iRFirst == -1)
 
796
    {
 
797
        // The remaining polygon is convex.
 
798
        rkV.IsEar = true;
 
799
        return true;
 
800
    }
 
801
 
 
802
    // Search the reflex vertices and test if any are in the triangle
 
803
    // <V[prev],V[curr],V[next]>.
 
804
    int iPrev = V(rkV.VPrev).Index;
 
805
    int iCurr = rkV.Index;
 
806
    int iNext = V(rkV.VNext).Index;
 
807
    rkV.IsEar = true;
 
808
    for (int j = m_iRFirst; j != -1; j = V(j).SNext)
 
809
    {
 
810
        // Check if the test vertex is already one of the triangle vertices.
 
811
        if (j == rkV.VPrev || j == i || j == rkV.VNext)
 
812
        {
 
813
            continue;
 
814
        }
 
815
 
 
816
        // V[j] has been ruled out as one of the original vertices of the
 
817
        // triangle <V[prev],V[curr],V[next]>.  When triangulating polygons
 
818
        // with holes, V[j] might be a duplicated vertex, in which case it
 
819
        // does not affect the earness of V[curr].
 
820
        int iTest = V(j).Index;
 
821
        if (m_kSPositions[iTest] == m_kSPositions[iPrev]
 
822
        ||  m_kSPositions[iTest] == m_kSPositions[iCurr]
 
823
        ||  m_kSPositions[iTest] == m_kSPositions[iNext])
 
824
        {
 
825
            continue;
 
826
        }
 
827
 
 
828
        // Test if the vertex is inside or on the triangle.  When it is, it
 
829
        // causes V[curr] not to be an ear.
 
830
        if (m_pkQuery->ToTriangle(iTest,iPrev,iCurr,iNext) <= 0)
 
831
        {
 
832
            rkV.IsEar = false;
 
833
            break;
 
834
        }
 
835
    }
 
836
 
 
837
    return rkV.IsEar;
 
838
}
 
839
//----------------------------------------------------------------------------
 
840
template <class Real>
 
841
void TriangulateEC<Real>::InsertAfterC (int i)
 
842
{
 
843
    if (m_iCFirst == -1)
 
844
    {
 
845
        // add first convex vertex
 
846
        m_iCFirst = i;
 
847
    }
 
848
    else
 
849
    {
 
850
        V(m_iCLast).SNext = i;
 
851
        V(i).SPrev = m_iCLast;
 
852
    }
 
853
    m_iCLast = i;
 
854
}
 
855
//----------------------------------------------------------------------------
 
856
template <class Real>
 
857
void TriangulateEC<Real>::InsertAfterR (int i)
 
858
{
 
859
    if (m_iRFirst == -1)
 
860
    {
 
861
        // add first reflex vertex
 
862
        m_iRFirst = i;
 
863
    }
 
864
    else
 
865
    {
 
866
        V(m_iRLast).SNext = i;
 
867
        V(i).SPrev = m_iRLast;
 
868
    }
 
869
    m_iRLast = i;
 
870
}
 
871
//----------------------------------------------------------------------------
 
872
template <class Real>
 
873
void TriangulateEC<Real>::InsertEndE (int i)
 
874
{
 
875
    if (m_iEFirst == -1)
 
876
    {
 
877
        // add first ear
 
878
        m_iEFirst = i;
 
879
        m_iELast = i;
 
880
    }
 
881
    V(m_iELast).ENext = i;
 
882
    V(i).EPrev = m_iELast;
 
883
    m_iELast = i;
 
884
}
 
885
//----------------------------------------------------------------------------
 
886
template <class Real>
 
887
void TriangulateEC<Real>::InsertAfterE (int i)
 
888
{
 
889
    Vertex& rkVFirst = V(m_iEFirst);
 
890
    int iCurrENext = rkVFirst.ENext;
 
891
    Vertex& rkV = V(i);
 
892
    rkV.EPrev = m_iEFirst;
 
893
    rkV.ENext = iCurrENext;
 
894
    rkVFirst.ENext = i;
 
895
    V(iCurrENext).EPrev = i;
 
896
}
 
897
//----------------------------------------------------------------------------
 
898
template <class Real>
 
899
void TriangulateEC<Real>::InsertBeforeE (int i)
 
900
{
 
901
    Vertex& rkVFirst = V(m_iEFirst);
 
902
    int iCurrEPrev = rkVFirst.EPrev;
 
903
    Vertex& rkV = V(i);
 
904
    rkV.EPrev = iCurrEPrev;
 
905
    rkV.ENext = m_iEFirst;
 
906
    rkVFirst.EPrev = i;
 
907
    V(iCurrEPrev).ENext = i;
 
908
}
 
909
//----------------------------------------------------------------------------
 
910
template <class Real>
 
911
void TriangulateEC<Real>::RemoveV (int i)
 
912
{
 
913
    int iCurrVPrev = V(i).VPrev;
 
914
    int iCurrVNext = V(i).VNext;
 
915
    V(iCurrVPrev).VNext = iCurrVNext;
 
916
    V(iCurrVNext).VPrev = iCurrVPrev;
 
917
}
 
918
//----------------------------------------------------------------------------
 
919
template <class Real>
 
920
int TriangulateEC<Real>::RemoveE (int i)
 
921
{
 
922
    int iCurrEPrev = V(i).EPrev;
 
923
    int iCurrENext = V(i).ENext;
 
924
    V(iCurrEPrev).ENext = iCurrENext;
 
925
    V(iCurrENext).EPrev = iCurrEPrev;
 
926
    return iCurrENext;
 
927
}
 
928
//----------------------------------------------------------------------------
 
929
template <class Real>
 
930
void TriangulateEC<Real>::RemoveR (int i)
 
931
{
 
932
    assert(m_iRFirst != -1 && m_iRLast != -1);
 
933
 
 
934
    if (i == m_iRFirst)
 
935
    {
 
936
        m_iRFirst = V(i).SNext;
 
937
        if (m_iRFirst != -1)
 
938
        {
 
939
            V(m_iRFirst).SPrev = -1;
 
940
        }
 
941
        V(i).SNext = -1;
 
942
    }
 
943
    else if (i == m_iRLast)
 
944
    {
 
945
        m_iRLast = V(i).SPrev;
 
946
        if (m_iRLast != -1)
 
947
        {
 
948
            V(m_iRLast).SNext = -1;
 
949
        }
 
950
        V(i).SPrev = -1;
 
951
    }
 
952
    else
 
953
    {
 
954
        int iCurrSPrev = V(i).SPrev;
 
955
        int iCurrSNext = V(i).SNext;
 
956
        V(iCurrSPrev).SNext = iCurrSNext;
 
957
        V(iCurrSNext).SPrev = iCurrSPrev;
 
958
        V(i).SNext = -1;
 
959
        V(i).SPrev = -1;
 
960
    }
 
961
}
 
962
//----------------------------------------------------------------------------
 
963
 
 
964
//----------------------------------------------------------------------------
 
965
// Tree support.
 
966
//----------------------------------------------------------------------------
 
967
template <class Real>
 
968
void TriangulateEC<Real>::Delete (Tree*& rpkRoot)
 
969
{
 
970
    if (rpkRoot)
 
971
    {
 
972
        std::queue<Tree*> kQueue;
 
973
        kQueue.push(rpkRoot);
 
974
 
 
975
        while (kQueue.size() > 0)
 
976
        {
 
977
            Tree* pkTree = kQueue.front();
 
978
            kQueue.pop();
 
979
            for (int i = 0; i < (int)pkTree->Child.size(); i++)
 
980
            {
 
981
               kQueue.push(pkTree->Child[i]);
 
982
            }
 
983
            WM4_DELETE pkTree;
 
984
        }
 
985
 
 
986
        rpkRoot = 0;
 
987
    }
 
988
}
 
989
//----------------------------------------------------------------------------
 
990
template <class Real>
 
991
int TriangulateEC<Real>::GetExtraElements (const Tree* pkTree)
 
992
{
 
993
    int iExtraElements = 0;
 
994
 
 
995
    std::queue<const Tree*> kQueue;
 
996
    kQueue.push(pkTree);
 
997
    while (kQueue.size() > 0)
 
998
    {
 
999
        const Tree* pkRoot = kQueue.front();
 
1000
        kQueue.pop();
 
1001
        int iNumChildren = (int)pkRoot->Child.size();
 
1002
        iExtraElements += 2*iNumChildren;
 
1003
 
 
1004
        for (int i = 0; i < iNumChildren; i++)
 
1005
        {
 
1006
            const Tree* pkChild = pkRoot->Child[i];
 
1007
            int iNumGrandChildren = (int)pkChild->Child.size();
 
1008
            for (int j = 0; j < iNumGrandChildren; j++)
 
1009
            {
 
1010
                kQueue.push(pkChild->Child[j]);
 
1011
            }
 
1012
        }
 
1013
    }
 
1014
 
 
1015
    return iExtraElements;
 
1016
}
 
1017
//----------------------------------------------------------------------------
 
1018
 
 
1019
//----------------------------------------------------------------------------
 
1020
// explicit instantiation
 
1021
//----------------------------------------------------------------------------
 
1022
template WM4_FOUNDATION_ITEM
 
1023
class TriangulateEC<float>;
 
1024
 
 
1025
template WM4_FOUNDATION_ITEM
 
1026
class TriangulateEC<double>;
 
1027
//----------------------------------------------------------------------------
 
1028
}