Polygon.cc

Go to the documentation of this file.

#include <esg/Definitions.h>
#include <esg/geometry/Polygon.h>
#include <esg/mesh/PolygonMesh.h>

using namespace esg;

Mesh* esg::Polygon::_mesh(int) const
{
    Mesh *tmp_polyg;
#if defined(_MSC_VER)
    unsigned *indices = new unsigned [_nVert];
#else
    unsigned indices [_nVert];
#endif

    for (unsigned i = 0; i < _nVert; i++) indices[i] = i;
    
    if (haveVertNormals()) {
        tmp_polyg = (Mesh*) new PolygonMesh(_vertices,
                                            indices,
                                            _nVert,
                                            _normals,
                                            indices,
                                            _normal,
                                            _fxyz);
#if defined(_MSC_VER)
        delete [] indices;
#endif
        return tmp_polyg;
    } else {
        tmp_polyg = (Mesh*) new PolygonMesh(_vertices,
                                            indices,
                                            _nVert,
                                            NULL,
                                            NULL,
                                            _normal,
                                            _fxyz);
#if defined(_MSC_VER)
        delete [] indices;
#endif
        return tmp_polyg;
    }
}
 
void esg::Polygon::_duplicate_attributes(const Geometry& src)
{
    Geometry::_duplicate_attributes(src);

    _nVert = ((Polygon&)src)._nVert;
    
    _vertices = new Vector3 [_nVert];
    for (unsigned i = 0; i < _nVert; i++)
        _vertices[i].set(((Polygon&)src)._vertices[i]);
    
    if (((Polygon&)src)._normals) {
        _normals = new Vector3 [_nVert];
        for (unsigned i = 0; i < _nVert; i++)
            _normals[i].set(((Polygon&)src)._normals[i]);
    } else
        _normals = NULL;

    if (((Polygon&)src)._uvCoords) {
        _uvCoords = new Vector2 [_nVert];
        for (unsigned i = 0; i < _nVert; i++)
            _uvCoords[i].set(((Polygon&)src)._uvCoords[i]);
    } else
        _uvCoords = NULL;
    
    _normal.set(((Polygon&)src)._normal);
    _normalFixed = ((Polygon&)src)._normalFixed;
    
    _fxyz  = ((Polygon&)src)._fxyz;

    _edgeProjection = ((Polygon&)src)._edgeProjection;

    _proj = NULL;
    _area = ((Polygon&)src)._area;
}

Vector3 esg::Polygon::_get_facet_normal() const
{
    Vector3 result; // = (0,0,0)
    result.cross(_get_vertex(1) - _get_vertex(0),
                 _get_vertex(_nVert-1) - _get_vertex(0));
    result.normalize();
    return result;
}

Vertex3 esg::Polygon::_get_vertex(unsigned index) const 
{
    return _vertices[index];
}

Vertex3 esg::Polygon::_get_vert_normal(unsigned index) const 
{
    return _normals[index];
}

Vertex2 esg::Polygon::_get_vert_uv_coord(unsigned index) const
{
    return _uvCoords[index];
}

void esg::Polygon::_set_edge_projection(void)
{
    // find projection plane
    Vector3 absNorm(_normal);
    absNorm.absolute();
    if (absNorm.x >= absNorm.y && absNorm.x >= absNorm.z)
        _edgeProjection = YZ;
    else if (absNorm.y >= absNorm.x && absNorm.y >= absNorm.z)
        _edgeProjection = XZ;
    else
        _edgeProjection = XY;
}

double (*esg::Polygon::_precompute_proj(void)) [6]
{
    double (*proj)[6] = new double [_nVert-2][6];
    double aux;
    Vector3 v0(_get_vertex(0));
    for (register unsigned i = 2; i < _nVert; i++) {
        switch (_edgeProjection) {
        case XY:
            proj[i-2][0] = _get_vertex(i-1).x - v0.x;
            proj[i-2][1] = _get_vertex(i-1).y - v0.y;
            proj[i-2][2] = _get_vertex(i).x   - v0.x;
            proj[i-2][3] = _get_vertex(i).y   - v0.y;
            break;
        case YZ:
            proj[i-2][0] = _get_vertex(i-1).y - v0.y;
            proj[i-2][1] = _get_vertex(i-1).z - v0.z;
            proj[i-2][2] = _get_vertex(i).y   - v0.y;
            proj[i-2][3] = _get_vertex(i).z   - v0.z;
            break;
        case XZ:
            proj[i-2][0] = _get_vertex(i-1).x - v0.x;
            proj[i-2][1] = _get_vertex(i-1).z - v0.z;
            proj[i-2][2] = _get_vertex(i).x   - v0.x;
            proj[i-2][3] = _get_vertex(i).z   - v0.z;
            break;
        case NONE_PROJ:
            delete [] proj;
            return NULL;
        }
        aux = proj[i-2][3] * proj[i-2][0] - proj[i-2][2] * proj[i-2][1];
        proj[i-2][4] = proj[i-2][0] / aux;
        proj[i-2][5] = proj[i-2][1] / aux;
    }

    return proj;
}

