Skip to content

Register

Register Initial

view_register_shift_info

For all shifts to imaging rounds from the reference round computed in the register_initial section of the pipeline, this plots the values of the shifts found and the score compared to the score_thresh.

For each round, there will be 3 plots:

  • y shift vs x shift for all tiles
  • z shift vs x shift for all tiles
  • score vs score_thresh for all tiles (a green score = score_thresh line is plotted in this).

In each case, the markers in the plots are numbers. These numbers indicate the tile the shift was found for. The number will be blue if score > score_thresh and red otherwise.

Parameters:

Name Type Description Default
nb Notebook

Notebook containing at least the register_initial page.

 required
outlier bool

If True, will plot nb.register_initial.shift_outlier instead of nb.register_initial.shift. In this case, only tiles for which the two are different are plotted for each round.

 False
Source code in coppafish/plot/register/shift.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
def view_register_shift_info(nb: Notebook, outlier: bool = False):
    """
    For all shifts to imaging rounds from the reference round computed in the `register_initial` section
    of the pipeline, this plots the values of the shifts found and the `score` compared to
    the `score_thresh`.

    For each round, there will be 3 plots:

    * y shift vs x shift for all tiles
    * z shift vs x shift for all tiles
    * `score` vs `score_thresh` for all tiles (a green score = score_thresh line is plotted in this).

    In each case, the markers in the plots are numbers.
    These numbers indicate the tile the shift was found for.
    The number will be blue if `score > score_thresh` and red otherwise.

    Args:
        nb: Notebook containing at least the `register_initial` page.
        outlier: If `True`, will plot `nb.register_initial.shift_outlier` instead of
            `nb.register_initial.shift`. In this case, only tiles for which
            the two are different are plotted for each round.
    """
    shift_info = {}
    if nb.basic_info.is_3d:
        ndim = 3
    else:
        ndim = 2
    for r in nb.basic_info.use_rounds:
        name = f'Round {r}'
        shift_info[name] = {}
        shift_info[name]['tile'] = nb.basic_info.use_tiles
        if outlier:
            shift_info[name]['shift'] = nb.register_initial.shift_outlier[nb.basic_info.use_tiles, r, :ndim]
            shift_info[name]['score'] = nb.register_initial.shift_score_outlier[nb.basic_info.use_tiles, r]
        else:
            shift_info[name]['shift'] = nb.register_initial.shift[nb.basic_info.use_tiles, r, :ndim]
            shift_info[name]['score'] = nb.register_initial.shift_score[nb.basic_info.use_tiles, r]
        shift_info[name]['score_thresh'] = nb.register_initial.shift_score_thresh[nb.basic_info.use_tiles, r]

    if outlier:
        title_start = "Outlier "
    else:
        title_start = ""
    shift_info_plot(shift_info, f"{title_start}Shifts found in register_initial part of pipeline "
                                f"from round {nb.basic_info.ref_round}, channel "
                                f"{nb.basic_info.ref_channel} to channel "
                                f"{nb.register_initial.shift_channel} for each round and tile")

Function to plot results of exhaustive search to find shift between ref_round/ref_channel and round r, channel c for tile t. This shift will then be used as the starting point when running point cloud registration to find affine transform. Useful for debugging the register_initial section of the pipeline.

Parameters:

Name Type Description Default
nb Notebook

Notebook containing results of the experiment. Must contain find_spots page.

 required
t int

tile interested in.

 required
r int

Want to find the shift between the reference round and this round.

 required
c Optional[int]

Want to find the shift between the reference channel and this channel. If None, config['shift_channel'] will be used, as it is in the pipeline.

 None
return_shift bool

If True, will return shift found and will not call plt.show() otherwise will return None.

 False

Returns:

Type Description
Optional[np.ndarray]

