def _ts_value(position, counts, bkg, model, C_0_map):
    """
    Compute TS value at a given pixel position using the approach described
    in Stewart (2009).

    Parameters
    ----------
    position : tuple
        Pixel position.
    counts : `~numpy.ndarray`
        Count map.
    bkg : `~numpy.ndarray`
        Background map.
    model : `~numpy.ndarray`
        Source model map.

    Returns
    -------
    TS : float
        TS value at the given pixel position.
    """
    extract_fn = _collect_wrapper(extract_large_array)
    truncate_fn = _collect_wrapper(extract_small_array)

    # Get data slices
    counts_slice = extract_fn(counts, model, position)
    bkg_slice = extract_fn(bkg, model, position)
    C_0_slice = extract_fn(C_0_map, model, position)
    model_slice = truncate_fn(model, counts, position)

    # Flattened Arrays
    counts_ = np.concatenate([t.flat for t in counts_slice])
    bkg_ = np.concatenate([t.flat for t in bkg_slice])
    model_ = np.concatenate([t.flat for t in model_slice])
    C_0_ = np.concatenate([t.flat for t in C_0_slice])
    C_0 = np.sum(C_0_)

    root_fn = _sum_wrapper(_f_cash_root)
    amplitude, niter = _root_amplitude_brentq(counts_, bkg_, model_,
                                              root_fn=_f_cash_root)

    if niter > MAX_NITER:
        print('Exceeded maximum number of function evaluations!')
        return np.nan, amplitude, niter

    with np.errstate(invalid='ignore', divide='ignore'):
        C_1 = f_cash_sum(amplitude, counts_, bkg_, model_)

    # Compute and return TS value
    return (C_0 - C_1) * np.sign(amplitude), amplitude, niter