How can we register three-dimensional points in space acquired from the same structure in different time-points? This is a problem I faced in my Ph.D. project where I was dealing with images of the same structure acquired at different time-points. For that project I was imaging dendritic spines, microscopic protrusions emanating from neurons, in alive mice (Figure 1).

Every spine was assigned an ID and tracked over all available time-points, and we were interested in spine positions on the dendrite.

Even though I made every effort to acquire the images under similar conditions their alignment was never going to be perfect from one time-point to the next. To overcome the problems that such misalignments may cause when considering changes in distance to a fixed point over time, say a cortical column’s center or the surface of the brain, I had to perform some sort of image registration.

Two sources of information that can be useful when registering images are the intensity profiles of the image and the geometry of some of the image’s features. For example, we could choose to calculate some correlation metric based on the intensity values of the image as a whole. Alternatively, we could take into account the relative positions of some structures in the image in order to generate the geometrical transformation. In our case, since spine morphology can change a lot from session to session, especially since the animal received externally-generated stimulation, I chose to focus on the second approach.

A very simple method to register images based on feature coordinates is to apply a linear transformation, which assumes that no warping can occur between imaging sessions, i.e. the transformation is uniform across all structures considered. The possible transformations under this assumption are translation, rotation, scaling, and possibly reflection. The solution to finding the optimal linear transformation based on the sum of squared errors criterion is called procrustes analysis. Below I provide an example of procrustes analysis implemented in MATLAB for the dendrite shown in Figure 1.

_{
(Please note: this post contains interactive visualizations through
plot.ly which are best viewed in landscape mode on small screens.)
}

## Implementation

First, let’s have a look at the data for this example:

```
data = readtable('data.csv'); % Load data
data = sortrows(data,'correlationID','ascend');
% Sort data based on correlation ID
disp(data(data.correlationID == 1, :))
% Select data for an example spine
```

```
session correlationID x0 x1 y0 y1 z
_______ _____________ ______ ______ ______ ______ __
1 1 7.1997 4.7421 33.728 33.421 13
2 1 8.0108 5.6118 32.435 32.26 13
3 1 8.1279 5.2061 32.497 32.387 13
4 1 9.4179 6.6599 34.94 35.295 17
5 1 11.554 8.9876 35.097 34.933 14
```

The above code loads the annotation data for the dendrite’s spines
and displays the data for a single spine.
The variable `session`

refers to the imaging day,
`correlationID`

to the spine identity,
`x0`

and `y0`

to the coordinates of the base of the spine,
`x1`

and `y1`

to the coordinates of the tip of the spine,
and `z`

to the slice corresponding to the plane
in which the spine had the brightest fluorescence.
All coordinates are in micrometers.

Let’s then plot the annotations for all spines from the first session.
To do that I have defined a function, `selectRows`

, that will simply
return all coordinates that match a logical index,
as I will be doing that a few times.

```
type('selectRows.m')
```

```
function [X, Y, Z] = selectRows(data, ind)
X = data.x0(ind);
Y = data.y0(ind);
Z = data.z(ind);
end
```

```
imatlab_export_fig('fig2plotly')
% Use the plotly graphics engine
% in a Jupyter notebook with the imatlab kernel
indSession1 = data.session == 1;
[X1, Y1, Z1] = selectRows(data, indSession1);
% Select coordinates from the first session for illustration
figure()
plot3(X1, Y1, Z1, 'ko', 'MarkerFaceColor', 'k')
```

We want to see what kind of displacements we are facing so I am going to add spines from the second session to the previous plot:

```
indSession2 = data.session == 2;
[X2, Y2, Z2] = selectRows(data, indSession2);
figure(); hold 'on';
plot3(X1, Y1, Z1, 'ko', 'MarkerFaceColor', 'k')
plot3(X2, Y2, Z2, 'ro', 'MarkerFaceColor', 'r')
legend('Session 1', 'Session 2')
```

We see that spines appear and disappear, and that part of the dendrite is shifted. Next, we will select those spines that are present in both the first and second sessions:

```
commonIDs = intersect(data.correlationID(indSession1), ...
data.correlationID(indSession2));
idInd = ismember(data.correlationID, commonIDs);
[sX1, sY1, sZ1] = selectRows(data, indSession1 & idInd);
[sX2, sY2, sZ2] = selectRows(data, indSession2 & idInd);
```

