Troubles With Image Resizing

Hello! I have a slice of a .nii file (essentially a 3D array) I'd like to print out. However, the dimensions of a voxel (a 3D pixel, so basically a cube) aren't exactly 1mm x 1mm x 1mm. They're around 1.2mm x 1.1mm x 1.1mm.

I took a slice of the 3D array, so it's just a 2D image. So, the dimensions of each pixel are ~1.2mm x 1.1mm. I'd like to resize the image such that one pixel is equal to one element of the array is equal to one mm. In theory, it doesn't seem too difficult. But I just can't figure this out.

Here's my code so far:

clear
clc
close all
% Load NIfTI image
info = niftiinfo('index_image.nii');
brain = niftiread(info);

% Extract a slice from the NIfTI image
slice = brain(:, :, 128);

% Create figure
hFig = figure;

% Display the slice with pixel-to-mm scaling
imagesc(slice);
axis off;
axis equal;
colormap(gray)

% Set figure size so that 1 pixel corresponds to 1mm on the printed figure
set(hFig, 'Units', 'centimeters', 'Position', [0 0 24 17.6]);
movegui(hFig, 'center');

% Set properties to control the output size
set(hFig, 'PaperUnits', 'centimeters');
set(hFig, 'PaperPosition', [0 0 24 17.6]);
set(hFig, 'PaperOrientation', 'landscape');

% Export to PDF and open file
print(hFig, '-dpdf', '-r0', 'out.pdf');
open('out.pdf'); 

The reason why use 24 and 17.6 is because my .nii file is 176 x 240 x 256 elements. I'd honestly like the code to work regardless of the size of the array, but even hard coding it into MATLAB isn't working. The code is compiling just fine, but the figure in the generated PDF doesn't have the dimensions I'd like it to have.

Would anyone be able to provide some insights? Thank you!

Answer 1

[Cris LaPierre]

If the Medical Imaging Toolbox is included in your license, I suggest using extractSlice to get the voxel information for your slice. The syntax is [X,position,spacings] = extractSlice(medVol,slice,direction)

imshow assumes equal spacing in all directions. When that is not the case, use the XData and YData inputs to specify the real world size of your image. Keep in mind that the first dimension of your array corresponds to Y, and the 2nd corresponds to X.

This code would display the image correctly.

brain = medicalVolume(info);
[img,position,spacings] = extractSlice(brain,128,'transverse') % replace direction accordingly
xMax = size(img, 2)*spacing(2);
yMax = size(img, 1)*spacing(1);
imshow(img, [], XData = [0 xMax], YData = [0 yMax])

Of course, if you don't have the Medical Imaging Toolbox, you could adapt the imshow code accordingly and still get the desired result.

[Cris LaPierre]

Another option is to use imwarp to resample your data to a 1x1x1 grid.

brain = medicalVolume('index_image.nii');
tform = intrinsicToWorldMapping(brain.VolumeGeometry);
worldBrain = imwarp(brain.Voxels, tform); % wholeBrain is resampled to voxelSpacing of 1x1x1

[Isabella]

After warping, how would I extract a slice and print it out?

Also, medicalVolume doesn't seem to preserve the file headers (which is why I was initally using niftiinfo), so the slice is all wonky in terms of orientation.

What does coordinate system does intrinsicToWorldMapping map to, exactly? I know the code can access the dimensions of the voxels via brain.VolumeGeometry, but where is it specified that I'd like the voxel to be remapped to have 1x1x1mm dimensions?

[Cris LaPierre]

The output of imwarp is an array, so you would need to use indexing (worldBrain(:,:,300)). Just be aware that your slice numbers will have changed to approximately

newSliceNums = ceil([176 240 256].*[1.2 1.1 1.1])

newSliceNums = 1x3

212 264 282

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Isabella]

brain = medicalVolume('index_image.nii');

tform = intrinsicToWorldMapping(brain.VolumeGeometry);

