OSGExtrusionGeometry.cpp

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
 *                                OpenSG                                     *
 *                                                                           *
 *                                                                           *
 *           Copyright (C) 2000,2001,2002 by the OpenSG Forum                *
 *                                                                           *
 *                            www.opensg.org                                 *
 *                                                                           *
 *   contact: dirk@opensg.org, gerrit.voss@vossg.org, jbehr@zgdv.de          *
 *                                                                           *
\*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*\
 *                                License                                    *
 *                                                                           *
 * This library is free software; you can redistribute it and/or modify it   *
 * under the terms of the GNU Library General Public License as published    *
 * by the Free Software Foundation, version 2.                               *
 *                                                                           *
 * This library is distributed in the hope that it will be useful, but       *
 * WITHOUT ANY WARRANTY; without even the implied warranty of                *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU         *
 * Library General Public License for more details.                          *
 *                                                                           *
 * You should have received a copy of the GNU Library General Public         *
 * License along with this library; if not, write to the Free Software       *
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 *
 *                                                                           *
\*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*\
 *                                Changes                                    *
 *                                                                           *
 *                                                                           *
 *                                                                           *
 *                                                                           *
 *                                                                           *
 *                                                                           *
\*---------------------------------------------------------------------------*/

//-----------------------------------------------------------------------------
//                                 Includes                                    
//-----------------------------------------------------------------------------


#include <memory>
#include <assert.h>

#include "OSGConfig.h"
#include <OSGLog.h>

#include "OSGExtrusionGeometry.h"

// #define DEBUG_VERTEX_NORMALS
// #define DEBUG_FACE_NORMALS


OSG_USING_NAMESPACE


/*---------------------------------------------------------------------*/
//-----------------------------------------------------------------------------
//                           Forward Declarations                              
//-----------------------------------------------------------------------------
template <class VectorTypeT>
void subdivide(const typename std::vector<VectorTypeT> &dataIn,
                     typename std::vector<VectorTypeT> *dataOut,
                     bool closed);

template <class VectorTypeT>
void subdivide(const typename std::vector<VectorTypeT> &dataIn,
                     typename std::vector<VectorTypeT> *dataOut,
                     UInt32 nTimes,
                     bool close);


//-----------------------------------------------------------------------------
// Constructor
//
// Author: afischle
//-----------------------------------------------------------------------------
ExtrusionSurface::ExtrusionSurface(const std::vector<Pnt2f> &crossSection,
                                   const std::vector<Quaternion> &orientation,
                                   const std::vector<Vec2f> &scale,
                                   const std::vector<Pnt3f> &spine,
                                   Real32 creaseAngle,
                                   bool beginCap,
                                   bool endCap,
                                   bool ccw,
                                   bool convex,
                                   bool buildNormal,
                                   bool buildTexCoord) :

    _creaseAngle(creaseAngle),
    _beginCap(beginCap),
    _endCap(endCap),
    _ccw(ccw),
    _convex(convex),
    _createNormals(buildNormal),
    _createTexCoords(buildTexCoord),
    _primitiveCount(0),
    _vertexCount(0),
    _totalVertexCount(0)
{
    // initialize cross section
    if(!crossSection.empty()) 
    {
        // copy argument
        _crossSection.assign(crossSection.begin(), crossSection.end());
    }
    else
    {
        // copy default
        FDEBUG(("OSGExtrusion: Empty cross section. Using default.\n"));
        for(UInt32 i = 0; i < DEF_N_CROSS_SECTION_POINTS; i++)
            _crossSection.push_back(Pnt2f(DEF_CROSS_SECTION[i]));
    }

    // initialize orientations
    if(!orientation.empty())
    {
        // copy argument
        _orientation.assign(orientation.begin(), orientation.end());
    }
    else
    {
        // copy default
        FDEBUG(("OSGExtrusion: Empty orientation. Using default.\n"));
        for(UInt32 i = 0; i < DEF_N_ORIENTATION_PARAMS; i++)
            _orientation.push_back(Quaternion(Vec3f(DEF_ORIENTATION[i]),
                                              DEF_ORIENTATION[i][3]));
    }
    
    // initialize scale parameters
    if(!scale.empty())
    {
        // copy argument
        _scale.assign(scale.begin(), scale.end());
    }
    else
    {
        // copy default
        FDEBUG(("OSGExtrusion: Empty scale parameter. Using default.\n"));
        for(UInt32 i = 0; i < DEF_N_SCALE_PARAMS; i++)
            _scale.push_back(Vec2f(DEF_SCALE[i]));
    }


    // initialize spine
    if(!spine.empty())
    {
        // copy argument
        _spine.assign(spine.begin(), spine.end());
    }
    else
    {
        // copy default
        FDEBUG(("OSGExtrusion: Empty spine. Using default.\n"));
        for(UInt32 i = 0; i < DEF_N_SPINE_POINTS; i++)
            _spine.push_back(Pnt3f(DEF_SPINE[i]));
    }
}