best_shift - float [shift_y, shift_x, shift_z]. Best shift found. shift_z is in units of z-pixels.
Source code in coppafish/plot/register/shift.py 
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
def view_register_search(nb: Notebook, t: int, r: int, c: Optional[int] = None,
                         return_shift: bool = False) -> Optional[np.ndarray]:
    """
    Function to plot results of exhaustive search to find shift between `ref_round/ref_channel` and
    round `r`, channel `c` for tile `t`. This shift will then be used as the starting point when running point cloud
    registration to find affine transform.
    Useful for debugging the `register_initial` section of the pipeline.

    Args:
        nb: Notebook containing results of the experiment. Must contain `find_spots` page.
        t: tile interested in.
        r: Want to find the shift between the reference round and this round.
        c: Want to find the shift between the reference channel and this channel. If `None`, `config['shift_channel']`
            will be used, as it is in the pipeline.
        return_shift: If True, will return shift found and will not call plt.show() otherwise will return None.

    Returns:
        `best_shift` - `float [shift_y, shift_x, shift_z]`.
            Best shift found. `shift_z` is in units of z-pixels.
    """
    config = nb.get_config()['register_initial']
    if c is None:
        c = config['shift_channel']
        if c is None:
            c = nb.basic_info.ref_channel
    if not np.isin(c, nb.basic_info.use_channels):
        raise ValueError(f"c should be in nb.basic_info.use_channels, but value given is\n"
                         f"{c} which is not in use_channels = {nb.basic_info.use_channels}.")

    coords = ['y', 'x', 'z']
    shifts = [{}]
    for i in range(len(coords)):
        shifts[0][coords[i]] = np.arange(config['shift_min'][i],
                                         config['shift_max'][i] +
                                         config['shift_step'][i] / 2, config['shift_step'][i]).astype(int)
    if not nb.basic_info.is_3d:
        config['nz_collapse'] = None
        config['shift_widen'][2] = 0  # so don't look for shifts in z direction
        config['shift_max_range'][2] = 0
        shifts[0]['z'] = np.array([0], dtype=int)
    shifts = shifts * nb.basic_info.n_rounds  # get one set of shifts for each round
    c_ref = nb.basic_info.ref_channel
    r_ref = nb.basic_info.ref_round
    # to convert z coordinate units to xy pixels when calculating distance to nearest neighbours
    z_scale = nb.basic_info.pixel_size_z / nb.basic_info.pixel_size_xy
    print(f'Finding shift between round {r_ref}, channel {c_ref} to round {r}, channel {c} for tile {t}')
    shift, shift_score, shift_score_thresh, debug_info = \
        compute_shift(spot_yxz(nb.find_spots.spot_details, t, r_ref, c_ref),
                      spot_yxz(nb.find_spots.spot_details, t, r, c),
                      config['shift_score_thresh'], config['shift_score_thresh_multiplier'],
                      config['shift_score_thresh_min_dist'], config['shift_score_thresh_max_dist'],
                      config['neighb_dist_thresh'], shifts[r]['y'], shifts[r]['x'], shifts[r]['z'],
                      config['shift_widen'], config['shift_max_range'], z_scale,
                      config['nz_collapse'], config['shift_step'][2])
    title = f'Shift between r={r_ref}, c={c_ref} and r={r}, c={c} for tile {t}. YXZ Shift = {shift}.'
    if return_shift:
        show = False
    else:
        show = True
    view_shifts(debug_info['shifts_2d'], debug_info['scores_2d'], debug_info['shifts_3d'],
                debug_info['scores_3d'], shift, debug_info['min_score_2d'], debug_info['shift_2d_initial'],
                shift_score_thresh, debug_info['shift_thresh'], config['shift_score_thresh_min_dist'],
                config['shift_score_thresh_max_dist'], title, show)
    if return_shift:
        return shift

Register

scale_box_plots

Function to plot distribution of chromatic aberration scaling amongst tiles for each round and channel. Want very similar values for a given channel across all tiles and rounds for each dimension. Also expect \(y\) and \(x\) scaling to be very similar. \(z\) scaling different due to unit conversion.

Parameters:

Name Type Description Default
nb Notebook

Notebook containing the register and register_debug NotebookPages.

 required
Source code in coppafish/plot/register/diagnostics.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def scale_box_plots(nb: Notebook):
    """
    Function to plot distribution of chromatic aberration scaling amongst tiles for each round and channel.
    Want very similar values for a given channel across all tiles and rounds for each dimension.
    Also expect $y$ and $x$ scaling to be very similar. $z$ scaling different due to unit conversion.

    Args:
        nb: *Notebook* containing the `register` and `register_debug` *NotebookPages*.
    """
    if nb.basic_info.is_3d:
        ndim = 3
        if np.ptp(nb.register.transform[np.ix_(nb.basic_info.use_tiles, nb.basic_info.use_rounds,
                                               nb.basic_info.use_channels)][:, :, :, 2, 2]) < 1e-5:
            ndim = 2
            warnings.warn("Not showing z-scaling as all are the same")
    else:
        ndim = 2

    fig, ax = plt.subplots(ndim, figsize=(10, 6), sharex=True)
    ax[0].get_shared_y_axes().join(ax[0], ax[1])
    fig.subplots_adjust(left=0.075, right=0.95, bottom=0.15)
    y_titles = ["Scaling - Y", "Scaling - X", "Scaling - Z"]
    n_use_channels = len(nb.basic_info.use_channels)
    # different boxplot color for each channel
    # Must be distinct from black and white
    channel_colors = distinctipy.get_colors(n_use_channels, [(0, 0, 0), (1, 1, 1)])
    for i in range(ndim):
        box_data = [nb.register.transform[nb.basic_info.use_tiles, r, c, i, i] for c in nb.basic_info.use_channels
                    for r in nb.basic_info.use_rounds]
        bp = ax[i].boxplot(box_data, notch=0, sym='+', patch_artist=True)
        leg_markers = []
        c = -1
        for j in range(len(box_data)):
            if j % n_use_channels == 0:
                c += 1
                leg_markers = leg_markers + [bp['boxes'][j]]
            bp['boxes'][j].set_facecolor(channel_colors[c])
        ax[i].set_ylabel(y_titles[i])

        if i == ndim-1:
            tick_labels = np.tile(nb.basic_info.use_rounds, n_use_channels).tolist()
            leg_labels = nb.basic_info.use_channels
            ax[i].set_xticks(np.arange(len(tick_labels)))
            ax[i].set_xticklabels(tick_labels)
            ax[i].legend(leg_markers, leg_labels, title='Channel')
            ax[i].set_xlabel('Round')
    ax[0].set_title('Boxplots showing distribution of scalings due to\nchromatic aberration amongst tiles for each '
                    'round and channel')
    plt.show()

