基于C++实现各种系统仿真模拟

本文主要介绍各种有关动画模拟的模型，直接进入正题

骨骼动画模拟

概念引入

对于网格体而言有不少实现动画的方式。直接对顶点进行操作也就是顶点动画，适用于一些比较简单的植物摆动、水面波动效果。此外，还有在两个网格之间进行插值的morphing动画；但它们本质上都是对顶点进行操作。

而在某些情况下，我们认为某些区域的顶点具有关联性，并且希望能够对其整体进行控制，比如整体移动腿部的顶点。为此，我们引入了骨骼动画系统。我们定义某个顶点可由一个或多个骨骼控制，每个骨骼对顶点的贡献有着不同的权重，对于同一个顶点而言，所有对其产生影响的骨骼的权重加起来应该为1。这样，我们就可以通过修改骨骼的平移、旋转、缩放数据，来控制顶点的变换。

对于骨骼而言，不同的骨骼存在父子关系。比如，当我们旋转手臂的时候，也会带动手部和手指的旋转。这些骨骼的父子关系之间形成了一棵骨骼“树”结构。

常用变换

骨骼动画中存在着大量坐标空间和坐标变换，搞懂这些空间或变换是理解骨骼动画最核心的部分。

在此之前，我们应该对坐标空间/姿态和坐标变换这两个概念有一定的了解。为了更方便描述这两个概念，我们可以用相机空间来举例。我们的相机位于空间中的某一点，我们认为它此时平移、旋转、缩放对应的矩阵就是相机空间，或者说相机的姿态矩阵。如名字所见，这反映的是相机当前的状态。

那么，如果我们希望将点/向量变换到相机空间呢？我们需要用到的变换矩阵和相机空间矩阵又是什么关系呢？此处，实际上可以得到一个一般性的结论：

● 如果希望将物体从世界空间变换到A空间，变换矩阵为A空间矩阵的逆矩阵；反之，从A空间变换到世界空间，变换矩阵即为A空间矩阵。

我们所谓的变换到某个空间，也就是找到当前物体在某个空间下的坐标表达。对于相机空间而言，我们以相机原点为例，它在相机空间的坐标应为(0,0,0)，这意味如果我们想要从世界空间的相机位置通过转换得到(0,0,0)这一结果，需要乘以相机姿态矩阵的逆矩阵，此时两个变换就刚好抵消得到0的结果。

对以上几个概念有了基本的了解后，我们开始引入骨骼变换中几个常见的空间/变换：绑定姿态、骨骼空间、offset变换矩阵、局部变换矩阵、全局变换矩阵、蒙皮变换矩阵等。

由于本文的重点是实现骨骼动画的渲染，因此，在此仅对相关的矩阵进行介绍。

BindPose

绑定姿态(bindPose) 也就是我们常说的T-Pose，这通常是由美术在建模软件中设定的，它定义了骨架的一种默认姿态，绑定矩阵存储了在这一姿态下，所有骨骼的变换数据。我们所有的骨骼动画变换都是在这一绑定姿态的基础上进行的。对于不同体型的角色而言，它们的高矮胖瘦是不同的，因而它们的绑定姿态也是不一样的。

Offset Matrix

根据上面的一般性结论，我们知道如果我们希望将物体从世界空间转换到绑定空间，我们需要乘以绑定姿态的逆矩阵。我们把这一矩阵称为Offset矩阵。这个矩阵在后面的计算会派上用场。

GlobalTransform & LocalTransform

还有两个比较重要的变换，一个是全局变换矩阵，一个是局部变换矩阵。它们反映了骨骼动画的每帧的变换。其中全局变换矩阵是是关节从根处变换到它最终所在位置所对应的变换矩阵，它通常用于骨骼动画渲染中，计算每个顶点的最终位置的时候；而拒不变换矩阵是指关节在父关节空间下的变换矩阵，它通常用于我们希望编辑某一关节的变换的时候，比如做出让头部抬起来这样的动作。

在最简单的情况下，已知一个关节的局部空间变换时，连乘它的所有父骨骼的局部空间变换矩阵，就能得到该关节的全局变换矩阵。某些格式实际上还可能有别的一些特别定义的矩阵，在实际实现的时候需要特别注意。

Mesh Transform

实际上确定每个顶点最终位置的是蒙皮变换矩阵。对于每个处于绑定姿态的顶点而言，它通过蒙皮矩阵变换得到最终位置。它和全局变换矩阵（GlobalTransform）的区别在于，全局变换矩阵是将顶点从root处变换到最终位置，而蒙皮矩阵是将顶点从绑定姿态变换到最终位置。因此，我们需要先将顶点变换到绑定空间（通过乘以offset矩阵），再变换到最终位置。

因此，蒙皮矩阵 = OffsetMatrix * GlobalTransform

此外，我们还需要考虑到一个顶点可能受到多个骨骼不同权重的影响。此处需要我们对所有蒙皮矩阵进行加权平均，得到最终的变换矩阵。

骨骼动画导入和绘制

我们可以把骨骼动画的加载分为三个部分。

第一部分，导入标准骨骼数据。对于不同的角色，如人和动物，我们有着不同的骨架，在这个过程中，我们导入多套骨架，包含每个骨架的名字，以及它对应的所有骨骼的名字，按顺序存储。我们得到了每个骨骼和其对应的索引下标，之后在连续地址存储骨骼矩阵的时候，我们将按此顺序存储（but,我贴出的代码只存了一份骨架）。

