///////////////////////////////////////////////////////////////////////////////
// TrefoilKnot.cpp
// ===============
// N-leaf Trefoil Knot geometry for OpenGL with (majorR, minorR, tubeRadius, n, sectors, sides)
// x = r * sin(a) - R * sin((n-1)*a)
// y = r * cos(a) + R * cos((n-1)*a)
// z = r * sin(n*a)
// where n = # of leaves
// 
// The minimum # of sectors is 18 and # of sides are 2.
// The minimum # of n is 2, and it becomes 2-leaf clover shape knot.
// - major radius(R): distance from the origin to the centre of the torus
// - minor radius(r): radius of the torus tube
// - tube radius: radius of the knot tube
// - n: # of leaves of the knot
// - sectors: # of sectors of the tube
// - sides: # of sides of the tube
// - smooth: smooth (default) or flat shading
// - up-axis: facing direction, X=1, Y=2, Z=3(default)
// - mode: tube generation mode, 0=projection, 1=transform
//
//  AUTHOR: Song Ho Ahn (song.ahn@gmail.com)
// CREATED: 2026-02-03
// UPDATED: 2026-02-03
///////////////////////////////////////////////////////////////////////////////

#ifdef _WIN32
#include <windows.h>    // include windows.h to avoid thousands of compile errors even though this class is not depending on Windows
#endif

#ifdef __APPLE__
#include <OpenGL/gl.h>
#else
#include <GL/gl.h>
#endif

#include <iostream>
#include <iomanip>
#include <cmath>
#include "TrefoilKnot.h"



// constants //////////////////////////////////////////////////////////////////
const int MIN_SECTOR_COUNT = 18;
const int MIN_SIDE_COUNT  = 2;
const int MIN_N_COUNT = 2;
const float PI = acos(-1.0f);
const float RAD2DEG = 180.0f / PI;
const float EPSILON = 0.017f;        // less than 1 degree, 1/180*PI = 0.0174



///////////////////////////////////////////////////////////////////////////////
// ctor
///////////////////////////////////////////////////////////////////////////////
TrefoilKnot::TrefoilKnot(float majorR, float minorR, float tubeRadius, int n, int sectors, int sides, bool smooth, int up, int mode) : interleavedStride(32)
{
    set(majorR, minorR, tubeRadius, n, sectors, sides, smooth, up, mode);
}



///////////////////////////////////////////////////////////////////////////////
// setters
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::set(float majorR, float minorR, float tubeR, int n, int sectors, int sides, bool smooth, int up, int mode)
{
    if(majorR > 0)
        this->majorRadius = majorR;
    if(minorR > 0)
        this->minorRadius = minorR;
    if(tubeR > 0)
        this->tubeRadius = tubeR;
    this->n = n;
    if(n < MIN_N_COUNT)
        this->n = MIN_N_COUNT;
    this->sectorCount = sectors;
    if(sectors < MIN_SECTOR_COUNT)
        this->sectorCount = MIN_SECTOR_COUNT;
    this->sideCount = sides;
    if(sides < MIN_SIDE_COUNT)
        this->sideCount = MIN_SIDE_COUNT;
    this->smooth = smooth;
    this->upAxis = up;
    if(up < 1 || up > 3)
        this->upAxis = 3;
    this->mode = mode;
    if(mode < 0 || mode > 1)
        this->mode = 0;

    if(smooth)
        buildVerticesSmooth();
    else
        buildVerticesFlat();
}

void TrefoilKnot::setMajorRadius(float majorRadius)
{
    if(majorRadius != this->majorRadius)
        set(majorRadius, minorRadius, tubeRadius, n, sectorCount, sideCount, smooth, upAxis, mode);
}

void TrefoilKnot::setMinorRadius(float minorRadius)
{
    if(minorRadius != this->minorRadius)
        set(majorRadius, minorRadius, tubeRadius, n, sectorCount, sideCount, smooth, upAxis, mode);
}

void TrefoilKnot::setTubeRadius(float tubeRadius)
{
    if(tubeRadius != this->tubeRadius)
        set(majorRadius, minorRadius, tubeRadius, n, sectorCount, sideCount, smooth, upAxis, mode);
}

void TrefoilKnot::setN(int n)
{
    if(n != this->n)
        set(majorRadius, minorRadius, tubeRadius, n, sectorCount, sideCount, smooth, upAxis, mode);
}

void TrefoilKnot::setSectorCount(int sectors)
{
    if(sectors != this->sectorCount)
        set(majorRadius, minorRadius, tubeRadius, n, sectors, sideCount, smooth, upAxis, mode);
}

