随机生成植物生长及舞动算法

几年前写过一套随机树木的生成算法，其中使用了分形和放样建模。那时候还不知道有speedtree这款软件，写的比较粗糙，最近看了speedtree的演示把原算法改进了一下，算是一个speedtree的简化版本。
重构主要是把原先使用递归函数实现的分形过程改为了非递归的，这样结构更清晰一些，尽量把树木生成分拆为若干独立步骤，减少每个步骤之间的耦合。使用骨骼动画模拟树木生长及随风舞动。

搜集的一些关于树木形态的数学规律：
分形：每条分枝就像一颗经过旋转的独立树。
达芬奇公式：当树木分岔时，树干的粗度等于同一高度树枝的总粗度，工程学的角度发现这样的结构最能抵抗强风。
一般树冠分为六种形态：长型、球型、椭圆型、瓶型，金字塔型和垂枝型。这跟顶端优势是否明显有关系，如果顶端优势明显，向两边生长就会抑制，形成高大的树，如松树；如果顶端优势不明显，树就会倾向于横向生长，如橡树。树冠的形态和生长地域息息相关，寒冷地区的针叶树，树冠比温带地区要小，而且呈金字塔状，这样的结构可以防止在下雪的时候侧枝积累过多的重量而折断。
分叉很多时候符合斐波纳契数列（Fibonacci Sequence），形如[1、1、2、3、5、8、13、21]，F0=0，F1=1，Fn=F(n-1)+F(n-2)。越往后面，相邻两个斐波纳契数的比（1/2， 2/3， 3/5， 5/8， 8/13）越近似等于0.618黄金分割。在发育成树枝的叶腋的排列中，记任意一个点为0然后往上数，直到另一个相同位置出现叶腋，其中的周数/叶腋数（1/2， 1/3， 2/5， 3/8， 8/13）很多都是斐波那契数的相邻两个。这种排列方式可以保证在空间上错落开，充分利用阳光。分叉黄金角为137.51度。

可调整分支参数：
       BranchType branchType; //分支类型
       BranchMeshType meshType; //模型类型
       RectF texRect ; //纹理
       int branchNum ; //当前层枝条数
       float length ; //枝干长度
       float radiusBottom ; //枝干底部半径
       float radiusTop ; //枝干顶部半径
       float sizeVariation ; // 大小异性
       float bend ; //枝干弯曲度
       float gravityW ; //全局重力（world）
       float gravityL ; //局部重力（local 平行父枝干方向）
       float gravityC ; //局部向心力（center垂直指向父枝干方向）适用于花瓣内收
       int pathStepNum ; //长度段数
       int shapeStepNum ; //切面段数
       bool axisXAlign ; //是否将x轴放到父枝的垂直面上，树叶(花瓣)的放样才正确
       float pitchAng ; //当前层枝条与父枝的夹角 pitch
       float pitchAngVariation; //pitchAng差异范围
       float rotAng ; //当前层枝条与父枝的旋转递增夹角 rot 一般为黄金角137.51，君子兰对称为180 十字花90
       //float rotAngVariation ;

       float minOffset ; //当前层枝条在父枝上分叉点高度偏移
       float maxOffset ; //.......最大
       float offsetRelax ; //分叉均匀度，0度时-分叉点集中到每节的起始处

       //叶子参数
       float leafConvex ; //叶子横向凹凸 (放样形状)
       String leafMesh ;

//风力动画
float windWeight ; //风动系数

       //生长动画
       float matureTime ; //成熟时间[0~1]，成熟后大小不变，成熟前线性长大并且被父枝的生长大小叠加缩放
       //生长期纹理融合 in treedef

参数调整技巧：
垂柳：设置全局重力参数gravityW
花：花瓣可以使用多层树叶代替，顶层和底层花瓣向心力不同、局部重力不同，模拟出花瓣的不同弧度。
球形树：设置分枝pitchAng=0、pitchAngVariation=180，在分枝的末端添加树叶pitchAng=90、pitchAngVariation=0。
草丛：可以使用一个微小的假主干，或者将主干数目设置为大于1。
树叶不一定是末端，也可以再分出次级树枝、树叶或花。