bool esg::Polygon::_separation_by_plane(const Geometry& geom) const
{
    Interval ext = geom.extent(_normal);
    if (ext.min <= _fxyz && _fxyz <= ext.max) return true;
    else return false;
}

bool esg::Polygon::_separation_by_edges(Geometry& geom) const
{
    PointEnv pointEnv1, pointEnv2;
    Vector3  ray1;
    Vector3  ray2;

    for (register unsigned i = 0; i < _nVert; i++) {
        Vector3 edge(_get_vertex((i+1)%_nVert) - _get_vertex(i));
        
        ray1.set(edge);
        ray1.normalize();
        ray2.set(ray1);
        ray2.negate();
        
        geom.rayIntersection(&pointEnv1,
                             ENV_WANT_INTERSECTION,
                             _get_vertex(i),
                             ray1);
        geom.rayIntersection(&pointEnv2,
                             ENV_WANT_INTERSECTION,
                             _get_vertex(i),
                             ray2);
        
        if ((pointEnv1.mask & ENV_HAVE_INTERSECTION) && // vertex i is inside
            (pointEnv2.mask& ENV_HAVE_INTERSECTION)) return false;

        if ((pointEnv1.mask & ENV_HAVE_INTERSECTION) &&
            (ray1.dot(pointEnv1.intersection) <= edge.length())) return false;

        if ((pointEnv2.mask & ENV_HAVE_INTERSECTION) &&
            (ray2.dot(pointEnv2.intersection) <= edge.length())) return false;
    }

    return true;
}

bool esg::Polygon::_point_inside_polygon(const Vector3& v, Vector3* n, Vector2* uv)
{
    /*
     * Check if point of intersection is inside the polygon's area
     *
     * Algorithm reference:
     *    Glassner A. S., Graphics Gems, Academic Press Professional,
     *    Cambridge, 1990, p. 390-393,735
     *
     * Here we use a little bit faster implementation.
     */

    float   alpha = 0.0, beta = 0.0;
    float   ux = 0.0, vx = 0.0;
    Vector3 v0(_get_vertex(0));

    switch (_edgeProjection) {
    case XY: ux = v.x - v0.x; vx = v.y - v0.y; break;
    case YZ: ux = v.y - v0.y; vx = v.z - v0.z; break;
    case XZ: ux = v.x - v0.x; vx = v.z - v0.z; break;
    case NONE_PROJ: return false;
    }

    if (!_proj) _proj = _precompute_proj();
    
    register unsigned i = 0;
    register double*  p = &_proj[0][0]; // fast access to precomputed values
    do {
        if (*p != 0.) {
            beta = vx * *(p+4) - ux * *(p+5);
            if ((beta >= 0.) && (beta <= 1.)) {
                alpha = (ux - beta * *(p+2)) / *p;
                if ((alpha >= 0.) && ((alpha+beta) <= 1.)) break;
            }
        } else {
            beta = ux / *(p+2);
            if ((beta >= 0.) && (beta <= 1.)) {
                alpha = (vx - beta * *(p+3)) / *(p+1);
                if ((alpha >= 0.) && ((alpha+beta) <= 1.)) break;
            }
        }
        if (++i >= _nVert-2) return false; // point is outside 
        p += 6;
    } while (true);

    // Tripple (vert[0], vert[i+1], vert[i+2]) determines triangle
    // including the point of intersection.

    if (n) { // interpolate normal
        n->set(_get_vert_normal(0)   * (1.0 - alpha - beta));
        n->add(_get_vert_normal(i+1) * alpha);
        n->add(_get_vert_normal(i+2) * beta);
    }

    if (uv) { // compute UV coord.  of the point of intersection
        uv->set(_get_vert_uv_coord(0)   * (1.0 - alpha - beta));
        uv->add(_get_vert_uv_coord(i+1) * alpha);
        uv->add(_get_vert_uv_coord(i+2) * beta);
    }

    return true;
}

/* this edge to edge test is based on Franlin Antonio's gem:
   "Faster Line Segment Intersection", in Graphics Gems III,
   pp. 199-202 */
