Dealing with pixels clipped at 0 for SNR calculation

[ethanbb]

Recently I've been trying to figure out why the SNR criterion isn't working well for my recordings. One thing I notice pretty often is that ROIs in dark parts of the image have unusually high SNRs, like 8 and above. At the same time, in the good parts of the image, higher SNR is better (and there may be some overlap between the "BS" and "good" SNR ranges), so this is a problem if I want to set a single SNR threshold.

I believe this is happening because some parts of the image are dim enough that the signal is clipping at 0. This makes the raw distributions of pixels in these areas very skewed, which leads eventually to a very high SNR when estimating based on the below-mode RMS. I'm realizing now that this is probably ideally a problem to be fixed on the recording side... But in any case, has anyone dealt with this before and do you have any suggestions for how to alter the SNR calculation to make it more useful?

Illustration of what I mean (comparing a junk ROI to a real cell, although the real cell is also somewhat clipped):

[kushalkolar]

I have seen infinite SNR before, couldn't figure out a fix: kushalkolar/mesmerize-viz#38 (comment)

[ethanbb]

Yes, I saw that and I have also seen infinite SNR in some of my data. Off the top of my head, I think it could happen if the mode of a pixel also happens to be the minimum value. Perhaps some special casing is needed (broader than just for the infinite case?)

…