view_affine_shift_info

For all affine transforms to imaging rounds/channels from the reference round computed in the register section of the pipeline, this plots the values of the shifts, n_matches (number of neighbours found) and error (average distance between neighbours).

For each round and channel (channel is changed by scrolling with the mouse), there will be 3 plots:

  • y shift vs x shift for all tiles
  • z shift vs x shift for all tiles
  • n_matches vs error for all tiles

In each case, the markers in the plots are numbers. These numbers indicate the tile the shift was found for. The number will be blue if nb.register_debug.n_matches > nb.register_debug.n_matches_thresh and red otherwise.

Parameters:

Name Type Description Default
nb Notebook

Notebook containing at least the register page.

 required
c Optional[int]

If None, will give option to scroll with mouse to change channel. If specify c, will show just one channel with no scrolling.

 None
outlier bool

If True, will plot shifts from nb.register_debug.transform_outlier instead of nb.register.transform. In this case, only tiles for which nb.register_debug.failed == True are plotted for each round/channel.

 False
Source code in coppafish/plot/register/diagnostics.py 
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def __init__(self, nb: Notebook, c: Optional[int] = None, outlier: bool = False):
    """
    For all affine transforms to imaging rounds/channels from the reference round computed in the `register` section
    of the pipeline, this plots the values of the shifts, `n_matches` (number of neighbours found) and
    `error` (average distance between neighbours).

    For each round and channel (channel is changed by scrolling with the mouse), there will be 3 plots:

    * y shift vs x shift for all tiles
    * z shift vs x shift for all tiles
    * `n_matches` vs `error` for all tiles

    In each case, the markers in the plots are numbers.
    These numbers indicate the tile the shift was found for.
    The number will be blue if `nb.register_debug.n_matches > nb.register_debug.n_matches_thresh` and red otherwise.

    Args:
        nb: Notebook containing at least the `register` page.
        c: If None, will give option to scroll with mouse to change channel. If specify c, will show just
            one channel with no scrolling.
        outlier: If `True`, will plot shifts from `nb.register_debug.transform_outlier` instead of
            `nb.register.transform`. In this case, only tiles for which
            `nb.register_debug.failed == True` are plotted for each round/channel.
    """
    self.outlier = outlier
    self.nb = nb
    if c is None:
        if self.outlier:
            # only show channels for which there is an outlier shift
            self.channels = np.sort(np.unique(np.where(nb.register_debug.failed)[2]))
            if len(self.channels) == 0:
                raise ValueError(f"No outlier transforms were computed")
        else:
            self.channels = np.asarray(nb.basic_info.use_channels)
    else:
        self.channels = [c]
    self.n_channels = len(self.channels)
    self.c_ind = 0
    self.c = self.channels[self.c_ind]

    n_cols = len(nb.basic_info.use_rounds)
    if nb.basic_info.is_3d:
        n_rows = 3
    else:
        n_rows = 2
    self.fig, self.ax = plt.subplots(n_rows, n_cols, figsize=(15, 7))
    self.fig.subplots_adjust(hspace=0.4, bottom=0.08, left=0.06, right=0.97, top=0.9)
    self.shift_info = self.get_ax_lims(self.nb, self.channels, self.outlier)
    self.update()
    if self.n_channels > 1:
        self.fig.canvas.mpl_connect('scroll_event', self.z_scroll)
    plt.show()

ICP

view_icp

Function to plot results of iterative closest point to find affine transform between ref_round/ref_channel and round r, channel c for tile t. Useful for debugging the register section of the pipeline.

Parameters:

Name Type Description Default
nb Notebook

Notebook containing results of the experiment. Must contain find_spots page. If contains register_initial and/or register pages, then transform from these will be used.

 required
t int

tile interested in.

 required
r int

Want to find the transform between the reference round and this round.

 required
c int

Want to find the transform between the reference channel and this channel.

 required
