GitLab

Biomedical imaging is broadly concerned with the identification and characterization of some biological state, and/or determining the relationships between some intervention and the future evolution of that biological state. In general, the benefits of imaging come from the ability to derive inference based on the spatial relationships present within the image. Volumetric imaging has been a steadily progressing application of imaging, and this is particularly evident in such modalities as MRI and structured light sheet microscopy. Modern imaging frequently includes a temporal component, forming videos, or a 'channel' or 'color' component, which capture a particular subset of the total information available from the imaging subject. One can think of, for example, the use of multiple fluorescent filters while performing live-cell microscopy, each of which captures a particular absorption/emission band combination; or one can think of various T_1 and T_2 weighted images, as well as derived images such as the ADC or BOLD/TOLD.

Broadly, then, we can describe an image as some object I(x,y,z,c,t) which has spatial extent in the \mathcal{X}, \mathcal{Y}, and \mathcal{Z} dimensions and has a defined value for each channel \mathcal{C} at each point in space, and which evolves through time \mathcal{T}.

If we fix one or more of these variables , then we describe the resultant image as being parametrized by the particular values of those variables. For example, I(x,y,z;c,t) denotes a volumetric image at a particular point in time for a particular channel, and I(x,y,c,t;z) describes a multi-channel video taken as a slice through a particular z-location. Occasionally we will suppress the parameter indices when the context makes the notation unambiguous, i.e. I(x,y,z;t,c) \equiv I(x,y,z)

Frequently it will happen that multiple channels do not share the same underlying basis, such as when observing raw vs. derived images such as T_1 maps vs. ADC maps. In these cases, we take the direct sum of those bases, i.e. C = C_{T_1} \oplus C_{ADC}. This allows us to treat arbitrarily many distinct channels as subspaces of the overall C dimension. This treatment is reasonable, as a single sensor or transducer is generally unable to detect in multiple channels at once.

Images satisfy all the properties required in order for them to be treated as vectors in a vector space, and as they are finite dimensional (though the dimension may be very large), we may represent them as scalar arrays. So, for example, a volumetric image I(x,y,z;t,c) may be thought of as a three-dimensional parallelipiped, whereas a particular slice through that volume I(x,z;y,t,c) is a two-dimensional slab. Treating the image as a tensor, we may perform a reshaping operation to make the handling of the image more convenient. This is technically due to taking the tensor product of the multiple dimensions involved in the reshaping. For example, an image I(x,y,t,c) may be flattened to a two-dimensional image I(\lbrace{x,c}\rbrace,\lbrace{y,t}\rbrace), where the brackets indicate the tensor product of the enclosed dimensional variables, i.e. \lbrace a,b \rbrace \equiv a \otimes b. This may be repeated for each dimension or combination of dimensions, so that we can convert a 4-D tensor of size (N_x,N_y,N_t,N_c) into a 1-D vector of size ({N_x}\times{N_y}\times{N_t}\times{N_c},1).

Treating tensors as vectors enables us to correspondingly represent the linear transformations between tensors as matrices. Thus we can perform a linear transformation of an image I_1(x_1,y_1,z_1,c_1,t_1) \overset{H}{\rightarrow} I_2(x_2,y_2,z_2,c_2,t_2) simply by performing a matrix multiplication \operatorname{vec}(I_2) = H\operatorname{vec}(I_1). In general the two spaces the images belong to do not need to be the same, but the overall procedure may be treated as an isometry of a common linear space I({x_1}\oplus{x_2},{y_1}\oplus{y_2},{z_1}\oplus{z_2},{c_1}\oplus{c_2},{t_1}\oplus{t_2}), of which each image space is a subspace, i.e. I_1,I_2 \subset I. By casting the imaging problem in this way, the transformation between the two images becomes a submatrix of a larger, square matrix with dimension equal to the sum of the dimensions of the input and output spaces.

Interpolation is an example of this, where the data points themselves and the interpolated points are joined into one overall space, and the interpolation itself is then a map from one to the other. Similarly, the process of imaging an object may be represented as a transformation from the 'real' object space \mathcal{O} to the measurement space \mathcal{M}. If we represent the original transformation as

\tag{1} I_2 = HI_1

then we write the corresponding operator in the large space as

\tag{2}
\begin{bmatrix}
   I_2 \\
   \bold{0} 
\end{bmatrix}
=
\begin{bmatrix}
   \mathbb{I}_{I_2} & H \\
   \bold{0} & \bold{0}
\end{bmatrix}
\begin{bmatrix}
   \bold{0} \\
   I_1 
\end{bmatrix}

which we write in condensed form as

\tag{3} \tilde{I_2} = \tilde{H}\tilde{I_1}

where \mathbb{I}_{I_2} denotes the identity matrix of dimension equal to the space \mathcal{I_2} : If I_2 has 3 coordinates, then \mathbb{I}_{I_2} is the 3x3 identity matrix.

\tilde{H} enjoys some desirable properties, not least of which is that it is a projection matrix, that is \tilde{H}^2 = \tilde{H}. The transformation between spaces has been lifted to become a subspace projection of a larger space.

If we invert (1) so that

\tag{4} H^+I_2 = H^+HI_1 = I_1

the corresponding lifted transformation is

\tag{5}
\begin{bmatrix}
   \bold{0} \\
   I_1 
\end{bmatrix}
=
\begin{bmatrix}
   \bold{0} & \bold{0} \\
   H^+ & \mathbb{I}_{I_1}
\end{bmatrix}
\begin{bmatrix}
   I_2 \\
   \bold{0} 
\end{bmatrix}

We can also see that we can define a so-called 'completion' operator by

\tag{6}
\begin{bmatrix}
   I_2 \\
   I_1 
\end{bmatrix}
=
\begin{bmatrix}
   \mathbb{I}_{I_2} & \bold{0} \\
   H^+ & \mathbb{I}_{I_1}
\end{bmatrix}
\begin{bmatrix}
   I_2 \\
   \bold{0} 
\end{bmatrix}

Practical examples

In MSOT, we acquire pressure transducer data at a certain time, wavelength, and axial position. This transducer data we denote p(\vec{r_s},\vec{t_s}), where \vec{r_s} denotes the spatial location of the detecting transducer and \vec{t_s} denotes the sample times at which pressure data is acquired. Indexing N_s transducers by j \subset \lbrace\mathcal{J} = 1,2,...,N_s\rbrace, we have \overrightarrow{r_{s=j}} \equiv (x_{s=j},y_{s=j},z_{s=j}) so that each transducer has a defined location. Similarly, we know the sample times \vec{t_s}.

Our end goal is to determine what the corresponding multispectral image is. Restricting ourselves to two-dimensional spatial images for the time being, we seek to solve the following equation:

p(\vec{r_s},\vec{t_s})
=
H(\vec{r_s},\vec{t_s},\bar{x},\bar{y},\bar{z})
I(\bar{x},\bar{y},\bar{z})

or, in matrix-vector form,

\vec{p} 
\equiv
{\vec{p}(\vec{r_s},\vec{t_s})
=
H(\vec{r_s},\vec{t_s},\bar{x},\bar{y})
\vec{I}(\bar{x},\bar{y})} 
\equiv
H\vec{I}