________________________________ From: Kushal Kolar ***@***.***> Sent: Tuesday, January 28, 2025 4:18:13 AM To: flatironinstitute/CaImAn ***@***.***> Cc: Ethan Blackwood ***@***.***>; Author ***@***.***> Subject: Re: [flatironinstitute/CaImAn] Dealing with pixels clipped at 0 for SNR calculation (Discussion #1456) I have seen infinite SNR before, couldn't figure out a fix: kushalkolar/mesmerize-viz#38 (comment)<kushalkolar/mesmerize-viz#38 (comment)> — Reply to this email directly, view it on GitHub<#1456 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACEO4AMYAAQKGVZQ6FOKHSL2M3ZHLAVCNFSM6AAAAABV7JYMOGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCOJXHE2TCMY>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[ethanbb]

Should we not maybe model the noise as a Poisson distribution rather than Gaussian?? I can give it a try...

…

________________________________ From: Ethan Blackwood ***@***.***> Sent: Tuesday, January 28, 2025 5:20:12 AM To: flatironinstitute/CaImAn ***@***.***>; flatironinstitute/CaImAn ***@***.***> Cc: Author ***@***.***> Subject: Re: [flatironinstitute/CaImAn] Dealing with pixels clipped at 0 for SNR calculation (Discussion #1456) Yes, I saw that and I have also seen infinite SNR in some of my data. Off the top of my head, I think it could happen if the mode of a pixel also happens to be the minimum value. Perhaps some special casing is needed (broader than just for the infinite case?)

________________________________ From: Kushal Kolar ***@***.***> Sent: Tuesday, January 28, 2025 4:18:13 AM To: flatironinstitute/CaImAn ***@***.***> Cc: Ethan Blackwood ***@***.***>; Author ***@***.***> Subject: Re: [flatironinstitute/CaImAn] Dealing with pixels clipped at 0 for SNR calculation (Discussion #1456) I have seen infinite SNR before, couldn't figure out a fix: kushalkolar/mesmerize-viz#38 (comment)<kushalkolar/mesmerize-viz#38 (comment)> — Reply to this email directly, view it on GitHub<#1456 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ACEO4AMYAAQKGVZQ6FOKHSL2M3ZHLAVCNFSM6AAAAABV7JYMOGVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTCOJXHE2TCMY>. You are receiving this because you authored the thread.Message ID: ***@***.***>

[kushalkolar]

Good question and I looked into this a few months ago! it's more than just poisson, more factors than the sensor's photon detection events: https://ieeexplore.ieee.org/document/8327626

[kushalkolar]

See the "3) Mixed Poisson-Gaussian Noise:"

[ethanbb]

Interesting, thanks. However, I'm trying to figure out whether modeling readout noise as a Gaussian makes sense/does this not take into account saturation? I might need to look further into how the PMT amplifier circuit and/or acquisition code works to produce values from 0 to 2^16-1. Maybe the 0 value does not actually correspond to "no photons" but rather is the result of rescaling the mean "no photons" level + negative readout noise.

In any case, maybe something can be done to take the Poisson component into account (although probably the better thing to do is just record with higher power).

[ethanbb]

OK here's one idea I just tested that might work specifically for getting an SNR after CNMF is done. If we trust the separation of YrA and C, we can just fit a distribution to all values of YrA rather than just the negative ones. Then we would use the CDF of that distribution to get the SNR for C + YrA the same way it's currently done with a Gaussian distribution.

Based on that paper, the most correct distribution to model the noise may be a Gaussian convolved with a Poisson. But that would be very difficult to work with, so I'm trying just a Gamma distribution instead, and on my first 2 test cases it seems to work pretty well...

from scipy.stats import gamma

cell_ids = [1170, 976]
labels = ['Real cell', 'Noise ROI']
fig, axss = plt.subplots(len(cell_ids), 2)

for axs, cell_id, label in zip(axss, cell_ids, labels):
    # fit gamma distribution
    noise = estimates.YrA[cell_id, :]
    gamma_params = gamma.fit(noise)

    # plot PDF with model fit
    ps, bins = axs[0].hist(noise, 100, density=True, label='noise empirical PDF')[:2]
    axs[0].plot(bins, gamma.pdf(bins, *gamma_params), label='Gamma distribution fit', lw=2)
    axs[0].legend()
    axs[0].set_title(label)

    # plot CDF with model fit
    probs = ps * np.diff(bins)
    cum_probs = np.insert(np.cumsum(probs), 0, 0)
    axs[1].plot(bins, cum_probs, label='noise empirical CDF', lw=2)
    axs[1].plot(bins, gamma.cdf(bins, *gamma_params), label='Gamma distribution fit', lw=2)
    axs[1].legend()
    axs[1].set_title(label)
plt.show()

I'm now going to try writing a function for the full SNR calculation and see how it works/how long it takes to run on a full set of cells.

[ethanbb]

Promising I think! (I switched to fitting a loggamma distribution because I was getting some degenerate fits with tiny shape parameters, and I saw on Wikipedia the suggestion to use loggamma instead to get around this.)

[ethanbb]

OK, I plugged the gamma-based SNR into mesmerize-viz to see if it fixes classification of my pathologically skewed ROIs, and it seems to work great. It's part of my own repository rather than caiman right now, but here is the code in case you want to adapt it for use in caiman:

def compute_snr_gamma(C: np.ndarray, YrA: np.ndarray, use_loggamma=True,
                      remove_baseline=True, N=5, sigma_factor=3., dview=None) -> np.ndarray:
    """
    Compute an SNR measure based on fitting a gamma distribution to the residuals (YrA)
    of each component, then using the gamma ppf to evaluate exceptionality of events
    in the full traces (C + YrA).
    See evaluate_components, compute_event_exceptionality in components_evaluation.py
    
    Args:
        C: ndarray
            denoised traces (output of CNMF)
        
        YrA: ndarray
            residuals (output of CNMF)
        
        use_loggamma: bool
            whether to fit loggamma function instead of gamma (avoid underflow w/ small shape param)
        
        remove_baseline: bool
            whether to remove the baseline from YrA in a rolling fashion (8th percentile)
        
        N: int
            N number of consecutive events probability multiplied
        
        sigma_factor: float
            multiplicative factor for spread of gamma distribution (higher = higher criterion)
    """
    logger = logging.getLogger('caiman')

    if N == 0:
        # Without this, numpy ranged syntax does not work correctly, and also N=0 is conceptually incoherent
        raise Exception("FATAL: N=0 is not a valid value for compute_event_exceptionality()")
    
    if remove_baseline:
        logger.debug('Removing baseline from YrA')
        YrA = YrA - estimate_baseline(YrA, YrA.shape[1], slow_baseline=False)

    logger.debug('Computing event exceptionality with gamma distribution')
    args = zip(C, YrA, repeat(sigma_factor))
    logsf = np.empty_like(C)
    logsf_fn = get_loggamma_logsfs if use_loggamma else get_gamma_logsfs
    if dview is None:
        for out, arg_tuple in zip(logsf, args):
            out[:] = logsf_fn(arg_tuple)
    else:
        logger.info('SNR calculation in parallel')
        map_fn = dview.map if 'multiprocessing' in str(type(dview)) else dview.map_sync
        for out, res in zip(logsf, map_fn(logsf_fn, args)):
            out[:] = res
    
    # find most improbable sequence of N events using moving sum
    logsf_cum = np.cumsum(logsf, 1)
    logsf_seq = logsf_cum[:, N:] - logsf_cum[:, :-N]

    # select the minimum value of log-probability for each trace
    fitness = np.nanmin(logsf_seq, 1)
    # technically maybe weird to use normal distribution here, but
    # not sure how to deal with shape parameter if using gamma and 
    # I think it makes sense to use the same transformation to SNR values as regular SNR
    comp_SNR = -norm.ppf(np.exp(fitness / N))
    return comp_SNR


def get_gamma_logsfs(args: tuple[np.ndarray, np.ndarray, float]) -> np.ndarray:
    """
    Helper to compute the log of survival function (called "erf" in compute_event_exceptionality)
    for each sample in C_one + Yr_one based on fitting a gamma distribution to Yr_one.
    Inputs are vectors of C and YrA for a single ROI.
    """
    C_one, YrA_one, sigma_factor = args
    a, loc, scale = gamma.fit(YrA_one)
    scale *= sigma_factor  # artificially increase probabilities
    return gamma.logsf(C_one + YrA_one, a=a, loc=loc, scale=scale)


def get_loggamma_logsfs(args: tuple[np.ndarray, np.ndarray, float]) -> np.ndarray:
    """
    Like get_gamma_logsfs, but uses loggamma distribution instead.
    According to Wikipedia (https://en.wikipedia.org/wiki/Gamma_distribution#Caveat_for_small_shape_parameter),
    this can help avoid underflow.
    """
    C_one, YrA_one, sigma_factor = args
    # offset YrA so that it is in range, then take log
    # First scale by sigma_factor
    YrA_scaled = YrA_one * sigma_factor
    offset = 1 - np.min(YrA_scaled)
    YrA_offset = YrA_scaled + offset
    YrA_log = np.log(YrA_offset)

    # fit loggamma distribution
    params = loggamma.fit(YrA_log)

    # now apply same transformation (except scaling) to C + YrA and get survival function values
    CYrA_offset = YrA_one + C_one + offset
    CYrA_log = np.log(CYrA_offset)
    return loggamma.logsf(CYrA_log, *params)