Source code in coppafish/plot/register/icp.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
def view_icp(nb: Notebook, t: int, r: int, c: int):
    """
    Function to plot results of iterative closest point to find affine transform between
    `ref_round/ref_channel` and round `r`, channel `c` for tile `t`.
    Useful for debugging the `register` section of the pipeline.

    Args:
        nb: Notebook containing results of the experiment. Must contain `find_spots` page.
            If contains `register_initial` and/or `register` pages, then transform from these will be used.
        t: tile interested in.
        r: Want to find the transform between the reference round and this round.
        c: Want to find the transform between the reference channel and this channel.
    """
    config = nb.get_config()['register']
    if nb.basic_info.is_3d:
        neighb_dist_thresh = config['neighb_dist_thresh_3d']
    else:
        neighb_dist_thresh = config['neighb_dist_thresh_2d']
    z_scale = [1, 1, nb.basic_info.pixel_size_z / nb.basic_info.pixel_size_xy]
    point_clouds = []
    # 1st point cloud is imaging one as does not change
    point_clouds = point_clouds + [spot_yxz(nb.find_spots.spot_details, t, r, c)]
    # only keep isolated spots, those whose second neighbour is far away
    # Do this only for imaging point cloud as that is what is done in pipeline/register
    isolated = get_isolated_points(point_clouds[0] * z_scale, 2 * neighb_dist_thresh)
    point_clouds[0] = point_clouds[0][isolated]

    # 2nd is untransformed reference point cloud
    r_ref = nb.basic_info.ref_round
    c_ref = nb.basic_info.ref_channel
    point_clouds = point_clouds + [spot_yxz(nb.find_spots.spot_details, t, r_ref, c_ref)]
    z_scale = z_scale[2]

    # Add shifted reference point cloud
    if nb.has_page('register_initial'):
        shift = nb.register_initial.shift[t, r]
    else:
        shift = view_register_search(nb, t, r, return_shift=True)
    point_clouds = point_clouds + [point_clouds[1] + shift]

    # Add reference point cloud transformed by an affine transform
    transform_outlier = None
    if nb.has_page('register'):
        transform = nb.register.transform[t, r, c]
        if nb.has_page('register_debug'):
            # If particular tile/round/channel was found by regularised least squares
            transform_outlier = nb.register_debug.transform_outlier[t, r, c]
            if np.abs(transform_outlier).max() == 0:
                transform_outlier = None
    else:
        start_transform = np.eye(4, 3)  # no scaling just shift to start off icp
        start_transform[3] = shift * [1, 1, z_scale]
        transform = get_single_affine_transform(point_clouds[1], point_clouds[0], z_scale, z_scale,
                                                start_transform, neighb_dist_thresh, nb.basic_info.tile_centre,
                                                config['n_iter'])[0]

    if not nb.basic_info.is_3d:
        # use numpy not jax.numpy as reading in tiff is done in numpy.
        tile_sz = np.array([nb.basic_info.tile_sz, nb.basic_info.tile_sz, 1], dtype=np.int16)
    else:
        tile_sz = np.array([nb.basic_info.tile_sz, nb.basic_info.tile_sz, nb.basic_info.nz], dtype=np.int16)

    if transform_outlier is not None:
        point_clouds = point_clouds + [apply_transform(point_clouds[1], transform_outlier, nb.basic_info.tile_centre,
                                                       z_scale, tile_sz)[0]]

    point_clouds = point_clouds + [apply_transform(point_clouds[1], transform, nb.basic_info.tile_centre, z_scale,
                                                   tile_sz)[0]]
    pc_labels = [f'Imaging: r{r}, c{c}', f'Reference: r{r_ref}, c{c_ref}', f'Reference: r{r_ref}, c{c_ref} - Shift',
                 f'Reference: r{r_ref}, c{c_ref} - Affine']
    if transform_outlier is not None:
        pc_labels = pc_labels + [f'Reference: r{r_ref}, c{c_ref} - Regularized']
    view_point_clouds(point_clouds, pc_labels, neighb_dist_thresh, z_scale,
                      f'Transform of tile {t} to round {r}, channel {c}')
    plt.show()

view_icp_reg

Function to plot how regularisation changes the affine transform found through iterative closest point between ref_round/ref_channel and round \(r\), channel \(c\) for tile \(t\).

Useful for finding suitable values for config['register']['regularize_constant'] and config['register']['regularize_factor'].

Parameters:

Name Type Description Default
nb Notebook

Notebook containing results of the experiment. Must contain find_spots page. Must contain register_initial/register_debug pages if start_transform/reg_transform not specified.

 required
t int

tile interested in.

 required
r int

Want to find the transform between the reference round and this round.

 required
c int

Want to find the transform between the reference channel and this channel.

 required
reg_constant Optional[List]

int [n_reg] Constant used when doing regularized least squares. Will be a point cloud produced for affine transformed produced with each of these parameters. Value will be indicated by \(\lambda\) in legend/buttons. If None, will use config['register']['regularize_constant'].

 None
reg_factor Optional[List]

float [n_reg] Factor to boost rotation/scaling term in regularized least squares. Must be same length as reg_constant. Value will be indicated by \(\mu\) in legend/buttons. If None, will use config['register']['regularize_factor'].

The regularized term in the loss function for finding the transform is:

\(0.5\lambda (\mu D_{scale}^2 + D_{shift}^2)\)