//-----------------------------------------------------------------------------
// Returns true if the member fields describe a valid extrusion surface and
// false otherwise.
// 
// Author: afischle
//-----------------------------------------------------------------------------
bool ExtrusionSurface::verifyInput(void)
{
    if(_crossSection.size() < 2)
    {
        FWARNING(("OSGExtrusion:checkInput: Cross section is degenerate\n"));
        return false;
    }

    if(_spine.size() < 2)
    {
        FWARNING(("OSGExtrusion:checkInput: spine is degenerate\n"));
        return false;
    }

    return true;
}


//-----------------------------------------------------------------------------
// Determines the topology of the sweep surface and sets the topology member
// variables accordingly.
//
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::determineTopology(void)
{
    // compute vector differences between start and end points
    Vec3f sDiff = _spine.back() - _spine.front();
    Vec2f cDiff = _crossSection.back() - _crossSection.front();

    // determine whether the curves are closed
    _spineClosed = (sDiff.length() < Eps) ? true : false;
    _crossSectionClosed = (cDiff.length() < Eps) ? true : false;
    
    FDEBUG(("OSGExtrusion:checkInput: Topology: (Spine: %s / Cross-section: %s)\n", 
            (_spineClosed ? "Closed" : "Open"),
            (_crossSectionClosed ? "Closed" : "Open")));
}


//-----------------------------------------------------------------------------
// Refines the cross-section n times using a simple subdivision scheme.
// 
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::refineCrossSection(UInt32 nTimes)
{
    FDEBUG(("OSGExtrusion:refineCrossSection: Refinining %d times\n",nTimes));
    std::vector<Vec2f> output;

    // cast <Pnt2f> to <Vec2f> because we need the vector interface
    std::vector<Vec2f> & csVec2f ((std::vector<Vec2f>&)(_crossSection));
    subdivide(csVec2f, &output, nTimes, _crossSectionClosed);

    // swap output data into _crossSection
    output.swap(csVec2f);
}

//-----------------------------------------------------------------------------
// Refines the spine n times using a simple subdivision scheme.
// 
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::refineSpine(UInt32 nTimes)
{
    FDEBUG(("OSGExtrusion:refineSpine: Refinining %d times\n",nTimes));

    std::vector<Vec3f> output;

    // cast <Pnt3f> to <Vec3f> because we need the vector interface    
    std::vector<Vec3f> & spVec3f ((std::vector<Vec3f>&)(_spine));
    subdivide<Vec3f>(spVec3f, &output, nTimes, _spineClosed);

    // swap output data into _spine
    output.swap(spVec3f);
}


//-----------------------------------------------------------------------------
// Refines the scale parameters n times using a simple subdivision scheme.
// 
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::refineScale(UInt32 nTimes)
{
    FDEBUG(("OSGExtrusion:refineScale: Refinining %d times\n",nTimes));

    std::vector<Vec2f> output;
    subdivide<Vec2f>(_scale, &output, nTimes, _spineClosed);

    // swap output data into _scale
    output.swap(_scale);
}


//-----------------------------------------------------------------------------
// Refines the scale parameters n times using a simple subdivision scheme.
//
// Todo: Not implemented yet
//
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::refineOrientation(UInt32 nTimes)
{
    FFATAL(("OSGExtrusion:refineOrientation: Not Implemented\n"));
}


//-----------------------------------------------------------------------------
// Calculates the y-axes of the base transform from the cross section plane
// to the spine-aligned cross section plane (SCP). The y- and z-axes are
// expected to be already stored in the _transform vector.
//
// The computed x-axes are stored in the _transform vector.
// 
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::calcXAxes(void)
{
    std::vector<Matrix>::iterator tIt = _transform.begin();

    // Either of the two cases below implies the cross product of
    // the y- and the z-axis vanish (and so do all x-axes)
    if(!_spineCollinear && !_revolutionSurface)
    {
        for(; tIt != _transform.end(); ++tIt)
        {
            Vec3f xAxis = calcXAxis(tIt);
            (*tIt)[0].setValue(xAxis);
        }
    }
    else
    {
        (*tIt)[0].setValue(Vec3f(0.f, 0.f, 0.f));
    }
}


