Titano/Titano/Motion/MotionBlurRenderer.cs

using System.Runtime.CompilerServices;
using Titano.Imaging;

namespace Titano.Motion;

/// <summary>Parametri del motion blur sintetico, esposti nel pannello "Elaborazione immagini".</summary>
public sealed class MotionBlurSettings
{
    public bool Enabled { get; set; } = true;

    /// <summary>Shutter angle desiderato: 180° è la convenzione cinematografica.</summary>
    public double TargetShutterAngle { get; set; } = 180.0;

    /// <summary>Quota della sfocatura mancante effettivamente sintetizzata.</summary>
    public double Strength { get; set; } = 1.0;

    /// <summary>Limite superiore della scia, in pixel: protegge dalle stime di movimento errate.</summary>
    public double MaxBlurPixels { get; set; } = 48.0;

    /// <summary>Numero massimo di campioni per pixel lungo la scia.</summary>
    public int MaxSamples { get; set; } = 25;

    public MotionBlurSettings Clone() => (MotionBlurSettings)MemberwiseClone();
}

/// <summary>
/// Motion blur direzionale sintetico lungo il campo vettoriale.
///
/// Il fattore di sfocatura mancante deriva dalla composizione delle varianze: la scia
/// realmente incisa nel fotogramma (proporzionale allo shutter angle di scatto) e quella
/// sintetica si sommano in quadratura, quindi per raggiungere l'apertura obiettivo serve
/// una scia lunga √(target² − reale²) volte lo spostamento. Sommare linearmente
/// produrrebbe un'immagine sistematicamente troppo morbida.
///
/// Il filtro di ricostruzione è di tipo "gather": ogni campione contribuisce al pixel
/// centrale solo se la propria scia lo raggiunge davvero, così lo sfondo fermo non viene
/// trascinato dentro i soggetti in movimento.
/// </summary>
public static class MotionBlurRenderer
{
    /// <summary>
    /// Frazione di spostamento da sintetizzare per portare <paramref name="actualAngle"/>
    /// al valore <paramref name="targetAngle"/>. Zero se il fotogramma è già abbastanza mosso.
    /// </summary>
    public static double MissingBlurFactor(double actualAngle, double targetAngle, double strength)
    {
        double target = Math.Clamp(targetAngle, 0, 360) / 360.0;
        double actual = Math.Clamp(actualAngle, 0, 360) / 360.0;
        if (target <= actual) return 0;
        return Math.Sqrt(target * target - actual * actual) * Math.Clamp(strength, 0, 1);
    }

    /// <summary>
    /// Applica la sfocatura da <paramref name="source"/> a <paramref name="destination"/>.
    /// Restituisce la lunghezza media della scia effettivamente resa, in pixel.
    /// </summary>
    public static double Render(ImageBuffer source, ImageBuffer destination, MotionField field,
                                double missingFactor, MotionBlurSettings settings)
    {
        if (missingFactor <= 1e-4)
        {
            destination.CopyFrom(source);
            return 0;
        }

        int width = source.Width;
        int height = source.Height;
        float factor = (float)missingFactor;
        float maxBlur = (float)Math.Max(1.0, settings.MaxBlurPixels);
        int maxSamples = Math.Clamp(settings.MaxSamples, 3, 129);

        var src = source.Data;
        var dst = destination.Data;
        double lengthSum = 0;
        object sumLock = new();

        Parallel.For(0, height, () => 0.0, (y, _, localSum) =>
        {
            int rowBase = y * width * ImageBuffer.Channels;

            for (int x = 0; x < width; x++)
            {
                int index = rowBase + x * ImageBuffer.Channels;
                field.Sample(x, y, out float vx, out float vy);

                float magnitude = MathF.Sqrt(vx * vx + vy * vy);
                float length = Math.Min(magnitude * factor, maxBlur);
                localSum += length;

                if (length < 0.5f || magnitude < 1e-4f)
                {
                    dst[index] = src[index];
                    dst[index + 1] = src[index + 1];
                    dst[index + 2] = src[index + 2];
                    continue;
                }

                float dirX = vx / magnitude;
                float dirY = vy / magnitude;
                float half = length * 0.5f;

                int taps = Math.Min(maxSamples, Math.Max(3, (int)MathF.Ceiling(length) | 1));
                int side = taps / 2;
                float step = half / side;

                // Il pixel centrale è sempre presente: garantisce continuità con le zone ferme.
                float accR = src[index], accG = src[index + 1], accB = src[index + 2];
                float weightSum = 1f;

                for (int k = 1; k <= side; k++)
                {
                    float offset = k * step;

                    weightSum += Accumulate(src, field, width, height, x, y,
                                            dirX, dirY, offset, factor, maxBlur, ref accR, ref accG, ref accB);
                    weightSum += Accumulate(src, field, width, height, x, y,
                                            dirX, dirY, -offset, factor, maxBlur, ref accR, ref accG, ref accB);
                }

                float inv = 1f / weightSum;
                dst[index] = accR * inv;
                dst[index + 1] = accG * inv;
                dst[index + 2] = accB * inv;
            }

            return localSum;
        },
        localSum =>
        {
            lock (sumLock) lengthSum += localSum;
        });

        return lengthSum / Math.Max(1, (long)width * height);
    }

