source: GTP/trunk/Lib/Vis/Preprocessing/src/Polygon3.cpp @ 1420

Revision 1420, 12.8 KB checked in by mattausch, 18 years ago (diff)

corrected raycasting bug for triangles because of ill defined triangles

RevLine 
[235]1#include "Polygon3.h"
2#include "Mesh.h"
[372]3#include "AxisAlignedBox3.h"
4#include "Ray.h"
[503]5#include "Triangle3.h"
[372]6
[860]7
[863]8namespace GtpVisibilityPreprocessor {
[860]9
10
[372]11Polygon3::Polygon3():
[1344]12mMaterial(NULL),
13mParent(NULL),
14mPiercingRays(0)
[1076]15, mPlane(NULL)
[574]16{
17}
[372]18
19Polygon3::Polygon3(const VertexContainer &vertices):
[1344]20mMaterial(NULL),
21mParent(NULL),
22mPiercingRays(0)
[1076]23, mPlane(NULL)
[574]24{
25        mVertices.reserve(vertices.size());
26        mVertices = vertices;
27}
[372]28
[574]29
[372]30Polygon3::Polygon3(MeshInstance *parent):
[1344]31mMaterial(NULL),
32mParent(parent),
33mPiercingRays(0)
[1076]34, mPlane(NULL)
[372]35{}
36
[574]37
[372]38Polygon3::Polygon3(Face *face, Mesh *parentMesh):
[1344]39mMaterial(NULL),
40mParent(NULL),
41mPiercingRays(0),
42mPlane(NULL)
[372]43{       
[574]44        mVertices.reserve(face->mVertexIndices.size());
[1344]45        VertexIndexContainer::iterator it, it_end = face->mVertexIndices.end();
46
[1420]47        for (it = face->mVertexIndices.begin(); it != it_end; ++ it)
[372]48        {
49                mVertices.push_back(parentMesh->mVertices[*it]);
50        }
[574]51
52        mMaterial = parentMesh->mMaterial;
[372]53}
54
[574]55
[1404]56Polygon3::Polygon3(const Triangle3 &tri):
57mMaterial(NULL),
58mParent(NULL),
59mPiercingRays(0),
60mPlane(NULL)
61{       
62        mVertices.reserve(3);
63       
64        mVertices.push_back(tri.mVertices[0]);
65        mVertices.push_back(tri.mVertices[1]);
66        mVertices.push_back(tri.mVertices[2]);
[1418]67
[1419]68        VertexContainer::iterator it, it_end = mVertices.end();
69
[1420]70        //for (it = mVertices.begin(); it != it_end;  ++ it){cout << (*it) << " ";}cout<<endl<<endl;
[1418]71        //cout << "poly n: " << GetNormal() << " t normal: " << tri.GetNormal() << endl;
[1404]72}
73
74
[1076]75Plane3 Polygon3::GetSupportingPlane()// const
[372]76{
[1076]77#if 0
[372]78        return Plane3(mVertices[0], mVertices[1], mVertices[2]);
[1076]79#else
80        if (!mPlane)
[1344]81        {
[1076]82                mPlane = new Plane3(mVertices[0], mVertices[1], mVertices[2]);
[1344]83        }
84
[1076]85        return *mPlane;
86#endif
[372]87}
88
[574]89
[306]90Vector3 Polygon3::GetNormal() const
91{
[321]92    return Normalize(CrossProd(mVertices[2] - mVertices[1],
93                                                           mVertices[0] - mVertices[1]));
[372]94}
95
[1344]96
[372]97void Polygon3::Split(const Plane3 &partition,
98                                         Polygon3 &front,
[448]99                                         Polygon3 &back,
100                                         const float epsilon)
[372]101{
102        Vector3 ptA = mVertices.back();
[448]103        int sideA = partition.Side(ptA, epsilon);
[372]104        bool foundSplit = false;
[1344]105    Vector3 lastSplit;
106               
107        //-- find line - plane intersections
[448]108
[1344]109        VertexContainer::const_iterator it, it_end = mVertices.end();
110
111        for (it = mVertices.begin(); it != it_end; ++ it)
[372]112        {
113                Vector3 ptB = *it;
[448]114                int sideB = partition.Side(ptB, epsilon);
[372]115       
116                // vertices on different sides => split
[237]117            if (sideB > 0)
118                {
119                        if (sideA < 0)
[372]120                        {
121                                //-- plane - line intersection
122                                Vector3 splitPt = partition.FindIntersection(ptA, ptB);
123                       
124                                // test if split point not too close to previous split point
[448]125                                if (!foundSplit || (!EpsilonEqualV3(splitPt, lastSplit, epsilon)))
[372]126                                {
127                                        // add vertex to both polygons
128                                        front.mVertices.push_back(splitPt);
129                                        back.mVertices.push_back(splitPt);
130                                       
[436]131                                        lastSplit = splitPt;
[372]132                                        foundSplit = true;
133                                }
134                        }
135                        front.mVertices.push_back(ptB);
136                }
[237]137                else if (sideB < 0)
138                {
139                        if (sideA > 0)
[372]140                        {
141                                //-- plane - line intersection
142                                Vector3 splitPt = partition.FindIntersection(ptA, ptB);
143                                // test if split point not too close to other split point
144                                        // test if split point not too close to previous split point
[448]145                                if (!foundSplit || (!EpsilonEqualV3(splitPt, lastSplit, epsilon)))
[372]146                                {
147                                        // add vertex to both polygons
148                                        front.mVertices.push_back(splitPt);
149                                        back.mVertices.push_back(splitPt);
150
[436]151                                        lastSplit = splitPt;
[372]152                                        foundSplit = true;
153                                }       
154                        }
155                        back.mVertices.push_back(ptB);
156                }
157                else
158                {
159                        // vertex on plane => add vertex to both polygons
160                        front.mVertices.push_back(ptB);
161                        back.mVertices.push_back(ptB);
162                }
163       
164                ptA = ptB;
165                sideA = sideB;
166        }
167}
168
[574]169
[372]170float Polygon3::GetArea() const
171{
[306]172        Vector3 v = CrossProd(mVertices.back(), mVertices.front());
173   
174    for (int i=0; i < mVertices.size() - 1; ++i)
175                v += CrossProd(mVertices[i], mVertices[i+1]);
176   
[544]177    return 0.5f * fabs(DotProd(GetNormal(), v));
[372]178}
[235]179
[574]180
[448]181int Polygon3::Side(const Plane3 &plane,
182                                   const float epsilon) const
[372]183{
[448]184        int classification = ClassifyPlane(plane, epsilon);
[372]185       
186        if (classification == BACK_SIDE)
187                return -1;
188        else if (classification == FRONT_SIDE)
189                return 1;
[538]190        // plane splits polygon
[372]191        return 0;
192}
193
[574]194
[448]195int Polygon3::ClassifyPlane(const Plane3 &plane,
196                                                        const float epsilon) const
[372]197{
[574]198        VertexContainer::const_iterator it, it_end = mVertices.end();
[372]199
200        bool onFrontSide = false;
201        bool onBackSide = false;
202       
203        int count = 0;
204        // find possible line-plane intersections
[574]205        for (it = mVertices.begin(); it != it_end; ++ it)
[372]206        {
[448]207                const int side = plane.Side(*it, epsilon);
[372]208               
209                if (side > 0)
[683]210                {
[372]211                        onFrontSide = true;
[683]212                }
[372]213                else if (side < 0)
[683]214                {
[372]215                        onBackSide = true;
[683]216                }
217
[372]218                //TODO: check if split goes through vertex
219                if (onFrontSide && onBackSide) // split
220                {
221                        return SPLIT;
222                }
223                // 3 vertices enough to decide coincident
224                else if (((++ count) >= 3) && !onFrontSide && !onBackSide)
225                {   
226                        return COINCIDENT;
227                }
228        }
229
230        if (onBackSide)
231        {
232                return BACK_SIDE;
233        }
234        else if (onFrontSide)
235        {
236                return FRONT_SIDE;
237        }
238       
239        return COINCIDENT; // plane and polygon are coincident
240}
241
242
243Vector3
244Polygon3::Center() const
245{
246        int i;
247        Vector3 sum = mVertices[0];
248        for (i=1; i < mVertices.size(); i++)
249                sum += mVertices[i];
250       
251        return sum/(float)i;
252}
253
254
255void
256Polygon3::Scale(const float scale)
257{
258        int i;
259        Vector3 center = Center();
260        for (i=0; i < mVertices.size(); i++) {
261                mVertices[i] = center + scale*(mVertices[i] - center);
262        }
263}
264
[1420]265
[448]266bool Polygon3::Valid(const float epsilon) const
[372]267{
268        if (mVertices.size() < 3)
269                return false;
270
[1420]271        // matt: removed for performance issues
272#if _DEBUG
273        if (0)
[372]274        {
[574]275                // check if area exceeds certain size
276                if (GetArea() < AREA_LIMIT)
277                {
278                        Debug << "area too small: " << GetArea() << endl;
279                        return false;
280                }
[372]281        }
[645]282       
283        if (1)
[574]284        {
285                Vector3 vtx = mVertices.back();
286                VertexContainer::const_iterator it, it_end = mVertices.end();
[372]287
[574]288                for (it = mVertices.begin(); it != it_end; ++it)
[372]289                {
[1420]290                        if (EpsilonEqualV3(vtx, *it, epsilon))
[574]291                        {
[1420]292                                //cout << "Malformed vertices:\n" << *this << endl;
[574]293                                return false;
294                        }
295                        vtx = *it;
[372]296                }
297        }
[574]298#endif
299
[372]300        return true;
301}
302
[574]303
[314]304// int_lineseg returns 1 if the given line segment intersects a 2D
305// ray travelling in the positive X direction.  This is used in the
306// Jordan curve computation for polygon intersection.
307inline int
308int_lineseg(float px,
309            float py,
310            float u1,
311            float v1,
312            float u2,
313            float v2)
314{
315        float t;
316        float ydiff;
317
318        u1 -= px; u2 -= px;       // translate line
319        v1 -= py; v2 -= py;
320
321        if ((v1 > 0 && v2 > 0) ||
322                (v1 < 0 && v2 < 0) ||
323                (u1 < 0 && u2 < 0))
324        return 0;
325
326        if (u1 > 0 && u2 > 0)
327        return 1;
328
329        ydiff = v2 - v1;
330
331        if (fabs(ydiff) < Limits::Small)
332        {         // denominator near 0
333                if (((fabs(v1) > Limits::Small) || (u1 > 0) || (u2 > 0)))
334                        return 0;
335                return 1;
336        }
337
338        t = -v1 / ydiff;                  // Compute parameter
339
340        return (u1 + t * (u2 - u1)) > 0;
341}
342
[574]343
[314]344int Polygon3::CastRay(const Ray &ray, float &t, const float nearestT)
345{
[1328]346        const Plane3 plane = GetSupportingPlane();
347        const float dot = DotProd(plane.mNormal, ray.GetDir());
[314]348
349        // Watch for near-zero denominator
350        // ONLY single sided polygons!!!!!
351        if (dot > -Limits::Small)
[1328]352        {
353                //  if (fabs(dot) < Limits::Small)
354                return Ray::NO_INTERSECTION;
355        }
[314]356
357        t = (-plane.mD - DotProd(plane.mNormal, ray.GetLoc())) / dot;
358
359        if (t <= Limits::Small)
[1328]360        {
361                return Ray::INTERSECTION_OUT_OF_LIMITS;
362        }
[314]363
[1328]364        if (t >= nearestT)
365        {
366                return Ray::INTERSECTION_OUT_OF_LIMITS; // no intersection was found
[314]367        }
368
369        int count = 0;
370        float u, v, u1, v1, u2, v2;
371        int i;
372
[1328]373        const int paxis = plane.mNormal.DrivingAxis();
[314]374
375        // Project the intersection point onto the coordinate plane
376        // specified by which.
377        ray.Extrap(t).ExtractVerts(&u, &v, paxis);
378
[1328]379        const int size = (int)mVertices.size();
[314]380
381        mVertices.back().ExtractVerts(&u1, &v1, paxis );
382
383        if (0 && size <= 4)
384        {
385                // assume a convex face
386                for (i = 0; i < size; i++)
387                {
388                mVertices[i].ExtractVerts(&u2, &v2, paxis);
[1328]389           
[314]390                        // line u1, v1, u2, v2
[1328]391                        if ((v2 - v1) * (u1 - u) > (u2 - u1)*(v1 - v))
[314]392                                return Ray::NO_INTERSECTION;
393
394                        u1 = u2;
395                        v1 = v2;
[1328]396                }
397               
398                return Ray::INTERSECTION;
[314]399        }
400
401        // We're stuck with the Jordan curve computation.  Count number
402        // of intersections between the line segments the polygon comprises
403        // with a ray originating at the point of intersection and
404        // travelling in the positive X direction.
405        for (i = 0; i < size; i++)
406        {
407                mVertices[i].ExtractVerts(&u2, &v2, paxis);
408
409                count += (int_lineseg(u, v, u1, v1, u2, v2) != 0);
410
411                u1 = u2;
412                v1 = v2;
413        }
414
415        // We hit polygon if number of intersections is odd.
416        return (count & 1) ? Ray::INTERSECTION : Ray::NO_INTERSECTION;
417}
418
[574]419
[372]420void Polygon3::InheritRays(Polygon3 &front_piece,
421                                               Polygon3 &back_piece) const
422{
423        if (mPiercingRays.empty())
424                return;
425
426        RayContainer::const_iterator it,
427                it_end = mPiercingRays.end();
428
429        for (it = mPiercingRays.begin(); it != it_end; ++ it)
430        {
431                switch((*it)->GetId())
432                {
433                case Ray::BACK:
434                        back_piece.mPiercingRays.push_back(*it);
435                        break;
436                case Ray::FRONT:
437                        front_piece.mPiercingRays.push_back(*it);
438                        break;
439                case Ray::FRONT_BACK:
440                        back_piece.mPiercingRays.push_back(*it);
441                        break;
442                case Ray::BACK_FRONT:
443                        front_piece.mPiercingRays.push_back(*it);
444                        break;
445                default:
446                        break;
447                }
448        }
449}
450
[574]451
[448]452int Polygon3::ClassifyPlane(const PolygonContainer &polys,
453                                                        const Plane3 &plane,
454                                                        const float epsilon)
[372]455{
[683]456        bool onFrontSide = false;
457        bool onBackSide = false;
[372]458        PolygonContainer::const_iterator it;
459
460        // find intersections
461        for (it = polys.begin(); it != polys.end(); ++ it)
462        {
[448]463        const int cf = (*it)->ClassifyPlane(plane, epsilon);
[372]464               
465        if (cf == FRONT_SIDE)
466                        onFrontSide = true;
467                else if (cf == BACK_SIDE)
468                        onBackSide = true;
469               
470                if ((cf == SPLIT) || (cf == COINCIDENT) || (onFrontSide && onBackSide))
471                {
472                        return SPLIT;
473                }
474        }
475
476        if (onBackSide)
477        {
478                return BACK_SIDE;
479        }
480        else if (onFrontSide)
481        {
482                return FRONT_SIDE;
483        }
484       
485        return SPLIT;
486}
487
[574]488
[375]489int Polygon3::ParentObjectsSize(const PolygonContainer &polys)
[372]490{
491        int count = 0;
492
493        PolygonContainer::const_iterator it, it_end = polys.end();
494
495        MeshInstance::NewMail();
496
497        for (it = polys.begin(); it != it_end; ++ it)
498        {
499                if ((*it)->mParent && !(*it)->mParent->Mailed())
500                {
501                        ++ count;
502                        (*it)->mParent->Mail();
503                }
504        }
505        return count;
[384]506}
507
[574]508
[384]509float Polygon3::GetArea(const PolygonContainer &cell)
510{
511        float area = 0;
512        PolygonContainer::const_iterator pit;
513
514        for (pit = cell.begin(); pit != cell.end(); ++ pit)
[448]515        {
[384]516                area += (*pit)->GetArea();
[448]517        }
[396]518
[384]519        return area;
[396]520}
521
522
523Polygon3 *Polygon3::CreateReversePolygon() const
524{
525        Polygon3 *revPoly = new Polygon3();
526
527        VertexContainer::const_reverse_iterator rit,
528                        rit_end = mVertices.rend();
529
530        for(rit = mVertices.rbegin(); rit != rit_end; ++ rit)
[448]531                revPoly->mVertices.push_back(*rit);
532
533        return revPoly;
[503]534}
535
536
[840]537void Polygon3::Triangulate(vector<Triangle3> &triangles) const
[503]538{
539        int i = 1;
540        int j = 0;
[508]541        int k = (int)mVertices.size() - 1;
[503]542        int count = 0;
543
544        while (i < k)
545        {
[1418]546                if (0) // matt: I don't know why i made it like this because the normals are wrong!!
547                        triangles.push_back(Triangle3(mVertices[i], mVertices[j], mVertices[k]));
548                else
549                        triangles.push_back(Triangle3(mVertices[k], mVertices[j], mVertices[i]));
[503]550
551                if ((count ++) % 2)
552                        j = i ++;
553                else
554                        j = k --;
555        }
556}
[508]557
558
[840]559void Polygon3::Triangulate(VertexIndexContainer &indices) const
[508]560{
561        int i = 1;
562        int j = 0;
563        int k = (int)mVertices.size() - 1;
564
565        int count = 0;
566
567        while (i < k)
568        {
569                indices.push_back(k);
570                indices.push_back(j);
571                indices.push_back(i);
572
573                if ((count ++) % 2)
[1328]574                {
[508]575                        j = i ++;
[1328]576                }
[508]577                else
[1328]578                {
[508]579                        j = k --;
[1328]580                }
[508]581        }
582}
583
584
[840]585void IncludePolyInMesh(const Polygon3 &poly, Mesh &mesh)
[508]586{
587        const int n = (int)mesh.mVertices.size();
588
[1420]589
[1419]590        /////////////
591    //-- add the vertices to the mesh
592
[840]593    VertexContainer::const_iterator vit, vit_end = poly.mVertices.end();
594        for (vit = poly.mVertices.begin(); vit != vit_end; ++ vit)
[508]595        {
596                mesh.mVertices.push_back(*vit);
597        }
598
599        // one quad => no triangulation necessary
[840]600        if (poly.mVertices.size() == 4)
[508]601        {
602                mesh.AddFace(new Face(n, n + 1, n + 2, n + 3));
[840]603                return;
[508]604        }
[840]605       
606        VertexIndexContainer indices;
607        poly.Triangulate(indices);
[508]608
[1420]609        //if (indices.size() < 3) return; // something is wrong
610       
[840]611        // add indices of triangle strip
[1420]612        for (int i = 0; i < (int)indices.size(); i += 3)
[840]613        {
[1420]614                Face *face = new Face(indices[i] + n,
615                                                          indices[i + 1] + n,
616                                                          indices[i + 2] + n);
[840]617                mesh.AddFace(face);
[508]618        }
[860]619}
620
[1328]621
622AxisAlignedBox3 Polygon3::GetBoundingBox() const
623{
624        VertexContainer::const_iterator vit, vit_end = mVertices.end();
625       
626        AxisAlignedBox3 box;
627        box.Initialize();
628
629        for (vit = mVertices.begin(); vit != vit_end; ++ vit)
630        {
631                box.Include(*vit);
632        }
633
634        return box;
635}
636
637
[508]638}
Note: See TracBrowser for help on using the repository browser.