PhysarumSimulation.compute

// Each #kernel tells which function to compile; you can have many kernels
#pragma kernel InitParticles
#pragma kernel InitInput
#pragma kernel Simulate
#pragma kernel Reset
#pragma kernel DiffuseDecay

// Create a RenderTexture with enableRandomWrite flag and set it
// with cs.SetTexture
RWTexture2D<float4> _resultTexture;

// texture of input patterns; use r channel
Texture2D<float4> _inputTexture;

// use texture for input or use noise instead
bool _useInputTexture;

// number of particles
int _particleNum;

// particle to move around
struct Particle
{
    float2 position;
    float2 direction;
};

// buffer for particles
RWStructuredBuffer<Particle> _particles;

// rotation matrix for sensor right
float3x3 _sensorRRot;

// rotation matrix for sensor left
float3x3 _sensorLRot;

// rotation matrix for turn right
float3x3 _turnRRot;

// rotation matrix for turn left
float3x3 _turnLRot;

// sensor angle in degrees
float _sensorAngle0;
float _sensorAngle1;

// turn angle in degrees
float _turnAngle0;
float _turnAngle1;

// move distance of particles per frame
float _moveDistance0;
float _moveDistance1;

// sensor distance of particles
float _sensorDistance0;
float _sensorDistance1;

// deposit value
float _depositValue0;
float _depositValue1;

// decay of the input
float _decayFactor0;
float _decayFactor1;

// texture resolution
uint2 _resultTexRes;

/* noise functions */
// get float random value
float random(float3 value, float3 dotDir)
{
    float3 smallV = sin(value);
    float random = dot(smallV, dotDir);
    random = frac(sin(random) * 123574.43212);
    return random;
}

// get 3d vector random value
float3 random3d(float3 value)
{
    return float3 (random(value, float3(12.898, 68.54, 37.7298)),
        random(value, float3(39.898, 26.54, 85.7238)),
        random(value, float3(76.898, 12.54, 8.6788)));
}

// get noise
float3 noise3d(float3 value, float scale, float offset)
{
    value *= scale;
    value.x += offset * 5;
    float3 interp = frac(value);
    interp = smoothstep(0.0, 1.0, interp);

    float3 ZValues[2];
    for (int z = 0; z <= 1; z++)
    {
        float3 YValues[2];
        for (int y = 0; y <= 1; y++)
        {
            float3 XValues[2];
            for (int x = 0; x <= 1; x++)
            {
                float3 cell = floor(value) + float3(x, y, z);
                XValues[x] = random3d(cell);
            }
            YValues[y] = lerp(XValues[0], XValues[1], interp.x);
        }
        ZValues[z] = lerp(YValues[0], YValues[1], interp.y);
    }

    float3 noise = -1.0 + 2.0 * lerp(ZValues[0], ZValues[1], interp.z);
    return noise;
}

// wrap position by the texture resolution
float2 WrapPosition(float2 pos)
{
    float wrappedPos = pos;
    pos.x = pos.x > _resultTexRes.x ? pos.x - _resultTexRes.x : pos.x;
    pos.x = pos.x < 0 ? pos.x + _resultTexRes.x : pos.x;
    pos.y = pos.y > _resultTexRes.y ? pos.y - _resultTexRes.y : pos.y;
    pos.y = pos.y < 0 ? pos.y + _resultTexRes.y : pos.y;

    return pos;
}

// initialize rotation matrix
float3x3 GetRotMatrix(float degree)
{
    float3x3 newRotMatrix = 
    {
        cos(radians(degree)), -sin(radians(degree)), 0,
        sin(radians(degree)), cos(radians(degree)), 0,
        0, 0, 1
    };
    return newRotMatrix;
}

// initialize particles for simulation
[numthreads(8, 1, 1)]
void InitParticles(uint3 id : SV_DispatchThreadID)
{
    Particle newParticle;

    // get resolution of the _resultTexture
    _resultTexture.GetDimensions(_resultTexRes.x, _resultTexRes.y);
    
    float3 baseNoise = noise3d((float3)id, 1, (float)id.x) * 0.5 + 0.5;

    // set initial position
    newParticle.position = baseNoise.xy * _resultTexRes;
    newParticle.position = WrapPosition(newParticle.position);

    // set initial direction
    newParticle.direction = float2(lerp(-1, 1, baseNoise.z), lerp(-1, 1, baseNoise.y));
    newParticle.direction = normalize(newParticle.direction);

    _particles[id.x] = newParticle;

}

// initialize input grid for simulation
[numthreads(8, 8, 1)]
void InitInput(uint3 id : SV_DispatchThreadID)
{
    //float3 baseNoise = noise3d((float3)id, 0.05, 0) * 0.5 + 0.5;
    //_resultTexture[id.xy] = float4(baseNoise.xxx, 1);
    _resultTexture[id.xy] = float4(0, 0,0, 1);

}