void TrefoilKnot::setSideCount(int sides)
{
    if(sides != this->sideCount)
        set(majorRadius, minorRadius, tubeRadius, n, sectorCount, sides, smooth, upAxis, mode);
}

void TrefoilKnot::setSmooth(bool smooth)
{
    if(this->smooth == smooth)
        return;

    this->smooth = smooth;
    if(smooth)
        buildVerticesSmooth();
    else
        buildVerticesFlat();
}

void TrefoilKnot::setUpAxis(int up)
{
    if(this->upAxis == up || up < 1 || up > 3)
        return;

    changeUpAxis(this->upAxis, up);
    this->upAxis = up;
}

void TrefoilKnot::setMode(int mode)
{
    if(this->mode == mode || mode < 0 || mode > 1)
        return;

    // regenerate tube
    if(smooth)
        buildVerticesSmooth();
    else
        buildVerticesFlat();
}



///////////////////////////////////////////////////////////////////////////////
// flip the face normals to opposite directions
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::reverseNormals()
{
    std::size_t i, j;
    std::size_t count = normals.size();
    for(i = 0, j = 3; i < count; i+=3, j+=8)
    {
        normals[i]   *= -1;
        normals[i+1] *= -1;
        normals[i+2] *= -1;

        // update interleaved array
        interleavedVertices[j]   = normals[i];
        interleavedVertices[j+1] = normals[i+1];
        interleavedVertices[j+2] = normals[i+2];
    }

    // also reverse triangle windings
    unsigned int tmp;
    count = indices.size();
    for(i = 0; i < count; i+=3)
    {
        tmp = indices[i];
        indices[i]   = indices[i+2];
        indices[i+2] = tmp;
    }
}



///////////////////////////////////////////////////////////////////////////////
// print itself
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::printSelf() const
{
    std::cout << "===== TrefoilKnot =====\n"
              << "  Major Radius: " << majorRadius << "\n"
              << "  Minor Radius: " << minorRadius << "\n"
              << "   Tube Radius: " << tubeRadius << "\n"
              << "  Leaves Count: " << n << "\n"
              << "  Sector Count: " << sectorCount << "\n"
              << "    Side Count: " << sideCount << "\n"
              << "Smooth Shading: " << (smooth ? "true" : "false") << "\n"
              << "       Up Axis: " << (upAxis == 1 ? "X" : (upAxis == 2 ? "Y" : "Z")) << "\n"
              << "     Tube Mode: " << (mode == 0 ? "Project" : "Transform") << "\n"
              << "Triangle Count: " << getTriangleCount() << "\n"
              << "   Index Count: " << getIndexCount() << "\n"
              << "  Vertex Count: " << getVertexCount() << "\n"
              << "  Normal Count: " << getNormalCount() << "\n"
              << "TexCoord Count: " << getTexCoordCount() << std::endl;
}



///////////////////////////////////////////////////////////////////////////////
// draw a trefoil knot in VertexArray mode
// OpenGL RC must be set before calling it
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::draw() const
{
    // interleaved array
    glEnableClientState(GL_VERTEX_ARRAY);
    glEnableClientState(GL_NORMAL_ARRAY);
    glEnableClientState(GL_TEXTURE_COORD_ARRAY);
    glVertexPointer(3, GL_FLOAT, interleavedStride, &interleavedVertices[0]);
    glNormalPointer(GL_FLOAT, interleavedStride, &interleavedVertices[3]);
    glTexCoordPointer(2, GL_FLOAT, interleavedStride, &interleavedVertices[6]);

    glDrawElements(GL_TRIANGLES, (unsigned int)indices.size(), GL_UNSIGNED_INT, indices.data());

    glDisableClientState(GL_VERTEX_ARRAY);
    glDisableClientState(GL_NORMAL_ARRAY);
    glDisableClientState(GL_TEXTURE_COORD_ARRAY);
}



///////////////////////////////////////////////////////////////////////////////
// draw lines only
// the caller must set the line width before call this
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::drawLines(const float lineColor[4]) const
{
    // set line colour
    glColor4fv(lineColor);
    glMaterialfv(GL_FRONT, GL_DIFFUSE, lineColor);

    // draw lines with VA
    glDisable(GL_LIGHTING);
    glDisable(GL_TEXTURE_2D);
    glEnableClientState(GL_VERTEX_ARRAY);
    glVertexPointer(3, GL_FLOAT, 0, vertices.data());

    glDrawElements(GL_LINES, (unsigned int)lineIndices.size(), GL_UNSIGNED_INT, lineIndices.data());

    glDisableClientState(GL_VERTEX_ARRAY);
    glEnable(GL_LIGHTING);
    glEnable(GL_TEXTURE_2D);
}