We are now ready to apply the procrustes analysis:

```
[d, XYZ, C] = procrustes([sX1, sY1, sZ1], [sX2, sY2, sZ2]);
```

The sum-of-squared-errors criterion is defined as:

$$ d = \sum _ {i=1 }^{n} \sum _ {j=1}^{p} (M _ {i,j} - N _ {i,j}) ^ 2 \enspace \text, \tag{1} $$

where $n$ the number of spines used for the registration,
$p$ the number of coordinates per spine $\{ x, y, z \}$, i.e. $p = 3$,
$M$ the coordinates of spines in the first session,
and $N$ the new set of coordinates for the second session after transformation. The function `procrustes`

reports $d’$, stored in the above call in variable `d`

, which is $d$ standarized by the sum of square errors of a centered version of $M$:

$$ d’ = \frac {d} {\sum \sum\limits _ {i = 1}^{n} {(M _ {i} - \overline{M})} ^ 2} \enspace \text, \tag{2} $$

where $\overline{M}$ the average $\{x, y, z \}$ coordinates across all spines of the first session.

In our example $d’$ was pretty low:

```
disp(d)
```

```
0.0145
```

which indicates that the registration was succesful.

The variable `XYZ`

contains the transformed coordinates for the included spines of the second session.
We are not going to be using it since we are interested in applying the transformation to all spines of the second session, not just those that were also present in the first session.

The variable `C`

contains the components for the tranformation:

```
disp(C)
```

```
T: [3×3 double]
b: 0.9974
c: [42×3 double]
```

`T`

corresponds to the rotation and reflection component, `b`

to the scale component, and `c`

to the translation component. They are chosen so as to minimize $d’$ while satisfying the following equation:
$$ N = b*M’*T + c \tag{3}
$$

By simply applying this transformation to all spines of the second session we obtain the new coordinates:

```
N = C.b * [X2, Y2, Z2] * C.T + C.c(1,:);
```

Let’s now plot again the coordinates for spines from the first session and our new coordinates for spines from the second session, and see how we did:

```
figure(); hold on;
plot3(X1, Y1, Z1, 'ko', 'MarkerFaceColor', 'k')
plot3(X2, Y2, Z2, 'ro', 'MarkerFaceColor', 'r')
plot3(N(:,1), N(:,2), N(:,3), 'bo', 'MarkerFaceColor', 'b')
legend('Original Session 1', 'Original Session 2', 'Transformed Session 2')
```

You may want to use your mouse or finger to navigate through the image and click on the legend of the figure above to toggle viewing the corresponding data. We now simply need to write a function that runs all the above for all sessions.

```
type('applyProcrustes.m')
```

```
function [N] = applyProcrustes(data)
nSession = max(data.session);
sessVector = num2cell(2:nSession);
N = cellfun(@applyProcrustesSingleSession, ...
repmat({data}, length(sessVector), 1), sessVector', ...
'UniformOutput', 0);
function [N] = applyProcrustesSingleSession(data, iS)
indSession1 = data.session == 1;
indSession2 = data.session == iS;
[X2, Y2, Z2] = selectRows(data, indSession2);
commonIDs = intersect(data.correlationID(indSession1), ...
data.correlationID(indSession2));
idInd = ismember(data.correlationID, commonIDs);
[sX1, sY1, sZ1] = selectRows(data, indSession1 & idInd);
[sX2, sY2, sZ2] = selectRows(data, indSession2 & idInd);
[~, ~, C] = procrustes([sX1, sY1, sZ1], [sX2, sY2, sZ2]);
N = C.b * [X2, Y2, Z2] * C.T + C.c(1,:);
end
end
```

All we have done is to wrap our previous commands in a sub-function
which accepts the session variable in a generalized way,
and pass that function to `cellfun`

that in turn applies it to data for all sessions.
The usage of `cellfun`

over a conventional for-loop is advantageous in this case,
as it provides an order of magnitude gain in speed,
something that can make a difference if the function is run over many datasets.

One thing that we could have done differently is, instead of arbitrarily choosing one session (the first in our case) as the reference, we could take into account information from all sessions. For example, in generalized procrustes analysis the reference is chosen as an optimally determined “mean” shape. Such functionality is already implemented in the R package shapes.

## Comments

No comments.

## Submit comment