#define ESG_EDGE_EDGE_TEST(V0,U0,U1)                  \
  switch (_edgeProjection) {                          \
  case Polygon::XY: Bx=U0.x-U1.x;By=U0.y-U1.y;Cx=V0.x-U0.x;Cy=V0.y-U0.y;break;\
  case Polygon::YZ: Bx=U0.y-U1.y;By=U0.z-U1.z;Cx=V0.y-U0.y;Cy=V0.z-U0.z;break;\
  case Polygon::XZ: Bx=U0.x-U1.x;By=U0.z-U1.z;Cx=V0.x-U0.x;Cy=V0.z-U0.z;break;\
  case Polygon::NONE_PROJ: return false;              \
  }                                                   \
  f = Ay*Bx-Ax*By;                                    \
  d = By*Cx-Bx*Cy;                                    \
  if((f>0 && d>=0 && d<=f) || (f<0 && d<=0 && d>=f))  \
  {                                                   \
    e=Ax*Cy-Ay*Cx;                                    \
    if(f>0)                                           \
    {                                                 \
      if(e>=0 && e<=f) return true;                   \
    }                                                 \
    else                                              \
    {                                                 \
      if(e<=0 && e>=f) return true;                   \
    }                                                 \
  }

#define ESG_EDGE_AGAINST_TRI_EDGES(V0,V1,U0,U1,U2)             \
{                                                              \
  double Ax=0.0,Ay=0.0,Bx=0.0,By=0.0,Cx=0.0,Cy=0.0,e,d,f;      \
  switch (p._edgeProjection) {                                 \
  case Polygon::XY: Ax = V1.x - V0.x; Ay = V1.y - V0.y; break; \
  case Polygon::YZ: Ax = V1.y - V0.y; Ay = V1.z - V0.z; break; \
  case Polygon::XZ: Ax = V1.x - V0.x; Ay = V1.z - V0.z; break; \
  case Polygon::NONE_PROJ: return false;                       \
  }                                                            \
  /* test edge U0,U1 against V0,V1 */                          \
  ESG_EDGE_EDGE_TEST(V0,U0,U1);                                \
  /* test edge U1,U2 against V0,V1 */                          \
  ESG_EDGE_EDGE_TEST(V0,U1,U2);                                \
  /* test edge U2,U1 against V0,V1 */                          \
  ESG_EDGE_EDGE_TEST(V0,U2,U0);                                \
}

bool esg::Polygon::_coplanar_tri_collision(Polygon&       p,
                                      const Vector3& v0,
                                      const Vector3& v1,
                                      const Vector3& v2,
                                      const Vector3& u0,
                                      const Vector3& u1,
                                      const Vector3& u2)
{
    ESG_EDGE_AGAINST_TRI_EDGES(v0,v1,u0,u1,u2);  
    ESG_EDGE_AGAINST_TRI_EDGES(v1,v2,u0,u1,u2);
    ESG_EDGE_AGAINST_TRI_EDGES(v2,v0,u0,u1,u2);
    
    /* finally, test if tri1 is totally contained in tri2 or vice versa */
    if (_point_inside_polygon(v0, NULL, NULL)) return true;
    if (p._point_inside_polygon(u0, NULL, NULL)) return true;
    return false;
}

#undef ESG_EDGE_EDGE_TEST
#undef ESG_AGAINST_TRI_EDGES


#define ESG_COMPUTE_INTERVALS(VV0,VV1,VV2,D0,D1,D2,D0D1,D0D2,A,B,C,X0,X1) \
{ \
        if(D0D1 > 0.0f) { \
            /* here we know that D0D2<=0.0 */ \
            /* that is D0, D1 are on the same side, D2 on the other */ \
            /* or on the plane */ \
            A=VV2; B=(VV0-VV2)*D2; C=(VV1-VV2)*D2; X0=D2-D0; X1=D2-D1; \
        } else if(D0D2>0.0f) { \
            /* here we know that d0d1<=0.0 */ \
            A=VV1; B=(VV0-VV1)*D1; C=(VV2-VV1)*D1; X0=D1-D0; X1=D1-D2; \
        } else if(D1*D2>0.0f || D0!=0.0f) { \
            /* here we know that d0d1<=0.0 or that D0!=0.0 */ \
            A=VV0; B=(VV1-VV0)*D0; C=(VV2-VV0)*D0; X0=D0-D1; X1=D0-D2; \
        } else if(D1!=0.0f) { \
            A=VV1; B=(VV0-VV1)*D1; C=(VV2-VV1)*D1; X0=D1-D0; X1=D1-D2; \
        } else if(D2!=0.0f) { \
            A=VV2; B=(VV0-VV2)*D2; C=(VV1-VV2)*D2; X0=D2-D0; X1=D2-D1; \
        } else { \
            /* triangles are coplanar */ \
            return _coplanar_tri_collision(p,                             \
                                           p2Vertex0,p2Vertex1,p2Vertex2, \
                                           p2Vertex0,p2Vertex1,p2Vertex2);\
        } \
}

