Chromatic Aberration Correction

[rfor10]

Nonrigid Image Transformation

Chromatic aberration causes light of different wavelengths to be imaged at different XY pixels in a systematic manner. This can be corrected after imaging with high accuracy, but to apply the correction requires non-rigid transformation of the image.

We correct this in the Wollman lab by imaging multicolored diffraction limited beads (100nm Tetraspeck) in all relevant wavelengths for MERFISH experiments. The distortion field is robustly estimated from nearest neighbor (between two channels) pairs of gaussian fitted centroids of these beads for each pair of image channels.

Candidate Solution

scipy.interpolate.interp2d

Interpolate images with coordinates modified from pixels to chromatic aberration shifted pixels back onto regular pixel grid. As a bonus one can also register translational stage drift/precision errors easily with the same computation.

def interp_warp(img, x, y):
    """
    Apply chromatic abberation shifts to images.

    Parameters
    ----------
    img : ndarray
    x : array
      Shifted pixel coordinates same size as img.shape[0]
    y : array
      Shifted pixel coordinates same size as img.shape[1]

    Returns
    -------
    nimg : ndarray
      Same size as img but interpolated from x, y onto 0, 1, ...
    """
    i2 = interpolate.interp2d(x, y, img)
    nimg = i2(range(img.shape[0]), range(img.shape[1]))
    return nimg

[dganguli]

Hi @rfor10 thanks for the issue! So far we have considered chromatic aberration correction out of scope of Starfish, but we are now starting to re-consider :)

For now, I have more questions than answers to make sure I both understand the problem and your proposed solution. If you help me work through this then we can figure out a plan for what to do next!