Where:

  • \(D_{scale}^2\) is the squared distance between transform[:3, :] and reg_transform[:3, :]. I.e. the squared difference of the scaling/rotation part of the transform from the target.
  • \(D_{shift}^2\) is the squared distance between transform[3] and reg_transform[3]. I.e. the squared difference of the shift part of the transform from the target.
  • \(\lambda\) is reg_constant - the larger it is, the smaller \(D_{scale}^2\) and \(D_{shift}^2\).
  • \(\mu\) is reg_factor - the larger it is, the smaller \(D_{scale}^2\).
 None
reg_transform Optional[np.ndarray]

Transform to regularize towards i.e. the expected transform. If not specified, will use average transformation based on nb.register_debug.av_scaling[c] and nb.register_debug.av_shifts[t, r] with no rotation. This is the same as what is used in coppafish.register.icp.

 None
start_transform Optional[np.ndarray]

Initial transform to use as starting point for ICP. If None, will use nb.register_initial.shift[t, r] for the non-regularized case and reg_transform for the regularized case to match method used in coppafish.register.icp.

 None
plot_residual bool

If True, will run plot_reg_residual as well.

 False
Source code in coppafish/plot/register/icp.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def view_icp_reg(nb: Notebook, t: int, r: int, c: int, reg_constant: Optional[List] = None,
                 reg_factor: Optional[List] = None, reg_transform: Optional[np.ndarray] = None,
                 start_transform: Optional[np.ndarray] = None, plot_residual: bool = False):
    """
    Function to plot how regularisation changes the affine transform found through iterative closest point between
    `ref_round/ref_channel` and round $r$, channel $c$ for tile $t$.

    Useful for finding suitable values for `config['register']['regularize_constant']`
    and `config['register']['regularize_factor']`.


    Args:
        nb: Notebook containing results of the experiment. Must contain `find_spots` page.
            Must contain `register_initial`/`register_debug` pages if `start_transform`/`reg_transform` not specified.
        t: tile interested in.
        r: Want to find the transform between the reference round and this round.
        c: Want to find the transform between the reference channel and this channel.
        reg_constant: `int [n_reg]`
            Constant used when doing regularized least squares.
            Will be a point cloud produced for affine transformed produced with each of these parameters.
            Value will be indicated by $\lambda$ in legend/buttons.
            If `None`, will use `config['register']['regularize_constant']`.
        reg_factor: `float [n_reg]`
            Factor to boost rotation/scaling term in regularized least squares.
            Must be same length as `reg_constant`.
            Value will be indicated by $\mu$ in legend/buttons.
            If `None`, will use `config['register']['regularize_factor']`.

            The regularized term in the loss function for finding the transform is:

            $0.5\lambda (\mu D_{scale}^2 + D_{shift}^2)$

            Where:

            * $D_{scale}^2$ is the squared distance between `transform[:3, :]` and `reg_transform[:3, :]`.
            I.e. the squared difference of the scaling/rotation part of the transform from the target.
            * $D_{shift}^2$ is the squared distance between `transform[3]` and `reg_transform[3]`.
            I.e. the squared difference of the shift part of the transform from the target.
            * $\lambda$ is `reg_constant` - the larger it is, the smaller $D_{scale}^2$ and $D_{shift}^2$.
            * $\mu$ is `reg_factor` - the larger it is, the smaller $D_{scale}^2$.
        reg_transform: Transform to regularize towards i.e. the expected transform.
            If not specified, will use average transformation based on `nb.register_debug.av_scaling[c]`
            and `nb.register_debug.av_shifts[t, r]` with no rotation.
            This is the same as what is used in `coppafish.register.icp`.
        start_transform: Initial transform to use as starting point for ICP.
            If `None`, will use `nb.register_initial.shift[t, r]` for the non-regularized case
            and `reg_transform` for the regularized case to match method used in `coppafish.register.icp`.
        plot_residual: If `True`, will run plot_reg_residual as well.
    """
    config = nb.get_config()['register']
    n_iter = config['n_iter']
    if nb.basic_info.is_3d:
        neighb_dist_thresh = config['neighb_dist_thresh_3d']
    else:
        neighb_dist_thresh = config['neighb_dist_thresh_2d']
    z_scale = [1, 1, nb.basic_info.pixel_size_z / nb.basic_info.pixel_size_xy]
    if not nb.basic_info.is_3d:
        # use numpy not jax.numpy as reading in tiff is done in numpy.
        tile_sz = np.array([nb.basic_info.tile_sz, nb.basic_info.tile_sz, 1], dtype=np.int16)
    else:
        tile_sz = np.array([nb.basic_info.tile_sz, nb.basic_info.tile_sz, nb.basic_info.nz], dtype=np.int16)
    point_clouds = []
    # 1st point cloud is imaging one as does not change
    point_clouds = point_clouds + [spot_yxz(nb.find_spots.spot_details, t, r, c)]
    # only keep isolated spots, those whose second neighbour is far away
    # Do this only for imaging point cloud as that is what is done in pipeline/register
    isolated = get_isolated_points(point_clouds[0] * z_scale, 2 * neighb_dist_thresh)
    point_clouds[0] = point_clouds[0][isolated]

    # 2nd is untransformed reference point cloud
    r_ref = nb.basic_info.ref_round
    c_ref = nb.basic_info.ref_channel
    point_clouds = point_clouds + [spot_yxz(nb.find_spots.spot_details, t, r_ref, c_ref)]
    z_scale = z_scale[2]

    # 3rd is ref point cloud transformed according to affine transform with no regularisation
    if start_transform is None:
        # no scaling just shift to start off icp as used in pipeline if no start_transform given
        shift = nb.register_initial.shift[t, r]
        start_transform = np.eye(4, 3)
        start_transform[3] = shift * [1, 1, z_scale]
    transform_no_reg = get_single_affine_transform(point_clouds[1], point_clouds[0], z_scale, z_scale,
                                                   start_transform, neighb_dist_thresh, nb.basic_info.tile_centre,
                                                   n_iter)[0]
    point_clouds = point_clouds + [apply_transform(point_clouds[1], transform_no_reg, nb.basic_info.tile_centre, z_scale,
                                                   tile_sz)[0]]

    # 4th is ref point cloud transformed according to target reg transform
    if reg_transform is None:
        # If reg transform not specified, use av scaling for c and av shift for t, r.
        # Also set start_transform to reg_transform if not specified
        # This is same as what is done in `coppafish.register.base.icp` for t/r/c where regularisation required.
        reg_transform = np.eye(4, 3) * nb.register_debug.av_scaling[c]
        reg_transform[3] = nb.register_debug.av_shifts[t, r]
        start_transform = reg_transform.copy()
    point_clouds = point_clouds + [apply_transform(point_clouds[1], reg_transform, nb.basic_info.tile_centre, z_scale,
                                                   tile_sz)[0]]
    pc_labels = [f'Imaging: r{r}, c{c}', f'Reference: r{r_ref}, c{c_ref}',
                 r'$\lambda=0$', r'$\lambda=\infty$']

    # Now add ref point cloud transformed according to reg transform found with all reg params given.
    # If None given, use default params.
    if reg_constant is None:
        reg_constant = [config['regularize_constant']]
    if reg_factor is None:
        reg_factor = [config['regularize_factor']]

    n_reg = len(reg_constant)
    if n_reg != len(reg_factor):
        raise ValueError(f"reg_constant and reg_factor need to be same size but they "
                         f"have {n_reg} and {len(reg_factor)} values respectively")
    transforms = np.zeros((n_reg, 4, 3))
    for i in range(n_reg):
        # Deviation in scale/rotation is much less than permitted deviation in shift so boost scale reg constant.
        reg_constant_scale = np.sqrt(0.5 * reg_constant[i] * reg_factor[i])
        reg_constant_shift = np.sqrt(0.5 * reg_constant[i])
        transforms[i] = \
            get_single_affine_transform(point_clouds[1], point_clouds[0], z_scale, z_scale, start_transform,
                                        neighb_dist_thresh, nb.basic_info.tile_centre, n_iter,
                                        reg_constant_scale, reg_constant_shift, reg_transform)[0]
        point_clouds += [apply_transform(point_clouds[1], transforms[i], nb.basic_info.tile_centre,
                                         z_scale, tile_sz)[0]]
        if reg_constant[i] > 1000:
            pc_labels += [r'$\lambda={:.0E},\mu={:.0E}$'.format(int(reg_constant[i]), int(reg_factor[i]))]
        else:
            pc_labels += [r'$\lambda={:.0f},\mu={:.0E}$'.format(int(reg_constant[i]), int(reg_factor[i]))]

    vpc = view_point_clouds(point_clouds, pc_labels, neighb_dist_thresh, z_scale,
                            f'Regularized transform of tile {t} to round {r}, channel {c}')
    n_matches_reg_target = np.sum(vpc.neighb[3] >= 0)

    if plot_residual and n_reg > 1:
        plot_reg_residual(reg_transform, transforms, reg_constant, reg_factor, transform_no_reg, n_matches_reg_target)
    else:
        plt.show()