///////////////////////////////////////////////////////////////////////////////
// draw a trefoil knot surfaces and lines on top of it
// the caller must set the line width before call this
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::drawWithLines(const float lineColor[4]) const
{
    glEnable(GL_POLYGON_OFFSET_FILL);
    glPolygonOffset(1.0, 1.0f); // move polygon backward
    this->draw();
    glDisable(GL_POLYGON_OFFSET_FILL);

    // draw lines with VA
    drawLines(lineColor);
}



///////////////////////////////////////////////////////////////////////////////
// draw only path of trefoil knot
// the caller must set the line width before call this
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::drawPath(const float lineColor[4]) const
{
    // set line colour
    glColor4fv(lineColor);
    glMaterialfv(GL_FRONT, GL_DIFFUSE, lineColor);

    // draw lines with VA
    glDisable(GL_LIGHTING);
    glDisable(GL_TEXTURE_2D);
    glEnableClientState(GL_VERTEX_ARRAY);
    glVertexPointer(3, GL_FLOAT, 0, path.data());

    glDrawElements(GL_LINES, (unsigned int)pathIndices.size(), GL_UNSIGNED_INT, pathIndices.data());

    glDisableClientState(GL_VERTEX_ARRAY);
    glEnable(GL_LIGHTING);
    glEnable(GL_TEXTURE_2D);
}



///////////////////////////////////////////////////////////////////////////////
// draw contoures' first points for debug purpose
// the caller must set the point size before call this
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::drawFirstPoints(const float pointColor[4]) const
{
    // set point colour
    glColor4fv(pointColor);
    glMaterialfv(GL_FRONT, GL_DIFFUSE, pointColor);

    // draw points
    glDisable(GL_LIGHTING);
    glDisable(GL_TEXTURE_2D);
    glEnableClientState(GL_VERTEX_ARRAY);
    glVertexPointer(3, GL_FLOAT, 0, vertices.data());

    int step = sideCount * 6;
    int i;
    for(i = 0; i < (int)indices.size(); i += step)
    {
        glDrawArrays(GL_POINTS, indices[i], 1);
    }
    // last contour
    glDrawArrays(GL_POINTS, indices[i-1], 1);

    glDisableClientState(GL_VERTEX_ARRAY);
    glEnable(GL_LIGHTING);
    glEnable(GL_TEXTURE_2D);
}




///////////////////////////////////////////////////////////////////////////////
// dealloc vectors
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::clearArrays()
{
    std::vector<float>().swap(vertices);
    std::vector<float>().swap(normals);
    std::vector<float>().swap(texCoords);
    std::vector<unsigned int>().swap(indices);
    std::vector<unsigned int>().swap(lineIndices);
    std::vector<vec3>().swap(path);
    std::vector<vec3>().swap(contour);
    std::vector<vec3>().swap(contourNormal);
    std::vector<unsigned int>().swap(pathIndices);
    std::vector<vec3>().swap(pathDirections);
    std::vector<vec3>().swap(firstVectors);
}



///////////////////////////////////////////////////////////////////////////////
// generate path of N-leaf trefoil knot with (sectorCount+1) points
// x = r * sin(a) - R * sin((N-1)*a)
// y = r * cos(a) + R * cos((N-1)*a)
// z = r * sin(N*a)
// where 0 <= a <= 360
//       N = # of leaves
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::buildPath()
{
    TrefoilKnot::vec3 v;
    float angle;
    float sectorStep = 2 * PI / sectorCount;

    for(int i = 0; i <= sectorCount; ++i)
    {
        angle = i * sectorStep;                 // starting from 0 to 2pi
        v.x = minorRadius * sinf(angle) - majorRadius * sinf((n-1)*angle);
        v.y = minorRadius * cosf(angle) + majorRadius * cosf((n-1)*angle);
        v.z = minorRadius * sinf(n*angle);
        path.push_back(v);

        // build indices array
        if(i != sectorCount)
        {
            pathIndices.push_back(i);
            pathIndices.push_back(i+1);
        }
    }

    // make sure the first and last are same
    path[path.size()-1] = path[0];

    // remember direction vectors of the contour plane
    TrefoilKnot::vec3 dir1, dir2, dir;
    for(int i = 0; i <= sectorCount; ++i)
    {
        if(i == 0)
        {
            dir1 = path[path.size()-1] - path[path.size()-2];
            dir2 = path[1] - path[0];
        }
        else if(i == sectorCount)
        {
            dir1 = path[i] - path[i-1];
            dir2 = path[1] - path[0];
        }
        else
        {
            dir1 = path[i] - path[i-1];
            dir2 = path[i+1] - path[i];
        }
        dir = (dir1 + dir2).normalize();
        pathDirections.push_back(dir);
    }
}



