.. _transformations: ========================== Transformation use cases ========================== Use cases for defining and using transforms on images. We should be very careful to only use the terms ``x, y, z`` to refer to physical space. For voxels, we should use ``i, j, k``, or ``i', j', k'`` (i prime, j prime k prime). I have an image *Img*. Image Orientation ----------------- I would like to know what the voxel sizes are. I would like to determine whether it was acquired axially, coronally or sagittally. What is the brain orientation in relation to the voxels? Has it been acquired at an oblique angle? What are the voxel dimensions?:: img = load_image(file) cm = img.coordmap print cm input_coords axis_i: axis_j: axis_k: effective pixel dimensions axis_i: 4mm axis_j: 2mm axis_k: 2mm input/output mapping x y z ------------ i| 90 90 0 j| 90 0 90 k| 180 90 90 input axis_i maps exactly to output axis_z input axis_j maps exactly to output axis_y input axis_k maps exactly to output axis_x flipped 180 output_coords axis0: Left -> Right axis1: Posterior -> Anterior axis2: Inferior -> Superior In the case of a mapping that does not exactly align the input and output axes, something like:: ... input/output mapping input axis0 maps closest to output axis2 input axis1 maps closest to output axis1 input axis2 maps closest to output axis0 ... If the best matching axis is reversed compared to input axis:: ... input axis0 maps [closest|exactly] to negative output axis2 and so on. Creating transformations / coordinate maps ------------------------------------------- I have an array *pixelarray* that represents voxels in an image and have a matrix/transform *mat* which represents the relation between the voxel coordinates and the coordinates in scanner space (world coordinates). I want to associate the array with the matrix:: img = load_image(infile) pixelarray = np.asarray(img) (*pixelarray* is an array and does not have a coordinate map.):: pixelarray.shape (40,256,256) So, now I have some arbitrary transformation matrix:: mat = np.zeros((4,4)) mat[0,2] = 2 # giving x mm scaling mat[1,1] = 2 # giving y mm scaling mat[2,0] = 4 # giving z mm scaling mat[3,3] = 1 # because it must be so # Note inverse diagonal for zyx->xyz coordinate flip I want to make an ``Image`` with these two:: coordmap = voxel2mm(pixelarray.shape, mat) img = Image(pixelarray, coordmap) The ``voxel2mm`` function allows separation of the image *array* from the size of the array, e.g.:: coordmap = voxel2mm((40,256,256), mat) We could have another way of constructing image which allows passing of *mat* directly:: img = Image(pixelarray, mat=mat) or:: img = Image.from_data_and_mat(pixelarray, mat) but there should be "only one (obvious) way to do it". Composing transforms '''''''''''''''''''' I have two images, *img1* and *img2*. Each image has a voxel-to-world transform associated with it. (The "world" for these two transforms could be similar or even identical in the case of an fmri series.) I would like to get from voxel coordinates in *img1* to voxel coordinates in *img2*, for resampling:: imgA = load_image(infile_A) vx2mmA = imgA.coordmap imgB = load_image(infile_B) vx2mmB = imgB.coordmap mm2vxB = vx2mmB.inverse # I want to first apply transform implied in # cmA, then the inverse of transform implied in # cmB. If these are matrices then this would be # np.dot(mm2vxB, vx2mmA) voxA_to_voxB = mm2vxB.composewith(vx2mmA) The (matrix) multiply version of this syntax would be:: voxA_to_voxB = mm2vxB * vx2mmA Composition should be of form ``Second.composewith(First)`` - as in ``voxA_to_voxB = mm2vxB.composewith(vx2mmA)`` above. The alternative is ``First.composewith(Second)``, as in ``voxA_to_voxB = vx2mmA.composewith(mm2vxB)``. We choose ``Second.composewith(First)`` on the basis that people need to understand the mathematics of function composition to some degree - see wikipedia_function_composition_. .. _wikipedia_function_composition: http://en.wikipedia.org/wiki/Function_composition Real world to real world transform '''''''''''''''''''''''''''''''''' We remind each other that a mapping is a function (callable) that takes coordinates as input and returns coordinates as output. So, if *M* is a mapping then:: [i',j',k'] = M(i, j, k) where the *i, j, k* tuple is a coordinate, and the *i', j', k'* tuple is a transformed coordinate. Let us imagine we have somehow come by a mapping *T* that relates a coordinate in a world space (mm) to other coordinates in a world space. A registration may return such a real-world to real-world mapping. Let us say that *V* is a useful mapping matching the voxel coordinates in *img1* to voxel coordinates in *img2*. If *img1* has a voxel to mm mapping *M1* and *img2* has a mm to voxel mapping of *inv_M2*, as in the previous example (repeated here):: imgA = load_image(infile_A) vx2mmA = imgA.coordmap imgB = load_image(infile_B) vx2mmB = imgB.coordmap mm2vxB = vx2mmB.inverse then the registration may return the some coordinate map, *T* such that the intended mapping *V* from voxels in *img1* to voxels in *img2* is:: mm2vxB_map = mm2vxB.mapping vx2mmA_map = vx2mmA.mapping V = mm2vxB_map.composewith(T.composedwith(vx2mmA_map)) To support this, there should be a CoordinateMap constructor that looks like this:: T_coordmap = mm2mm(T) where *T* is a mapping, so that:: V_coordmap = mm2vxB.composewith(T_coordmap.composedwith(vx2mmA)) I have done a coregistration between two images, *img1* and *img2*. This has given me a voxel-to-voxel transformation and I want to store this transformation in such a way that I can use this transform to resample *img1* to *img2*. :ref:`resampling` I have done a coregistration between two images, *img1* and *img2*. I may want this to give me a worldA-to-worldB transformation, where worldA is the world of voxel-to-world for *img1*, and worldB is the world of voxel-to-world of *img2*. My *img1* has a voxel to world transformation. This transformation may (for example) have come from the scanner that acquired the image - so telling me how the voxel positions in *img1* correspond to physical coordinates in terms of the magnet isocenter and millimeters in terms of the primary gradient orientations (x, y and z). I have the same for *img2*. For example, I might choose to display this image resampled so each voxel is a 1mm cube. Now I have these transformations: ST(*img1*-V2W), and ST(*img2*-V2W) (where ST is *scanner transform* as above, and *V2W* is voxel to world). I have now done a coregistration between *img1* and *img2* (somehow) - giving me, in addition to *img1* and *img2*, a transformation that registers *img1* and *img2*. Let's call this transformation V2V(*img1*, *img2*), where V2V is voxel-to-voxel. In actuality *img2* can be an array of images, such as series of fMRI images and I want to align all the *img2* series to *img1* and then take these voxel-to-voxel aligned images (the *img1* and *img2* array) and remap them to the world space (voxel-to-world). Since remapping is an interpolation operation I can generate errors in the resampled pixel values. If I do more than one resampling, error will accumulate. I want to do only a single resampling. To avoid the errors associated with resampling I will build a *composite transformation* that will chain the separate voxel-to-voxel and voxel-to-world transformations into a single transformation function (such as an affine matrix that is the result of multiplying the several affine matrices together). With this single *composite transformatio* I now resample *img1* and *img2* and put them into the world coordinate system from which I can make measurements.