plot_reg_residual

This shows how changing the regularization parameters affect how close the affine transform is to that which it was being regularized towards. E.g. it should show that the larger reg_constant, the smaller the difference (y-axis values in the plots).

There will be up to 4 plots, in each, the different colors refer to the different reg_constant/reg_factor combinations and the smaller the y-axis value, the closer the transform is to reg_transform. The different axis variables in the plot are explained in the reg_factor variable description below.

Parameters:

Name Type Description Default
reg_transform np.ndarray

float [4 x 3] Transform which was regularized towards i.e. the expected transform.

 required
transforms_plot List[np.ndarray]

[n_reg]. transforms_plot[i] is the [4 x 3] transform found with regularization parameters reg_constant[i] and reg_factor[i].

 required
reg_constant List

int [n_reg] Constant used when doing regularized least squares. Value will be indicated by \(\lambda\) on x-axis.

 required
reg_factor List

float [n_reg] Factor to boost rotation/scaling term in regularized least squares. Must be same length as reg_constant. Value will be indicated by \(\mu\) on x-axis.

The regularized term in the loss function for finding the transform is:

\(0.5\lambda (\mu D_{scale}^2 + D_{shift}^2)\)

Where:

  • \(D_{scale}^2\) is the squared distance between transform[:3, :] and reg_transform[:3, :]. I.e. the squared difference of the scaling/rotation part of the transform from the target.
  • \(D_{shift}^2\) is the squared distance between transform[3] and reg_transform[3]. I.e. the squared difference of the shift part of the transform from the target.
  • \(\lambda\) is reg_constant - the larger it is, the smaller \(D_{scale}^2\) and \(D_{shift}^2\).
  • \(\mu\) is reg_factor - the larger it is, the smaller \(D_{scale}^2\).
 required