算法大部分代码集中在分枝路径的计算和骨骼动画模型的生成，大体的步骤：
1、计算枝条个数（因为不想用动态数组，所以先计算个数分配空间）
2、构造枝条结构
3、计算骨骼个数
4、构造骨骼结构
5、计算枝条曲线、放样形状。这一步代码较多，路径受弯曲度影响分形叠加曲线，路径受重力、局部重力、向心力等影响。放样形状可以是闭合曲线（枝条）或弯曲线（面片，十字面片）。
6、计算骨骼矩阵
7、生成骨骼动画驱动的树
其中只有枝条计算和骨骼计算有一些穿插耦合，因为枝条曲线的计算需要用到父骨骼的矩阵，骨骼矩阵又需要根据自身枝条曲线计算，不好完全独立开。

数据结构：


enum BranchType
{BT_Trunk , //树干BT_Branch, //树枝BT_Leaves, //树叶BT_Root  , //树根BT_Flower, //花、果实
};
const char* BranchTypeToString(int enumeration);
bool StringToBranchType(const char* str,int& type);enum BranchMeshType   //模型类型
{MT_Strip       , //线放样MT_Cross       , //十字放样MT_Cylinder    , //圆放样MT_CylinderCap , //圆放样封顶MT_Billboard   , //MT_Sphere      , //MT_Mesh        , //
};
const char* BranchMeshTypeToString(int enumeration);
bool StringToBranchMeshType(const char* str,int& type);//随机树配置
class TreeDef
{
public:
#define  MaxDepth        4
#define  MaxShapeStep    16
#define  MaxPathStep     8
#define  MaxSubBranch    50	#define  MaxSubBranchDef 3
#define  MaxBranchDef    32 //128#define  MaxGrowTex      4enum Method{TT_StaticAndUV,  //静态树 纹理动画TT_AnimBone   ,  //绑定骨骼动画};//每一级分枝class BranchDef{public:char  name[64]         ;int   ID               ;int   parentID         ;int   depth            ;  //层深度//int   seed             ;  //单独seedBranchType     branchType;  //分支类型BranchMeshType meshType; RectF texRect          ;  //纹理int   branchNum        ;  //当前层枝条数//todo枝长曲线  lengthBottom lengthCen lengthTopfloat length           ;  //枝干长度float radiusBottom     ;  //枝干底部半径float radiusTop        ;  //枝干顶部半径float sizeVariation    ;  //   float bend             ;  //枝干弯曲度（1-笔直度） 弯曲度和俯仰角效果有重叠？ 点在和路径垂直的横截面上偏移 起点无偏移 末端在横截面的圆环上float gravityW         ;  //全局重力（world） 路径受重力的影响float gravityL         ;  //局部重力（local平行父枝干方向）float gravityC         ;  //局部向心力（center垂直指向父枝干方向） 适用于花瓣内收int   pathStepNum      ;  //长度段数int   shapeStepNum     ;  //切面段数bool  axisXAlign       ;  //将x轴放到父枝的垂直面上，树叶(花瓣)的放样才正确//float pitchAng         ;  //当前层枝条与父枝的夹角 pitchfloat pitchAngVariation;  //pitchAng差异范围 float rotAng           ;  //当前层枝条与父枝的旋转递增夹角 rot 一般为黄金角137.51，君子兰对称为180 十字花90//float rotAngVariation  ; float minOffset        ;  //当前层枝条在父枝上分叉点高度偏移 float maxOffset        ;  float offsetRelax      ;  //分叉均匀度，0度时-分叉点集中到每节的起始处//叶子参数float    leafConvex    ; //叶子横向凹凸 (放样形状)String   leafMesh      ;//风力动画float windWeight       ; //风动系数float windWeightSq     ; //修正风力权重呈平方增长，且避免段数越大弯曲累积越多//生长动画float matureTime       ; //成熟时间[0~1]，成熟后大小不变，成熟前线性长大并且被父枝的生长大小叠加缩放//生长期纹理融合 in treedefBranchDef* parent;int        subDefNum;BranchDef* subDef[MaxSubBranchDef];int   numW;RectF uiRect; //for diagram edit};TreeDef();void RandGen();void RandChange();void Load(const char* fileName);void Save(const char* fileName);BranchDef* GetBranchDef(int ID);void GenTopology();void UpdateAllItems(ListItemUpdater& updater);public:String m_fileName;String m_texture ;//生长期纹理融合 int    m_growTexNum;String m_growTexture[MaxGrowTex];float  m_growTexTime[MaxGrowTex]; //在此生长时间点开始改变材质Method m_method  ;int    m_seed    ;  BranchDef* m_trunkDef;  //主干，规定只有一个trunkdef，其它未连接parentID的分枝不被构造。int        m_branchDefNum;BranchDef  m_branchDef[MaxBranchDef];//mesh//RendSys::MovieClip* m_movieLib;
};//随机树模型生成及动画更新器 
class TreeProducer
{
public:class TBone;class Branch;TreeProducer();~TreeProducer();void Free();void EditRender();void GenRandTree(Mesh* mesh,Skeleton* skeleton,TreeDef* def,int randSeed=-1);//private://生成静态树：树叶和树枝纹理drawcall分开，树叶可以单独uv动画void GenRandTreeStaticAndUV(Mesh* mesh,TreeDef* def,int randSeed=-1);	//生成骨骼动画驱动的树：骨骼动画为了加速后期可以塌陷到morph，当skeleton为空时生成静态树void GenRandTreeAnimBone (Mesh* mesh,Skeleton* skeleton,TreeDef* def,int randSeed=-1);	void CalTBoneMatrix      (Branch* branch);void GenRendMeshSkeleton (Mesh* mesh,Skeleton* skeleton);	void RefreshRendSkeleton   (Skeleton* skeleton,int frame,int offset);void RefreshRendMixSkeleton(MixSkeleton* mixSkeleton,int frame,int offset);//生长 0~1void RefreshGrow(float life);//风的角度和风力，如果叠加生长动画，则RefreshGrow调用要在RefreshWind前面void RefreshWind(TreeWind* wind);class Branch{public :Branch();~Branch();public:int   indexG;  //在全局分枝中的索引int   indexL;  //在父枝分枝中的索引int   indexLT; //在父枝同类型分枝中的索引TreeDef::BranchDef* def; //TBone* bones; //arr[path step num+1]float lengthScale  ;  //枝干长度半径缩放系数float stemOffset   ;  //在父枝干上的分叉位置[0~parentDef->pathStepNum]float stemRot      ;  //分支方向float stemAngle    ;  //pitch//枝成熟期原始造型vec3  shapePos   [MaxShapeStep+1] ;   //放样形状坐标    vec3  shapeNormal[MaxShapeStep+1] ;   //放样形状法线  vec3  pathPos    [MaxPathStep+1 ] ;   //放样路径坐标    int      subBranchNum; //multi typeBranch** subBranchs;Branch*  parent;TexturePtr m_texture;};class TBone{public :TBone();~TBone();public:int     index      ; //indexWBranch* branch     ; //TBone*  parent     ;int     subBoneNum ;TBone** subBones   ;	mat4    matL       ; //原始局部矩阵mat4    matLGrow   ; //经过生长动画模拟后的局部矩阵，不叠加风力，若未执行生长模拟则等于matLmat4    matLWind   ; //经过生长动画模拟+风力后的局部矩阵，使用matLGrow作为原始局部矩阵叠加风力，若未执行生长模拟则等于matLGrowmat4    matW       ; //世界矩阵vec3    axisYW     ;风动//vec3  windRotSpeed; //角速度 此处使用简化的快速模拟，真实的物理模拟对每根骨骼创建刚体速度太慢//生长float   matureTime ; //成熟时间[0~1]，大小随时间呈ln曲线变化，成熟时达到0.8f，成熟后大小微变，成熟前线性长大并且被父枝的生长大小叠加缩放float   scaleW     ;};//private:
public :int      m_seed;Branch*  m_branchs  ; //主干chunk可以有多个int      m_branchNum;	TBone*   m_bones;int      m_boneNum;
};

