Read, Process, and Write 3-D Medical Images

This example shows how to import and display volumetric medical image data, apply a smoothing filter to the image, and write the processed image to a new file. You can use the medicalVolume object to import image data and spatial information about a 3-D medical image in one object. This example also shows how you can use the properties and object functions of a medicalVolume object to work with image volumes in the intrinsic and patient coordinate systems.

Download Image Volume Data

This example uses one chest CT volume saved as a directory of DICOM files. The volume is part of a data set containing three CT scans. The size of the entire data set is approximately 81 MB. Download the data set from the MathWorks® website, then unzip the folder.

zipFile = matlab.internal.examples.downloadSupportFile("medical","MedicalVolumeDICOMData.zip");
filepath = fileparts(zipFile);
unzip(zipFile,filepath)
dataFolder = fullfile(filepath,"MedicalVolumeDICOMData","LungCT01");

Read Image File

The medicalVolume object imports data from the DICOM, NIfTI, and NRRD medical image file formats. DICOM volumes can be stored as a single file or as a directory containing individual files for each 2-D slice. The medicalVolume object automatically detects the file format and extracts the image data, spatial information, and modality from the file metadata. For this example, specify the data source as the download directory of the chest CT scan.

medVol = medicalVolume(dataFolder)

medVol = 
  medicalVolume with properties:

                 Voxels: [512×512×88 int16]
         VolumeGeometry: [1×1 medicalref3d]
           SpatialUnits: "mm"
            Orientation: "transverse"
           VoxelSpacing: [0.7285 0.7285 2.5000]
           NormalVector: [0 0 1]
       NumCoronalSlices: 512
      NumSagittalSlices: 512
    NumTransverseSlices: 88
           PlaneMapping: ["sagittal"    "coronal"    "transverse"]
               Modality: "CT"
          WindowCenters: [88×1 double]
           WindowWidths: [88×1 double]

The Voxels property contains the image intensity values. If the file metadata specifies the rescale intercept and slope, medicalVolume automatically rescales the voxel intensities to the specified units. In this example, the CT intensity values are rescaled to Hounsfield units.

vol = medVol.Voxels;

The VolumeGeometry property contains a medicalref3d object that defines the spatial referencing for the image volume, including the mapping between the intrinsic and patient coordinate systems. The intrinsic coordinate system is defined by the rows, columns, and slices of the Voxels array, with coordinates in voxel units. The patient coordinate system is defined relative to the anatomical axes of the patient, with real-world units such as millimeters.

R = medVol.VolumeGeometry

R = 
  medicalref3d with properties:

                 VolumeSize: [512 512 88]
                   Position: [88×3 double]
             VoxelDistances: {[88×3 double]  [88×3 double]  [88×3 double]}
    PatientCoordinateSystem: "LPS+"
               PixelSpacing: [88×2 double]
                   IsAffine: 1
              IsAxesAligned: 1
                    IsMixed: 0

Display Center Slice in Intrinsic Coordinates

Extract the center transverse slice from the image volume and display it using intrinsic image coordinates.

Calculate the index of the center slice as half of the number of transverse slices.

sliceTransverse = round(medVol.NumTransverseSlices/2);

Extract the voxel data for the center slice by using the extractSlice object function. The function orients the extracted slice imTransverse to match the standard radiological orientation. For example, when you display an extracted transverse slice, the anterior direction always points toward the top of the image, and the left anatomical direction points toward the right side of the image.

[imTransverse,positionTransverse,spacingsTransverse] = extractSlice(medVol,sliceTransverse,"transverse");

If the file metadata specifies a recommended display window to improve the contrast of the displayed image, the WindowCenters and WindowWidths properties of the medicalVolume object specify the center and width of the window, respectively. Calculate the minimum and maximum of the recommended display window for medVol.

windowMin = medVol.WindowCenters(1)-0.5*medVol.WindowWidths(1);
windowMax = medVol.WindowCenters(1)+0.5*medVol.WindowWidths(1);

Display the extracted slice. By default, the imshow function displays the image in intrinsic coordinates, in voxel units.

imshow(imTransverse,[windowMin windowMax]);
title("Center Transverse Slice, Intrinsic Coordinates")
axis on
xlabel("Right \rightarrow Left [voxels]")
ylabel("Posterior \rightarrow Anterior [voxels]")

Display Center Slice in Patient Coordinates

Display the center transverse slice in real-world units using patient coordinates.

Define the spatial referencing for the extracted slice. Get the limits of the slice in the patient coordinate system by using the sliceLimits object function.

[XLim,YLim,ZLim] = sliceLimits(medVol,sliceTransverse,"transverse");

Create an imref2d object, RTransverse, that specifies the limits of the slice in patient coordinates. The PlaneMapping property of medVol maps between the transverse, coronal, and sagittal planes and the xyz-axes of the patient coordinate system. "transverse" is the third element of PlaneMapping, meaning transverse slices are normal to the z-axis and parallel to the xy-plane. Therefore, XLim and YLim define the transverse slice limits.

RTransverse = imref2d(size(imTransverse),XLim,YLim);

Display the transverse slice, specifying the imref2d object to plot the image in patient coordinates.

imshow(imTransverse,RTransverse,[windowMin windowMax])
title("Center Transverse Slice, Patient Coordinates")
xlabel("Right \rightarrow Left [mm]")
ylabel("Posterior \rightarrow Anterior [mm]")

Smooth Voxel Intensity Data

Smooth the image with a 3-D Gaussian filter. Applying a Gaussian filter is one approach for reducing noise in medical images.

sigma = 2;
volSmooth = imgaussfilt3(vol,sigma);

Create a new medicalVolume object that contains the smoothed voxel intensities and preserves the spatial referencing of the original file. Create a copy of the original object medVol and set the Voxels property of the new object, medVolSmooth, to the smoothed image data.

medVolSmooth = medVol;
medVolSmooth.Voxels = volSmooth;

Extract and display the center slice of the smoothed voxel data. Use the same spatial referencing and display window information as the original object to plot the image in patient coordinates.

[imSmooth,position,spacings] = extractSlice(medVolSmooth,sliceTransverse,"transverse");
figure
imshow(imSmooth,RTransverse,[windowMin windowMax])
title("Center Transverse Slice, Smoothed")
xlabel("Right \rightarrow Left [mm]")
ylabel("Posterior \rightarrow Anterior [mm]")

Write Processed Data to New NIfTI File

Write the smoothed image data to a new NIfTI file by using the write object function. The function supports writing medical volume data in only the NIfTI file format.

niftiFilename = "LungCT01_smoothed.nii";
write(medVolSmooth,niftiFilename)

Read the new file using medicalVolume. Because the NIfTI file format does not contain metadata related to modality or display windows, the Modality property value is "unknown" and the WindowCenters and WindowWidths properties are empty.

medVolNIfTI = medicalVolume(niftiFilename);