// actual physarum simulation
[numthreads(8,1,1)]
void Simulate(uint3 id : SV_DispatchThreadID)
{
    // separate parameters by the optional input
    float sensorAngle;
    float turnAngle;
    float sensorDistance;
    float moveDistance;
    float depositValue;
    if (_useInputTexture && _inputTexture[(uint2)_particles[id.x].position].r > 0.5)
    {
        sensorAngle = _sensorAngle1;
        turnAngle = _turnAngle1;
        sensorDistance = _sensorDistance1;
        moveDistance = _moveDistance1;
        depositValue = _depositValue1;
    }
    else 
    {
        sensorAngle = _sensorAngle0;
        turnAngle = _turnAngle0;
        sensorDistance = _sensorDistance0;
        moveDistance = _moveDistance0;
        depositValue = _depositValue0;
    }
    
    // calculate sensor position
    float2 sensorCPos = _particles[id.x].position + _particles[id.x].direction * sensorDistance;
    float2 sensorLPos = _particles[id.x].position + mul(_sensorLRot, _particles[id.x].direction) * sensorDistance;
    float2 sensorRPos = _particles[id.x].position + mul(_sensorRRot, _particles[id.x].direction) * sensorDistance;

    // wrap sensor position
    sensorCPos = WrapPosition(sensorCPos);
    sensorLPos = WrapPosition(sensorLPos);
    sensorRPos = WrapPosition(sensorRPos);

    // read texture value
    float centerValue = _resultTexture[uint2(sensorCPos)].r;
    float leftValue = _resultTexture[uint2(sensorLPos)].r;
    float rightValue = _resultTexture[uint2(sensorRPos)].r;

    // initialize rotation matrix
    _sensorRRot = GetRotMatrix(sensorAngle);
    _sensorLRot = GetRotMatrix(-sensorAngle);
    _turnRRot = GetRotMatrix(turnAngle);
    _turnLRot = GetRotMatrix(-turnAngle);
    
    // calculate turn angle
    float2 newDir = _particles[id.x].direction;
    if (rightValue > centerValue && leftValue > centerValue)
    {
        newDir = random(float3(_particles[id.x].position, 0), float3(_particles[id.x].direction, 0)) < 0.5 ?
            mul(_turnRRot, _particles[id.x].direction) : mul(_turnLRot, _particles[id.x].direction);
    }
    else if (leftValue > centerValue && leftValue > rightValue)
    {
        // in case the leftValue is the largest
        newDir = mul(_turnLRot, _particles[id.x].direction);
    }
    else if (rightValue > leftValue && rightValue > centerValue)
    {
        // in case the rightValue is the largest
        newDir = mul(_turnRRot, _particles[id.x].direction);
    }
    else
    {
        newDir = mul(_turnRRot, _particles[id.x].direction);
    }
    _particles[id.x].direction = newDir;

    // set new position of particles
    float2 newPos = _particles[id.x].position + _particles[id.x].direction * moveDistance;
    newPos = WrapPosition(newPos);
    _particles[id.x].position = newPos;

    // deposit
    uint2 depositCoord = (uint2)newPos;
    _resultTexture[depositCoord] += depositValue;
    _resultTexture[depositCoord] = clamp(_resultTexture[depositCoord], 0, 1);
}

// decay texture paint result
[numthreads(8, 8, 1)]
void DiffuseDecay(uint3 id : SV_DispatchThreadID)
{
    // separate parameters by the optional input
    float decayFactor;
    if (_useInputTexture && _inputTexture[id.xy].r > 0.5)
    {
        decayFactor = _decayFactor1;
    }
    else
    {
        decayFactor = _decayFactor0;
    }
    
    // diffuse result texture
    // define sample range in pixels neareby the target pixel
    int2 sampleRage = int2(-5, 5);

    // sum of the color value
    float4 colorSum = float4(0, 0, 0, 0);
    
    // number of pixels to add
    int count = 0;

    // sample pixel color
    for (int x = sampleRage.x; x <= sampleRage.y; x++)
    {
        for (int y = sampleRage.x; y < sampleRage.y; y++)
        {
            uint2 coord = uint2((int)id.x + x, (int)id.y + y);

            // execute in the range of the resolution of _resultTexture
            //if (coord.x < 0 || coord.x >_resultTexRes.x || coord.y < 0 || coord.y >_resultTexRes.y) { continue; }

            colorSum += _resultTexture[coord];
            count++;
        }
    }
    
    colorSum /= (float)count;

    // apply decay
    colorSum /= decayFactor;
    _resultTexture[id.xy] = float4(colorSum.xyz, 1);
}

// reset simulation
[numthreads(8, 8, 1)]
void Reset(uint3 id : SV_DispatchThreadID)
{
    _resultTexture[id.xy] = float4(0,0,0,1);
    //float3 baseNoise = (noise3d((float3)id, 0.05, 0) + 1) / 2;
    //_resultTexture[id.xy] = float4(baseNoise.xxx, 1);
}