//-----------------------------------------------------------------------------
// Calculates the y-axes of the base transform from the cross-section plane
// to the spine-aligned cross section plane (SCP). If no y-axis could be
// defined it sets the _revolutionSurface flag to true (and false otherwise.)
//
// The computed y-axes are stored in the _transform vector.
//
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::calcYAxes(void)
{
    bool yAxisDefined = false;

    std::vector<Pnt3f>::const_iterator  spineIt = _spine.begin();
    std::vector<Matrix>::iterator tIt = _transform.begin();
    for(; spineIt != _spine.end(); ++spineIt, ++tIt)
    {
        Vec3f yAxis = calcNonUnitYAxis(spineIt);
        if(yAxis.length() > Eps) // we have found a valid tangent
        {
            yAxis.normalize();
            if(!yAxisDefined)
            {
                FDEBUG(("OSGExtrusion:calcYAxes: yAxis is defined\n"));
                yAxisDefined = true;

                // assign found y-axis to all transforms
                // up to and including the current one
                std::vector<Matrix>::iterator yIt;
                for(yIt = _transform.begin(); yIt != (tIt + 1); ++yIt)
                    (*yIt)[1].setValue(yAxis);
            }
            else 
            {
                // assign found y-axis to the current transform only
                (*tIt)[1].setValue(yAxis);
            }
        }
        else if(yAxisDefined)
        {
            // reuse y-axis determined during last step
            (*tIt)[1] = (*(tIt - 1))[1];
        }
    }

    // all spine points coincided => surface of revolution
    _revolutionSurface = !yAxisDefined;
}

//-----------------------------------------------------------------------------
// Calculates the z-axes of the base transform from the cross-section plane (CP)
// to the spine-aligned cross section plane (SCP). If no z-axis could be
// defined it sets the _spineCollinear flag to true (and false otherwise.)
//
// The computed z-axes are stored in the _transform vector.
//
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::calcZAxes(void)
{
    bool zAxisDefined = false;
    
    std::vector<Pnt3f>::const_iterator  spineIt = _spine.begin();
    std::vector<Matrix>::iterator tIt = _transform.begin();
    for(; spineIt != _spine.end(); ++spineIt, ++tIt)
    {
        Vec3f zAxis = calcNonUnitZAxis(spineIt);
        if(zAxis.length() > Eps)
        {
            zAxis.normalize();
            if(!zAxisDefined)
            {
                FDEBUG(("OSGExtrusion:calcZAxes: zAxis is defined\n"));
                zAxisDefined = true;
                
                // set the current z-Axis for the current transform
                (*tIt)[2].setValue(zAxis);

                // assign found z-axis to all previous transforms
                std::vector<Matrix>::iterator zIt;
                for(zIt = _transform.begin(); zIt != tIt; ++zIt)
                {
                    (*zIt)[2].setValue(zAxis);
                }
            }
            else // assign found z-axis to this basis _only_
            {
                (*tIt)[2].setValue(zAxis);

                // as the zAxis is defined tIt - 1 contains a valid value 
                Vec3f prevZAxis((*(tIt - 1))[2][0],
                                (*(tIt - 1))[2][1],
                                (*(tIt - 1))[2][2]);

                // check whether the z-axis needs to be flipped
                if(prevZAxis.dot(zAxis) + Eps < 0.f)
                {
                    FDEBUG(("OSGExtrusion:calcZAxes: Flipping z-axis\n"));
                    zAxis.negate();
                }

                // assign zAxis to current transform
                (*tIt)[2].setValue(zAxis);
            }
        }
        else if(zAxisDefined)
        {
            // reuse previous value
            (*tIt)[2] = (*(tIt - 1))[2];
        }
    } // end for

    // All cross products vanished => spine has no non-collinear part
    // and is thus composed of collinear parts and revolutions only
    _spineCollinear = !zAxisDefined;
}