///////////////////////////////////////////////////////////////////////////////
// generate circle on XY plane with tube radius
// x = r * cos(v)
// y = r * sin(v)
// z = 0
// where 0 <= v <= 360
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::buildContour()
{
    vec3 v;
    float angle;
    float sideStep = 2 * PI / sideCount;

    v.z = 0;
    for(int i = 0; i <= sideCount; ++i)
    {
        angle = i * sideStep;           // starting from pi to -pi
        v.x = tubeRadius * cosf(angle);
        v.y = tubeRadius * sinf(angle);

        contour.push_back(v);
        contourNormal.push_back(v.normalize());
    }
}



///////////////////////////////////////////////////////////////////////////////
// transform the contour at the first path point
// it will replace the first contour vertices
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::transformFirstContour()
{
    vec3 t = path[0];                           // translation
    vec3 f = pathDirections[0];                 // forward vector
    vec3 u(0, 0, 1);                            // tmp up vector
    vec3 l = u.cross(f).normalize();            // left vector
    u = f.cross(l);                             // recompute up vector

    // transform contour to the target(to) path point
    vec3 c, v;
    for(int i = 0; i <= sideCount; ++i)
    {
        c = contour[i];
        v.x = l.x*c.x + u.x*c.y + f.x*c.z + t.x;
        v.y = l.y*c.x + u.y*c.y + f.y*c.z + t.y;
        v.z = l.z*c.x + u.z*c.y + f.z*c.z + t.z;
        contour[i] = v;
        // update normal
        contourNormal[i] = (v - t).normalize();
    }
}



///////////////////////////////////////////////////////////////////////////////
// transform contour from the from-point to the to-point of the knot path
///////////////////////////////////////////////////////////////////////////////
std::vector<TrefoilKnot::vec3> TrefoilKnot::transformContour(int fromIndex, int toIndex)
{
    vec3 c = contour[0];
    vec3 t = path[toIndex];                     // translation
    vec3 f = pathDirections[toIndex];           // forward vector
    vec3 u(0,0,1);                              // tmp up vector
    vec3 l = u.cross(f).normalize();            // left vector
    u = f.cross(l);                             // recompute up vector

    // transform contour to the target(to) path point
    vec3 v;
    std::vector<vec3> toContour;
    for(int i = 0; i <= sideCount; ++i)
    {
        c = contour[i];
        v.x = l.x*c.x + u.x*c.y + f.x*c.z + t.x;
        v.y = l.y*c.x + u.y*c.y + f.y*c.z + t.y;
        v.z = l.z*c.x + u.z*c.y + f.z*c.z + t.z;
        toContour.push_back(v);
    }
    return toContour;
}



///////////////////////////////////////////////////////////////////////////////
// project contour from the from-point to the to-point of the knot path
///////////////////////////////////////////////////////////////////////////////
std::vector<TrefoilKnot::vec3> TrefoilKnot::projectContour(std::vector<TrefoilKnot::vec3>& fromContour, int fromIndex, int toIndex)
{
    // define plane: n.p + d = 0, at to-point
    vec3 p = path[toIndex];     // a point on the plane, (x0,y0,z0)
    vec3 lineDir = path[toIndex] - path[fromIndex];
    vec3 n = pathDirections[toIndex]; // normal of plane, (a,b,c)
    float d = -n.dot(p);            // d = -(a*x0 + b*y0 + c*z0)
    float dot1 = n.dot(lineDir);    // a*x + b*y + c*z

    std::vector<vec3> toContour;
    for(int i = 0; i <= sideCount; ++i)
    {
        // define line: c + v * t
        // find t = -(a*x0 + b*y0 + c*z0 + d) / (a*x + b*y + c*z)
        vec3 c = fromContour[i];
        float dot2 = n.dot(c);  // a*x0 + b*y0 + c*z0
        float t = -(dot2 + d) / dot1;
        c += lineDir * t;  // find the intersect point on the line using t
        toContour.push_back(c);
    }

    return toContour;
}



