Affine stitching algorithm

Affine stitching algorithm with iterative prediction

Support for affine models is experimental and is available in this branch: https://github.com/saalfeldlab/stitching-spark/tree/affine-stitching-with-uncertainty

The stitching algorithm takes several tile configuration files as an input, one for each channel. It is assumed that the tiles are aligned between channels, such that the tiles at the same grid position have the same index, stage coordinates, and affine transformation if any. All tiles are expected to be of the same size.

The core of the stitching algorithm is pairwise phase-correlation matching (based on Fiji stitching plugin) and global optimization (from mpicbg library).

These two main steps (phase correlation and global optimization) are performed iteratively in order to improve the resulting stitching quality and minimize the error while retaining as many tiles as possible.

Each iteration yields a solution to the stitching problem: a set of tiles with an estimated affine transformation for each of them along with estimation error statistics. The first iteration is intended to find a well-stitched subset of tiles, even though it is possible that many tiles are excluded from this set. Starting from the second iteration, the tiles are rematched where the matching is constrained by the predicted position and estimated uncertainty based on the previous solution. Two rematching modes are supported:

• full: all pairs are rematched
• incremental: only those pairs are rematched that were not included in the solution on the previous iteration. This optimization allows to save computation time but may give less satisfactory results in certain situations.

In context of pairwise matching, a pair of tiles is considered ordered where the first tile is fixed and the second is moving. Fixed tile box position is denoted as F, and moving tile box position is denoted as M (both represent local positions inside corresponding tiles).

---

First iteration:

1. Each tile is subdivided into 2x2x2 boxes. Pairwise matching will be done between these boxes and not the full tiles, this is to ensure that each tile has enough control points for affine model estimation.

2. Find pairs of overlapping boxes. For that, each tile box is transformed into the world space by either transforming its middle point (local coordinates in the tile space):

   • using the predefined affine transformation for the tile, or
   • using the stage coordinates of the tile

   Then, the intersection interval is calculated for each pair of boxes, and the pair is considered overlapping if:

   • the intersection interval is non-empty, and
   • at most one side of the intersection interval is shorter than half the length of the tile box side along that dimension (optional, enabled by default). This is expected to select only X/Y/Z adjacent pairs but not the corner ones.

3. Render both ROI images (tile boxes) in the relative coordinate space of the fixed box:

   • load tile image
   • apply flatfield correction (optional if flatfields exist in the same folder as the configuration file)
   • crop the ROI part (tile box)
   • transform into coordinate space of the fixed box. For the fixed box, it is just an identity transformation. For the moving box, the transformation is as follows:
     • moving tile box -> moving tile -> world -> fixed tile -> fixed box
   • average using all channels
   • apply Gaussian blur with given sigma for XY, and adjusted sigma for Z with respect to the pixel resolution (default is 2.0)

4. Perform pairwise matching between ROIs:

   • compute phase correlation
   • select best N peaks (default is 50)
   • test possible positions for each of them and pick the one with the highest cross correlation

5. Having the offset between the bounding box of the transformed moving tile box and the fixed box, find the global offset between the tile boxes:

   • find new transformation for the moving tile: (?)
     • find the offset between the position of the transformed moving tile box and its bounding box in the fixed box space
     • map the new offset position of the transformed moving tile box into the world space
     • calculate the offset between the new and the old positions of the moving tile box in the world space
     • calculate new translational component by adding the offset to the old translational component
     • create new transformation for the moving tile by merging the new translational component and the existing linear component
   • undo the linear component of the fixed tile transformation and map F into this coordinate space (F')
   • undo the linear component of the new moving tile transformation and map M into this coordinate space (M')
   • the new offset between the tile boxes will be M' - F'

6. Global optimization. Perform grid search on threshold value pairs (min.cross correlation; min.variance) within certain range and run the following for each combination:

   • create an identity affine model for each tile
   • filter pairwise matches using given threshold values
   • for each overlapping pair, calculate the full tile offset: M' - F' + F - M
   • for each overlapping pair, add two point matches:
     • middle point of the fixed tile box -> moving tile
     • middle point of the moving tile box -> fixed tile
   • randomly shift the points to ensure that point matches for any tile do not lie on the same plane
   • find and retain only the largest component of the tile configuration graph
   • replace some of the models to the translational model if there are not enough point matches to estimate the affine model
   • prealign and optimize the configuration using translational models for all tiles
   • optimize the configuration using affine models (for tiles with enough point matches)
   • collect min/max/avg displacement statistics
   • choose the resulting configuration with the largest number of tiles that has max.error below certain value

---

Subsequent iterations:

1. Create search radius estimator:

   • build a KD-tree for stitched tiles from the previous iteration based on their stage coordinates
   • set the neighborhood window size to 3x tile size in all dimensions

2. Estimate expected affine transformations for tiles (only for excluded tiles if the incremental mode is ON):

   • for each tile select stitched tiles from its neighborhood
   • average linear components of their affine transformations (L)
   • estimate translation by applying the expected offset to the stage position (T)
   • combine them in an 'inverse' way: A = LT (because the translational component is estimated in the 'offset' space)

3. Subdivide the tiles and find overlapping pairs (same as on the first iteration).

4. For each overlapping pair, estimate search radius:

   • estimate search radius for both tiles boxes in the pair:
     • select stitched tiles within the neighborhood window around the current tile
     • for each selected tile, get a pair (stage position; transformed position by undoing linear component)
     • calculate average offset
     • calculate covariance matrix
     • find its eigenvalues/eigenvectors which represents an error ellipse
   • estimate a combined search radius for the moving tile box (zero-centered):
     • calculate combined covariance matrix as a sum of the covariance matrices for individual tile boxes
     • find its eigenvalues/eigenvectors which represents an error ellipse

5. Build the transformation to map the error ellipse into the fixed box space:

   • moving box -> moving tile -> position in fixed box space -> bounding box position in fixed box space

6. For each overlapping pair, render both ROI images in the fixed box space (same as on the first iteration)

7. Perform pairwise matching between ROIs:

   • compute phase correlation (same as on the first iteration)
   • select best N peaks that fall within the search radius (default is 50)
   • test possible positions for each of them and pick the one with the highest cross correlation that falls within the search radius

   Whether a peak falls within the search radius is validated as follows:

   • test offset is the offset between the fixed tile box and the bounding box of the transformed moving box
   • the peak position is mapped to the unit sphere by applying the inverse transformation of the error ellipse
   • check if the mapped point falls within the unit sphere

8. Find the global offset between the tile boxes (same as on the first iteration)

9. Global optimization (same as on the first iteration)