//-----------------------------------------------------------------------------
// Calculates the affine basis transformations from the cross-section plane
// (i.e. the xz-plane) to the spine-aligne cross-section plane (SCP) and
// determines whether:
//
// 1. the surface has a collinear spine,
// 2. the surface is a pure surface of revolution
// 3. the surface is not of one of the above special types.
//
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::calcTransforms(void)
{
    _transform.clear();

    // initialize the transformations with Id
    _transform.resize(_spine.size());

    /* X3D Spec:
     * The SCP is computed by first computing its Y-axis and Z-axis, then 
     * taking the cross product of these to determine the X-axis. These three 
     * axes are then used to determine the rotation value needed to rotate the
     * Y=0 plane to the SCP. This results in a plane that is the approximate 
     * tangent of the spine at each point, as shown in Figure 13.5.
     */
    calcZAxes(); // sets bool _spineCollinear
    calcYAxes(); // sets bool _revolutionSurface

    // We need to make a difference between surfaces having at least
    // one non-trivial collinear part in their spine and surfaces created
    // by pure revolution of the crossection.
    if(_revolutionSurface)
    {
        _spineCollinear = false;
    }
    calcXAxes(); // uses _spineCollinear and _revolutionSurface

    /* X3D Spec:
     *
     * If the entire spine is collinear, the SCP is computed by finding the 
     * rotation of a vector along the positive Y-axis (v1) to the vector formed
     * by the spine points (v2). The Y=0 plane is then rotated by this value.
     *
     * If two points are coincident, they both have the same SCP. If each point
     * has a different orientation value, then the surface is constructed by 
     * connecting edges of the cross-sections as normal. This is useful in 
     * creating revolved surfaces.
     */
    Matrix rotation; // m = Id and thus ok for revolution case
    if(_spineCollinear && !_revolutionSurface)
    {
        // calculate the rotation that rotates the y unit vector to the
        // direction of the spine.
        Vec3f v1(0.f, 1.f, 0.f);
        Vec3f v2(_transform[0][1]);
            
        Quaternion q;
        q.setValue(v1, v2);
        q.getValue(rotation);
    }

    std::vector<Matrix>::iterator tIt             = _transform.begin();
    std::vector<Pnt3f>::const_iterator spineIt    = _spine.begin();
    std::vector<Vec2f>::const_iterator scaleIt    = _scale.begin(); 
    std::vector<Quaternion>::const_iterator rotIt = _orientation.begin();

    // calculate an affine transformation for the cross-section at
    // each spine point
    for(; spineIt != _spine.end(); ++spineIt, ++tIt)
    {
        // special cases
        if(_spineCollinear || _revolutionSurface)
        {
            // defines the x,y and z-axes in this case
            *tIt = rotation;
        }        

        // setup transform of the CP
        Matrix cSPlaneTransform;

        // Scale the x- and z-axes
        cSPlaneTransform[0] *= (*scaleIt)[0];
        cSPlaneTransform[2] *= (*scaleIt)[1];

        // apply orientation
        Matrix orientationMatrix;
        rotIt->getValue(orientationMatrix);
        cSPlaneTransform.multLeft(orientationMatrix);

        tIt->mult(cSPlaneTransform);

        // affine translation to the spine point
        tIt->setTranslate(*spineIt);
        
        // The affine transformation *tIt should now
        //
        //   1. scale the CP
        //   2. apply the orientation to the CP        
        //   2. translate the origin to spine[i]
        //
        // in exactly this order.

        /* X3D Spec quotation:
         *
         * If the number of scale or orientation values is greater than the
         * number of spine points, the excess values are ignored. If they 
         * contain one value, it is applied at all spine points. The 
         * results are undefined if the number of scale or orientation 
         * values is greater than one but less than the number of spine 
         * points. The scale values shall be positive.
         */
        
        // always use last available value
        if(rotIt + 1   != _orientation.end()) { ++rotIt;   }
        if(scaleIt + 1 != _scale.end())       { ++scaleIt; }
    }

    FDEBUG(("OSGExtrusion:calcTransforms: spine collinear: %s\n", 
            (_spineCollinear ? "true" : "false")));

    FDEBUG(("OSGExtrusion:calcTransforms: revolution surface: %s\n",
            (_revolutionSurface ? "true" : "false")));
}


//-----------------------------------------------------------------------------
// Initializes the vertex grid with _spine.size() rows and _crossSection.size()
// columns.
//
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::initGrid(void)
{
    FDEBUG(("OSGExtrusion:initGrid: Initializing vertex grid\n"));
    
    _grid.clear();
    _grid.resize(_spine.size());

    VertexGrid::iterator gridRowIt = _grid.begin();
    for(Int32 i = 0; i < _spine.size(); i++, ++gridRowIt)
    {
        gridRowIt->resize(_crossSection.size());
    }
}


//-----------------------------------------------------------------------------
// Calculates the vertex positions of the sweep surface by transforming the 
// cross-section using the _transforms member field.
//
// The calculated positions are stored in the position components of the
// vertices in the _grid.
//
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::calcSweepSurfacePositions(void)
{
    FDEBUG(("OSGExtrusion:calcSweepSurfacePositions: Calculation vertex positions\n"));

    if(_grid.size() < _spine.size())
    {
        FFATAL(("OSGExtrusion:calcSweepSurfacePositions: Grid not initialized\n"));
        return;
    }

    _posMap.clear();

    // compute transformation matrices
    calcTransforms();

    if(_grid.size() < _spine.size()) // grid not yet or incorrectly initialized
    {
        initGrid();
    }

    std::vector<Matrix>::iterator transIt = _transform.begin();
    VertexGrid::iterator gridRowIt = _grid.begin();

    // iterate over the spine points
    for(Int32 i = 0; i < _spine.size(); i++, ++gridRowIt)
    {
        if(transIt->det3() < Eps) // degenerate transform
        {
            FWARNING(("OSGExtrusion:calcSurfacePositions: Degenerate transformation matrix\n"));
            SWARNING << "Degenerate transform[" << i << "]:\n "
                     << *transIt << std::endl;
        }

        std::vector<Vertex>::iterator gridColIt = gridRowIt->begin();
        std::vector<Pnt2f>::iterator crossIt;

        // iterate over the cross-section curve
        for(crossIt = _crossSection.begin(); crossIt != _crossSection.end();
            ++crossIt, ++gridColIt)
        {
            Pnt3f p = Pnt3f((*crossIt)[0], 0.f, (*crossIt)[1]);

            // apply "CP to SCP"-transform
            transIt->multMatrixPnt(p, gridColIt->position);
        }

        // always use last defined transformation matrix
        if(transIt + 1 != _transform.end()) { ++transIt; }
    }
}


