Skip to content

Estimation API Reference

This page documents the parameter estimation functions.

compute_conversion_gain

compute_conversion_gain(movie, count_weight_gamma=0.2)

Calculate photon sensitivity and zero level from temporal variance analysis.

This function estimates camera parameters by fitting the noise transfer function from temporal variance. It uses HuberRegressor to robustly fit the relationship between mean signal and variance.

Parameters:

Name Type Description Default
movie ndarray

A movie in the format (time, height, width).

required
count_weight_gamma float

Weighting exponent for pixel counts in regression, by default 0.2. - 0.0: weigh each intensity level equally - 1.0: weigh each intensity in proportion to pixel counts

0.2

Returns:

Type Description
dict

ConversionGainConfig

Raises:

Type Description
AssertionError

If movie is not 3-dimensional or if insufficient intensity range is present.

Source code in src/anscombe_transform/estimate.py
def compute_conversion_gain(movie: np.array, count_weight_gamma: float = 0.2) -> ConversionGainConfig:
    """
    Calculate photon sensitivity and zero level from temporal variance analysis.

    This function estimates camera parameters by fitting the noise transfer function
    from temporal variance. It uses HuberRegressor to robustly fit the relationship
    between mean signal and variance.

    Parameters
    ----------
    movie : np.ndarray
        A movie in the format (time, height, width).
    count_weight_gamma : float, optional
        Weighting exponent for pixel counts in regression, by default 0.2.
        - 0.0: weigh each intensity level equally
        - 1.0: weigh each intensity in proportion to pixel counts

    Returns
    -------
    dict
        `ConversionGainConfig`

    Raises
    ------
    AssertionError
        If movie is not 3-dimensional or if insufficient intensity range is present.
    """
    assert movie.ndim == 3, (
        f"Thee dimensions (Time x Height x Width) of grayscale movie expected, got {movie.ndim} dimensions"
    )

    # assume that negative values are due to noise
    movie = np.maximum(0, movie.astype(np.int32, copy=False))
    intensity = (movie[:-1, :, :] + movie[1:, :, :] + 1) // 2
    difference = movie[:-1, :, :].astype(np.float32) - movie[1:, :, :]

    select = intensity > 0  # discard non-positive values
    intensity = intensity[select]
    difference = difference[select]

    counts = np.bincount(intensity.flatten())
    bins = _longest_run(
        counts > 0.01 * counts.mean()
    )  # consider only bins with at least 1% of mean counts
    bins = slice(max(bins.stop * 3 // 100, bins.start), bins.stop)
    assert bins.stop - bins.start > 100, (
        "The image does not have a sufficient range of intensities to compute the noise transfer function."
    )

    counts = counts[bins]
    idx = (intensity >= bins.start) & (intensity < bins.stop)
    variance = (
        np.bincount(
            intensity[idx] - bins.start,
            weights=(difference[idx] ** 2) / 2,
        )
        / counts
    )
    model = Regressor()
    model.fit(np.c_[bins], variance, counts**count_weight_gamma)
    conversion_gain = model.coef_[0]
    zero_level = -model.intercept_ / model.coef_[0]

    return {
        "model": model,
        "counts":  counts,
        "min_intensity": bins.start,
        "max_intensity" : bins.stop,
        "variance": variance,
        "conversion_gain": conversion_gain,
        "zero_level": zero_level,
    }