1. Do I understand the problem setup correctly? Let's say I have 3 color channels. I'm imaging multi-spectral beads, and I obtain 3 images of the beads, one per each acquisition at the relevant wavelength. I know the beads in each image should be in the same spatial position (since they inherently aren't moving spatially) but due to chromatic aberration, they are unfortunately in different positions in different channels.

2. How do you learn the distortion field? Naively, I can find the spatial positions of all the beads in all my channels independently (via Gaussian spot fitting, for example). Next, I can pick an anchor channel, lets say channel 1. Then, for channel 2 I can learn a mapping from channel 2->1 using, for example, a point-cloud registration algorithm. I can repeat this exercise to learn a mapping from channel 3->1.

3. OK, I've learned my distortion field, but now what? Naively, at this point, my usage of the point cloud-registration algorithm has returned to me two functions, f_2() and f_3(), which I can apply to point sets from channel 2 and 3 respectively to make them correspond spatially to the point set from channel 1.

4. OK, but I'm not interested in simply applying f()'s to point sets from the multi-spectral beads right? Rather, what I want is to apply f()'s to my hybridization images. I don't really care about aligning my beads across channels, I care about aligning my RNA molecules across channels.

5. OK, so how am I going to apply my learned f()'s to my hybridization images? Well, scipy.interpolate.interp2d could help me here. I think? I'm a little stuck here. But I think the intuition is something like, I can leverage interpolation of my shifted points to learn how to shift everything on a pixel grid. Perhaps, given the language I've described here, could you describe to me what exactly x, y, and img are?

Am I on the right track here?

[dganguli]

Follow up question, I think chromatic aberration separates nicely into two problems: 1. Learning the correction factor and 2. Applying the correction factor. Your issue only addresses 2 right? Should Starfish also implement a solution to 1?

[rfor10]

1. Do I understand the problem setup correctly? Let's say I have 3 color channels. I'm imaging multi-spectral beads, and I obtain 3 images of the beads, one per each acquisition at the relevant wavelength. I know the beads in each image should be in the same spatial position (since they inherently aren't moving spatially) but due to chromatic aberration, they are unfortunately in different positions in different channels.

Exactly

2. How do you learn the distortion field? Naively, I can find the spatial positions of all the beads in all my channels independently (via Gaussian spot fitting, for example). Next, I can pick an anchor channel, lets say channel 1. Then, for channel 2 I can learn a mapping from channel 2->1 using, for example, a point-cloud registration algorithm. I can repeat this exercise to learn a mapping from channel 3->1.

Yes conceptually - what I have done is just pair beads based on nearest neighbor to find delta xy for each bead from many FOVs yielding a sparse matrix of offset(imageCoordinateX, imageCoordinateY). With our optics the offset of X is relatively independent of Y and vice versa. So ultimately the output of this step are two vectors offsetX(imageCoordinateX) and offsetY(imageCoordinateY) which are interpolated from the beads to fill missing values. Perhaps a more general solution would be needed to account for complex relationship in XY of the distortion.

Note - this distortion field with our optics is simple and smooth, but not centered coincident with the center of the image. I.e. center of distortion != 1024, 1024 in our (2048, 2048) images.

3. OK, I've learned my distortion field, but now what? Naively, at this point, my usage of the point cloud-registration algorithm has returned to me two functions, f_2() and f_3(), which I can apply to point sets from channel 2 and 3 respectively to make them correspond spatially to the point set from channel 1.

Yes, one could imagine different approaches (correct all to some pseudo reference such as mean between all). We've been applying the correction to map onto a reference.

4. OK, but I'm not interested in simply applying f()'s to point sets from the multi-spectral beads right? Rather, what I want is to apply f()'s to my hybridization images. I don't really care about aligning my beads across channels, I care about aligning my RNA molecules across channels.

Yes - I mean they can be applied to check visually or quantitatively how well the correction is working, but that is not the end product/goal.

5. OK, so how am I going to apply my learned f()'s to my hybridization images? Well, scipy.interpolate.interp2d could help me here. I think? I'm a little stuck here. But I think the intuition is something like, I can leverage interpolation of my shifted points to learn how to shift everything on a pixel grid. Perhaps, given the language I've described here, could you describe to me what exactly x, y, and img are?

Am I on the right track here?

I think that is the idea. For example, our images are (2048, 2048) pixels. x = range(2048) + offSetX(range(2048)) + stageTranslationX (optional) y = range(2048) + offSetY(range(2048)) + stageTranslationY (optional) img = the spot intensity image one wishes to chromatically/translationally correct

So the intuition is that we are redefining the coordinates of this image to be aligned with the reference, and then interpolating back onto the 0 to 2048 coordinate values that the reference exists in naturally.

[rfor10]

Follow up question, I think chromatic aberration separates nicely into two problems: 1. Learning the correction factor and 2. Applying the correction factor. Your issue only addresses 2 right? Should Starfish also implement a solution to 1?

I think that for the near term any lab that wants to try an approach like this will have to be able to fit a bunch of beads and supply that as input (sort of an acid test). Long term for the democratization of the technique it would be good to have some stock routines available for fitting. It is mostly a one-time thing to fit the chromatic aberration, but has to be repeated whenever changes are made to the microscope.

After some more thought - I think starfish can probably implement three routines with the first two being lower priority. The priority statement is based on users of starfish who want to do chromatic aberration correction must already have that data available and it is relatively simple to do. However, if they want to use the starfish framework to process their data they need the ability to apply the correction.

1. Find candidate beads in each image. Essentially can be fitting gaussians to local maxima and thresholding on spot widths and amplitudes.

2. Fit a distortion model from sets of paired bead candidates. Reject local outliers from the bead pairs. E.g. one pair has a significantly different shift than its neighbors. Then fit whatever the input is for the chromatic aberration correction function.

3. Apply chromatic aberration model to unwarp/align images.

[dganguli]

Ah this clarifies things immensely, thank you!

Regarding this statement:

what I have done is just pair beads based on nearest neighbor to find delta xy for each bead from many FOVs yielding a sparse matrix of offset(imageCoordinateX, imageCoordinateY)

At the end of this operation, is the dimensionality of offset(imageCoordinateX, imageCoordinateY) = [num_beads+total, 2] ?

If that's true, then I'm confused by this:

x = range(2048) + offSetX(range(2048))

Which I think you resolve by saying this

which are interpolated from the beads to fill missing values

I *think what's going on is that your offSetX is a function, f() that you've learned from running a 1-D interpolation between x_reference_bead_position and x_reference_bead_position+delta_x such that x_reference_bead_position = f(x_reference_bead_position+delta_x). Now you can apply this function on a something of dimension range(2048), as implied above, such that it works not just for bead locations, but all possible pixel coordinate values.

Is that right? Sorry for my obtuseness.

[rfor10]

Ah this clarifies things immensely, thank you!

Regarding this statement:

what I have done is just pair beads based on nearest neighbor to find delta xy for each bead from many FOVs yielding a sparse matrix of offset(imageCoordinateX, imageCoordinateY)

At the end of this operation, is the dimensionality of offset(imageCoordinateX, imageCoordinateY) = [num_beads+total, 2] ?

If that's true, then I'm confused by this:

x = range(2048) + offSetX(range(2048))

Which I think you resolve by saying this

which are interpolated from the beads to fill missing values

I *think what's going on is that your offSetX is a function, f() that you've learned from running a 1-D interpolation between x_reference_bead_position and x_reference_bead_position+delta_x such that x_reference_bead_position = f(x_reference_bead_position+delta_x). Now you can apply this function on a something of dimension range(2048), as implied above, such that it works not just for bead locations, but all possible pixel coordinate values.

Is that right? Sorry for my obtuseness.

I think you have it figured out. Sorry for not writing it in a more clear manner. I will restate just to try to make sure we're on the same page.

the dimensionality of offset(imageCoordinateX, imageCoordinateY) = [num_beads+total, 2] I am not sure what the total is here, but it would be [num_beads, 2], and yes I am then interpolating these offsets to find two offset vectors (num_x_pixels, 1) and (num_y_pixels, 1). In practice we store these and rather than x = range(2048) + offSetX(range(2048)) it is actually just x = range(2048) + savedOffsetVectorX.