//-----------------------------------------------------------------------------
// Calculates the texture coordinates of the sweep surface. The texture coords
// are obtained by measuring the lengths along the row- and column-curves
// of the grid and renormalizing them to the usual [0, 1]-range. The grid
// is expected to already hold the position data of the sweep surface.
//
// The calculated texture coordinates are stored in the _grid.
//
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::calcSweepSurfaceTexCoords(void)
{
    FDEBUG(("OSGExtrusion:calcSweepSurfaceTexCoords: Calculating texture coordinates\n"));
    if(_grid.size() < _spine.size())
    {
        FFATAL(("OSGExtrusion:calcSweepSurfaceTexCoords: Grid not initialized\n"));
        return;
    }
    _texCoordMap.clear();
    
    /* X3D Spec:
     *
     *  Textures are mapped so that the coordinates range in the U direction 
     *  from 0 to 1 along the crossSection curve (with 0 corresponding to the 
     *  first point in crossSection and 1 to the last) and in the V direction
     *  from 0 to 1 along the spine curve (with 0 corresponding to the first 
     *  listed spine point and 1 to the last) 
     */

    // accumulate distances along u and v directions
    for(UInt32 i = 0; i < _spine.size(); ++i)
    {
        for(UInt32 j = 0; j < _crossSection.size(); ++j)
        {
            Real32 u,v;
            UInt32 iPrev = (i == 0) ? 0 : i - 1;
            UInt32 jPrev = (j == 0) ? 0 : j - 1;

            // contains length of i-th column curve up to the (j-1)-th point
            u = _grid[i][jPrev].texCoord[0];
            // add distance to j-th point
            u += _grid[i][j].position.dist(_grid[i][jPrev].position);

            // contains length of j-th row curve up to the (i-1)-th point
            v = _grid[iPrev][j].texCoord[1];
            // add distance to i-th point
            v += _grid[i][j].position.dist(_grid[iPrev][j].position);

            _grid[i][j].texCoord[0] = u;
            _grid[i][j].texCoord[1] = v;
        }
    }

    // renormalize distances to [0,1]x[0,1]
    for(UInt32 i = 0; i < _spine.size(); ++i)
    {
        for(UInt32 j = 0; j < _crossSection.size(); ++j)
        {
            // length is maximal in the row/column curve end points
            Real32 uMax = _grid[i][_crossSection.size() - 1].texCoord[0];
            Real32 vMax = _grid[_spine.size() - 1][j].texCoord[1];

            if(uMax > Eps)
                _grid[i][j].texCoord[0] /= uMax;

            if(vMax > Eps)
                _grid[i][j].texCoord[1] /= vMax;
        }
    }
}


void ExtrusionSurface::calcSweepSurfaceFaceNormals(void)
{
    FDEBUG(("OSGExtrusion:calcSweepSurfaceFaceNormals: Calculating sweep surface face normals\n"));

    if(_grid.size() < _spine.size())
    {
        FFATAL(("OSGExtrusion:calcSweepSurfaceFaceNormals: Grid not initialized\n"));
        return;
    }
    
    // calculate face normals
    for(UInt32 i = 0; i < _spine.size(); ++i)
    {
        UInt32 iPrev, iNext;

        // determine adjacent row indices
        if(_spineClosed)
        {
            // first point == last point
            iPrev = (i != 0) ? i - 1 : _spine.size() - 2;
            iNext = (i != _spine.size() - 1) ? i + 1 : 1;
        }
        else
        {
            iPrev = (i != 0) ? i - 1 : 0;
            iNext = (i != _spine.size() - 1) ? i + 1 : i;
        }

        for(UInt32 j = 0; j < _crossSection.size(); ++j)
        {
            Vec3f n;
            UInt32 jPrev, jNext;

            // determine adjacent column indices
            if(_crossSectionClosed)
            {
                // first point == last point
                jPrev = (j != 0) ? j - 1 : _crossSection.size() - 2;
                jNext = (j != _crossSection.size() - 1) ? j + 1 : 1;
            }
            else
            {
                jPrev = (j != 0) ? j - 1 : 0;
                jNext = (j != _crossSection.size() - 1) ? j + 1 : j;
            }

            // upper right face normal (1st quadrant)
            _grid[i][j].adjFaceNormals[0] =
                calcQuadFaceNormal(_grid[i][j].position,
                                   _grid[i][jNext].position,
                                   _grid[iNext][jNext].position,
                                   _grid[iNext][j].position);


            // upper left face normal (2nd quadrant)
            _grid[i][j].adjFaceNormals[1] = 
                calcQuadFaceNormal(_grid[i][j].position,
                                   _grid[iNext][j].position,
                                   _grid[iNext][jPrev].position,
                                   _grid[i][jPrev].position);


            // lower left face normal (3nd quadrant)
            _grid[i][j].adjFaceNormals[2] = 
                calcQuadFaceNormal(_grid[i][j].position,
                                   _grid[i][jPrev].position,
                                   _grid[iPrev][jPrev].position,
                                   _grid[iPrev][j].position);


            // lower right face normal (4th quadrant)
            _grid[i][j].adjFaceNormals[3] = 
                calcQuadFaceNormal(_grid[i][j].position,
                                   _grid[iPrev][j].position,
                                   _grid[iPrev][jNext].position,
                                   _grid[i][jNext].position);
        }
    }
}



