I have been trying to write the lucas-kanade algorithm to use inside unity. I have written the algorithm in .compute file(had to take help of AI, as I am not at all familiar with shader languages). The code seems to be working in real time, the only issue I'm facing is that it doesn't work when the camera moves a little fast. Like if I move the camera slowly, it works almost perfectly but as soon as I increase the speed of the movement, it gets lost.
I have tried looking into getting the exact implementation from scratch, without opencv, but I couldn't find any. Here's my compute shader code:

// Sparse Pyramidal Lucas-Kanade Optical Flow
// Production-quality implementation with Gaussian pyramid and weighted windows

#pragma kernel KToLuma
#pragma kernel KGaussianH
#pragma kernel KGaussianV
#pragma kernel KDownsample2
#pragma kernel KComputeGradients
#pragma kernel KInitPoints
#pragma kernel KScalePoints
#pragma kernel KSparseLK


int _Width;
int _Height;


int _Radius;        // window radius (e.g., 7 => 15x15)
int _Iters;         // iterations per level (e.g., 10)
float _Eps;         // convergence eps (e.g., 1e-3)
float _MaxStep;     // clamp per-iter update (e.g., 1.0)

// Pyramid scale for p0 (base points): p0_level = p0_full * _P0Scale
float _P0Scale;

// Gaussian window sigma (typically radius / 2)
float _WinSigma;

// --------------------
// Keypoints
// --------------------
int _Count;
StructuredBuffer _Pts0Full;     // original points in FULL-res pixel coords (frame1)
StructuredBuffer _Velocity;     // predicted motion (full-res pixels) for initial guess
RWStructuredBuffer _Pts1Level;  // current estimate in CURRENT level coords (updated in place)
RWStructuredBuffer _Status;       // 1 = ok (never set to 0)

Texture2D _BlurSrc;
RWTexture2D _BlurDst;
Texture2D _SrcGray;
RWTexture2D _DstGraySmall;
Texture2D _I1;
Texture2D _I2;
RWTexture2D _GradI1;        
Texture2D _GradI1Tex;       
SamplerState sampler_linear_clamp;
static const float GAUSS5[5] = { 0.0625, 0.25, 0.375, 0.25, 0.0625 };
float SampleGray(Texture2D tex, float2 p)
{
    float2 uv = (p + 0.5) / float2(_Width, _Height);
    return tex.SampleLevel(sampler_linear_clamp, uv, 0);
}

float2 SampleGrad(float2 p)
{
    float2 uv = (p + 0.5) / float2(_Width, _Height);
    return _GradI1Tex.SampleLevel(sampler_linear_clamp, uv, 0).xy;
}


float GaussianWeight(float x, float y, float sigma)
{
    float sigma2 = sigma * sigma;
    return exp(-0.5 * (x * x + y * y) / sigma2);
}

[numthreads(8,8,1)]
void KToLuma(uint3 id : SV_DispatchThreadID)
{
    if (id.x >= (uint)_Width || id.y >= (uint)_Height) return;
    float3 rgb = _SrcColor.Load(int3(id.xy, 0)).rgb;
    _DstGray[id.xy] = dot(rgb, float3(0.2126, 0.7152, 0.0722));
}


[numthreads(8,8,1)]
void KGaussianH(uint3 id : SV_DispatchThreadID)
{
    if (id.x >= (uint)_Width || id.y >= (uint)_Height) return;
    
    int2 p = int2(id.xy);
    float sum = 0;
    
    [unroll]
    for (int i = -2; i <= 2; i++)
    {
        int2 sp = p + int2(i, 0);
        // Clamp to bounds
        sp.x = clamp(sp.x, 0, _Width - 1);
        sum += _BlurSrc.Load(int3(sp, 0)) * GAUSS5[i + 2];
    }
    
    _BlurDst[p] = sum;
}

// --------------------
// Kernel: Gaussian blur - Vertical pass (5-tap)
// --------------------
[numthreads(8,8,1)]
void KGaussianV(uint3 id : SV_DispatchThreadID)
{
    if (id.x >= (uint)_Width || id.y >= (uint)_Height) return;
    
    int2 p = int2(id.xy);
    float sum = 0;
    
    [unroll]
    for (int i = -2; i <= 2; i++)
    {
        int2 sp = p + int2(0, i);
        // Clamp to bounds
        sp.y = clamp(sp.y, 0, _Height - 1);
        sum += _BlurSrc.Load(int3(sp, 0)) * GAUSS5[i + 2];
    }
    
    _BlurDst[p] = sum;
}

// --------------------
// Kernel: Downsample by 2x (after Gaussian blur)
// Simple 2x2 average since blur already applied
// --------------------
[numthreads(8,8,1)]
void KDownsample2(uint3 id : SV_DispatchThreadID)
{
    if (id.x >= (uint)_Width || id.y >= (uint)_Height) return;

    int2 p2 = int2(id.xy);
    int2 p = p2 * 2;

    float a = _SrcGray.Load(int3(p + int2(0,0), 0));
    float b = _SrcGray.Load(int3(p + int2(1,0), 0));
    float c = _SrcGray.Load(int3(p + int2(0,1), 0));
    float d = _SrcGray.Load(int3(p + int2(1,1), 0));
    _DstGraySmall[p2] = 0.25 * (a + b + c + d);
}