生长动画算法：

很短也许不是最好的，但是比较快速，主要是逐级改变骨骼缩放。

void TreeProducer::RefreshGrow(float life)
{Clamp(life,0.0f,1.0f);mat4  matScale;float scaleL;float scaleL2; //y轴非线性缩放float matureTime;TBone* bone = m_bones;for(int i = 0; i < m_boneNum; ++i,++bone){if(bone->parent){matureTime = bone->branch->def->matureTime;//scaleW-time 为过三点的折线if(matureTime==1)bone->scaleW = 0.8f * life / matureTime;else if(life >= matureTime || matureTime==0)bone->scaleW = 0.8f + 0.2f*(life-matureTime)/(1-matureTime);else bone->scaleW = 0.8f * life / matureTime;if (bone->parent->scaleW<_EPSILON){scaleL  = 0;scaleL2 = 0;}else{scaleL  = bone->scaleW/bone->parent->scaleW;scaleL2 = bone->scaleW/sqrt(bone->parent->scaleW);}//优化展开matScale.FromScale(scaleL,scaleL2,scaleL);bone->matLGrow = matScale * bone->matL; bone->matW = bone->parent->matW * bone->matLGrow;}else{bone->scaleW = 1;}}
}

随风舞动动画模拟：

此处使用简化的快速模拟。真实的物理模拟需要对每根骨骼创建刚体，模拟风的作用力，速度太慢。如果同时开启生长模拟和风力模拟，则风力模拟需要放在生长模拟之后计算。