///////////////////////////////////////////////////////////////////////////////
// build vertices of trefoil knot with smooth shading using parametric equation
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::buildVerticesSmooth()
{
    // clear memory of prev arrays
    clearArrays();

    this->buildPath();
    this->buildContour();

    if(this->mode == 0) // projection mode
    {
        // add first contour
        this->transformFirstContour();
        std::vector<vec3> circle = contour;
        vec3 p = path[0];
        for(int i = 0; i <= sideCount; ++i)
        {
            vec3 c = circle[i];
            vec3 n = (c - p).normalize();
            addVertex(c.x, c.y, c.z);
            addNormal(n.x, n.y, n.z);
            addTexCoord(0, (float)i/sideCount);
        }
        // add rest of contours
        for(int i = 0; i < sectorCount; ++i)
        {
            vec3 p = path[i+1];
            circle = projectContour(circle, i, i+1);
            for(int j = 0; j <= sideCount; ++j)
            {
                vec3 c = circle[j];
                vec3 n = (c - p).normalize();
                addVertex(c.x, c.y, c.z);
                addNormal(n.x, n.y, n.z);
                addTexCoord((float)(i+1)/sectorCount, (float)j/sideCount);
            }
        }
        // first and last contour vertices are not aligned
        // rotate each contour with delta angle to align contours after loop
        alignContoursSmooth();
    }
    else // transform mode
    {
        // add contours
        for(int i = -1; i < sectorCount; ++i)
        {
            std::vector<vec3> circle = transformContour(i, i+1);
            for(int j = 0; j <= sideCount; ++j)
            {
                vec3 c = circle[j];
                vec3 p = path[0];
                vec3 n = (c - p).normalize();
                addVertex(c.x, c.y, c.z);
                addNormal(n.x, n.y, n.z);
                addTexCoord((float)(i+1)/sectorCount, (float)j/sideCount);
            }
        }
    }

    // indices
    //  k1--k2
    //  |  / |
    //  | /  |
    //  k1+1-k2+1
    unsigned int k1, k2;
    for(int i = 0; i < sectorCount; ++i)
    {
        k1 = i * (sideCount + 1);       // beginning of current sector
        k2 = k1 + sideCount + 1;        // beginning of next sector

        for(int j = 0; j < sideCount; ++j, ++k1, ++k2)
        {
            // 2 triangles per sector
            addIndices(k1, k1+1, k2);   // k1---k+1---k2
            addIndices(k2, k1+1, k2+1); // k2---k1+1--k2+1

            // vertical lines for all sides
            lineIndices.push_back(k1);
            lineIndices.push_back(k1+1);
            // horizontal lines for all sides
            lineIndices.push_back(k1);
            lineIndices.push_back(k2);
        }
    }

    // generate interleaved vertex array as well
    buildInterleavedVertices();

    // change up axis from Z-axis to the given
    if(this->upAxis != 3)
        changeUpAxis(3, this->upAxis);
}



///////////////////////////////////////////////////////////////////////////////
// generate vertices with flat shading
// each triangle is independent (no shared vertices)
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::buildVerticesFlat()
{
    // compute all vertices first, each vertex contains (x,y,z,s,t) except normal
    Vertex vertex;
    std::vector<Vertex> tmpVertices;

    // clear memory of prev arrays
    clearArrays();

    this->buildPath();
    this->buildContour();

    if(this->mode == 0) // projection mode
    {
        // add first contour
        this->transformFirstContour();
        std::vector<vec3> circle = contour;
        for(int i = 0; i <= sideCount; ++i)
        {
            vec3 c = circle[i];
            vertex.x = c.x;
            vertex.y = c.y;
            vertex.z = c.z;
            vertex.s = 0;
            vertex.t = (float)i/sideCount;
            tmpVertices.push_back(vertex);
        }
        // add rest of contours
        for(int i = 0; i < sectorCount; ++i)
        {
            circle = projectContour(circle, i, i+1);
            for(int j = 0; j <= sideCount; ++j)
            {
                vec3 c = circle[j];
                vertex.x = c.x;
                vertex.y = c.y;
                vertex.z = c.z;
                vertex.s = (float)(i+1)/sectorCount;
                vertex.t = (float)j/sideCount;
                tmpVertices.push_back(vertex);
            }
        }
        // first and last contour vertices are not aligned
        // rotate each contour with delta angle to align contours after loop
        alignContoursFlat(tmpVertices);
    }
    else // transform mode
    {
        for(int i = -1; i < sectorCount; ++i)
        {
            std::vector<vec3> circle = transformContour(i, i+1);
            for(int j = 0; j <= sideCount; ++j)
            {
                vec3 c = circle[j];
                vertex.x = c.x;
                vertex.y = c.y;
                vertex.z = c.z;
                vertex.s = (float)(i+1)/sectorCount;
                vertex.t = (float)j/sideCount;
                tmpVertices.push_back(vertex);
            }
        }
    }

    Vertex v1, v2, v3, v4;                          // 4 vertex positions and tex coords
    int i, j, k, vi1, vi2;
    int index = 0;                                  // index for vertex
    for(i = 0; i < sectorCount; ++i)
    {
        vi1 = i * (sideCount + 1);                // index of tmpVertices
        vi2 = (i + 1) * (sideCount + 1);

        for(j = 0; j < sideCount; ++j, ++vi1, ++vi2)
        {
            // get 4 vertices per side
            //  v1--v3
            //  |    |
            //  v2--v4
            v1 = tmpVertices[vi1];
            v2 = tmpVertices[vi1 + 1];
            v3 = tmpVertices[vi2];
            v4 = tmpVertices[vi2 + 1];

            // store 2 triangles (quad) per side
            // put quad vertices: v1-v2-v3-v4
            addVertex(v1.x, v1.y, v1.z);
            addVertex(v2.x, v2.y, v2.z);
            addVertex(v3.x, v3.y, v3.z);
            addVertex(v4.x, v4.y, v4.z);

            // put tex coords of quad
            addTexCoord(v1.s, v1.t);
            addTexCoord(v2.s, v2.t);
            addTexCoord(v3.s, v3.t);
            addTexCoord(v4.s, v4.t);

            // put normal
            vec3 n = computeFaceNormal(v1.x,v1.y,v1.z, v2.x,v2.y,v2.z, v3.x,v3.y,v3.z);
            for(k = 0; k < 4; ++k)  // same normals for 4 vertices
            {
                addNormal(n.x, n.y, n.z);
            }

            // put indices of quad (2 triangles)
            addIndices(index, index+1, index+2);
            addIndices(index+2, index+1, index+3);

            // indices for lines
            lineIndices.push_back(index);
            lineIndices.push_back(index+1);
            lineIndices.push_back(index);
            lineIndices.push_back(index+2);

            index += 4;     // for next
        }
    }

    // generate interleaved vertex array as well
    buildInterleavedVertices();

    // change up axis from Z-axis to the given
    if(this->upAxis != 3)
        changeUpAxis(3, this->upAxis);
}