// --------------------
// Kernel: Precompute gradients for I1 (Scharr-like for better accuracy)
// --------------------
[numthreads(8,8,1)]
void KComputeGradients(uint3 id : SV_DispatchThreadID)
{
    if (id.x >= (uint)_Width || id.y >= (uint)_Height) return;

    int2 p = int2(id.xy);
    
    // Scharr gradient (more accurate than simple central difference)
    // Horizontal: [-3 0 3; -10 0 10; -3 0 3] / 32
    // Vertical:   [-3 -10 -3; 0 0 0; 3 10 3] / 32
    
    float tl = _I1.Load(int3(p + int2(-1,-1), 0));
    float tc = _I1.Load(int3(p + int2( 0,-1), 0));
    float tr = _I1.Load(int3(p + int2( 1,-1), 0));
    float ml = _I1.Load(int3(p + int2(-1, 0), 0));
    float mr = _I1.Load(int3(p + int2( 1, 0), 0));
    float bl = _I1.Load(int3(p + int2(-1, 1), 0));
    float bc = _I1.Load(int3(p + int2( 0, 1), 0));
    float br = _I1.Load(int3(p + int2( 1, 1), 0));
    
    // Scharr horizontal gradient (dI/dx)
    float Ix = (3.0 * (tr - tl) + 10.0 * (mr - ml) + 3.0 * (br - bl)) / 32.0;
    
    // Scharr vertical gradient (dI/dy)
    float Iy = (3.0 * (bl - tl) + 10.0 * (bc - tc) + 3.0 * (br - tr)) / 32.0;
    
    _GradI1[p] = float4(Ix, Iy, 0, 0);
}

// --------------------
// Kernel: Initialize points at coarsest pyramid level with motion prediction
// --------------------
[numthreads(256,1,1)]
void KInitPoints(uint3 id : SV_DispatchThreadID)
{
    uint i = id.x;
    if (i >= (uint)_Count) return;

    _Status[i] = 1;
    // Use predicted position: (p0 + velocity) scaled to this level
    _Pts1Level[i] = (_Pts0Full[i] + _Velocity[i]) * _P0Scale;
}

// --------------------
// Kernel: Scale points by 2x when moving to finer level
// --------------------
[numthreads(256,1,1)]
void KScalePoints(uint3 id : SV_DispatchThreadID)
{
    uint i = id.x;
    if (i >= (uint)_Count) return;
    _Pts1Level[i] *= 2.0;
}

// --------------------
// Kernel: Sparse Lucas-Kanade iteration at current pyramid level
// Optimized: Structure tensor computed once, Gaussian window weighting
// --------------------
[numthreads(256,1,1)]
void KSparseLK(uint3 id : SV_DispatchThreadID)
{
    uint i = id.x;
    if (i >= (uint)_Count) return;

    float2 p0 = _Pts0Full[i] * _P0Scale;
    float2 p1 = _Pts1Level[i];
    float2 d = p1 - p0;

    // Gaussian sigma for window weighting
    float sigma = _WinSigma;
    
    // ---------- Compute structure tensor ONCE (doesn't depend on d) ----------
    float A = 0, B = 0, C = 0;
    
    [loop]
    for (int yy = -_Radius; yy <= _Radius; yy++)
    {
        [loop]
        for (int xx = -_Radius; xx <= _Radius; xx++)
        {
            float2 q = p0 + float2(xx, yy);
            float2 g = SampleGrad(q);
            float Ix = g.x;
            float Iy = g.y;
            
            // Gaussian weight
            float w = GaussianWeight((float)xx, (float)yy, sigma);
            
            A += w * Ix * Ix;
            B += w * Ix * Iy;
            C += w * Iy * Iy;
        }
    }
    
    float det = A * C - B * B;
    
    // Check if structure tensor is invertible
    float invDet = (abs(det) > 1e-12) ? (1.0 / det) : 0.0;
    
    // ---------- Iterative refinement (only mismatch vector changes) ----------
    [loop]
    for (int it = 0; it < _Iters; it++)
    {
        float bx = 0, by = 0;

        // Compute mismatch vector with Gaussian weighting
        [loop]
        for (int yy = -_Radius; yy <= _Radius; yy++)
        {
            [loop]
            for (int xx = -_Radius; xx <= _Radius; xx++)
            {
                float2 q  = p0 + float2(xx, yy);
                float2 q2 = q + d;

                float I1v = SampleGray(_I1, q);
                float I2v = SampleGray(_I2, q2);
                float It  = I2v - I1v;

                float2 g = SampleGrad(q);
                float Ix = g.x;
                float Iy = g.y;
                
                // Gaussian weight
                float w = GaussianWeight((float)xx, (float)yy, sigma);

                bx += w * Ix * It;
                by += w * Iy * It;
            }
        }

        // Solve: [A B; B C] * dd = -[bx by]
        float2 dd;
        dd.x = (-C * bx + B * by) * invDet;
        dd.y = ( B * bx - A * by) * invDet;

        // Clamp per-iter step
        float m = length(dd);
        if (m > _MaxStep) dd *= (_MaxStep / m);

        d += dd;

        // Convergence
        if (abs(dd.x) + abs(dd.y) < _Eps)
            break;
    }

    _Pts1Level[i] = p0 + d;
}

Can someone help me please!