6. astec_drift

This method aims at estimating the drift (motion) that may exists between any two successive acquisitions. It has been developed to handle embryo motion, that may be due to the stage rotation (see section Section 5.2). This motion may then occur between the two stacks (stack #0 and stack #1) of the same time point, or between two time points.

Given the fusion of one stack (based on two views, ie left and right), this method compute the co-registration between any two successive time points. It somehow acts similarly to astec_intraregistration, but will test a number of initial rotations (spanned over the space of rotations) to find the transformations that best superimposes two successive acquisitions: astec_intraregistration assumed that two successive acquisitions are spatially close, so that the identity will be sufficient as the only initial transformation to retrieve the searched one.

To do so, a first co-registration is done between two successive images, the found transformation is assessed by an ad-hoc score, and, based on this score, it is decided to test other initial transformations if necessary.

Section 4.3 presents a tutorial for a drift-compensated fusion.

A comprehensive list of drift parameters can be found in section Section 16.11.

6.1. Registration score

The score of the transformation \(T_{t \leftarrow t+1}\) (that allows to resample \(I_{t}\) into the frame than \(I_{t+1}\)) is computed by

\[S(T_{t \leftarrow t+1}) = \max \left( I_{t} \circ T_{t \leftarrow t+1}, I_{t+1} \oplus B(r) \right) - I_{t+1} \oplus B(r)\]

For a perfect registration (ie the right transformation), the membranes of \(I_{t} \circ T_{t \leftarrow t+1}\) (\(I_{t}\) resampled into the frame of \(I_{t+1}\)) should superimpose on those of \(I_{t+1}\), thus \(\max \left( I_{t} \circ T_{t \leftarrow t+1}, I_{t+1} \right)\) should be close to \(I_{t+1}\), and the difference should be small.

To handle image differences and registration imperfections, a dilated version of \(I_{t+1}\) is used. The size \(r\) of the structuring element \(B(r)\) (parameter dilation_radius of astec_drift parameters, see section 16.11) is related to the range of expected spatial differences between \(I_{t} \circ T_{t \leftarrow t+1}\) and \(I_{t+1}\).

Figure 6.1 illustrates the score principle. On the left, when \(I_{t} \circ T_{t \leftarrow t+1}\) is below \(I_{t+1} \oplus B(r)\), the score is zero. When \(I_{t} \circ T_{t \leftarrow t+1}\) is above \(I_{t+1} \oplus B(r)\), the score is high when the membrane (the high values of \(I_{t} \circ T_{t \leftarrow t+1}\)) is far from a membrane of \(I_{t+1}\).

_images/drift-score-scheme.png

Fig. 6.1 Red: \(I_{t+1}\); dashed red: \(I_{t+1} \oplus B(r)\); green: \(I_{t} \circ T_{t \leftarrow t+1}\).

\(S(T_{t \leftarrow t+1})\) above defines an image (there is one value per voxel). To end up with a scalar value, the average of the higher values of \(I_{t} \circ T_{t \leftarrow t+1}\) (where membranes are likely to be) is computed. The percentile of points used for the average is defined by the parameter threshold_quantile (see section 16.11).

6.2. Drift prerequisite

It is assumed that the stack for which drifts have to be computed is already fused (see section 5) as exemplified in the in the /path/to/experiment/ directory as depicted below.

$ /path/to/experiment/
├── FUSE/
│  └── FUSE_stack0/
│     ├── 241010-Zidane_fuse_t000.mha
│     ├── 241010-Zidane_fuse_t001.mha
.  .  .

Note

Since this fusion will only serve to estimate the drifts, it can be computed at a low resolution (eg target-resolution = 0.60).

6.3. Drift / output data

The following parameter file allows to compute the transformations between the couple of successive time points.

PATH_EMBRYO = "."
EN = "241010-Zidane"
begin = 0
end = 199

EXP_FUSE = 'stack0'

movie_fusion_images = True
xy_movie_fusion_images = True
yz_movie_fusion_images = True

resolution = 0.60
template_type = 'FUSION'
template_threshold = 140

EXP_DRIFT = 'stack0'
rotation_sphere_radius = 4.2

Most parameters (see section 16.11) are the same than the ones of astec_intraregistration (see section 16.12 and section 11) since the principle of astec_drift is similar to that of astec_intraregistration.

After computation, a DRIFT directory is created.

/path/to/experiment/
├── ...
├── DRIFT/
│  └── DRIFT_<EXP_DRIFT>/
│     ├── CORRECTED0-CO-REGISTERED/
│     │  ├── 241010-Zidane_flo(t)_in_ref(t+1).mha
│     │  .
│     ├── CORRECTED0-CO-SCORE/
│     │  ├── 241010-Zidane_score_flo(t)_ref(t+).py
│     │  .
│     ├── CORRECTED0-CO-TRSFS/
│     │  ├── 241010-Zidane_drift_flo(t)_ref(t+1).trsf
│     │  .
│     ├── ITER0-CO-REGISTERED/
│     │  ├── 241010-Zidane_flo(t)_in_ref(t+1).mha
│     │  .
│     ├── ITER0-CO-SCORE/
│     │  ├── 241010-Zidane_score_flo(t)_ref(t+).py
│     │  .
│     │  └── figure_iter0_coregistration_analyze.py
│     ├── ITER0-CO-TRSFS
│     │  ├── 241010-Zidane_drift_flo(t)_ref(t+1).trsf
│     │  .
│     ├── ITER0-FUSE
│     │  └── FUSE_stack0/
│     │     ├── 241010-Zidane_drift_fuse_t(t).mha.gz
│     │     .
│     ├── ITER0-MOVIES_t007-010
│     │  └── FUSE/
│     │     └── FUSE_stack0/
│     │        ├── 241010-Zidane_drift_fuse_t000-199_xy0281.mha.gz
│     │        .
│     ├── ITER0-TRSFS_t000-199
│     │  ├── 241010-Zidane_drift_t(t).trsf
│     │  .
│     ├── ITER1-CO-REGISTERED
│     ├── ITER1-CO-SCORE
│     ├── ITER1-CO-TRSFS
│     ├── ITER1-FUSE
│     ├── ITER1-MOVIES_t007-010
│     ├── ITER1-TRSFS_t007-010
│     └── LOGS

Subdirectories prefixed by ITER<i> are results issued from the ith iteration. These directories are similar to the ones generated by astec_intraregistration

  • ITER<i>-CO-TRSFS/ contains the (rigid) transformations 241010-Zidane_drift_flo(t)_ref(t+1).trsf that allow to resampled the image at time t onto the image at time t+1 (see also section section 11.5).

  • ITER<i>-CO-REGISTERED/ contains the images at time t resampled onto the image at time t+1. It is kept for the score computation, and for visual assessment of the registration quality (to be compared for the image at time t+1 from the fusion subdirection, ie FUSE/FUSE_<EXP_FUSE>). Once the drift is found satisfactory, these sub-directories can be deleted.

  • ITER<i>-CO-SCORE/ contains the coregistration scores, each of them being saved into a single. More interestingly, a file figure_iter0_coregistration_analyze.py is also generated, that allows to generate a figure (see figure 6.2), where the co-registration scores can be assessed. In addition, a list of corrections (co-registration to be re-calculated) is proposed based on a threshold computed by Otsu’s method.

  • ITER<i>-TRSFS_t<begin>-<end> contains the transformations 241010-Zidane_drift_t(t).trsf that allow to resampled the image at time t onto a common template (see also section 11.6). The still image is defined by the parameter reference_index.

    Important

    This subdirectory contains the transformations that will be used for drift-compensated fusion.

  • ITER<i>-FUSE/ contains the fusion images after resampling by the transformations of ITER<i>-TRSFS_t<begin>-<end>. Visually comparing any two successive images is also a means to visually assess their co-registration (note that all images here have the same geometry) (see also section 11.7).

  • ITER<i>-MOVIES_t<begin>-<end> contains 2D+t movies built from ITER<i>-FUSE/. Like the above figure, these movies are a precious means to assess whether the co-resgitration were suceessful (see also section 11.8).

_images/drift-score-example-iter0.png

Fig. 6.2 Figure generated by execution of figure_iter0_coregistration_analyze.py obtained at the initial calculation of astec_drift. Top: scores with respect to time; bottom: rotation angle with respect to time. Otsu’s method is used to compute a threshold to propose a list of corrections (ie co-registrations to be re-computed).

Note

The figures generated by figure_iter(i)_coregistration_analyze.py as well as the movies in ITER<i>-MOVIES_t<begin>-<end> allows to visually assess whether the co-registrations perform well. They have to serve as cues to decide whether the drift have been successfully corrected along the whole sequence, or whether an other iteration of corrections have to be done.

6.4. Drift / principle

To the opposite of the other modules of astec, astec_drift can be run iteratively within the same experiment (with the same value of EXP_DRIFT)

At the first run of astec_drift, the ITER0-* subdirectories are created, to initialize the drift computation process.

Then, the first run and the other runs have the same behavior. Based on the results of the previous (say the ith) iteration (in ITER<i>-* subdirectories), some co-registrations are recomputed (results will be stored in CORRECTED(i)-* subdirectories), then these corrections and the unchanged are pulled together in the ITER(i+1)-* subdirectories.

Note

The transformations used for drift-compensated fusion are the most recent ones, ie the ones associated with the latest iteration.

6.5. Drift / selection of corrections

At the first run, a threshold is computed from the scores collected at initialization (in the ITER0-CO-SCORE/ sub-directory), as illustrated in figure 6.2. Assessment of the computed corrections can be then either done with the same figure (in the ITER1-CO-SCORE/ sub-directory then) or by the movies (in the ITER1-MOVIES_t<begin>-<end>/ sub-directory then) generated at the next iteration.

_images/drift-score-example-iter1.png

Fig. 6.3 Figure generated by execution of figure_iter1_coregistration_analyze.py obtained after a first round of corrections. Top: scores with respect to time; bottom: rotation angle with respect to time. Otsu’s method is used to compute a threshold to propose a list of corrections (ie co-registrations to be re-computed) for the next iteration.

After visual assessment of the re-computed co-registrations (at a given iteration), a number of additional corrections can be identified for a next iteration. Using then the threshold computed from the collected scores may not be optimal, since too many couples of co-registrations may be selected. From figure 6.3, it can be seen that the threshold computed from results of the first iteration will be 3.35, and will thus select couples of successive images that arealready correctly co-registered.

Selection can then be done (see section 16.11 for details)

  • either by specifying a threshold (with the parameter score_threshold)

  • or by specyfing the time points to be corrected (with parameters corrections_to_be_done or corrections_to_be_added) acquisition time points to be corrected (even if the corresponding score is below the threshold), in addition to the time points above the computed or given threshold. For ‘i’ in the list, it means that the ith image is not correctly registered with the (i+1)th one.

6.6. Drift / computation of corrections

The registration procedure is based on an iterative minimization [OurselinRPA00], and thus required an initial transformation to start with, usually chosen as the identity. Because of the tesselated nature of embryo, two images are unlikely to be correctly co-registered if they are too far apart (in terms of transformation).

Identifying the correct transformation comes then to minimize the co-registration score over the set of possible transformations. Since the registration procedure [OurselinRPA00] is able to optimize the registration locally (from an initial transformation), it is sufficient to discretize the rotation space to define a set of initial transformations.

The 3D rotation space is defined by a 3D sphere (of radius \(pi\)), and the set of initial rotations is defined by a 3D discrete sphere (whose radius is defined by the rotation_sphere_radius parameter). Choosing a smaller value of rotation_sphere_radius will decrease the number of rotations to be tested, thus decreases the computational cost (at the risk of missing the right convergence basin).

6.7. Drift / inter-stack registration

astec_drift is first designed to compute the transformations between any two successive fused images of a given stack (obtained by fusion of the left and right images of a given angle). These transformations can be then used for a drift-compensated fusion. However, it assumes that the two drift-compensated series of stacks are already (almost) co-registered.

The two stacks may have to be co-registered. Since the drift compensation aimed at stabilizing the stack time series, co-registering the two stacks comes to co-register one time point of the two stacks.

Recall that the transformations in ITER<i>-TRSFS_t<begin>-<end>/ are defined with respect to a still/reference image defined by the parameter reference_index.

Important

The parameter reference_index has to be the same for the two stacks to be processed. If not specified in the parameter file, it is set to the begin (first time point).

To co-register the two stacks, it is sufficient to specify the two fusion directories (one per stack) in the <EXP_FUSE> parameter, as in

PATH_EMBRYO = "."
EN = "241010-Zidane"
begin = 0
end = 0

EXP_FUSE = ['stack1', 'stack0']
score_threshold = 11.00
...

  • It is mandatory to specify the two directories in the second stack (stack #1) drift parameter file, and to put the fusion directory of the first stack (stack #0) after the the fusion directory of the second stack.

  • end has been set to the same value than begin. This way, only the stack-to-stack drift will be computed. If end has a larger value (say end = 199), a new iteration of intra-stack drift estimation (for the second stack, ie stack #1) will be performed.

  • score_threshold can be set to the threshold computed at the first run. It will allow an earlier stop when testing the set of initial rotations.

This will create a new subdirectory named CO-STACK/ into DRIFT_<EXP_DRIFT>/ containing the fused image of the second stack (stack #1) at time reference_index resampled in the frame of the corresponding image of the first stack (stack #0), so that the co-registration can be visually assessed, as well as the transformation for this resampling, which will be used as the transformation between the two stacks.

/path/to/experiment/
├── ...
├── DRIFT/
│  └── DRIFT_<EXP_DRIFT>/
│     ├── CO-STACK/
│     │  ├── 241010-Zidane_fuse_t003_costack.mha
│     │  └── 241010-Zidane_fuse_t003_costack.trsf

6.8. Drift / disk cleaning

Once the computed drift is satisfactory, which can be assessed by visual inspection of the 2D+t movies in the last ITER<i>-MOVIES_t<begin>-<end> directory, some directories or file can be removed to spare room on the physycal device.

  • The ITER<i>-CO-REGISTERED directories that contain the images at time t resampled onto the image at time t+1 for the ith iteration (see section 6.3).

  • The ITER<i>-FUSE/ directories that contain the fusion images after resampling by the transformations of ITER<i>-TRSFS_t<begin>-<end> (see section 6.3). These images serve to generate the 2D+t movies.

  • In addition, the image file template_t<begin>-<end>.mha, that serves as geometry template to generate the images of ITER<i>-FUSE/ can be at least compressed,