#define ESG_SORT(a,b) if (a>b) { double c; c=a; a=b; b=c; }

int esg::Polygon::_triTriCollision(Polygon& p)
{
    if (_nVert < 3 || p.numVertices() < 3) return -1;
    
    double  du0,du1,du2,dv0,dv1,dv2;
    Vector3 D;
    Vertex3 p1Vertex0(_get_vertex(0));
    Vertex3 p1Vertex1(_get_vertex(1));
    Vertex3 p1Vertex2(_get_vertex(2));
    Vertex3 p2Vertex0(p._get_vertex(0));
    Vertex3 p2Vertex1(p._get_vertex(1));
    Vertex3 p2Vertex2(p._get_vertex(2));
    double  isect1[2], isect2[2];
    double  du0du1,du0du2,dv0dv1,dv0dv2;
    double  vp0,vp1,vp2;
    double  up0,up1,up2;

    // put U0,U1,U2 into plane equation 1 to compute signed
    // distances to the plane
    du0 = p._normal.dot(p1Vertex0) + p._fxyz;
    du1 = p._normal.dot(p1Vertex1) + p._fxyz;
    du2 = p._normal.dot(p1Vertex2) + p._fxyz;

    // coplanarity robustness check
    if (fabs(du0) < Geometry::EPS) du0 = 0.0;
    if (fabs(du1) < Geometry::EPS) du1 = 0.0;
    if (fabs(du2) < Geometry::EPS) du2 = 0.0;
    
    du0du1 = du0 * du1;
    du0du2 = du0 * du2;

    // same sign on all of them + not equal 0 => no intersection occurs
    if (du0du1 > 0.0 && du0du2 > 0.0) return false;

    // put V0,V1,V2 into plane equation 2 
    dv0 = _normal.dot(p2Vertex0) + _fxyz;
    dv1 = _normal.dot(p2Vertex1) + _fxyz;
    dv2 = _normal.dot(p2Vertex2) + _fxyz;

    if (fabs(dv0) < Geometry::EPS) dv0 = 0.0;
    if (fabs(dv1) < Geometry::EPS) dv1 = 0.0;
    if (fabs(dv2) < Geometry::EPS) dv2 = 0.0;

    dv0dv1 = dv0 * dv1;
    dv0dv2 = dv0 * dv2;

    // same sign on all of them + not equal 0 => no intersection occurs
    if (dv0dv1 > 0.0 && dv0dv2 > 0.0) return false;

    // compute direction of intersection line
    D.cross(p._normal, _normal);

    // compute and index to the largest component of D
    D.absolute();
    if (D.x > D.y) {
        if (D.x > D.z) {
            vp0 = p2Vertex0.x;
            vp1 = p2Vertex1.x;
            vp2 = p2Vertex2.x;
            up0 = p1Vertex0.x;
            up1 = p1Vertex1.x;
            up2 = p1Vertex2.x;
        } else {
            vp0 = p2Vertex0.z;
            vp1 = p2Vertex1.z;
            vp2 = p2Vertex2.z;
            up0 = p1Vertex0.z;
            up1 = p1Vertex1.z;
            up2 = p1Vertex2.z;
        }
    } else
        if (D.y > D.z) {
            vp0 = p2Vertex0.y;
            vp1 = p2Vertex1.y;
            vp2 = p2Vertex2.y;
            up0 = p1Vertex0.y;
            up1 = p1Vertex1.y;
            up2 = p1Vertex2.y;
        } else {
            vp0 = p2Vertex0.z;
            vp1 = p2Vertex1.z;
            vp2 = p2Vertex2.z;
            up0 = p1Vertex0.z;
            up1 = p1Vertex1.z;
            up2 = p1Vertex2.z;
        }
    
    // compute interval for triangle 1 
    double a,b,c,x0,x1;
    ESG_COMPUTE_INTERVALS(vp0,vp1,vp2,dv0,dv1,dv2,dv0dv1,dv0dv2,a,b,c,x0,x1);

    // compute interval for triangle 2 
    double d,e,f,y0,y1;
    ESG_COMPUTE_INTERVALS(up0,up1,up2,du0,du1,du2,du0du1,du0du2,d,e,f,y0,y1);

    double xx,yy,xxyy,tmp;
    xx   = x0 * x1;
    yy   = y0 * y1;
    xxyy = xx * yy;

    tmp       = a * xxyy;
    isect1[0] = tmp + b * x1 * yy;
    isect1[1] = tmp + c * x0 * yy;

    tmp       = d * xxyy;
    isect2[0] = tmp + e * xx * y1;
    isect2[1] = tmp + f * xx * y0;

    ESG_SORT(isect1[0],isect1[1]);
    ESG_SORT(isect2[0],isect2[1]);

    return ((isect1[1] < isect2[0] || isect2[1] < isect1[0]) ? false : true);
}