worldBrain = imwarp(brain.Voxels, tform);

I ran this code you sent, however, the dimensions of worldBrain are 293 x 272 x 301 voxels.

[Cris LaPierre]

That's ok. The point I was making is that the number of slices was going to change.

[Cris LaPierre]

unzip('anat.nii.zip')
brain = niftiread('anat.nii.gz');
info = niftiinfo('anat.nii.gz')

info = struct with fields:
                 Filename: '/users/mss.system.5vsm6/anat.nii.gz'
              Filemoddate: '20-Oct-2022 00:36:26'
                 Filesize: 2408857
                  Version: 'NIfTI1'
              Description: 'FreeSurfer Mar 11 2010'
                ImageSize: [256 256 128]
          PixelDimensions: [0.9375 0.9375 1.5000]
                 Datatype: 'int16'
             BitsPerPixel: 16
               SpaceUnits: 'Millimeter'
                TimeUnits: 'Second'
           AdditiveOffset: 0
    MultiplicativeScaling: 0
               TimeOffset: 0
                SliceCode: 'Unknown'
       FrequencyDimension: 0
           PhaseDimension: 0
         SpatialDimension: 0
    DisplayIntensityRange: [0 0]
            TransformName: 'Sform'
                Transform: [1x1 affine3d]
                  Qfactor: -1
                      raw: [1x1 struct]

img = squeeze(brain(120,:,:));
spacing = info.PixelDimensions(2:3)

spacing = 1x2

0.9375 1.5000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

tiledlayout(1,2)
nexttile
imshow(img,[])
nexttile
xMax = size(img,2)*spacing(2);
yMax = size(img,1)*spacing(1);
imshow(img, [], XData = [0 xMax], YData = [0 yMax])

Now recreate the image using imwarp. I can't use the medical imaging toolbox here, so here is the same image using the older approach you were attempting. I manually determined the corresponding slice number in the warped array.

tform = info.Transform;
newBrain = imwarp(brain,tform);

newimg = squeeze(newBrain(:,:,128));
newimg = permute(newimg,[2 1]);
figure
tiledlayout(1,2)
nexttile
imshow(img,[])
nexttile
imshow(newimg,[])
axis xy

The permute is because, in MATLAB, the first dimension of the image array corresponds to Y and the sencond dimention to X. Additionally, the Y axis is reversed in image axes, so axis xy changes that so the image is displayed as expected.

When the volume is warped to real-world coordinates, the dimensions get rearranged, which is why sagittal is now the 3rd dimension instead of the first. As a sanity check, compare the image sizes

% original - real-world size
newSliceNums = ceil(info.PixelDimensions.*info.ImageSize)

newSliceNums = 1x3

240 240 192

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% warped - transformed to a 1x1x1 grid
size(newBrain)

ans = 1x3

192 240 240

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

When using the Medical Imaging Toolbox, this would be the code

b = medicalVolume("anat.nii.gz");
[i,~,spacing] = extractSlice(b,120,'sagittal')

tiledlayout(1,2)
nexttile
imshow(i,[])
nexttile
xMax = size(i,2)*spacing(2);
yMax = size(i,1)*spacing(1);
imshow(i, [], XData = [0 xMax], YData = [0 yMax])

T = intrinsicToWorldMapping(b.VolumeGeometry)
bwarp = imwarp(b.Voxels,T);

newimg = squeeze(bwarp(:,128,:));
newimg = permute(newimg,[2 1]);
figure
tiledlayout(1,2)
nexttile
imshow(img,[])
nexttile
imshow(newimg,[])
axis xy

[Isabella]

Thank you so, so much for all the help. After messing around some more with the Image Processing Toolkit, I ended up using one of MATLAB's custom functions, preProcessBrainCANDIData, to resample the data.

Answer 2

[Image Analyst]

Maybe use interpn to resample the array along the "longer pixel" direction to have more samples, like the number of samples it would have if it were 1 mm instead of 1.2 mm.