//-----------------------------------------------------------------------------
// 
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::storeSweepSurfaceWithoutNormals(GeoIndicesUI32Ptr  indicesPtr,
                                                       GeoPLengthsUI32Ptr lensPtr,
                                                       GeoPTypesUI8Ptr    typesPtr)
{
    assert(indicesPtr != NullFC);
    assert(lensPtr    != NullFC);
    assert(typesPtr   != NullFC);

    for(UInt32 i = 0; i < _spine.size() - 1; i++)
    {
        for(UInt32 j = 0; j < _crossSection.size(); j++)
        {
            Vertex v;
            // begin new triangle strip
            v = _grid[i+1][j];

            // store position (and texcoord if available)
            storeVertex(v, indicesPtr);

            v = _grid[i][j];
            // store position (and texcoord if available)
            storeVertex(v, indicesPtr);
        }
        // store row strip
        storePrimitive(GL_TRIANGLE_STRIP, lensPtr, typesPtr);
    }
}

void ExtrusionSurface::storeSweepSurfaceWithNormals(GeoIndicesUI32Ptr  indicesPtr,
                                                    GeoPLengthsUI32Ptr lensPtr,
                                                    GeoPTypesUI8Ptr    typesPtr)
{
    assert(indicesPtr != NullFC);
    assert(lensPtr    != NullFC);
    assert(typesPtr   != NullFC);

    // we need to calculate the face normals first
    calcSweepSurfaceFaceNormals();
 
    if(_ccw)
    {
        // The current face is the face in quadrant 1 of the vertex
        // grid[i,j]. Stripes are generated in positive j direction
        // along the cross-section.
        UInt32 spineEnd = _spine.size() - 1;
        UInt32 cSEnd = _crossSection.size() - 1;
        for(UInt32 i = 0; i < spineEnd; i++)
        {
            // upper left vertex of current quad face
            // current face is located in 4th quadrant of vertex
            calcVertexNormal(&_grid[i+1][0], 3);
            storeVertex(_grid[i+1][0], indicesPtr);

            // lower left vertex of current quad face
            // current face is located in 1st quadrant of vertex            
            calcVertexNormal(&_grid[i][0], 0);
            storeVertex(_grid[i][0], indicesPtr);

            Vec3f fN = _grid[i][0].adjFaceNormals[0];
            for(UInt32 j = 0; j < cSEnd; j++)
            {
                // The next face is the face in quadrant 1 of the vertex
                // grid[i,j+1]                
                Vec3f fNNext = _grid[i][j+1].adjFaceNormals[0];

                // if the creaseAngle is too small the strip has to be split,
                // because the normals are not shared along the common edge
                // of the current and the next face.
                // if j + 1 == csEnd the current quad is the last one
                if((fN.enclosedAngle(fNNext) + Eps > _creaseAngle - Eps)
                   || (j + 1 == cSEnd))
                {
                    // upper right vertex of current quad face
                    // current face is located in 3rd quadrant of vertex
                    calcVertexNormal(&_grid[i+1][j+1], 2);
                    storeVertex(_grid[i+1][j+1], indicesPtr);

                    // upper right vertex of current quad face
                    // current face is located in 2nd quadrant of vertex
                    calcVertexNormal(&_grid[i][j+1], 1);
                    storeVertex(_grid[i][j+1], indicesPtr);
                    
                    // terminate strip
                    storePrimitive(GL_TRIANGLE_STRIP, lensPtr, typesPtr);
                }

                // Either the strip was split above or not. In both cases
                // continue/begin new strip if there's still a quad left
                if(j + 1 < cSEnd)
                {
                    
                    // upper left vertex of _next_ quad face
                    // next face is located in 4th quadrant of vertex
                    calcVertexNormal(&_grid[i+1][j+1], 3); 
                    storeVertex(_grid[i+1][j+1], indicesPtr);
                    
                    // lower left vertex of _next_ quad face
                    // next face is located in 1st quadrant of vertex
                    calcVertexNormal(&_grid[i][j+1], 0);
                    storeVertex(_grid[i][j+1], indicesPtr);
                }

                // update face normal
                fN = fNNext;                
            }
        }
    }
    else // cw case
    {
        // The current face is the face in quadrant 1 of the vertex
        // grid[i,j-1]. Stripes are generated in negative j direction
        // along the cross-section.
        for(UInt32 i = 0; i < _spine.size() - 1; i++)
        {
            UInt32 cSEnd = _crossSection.size() - 1;

            // begin new triangle strip
            calcVertexNormal(&_grid[i+1][cSEnd], 2);
            storeVertex(_grid[i+1][cSEnd], indicesPtr);

            calcVertexNormal(&_grid[i][cSEnd], 1);
            storeVertex(_grid[i][cSEnd], indicesPtr);

            Vec3f fN = _grid[i][cSEnd].adjFaceNormals[0];            
            for(UInt32 j = cSEnd; j > 0; j--)
            {
                Vec3f fNNext = _grid[i][j - 1].adjFaceNormals[0];
                if((fN.enclosedAngle(fNNext) + Eps > _creaseAngle - Eps)
                   || (j == 1))
                {
                    calcVertexNormal(&_grid[i + 1][j - 1], 3);
                    storeVertex(_grid[i + 1][j - 1], indicesPtr);
                    
                    calcVertexNormal(&_grid[i][j - 1], 0);
                    storeVertex(_grid[i][j - 1], indicesPtr);
                    
                    // terminate strip
                    storePrimitive(GL_TRIANGLE_STRIP, lensPtr, typesPtr);
                }

                if(j > 1) // still one quad to go
                {
                    calcVertexNormal(&_grid[i + 1][j - 1], 2);
                    storeVertex(_grid[i + 1][j - 1], indicesPtr);
                
                    calcVertexNormal(&_grid[i][j - 1], 1);
                    storeVertex(_grid[i][j - 1], indicesPtr);
                }
                // update face normal
                fN = fNNext;                
            }
        }
    }
}