///////////////////////////////////////////////////////////////////////////////
// align contour vertices, so first and last contours are aligned
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::alignContoursSmooth()
{
    // compute angle between first contour and last contour
    int lastIndex = sectorCount * (sideCount + 1) * 3;
    vec3 f1 = vec3(vertices[0], vertices[1], vertices[2]) - path[0];
    f1.normalize();
    vec3 f2 = vec3(vertices[lastIndex], vertices[lastIndex+1], vertices[lastIndex+2]) - path[0];
    f2.normalize();
    vec3 d = pathDirections[0];
    float dot = f1.dot(f2);
    vec3 cross = f1.cross(f2);
    float det = d.dot(cross);
    float a = atan2(-det, -dot) + PI; // angle range: 0 ~ 2PI
    //std::cout << "angle: " << a << ", " << a * RAD2DEG << std::endl;

    // if angle is small no need to rotate
    if(a > -EPSILON && a < EPSILON)
    {
        // make sure the last vertices are same as the first
        for(int i = 0, ii = lastIndex; i <= (sideCount*3); i+=3, ii += 3)
        {
            vertices[ii]   = vertices[i];
            vertices[ii+1] = vertices[i+1];
            vertices[ii+2] = vertices[i+2];
            normals[ii]   = normals[i];
            normals[ii+1] = normals[i+1];
            normals[ii+2] = normals[i+2];
        }
        return;
    }

    // find delta angle per each path
    float deltaAngle = 0;
    if(a <= PI)
        deltaAngle = -a / sectorCount;
    else
        deltaAngle = (2 * PI - a) / sectorCount;
    //std::cout << "delta angle: " << deltaAngle << ", " << (a -PI) << std::endl;

    a = 0;
    for(int i = 1; i < sectorCount; ++i)
    {
        d = pathDirections[i];
        vec3 center = path[i];
        a += deltaAngle;
        float c = cosf(a);
        float s = sinf(a);
        //std::cout << "a: " << a << ", " << a * RAD2DEG << std::endl;
        for(int j = 0, k = (sideCount + 1) * i * 3; j <= sideCount; ++j, k += 3)
        {
            //int k = (sideCount + 1) * i * 3 + (j * 3);
            vec3 v = vec3(vertices[k], vertices[k+1], vertices[k+2]);
            vec3 n = vec3(normals[k], normals[k+1], normals[k+2]);
            // rotate with rodrigues formula
            v -= center;
            v = d * (1 - c) * v.dot(d) + v * c + d.cross(v) * s;
            v += center;
            vertices[k]   = v.x;
            vertices[k+1] = v.y;
            vertices[k+2] = v.z;
            n = d * (1 - c) * n.dot(d) + n * c + d.cross(n) * s;
            normals[k]   = n.x;
            normals[k+1] = n.y;
            normals[k+2] = n.z;
        }
    }

    // snap the last contour with last contour
    for(int i = 0, ii = lastIndex; i <= (sideCount*3); i += 3 , ii += 3)
    {
        vertices[ii]   = vertices[i];
        vertices[ii+1] = vertices[i+1];
        vertices[ii+2] = vertices[i+2];
        normals[ii]   = normals[i];
        normals[ii+1] = normals[i+1];
        normals[ii+2] = normals[i+2];
    }
}



