Exploring Slices from a 3-Dimensional MRI Data Set

This example shows how to explore a volume of data by extracting slices through a three-dimensional MRI data set using imtransform and tformarray functions.

Step 1: Load and View Horizontal MRI

This example uses the MRI data set that comes with MATLAB® and that is used in the help examples for both montage and immovie. Loading mri.mat adds two variables to the workspace: D (128-by-128-by-1-by-27, class uint8) and a grayscale colormap, map (89-by-3, class double).

D comprises 27 128-by-128 horizontal slices from an MRI data scan of a human cranium. Values in D range from 0 through 88, so the colormap is needed to generate a figure with a useful visual range. The dimensionality of D makes it compatible with montage. The first two dimensions are spatial. The third dimension is the color dimension, with size 1 because it indexes into the colormap. (size(D,3) would be 3 for an RGB image sequence.) The fourth dimension is temporal (as with any image sequence), but in this particular case it is also spatial. So there are three spatial dimensions in D and we can use imtransform or tformarray to convert the horizontal slices to sagittal slices (showing the view from the side of the head) or coronal (frontal) slices (showing the view from the front or back of the head).

The spatial dimensions of D are ordered as follows:

• Dimension 1: Front to back of head (rostral/anterior to caudal/posterior)

• Dimension 2: Left to right of head

• Dimension 4: Bottom to top of head (inferior to superior).

An important factor is that the sampling intervals are not the same along the three dimensions: samples along the vertical dimension (4) are spaced 2.5 times more widely than along the horizontal dimensions.

Load the MRI data set and view the 27 horizontal slices as a montage.

load mri;
montage(D,map)
title("Horizontal Slices")

Step 2: Extract Sagittal Slice from Horizontal Slices Using imtransform

We can construct a midsagittal slice from the MRI data by taking a subset of D and transforming it to account for the different sampling intervals and the spatial orientation of the dimensions of D.

The following statement extracts all the data needed for a midsagittal slice.

M1 = D(:,64,:,:); size(M1)

ans = 1×4

   128     1     1    27

However we cannot view M1 as an image because it is 128-by-1-by-1-by-27. reshape (or squeeze) can convert M1 into a 128-by-27 image that is viewable with imshow.

M2 = reshape(M1,[128 27]); 
size(M2)

ans = 1×2

   128    27

imshow(M2,map)
title("Sagittal - Raw Data")

The dimensions in M2 are ordered as follows:

• Dimension 1: Front to back of head (rostral to caudal)

• Dimension 2: Bottom to top of head (inferior to superior).

We can obtain a much more satisfying view by transforming M2 to change its orientation and increase the sampling along the vertical (inferior-superior) dimension by a factor of 2.5—making the sampling interval equal in all three spatial dimensions. We could do this in steps starting with a transpose, but the following affine transformation enables a single-step transformation and more economical use of memory.

T0 = maketform("affine",[0 -2.5; 1 0; 0 0]);

The upper 2-by-2 block of the matrix passed to maketform, [0 -2.5;1 0], combines the rotation and scaling. After transformation we have:

• Dimension 1: Top to bottom of head (superior to inferior).

• Dimension 2: Front to back of head (rostral to caudal)

The call

imtransform(M2,T0,"cubic")

would suffice to apply T to M2 and provide good resolution while interpolating along the top to bottom direction. However, there is no need for cubic interpolation in the front to back direction, since no resampling will occur along (output) dimension 2. Therefore we specify nearest-neighbor resampling in this dimension, with greater efficiency and identical results.

R2 = makeresampler(["cubic","nearest"],"fill");
M3 = imtransform(M2,T0,R2);  
imshow(M3,map)
title("Sagittal - imtransform")

Step 3: Extract Sagittal Slice from the Horizontal Slices Using tformarray

In this step we obtain the same result as step 2, but use tformarray to go from three spatial dimensions to two in a single operation. Step 2 does start with an array having three spatial dimensions and end with an array having two spatial dimensions, but intermediate two-dimensional images (M1 and M2) pave the way for the call to imtransform that creates M3. These intermediate images are not necessary if we use tformarray instead of imtransform. imtransform is very convenient for 2-D to 2-D transformations, but tformarray supports N-D to M-D transformations, where M need not equal N.