transform_no_reg Optional[np.ndarray]

float [4 x 3]. Transform found with no regularization. If given, will show white line labelled by \(\lambda=0\) with y-axis value indicating value with no regularization.

 None
n_matches Optional[int]

Number of nearest neighbours found for reg_transform. If given, will show white line labelled by \(n_{matches}\) with x-axis value equal to this on the plots where \(\lambda\) is the x-axis variable. This is because we expect the regularization should have more of an effect if reg_constant > n_matches, i.e. the y-axis variable should go towards zero at x-axis values beyond this line.

 None
Source code in coppafish/plot/register/icp.py 
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
def plot_reg_residual(reg_transform: np.ndarray, transforms_plot: List[np.ndarray],
                      reg_constant: List, reg_factor: List, transform_no_reg: Optional[np.ndarray] = None,
                      n_matches: Optional[int] = None):
    """
    This shows how changing the regularization parameters affect how close the affine transform
    is to that which it was being regularized towards. E.g. it should show that
    the larger `reg_constant`, the smaller the difference (y-axis values in the plots).

    There will be up to 4 plots, in each, the different colors refer to the different
    `reg_constant`/`reg_factor` combinations and the smaller the y-axis value,
    the closer the transform is to `reg_transform`. The different axis variables in the plot are
    explained in the `reg_factor` variable description below.

    Args:
        reg_transform: `float [4 x 3]`
            Transform which was regularized towards i.e. the expected transform.
        transforms_plot: `[n_reg]`.
            `transforms_plot[i]` is the `[4 x 3]` transform found with regularization parameters
            `reg_constant[i]` and `reg_factor[i]`.
        reg_constant: `int [n_reg]`
            Constant used when doing regularized least squares.
            Value will be indicated by $\lambda$ on x-axis.
        reg_factor: `float [n_reg]`
            Factor to boost rotation/scaling term in regularized least squares.
            Must be same length as `reg_constant`.
            Value will be indicated by $\mu$ on x-axis.

            The regularized term in the loss function for finding the transform is:

            $0.5\lambda (\mu D_{scale}^2 + D_{shift}^2)$

            Where:

            * $D_{scale}^2$ is the squared distance between `transform[:3, :]` and `reg_transform[:3, :]`.
            I.e. the squared difference of the scaling/rotation part of the transform from the target.
            * $D_{shift}^2$ is the squared distance between `transform[3]` and `reg_transform[3]`.
            I.e. the squared difference of the shift part of the transform from the target.
            * $\lambda$ is `reg_constant` - the larger it is, the smaller $D_{scale}^2$ and $D_{shift}^2$.
            * $\mu$ is `reg_factor` - the larger it is, the smaller $D_{scale}^2$.
        transform_no_reg: `float [4 x 3]`.
            Transform found with no regularization. If given, will show white line
            labelled by $\lambda=0$ with y-axis value indicating value with no regularization.
        n_matches: Number of nearest neighbours found for `reg_transform`.
            If given, will show white line labelled by $n_{matches}$ with x-axis value
            equal to this on the plots where $\lambda$ is the x-axis variable.
            This is because we expect the regularization should have more of an effect if `reg_constant > n_matches`,
            i.e. the y-axis variable should go towards zero at x-axis values beyond this line.
    """
    n_reg = len(reg_constant)
    col_info = {}
    col_ind = 0
    if len(np.unique(reg_constant)) > 1:
        col_info[col_ind] = {}
        if len(np.unique(reg_factor)) == 1:
            col_info[col_ind]['title'] = r"Varying $\lambda$ ($\mu = {:.0E}$)".format(int(reg_factor[0]))
        else:
            col_info[col_ind]['title'] = f"Varying $\lambda$"
        col_info[col_ind]['x_label'] = r"$\lambda$"
        col_info[col_ind]['x'] = reg_constant
        col_info[col_ind]['x_lims'] = [int(0.9 * np.min(reg_constant)), int(1.1 * np.max(reg_constant))]
        if n_matches is not None:
            col_info[col_ind]['x_lims'] = [int(0.9 * np.min(list(reg_constant) + [n_matches])),
                                           int(1.1 * np.max(list(reg_constant) + [n_matches]))]
        col_info[col_ind]['log_x'] = np.ptp(col_info[col_ind]['x_lims']) > 100
        col_info[col_ind]['n_matches'] = n_matches
        col_ind += 1
    if len(np.unique(reg_factor)) > 1:
        col_info[col_ind] = {}
        if len(np.unique(reg_constant)) == 1:
            col_info[col_ind]['title'] = rf"Varying $\mu$ ($\lambda = {int(reg_constant[0])}$)"
        else:
            col_info[col_ind]['title'] = f"Varying $\mu$"
        col_info[col_ind]['x_label'] = r"$\mu$"
        col_info[col_ind]['x'] = reg_factor
        col_info[col_ind]['x_lims'] = [int(0.9 * np.min(reg_factor)), int(1.1 * np.max(reg_factor))]
        col_info[col_ind]['log_x'] = np.ptp(col_info[col_ind]['x_lims']) > 100
        col_info[col_ind]['n_matches'] = None

    if n_reg != len(reg_factor):
        raise ValueError(f"reg_constant and reg_factor need to be same size but they "
                         f"have {n_reg} and {len(reg_factor)} values respectively")
    if n_reg != len(transforms_plot):
        raise ValueError(f"reg_constant and transforms_plot need to be same size but they "
                         f"have {n_reg} and {len(transforms_plot)} values respectively")
    if len(col_info) == 0:
        raise ValueError("Not enough data to plot")

    y0 = [np.linalg.norm(transforms_plot[i][:3, :]-reg_transform[:3, :]) for i in range(n_reg)]
    y1 = [np.linalg.norm(transforms_plot[i][3] - reg_transform[3]) for i in range(n_reg)]
    if transform_no_reg is not None:
        y0_no_reg = np.linalg.norm(transform_no_reg[:3, :]-reg_transform[:3,:])
        y1_no_reg = np.linalg.norm(transform_no_reg[3] - reg_transform[3])
        y0_lim = [np.min(y0 + [y0_no_reg]) - 0.0002, np.max(y0 + [y0_no_reg]) + 0.0002]
        y1_lim = [np.min(y1 + [y1_no_reg]) - 0.2, np.max(y1 + [y1_no_reg]) + 0.2]
    else:
        y0_lim = [np.min(y0) - 0.0002, np.max(y0) + 0.0002]
        y1_lim = [np.min(y1) - 0.2, np.max(y1) + 0.2]
    y0_lim = np.clip(y0_lim, 0, np.inf)
    y1_lim = np.clip(y1_lim, 0, np.inf)
    colors = distinctipy.get_colors(n_reg, [(0, 0, 0), (1, 1, 1)])
    fig, ax = plt.subplots(2, len(col_info), figsize=(15, 7))
    if len(col_info) == 1:
        ax = ax[:, np.newaxis]
    for i in range(len(col_info)):
        if transform_no_reg is not None:
            ax[0, i].plot(col_info[i]['x_lims'], [y0_no_reg, y0_no_reg], 'w:')
            ax[0, i].text(np.mean(col_info[i]['x_lims']), y0_no_reg, r"$\lambda = 0$",
                          va='bottom', ha="center", color='w')
            ax[1, i].plot(col_info[i]['x_lims'], [y1_no_reg, y1_no_reg], 'w:')
            ax[1, i].text(np.mean(col_info[i]['x_lims']), y1_no_reg, r"$\lambda = 0$",
                          va='bottom', ha="center", color='w')
        if col_info[i]['n_matches'] is not None:
            ax[0, i].plot([n_matches, n_matches], y0_lim, 'w:')
            ax[0, i].text(n_matches, np.percentile(y0_lim, 10), r"$n_{matches}$",
                          va='bottom', ha="right", color='w', rotation=90)
            ax[1, i].plot([n_matches, n_matches], y1_lim, 'w:')
            ax[1, i].text(n_matches, np.percentile(y1_lim, 10), r"$n_{matches}$",
                          va='bottom', ha="right", color='w', rotation=90)
        ax[0, i].scatter(col_info[i]['x'], y0, color=colors)
        ax[1, i].scatter(col_info[i]['x'], y1, color=colors)
        ax[0, i].get_shared_x_axes().join(ax[0, i], ax[1, i])
        ax[1, i].set_xlabel(col_info[i]['x_label'])
        ax[0, i].set_title(col_info[i]['title'])
        if col_info[i]['log_x']:
            ax[1, i].set_xscale('log')
        ax[1, i].set_xlim(col_info[i]['x_lims'])
        if i == 1:
            ax[0, 0].get_shared_y_axes().join(ax[0, 0], ax[0, 1])
            ax[1, 0].get_shared_y_axes().join(ax[1, 0], ax[1, 1])
        ax[0, 0].set_ylim(y0_lim)
        ax[0, 0].set_ylabel("$D_{scale}$")
        ax[1, 0].set_ylim(y1_lim)
        ax[1, 0].set_ylabel("$D_{shift}$")
    fig.suptitle("How varying regularization parameters affects how similar transform found is to target transform")
    plt.show()