    /// <summary>
    /// Aggiunge il contributo del campione spostato di <paramref name="offset"/> lungo la scia.
    /// Il peso vale 1 solo se il movimento proprio del campione arriva a coprire il pixel centrale.
    /// </summary>
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static float Accumulate(float[] src, MotionField field, int width, int height,
                                    int x, int y, float dirX, float dirY, float offset,
                                    float factor, float maxBlur,
                                    ref float accR, ref float accG, ref float accB)
    {
        float sx = x + dirX * offset;
        float sy = y + dirY * offset;
        if (sx < 0 || sy < 0 || sx > width - 1 || sy > height - 1) return 0f;

        field.Sample(sx, sy, out float tvx, out float tvy);
        float projection = MathF.Abs(tvx * dirX + tvy * dirY) * factor * 0.5f;
        projection = MathF.Min(projection, maxBlur * 0.5f);

        // Transizione morbida su un pixel: evita i gradini sui bordi della scia.
        float weight = projection - MathF.Abs(offset) + 1f;
        if (weight <= 0f) return 0f;
        if (weight > 1f) weight = 1f;

        SampleBilinear(src, width, height, sx, sy, out float r, out float g, out float b);
        accR += r * weight;
        accG += g * weight;
        accB += b * weight;
        return weight;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static void SampleBilinear(float[] data, int width, int height, float x, float y,
                                        out float r, out float g, out float b)
    {
        if (x < 0) x = 0; else if (x > width - 1) x = width - 1;
        if (y < 0) y = 0; else if (y > height - 1) y = height - 1;

        int x0 = (int)x;
        int y0 = (int)y;
        int x1 = x0 + 1 < width ? x0 + 1 : x0;
        int y1 = y0 + 1 < height ? y0 + 1 : y0;
        float fx = x - x0;
        float fy = y - y0;

        int i00 = (y0 * width + x0) * ImageBuffer.Channels;
        int i10 = (y0 * width + x1) * ImageBuffer.Channels;
        int i01 = (y1 * width + x0) * ImageBuffer.Channels;
        int i11 = (y1 * width + x1) * ImageBuffer.Channels;

        float w00 = (1 - fx) * (1 - fy);
        float w10 = fx * (1 - fy);
        float w01 = (1 - fx) * fy;
        float w11 = fx * fy;

        r = data[i00] * w00 + data[i10] * w10 + data[i01] * w01 + data[i11] * w11;
        g = data[i00 + 1] * w00 + data[i10 + 1] * w10 + data[i01 + 1] * w01 + data[i11 + 1] * w11;
        b = data[i00 + 2] * w00 + data[i10 + 2] * w10 + data[i01 + 2] * w01 + data[i11 + 2] * w11;
    }
}