#undef ESG_COMPUTE_INTERVALS
#undef ESG_SORT


void esg::Polygon::_rotateX(float a)
{
    Matrix3d trMat;
    trMat.rotX(a);
    for (unsigned i = 0; i < _nVert; i++) {
        trMat.transform(_vertices[i]);
        if (haveVertNormals()) trMat.transform(_normals[i]);
    }
    trMat.transform(_normal);
    _set_edge_projection();
    if (_proj) { delete [] _proj; _proj = NULL; }
    _proj = _precompute_proj();
    _fxyz = - _normal.dot(_get_vertex(0));
}

void esg::Polygon::_rotateY(float a)
{
    Matrix3d trMat;
    trMat.rotY(a);
    for (unsigned i = 0; i < _nVert; i++) {
        trMat.transform(_vertices[i]);
        if (haveVertNormals()) trMat.transform(_normals[i]);
    }
    trMat.transform(_normal);
    _set_edge_projection();
    if (_proj) { delete [] _proj; _proj = NULL; }
    _fxyz = - _normal.dot(_get_vertex(0));
}

void esg::Polygon::_rotateZ(float a)
{
    Matrix3d trMat;
    trMat.rotZ(a);
    for (unsigned i = 0; i < _nVert; i++) {
        trMat.transform(_vertices[i]);
        if (haveVertNormals()) trMat.transform(_normals[i]);
    }
    trMat.transform(_normal);
    _set_edge_projection();
    if (_proj) { delete [] _proj; _proj = NULL; }
    _fxyz = - _normal.dot(_get_vertex(0));
}

void esg::Polygon::_rotate(float a, const Vector3& axis)
{
    Matrix4d trMat;
    Matrix3d rotMat;
    trMat.rotationGL(a, axis);
    trMat.get(rotMat);
    for (unsigned i = 0; i < _nVert; i++) {
        trMat.transform(_vertices[i]);
        if (haveVertNormals()) rotMat.transform(_normals[i]);
    }
    rotMat.transform(_normal);
    _set_edge_projection();
    if (_proj) { delete [] _proj; _proj = NULL; }
    _fxyz = - _normal.dot(_get_vertex(0));
}

void esg::Polygon::_rotate(const Matrix3& rotMat)
{
    for (unsigned i = 0; i < _nVert; i++) {
        rotMat.transform(_vertices[i]);
        if (haveVertNormals()) rotMat.transform(_normals[i]);
    }
    rotMat.transform(_normal);
    _set_edge_projection();
    if (_proj) { delete [] _proj; _proj = NULL; }
    _fxyz = - _normal.dot(_get_vertex(0));
}

void esg::Polygon::_translate(float x, float y, float z)
{
    for (unsigned i = 0; i < _nVert; i++) {
        _vertices[i].add(Vector3d(x,y,z));
    }
    _fxyz = - _normal.dot(_get_vertex(0));
}

void esg::Polygon::_transform(const Matrix4& trMat)
{
    Matrix3d rotMat;
    trMat.get(rotMat);
    for (unsigned i = 0; i < _nVert; i++) {
        trMat.transform(_vertices[i]);
        if (haveVertNormals()) rotMat.transform(_normals[i]);
    }
    rotMat.transform(_normal);
    _set_edge_projection();
    if (_proj) { delete [] _proj; _proj = NULL; }
    _fxyz = - _normal.dot(_get_vertex(0));
}

void esg::Polygon::_scale(float s)
{
    for (unsigned i = 0; i < _nVert; i++) _vertices[i].scale(s);
    if (_proj) { delete [] _proj; _proj = NULL; }
    //_proj = _precompute_proj();
    _fxyz *= s;
}



//------ public ------

esg::Polygon::Polygon()
{
    _vertices  = NULL;
    _normals   = NULL;
    _uvCoords  = NULL;
    _nVert     = 0;
    _fxyz      = 0;
    _area      = 0;
    _proj      = NULL;    
}

esg::Polygon::Polygon(const Vector3* va,
                      const Vector3* na,
                      const Vector2* uva,
                      unsigned       n)
{
    _vertices  = NULL;
    _normals   = NULL;
    _uvCoords  = NULL;
    _nVert     = n;
    _proj      = NULL;

    (setVertices(va, n) && setNormals(na) && setUVCoords(uva));
}