需要导入的主要为骨骼名字和对应的绑定姿态，索引id可在导入的过程中生成。

此处我们可以借助于assimp库，不过在此之前，需要对模型加载有一定了解(参考https://blog.csdn.net/ZJU_fish1996/article/details/90143844）。在获得了aiMesh的基础上，我们可以继续提取骨骼数据，同时存储offset矩阵：

struct Bone
{
    string m_name;
 
    QMatrix4x4 m_invBindPose;
    Bone(const string& name, const QMatrix4x4& pose)
        :m_name(name),  m_invBindPose(pose) { }
};
12345678

void GeometryEngine::processBone(const aiMesh* pMesh, vector<BoneVertexData>& vertices, bool bAdd)
{
    for (uint i = 0 ; i < pMesh->mNumBones ; i++)
    {
        string BoneName(pMesh->mBones[i]->mName.data);
        QMatrix4x4 dst;
        TransformMatrix(pMesh->mBones[i]->mOffsetMatrix, dst);
        CAnimationEngine::Inst()->AddBone(new Bone(BoneName, dst));
    }
}
12345678910

第二部分是导入带绑定骨骼的模型，但不包含动画。此时我们不仅获取每个点的法线、位置等信息，还存储了影响每个顶点的骨骼id和对应的权重，以生成带骨骼模型的vao/vbo相关数据。默认情况下，我们认为最多有四个骨骼影响同一顶点。具体而言，我们把顶点格式定义如下：

struct BoneVertexData
{
    QVector3D position;
    QVector3D normal;
    QVector3D tangent;
    QVector2D texcoord;
    QVector4D boneWeight;
    QVector4D boneIndex;
};

（注：骨骼索引应该存为int类型，由于我本地在向顶点着色器传int类型数据出现了一些问题，我将其存储为了浮点数，并在着色器中转换为ivec4）

相比起静态物体，带绑骨的模型多存储了boneWeight和boneIndex这两个信息，我们在读入了其它数据之后（见导入模型代码），再读入绑骨信息：

void GeometryEngine::processBone(const aiMesh* pMesh, vector<BoneVertexData>& vertices, bool bAdd)
{
    for (uint i = 0 ; i < pMesh->mNumBones ; i++)
    {
        string BoneName(pMesh->mBones[i]->mName.data);
        int boneIdx = CAnimationEngine::Inst()->GetBoneIndex(BoneName);
 
        if(boneIdx == -1) continue;
 
        for (uint j = 0 ; j < pMesh->mBones[i]->mNumWeights; j++)
        {
            size_t VertexID = pMesh->mBones[i]->mWeights[j].mVertexId;
            float Weight = pMesh->mBones[i]->mWeights[j].mWeight;
            if(VertexID >= vertices.size())
            {
                qDebug() << "err larger " << VertexID;
            }
            for(int k = 0;k < 4;k++)
            {
                if(vertices[VertexID].boneIndex[k] == INVALID_BONE_IDX)
                {
                    vertices[VertexID].boneIndex[k] = boneIdx;
                    vertices[VertexID].boneWeight[k] = Weight;
                    break;
                }
            }
        }
    }
}

第三部分是解析动画，并应用于对应的带绑骨的模型（原则上动画中的骨骼信息应该和应用的骨架匹配）。此处我们需要做的是存储每一帧，每个骨骼的蒙皮变换矩阵（在这里我们不考虑动画信息的压缩算法）。解析动画处使用了fbx自带的sdk（assimp也是可行的，但是需要自己解算全局变换矩阵，具体可以参考https://blog.csdn.net/ZJU_fish1996/article/details/52450008）：

CFbxImporter::CFbxImporter()
{
    InitSdk();
}
 
 
QMatrix4x4 TransformToQMatrix(const FbxAMatrix& mat)
{
    QMatrix4x4 res;
    res.setRow(0,{float(mat.Get(0,0)),float(mat.Get(0,1)),float(mat.Get(0,2)),float(mat.Get(0,3))});
    res.setRow(1,{float(mat.Get(1,0)),float(mat.Get(1,1)),float(mat.Get(1,2)),float(mat.Get(1,3))});
    res.setRow(2,{float(mat.Get(2,0)),float(mat.Get(2,1)),float(mat.Get(2,2)),float(mat.Get(2,3))});
    res.setRow(3,{float(mat.Get(3,0)),float(mat.Get(3,1)),float(mat.Get(3,2)),float(mat.Get(3,3))});
    res.setRow(3,{float(mat.Get(3,0)) * 0.01f,float(mat.Get(3,1)) * 0.01f,float(mat.Get(3,2)) * 0.01f,float(mat.Get(3,3))});
    res = res.transposed();
    return res;
}
 
bool CFbxImporter::InitSdk()
{
    pManager = FbxManager::Create();
    if( !pManager )
    {
        return false;
    }
 
    FbxIOSettings* ios = FbxIOSettings::Create(pManager, IOSROOT);
    pManager->SetIOSettings(ios);
 
    FbxString lPath = FbxGetApplicationDirectory();
    pManager->LoadPluginsDirectory(lPath.Buffer());
 
    pScene = FbxScene::Create(pManager, "Scene");
    if( !pScene )
    {
        return false;
    }
    qDebug() << "success init sdk " ;
    return true;
}
 
void CFbxImporter::ProcessAnimation()
{
    FbxAnimStack* pAnimStack = pScene->GetCurrentAnimationStack();
 
    FbxTime timePerFrame;
    timePerFrame.SetTime(0, 0, 0, 1, 0, pScene->GetGlobalSettings().GetTimeMode());
 
    const FbxTimeSpan animTimeSpan = pAnimStack->GetLocalTimeSpan();
    const FbxTime startTime = animTimeSpan.GetStart();
    const FbxTime endTime = animTimeSpan.GetStop();
    unsigned int frameNum = static_cast<unsigned int>(animTimeSpan.GetDuration().GetFrameCount(pScene->GetGlobalSettings().GetTimeMode())) + 1;
    int boneNum = CAnimationEngine::Inst()->GetBoneNum();
 
    if(pCurrentAnimator)
    {
        pCurrentAnimator->Init(boneNum, frameNum);
    }
 
    for (FbxNode* pNode : vecNodes)
    {
        unsigned int numFrames = 0;
        int boneIdx = CAnimationEngine::Inst()->GetBoneIndex(pNode->GetName());
        if(boneIdx == -1)
        {
            continue;
        }
 
        FbxAMatrix kBindPose = pNode->EvaluateGlobalTransform();
        const QMatrix4x4& bindPose = TransformToQMatrix(kBindPose);
        const QMatrix4x4& invBindPose = bindPose.inverted();
        // 这里直接提取了bindPose并计算offset矩阵，可以取我们之前预先算好的
        for (FbxTime time = startTime; time <= endTime; time += timePerFrame, ++numFrames)
        {
            FbxAMatrix kMatGlobal = pNode->EvaluateGlobalTransform(time);//transform of bone in world space at time t
            const QMatrix4x4& globalMat = TransformToQMatrix(kMatGlobal);
 
            if(pCurrentAnimator)
            {
                QMatrix4x4 renderMatrix = globalMat * invBindPose;
                pCurrentAnimator->AddFrame(numFrames, boneIdx, renderMatrix);
            }
        }
    }
 
}
bool CFbxImporter::LoadFbx(const string& pFilename)
{
    vecNodes.clear();
    pCurrentAnimator = &CAnimationEngine::Inst()->CreateAnimator(pFilename);
    FbxImporter* lImporter = FbxImporter::Create(pManager,"");
 
    const bool lImportStatus = lImporter->Initialize(pFilename.c_str(), -1, pManager->GetIOSettings());
 
    if( !lImportStatus )
    {
        FbxString error = lImporter->GetStatus().GetErrorString();
        qDebug() << "load false";
        return false;
    }
 
    if(!lImporter->Import(pScene))
    {
        lImporter->Destroy();
        return false;
    }
 
    lImporter->Destroy();
 
    ProcessNode(pScene->GetRootNode());
    ProcessAnimation();
    return true;
}

首先出于简单考虑，我们的动作没有经过任何曲线压缩处理，而是直接存储了每帧的所有骨骼蒙皮数据。

蒙皮的这个过程可以在CPU或GPU中完成，不同引擎出于不同的考虑会有自己的实现。此处我们暂且放到GPU中实现，这意味着需要我们向顶点着色器发送所有骨骼的蒙皮变换矩阵，顶点根据自己相关的骨骼下标去查找对应的矩阵，并根据权重进行加权平均，得到最终的变换矩阵。

以上是最为普通的线性蒙皮计算方式。在实际应用中可能会遇到肩膀之类的地方变得非常细长的问题，此处需要我们考虑新的蒙皮计算来规避这一问题，具体内容可以尝试查阅相关资料。

void main()
{
    vec4 position = vec4(0,0,0,1);
 
    if(HasAnim)
    {
        ivec4 index = ivec4(a_boneindex);
        if(index.x < 0 || index.x >= BONE_NUM)
        {
            position = a_position;
        }
        else
        {
            mat4 BoneTransform;
            BoneTransform =   Bones[index.x] * a_boneweight.x;
            BoneTransform +=  Bones[index.y] * a_boneweight.y;
            BoneTransform +=  Bones[index.z] * a_boneweight.z;
            BoneTransform +=  Bones[index.w] * a_boneweight.w;
            position = BoneTransform * a_position;
            v_tangent = mat3(BoneTransform) * a_tangent;
            v_normal = mat3(BoneTransform) * a_normal; // 有非等比缩放需为逆转置
        }
    }
    else
    {
        position = a_position;
        v_normal = a_normal;
        v_tangent = a_tangent;
    }
 
    mat3 M1 = mat3(IT_ModelMatrix[0].xyz, IT_ModelMatrix[1].xyz, IT_ModelMatrix[2].xyz);
    v_normal = normalize(M1 * v_normal);
 
    mat3 M2 = mat3(ModelMatrix[0].xyz, ModelMatrix[1].xyz, ModelMatrix[2].xyz);
    v_tangent = normalize(M2 * v_tangent);
 
    v_texcoord = a_texcoord;
 
    gl_Position = ModelMatrix * position;
    gl_Position = ViewMatrix  * gl_Position;
    gl_Position = ProjectMatrix * gl_Position;
}
123456789101112131415161718192021222324252627282930313233343536373839404142

在GPU中，我们做一个简单的线性混合，此时还要注意处理法线相关数据的变换，避免做动作时光照错误。

动画管理实现

为了控制骨骼动画的播放、更新，我们需要一个动画管理的类。我们默认每个角色同一时间只能播放一个动作。

首先我们需要一个存储动画信息的类，也就是导入动画时我们用到的pCurrentAnimator变量。使用一个二维vector存储每帧每骨骼的蒙皮矩阵。

class CAnimator
{
private:
    vector<vector<QMatrix4x4>> m_vecFrames;
public:
    QMatrix4x4& GetFrame(int time, unsigned int boneIdx) { return m_vecFrames[time][boneIdx]; }
    vector<QMatrix4x4>& GetFrame(int time) { return m_vecFrames[time]; }
 
    void AddFrame(int time, unsigned int boneIdx, const QMatrix4x4& frame)
    {
        m_vecFrames[time][boneIdx] = frame;
    }
    int GetFrameNum() { return static_cast<int>(m_vecFrames.size()); }
    void Init(unsigned int boneNum, unsigned int frameNum);
};
123456789101112131415

其次，为了播放每个动作，我们需要定义我们如何播放这个动作，比如是否循环、混合时间、回调等信息。

struct SEvent
{
    float                           m_time = 0;
    float                           m_lastTime = 0;
    string                          m_path;
    bool                            m_bLoop = true;
    int                             m_blendFrame = 10;
    unordered_map<int,function<void()>> m_callbacks;
    vector<QMatrix4x4>              m_cachePose; // 和上一动作混合时，在事件中缓存一下上一动作姿态
    SEvent() { }
    SEvent(const string& path, bool bLoop, int blendFrame)
        : m_path(path), m_bLoop(bLoop), m_blendFrame(blendFrame) { }
};
12345678910111213

最后是我们最终控制播放的动画类，执行播放、更新、管理相关操作。

#define FRAME_PER_MS 0.03f
class CAnimationEngine
{
private:
    CAnimationEngine() {  }
 
    static CAnimationEngine* m_inst;
    vector<Bone*> m_bones;
    unordered_map<string, int> m_mapBonesIndex;
    unordered_map<string, CAnimator> m_animators;
    unordered_map<Object*, SEvent> m_events;
 
public:
 
    static CAnimationEngine* Inst()
    {
        if(!m_inst) m_inst = new CAnimationEngine();
        return m_inst;
    }
 
    void        Init();
    void        PlayAnimation(Object* obj, const string& path, bool bLoop = true, int blendFrame = 10);
    bool        UpdateAnimation(Object* obj, QOpenGLShaderProgram* program);
    bool        HasAnimator(const string& name) { return m_animators.find(name) != m_animators.end(); }
    CAnimator&  GetAnimator(const string& name) { return m_animators[name]; }
    CAnimator&  CreateAnimator(const string& name) { m_animators[name] = CAnimator(); return m_animators[name]; }
 
    void        AddBone(Bone* bone) {  m_bones.push_back(bone); m_mapBonesIndex[bone->m_name] = m_bones.size() - 1;}
    int         GetBoneNum() { return static_cast<int>(m_bones.size());}
    int         GetBoneIndex(const string& name) { return m_mapBonesIndex.find(name) == m_mapBonesIndex.end() ? -1 : m_mapBonesIndex[name]; }
    Bone*       GetBones(int i) { return m_bones[static_cast<size_t>(i)];}
};

播放动作实际上就是给角色添加新的动作指令，并且覆盖掉旧的数据。

void CAnimationEngine::PlayAnimation(Object* obj, const string& path, bool bLoop, int blendFrame)
{
    if(m_animators.find(path) == m_animators.end())
    {
        return;
    }
    int frame；
    string oldPath;
    if(m_events.find(obj) != m_events.end() && m_animators.find(m_events[obj].m_path) != m_animators.end())
    {
        SEvent& event = m_events[obj];
        if(blendFrame)
        {
            oldPath = event.m_path;
            frame = static_cast<int>(event.m_time * FRAME_PER_MS);
        }
    }
    m_events[obj] = SEvent(path, bLoop, blendFrame);
    if(!oldPath.empty())
    {
        CAnimator& animator = m_animators[oldPath];
        m_events[obj].m_cachePose = animator.GetFrame(frame);
    }
}

此处执行的更新动画也就是根据计算得到的当前帧数取得对应蒙皮矩阵，然后根据传入的shaderProgram向着色器发送蒙皮矩阵信息。我们采样的时间可能在两帧之间，此时可以选择在前后两帧之间按时间插值，也可以选择我这样偷懒的方法，在两者之间直接取一个作为结果：

bool CAnimationEngine::UpdateAnimation(Object* obj, QOpenGLShaderProgram* program)
{
    if(m_events.find(obj) == m_events.end())
    {
        return false;
    }
    SEvent& event = m_events[obj];
    if(m_animators.find(event.m_path) == m_animators.end())
    {
        return false;
    }
    CAnimator& animator = m_animators[event.m_path];
    if(animator.GetFrameNum() == 0)
    {
        return false;
    }
 
    float curTime = RenderCommon::Inst()->GetMsTime();
    event.m_time += curTime - event.m_lastTime;
    int frame = static_cast<int>(event.m_time * FRAME_PER_MS); // 30/1000帧每毫秒
 
    if(frame >= animator.GetFrameNum())
    {
        if(event.m_bLoop)
        {
            event.m_time = 0;
            frame = 0;
        }
        else
        {
            frame = animator.GetFrameNum() - 1;
        }
    }
 
    int size = static_cast<int>(m_bones.size());
    int location = program->uniformLocation("Bones");
    if(event.m_blendFrame > 0 && frame <= event.m_blendFrame && event.m_cachePose.size() > 0)
    {
        // ...
        // 此处做动作融合，代码略，主要是根据cachePose和currentPose先Decompose反解出位移旋转
        // 缩放信息，然后分别插值，计算得到新的变换矩阵
        // program->setUniformValueArray(location,final.data(), size);
 
    }
    else
    {
        program->setUniformValueArray(location,animator.GetFrame(frame).data(), size);
    }
 
    if(callbacks.find(frame) != callbacks.end())
    {
        callbacks[frame]();
    }
    event.m_lastTime = curTime;
 
    return true;
}

三维小球碰撞模拟

添加小球

首先，为了能够无限添加小球，采用链表结构，并定义小球结构体，其中包含小球的各个物理属性。

struct ball {

glm::vec3 position; //球心坐标

glm::vec3 speed; //速度矢量

glm::vec3 color；//可有可无

float r; //小球半径

float m; //小球质量

struct ball* next；

};

在渲染循环里面加上p->position += p->speed * deltaTime;实现小球移动。deltatime为渲染的时间间隔。

然后就是简单的循环，用来筛选发生碰撞的小球。

struct ball* p = head;

struct ball* q;

while (p != NULL) {

q = p->next;

while ((q != NULL)) {

if ((veclength(p->position - q->position) <= (p->r + q->r))) {

}

q = q->next;

}

p = p->next;

}

接着最关键的就是发生碰撞的两个小球的代码了。

碰撞

弹力是关键。小球碰撞速度的改变，终究还是因为他们之间的相互作用力，给了他们加速度。

于是我在小球属性里面添加了加速度glm::vec3 a，并且在渲染循环里面添加了p->speed += p->a * deltaTime;当小球发生碰撞时，根据质量反比，赋给他们分离的加速度99999.0f/m。

if ((veclength(p->position - q->position) <= (p->r + q->r))) {

p->a = glm::normalize(p->position - q->position) * 999999.0f / p->m;

q->a = glm::normalize(q->position - p->position) * 999999.0f / q->m;

}

此外，由于自然界不存在完全弹性碰撞，因此加了一个随速度阻尼，保证熵增。

最终代码如下：

void BALLMOVE() {

struct ball* p = head;

struct ball* q;

while (p != NULL) {

p->speed += p->a * deltaTime;

p->speed += 30.0f * deltaTime * glm::vec3(0, -1, 0);//这是重力加速度

p->position += p->speed * deltaTime;

q = p->next;

p->a = glm::vec3(0, 0, 0);

while ((q != NULL)) {

if ((veclength(p->position - q->position) <= (p->r + q->r))) {

float k = fabs(veclength(p->position - q->position) - (p->r + q->r));

float kn = pow(7,k*k)-1;

p->a += kn*(glm::normalize(p->position - q->position) * 999999.0f - (p->speed - q->speed) * 100000.0f) / p->m;

q->a += kn*(glm::normalize(q->position - p->position) * 999999.0f - (q->speed - p->speed)*100000.0f) / q->m;

}

q = q->next;

}

p = p->next;

}

}

将这个函数插入到渲染循环里面，就可以实现小球碰撞。

质点弹簧系统仿真

质点弹簧模型

弹簧质点模型是利用牛顿运动定律来模拟物体变形的方法，把薄膜看作一张由质点排列组成的网格，质点的位置代表膜上一点的空间位置，质点之间用无质量的、自然长度大于零的弹簧连接。质点弹簧模型中，弹簧的类型有以下三种：

结构弹簧（Structure Spring）：连接紧密相邻的横向和纵向质点，提供薄膜的拉伸刚度，固定整体结构；

剪切弹簧（Shear Spring）：连接在同一对角线上的相邻质点，模拟薄膜倾斜方向的作用力，防止剪切变形错位等；

弯曲弹簧（Bending Spring）：连接横向和纵向隔一个质点的两质点，提供弯曲刚度，薄膜折叠时弯折边缘圆滑。

具体如下图所示

一个常见的模拟就是布料模拟，接下来做一个简单的介绍，具体实现不做展开

我们将布料模型构建为一个包含n个顶点的三角形网格模型。对于顶点来说，其位置为

那么整个布料系统的几何状态向量

同样对于力

i点处的纹理坐标为（ui,vi）

计算时，需要将上述所有力汇集到力的状态向量f中。可知对于质点i，x=fi/mi

mi为质点的质量。为所有包含该质点三角形质量的三分之一。定义一个质量对角矩阵

M∈R3n×3n为

diagM=（m1,m1,m1,m2,m2,m2,…,mn,mn,mn），则对于系统的加速度状态可表达为

由上文可知，对布料系统的模拟就是对该微分方程的求解。将该微分方程转换为2阶微分系统：

然而此法需要步长h越小越好，否则可能出现刚度（stiff）问题，导致计算的结果不收敛，如下图所示

故我们使用过另一种反向欧拉法来计算，令

欧拉法的计算只依赖于t0时的状态，而反向欧拉法则需要额外依赖终点状态，此方程为非线性方程，无法直接求解，可以使用泰勒展开来求近似值，令

后两项代表的是整个系统对质点x作用的力，则

力的分析

由上文可知，施加在布料上的力为

给定内力条件向量C(x)，系统的能量定义为

k为该系统刚度（stiffness），质点i受到的内力为

定义微分矩阵

可知K是一个对称矩阵。由于内列条件C(x)不依赖于速度v，即

拉伸力

由于网格顶点有位置属性xi, 纹理坐标（ui,vi）。令函数w(u,v)为将纹理坐标转换为世界坐标的转换函数。布料的在坐标uv处的拉伸力可表达为偏导数

若布料处于松弛状态，wu=wv=1，可以用于测量沿u,v方向上的拉伸力。可以基于三角形来对拉伸力进行分析，对于三角形xi,xj,xk来说，定义

a为三角面的面积，bu，bv一般设置 为1

剪切力

阻尼力

算法简介

简单介绍一下算法部分cpp的思路

首先初始化部分需要完成 Key tempEdge

调用gpu

在例如uploadToGPU()中完成 create vbos以及create indexBuffer

同理需要更新位置以及法向量的缓存vbo，即 update vbo

对于力的系统来说，同样需要进行更新

首先需要进行check constrains

接下来依次按照文章介绍计算 stretch force，shear forces，shear damping，damper force

计算完成之后则判断是否添加这个力，例如使用一个图形化界面中的按钮来控制，则可以如下步骤

add external force，add damping force

添加之后更新速度和位置

update velocity，update position

除此之外或许还需要用到碰撞系统

collision

最后不要忘记标准化法向量并且更新

normalize normal，update normals

在此概括一下

模拟开始时，每个弹簧都处于原始状态。根据胡克定律，弹簧的形变大小正比与施加的力。阻尼器可以施加一个与弹簧形变方向相反的力，可以使弹簧的震荡逐渐衰减。

粒子和弹簧的力可以对时间进行积分，以改变每个粒子的状态。对整个对力和位置的计算过程不断迭代，以实现布料的动画模拟。

水波纹模拟

glsl水波纹效果，可使用shadertoy直接运行。sin和iTime配合得到水波纹mask，通过mask流动变化得到水波纹效果。

脚本1

#iChannel0 "file://./bg0.jpg"

// 水波纹中心点

const vec2 center = vec2(0.5, 0.5);

// 水波纹速度

const float speed = 2.5;

// 水波纹强度

const float intensity = 36.0;

void mainImage(out vec4 fragColor, in vec2 fragCoord) {

// 计算uv

vec2 uv = fragCoord.xy / iResolution.xy;

// 计算center到各个uv的值

float dst = distance(uv, center);

// fragColor = vec4(vec3(dst), 1.0);

// 计算水波纹

// sin将对这个范围限制到[-1, 1]

// dst范围是[0, sqrt(2)]

// intensity可以理解为波纹数

// iTime的引入让一个点可以从[-1, 1]的变化

// speed是速度

// 可以理解为dst*intensity是初始波纹, iTime的引入让波纹动了起来

float wave = sin(sqrt(dst)*intensity-iTime*speed);

// fragColor = vec4(vec3(wave), 1.0);

// 接下来使用da和偏移的db来mix, 也就是水波纹白色区域取db, 黑色区域取da

vec3 da = texture2D(iChannel0, uv).rgb;

vec3 db = texture2D(iChannel0, uv-vec2(0.01)).rgb;

vec3 res = mix(da, db, wave*wave); // wave*wave范围为[0, 1]

fragColor = vec4(res, 1.0);

}

脚本2

#iChannel1 "file://./bg0.jpg"

#iChannel2 "file://./mask.png"

const float intensity = 0.02;

const float speed = 10.0;

void mainImage(out vec4 fragColor, in vec2 fragCoord) {

vec2 uv1 = fragCoord.xy / iResolution.xy;

// mask的uv流动

vec2 uv2 = uv1 + vec2(0.0, iTime * speed / 100.0);

// 由于vs中不能对iChannel2进行repeat, 所以进行一个求余操作

uv2 = mod(uv2, 1.0);

// 得到mask的流动

vec2 wave = texture2D(iChannel2, uv2).xy;

// wave范围是[0, 1], intensity限制了范围, 错位得取像素

uv1 = uv1 + wave*intensity;

fragColor = texture(iChannel1, uv1);

}

脚本3

#iChannel0 "file://./bg0.jpg"

const float speed = 1.0;

const float intensity = 6.0;

float wave(vec2 pos, float t, float freq, float numWaves, vec2 center) {

float d = length(pos - center);

d = log(1.0 + exp(d));

return 1.0/(1.0+20.0*d*d) * sin(2.0*3.1415*(-numWaves*d + t*freq));

}

float height(vec2 pos, float t) {

float w;

w = wave(pos, t, speed, intensity, vec2(0.5, -0.5));

w += wave(pos, t, speed, intensity, -vec2(0.5, -0.5));

return w;

}

vec2 normal(vec2 pos, float t) {

return vec2(height(pos - vec2(0.01, 0), t) - height(pos, t),

height(pos - vec2(0, 0.01), t) - height(pos, t));

}

void mainImage(out vec4 fragColor, in vec2 fragCoord) {

vec2 uv = fragCoord.xy / iResolution.xy;

vec2 uvn = 2.0*uv - vec2(1.0);

uv += normal(uvn, iTime);

fragColor = texture2D(iChannel0, uv);

}

天体模拟

一些概念

OpenGL在光照模拟中将光线分为辐射光、环境光、漫反射光和镜面反射光4 种独立的成分。

OpenGL中，材料的颜色决定于其对入射光中的红、绿、蓝各成分的反射比例。

材料的设置参数中也包括对R、G、B值的设定，但与光照的参数设置相比，两者的含义是不同的。对材料来说， R、G、B值分别对应材料对各种颜色光的反射比例。

与光线的情况相似，材料也具有辐射色、环境色、漫反射色和镜面反射色。为了模拟场景中的发光体，可以设定材料的辐射光特性；而环境色、漫反射色和镜面反射色则反映材料对环境光、漫反射光和镜面反射光的反射系数。

OpenGl搭建

VS中搭建OpenGL所需的头文件及库文件：

云盘链接：https://pan.baidu.com/s/1fFsic53et67IjIHdV-eUvg 密码：unyr

搭建OpenGL所需文件的方法：

① 找到并打开Visual Studio的安装目录，比如：

C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Tools\MSVC\14.14.26428

② 打开\include文件夹，把资源里的“GL”文件夹整个复制到里面。

③ 回到①步骤，再打开lib\x86，把资源“LIB”文件夹里的5个.lib文件复制到lib文件夹里面。

④ 最后把资源“C盘SysWOW64或system32目录下的文件”文件夹里的两个文件复制到以下目录：

C:\Windows\system32文件夹内（如果你的电脑是32位系统请选这个）

或‪C:\Windows\SysWOW64(64位系统请选这个）

使用vs插件：

① 文件 -- 新建 -- 项目 -- Visual C++ -- 空项目；

② 项目 -- 管理Nuget程序包 -- 浏览 -- (搜索栏输入)NupenGL；

③ 选择第一个，安装即可完成配置。

原文链接：https://blog.csdn.net/weixin_41918712/article/details/82502347

天体源码

#include <GL/glut.h>

#include <stdlib.h>

static int year = 0, day = 0, moon = 0;

void init(void)

{

GLfloat sun_light_position[] = { 1.0f, 1.0f, 1.0f, 0.0f }; //表示光源所在的位置,（X, Y, Z, W）

GLfloat sun_light_ambient[] = { 0.0f, 0.0f, 0.0f, 1.0f }; //环境光分量

GLfloat sun_light_diffuse[] = { 1.0f, 1.0f, 1.0f, 1.0f }; //漫反射分量

GLfloat sun_light_specular[] = { 1.0f, 1.0f, 1.0f, 1.0f }; //镜面反射分量

glLightfv(GL_LIGHT0, GL_POSITION, sun_light_position); //指定第0号光源的位置

glLightfv(GL_LIGHT0, GL_AMBIENT, sun_light_ambient); //GL_AMBIENT表示各种光线照射到该材质上，

//经过很多次反射后最终遗留在环境中的光线强度（颜色）

glLightfv(GL_LIGHT0, GL_DIFFUSE, sun_light_diffuse); //漫反射后

glLightfv(GL_LIGHT0, GL_SPECULAR, sun_light_specular);//镜面反射后

glEnable(GL_LIGHT0); // 激活指定光源，使用第0号光照

glEnable(GL_LIGHTING); //激活光照功能

glEnable(GL_DEPTH_TEST); //这句是启用深度测试，这样，在后面的物体会被挡着

glClearColor(0.0, 0.0, 0.0, 0.0);

glShadeModel(GL_SMOOTH); //设置光滑着色模型

}

void display(void) {

glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //清空颜色和深度缓冲

glColor3f(1.0, 1.0, 1.0);

// 定义太阳的材质并绘制太阳

GLfloat sun_mat_ambient[] = { 0.5f, 0.5f, 0.0f, 1.0f }; //材料环境光颜色

GLfloat sun_mat_diffuse[] = { 0.5f, 0.0f, 0.5f, 1.0f }; //材料漫反射光颜色

GLfloat sun_mat_specular[] = { 0.0f, 0.0f, 0.0f, 1.0f }; //材料镜面反射光颜色

GLfloat sun_mat_emission[] = { 0.0f, 0.0f, 0.0f, 1.0f }; //材料辐射光颜色

GLfloat sun_mat_shininess = 60.0f; // 材料光滑指数

glMaterialfv(GL_FRONT, GL_AMBIENT, sun_mat_ambient); //定义材料的前面采用环境光颜色

glMaterialfv(GL_FRONT, GL_DIFFUSE, sun_mat_diffuse); //材料的前面为漫反射光颜色

glMaterialfv(GL_FRONT, GL_SPECULAR, sun_mat_specular); //定义材料的前面为镜面反射光颜色

glMaterialfv(GL_FRONT, GL_EMISSION, sun_mat_emission); //定义材料的前面辐射光颜色

glMaterialf(GL_FRONT, GL_SHININESS, sun_mat_shininess); //材料的前面的光滑指数

glPushMatrix(); //a 此处的push是为了表明堆栈当前状态

glutSolidSphere(1.0, 20, 16); /* draw sun */ //参数分别是半径、以Z轴上线段为直径分布的圆周线的条数（类似经线）、围绕在Z轴周围的线的条数（类似纬线）

// 定义地球的材质并绘制地球

GLfloat earth_mat_ambient[] = { 0.0f, 0.0f, 0.5f, 1.0f }; //材料环境光颜色，偏蓝色

GLfloat earth_mat_diffuse[] = { 0.0f, 0.0f, 0.5f, 1.0f }; //材料漫反射光为偏蓝色的光源

GLfloat earth_mat_specular[] = { 0.0f, 0.0f, 1.0f, 1.0f }; //材料镜面反射光为蓝色的光源

GLfloat earth_mat_emission[] = { 0.0f, 0.0f, 0.0f, 1.0f }; //材料辐射光颜色，白光

GLfloat earth_mat_shininess = 30.0f;

glMaterialfv(GL_FRONT, GL_AMBIENT, earth_mat_ambient);

glMaterialfv(GL_FRONT, GL_DIFFUSE, earth_mat_diffuse);

glMaterialfv(GL_FRONT, GL_SPECULAR, earth_mat_specular);

glMaterialfv(GL_FRONT, GL_EMISSION, earth_mat_emission);

glMaterialf(GL_FRONT, GL_SHININESS, earth_mat_shininess);

glRotatef((GLfloat)year, 0.0, 1.0, 0.0); //对局部坐标系统进行旋转（局部坐标系统最初与全局坐标系统是一致的），让地球绕太阳旋转

glTranslatef(2.0, 0.0, 0.0); //把局部坐标系统平移到地球轨道上的一个位置

glRotatef((GLfloat)day, 0.0, 1.0, 0.0); //使局部坐标轴进行旋转，让地球绕局部坐标系统y轴旋转(即绕自身旋转)

glutSolidSphere(0.3, 10, 8); /* draw smaller planet */

//自己

glPushMatrix();

GLfloat moon_mat_ambient[] = { 0.0f, 0.5f, 0.5f, 1.0f }; //材料环境光颜色

GLfloat moon_mat_diffuse[] = { 0.0f, 0.5f, 0.5f, 1.0f }; //材料漫反射光

GLfloat moon_mat_specular[] = { 0.0f, 1.0f, 1.0f, 1.0f }; //材料镜面反射光

GLfloat moon_mat_emission[] = { 0.0f, 0.0f, 0.0f, 1.0f }; //材料辐射光颜色，白光

GLfloat moon_mat_shininess = 15.0f;

glMaterialfv(GL_FRONT, GL_AMBIENT, moon_mat_ambient);

glMaterialfv(GL_FRONT, GL_DIFFUSE, moon_mat_diffuse);

glMaterialfv(GL_FRONT, GL_SPECULAR, moon_mat_specular);

glMaterialfv(GL_FRONT, GL_EMISSION, moon_mat_emission);

glMaterialf(GL_FRONT, GL_SHININESS, moon_mat_shininess);

glRotatef((GLfloat)day, 0.0, 1.0, 0.0); //对局部坐标系统进行旋转（局部坐标系统最初与全局坐标系统是一致的），让地球绕太阳旋转

glTranslatef(0.5, 0.0, 0.0); //把局部坐标系统平移到地球轨道上的一个位置

glRotatef((GLfloat)day, 0.0, 1.0, 0.0); //使局部坐标轴进行旋转，让地球绕局部坐标系统y轴旋转(即绕自身旋转)

glutSolidSphere(0.1, 8, 7); /* draw smaller planet */

glPopMatrix();

//自己

//glPopMatrix(); // 这个pop使得堆栈状态回到b状态

glPopMatrix(); // 这个pop使得堆栈状态回到a状

glutSwapBuffers(); //交换前台和后台缓冲区

}

void reshape(int w, int h) {

glViewport(0, 0, (GLsizei)w, (GLsizei)h);

glMatrixMode(GL_PROJECTION); //设置投影变换

glLoadIdentity();

gluPerspective(60.0, (GLfloat)w / (GLfloat)h, 1.0, 20.0); //设置透视投影，视角大小及宽高比来确定沿视线方向的棱锥，并通过指定近、远剪切面与视点间的距离来截断棱锥，得到取景体。

// glOrtho(-3.0, 3.0, -3.0, 3.0, -10, 10); //设置正交投影，取景体是一个各面均为矩形的六面体。前4个参数是坐标，后两个是距离

glMatrixMode(GL_MODELVIEW); //设置视图变换

glLoadIdentity();

gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);

}

void keyboard(unsigned char key, int x, int y)

{

switch (key)

{

case 'd': //逆时针自转

case 'D':

day = (day + 10) % 360;

moon = (moon + 5) % 360;

glutPostRedisplay(); //标记当前窗口需要重新绘制

break;

case 'a': //顺时针自转

case 'A':

day = (day - 10) % 360;

glutPostRedisplay();

break;

case 'w': //逆时针公转

case 'W':

year = (year + 5) % 360;

day = (day + 10) % 360;

moon = (moon + 5) % 360;

glutPostRedisplay();

break;

case 's': //顺时针公转

case 'S':

year = (year - 5) % 360;

moon = (moon - 10) % 360;

glutPostRedisplay();

break;

case 27:

exit(0);

break;

default:

break;

}

}

int main(int argc, char** argv)

{

glutInit(&argc, argv);

glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);

glutInitWindowSize(800, 600);

glutInitWindowPosition(100, 100);

glutCreateWindow("Sun & Earth");

init();

glutDisplayFunc(display); //设置display函数，当需要进行画图时，这个函数就会被调用

glutReshapeFunc(reshape); //调用此函数，重新建立用作新渲染画布的矩形区域

glutKeyboardFunc(keyboard); //这个函数是告诉窗口系统，哪一个函数将会被调用来处理普通按键消息

glutMainLoop(); //进入 GLUT 消息循环，开始执行程序，显示窗口，并且等待窗口关闭才会返回

return 0;

}

效果图

太阳&地球&月亮 如下图所示：

功能描述

键盘A & a：地球顺时针自转； 键盘D & d：地球逆时针自转；

键盘S & s：地球顺时针公转； 键盘W & w：地球逆时针公转；