Through its TDIMS_A argument, tformarray allows us to define a permutation for the input array. Since we want to create an image with:

• Dimension 1: Superior to inferior (original dimension 4, reversed)

• Dimension 2: Caudal to rostral (original dimension 1)

and extract just a single sagittal plane via the original dimension 2, we specify tdims_a = [4 1 2]. We create a tform via composition starting with a 2-D affine transformation T1 that scales the (new) dimension 1 by a factor of -2.5 and adds a shift of 68.5 to keep the array coordinates positive. The second part of the composite is a custom transformation T2 that extracts the 64th sagittal plane using a very simple INVERSE_FCN. Note that T2 and Tc take a 3-D input to a 2-D input.

T1 = maketform("affine",[-2.5 0; 0 1; 68.5 0]);  
inverseFcn = @(X,t) [X repmat(t.tdata,[size(X,1) 1])];
T2 = maketform("custom",3,2,[],inverseFcn,64);
Tc = maketform("composite",T1,T2);

We use the same approach to resampling as before, but include a third dimension.

R3 = makeresampler(["cubic","nearest","nearest"],"fill");

tformarray transforms the three spatial dimensions of D to a 2-D output in a single step. Our output image is 66-by-128, with the original 27 planes expanding to 66 in the vertical (inferior-superior) direction.

M4 = tformarray(D,Tc,R3,[4 1 2],[1 2],[66 128],[],0);

The result is identical to the previous output of imtransform.

imshow(M4,map)
title("Sagittal - TFORMARRAY")

Step 4: Create and Display Sagittal Slices

We create a 4-D array (the third dimension is the color dimension) that can be used to generate an image sequence that goes from left to right, starts 30 planes in, skips every other plane, and has 35 frames in total. The transformed array has:

• Dimension 1: Top to bottom (superior to inferior)

• Dimension 2: Front to back (rostral to caudal)

• Dimension 4: Left to right.

As in the previous step, we permute the input array using TDIMS_A = [4 1 2], again flipping and rescaling/resampling the vertical dimension. Our affine transformation is the same as the T1 above, except that we add a third dimension with a (3,3) element of 0.5 and (4,3) element of -14 chosen to map 30, 32, ... 98 to 1, 2, ..., 35. This centers our 35 frames on the mid-sagittal slice.

T3 = maketform("affine",[-2.5 0 0; 0 1 0; 0 0 0.5; 68.5 0 -14]);

In our call to tformarray, TSIZE_B = [66 128 35] now includes the 35 frames in the 4th, left-to-right dimension (which is the third transform dimension). The resampler remains the same.

S = tformarray(D,T3,R3,[4 1 2],[1 2 4],[66 128 35],[],0);

View the sagittal slices as a montage (padding the array slightly to separate the elements of the montage).

S2 = padarray(S,[6 0 0 0],0,"both");
montage(S2,map)
title("Sagittal Slices")

Step 5: Create and Display Coronal Slices

Constructing coronal slices is almost the same as constructing sagittal slices. We change TDIMS_A from [4 1 2] to [4 2 1]. We create a series of 45 frames, starting 8 planes in and moving from back to front, skipping every other frame. The dimensions of the output array are ordered as follows:

• Dimension 1: Top to bottom (superior to inferior)

• Dimension 2: Left to right

• Dimension 4: Back to front (caudal to rostral).

T4 = maketform("affine",[-2.5 0 0; 0 1 0; 0 0 -0.5; 68.5 0 61]);

In our call to tformarray, TSIZE_B = [66 128 48] specifies the vertical, side-to-side, and front-to-back dimensions, respectively. The resampler remains the same.

C = tformarray(D,T4,R3,[4 2 1],[1 2 4],[66 128 45],[],0);

Note that all array permutations and flips in steps 3, 4, and 5 were handled as part of the tformarray operation.

View the coronal slices as a montage (padding the array slightly to separate the elements of the montage).

C2 = padarray(C,[6 0 0 0],0,"both");
montage(C2,map)
title("Coronal Slices")