glSpline/src/bspline.cpp

#include "bspline.h"
#include <cmath>

BSpline::BSpline(
        ESpline typeSpline,
        int k,
        EPolygone typePolygone,
        EVecteurNodal typeVecteur,
        int nbPtsCtrlX,
        int nbPtsCtrlZ,
        float stepX,
        float stepZ)
{
    // Simple vérification
    switch (typeSpline) {
    case ESpline::SPLINE1D:
        assert(k <= nbPtsCtrlX);
        assert(typePolygone != EPolygone::GAUSSIEN3D);
        assert(typePolygone != EPolygone::FIXE3D);
        assert(typePolygone != EPolygone::FIXE3D);
        break;
    case ESpline::SPLINE2D:
        assert(k <= nbPtsCtrlX && k <= nbPtsCtrlZ);
        assert(typePolygone != EPolygone::FIXE2D);
        assert(typePolygone != EPolygone::RANDOM2D);
        break;
    }

    // Initialisation des valeurs
    _k = k;
    _stepX = stepX;
    _stepZ = stepZ;
    srand(time(NULL));

    // Génération du polygone de contrôle
    construirePolygone(typePolygone, nbPtsCtrlX, nbPtsCtrlZ);

    // Calcul du vecteur nodal de nbPtsCtrl valeurs
    construireVecteurNodal(typeVecteur);

    // Appel de la construction de la B-Spline
    switch (typeSpline) {
    case ESpline::SPLINE1D:
        calculerSpline1D();
        break;
    case ESpline::SPLINE2D:
        calculerSpline2D();
        break;
    default:
        std::cerr << "Type de spline inconnu" << std::endl;
        break;
    }

}

void BSpline::calculerSpline1D()
{
    // Init
    int dec = 0, m = _k - 1;
    uint i = 0;
    double u = 0;
    std::vector<glm::vec3> vecteursPointsControle (_nbPtsCtrlX);

    // Parcours de u
    for(u = _vecteurNodal[m]; u <= _vecteurNodal[_nbPtsCtrlX]; u += _stepX){
        dec = 0;
        // Calcul du décalage
        for(int i = _k; u > _vecteurNodal[i]; ++i)
            ++dec;

        // Récupération des k points de contrôles nécessaires à la floraison
        for(int i = 0; i < _k; i++){
            vecteursPointsControle[i] = _pointsDeControle[dec + i];
        }

        // Appel de la floraison récursive
        glm::vec3 newVertex = floraison(u, dec, _k, vecteursPointsControle);

        // Sauvegardes des paramètres de la B-spline
        _vertices.insert(_vertices.end(), {newVertex.x, newVertex.y, newVertex.z});
        _normals.insert(_normals.end(), {1, -1, -1});
        _indices.insert(_indices.end(), {i, i, ++i});
    }

    // Enlever la ligne d'indices en trop
    for(int i = 0; i < 3; ++i)
        _indices.pop_back();

    // Ajouter le polygone de controle aux vecteurs
    for(unsigned int i = 0; i < _pointsDeControle.size(); ++i) {
        _indices.insert(_indices.end(), {static_cast<unsigned int>(_vertices.size())/3, static_cast<unsigned int>(_vertices.size()/3), static_cast<unsigned int>(_vertices.size())/3+1});
        _vertices.insert(_vertices.end(), {_pointsDeControle[i].x, _pointsDeControle[i].y, _pointsDeControle[i].z});
        _normals.insert(_normals.end(), {0, 0, 0});
    }

    // Enlever la ligne d'indices en trop
    for(int i = 0; i < 3; ++i)
        _indices.pop_back();
}


