Read, Process, and Write 3-D Medical Images

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.

Run this code to download the data set from the MathWorks® website and unzip the folder.

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

The dataFolder folder contains the downloaded and unzipped data.

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.

intensities = 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 Volume as Slices

Medical Imaging Toolbox™ provides several options for visualizing medical image data. For details, see Choose Approach for Medical Image Visualization. For this example, display the transverse slices of the CT volume by creating a sliceViewer object. By default, the viewer opens to display the center slice. Use the slider to navigate through the volume.

sliceViewer(medVol,Parent=figure)
title("CT Volume, Transverse Slices")

Display the Volume in 3-D

Display the CT volume as a 3-D object by using the volshow function. The volshow function uses the spatial details in medVol to set the Transformation property of the output Volume object and display the volume in the patient coordinate system.

You can customize the volume display by setting properties of the Volume object. Specify a custom transparency map and colormap that highlight the rib cage. The alpha and color values are based on the CT-bone rendering style from the Medical Image Labeler app. The intensity values have been tuned for this volume using trial and error.

alpha = [0 0 0.72 1.0];
color = [0 0 0; 186 65 77; 231 208 141; 255 255 255] ./ 255;
intensity = [-3024 100 400 1499];
queryPoints = linspace(min(intensity),max(intensity),256);
alphamap = interp1(intensity,alpha,queryPoints)';
colormap = interp1(intensity,color,queryPoints);

Display the volume with the custom display settings. Specify the cinematic rendering style, which displays the volume with realistic lighting and shadows.

vol = volshow(medVol,Colormap=colormap,Alphamap=alphamap,RenderingStyle="CinematicRendering");

Pause to apply all of the cinematic rendering iterations before updating the display in Live Editor.

pause(3.5)
drawnow

Optionally, you can clean up the viewer window by using the 3-D Scissors tool, , to remove the patient bed. For a detailed example, see Remove Objects from Volume Display Using 3-D Scissors. This image shows the final volume display after removing the patient bed and rotating to view the anterior plane.

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;
intensitiesSmooth = imgaussfilt3(intensities,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 = intensitiesSmooth;

Display the smoothed voxel data.

figure
sliceViewer(medVolSmooth,Parent=figure)
title("Smoothed CT Volume, Transverse Slices")

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);