esg::Polygon::Polygon(const Vector3* va,
                      const Vector3* na,
                      const Vector2* uva,
                      unsigned       n,
                      const Vector3& norm,
                      bool           fixedn)
{
    _vertices  = NULL;
    _normals   = NULL;
    _uvCoords  = NULL;
    _nVert     = n;
    _proj      = NULL;

    if (setVertices(va, n) && setNormals(na) && setUVCoords(uva))
        setFacetNormal(norm, fixedn);
}

esg::Polygon::~Polygon()
{
    if (_vertices) delete [] _vertices;
    if (_normals)  delete [] _normals;
    if (_uvCoords) delete [] _uvCoords;
    if (_proj)     delete [] _proj;
}

void esg::Polygon::rayIntersection(PointEnv* pPE,
                              int            mask, 
                              const Vector3& origin,
                              const Vector3& direction,
                              float          maxDist)
{
    if (!pPE) return;
    
    pPE->mask = ENV_HAVE_NOTHING;

    if (_nVert < 3) return;

    /*
     * get point of intersection with support plane
     */

    double&  np1 = pPE->nd;
    double&  no  = pPE->no;
    double&  t   = pPE->distance;
    Vector3& v   = pPE->intersection;

    np1 = _normal.dot(direction);
    if (fabs(np1) < Geometry::EPS) return; // Parallel
    no = _normal.dot(origin);
    t = (no + _fxyz) / -np1;
    if (t < Geometry::EPS || t > maxDist) return; // Directs outward or too far

    v.set(direction); // point of intersection itself
    v.scaleAdd(t, origin);


    /*
     * Set output
     */

    Vector3* n  = (((mask&ENV_WANT_NORMAL)  && haveVertNormals())
                   ? &(pPE->normal): NULL);
    Vector2* uv = (((mask&ENV_WANT_UV_COORD) && haveVertUVCoords())
                   ? &(pPE->uvCoord) : NULL);

    if (!_point_inside_polygon(v, n, uv)) return;

    if (n) {
        pPE->normalOrientation = PointEnv::OUTWARDS_NORMAL;
        pPE->mask |= ENV_HAVE_NORMAL;
    } else {
        if (mask & ENV_WANT_NORMAL) {
            // intersection is found but the normal vector
            // cannot be interpolated  => set it to the facet normal
            pPE->normal.set(_normal);
            if (_normalFixed) pPE->normalOrientation = PointEnv::OUTWARDS_NORMAL;
            else              pPE->normalOrientation = PointEnv::RANDOM_NORMAL;
            pPE->mask |= ENV_HAVE_NORMAL|ENV_HAVE_N_DOT_O|ENV_HAVE_N_DOT_D;
        }
    }

    pPE->mask |= ((uv) ?
                  ENV_HAVE_INTERFERENCE|ENV_HAVE_N_DOT_O|
                  ENV_HAVE_INTERSECTION|ENV_HAVE_DISTANCE|ENV_HAVE_UV_COORD :
                  ENV_HAVE_INTERFERENCE|ENV_HAVE_N_DOT_O|
                  ENV_HAVE_INTERSECTION|ENV_HAVE_DISTANCE);
}

bool esg::Polygon::mapToUV(const Vector3& v, Vector2& uv)
{
    return _point_inside_polygon(v, NULL, &uv);
}