//----------------------------------------------------------------------
// Returns true if the origin lies in the left half plane of the
// oriented line AB.
//
// Author: afischle
//----------------------------------------------------------------------
inline bool ExtrusionSurface::isLeft(const Pnt2f &a, const Pnt2f &b)
{
    Vec2f d = b - a;
    d.normalize();
    // compute the normal pointing into the _right_ half plane of the oriented
    // line
    Vec2f n = Vec2f(d.y(), -d.x());
    // the dot product is >0 when two vectors point into the same direction
    return n.dot(a) > 0.f;
}


//-----------------------------------------------------------------------------
// Computes the minimum of the y-components of all points of the contour.
//
// Author: afischle
//-----------------------------------------------------------------------------
Real32 ExtrusionSurface::computeMinYAbs(const std::vector<Pnt2f> &contour,
                                        Real32 alpha)
{
    const Real32 sinAlpha = osgsin(alpha);
    const Real32 cosAlpha = osgcos(alpha);

    Real32 minYAbs = 1.f;
    for(std::vector<Pnt2f>::const_iterator it = contour.begin();
        it != contour.end(); ++it)
    {
        Real32 y = sinAlpha * it->x() + cosAlpha * it->y();
        if(osgabs(y) < minYAbs)
            minYAbs = y;
    }

    return minYAbs;
}

//-----------------------------------------------------------------------------
// Computes a rotation of the curve such that the minimum of the y-components
// of all of its points is the maximal for alpha = n * (2 * Pi / nAngles),
// where 0 <= n <= nAngles.
//
// Author: afischle
//-----------------------------------------------------------------------------
Real32 ExtrusionSurface::calcBetterRotationAngle(const std::vector<Pnt2f> &contour,
                                                 UInt32 nAngles)
{
    Real32 maxMinY = 0.f;
    Real32 bestAlpha = 0.f;

    Real32 alphaInc = 2 * Pi / nAngles;
    for(Real32 alpha = 0.f; alpha < 2 * Pi; alpha += alphaInc)
    {
        Real32 minYAlpha = computeMinYAbs(contour, alpha);
        if(minYAlpha > maxMinY)
        {
            maxMinY = minYAlpha;
            bestAlpha = alpha;
        }
    }

    return bestAlpha;
}


