Seeking Clarification on Calculation Details in local_correlations

[rcalfredson]

Hi Caiman developers and community,

I hope this message finds you all well! Firstly, I want to extend my appreciation to all the contributors of the Caiman project. Having discovered the project late last week, I've spent the whole week working with the tools and studying the theory supporting them, which has been both helpful and interesting.

Recently, while delving into the summary_images.local_correlations function, I noticed something intriguing about how the correlation image is calculated, especially when diagonal neighbors are considered, and wanted to hear your thoughts to ensure I'm on the right track.

In the process of calculating correlations with diagonal neighbors, two matrices, $\rho_{d1}$ and $\rho_{d2}$, are used, with each one intended to capture correlations in two distinct diagonal directions via matrix offsets when updating the accumulator.

Here's a succinct breakdown of my observation and queries:

Matrix $\rho_{d1}$

1. Intended to calculate the correlation of each point with its upper-right neighbor, so $\rho_{d1}[0,0]$ contains the correlation between pixels at $(1,0)$ and $(0,1)$ in a given image/video
2. Offsets applied during the accumulator update are:
   a. $(0,1)$, correlating $(i,j)$ with $(i+1,j-1)$.
   b. $(1,1)$, which appears to correlate each point's top and left neighbors, i.e., $(i−1,j)$ with $(i,j−1)$. I’m wondering if an offset of $(1,0)$ would align more accurately with the directional intent, correlating $(i,j)$ with $(i+1,j−1)$?

Matrix $\rho_{d2}$

1. Intended to calculate the correlation of each point with its lower-right neighbor, so $\rho_{d2}[0,0]$ contains the correlation between pixels at $(0,0)$ and $(1,1)$ in a given image/video
2. Offsets applied during the accumulator update are:
   a. $(0,0)$, correlating $(i,j)$ with $(i+1,j+1)$
   b. $(1,0)$, aiming to correlate each point's bottom and left neighbors, i.e., $(i -1,j)$ with $(i,j+1)$. Might it be more pertinent to use $(1,1)$ here to correlate $(i,j)$ with $(i−1,j−1)$?

In summary, with the current method, for a point $(i,j)$, the pairs of points being correlated and added to the accumulator appear to be:

• $(i,j)$ and $(i+1,j-1)$
• $(i,j)$ and $(i+1,j+1)$
• $(i−1,j)$ and $(i,j−1)$
• $(i-1,j)$ and $(i,j+1)$

In essence, whilst the existing schema does correlate the reference point with its upper-right and lower-right neighbors, it seems to simultaneously correlate two pairs of its immediate neighbors (top-and-left pair and bottom-and-left pair) rather than the reference point with all four diagonal directions, which, if my understanding serves me correctly, was the intent behind the implementation.

Please do enlighten me if there's something I might have misunderstood or overlooked. I'm eager to learn more about the inner workings of Caiman and to contribute in any way that I can.

Thank you so much for your time and for fostering such a collaborative environment within the Caiman community. I'm looking forward to your insights.

Sincerely,

Robert Alfredson

[EricThomson]

@rcalfredson thanks for the question. I have to admit I haven't dug into the details of this function. It definitely seems worth taking some time to work through it as you are doing. At best we should better document this function.

I would need to take some time to work through this before I have a reply. @j-friedrich not sure if you have any thoughts here.

I should add that local_correlations_fft() is more vetted code (for instance for 1p initialization it is what we use), so I would recommend focusing more on that than local_correlations() if you want the most "canonical" up-to-date working code (though looking through the code based, both functions are still used, so it is worth getting this sorted).

[j-friedrich]

@rcalfredson you are correct, $\rho_{d1}$ and $\rho_{d2}$ need to be swapped in lines 257 and 259. Thanks for the heads up. I personally use a PyTorch version of this function (it's ~10x faster on GPU than either of the CaImAn implementations) but haven't yet noticed the bug. If you run the corrected local_correlations and local_correlations_fft you'll see that the results now agree within numerical precision. (Someone with time on their hands could add this as a unit test)

Here the corrected code:

def local_correlations(
    mov: np.ndarray,
    eight_neighbours: bool = True,
    device: str = "cuda" if torch.cuda.is_available() else "cpu",
) -> np.ndarray:
    """Computes the correlation image for the input dataset mov

    Parameters
    ----------
    mov: ndarray
        Input movie data in 3D format.
    eight_neighbours: bool
        Use 8 neighbors if true, and 4 if false.
    device: str
        'cuda' or 'cpu', default is 'cuda' if GPU is available.

    Returns
    -------
    rho: ndarray
        Cross-correlation with adjacent pixels.
    """
    Y = torch.tensor(mov, dtype=torch.float32, device=device)
    rho = torch.zeros(Y.shape[1:], device=device)
    w_mov = (Y - torch.mean(Y, axis=0)) / (
        torch.std(Y, axis=0, correction=0) + torch.finfo(torch.float32).eps
    )

    rho_h = torch.mean(
        torch.multiply(w_mov[:, :-1, :], w_mov[:, 1:, :]), axis=0
    )
    rho_w = torch.mean(
        torch.multiply(w_mov[:, :, :-1], w_mov[:, :, 1:]), axis=0
    )

    rho[:-1, :] += rho_h
    rho[1:, :] += rho_h
    rho[:, :-1] += rho_w
    rho[:, 1:] += rho_w

    if eight_neighbours:
        rho_d1 = torch.mean(
            torch.multiply(w_mov[:, 1:, :-1], w_mov[:, :-1, 1:]), axis=0
        )
        rho_d2 = torch.mean(
            torch.multiply(w_mov[:, :-1, :-1], w_mov[:, 1:, 1:,]), axis=0
        )

        rho[1:, :-1] += rho_d1
        rho[:-1, 1:] += rho_d1
        rho[:-1, :-1] += rho_d2
        rho[1:, 1:] += rho_d2

        neighbors = 8 * torch.ones(Y.shape[1:3], device=device)
        neighbors[0, :] -= 3
        neighbors[-1, :] -= 3
        neighbors[:, 0] -= 3
        neighbors[:, -1] -= 3
        neighbors[0, 0] += 1
        neighbors[-1, -1] += 1
        neighbors[-1, 0] += 1
        neighbors[0, -1] += 1
    else:
        neighbors = 4 * torch.ones(Y.shape[1:3], device=device)
        neighbors[0, :] -= 1
        neighbors[-1, :] -= 1
        neighbors[:, 0] -= 1
        neighbors[:, -1] -= 1

    rho /= neighbors

    return rho.cpu().numpy()

[EricThomson]

Thanks a lot for noticing this @rcalfredson and for checking in on it @j-friedrich . I will do a PR on it, and add a test!

[rcalfredson]

Thank you both for being so responsive to my message- it’s much appreciated! I will also try the FFT variant as well of the local correlation calculation- had missed that one before.