validateTransientShape.m

function [ok, info] = validateTransientShape(s, pk, fps, sigma, P, prominence, sd)

%% VALIDATE THAT A DETECTED PEAK HAS CALCIUM-TRANSIENT SHAPE
% Checks a candidate peak for the characteristic FAST RISE / SLOW DECAY of a
% real calcium transient, rejecting slow drift, slow humps, step artifacts, and
% noise. This is the biophysical prior: an action-potential-coupled calcium
% transient rises fast (Ca influx + indicator on-kinetics) and decays slowly
% (indicator off-kinetics + Ca extrusion).
%
% INPUTS:
%   s          - one cell's smoothed (detection) signal, a row vector
%   pk         - frame index of the candidate peak in s
%   fps        - frames per second
%   sigma      - robust noise sigma (accepted for backward compatibility; unused)
%   P          - struct of thresholds (the two tunable knobs):
%       .maxRiseTimeSeconds  the 20%->80% rise time must be FASTER than this.
%                            Lower = stricter. Raise it if real transients are
%                            being rejected for "rise too slow".
%       .minDecayToRiseRatio the decay must be at least this many times slower
%                            than the rise (transient asymmetry). Higher =
%                            stricter. Lower if fast/compound events are wrongly
%                            rejected as "not asymmetric".
%   prominence - the peak's MinPeakProminence (from findpeaks). REQUIRED: it is
%                the amplitude reference. Using the prominence (peak height above
%                its bounding valleys) instead of the nearest local minimum is
%                what stops a broad slow hump from being scored as a 1-frame
%                micro-transient sitting on a crest wiggle.
%   sd         - (optional) precomputed Savitzky-Golay-smoothed derivative of s
%                (per second). Recomputed internally if omitted.
%
% OUTPUTS:
%   ok   - true if the event looks like a real transient
%   info - struct with measured features (riseTime, decayTime, amplitude)
%          and a 'reason' string (useful for the preview/QC tool)

    ok = false;
    info = struct('riseTime', NaN, 'decayTime', NaN, 'amplitude', NaN, 'reason', '');

    T = numel(s);
    if pk < 2 || pk > T - 1 || ~isfinite(fps) || fps <= 0
        info.reason = 'invalid peak/fps';
        return;
    end
    if nargin < 6 || isempty(prominence) || ~isfinite(prominence) || prominence <= 0
        info.reason = 'no prominence';
        return;
    end
    if nargin < 7 || isempty(sd)
        sd = smoothedDerivative(s, fps);
    end

    % --- Amplitude & baseline from the PEAK PROMINENCE (robust, non-local) ---
    % The prominence is the peak's height above the higher of its two bounding
    % valleys, so baseline = peak - prominence is a noise-immune foot. It does
    % NOT collapse onto a crest noise-wiggle the way a nearest-local-min walk-back
    % does, so a broad slow hump keeps its true (large) amplitude and therefore a
    % slow, correctly-measured rise.
    amp = prominence;
    baseline = s(pk) - amp;
    info.amplitude = amp;

    level20 = baseline + 0.20 * amp;
    level80 = baseline + 0.80 * amp;

    % --- Rise time = 20%->80% fractional rise on the rising edge (nearest
    %     crossings to the peak). Ignores the slow exponential tails that heavy
    %     smoothing adds, so it reflects the true fast rise. ---
    i = pk;
    while i > 1 && s(i) > level80, i = i - 1; end
    t80 = i;
    while i > 1 && s(i) > level20, i = i - 1; end
    t20 = i;
    info.riseTime = max(1, t80 - t20) / fps;
    if info.riseTime > P.maxRiseTimeSeconds
        info.reason = 'rise too slow';
        return;
    end

    % --- Decay end = the next robust trough after the peak, found from the
    %     SG-smoothed derivative's negative->positive crossing (the same method
    %     findCalciumSpikes uses for decay windows). Bounding the decay at the
    %     valley means an interrupted decay -- a following transient rising before
    %     the signal returns to 20% -- is measured as a partial decay instead of
    %     being wrongly rejected as a step/plateau. ---
    decayEnd = followingTrough(sd, pk, fps);
    if isnan(decayEnd), decayEnd = T; end

    d80 = NaN; d20 = NaN;
    for f = pk : decayEnd
        if isnan(d80) && s(f) <= level80, d80 = f; end
        if s(f) <= level20, d20 = f; break; end
    end
    if isnan(d80), d80 = pk; end
    if isnan(d20)
        % Never returned to 20% within the decay window.
        if decayEnd <= pk + 1
            info.reason = 'no decay (step/plateau)';   % never turned over at all
            return;
        end
        info.decayTime = max(1, decayEnd - d80) / fps;  % partial decay to the valley
    else
        info.decayTime = max(1, d20 - d80) / fps;
    end

    % --- Asymmetry: the decay must be slower than the rise ---
    if info.decayTime < P.minDecayToRiseRatio * info.riseTime
        info.reason = 'not asymmetric (decay too fast)';
        return;
    end

    ok = true;
    info.reason = 'valid';

end

% ================= local helpers =================

function sd = smoothedDerivative(s, fps)
    d = diff(s) * fps;
    w = 21;                          % ~matches findCalciumSpikes
    if numel(d) < w
        w = numel(d) - (1 - mod(numel(d), 2));   % largest odd <= length
    end
    if w >= 5
        sd = sgolayfilt(d, 3, w);
    else
        sd = d;
    end
end

function e = followingTrough(sd, pk, fps)
    % First robust local minimum after the peak: three descending derivative
    % frames followed by three non-descending ones (negative->positive crossing).
    e = NaN;
    ss = max(1, round(0.5 * fps));   % skip the immediate post-peak neighbourhood
    for f = (pk + ss) : (numel(sd) - 3)
        if all(sd(f-2:f) < 0) && all(sd(f+1:f+3) >= 0)
            e = f + 1;
            return;
        end
    end
end