Change coordinate system

An useful visualization tool in bmiptools is ChangeCoordinateSystem, which allows to change the coordinate system on which the stack is visualized. Its structure is similar to a plugin, but it has no transformation dictionary. After the initialization, this tool takes a given stack \(S(k,j,i)\) evaluated on a 3d grid with coordinates \((k,j,i)\) and return another stack \(S'(r,s,t)\) evaluated on a different 3d grid with coordinates \((r,s,t)\), related to the previous coordinate by a change of reference frame.

Attention

This tool can be used via Python API only.

Preliminary facts

Given a stack \(S\) the value assumed in the voxel \((k,j,i)\) has to be understood as the value of some physical quantity describing the sample in the point \((k\Delta Z,j\Delta Y,i \Delta X)\) of the 3d space, where \((\Delta Z,\Delta Y,\Delta X)\) are the voxel size. All the points \((k,j,i)\) associated to the various voxels of the stack form a 3d grid, which is typically cartesian, i.e. expressed in a cartesian reference frame. For visualization purpose one may want to see the stack using a different coordinate system.

Let \((z,y,x)\) and :(s,t,u): be two coordinate system related by the invertible function \(f\) such that

\[\begin{split}(s,t,u) &= f(z,y,x) \\ (z,y,x) &= f^{-1}(s,t,u).\end{split}\]

Given a \(K \times J \times I\) stack \(S(k,j,i)\) evaluated in the \((z,y,x)\) coordinates with voxel size given by \((\Delta z, \Delta y, \Delta x)\), it can be evaluated on the \((s,t,u)\) coordinate system by following this procedure.

Define the new grid in the \((s,t,u)\) coordinate system. Practically this means to specify two of the following 3 quantities: number of points for each coordinate, coordinate range, and step size for each coordinate. By specifying two of them, the third follow. Indeed:
- specifying for the coordinate \(s\) the number of points \(N_s\) and its range \([S_{min},S_{max}]\), then the step step is given by \(\Delta s = (S_{max}-S_{min})/N_s\). Similar considerations hold for the other two coordinates.
- specifying for the coordinate \(s\) the number of points \(N_s\) and its step size \(\Delta s\), then the coordinate range is \([0,N_s \Delta s]\) which can be translated by specifying an offset value. Similar considerations hold for the other two coordinates.
- specifying for the coordinate \(s\) its range \([S_{min},S_{max}]\) and its step size \(\Delta s\), then the number of points is given by \(N_s = (S_{max}-S_{min})/\Delta s\). Similar considerations hold for the other two coordinates.
The new grid will contain \(N_s \times N_t \times N_u\) points \((l,m,n)\) corresponding to the point of the 3d space having \((s,t,u)\)-coordinates \((l\Delta s, m\Delta t, u\Delta u)\).
Interpolate the stack \(S\) on the grid defined in the \((z,y,x)\) coordinate system, obtaining the function \(\tilde{S}(k,j,i)\).
Evaluate the interpolated stack on the new grid, by using the inverse function \(f^{-1}\). More precisely, given a point in the new grid \((l,m,n)\) one computes the corresponding point in the 3d space \((s_l,t_m,u_n) = (l\Delta s, m\Delta t, n\Delta u)\). Then one compute the corresponding point in the grid expressed in the old coordinate system \((z,y,x)\), this can be done first by computing the point

\[(z_{(l,m,n)},y_{(l,m,n)},x_{(l,m,n)}) = f^{-1}(l\Delta s, m\Delta t, n\Delta u),\]

and later diving each coordinate by the corresponding step size, obtaining

\[(k_{(l,m,n)},j_{(l,m,n)},i_{(l,m,n)}) = \left(\frac{z_{(l,m,n)}}{\Delta z},\frac{y_{(l,m,n)}}{\Delta y},\frac{x_{(l,m,n)}}{\Delta x}\right).\]

At this point the stack in the expressed in the new coordinate system \(S'\), is defined as

\[S'(l,m,n) = \tilde{S}(k_{(l,m,n)},j_{(l,m,n)},i_{(l,m,n)}).\]

Spherical and cylindrical may be particularly useful and will be briefly described below.

Spherical coordinates

The function \(f\) mapping the cartesian coordinates \((z,y,x)\) into spherical coordinates \((r,\theta,\phi)\) is defined as follow

\[\begin{split}\left(\begin{array}{c} r \\ \theta \\ \phi \\ \end{array}\right) = \left(\begin{array}{c} \sqrt{x^2+y^2+z^2} \\ \mbox{arccos} \left(\frac{z}{r}\right) \\ \mbox{arctan2} \left(y,x\right) \\ \end{array}\right)\end{split}\]

where \(\mbox{arctan2}\) is the “2-argument arctangent”, and the inverse transformation is

\[\begin{split}\left(\begin{array}{c} z \\ y \\ x \\ \end{array}\right) = \left(\begin{array}{c} r\cos(\theta) \\ r\sin(\theta)\sin(\phi) \\ r\sin(\theta)\cos(\phi) \\ \end{array}\right)\end{split}\]

In bmiptools \(f\) and \(f^{-1}\) corresponds to the function cartesian_to_spherical and spherical_to_cartesian.

Cylindrical coordinates

The function \(f\) mapping the cartesian coordinates \((z,y,x)\) into spherical coordinates \((r,\theta,z)\) is defined as follow

\[\begin{split}\left(\begin{array}{c} r \\ \theta \\ z \\ \end{array}\right) = \left(\begin{array}{c} \sqrt{x^2+y^2} \\ \mbox{arctan2} \left(y,x\right) \\ z \\ \end{array}\right)\end{split}\]

where \(\mbox{arctan2}\) is the “2-argument arctangent”, and the inverse transformation is

\[\begin{split}\left(\begin{array}{c} z \\ y \\ x \\ \end{array}\right) = \left(\begin{array}{c} z \\ r\sin(\theta) \\ r\cos(\theta) \\ \end{array}\right)\end{split}\]

In bmiptools \(f\) and \(f^{-1}\) corresponds to the function cartesian_to_cylindrical and cylindrical_to_cartesian.

Tool usage

To use this tool the first thing to is to initialize it. During the initialization the following information need to be declared:

reference_frame_origin: it is the position of the origin in the new reference frame, and have to be specified as a list, or tuple, or a numpy array.
xyz_to_XYZ_inv_map: this is a python function with 3 input representing the inverse mapping between the reference frame ‘xyz’ and the new reference frame ‘XYZ’ (i.e. is the function \(f^{-1}\) above).
xyz_to_XYZ_specs: this is a dictionary used to define the new reference frame, which need to have the following structure:
- new_shape, a tuple specifying the shape of the stack in the new reference frame, i.e. the number of points for each new coordinate of the new grid on which the stack will be define after the coordinate change.
- X_bounds, a list with two arguments: the minimum and maximum value of the new X coordinate.
- Y_bounds, a list with two arguments: the minimum and maximum value of the new Y coordinate.
- Z_bounds, a list with two arguments: the minimum and maximum value of the new Z coordinate.
- XYZ_ordering, is used to specify the ordering of the new coordinates. If None the order is the ‘ZYX’ otherwise on can specify a list with the new ordering, e.g. for ‘XYZ’ use [2,1,0], for ‘YXZ’ use [1,2,0], while for ‘XZY’ use [2,0,1].
use_xyz_ordering: if True, the default ordering of the axis in the stack (which is ‘zyx’) is converted to the cartesian ordering (i.e. ‘xyz’).

Therefore to initialize this tool the code below can be used, which assume the initial coordinate system as cartesian, the new coordinate system as spherical.

import numpy as np
from bmiptools.visualization.geometric.change_coordinate_system import ChangeCoordinateSystem
from bmiptools.core.math_utils import spherical_to_cartesian


ccs = ChangeCoordinateSystem(reference_frame_origin = [100,100,100],
                             xyz_to_XYZ_inv_map = spherical_to_cartesian,
                             xyz_to_XYZ_specs = {'new_shape': (100,180,360),
                                                 'X_bounds': [0,100],
                                                 'Y_bounds': [0,np.pi],
                                                 'Z_bounds': [0,2*np.pi]},
                             use_xyz_ordering = True)

Once the initialization is done, the application to this tool to a stack stack can be done as for the plugin with the method .transform.

from bmiptools.visualization.graphic_tools.basic_graphic_tools import Basic2d as b2d

# load/ create a stack.
# stack = ...

# See slice 0 before the change of coordinate system.
# b2d.show_image(stack[0])

ccs.transform(stack)

# See slice 0 after the change of coordinate system.
b2d.show_image(stack[0])