海飞丝头发的研究和实现

先展示成果：

https://v.qq.com/x/page/z0853wtw21a.html

很早就看了milo大佬的爱丽丝的海飞丝，那个时候完全看不懂。最近再去看了一遍，还是看不懂。打算自己研究看看。

首先查了gpugems2里面关于海飞丝的简要实现方案：



打算就先从这个基本流程入手。

首先是准备了一个头的模型，和主要的头发束模型，然后根据头皮的顶点，去生成头发。头发直接用9个顶点连成一条线。先在编辑器中把这些线画出来，方便查看效果。

private void OnDrawGizmos()
{
    if (!DebugDraw || GetVertices() == null || !ValidateImpl(false))
        return;

    var scalpToWorld = ScalpProvider.ToWorldMatrix;
    var vertices = GetVertices();

    for (var i = 1; i < vertices.Count; i++)
    {
        if (i % Segments == 0)
            continue;

        var vertex1 = scalpToWorld.MultiplyPoint3x4(vertices[i - 1]);
        var vertex2 = scalpToWorld.MultiplyPoint3x4(vertices);

        Gizmos.DrawLine(vertex1, vertex2);
    }

    var worldBounds = GetBounds();
    Gizmos.DrawWireCube(worldBounds.center, worldBounds.size);
}
大概是这样子的：



按照gpugems2里面的描述，我们要把每一个顶点想象成一个小珍珠，用Verlet积分来计算粒子的运动，并且增加头发顶点之间的约束，即头发长度会倾向于保持恒定避免拉伸。

头发顶点部分就是对应的网格顶点，而头发的连线部分则需要单独定义。

private void UpdateBodies(RenderParticle[] renderParticles)
{
    var renderSettings = settings.RenderSettings;
    var sizeY = settings.StandsSettings.Provider.GetSegmentsNum();
    //sizeY是单根头发的段数，i是第几根头发，y是这个头发里的第几段,t是百分比
    for (var i = 0; i < renderParticles.Length; i++)
    {
         var x = i / sizeY;
         var y = i % sizeY;
         var t = (float)y / sizeY;
         //下面是渲染的三个参数，后面渲染部分在解释
         var data = new RenderParticle
         {
             Color = ColorToVector(renderSettings.ColorProvider.GetColor(settings, x, y, sizeY)),
             Interpolation = Mathf.Clamp01(renderSettings.InterpolationCurve.Evaluate(t)),
             WavinessScale = Mathf.Clamp01(renderSettings.WavinessScaleCurve.Evaluate(t)) * renderSettings.WavinessScale,
             WavinessFrequency = Mathf.Clamp01(renderSettings.WavinessFrequencyCurve.Evaluate(t)) * renderSettings.WavinessFrequency,
         };

         renderParticles = data;
         }
}
曲面细分的部分需要利用到computershader先插值出顶点。

half3 CurveDirrection(half3 axis, half2 uv, half amplitude, half frequency)
{
//波动，uv.x就是当前所占百分比，uv.y是第几根，这个算法看上去比较难理解，其实就是每个轴，以圆形进行运动。大波浪的感觉。
        half angle = uv.x*frequency + uv.y;

        half c = cos(angle);
        half s = sin(angle);

        half3 vecX = half3(0, c, s);
        half3 vecY = half3(c, 0, s);
        half3 vecZ = half3(c, s, 0);

        half3 vec = normalize(vecX*axis.x + vecY*axis.y + vecZ*axis.z);

        return vec*amplitude;  
}


RenderParticle GetSplineBodyData(int x, half t, uint sizeY)
{
//渲染部分因为不需要贝塞尔插值，只需要当前点和下一个点lerp一下就行了
        int sizeYm1 = sizeY - 1;
        int y = (int)(t*sizeY);
        half tStep = 1.0f / sizeY;
        half localT = (t % tStep) * sizeY;

        int startI = x*sizeY;

        int y1 = min(y, sizeYm1);
        int y2 = min(y + 1, sizeYm1);

        RenderParticle b1 = renderParticles[startI + y1];
        RenderParticle b2 = renderParticles[startI + y2];

        RenderParticle b;
        b.color = lerp(b1.color, b2.color, localT);
        b.interpolation = lerp(b1.interpolation, b2.interpolation, localT);
        b.wavinessScale = lerp(b1.wavinessScale, b2.wavinessScale, localT);
        b.wavinessFrequency = lerp(b1.wavinessFrequency, b2.wavinessFrequency, localT);

        return b;
}