//-----------------------------------------------------------------------------
// Translates the given curve such that the specified reference point is
// translated to the origin and rescales it to the unit square.
//
// Author: afischle
//-----------------------------------------------------------------------------
bool ExtrusionSurface::calcOptimizedContour(const std::vector<Pnt2f> &contour,
                                            const Pnt2f &point,
                                            std::vector<Pnt2f> *optimizedContour)
{
    // copy the input
    std::vector<Pnt2f> coords(contour);

    // translate the reference point to the origin
    // and compute a scale factor
    Real32 maxLen = std::numeric_limits<Real32>::min();
    Real32 minLen = std::numeric_limits<Real32>::max();

    for(std::vector<Pnt2f>::iterator it = coords.begin();
            it != coords.end(); ++it)
    {
        // translate the problem to the origin
        *it -= Vec2f(point);
        Real32 len = Vec2f(*it).length();

        if(len > maxLen)  maxLen = len;
        if(len < minLen)  minLen = len;
    }

    if(maxLen < Eps)
    {
        FWARNING(("OSGExtrusion:calcOptimizedContour: Contour completely degenerate\n"));
        return false;
    }

    if(minLen < Eps)
    {
        FWARNING(("OSGExtrusion:calcOptimizedContour: Reference point on curve\n"));
        return false;
    }

    // rescale the contour to the unit disk and determine the minimal
    // absolute y coordinate
    Real32 minAbsY = 1.f;
    for(std::vector<Pnt2f>::iterator it = coords.begin();
            it != coords.end(); ++it)
    {
        *it /= maxLen;
        // osgabs(it->y()) <= 1.f after rescale
        if(osgabs(it->y()) < minAbsY)
            minAbsY = osgabs(it->y());
    }

    optimizedContour->clear();
    optimizedContour->reserve(coords.size());
    
    if(minAbsY > Eps) // we can probably determine all crossings
    {
        optimizedContour->assign(coords.begin(), coords.end());
    }
    else // try to rotate the contour
    {
        // we need to rotate the contour to get a reasonable result
        FDEBUG(("OSGExtrusion:calcOptimizedContour: Contour needs rotation\n"));
        // compute the best rotation angle for 64 equal partitions 
        // of the circle
        Real32 alpha = calcBetterRotationAngle(contour, 64);
        FDEBUG(("OSGExtrusion:calcOptimizedContour: Rotation angle is: %f\n",
                alpha));

        Real32 sinAlpha = osgsin(alpha);
        Real32 cosAlpha = osgcos(alpha);

        // rotate the contour and copy it into the output vector
        for(std::vector<Pnt2f>::iterator it = coords.begin();
            it != coords.end(); ++it)
        {
            Pnt2f p;
            p[0] = cosAlpha * it->x() - sinAlpha * it->y();
            p[1] = sinAlpha * it->x() + cosAlpha * it->y();

            optimizedContour->push_back(p);
        }
    }

    // optimization successful
    return true;
}


//----------------------------------------------------------------------
// Returns the winding number of the given contour with respect to the
// specified point by counting how often the infinite line f(x) = point.y
// was crossed by the curve.
// 
// Author: afischle
//----------------------------------------------------------------------
Int32 ExtrusionSurface::calcWindingNumber(const std::vector<Pnt2f> &contour, 
                                          const Pnt2f &point)
{
    Int32 windingNumber = 0;
    std::vector<Pnt2f> coords;

    // translate contour to origin and rotate it to a better position
    calcOptimizedContour(contour, point, &coords);

    const Pnt2f *prevPoint = &(coords[coords.size() - 1]);
    for(UInt32 i = 0; i < coords.size(); ++i)
    {
        const Pnt2f *curPoint = &(coords[i]);
        if(prevPoint->y() < 0.f)
        {
            // upward crossing
            if(curPoint->y() > 0.f)
            {
                if(isLeft(*prevPoint, *curPoint) == true)
                    windingNumber++;
            }
        }
        else // prevPoint->y() > 0.f
        {
            // downward crossing
            if(curPoint->y() < 0.f)
            {
                if(isLeft(*prevPoint, *curPoint) == false)
                    windingNumber--;
            }
        }
        prevPoint = curPoint;
    }

    return windingNumber;
}


//-----------------------------------------------------------------------------
// Calculates a convex cap surface bounded by the row-curve with index
// spineIndex in the grid. If it is a begin cap its normal is reversed.
//
// The positions, normals and texture coordinates of the generated primitives
// are inserted into maps which map the respective properties to an integer
// index. The indices so obtained are stored in the indices field container.
// The lengths and types of the generated primitives are directly stored in
// the respective field containers.
//
// Author: afischle
//-----------------------------------------------------------------------------
void ExtrusionSurface::storeConvexCap(UInt32 spineIndex, 
                                      bool isBeginCap, 
                                      GeoIndicesUI32Ptr  indicesPtr,
                                      GeoPLengthsUI32Ptr lensPtr,
                                      GeoPTypesUI8Ptr    typesPtr)
{
    assert(indicesPtr != NullFC);
    assert(lensPtr    != NullFC);
    assert(typesPtr   != NullFC);

    if(!_convex)
    {