void BSpline::calculerSpline2D()
{
    // Init
    int decX = 0, decZ = 0, m = _k - 1, nbPtsMaillageX = 0, nbPtsMaillageZ = 0;
    uint i = 0, n = 0;
    double u = 0, v = 0;
    //std::vector<glm::vec3> vecteursPointsControle (_k);
    std::vector<glm::vec3> generatrice;
    std::vector<std::vector<glm::vec3>> pointsInfluents;
    std::vector<glm::vec3> vecteurHorizontal;

    // Évaluation de p(u, v), revient à évaluer la génératrice sur les directrices
    //Parcours de u
    for(u = _vecteurNodal[m]; u <= _vecteurNodal[_nbPtsCtrlX]; u += _stepX){
        // Parcours de v
        for(v = _vecteurNodal[m]; v <= _vecteurNodal[_nbPtsCtrlZ]; v += _stepZ){
            // Réinitialisation
            vecteurHorizontal.clear();
            pointsInfluents.clear();
            generatrice.clear();
            decX = 0;
            decZ = 0;

            // Calcul du décalage
            for(int i = _k; u > _vecteurNodal[i]; ++i)
                ++decX;
            for(int j = _k; v > _vecteurNodal[j]; ++j)
                ++decZ;

            // Calcul des points influents
            pointsInfluents.clear();
            for(int i = 0; i < _k; ++i) {
                vecteurHorizontal.clear();
                for(int j = 0; j < _k; ++j) {
                    glm::vec3 pt = _pointsDeControle[((decX + i) * _nbPtsCtrlZ) + (decZ + j)];
                    vecteurHorizontal.push_back(pt);
                }
                pointsInfluents.push_back(vecteurHorizontal);
            }

            // m+1 directrices floraisons
            for(int i = 0; i < _k; ++i) {
                generatrice.push_back(floraison(v, decZ, _k, pointsInfluents[i]));
            }


            // Calcul de b-spline sur la génératrice (locale)
            // Appel de la floraison récursive
            glm::vec3 newVertex = floraison(u, decX, _k, generatrice);

            // Sauvegardes des paramètres de la B-spline

            _vertices.insert(_vertices.end(), {newVertex.x, newVertex.y, newVertex.z});
            //_normals.insert(_normals.end(), {1, -1, -1});
            //_indices.insert(_indices.end(), {i, i, ++i});


            // Incrémenter taille en z
            if(u == _vecteurNodal[m]) nbPtsMaillageZ++;
        }

        nbPtsMaillageX++;
    }

    std::cout << nbPtsMaillageX * nbPtsMaillageZ << " INDICES CALCULÉS" << std::endl;

    // Maillage
    for(unsigned int i = 0; i < _vertices.size() / 3; ++i) {
        // Indices
        int x = i / nbPtsMaillageZ;
        int z = i % nbPtsMaillageX;
        // Sommets
        uint sommetA = x * nbPtsMaillageX + z;
        uint sommetB = x * nbPtsMaillageX + z+1;
        uint sommetC = (x+1) * nbPtsMaillageX + z;
        uint sommetD = (x+1) * nbPtsMaillageX + z+1;
        // Triangles
        if (x != nbPtsMaillageX - 1 && z != nbPtsMaillageZ - 1) {
            // Faces
            _indices.insert(_indices.end(), {sommetA, sommetD, sommetB}); // Triangle 1
            _indices.insert(_indices.end(), {sommetA, sommetC, sommetD}); // Triangle 2
        }
    }

    std::cout << _vertices.size() / 3 * 2 << " TRIANGLES CALCULÉS" << std::endl;

    // Normales
    _normals.assign(_vertices.size(), 0.0f);

    for(int i = 0; i < _indices.size(); i+=3) {
        // Récupération des sommets
        uint indiceSommetA = _indices[i];
        uint indiceSommetB = _indices[i+1];
        uint indiceSommetC = _indices[i+2];
        glm::vec3 sommetA = glm::vec3(_vertices[indiceSommetA*3], _vertices[indiceSommetA*3+1], _vertices[indiceSommetA*3+2]);
        glm::vec3 sommetB = glm::vec3(_vertices[indiceSommetB*3], _vertices[indiceSommetB*3+1], _vertices[indiceSommetB*3+2]);
        glm::vec3 sommetC = glm::vec3(_vertices[indiceSommetC*3], _vertices[indiceSommetC*3+1], _vertices[indiceSommetC*3+2]);

        // Calcul de la normale d'un triangle TODO: vérifier orientation normale
        glm::vec3 normale = glm::cross(sommetB - sommetA, sommetC - sommetA);

        // Ajout de la normale de la face aux normales des sommets participants
        _normals[indiceSommetA*3] += normale.x; _normals[indiceSommetA*3+1] += normale.y; _normals[indiceSommetA*3+2] += normale.z;
        _normals[indiceSommetB*3] += normale.x; _normals[indiceSommetB*3+1] += normale.y; _normals[indiceSommetB*3+2] += normale.z;
        _normals[indiceSommetC*3] += normale.x; _normals[indiceSommetC*3+1] += normale.y; _normals[indiceSommetC*3+2] += normale.z;
    }

    // Normalisation des normales en chaque sommet
    for(int i = 0; i < _vertices.size(); i+=3) {
        glm::vec3 normale = glm::normalize(glm::vec3(_normals[i], _normals[i+1], _normals[i+2]));
        _normals.at(i) = normale.x;
        _normals.at(i+1) = normale.y;
        _normals.at(i+2) = normale.z;
    }

    std::cout << _normals.size()/3 << " NORMALES CALCULÉES" << std::endl;

    // Ajouter le polygone de controle aux vecteurs
    int offset = static_cast<unsigned int>(_vertices.size())/3;

    for(unsigned int i = 0; i < _pointsDeControle.size(); ++i) {
        // Indices
        int x = i / _nbPtsCtrlZ;
        int z = i % _nbPtsCtrlX;
        //std::cout << "z: " << z << "; x: " << x << std::endl;

        // Sommets
        uint sommetA = x * _nbPtsCtrlX + z       + offset;
        uint sommetB = x * _nbPtsCtrlX + z+1     + offset;
        uint sommetC = (x+1) * _nbPtsCtrlX + z   + offset;
        //uint sommetD = (x+1) * _nbPtsCtrlX + z+1 + offset;

        // Triangles
        if (x != _nbPtsCtrlX - 1 && z != _nbPtsCtrlZ - 1) {
            // Faces
            //_indices.insert(_indices.end(), {sommetA, sommetD, sommetB}); // Triangle 1
            //_indices.insert(_indices.end(), {sommetA, sommetC, sommetD}); // Triangle 2
        }

        // Arêtes
        if (x != _nbPtsCtrlX - 1)
            _indices.insert(_indices.end(), {sommetA, sommetC, sommetA}); // Arête Horizontale
        if (z != _nbPtsCtrlZ - 1)
            _indices.insert(_indices.end(), {sommetA, sommetB, sommetA}); // Arête Verticale
        _vertices.insert(_vertices.end(), {_pointsDeControle[i].x, _pointsDeControle[i].y, _pointsDeControle[i].z});
        _normals.insert(_normals.end(), {0, 0, 0});
    }
    std::cout << "POLYGONE AJOUTÉ" << std::endl;
}