float3 GetSplinePoint(int x, float t, uint sizeY)
{
//sizeYm1就是八个间隔，y就是整根头发对应下的这个部分的值。这里要找到原始的珍珠对应新的曲面细分下的百分比。localT就是把0-1这个区间划分成九个0-1，startI就是找到新的曲面细分对应原始的顶点的位置。例如x=1，那么起始点就是9，y0就是y-1,y1就是y，y2就是y+1，这里还处理的范围，防止越界。找到对应的三个点之后，就可以做贝塞尔曲线的插值。插值的部分理论上就是这个点位于这三个点中的具体位置。所以t取余1/sizeY之后，还要乘以sizeY，就是正确的插值位置。
        int sizeYm1 = sizeY - 1;
        int y = (uint)(t*sizeY);
        half tStep = 1.0f / sizeY;
        half localT = (t % tStep) * sizeY;

        int startI = x*sizeY;

        int y0 = max(0, y - 1);
        int y1 = min(y, sizeYm1);
        int y2 = min(y + 1, sizeYm1);

        float3 p0 = particles[startI + y0].position;
        float3 p1 = particles[startI + y1].position;
        float3 p2 = particles[startI + y2].position;

        float3 cPoint1 = (p0 + p1)*0.5f;
        float3 cPoint2 = (p1 + p2)*0.5f;

        return GetBezierPoint(cPoint1, p1, cPoint2, localT);
}


[numthreads(THREADS,1,1)]
void CSTesselate (uint3 id : SV_DispatchThreadID)
{
//id是这个线程在所有线程组中的位置，因为定义的yz都是1，所以其实就是id.x就是第几个线程，而一个顶点对应一个线程。tessSegments是曲面细分出的顶点数，这个可以优化成离摄像机越远，顶点数就越少。
那么y在这里就是在一个曲面细分小段里的余数，x就是第几个曲面细分小段。t就是他所占曲面细分的百分比
segments是每根头发的段数。工程里是9.
        uint y = id.x % tessSegments;
        uint x = id.x / tessSegments;

        float t = y / (float)tessSegments;
        float tessStep = 1.0/tessSegments;

    float3 tessPosition =  GetSplinePoint(x, saturate(t), segments);

        RenderParticle renderParticle = GetSplineBodyData(x, saturate(t), segments);
       
        float3 curve = CurveDirrection(normalize(wavinessAxis), half2(t, x), renderParticle.wavinessScale, renderParticle.wavinessFrequency);

        TessRenderParticle tessParticle;

        tessParticle.position = tessPosition + curve;
        tessParticle.tangent = float3(0,0,0);
        tessParticle.color = renderParticle.color;
        tessParticle.interpolation = renderParticle.interpolation;
        tessRenderParticles[id.x] = tessParticle;

        AllMemoryBarrierWithGroupSync();

        int sign = y == 0 ? -1 : 1;
        tessRenderParticles[id.x].tangent = normalize(tessParticle.position - tessRenderParticles[id.x - sign].position)*sign;
}
而着色器还需要头发索引。

