零度冰山

自建一个类似Houdini的intersect()函数。

定义核心函数

判断某个位置朝特定方向是否投射到了指定三角面上

bool IntersectRayTriangle(Vector3 originPosition, Vector3 direction, Point A, Point B, Point C, out float distance)
{
    distance = 0f;

    Vector3 e1 = B.position - A.position;
    Vector3 e2 = C.position - A.position;
    Vector3 h = Vector3.Cross(direction, e2);
    float a = Vector3.Dot(e1, h);

    if (a > -Mathf.Epsilon && a < Mathf.Epsilon)
        return false;

    float f = 1.0f / a;
    Vector3 s = originPosition - A.position;
    float u = f * Vector3.Dot(s, h);

    if (u < 0.0f || u > 1.0f)
        return false;

    Vector3 q = Vector3.Cross(s, e1);
    float v = f * Vector3.Dot(direction, q);

    if (v < 0.0f || u + v > 1.0f)
        return false;

    distance = f * Vector3.Dot(e2, q);

    if (distance > Mathf.Epsilon)
        return true;

    return false;
}

主函数

代入List<Point> meshPoints，其为代表三角面mesh的点集。

distance代表最大投射距离。

输出hitCount，如果是偶数说明在某个封闭Mesh内。

加了两层多线程优化。其一是CPU多线程，其二是GPU多线程。

1. 主函数

void RayCastWithMeshPoints(Vector3 originPosition, Vector3 direction, float distance, List<Point> meshPoints, ref int hitCount)
{
    hitCount = 0;
    direction.Normalize();
    //int handle = 0;
    //while (handle <= meshPoints.Count - 3)  //剪枝，去掉射线背面的triangle，感觉没太大必要，因为时间复杂度可能还比去运行下投射函数到return的还高。
    //{
    //    Vector3 pos0 = meshPoints[handle].position, pos1 = meshPoints[handle + 1].position, pos2 = meshPoints[handle + 2].position;
    //    Vector3 dir0 = (pos0 - originPosition).normalized, dir1 = (pos1 - originPosition).normalized, dir2 = (pos2 - originPosition).normalized;

    //    if (Vector3.Dot(direction, dir0) < 0.0f && Vector3.Dot(direction, dir1) < 0.0f && Vector3.Dot(direction, dir2) < 0.0f)
    //    {
    //        meshPoints.RemoveRange(handle, 3);
    //        continue;
    //    }
    //    handle += 3;
    //}
    //if (meshPoints.Count == 0) return;


    ///<summary>
    ///单线程跑
    /// </summary>
    if (meshPoints.Count < Environment.ProcessorCount)
    {
        for (int i = 0; i < meshPoints.Count; i += 3)
        {
            if (IntersectRayTriangle(originPosition, direction, meshPoints[i], meshPoints[i + 1], meshPoints[i + 2], out float intersectionDistance))
                if (intersectionDistance <= distance)
                    hitCount++;
        }
        return;
    }

    ///<summary>
    ///CPU多线程跑
    /// </summary>
    if (meshPoints.Count < 1024)
    {
        ConcurrentBag<char> results = new ConcurrentBag<char>();
        Task[] tasks = new Task[Environment.ProcessorCount];
        for (int i = 0; i < tasks.Length; i++)
        {
            int cur_i = i;
            tasks[i] = Task.Run(() =>
            {
                int triangleCount = meshPoints.Count / 3;
                int start = cur_i * Mathf.CeilToInt(triangleCount * 1f / tasks.Length);
                int end = start + Mathf.CeilToInt(triangleCount * 1f / tasks.Length);
                for (int j = start; j < end && j < triangleCount; j++)
                {
                    int index = j * 3;
                    if (IntersectRayTriangle(originPosition, direction, meshPoints[index], meshPoints[index + 1], meshPoints[index + 2], out float intersectionDistance))
                        if (intersectionDistance <= distance)
                            results.Add('H');
                }
            });
        }
        Task.WaitAll(tasks);
        hitCount = results.Count;
        return;
    }

    ///<summary>
    ///GPU多线程跑
    /// </summary>
    {
        //AssetDatabase.FindAssets("GeneralFunctionsGPUKernels.compute")[0];
        int ThreadCount = 1024;
        ComputeShader computeShader = AssetDatabase.LoadAssetAtPath<ComputeShader>("Assets/Plugins/GlightUtilities/PCGDevKit/Utilities/GeneralFunctionsGPUKernels.compute"); ;
        ComputeBuffer originPosition_inputBuffer, direction_inputBuffer, meshPoints_inputBuffer, triangleCount_inputBuffer, hitDistance_outputBuffer;

        originPosition_inputBuffer = new ComputeBuffer(ThreadCount, sizeof(float) * 3);
        direction_inputBuffer = new ComputeBuffer(ThreadCount, sizeof(float) * 3);
        meshPoints_inputBuffer = new ComputeBuffer(meshPoints.Count, sizeof(float) * 3);
        triangleCount_inputBuffer = new ComputeBuffer(ThreadCount, sizeof(int));
        hitDistance_outputBuffer = new ComputeBuffer(meshPoints.Count / 3, sizeof(float));

        float3[] _originPosition = new float3[ThreadCount], _direction = new float3[ThreadCount], _meshPoints = new float3[meshPoints.Count];
        int[] _triangleCount = new int[ThreadCount];
        for (int i = 0; i < ThreadCount; i++)
        {
            _originPosition[i] = originPosition;
            _direction[i] = direction;
            _triangleCount[i] = meshPoints.Count / 3;
        }
        for (int i = 0; i < meshPoints.Count; i++)
            _meshPoints[i] = meshPoints[i].position;


        originPosition_inputBuffer.SetData(_originPosition);
        direction_inputBuffer.SetData(_direction);
        meshPoints_inputBuffer.SetData(_meshPoints);
        triangleCount_inputBuffer.SetData(_triangleCount);

        int kernel = computeShader.FindKernel("IntersectRayTriangle");
        computeShader.SetBuffer(kernel, "originPosition", originPosition_inputBuffer);
        computeShader.SetBuffer(kernel, "direction", direction_inputBuffer);
        computeShader.SetBuffer(kernel, "points", meshPoints_inputBuffer);
        computeShader.SetBuffer(kernel, "triangleCount", triangleCount_inputBuffer);
        computeShader.SetBuffer(kernel, "hitDistance", hitDistance_outputBuffer);

        computeShader.Dispatch(kernel, ThreadCount, 1, 1);
        float[] hitDistance = new float[meshPoints.Count / 3];

        hitDistance_outputBuffer.GetData(hitDistance);

        for (int i = 0; i < meshPoints.Count / 3; i++)
        {
            if (hitDistance[i] > Mathf.Epsilon && hitDistance[i] <= distance)
                hitCount++;
        }

        direction_inputBuffer.Release();
        meshPoints_inputBuffer.Release();
        originPosition_inputBuffer.Release();
        triangleCount_inputBuffer.Release();
        hitDistance_outputBuffer.Release();
    }
}