glm::vec3 BSpline::floraison(float u, int dec, int k, std::vector<glm::vec3> vecteursPointsControle)
{
    if(k <= 1)
        return vecteursPointsControle.front();

    for(int i = 0; i < k - 1; ++i) {
        float max = _vecteurNodal[dec + k + i];
        float min = _vecteurNodal[dec + 1 + i];
        vecteursPointsControle[i] =
                ((( max - u) / (max - min)) * vecteursPointsControle[i]) +
                ((( u - min) / (max - min)) * vecteursPointsControle[i+1]);
    }

    return floraison(u, dec+1, k-1, vecteursPointsControle);
}

void BSpline::construirePolygone(EPolygone typePolygone, int nbPtsCtrlX, int nbPtsCtrlZ)
{
    // Intervalle pour afficher dans le repère caméra
    float xMin = -3, xMax =  3;
    float yMin = -2, yMax =  2;
    float zMin = -3, zMax = -6;

    switch (typePolygone) {
    case EPolygone::FIXE2D:
        _pointsDeControle.push_back(glm::vec3(-2, 0, -5));
        _pointsDeControle.push_back(glm::vec3(-1, 0, -5));
        _pointsDeControle.push_back(glm::vec3(-1, 1, -5));
        _pointsDeControle.push_back(glm::vec3( 1, 1, -5));
        _pointsDeControle.push_back(glm::vec3( 1, 0, -5));
        _pointsDeControle.push_back(glm::vec3( 2, 0, -5));
        _nbPtsCtrlX = 6;
        _nbPtsCtrlZ = 0;
        break;

    case EPolygone::RANDOM2D: {
        float x, y, z;
        for(int i = 0; i < nbPtsCtrlX; ++i) {
            //x = xMin + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(xMax-xMin)));
            x = (i * xMin + (nbPtsCtrlX - i) * xMax) / (xMax - xMin);
            y = yMin + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(yMax-yMin)));
            z = zMin + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(zMax-zMin)));
            _pointsDeControle.push_back(glm::vec3(x, y, z));
        }
        _nbPtsCtrlX = _pointsDeControle.size();
        _nbPtsCtrlZ = 0;
        break;
    }

    case EPolygone::GAUSSIEN3D: {
        float sigma = 0.3f;
        float r, s = 2.0f * sigma * sigma;
        for(int x = - (nbPtsCtrlX / 2) + 1; x < nbPtsCtrlX / 2; ++x) {
            for(int z = - (nbPtsCtrlZ / 2) + 1; z < nbPtsCtrlZ / 2; ++z) {
                r = sqrt(x * x + z * z);
                _pointsDeControle.push_back(glm::vec3( x, (exp(-(r * r) / s)) / (M_PI * s), z));
            }
        }
        _nbPtsCtrlX = nbPtsCtrlX - 1 - (nbPtsCtrlX % 2);
        _nbPtsCtrlZ = nbPtsCtrlZ - 1 - (nbPtsCtrlZ % 2);
        break;
    }

    case EPolygone::FIXE3D: {
        float x, y, z;
        for(int i = 0; i < nbPtsCtrlX; ++i) {
            for(int j = 0; j < nbPtsCtrlZ; ++j) {
                x = (i * xMin + (nbPtsCtrlX - i) * xMax) / (xMax - xMin);
                y = fmod(x + z, 3);
                z = (j * zMin + (nbPtsCtrlZ - j) * zMax) / (zMin - zMax);
                _pointsDeControle.push_back(glm::vec3(x, y, z));
            }
        }
        _nbPtsCtrlX = nbPtsCtrlX;
        _nbPtsCtrlZ = nbPtsCtrlZ;
        break;
    }

    case EPolygone::RANDOM3D: {
        float x, y, z;
        for(int i = 0; i < nbPtsCtrlX; ++i) {
            for(int j = 0; j < nbPtsCtrlZ; ++j) {
                x = (i * xMin + (nbPtsCtrlX - i) * xMax) / (xMax - xMin);
                y = yMin + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(yMax-yMin))) / 2;
                z = (j * zMin + (nbPtsCtrlZ - j) * zMax) / (zMin - zMax);
                _pointsDeControle.push_back(glm::vec3(x, y, z));
            }
        }
        _nbPtsCtrlX = nbPtsCtrlX;
        _nbPtsCtrlZ = nbPtsCtrlZ;
        break;
    }

    default:
        std::cerr << "Type de polygone inconnu" << std::endl;
        break;
    }
}

void BSpline::construireVecteurNodal(EVecteurNodal typeVecteur)
{
    switch (typeVecteur) {
    case EVecteurNodal::UNIFORME:
        for(int i = 0; i <= _k + _nbPtsCtrlX; ++i) {
            _vecteurNodal.push_back(i);
        }
        break;
    case EVecteurNodal::OUVERTUNIFORME:
        for(int i = 0; i <= _k + _nbPtsCtrlX; ++i) {
            if (i < _k)
                _vecteurNodal.push_back(0);
            else if (i > _nbPtsCtrlX)
                _vecteurNodal.push_back(_nbPtsCtrlX - _k + 1);
            else
                _vecteurNodal.push_back(i - _k + 1);
            std::cout << _vecteurNodal.back() << " ";
        }
        break;
    case EVecteurNodal::QUELCONQUE: {
        int r = 0;
        _vecteurNodal.push_back(0);
        for(int i = 1; i <= _k + _nbPtsCtrlX; ++i) {
            r = rand() / (RAND_MAX/5);
            _vecteurNodal.push_back(_vecteurNodal.back() + r);
        }
        break;
    }
    default:
        std::cerr << "Type de vecteur nodal inconnu" << std::endl;
        break;
    }

}