private static List<int> ProcessIndicesForMesh(int startIndex, List<Vector3> scalpVertices, List<int> scalpIndices, List<Vector3> hairVertices, int segments, float accuracy = Accuracy)
{
     var hairIndices = new List<int>();
//遍历头皮的顶点，每三个构成一个三角形，对应的头发也要是3的倍数才行。遍历头发的顶点，将靠近头皮的部分一一对应，最终，每一个头皮的顶点都获得对应的一个头发的顶点的索引。
     for (var i = 0; i < scalpIndices.Count; i++)
     {
         var index = scalpIndices;
         var scalpVertex = scalpVertices[index];

         if (i % 3 == 0)
              FixNotCompletedPolygon(hairIndices);

         for (var j = 0; j < hairVertices.Count; j += segments)
         {
              var hairVertex = hairVertices[j];

              if ((hairVertex - scalpVertex).sqrMagnitude < accuracy*accuracy)
              {
                   hairIndices.Add(startIndex + j);
                   break;
              }
          }
      }

      FixNotCompletedPolygon(hairIndices);
      return hairIndices;
}

public static List<int> ProcessIndices(List<int> scalpIndices, List<Vector3> scalpVertices, List<List<Vector3>> hairVerticesGroups, int segments, float accuracy = Accuracy)
{
      var hairIndices = new List<int>();

      var grouStartIndex = 0;
      foreach (var hairVertices in hairVerticesGroups)
      {
           var groupIndices = ProcessIndicesForMesh(grouStartIndex, scalpVertices, scalpIndices, hairVertices, segments, accuracy);
           hairIndices.AddRange(groupIndices);

           grouStartIndex += hairVertices.Count;
       }
//拿到所有的索引之后，再除以9，就是每个顶点对应的头发第几根
       for (var i = 0; i < hairIndices.Count; i++)
       {
            hairIndices = hairIndices / segments;
       }

       return hairIndices;
}
接下来是插值，因为我们网格中定义的主发束就那么几百根，但人的头发要远远多于这些，所以我们需要使用插值来获得浓密的头发。这在gpugemes中也有描述：





//这是使用重心坐标插值，一开始使用了0.9,0.05,0.05的坐标，保持头发的稳定性，后面的部分就随机，总共生成了64个插值发束
private void Gen(bool forceUpdate = false)
        {
            var off = 0.1f;
            var n = 2;

            //oldN = n;
            var m = 1 - off;
            var mm = (1 - m) * 0.5f;

            barycentric.Reset();
            Split(new Vector3(m, mm, mm), new Vector3(mm, m, mm), new Vector3(mm, mm, m), n);

            while (barycentric.Count < MaxCount)
            {
                var k = GetRandomK();
                if (!barycentric.Contains(k))
                    barycentric.Add(GetRandomK());
            }

            settings.RuntimeData.Barycentrics.PushData();
        }

        private void Split(Vector3 b1, Vector3 b2, Vector3 b3, int steps)
        {
            steps--;

            TryAdd(b1);
            TryAdd(b2);
            TryAdd(b3);

            var n1 = (b1 + b2) * 0.5f;
            var n2 = (b2 + b3) * 0.5f;
            var n3 = (b3 + b1) * 0.5f;

            if (steps < 0)
                return;

            Split(b1, n1, n3, steps);
            Split(b2, n1, n2, steps);
            Split(b3, n2, n3, steps);
            Split(n1, n2, n3, steps);
        }

        private void TryAdd(Vector3 v)
        {
            if (!barycentric.Contains(v))
            {            
                barycentric.Add(v);
            }
        }


        private Vector3 GetRandomK()
        {
            var ka = Random.Range(0f, 1f);
            var kb = Random.Range(0f, 1f);

            if (ka + kb > 1)
            {
                ka = 1 - ka;
                kb = 1 - kb;
            }

            var kc = 1 - (ka + kb);
            return new Vector3(ka, kb, kc);
        }



---------------------
作者：yxriyin
来源：CSDN
原文：https://blog.csdn.net/yxriyin/article/details/88727377
版权声明：本文为博主原创文章，转载请附上博文链接！


准备工作基本完成，然后就可以开始搞着色了。先从曲面细分开始，曲面细分的基础知识可以看凯奥斯大佬的文章：

https://zhuanlan.zhihu.com/p/42550699