///////////////////////////////////////////////////////////////////////////////
// align contour vertices, so first and last contours are aligned
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::alignContoursFlat(std::vector<TrefoilKnot::Vertex>& vertices)
{
    // compute angle between first contour and last contour
    int lastIndex = sectorCount * (sideCount + 1);
    vec3 f1 = vec3(vertices[0].x, vertices[0].y, vertices[0].z) - path[0];
    f1.normalize();
    vec3 f2 = vec3(vertices[lastIndex].x, vertices[lastIndex].y, vertices[lastIndex].z) - path[0];
    f2.normalize();
    vec3 n = pathDirections[0];
    float dot = f1.dot(f2);
    vec3 cross = f1.cross(f2);
    float det = n.dot(cross);
    float a = atan2(-det, -dot) + PI; // angle range: 0 ~ 2PI
    //std::cout << "angle: " << a << ", " << a * RAD2DEG << std::endl;

    // if angle is small no need to rotate
    if(a > -EPSILON && a < EPSILON)
    {
        // make sure the last vertices are same as the first
        for(int i = 0; i <= sideCount; ++i)
        {
            vertices[lastIndex + i].x = vertices[i].x;
            vertices[lastIndex + i].y = vertices[i].y;
            vertices[lastIndex + i].z = vertices[i].z;
        }
        return;
    }

    // find delta angle per each path
    float deltaAngle = 0;
    if(a <= PI)
        deltaAngle = -a / sectorCount;
    else
        deltaAngle = (2 * PI - a) / sectorCount;
    //std::cout << "delta angle: " << deltaAngle << ", " << (a -PI) << std::endl;

    a = 0;
    for(int i = 1, ii = sideCount+1; i < sectorCount; ++i, ii += (sideCount+1))
    {
        n = pathDirections[i];
        vec3 center = path[i];
        a += deltaAngle;
        float c = cosf(a);
        float s = sinf(a);
        for(int j = 0; j <= sideCount; ++j)
        {
            vec3 p = vec3(vertices[ii + j].x, vertices[ii + j].y, vertices[ii + j].z);
            // rotate with rodrigues formula
            p -= center;
            p = n * (1 - c) * p.dot(n) + p * c + n.cross(p) * s;
            p += center;
            vertices[ii + j].x = p.x;
            vertices[ii + j].y = p.y;
            vertices[ii + j].z = p.z;
        }
    }

    // snap the last contour with last contour
    for(int i = 0; i <= sideCount; ++i)
    {
        vertices[lastIndex + i].x = vertices[i].x;
        vertices[lastIndex + i].y = vertices[i].y;
        vertices[lastIndex + i].z = vertices[i].z;
    }
}



///////////////////////////////////////////////////////////////////////////////
// generate interleaved vertices: V/N/T
// stride must be 32 bytes
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::buildInterleavedVertices()
{
    std::vector<float>().swap(interleavedVertices);

    std::size_t i, j;
    std::size_t count = vertices.size();
    for(i = 0, j = 0; i < count; i += 3, j += 2)
    {
        interleavedVertices.push_back(vertices[i]);
        interleavedVertices.push_back(vertices[i+1]);
        interleavedVertices.push_back(vertices[i+2]);

        interleavedVertices.push_back(normals[i]);
        interleavedVertices.push_back(normals[i+1]);
        interleavedVertices.push_back(normals[i+2]);

        interleavedVertices.push_back(texCoords[j]);
        interleavedVertices.push_back(texCoords[j+1]);
    }
}