void esg::Polygon::randomSample(int mask, PointEnv& pe, double* pdf)
{
    pe.mask = ENV_HAVE_NOTHING;
    if (!(mask & ENV_WANT_SURFACE_POINT|ENV_WANT_UV_COORD|ENV_WANT_NORMAL))
        return;

    // get rectangle packing the plygon
    Vector3  u(_get_vertex(1)); u.sub(_get_vertex(0)); u.normalize();
    Vector3  v; v.cross(u, _normal); v.normalize();
    Interval uExt, vExt;
    double   ud, vd, p0u = 0, p0v = 0;
    for (register unsigned i = 0; i < _nVert; i++) {
        ud = u.dot(_get_vertex(i));
        vd = v.dot(_get_vertex(i));
        if (i == 0) { p0u = ud; p0v = vd; }
        if (ud > uExt.max) uExt.max = ud;
        if (ud < uExt.min) uExt.min = ud;
        if (vd > vExt.max) vExt.max = vd;
        if (vd < vExt.min) vExt.min = vd;
    }

    Vector3* n  = ((((mask&ENV_WANT_NORMAL)  && haveVertNormals())
                    ? &(pe.normal) : NULL));
    Vector2* uv = ((((mask&ENV_WANT_UV_COORD) && haveVertUVCoords())
                    ? &(pe.uvCoord) : NULL));
    Vector3  auxU, auxV;

    // set p0 to the bottom-left corner of rectangle
    auxU.set(u);
    auxV.set(v);
    auxU.scale(uExt.min - p0u);
    auxV.scale(vExt.min - p0v);
    Vector3 p0(_get_vertex(0));
    p0.add(auxU);
    p0.add(auxV);

    ud = uExt.max - uExt.min;
    vd = vExt.max - vExt.min;

    do { // rejection sampling
        auxU.set(u);
        auxU.scale(ESG_DBL_RAND * ud);
        auxV.set(v);
        auxV.scale(ESG_DBL_RAND * vd);
        pe.intersection.set(p0);
        pe.intersection.add(auxU);
        pe.intersection.add(auxV);
    } while(!_point_inside_polygon(pe.intersection, n, uv));

    if (n) {
        pe.normalOrientation = PointEnv::OUTWARDS_NORMAL;
        pe.mask |= ENV_HAVE_NORMAL;
    } else {
        if (mask & ENV_WANT_NORMAL) {
            // intersection is found but the normal vector
            // cannot be interpolated  => set it to the facet normal
            pe.normal.set(_normal);
            if (_normalFixed) pe.normalOrientation = PointEnv::OUTWARDS_NORMAL;
            else              pe.normalOrientation = PointEnv::RANDOM_NORMAL;
            pe.mask |= ENV_HAVE_NORMAL;
        }
    }
    
    pe.mask |= ((uv) ?
                ENV_HAVE_INTERFERENCE|ENV_HAVE_INTERSECTION|ENV_HAVE_UV_COORD :
                ENV_HAVE_INTERFERENCE|ENV_HAVE_INTERSECTION);

    if (pdf) *pdf = 1./_area;
}

bool esg::Polygon::randomDirection(const Vector3& pov, Vector3&, double* pdf) 
{
    if (_normalFixed && _normal.dot(pov) < _fxyz) return false;
    fprintf(stderr,"Polygon::randomDirection(): Not implemented\n");
    return false;
}

Interval esg::Polygon::extent(const Vector3& direction) const
{
    float    vertExt;
    Interval retInt;
    
    retInt.min = MAXFLOAT;
    //  retInt.max = MINFLOAT;
    retInt.max = - MAXFLOAT;
    
    for (register unsigned i = 0; i < _nVert; i++) {
        vertExt = direction.dot(_get_vertex(i));
        if (retInt.min > vertExt) retInt.min = vertExt;
        if (retInt.max < vertExt) retInt.max = vertExt;
    }
    return retInt;
}


Geometry* esg::Polygon::clone(const Matrix4* pTrMat) const
{
    Polygon* pRet = new Polygon;
    pRet->_duplicate_attributes(*this);
    if (pTrMat) pRet->_transform(*pTrMat);
    return pRet;
}

Vertex3 esg::Polygon::centroid() const
{
    Vertex3 pog(0.0, 0.0, 0.0);
    for (register unsigned i = 0; i < _nVert; i++) {
        Vertex3 aux(_get_vertex(i));
        aux.scale(1.0/_nVert);
        pog.add(aux);
    }
    return pog;
}

double esg::Polygon::radius(const Vector3& centroid) const
{
    double  r = -MAXDOUBLE;
    double  auxR;
    Vector3 auxV;
    for (register unsigned i = 0; i < _nVert; i++) {
        auxV.set(_get_vertex(i));
        auxV.sub(centroid);
        auxR = auxV.length();
        if (auxR > r) r = auxR;
    }
    return r;
}

bool esg::Polygon::separation(Geometry& geom, Vector3* pDir) 
{
    try {
        int c = _triTriCollision(dynamic_cast<Polygon&>(geom));
        if (c >= 0) return (!c);
    } catch (bad_cast) {
        // geom is not a polygon => continue
    }

    // either geom is not polygon or triTriCollision failed:
    
    if (_separation_by_plane(geom)) {
        if (pDir) pDir->set(_normal);
        return true;
    }

    if (!_separation_by_edges(geom)) return false;
    
    return false; // don't know => return false
}

double esg::Polygon::distance(const Geometry& geom, Vector3* pDir) 
{
    return MINDOUBLE; 
}