domain是isoline，分割模式是integer，输出是线段，输出的控制点是3个，对于每3个一组的头发根节点，
输出沿着y轴生成等y值的水平线图元。具体可以参考https://www.jianshu.com/p/3d974e69f842
[domain("isoline")]
[partitioning("integer")]
[outputtopology("line")]
[outputcontrolpoints(3)]
[patchconstantfunc("HSConst")]
HS_OUTPUT HS(InputPatch<VS_OUTPUT, 3> ip, uint id : SV_OutputControlPointID)
{
        HS_OUTPUT output;
        output.id = ip[id].id;
        return output;
}

float3 GetBarycentric(float3 a, float3 b, float3 c, fixed3 k)
{
        return a*k.x + b*k.y + c*k.z;
}

fixed GetBarycentricFixed(fixed a, fixed b, fixed c, fixed3 k)
{
        return a*k.x + b*k.y + c*k.z;
}
//op是细分完毕之后的等y值线段
StepData GetPosition(OutputPatch<HS_OUTPUT, 3> op, fixed2 uv)
{
//根据uv.y*64，就可以知道插值在第几个重心坐标，而index就是uv.x乘以第几段。
        fixed3 barycentric = _Barycentrics[uv.y*64];
            
        half index = uv.x*_TessFactor.y;
        根据重心坐标，可以得到插值后的长度，再乘以Index，就是这个顶点离根节点的长度。
        half length = GetBarycentricFixed(_Length.x, _Length.y, _Length.z, barycentric);
        half length1 = length*index;
        //op的id就是对应的根节点的id，id乘以段数，加上length1，就可以得到珍珠的下标。
        ParticleData p1 = _Particles[op[0].id*_TessFactor.y + length1];
        ParticleData p2 = _Particles[op[1].id*_TessFactor.y + length1];
        ParticleData p3 = _Particles[op[2].id*_TessFactor.y + length1];
        //根据这三个下标，再和重心坐标进行插值，得到位置后，在根据参数进行最终插值。
        float3 position = GetBarycentric(p1.position, p2.position, p3.position, barycentric);
        position = lerp(position, p1.position, p1.interpolation);

        //切线也是按照这种方式处理
        float3 tangent = GetBarycentric(p1.tangent, p2.tangent, p3.tangent, barycentric);
        tangent = lerp(tangent, p1.tangent, p1.interpolation);        
        
        float3 color = GetBarycentric(p1.color, p2.color, p3.color, barycentric);
        color = lerp(color, p1.color, p1.interpolation);

        StepData data;
        data.position = position;
        data.tangent = tangent;
        data.color = color;
        return data;
}




[domain("isoline")]
DS_OUTPUT DS(HS_CONSTANT_OUTPUT input, OutputPatch<HS_OUTPUT, 3> op, float2 uv : SV_DomainLocation)
{

        DS_OUTPUT output;

        StepData step = GetPosition(op, uv);

        float4 lightData = LightData(step.position);
        float3 lightDir = lightData.xyz;
        half attenuation = lightData.w;

        float3 viewDir = ViewDir(step.position);

        //光照模型并没有使用gpugems2中的Marschner，而是使用了更简单的Kajiya-Kay Model
        //主要参考了https://blog.csdn.net/noahzuo/article/details/51162472
        //法线用了近似计算，然后与光线方向点积，菲尼尔根据法线和视线的点积，夹角越大就越亮，
        half3 psevdoNormal = normalize(step.position - _LightCenter);
        attenuation *= Diffuse(psevdoNormal, lightDir, _Diffuse) + Fresnel(psevdoNormal, viewDir, _FresnelPower)*_FresnelAtten;


        half shift = saturate(tex2Dlod(_ColorTex, half4(uv.yx, 0, 0)).r - 0.5);
        fixed thickness = 1 - pow(2, -10 * (1 - uv.x));//curve

        output.vertex = float4(step.position, 1);
        output.tangent = step.tangent;
        output.normal = cross(step.tangent, cross(lightDir, step.tangent));
        output.viewDir = viewDir;
        output.lightDir = lightDir;
        output.factor = half4(saturate(attenuation), shift, 0, 0);
        output.right = normalize(cross(step.tangent, output.viewDir))*thickness*_StandWidth;
        output.color = step.color;

        return output;                        
}
曲面细分做完了，接下来就是几何着色器。