///////////////////////////////////////////////////////////////////////////////
// transform vertex/normal (x,y,z) coords
// assume from/to values are validated: 1~3 and from != to
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::changeUpAxis(int from, int to)
{
    // initial transform matrix cols
    vec3 l(1, 0, 0);                    // x-axis (left)
    vec3 u(0, 1, 0);                    // y-axis (up)
    vec3 f(0, 0, 1);                    // z-axis (forward)

    // X -> Y
    if(from == 1 && to == 2)
    {
        l.x = 0.0f;  l.y = 1.0f;
        u.x = 1.0f;  l.y = 0.0f;
    }
    // X -> Z
    else if(from == 1 && to == 3)
    {
        l.x = 0.0f;  l.z = 1.0f;
        f.x = -1.0f; f.z = 0.0f;
    }
    // Y -> X
    else if(from == 2 && to == 1)
    {
        l.x = 0.0f;  l.y = -1.0f;
        u.x = 1.0f;  u.y = 0.0f;
    }
    // Y -> Z
    else if(from == 2 && to == 3)
    {
        u.y = 0.0f;  u.z = 1.0f;
        f.y = -1.0f; f.z = 0.0f;
    }
    //  Z -> X
    else if(from == 3 && to == 1)
    {
        l.x = 0.0f;  l.z = 1.0f;
        f.x = 1.0f;  f.z = 0.0f;
    }
    // Z -> Y
    else
    {
        u.y = 0.0f;  u.z = -1.0f;
        f.y = 1.0f;  f.z = 0.0f;
    }

    std::size_t i, j;
    std::size_t count = vertices.size();
    vec3 v, n;
    for(i = 0, j = 0; i < count; i += 3, j += 8)
    {
        // transform vertices
        v.x = vertices[i];
        v.y = vertices[i+1];
        v.z = vertices[i+2];
        vertices[i]   = l.x * v.x + u.x * v.y + f.x * v.z;  // x
        vertices[i+1] = l.y * v.x + u.y * v.y + f.y * v.z;  // y
        vertices[i+2] = l.z * v.x + u.z * v.y + f.z * v.z;  // z

        // transform normals
        n.x = normals[i];
        n.y = normals[i+1];
        n.z = normals[i+2];
        normals[i]   = l.x * n.x + u.x * n.y + f.x * n.z;   // nx
        normals[i+1] = l.y * n.x + u.y * n.y + f.y * n.z;   // ny
        normals[i+2] = l.z * n.x + u.z * n.y + f.z * n.z;   // nz

        // trnasform interleaved array
        interleavedVertices[j]   = vertices[i];
        interleavedVertices[j+1] = vertices[i+1];
        interleavedVertices[j+2] = vertices[i+2];
        interleavedVertices[j+3] = normals[i];
        interleavedVertices[j+4] = normals[i+1];
        interleavedVertices[j+5] = normals[i+2];
    }

    // change path as well
    for(i = 0, j = 0; i < path.size(); ++i)
    {
        vec3 p = path[i];
        v.x = l.x * p.x + u.x * p.y + f.x * p.z;  // x
        v.y = l.y * p.x + u.y * p.y + f.y * p.z;  // y
        v.z = l.z * p.x + u.z * p.y + f.z * p.z;  // z
        path[i] = v;
    }
}



///////////////////////////////////////////////////////////////////////////////
// add single vertex to array
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::addVertex(float x, float y, float z)
{
    vertices.push_back(x);
    vertices.push_back(y);
    vertices.push_back(z);
}



///////////////////////////////////////////////////////////////////////////////
// add single normal to array
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::addNormal(float nx, float ny, float nz)
{
    normals.push_back(nx);
    normals.push_back(ny);
    normals.push_back(nz);
}



///////////////////////////////////////////////////////////////////////////////
// add single texture coord to array
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::addTexCoord(float s, float t)
{
    texCoords.push_back(s);
    texCoords.push_back(t);
}



///////////////////////////////////////////////////////////////////////////////
// add 3 indices to array
///////////////////////////////////////////////////////////////////////////////
void TrefoilKnot::addIndices(unsigned int i1, unsigned int i2, unsigned int i3)
{
    indices.push_back(i1);
    indices.push_back(i2);
    indices.push_back(i3);
}



///////////////////////////////////////////////////////////////////////////////
// return face normal of a triangle v1-v2-v3
// if a triangle has no surface (normal length = 0), then return a zero vector
///////////////////////////////////////////////////////////////////////////////
TrefoilKnot::vec3 TrefoilKnot::computeFaceNormal(float x1, float y1, float z1,
                                             float x2, float y2, float z2,
                                             float x3, float y3, float z3)
{
    // find 2 edge vectors: v1-v2, v1-v3
    vec3 v1(x1, y1, z1);
    vec3 v2(x2, y2, z2);
    vec3 v3(x3, y3, z3);
    vec3 e1 = v2 - v1;
    vec3 e2 = v3 - v1;

    // cross product: e1 x e2
    vec3 n = e1.cross(e2).normalize();
    return n;
}