void esg::Polygon::dump(const char* intent, const char* tab)
{
    printf("%s geometry Polygon {\n", intent);
    printf("%s %s numVertices %d\n", intent, tab, _nVert);
    if (_normalFixed) 
        printf("%s %s normal %f %f %f\n", intent, tab, _normal.x, _normal.y, _normal.z);
    if (_vertices) {
        printf("%s %s vertex [\n", intent, tab);
        for (register unsigned i = 0; i < _nVert; i++)
            printf("%s %s %s %f %f %f\n", intent, tab, tab, _vertices[i].x, _vertices[i].y, _vertices[i].z);
        printf("%s %s ]\n", intent, tab);
    }
    if (haveVertNormals()) {
        printf("%s %s vertexNormal [\n", intent, tab);
        for (register unsigned i = 0; i < _nVert; i++)
            printf("%s %s %s %f %f %f\n", intent, tab, tab, _normals[i].x, _normals[i].y, _normals[i].z);
        printf("%s %s ]\n", intent, tab);
    }
    if (_uvCoords) {
        printf("%s %s uvCoord [\n", intent, tab);
        for (register unsigned i = 0; i < _nVert; i++)
            printf("%s %s %s %f %f\n", intent, tab, tab, _uvCoords[i].x, _uvCoords[i].y);
        printf("%s %s ]\n", intent, tab);
    }
    printf("%s }\n", intent);
}

Vector3 esg::Polygon::getFacetNormal() const
{
    return _normal;
}

bool esg::Polygon::haveFixedNormal() const
{
    return _normalFixed;
}

Vertex3 esg::Polygon::getVertex(unsigned index) const 
{
    return ((index < _nVert) ? _get_vertex(index) : Vertex3(0,0,0));
}

Vertex3 esg::Polygon::getVertNormal(unsigned index) const 
{
    return ((haveVertNormals() && index < _nVert) ?
            _get_vert_normal(index) :
            Vertex3(0,0,0));
}

bool esg::Polygon::haveVertNormals() const
{
    return (_normals != NULL);
}

Vertex2 esg::Polygon::getVertUVCoord(unsigned index) const 
{
    return ((haveVertUVCoords() && index < _nVert)
            ? _get_vert_uv_coord(index)
            : Vertex2(0,0));
}

bool esg::Polygon::haveVertUVCoords() const
{
    return (_uvCoords != NULL);
}

unsigned esg::Polygon::numVertices() const
{
    return _nVert;
}

double esg::Polygon::getFXYZ() const
{
    return _fxyz;
}

esg::Polygon::edgeProjection esg::Polygon::getEdgeProjection() const
{
    return _edgeProjection;
}

bool esg::Polygon::setVertices(const Vector3* va, unsigned n)
{
    if (_vertices) delete [] _vertices;
    if (_normals)  delete [] _normals;
    if (_uvCoords) delete [] _uvCoords;
    if (_proj)     delete [] _proj;

    _vertices  = NULL;
    _normals   = NULL;
    _uvCoords  = NULL;
    _nVert     = n;
    _proj      = NULL;

    if (_nVert > 2) {
        _vertices = new Vector3d [_nVert];
        for (unsigned i = 0; i < _nVert; i++) _vertices[i].set(va[i]);
        _normal.set(_get_facet_normal());
        _normalFixed = false;
        _fxyz = - _normal.dot(_get_vertex(0));
        _set_edge_projection();
    } else {
        cerr << "setVerices() error: Number of vertices is less then 3" << endl;
        return false;
    }

#ifndef WIN32
#warning "FIXME: Area of polygon"
#endif
    _area = 1;

    return true;
}

bool esg::Polygon::setNormals(const Vector3* na)
{
    if (_nVert < 3) {
        cerr << "setNormals() error: Number of vertices is less then 3" << endl;
        return false;
    }

    if (_normals)  {
        delete [] _normals;
        _normals   = NULL;
    }

    if (na) {
        _normals = new Vector3 [_nVert];
        for (unsigned i = 0; i < _nVert; i++) _normals[i].set(na[i]);
    }

    return true;
}

bool esg::Polygon::setUVCoords(const Vector2* uva)
{
    if (_nVert < 3) {
        cerr << "setUVCoords() error: Number of vertices is less then 3" << endl;
        return false;
    }

    if (_uvCoords)  {
        delete [] _uvCoords;
        _uvCoords = NULL;
    }

    if (uva) {
        _uvCoords = new Vector2 [_nVert];
        for (unsigned i = 0; i < _nVert; i++) _uvCoords[i].set(uva[i]);
    }

    return true;
}

void esg::Polygon::setFacetNormal(const Vector3& norm, bool fixedn)
{
    _normal.set(norm);
    _normalFixed = fixedn;
    _fxyz = - _normal.dot(_get_vertex(0));
    _set_edge_projection();
}


void esg::Polygon::__debug()
{
    //fprintf(stderr,"Polygon: normal: %f %f %f\n",
    //      _normal.x,_normal.y,_normal.z);
    fprintf(stderr,"Polygon vertices:\n");
    for (unsigned i = 0; i < _nVert; i++)
        fprintf(stderr,"%f %f %f\n",
                getVertex(i).x,getVertex(i).y,getVertex(i).z);
    fprintf(stderr,"Polygon's normal: %f %f %f\n",_normal.x,_normal.y,_normal.z);
}

---