GS_OUTPUT CopyToFragment(DS_OUTPUT v, float4 position)
{
//算出世界坐标和裁剪空间坐标
        float4 objectPosition = mul(unity_WorldToObject, position);
        float4 clipPosition = UnityObjectToClipPos(objectPosition);

        GS_OUTPUT output;

        output.pos = clipPosition;
        output.tangent = v.tangent;
        output.normal = v.normal;
        output.viewDir = v.viewDir;
        output.lightDir = v.lightDir;
        output.factor = v.factor;
        output.color = v.color;

        TRANSFER_VERTEX_TO_FRAGMENT(output);
        UNITY_TRANSFER_FOG(output, output.pos);
        
        return output;
}

void GS(line DS_OUTPUT p[2], inout TriangleStream<GS_OUTPUT> triStream)
{
        float4 v[4];
        v[0] = float4(p[0].vertex + p[0].right, 1);
        v[1] = float4(p[1].vertex + p[1].right, 1);
        v[2] = float4(p[0].vertex - p[0].right, 1);
        v[3] = float4(p[1].vertex - p[1].right, 1);

        triStream.Append(CopyToFragment(p[0], v[0]));
        triStream.Append(CopyToFragment(p[1], v[1]));
        triStream.Append(CopyToFragment(p[0], v[2]));
        triStream.Append(CopyToFragment(p[1], v[3]));
}
最后是片段着色器

//这事各项异性高光，将切线沿着法线进行一定程度的移动，可以展现头发的高光
half3 ShiftTangent(half3 tangent, half3 normal, half shift)
{
        return normalize(tangent + shift * normal);
}

这就是标准的Kajiya-Kay Model公式
half Specular(half3 tangent, half3 viewDir, half3 lightDir, half exponent)
{
        half3 h = normalize(viewDir + lightDir);
        half dotTH = dot(tangent, h);
        half sinTH = sqrt(1.0 - dotTH * dotTH);
        half dirAtten = smoothstep(-1.0, 0.0, dotTH);
        return dirAtten * pow(sinTH, exponent);
}
fixed3 SpecularColor(GS_OUTPUT i, fixed shift, half width1, half width2, fixed3 color)
{
        half3 tangent1 = ShiftTangent(i.tangent, i.normal, i.factor.y - shift);
        half3 tangent2 = ShiftTangent(i.tangent, i.normal, i.factor.y + shift);

        half3 specular1 = Specular(tangent1, i.viewDir, i.lightDir, width1);
        half3 specular2 = Specular(tangent2, i.viewDir, i.lightDir, width2);

        return color * specular1*specular2;
}

float4 FS(GS_OUTPUT i) :SV_Target
{
        fixed3 lightColor = _LightColor0 * LIGHT_ATTENUATION(i)*i.factor.x;

        fixed3 ambient = UNITY_LIGHTMODEL_AMBIENT.rgb*i.color;
        fixed3 diffuse = Diffuse(i.normal, i.lightDir, _Diffuse)*i.color*lightColor;
        fixed3 specular = SpecularColor(i, _SpecularShift, _PrimarySpecular, _SecondarySpecular, _SpecularColor)*(max(lightColor, 0.35));

        fixed4 final = fixed4(diffuse + specular + ambient, 1);
        UNITY_APPLY_FOG(i.fogCoord, final);
        return final;
}
终于，头发有了基础的样子。要想人生过得去，带点绿。





接下来就要开始处理动力学的部分。

[numthreads(THREADS,1,1)]
void CSIntegrate (uint3 id : SV_DispatchThreadID)
{
    if(id.x >= particlesLength)
        return;

//重力加上风，移动过程中和以前的位置差，那么速度就是位置差加上加速度
        Particle particle = particles[id.x];

        float3 acceleration = (gravity + wind)*step;

        float3 difference = particle.position - particle.lastPosition;
        float3 velocity = difference*invDrag + acceleration;
        float3 nextPosition = particle.position + velocity;

        particle.lastPosition = particle.position;
        particle.position = nextPosition;
        
        particles[id.x] = particle;
}