//简化的快速风动模拟
class TreeWind
{
public:TreeWind();void Random();void Update();const vec3& GetWindDir();float GetWindForce();//参数float strength  ; //强度float turbulence; //湍流强度float frequency ; //湍流频率  周期=TWOPI/frequencyprivate://输出vec3    m_tree3DWindDir;float   m_tree3DWindForce;vec3    m_tree3DWindAngSpeed;
};void TreeWind::Update()
{float acumTime = G_Timer->GetAccumTime();float time     = G_Timer->GetStepTimeLimited();m_tree3DWindForce = (sin(acumTime*frequency)*0.7f+sin(acumTime*frequency*2.337f+0.123f)*0.2f//+sin(acumTime*frequency*5.537f+0.823f)*0.1f)*strength; //周期=TWOPI/frequencyvec3 angAcc;RandDir(angAcc);m_tree3DWindAngSpeed += angAcc*0.2f*time;m_tree3DWindAngSpeed.Normalize();mat4 mat;mat.FromAxisAngle(m_tree3DWindAngSpeed,RandRange(60.0f,100.0f)*turbulence*DEG2RAD*time);m_tree3DWindDir = mat*m_tree3DWindDir;
}void TreeProducer::RefreshWind(TreeWind* wind)
{const vec3& windDir = wind->GetWindDir();float windForce = wind->GetWindForce();float windWeight;vec3  axisY(0,1,0);mat4  rot;vec3  axis;TBone* bone = m_bones;TreeDef::BranchDef* branchDef;for(int i = 0; i < m_boneNum; ++i,++bone){branchDef  = bone->branch->def;if(bone->parent){windWeight = branchDef->windWeightSq;if (windWeight>_EPSILON){CrossVec3(axis,bone->axisYW,windDir);float angRad = (windForce+0.4f) * windWeight * LengthVec3(axis); //LengthVec3(axis)减少抖动rot.FromAxisAngle(axis,angRad);bone->matLWind = bone->matLGrow * rot;}else{bone->matLWind = bone->matLGrow;}bone->matW = bone->parent->matW * bone->matLWind;bone->axisYW = axisY;bone->matW.AffectNormal(bone->axisYW);}}
}