Titano/Titano/Analysis/ExposureProcessor.cs

using System.Numerics;
using Titano.Imaging;

namespace Titano.Analysis;

/// <summary>
/// Applicazione dei guadagni di esposizione sui canali colore.
///
/// L'operazione avviene in luce lineare (moltiplicazione pura, fisicamente corretta) e,
/// se richiesto, con una compressione dolce delle alte luci: la parte di segnale che
/// supererebbe il bianco viene ripiegata con una tangente iperbolica invece di essere
/// troncata, evitando le classiche macchie piatte e le derive di tinta sui bordi bruciati.
/// </summary>
public static class ExposureProcessor
{
    /// <summary>
    /// Moltiplica il buffer per i tre guadagni indicati, in-place.
    /// Restituisce la frazione stimata di campioni saturati dopo la correzione.
    /// </summary>
    public static double Apply(ImageBuffer frame, ReadOnlySpan<double> channelGain,
                               bool protectHighlights, double knee)
    {
        float gr = (float)channelGain[0];
        float gg = (float)channelGain[1];
        float gb = (float)channelGain[2];

        bool identity = Math.Abs(gr - 1f) < 1e-5f && Math.Abs(gg - 1f) < 1e-5f && Math.Abs(gb - 1f) < 1e-5f;
        bool rolloff = protectHighlights && Math.Max(gr, Math.Max(gg, gb)) > 1.0001f;

        var data = frame.Data;
        int count = frame.SampleCount;

        if (!identity)
        {
            if (rolloff) ApplyWithRolloff(data, count, gr, gg, gb, (float)Math.Clamp(knee, 0.05, 0.98));
            else ApplyLinear(data, count, gr, gg, gb);
        }

        return MeasureClipping(data, count);
    }

    /// <summary>
    /// Applica due terne di guadagni miscelate secondo la maschera delle regioni.
    ///
    /// La miscela avviene in scala logaritmica: fra un guadagno di 1× e uno di 4×, il valore
    /// intermedio è 2×, non 2,5×. È l'unica interpolazione che produce una transizione
    /// visivamente uniforme lungo la sfumatura fra cielo e paesaggio; in scala lineare il
    /// passaggio si addenserebbe verso il guadagno più alto e si vedrebbe come una banda.
    /// </summary>
    public static double ApplyRegional(ImageBuffer frame, RegionMask mask,
                                       ReadOnlySpan<double> gainHigh, ReadOnlySpan<double> gainLow,
                                       bool protectHighlights, double knee)
    {
        Span<float> logHigh = stackalloc float[ImageBuffer.Channels];
        Span<float> logLow = stackalloc float[ImageBuffer.Channels];
        float maxGain = 0;

        for (int c = 0; c < ImageBuffer.Channels; c++)
        {
            logHigh[c] = (float)Math.Log2(Math.Max(gainHigh[c], 1e-6));
            logLow[c] = (float)Math.Log2(Math.Max(gainLow[c], 1e-6));
            maxGain = Math.Max(maxGain, (float)Math.Max(gainHigh[c], gainLow[c]));
        }

        bool rolloff = protectHighlights && maxGain > 1.0001f;
        float kneeValue = (float)Math.Clamp(knee, 0.05, 0.98);
        float span = 1f - kneeValue;

        int width = frame.Width;
        int height = frame.Height;
        var data = frame.Data;
        float invWidth = width > 1 ? 1f / (width - 1) : 0f;
        float invHeight = height > 1 ? 1f / (height - 1) : 0f;

        // Il compilatore non lascia usare uno Span dentro una lambda: le tre coppie di
        // esponenti passano quindi per variabili locali, che è anche la forma più veloce.
        float hR = logHigh[0], hG = logHigh[1], hB = logHigh[2];
        float lR = logLow[0], lG = logLow[1], lB = logLow[2];

        Parallel.For(0, height, y =>
        {
            int rowBase = y * width * ImageBuffer.Channels;
            float v = y * invHeight;

            for (int x = 0; x < width; x++)
            {
                float weight = mask.Sample(x * invWidth, v);
                float gr = Exp2(lR + (hR - lR) * weight);
                float gg = Exp2(lG + (hG - lG) * weight);
                float gb = Exp2(lB + (hB - lB) * weight);

                int index = rowBase + x * ImageBuffer.Channels;
                data[index] = Compress(data[index] * gr, rolloff, kneeValue, span);
                data[index + 1] = Compress(data[index + 1] * gg, rolloff, kneeValue, span);
                data[index + 2] = Compress(data[index + 2] * gb, rolloff, kneeValue, span);
            }
        });

        return MeasureClipping(data, frame.SampleCount);
    }

    private static float Exp2(float value) => MathF.Pow(2f, value);

    private static float Compress(float value, bool rolloff, float knee, float span)
        => !rolloff || value <= knee ? value : knee + span * MathF.Tanh((value - knee) / span);

    private static void ApplyLinear(float[] data, int count, float gr, float gg, float gb)
    {
        int width = Vector<float>.Count;

        if (width >= 4 && count >= width * 3)
        {
            // I guadagni si ripetono ogni 3 campioni: con vettori di larghezza non multipla di 3
            // il pattern si richiude su 3 vettori consecutivi (lcm(3, width) / width == 3).
            var phases = BuildGainPhases(gr, gg, gb, width, out int phaseCount);
            int blocks = count / width;
            int i = 0;

            for (int block = 0; block < blocks; block++, i += width)
            {
                var v = new Vector<float>(data, i);
                (v * phases[block % phaseCount]).CopyTo(data, i);
            }

            for (; i < count; i++) data[i] *= GainFor(i, gr, gg, gb);
            return;
        }

        for (int i = 0; i < count; i++) data[i] *= GainFor(i, gr, gg, gb);
    }

    private static void ApplyWithRolloff(float[] data, int count, float gr, float gg, float gb, float knee)
    {
        float span = 1f - knee;
        for (int i = 0; i < count; i++)
        {
            float v = data[i] * GainFor(i, gr, gg, gb);
            data[i] = v <= knee ? v : knee + span * MathF.Tanh((v - knee) / span);
        }
    }

    private static Vector<float>[] BuildGainPhases(float gr, float gg, float gb, int width, out int phaseCount)
    {
        phaseCount = width % 3 == 0 ? 1 : 3;
        var phases = new Vector<float>[phaseCount];
        var scratch = new float[width];

        for (int phase = 0; phase < phaseCount; phase++)
        {
            for (int k = 0; k < width; k++)
            {
                scratch[k] = GainFor(phase * width + k, gr, gg, gb);
            }
            phases[phase] = new Vector<float>(scratch);
        }
        return phases;
    }

    private static float GainFor(int sampleIndex, float gr, float gg, float gb)
        => (sampleIndex % 3) switch { 0 => gr, 1 => gg, _ => gb };

    /// <summary>Stima del clipping su campionamento regolare: costo indipendente dalla risoluzione.</summary>
    private static double MeasureClipping(float[] data, int count)
    {
        if (count == 0) return 0;
        int step = Math.Max(3, (count / 3 / 200_000) * 3);
        int clipped = 0, samples = 0;

        for (int i = 0; i + 2 < count; i += step)
        {
            float luma = ColorSpace.Luminance(data[i], data[i + 1], data[i + 2]);
            if (luma >= 0.995f) clipped++;
            samples++;
        }
        return samples == 0 ? 0 : (double)clipped / samples;
    }
}