float3 DistanceJointSolveImpl(float3 position1, float3 position2, float distance)
{
//将两个位置差和初始位置差进行比较，就可以得到缩放的百分比。
        float3 relPosition = position1 - position2;
        float actualDistance = length(relPosition);

        float penetration = (distance - actualDistance) / actualDistance;
        return relPosition*penetration;
}

void DistanceJointsSolve(DistanceJoint joint)
{
//取到这两个珍珠点，得到矫正数，两个都进行矫正。这里因为是并行的，所以不能连续矫正，只能断开成两组group分别矫正。
        Particle particle1 = particles[joint.body1Id];
        Particle particle2 = particles[joint.body2Id];

        float3 correction = DistanceJointSolveImpl(particle1.position, particle2.position, joint.distance)*joint.elasticity*step*0.5f;

        particle1.position += correction;
        particle2.position -= correction;

        particles[joint.body1Id] = particle1;
        particles[joint.body2Id] = particle2;
}

[numthreads(THREADS,1,1)]
void CSDistanceJoints (uint3 id : SV_DispatchThreadID)
{
        int i = startGroup + id.x;
        if(i < startGroup + sizeGroup)
        {
                DistanceJointsSolve(distanceJoints);
        }
}

很早就看到文章说只是这样的约束无法处理迅速移动的情况，会导致头发被拉长，于是需要一个专门处理迅速移动的判断

float3 DistanceJointSolveImpl(float3 position1, float3 position2, float distance)
{
        float3 relPosition = position1 - position2;
        float actualDistance = length(relPosition);

        float penetration = (distance - actualDistance) / actualDistance;
        return relPosition*penetration;
}

//按照现在的距离和上次的距离的比较，得到的矫正全部应用到第二个点的位置。
void DistanceJointsSolve(uint i1, uint i2, float distance)
{
        Particle particle1 = particles[i1];
        Particle particle2 = particles[i2];

        float3 correction = DistanceJointSolveImpl(particle1.position, particle2.position, distance)*step;

        particle2.position -= correction;
        particle2.lastPosition -= correction*0.9;

        particles[i2] = particle2;
}

[numthreads(THREADS,1,1)]
void CSSplineJoints (uint3 id : SV_DispatchThreadID)
{
    if(id.x*segments >= pointJointsLength)
        return;
        
//从本地坐标转到世界坐标之后的两个点
        for(uint i = 1; i < segments; i++)
        {
                uint index = id.x*segments + i;

                PointJoint joint1 = pointJoints[index - 1];
                float4x4 m1 = transforms[joint1.matrixId];
                float3 guidePosition1 = mul(m1, float4(joint1.position, 1.0)).xyz;
        
                PointJoint joint2 = pointJoints[index];
                float4x4 m2 = transforms[joint2.matrixId];
                float3 guidePosition2 = mul(m2, float4(joint2.position, 1.0)).xyz;

                float distance = length(guidePosition2 - guidePosition1);
                DistanceJointsSolve(joint1.bodyId, joint2.bodyId, distance);
        }
}
这样子，基本的动力学算是完成了。还剩下最后的影子问题了。

本来直接用自带的影子，发现精度太低，那么只能用老办法单独绘制角色阴影了，因为这个以前做过，就暂时先不做了。

最后我想到了大佬Kvant Wig的那个骚气的男子，就把他开源工程里的那个模型下载下来，配置了一下，然后跑起来看看。

感觉还行，但由于头发中间有点镂空，感觉掉发一样，于是干脆把头皮搞成和头发一样的颜色，这样就好多了。

最终演示发现头发和碰撞体依然有穿插，再精确判断电脑就有点卡了，就适当增加了以下头发的半径。暂时先这样。

最后展示结果：











最后是插件地址：

https://assetstore.unity.com/packages/vfx/shaders/best-hair-142637
---------------------

有点骚 大佬

niu比