2. ComputeShader中

#pragma kernel IntersectRayTriangle
RWStructuredBuffer<uint> triangleCount; 
RWStructuredBuffer<float3> originPosition; 
RWStructuredBuffer<float3> direction;   
RWStructuredBuffer<float3> points;      
RWStructuredBuffer<float> hitDistance;  
[numthreads(1024, 1, 1)]
void IntersectRayTriangle(uint3 id : SV_DispatchThreadID)
{
    uint step = (uint)ceil(triangleCount[id.x] / 1024.0);
    uint start = id.x * step;
    uint end = start + step;
    for (uint i = start; i < end && i < triangleCount[id.x]; i++)
    {
        float3 p0 = points[i * 3 + 0];
        float3 p1 = points[i * 3 + 1];
        float3 p2 = points[i * 3 + 2];
        hitDistance[i] = -1.0;
        
        float3 e1 = p1 - p0;
        float3 e2 = p2 - p0;
        
        float3 h = cross(direction[id.x], e2);
        float a = dot(e1, h);

        if (a > -0.000000001 && a < 0.000000001)
            continue;

        float f = 1.0f / a;
        float3 s = originPosition[id.x] - p0;
        float u = f * dot(s, h);

        if (u < 0.0f || u > 1.0f)
            continue;

        float3 q = cross(s, e1);
        float v = f * dot(direction[id.x], q);

        if (v < 0.0f || u + v > 1.0f)
            continue;

        
        hitDistance[i] = f * dot(e2, q);
        
        if (hitDistance[i] < 0.000000001)
            hitDistance[i] = -1.0;
    }
    
}