How can I apply a Butterworth filter to a 3D DICOM image?

15 views (last 30 days)
Alex Doruyter
Alex Doruyter on 30 Jan 2018
Commented: Michal on 18 May 2022
I want to apply a Butterworth filter to a 3d image (single file, multiframe DICOM), in this case a perfusion SPECT of the lungs.
My script:
%Apply a Butterworth filter to a volumetric (3D) DICOM image.
clear all;
[FileName,PathName] = uigetfile('*.dcm', 'Select studies to filter','MultiSelect','on');
cutoff = 0.44;
order = 5;
[b,a] = butter(order,cutoff);
FileName = cellstr(FileName);
filenumber = size(FileName,1);
for f = 1:filenumber
location = char(fullfile(PathName,FileName(f)));
Y = dicomread(location);
metadata = dicominfo(location);
P = squeeze(Y);
F = fftn(P);
F = filter(b,a,F);
G = ifftn(F);
G = uint16(G);
G = real(G);
G = reshape(G,[64 64 1 36]);
metadata.StudyDescription = strcat('bw',metadata.StudyDescription);
dicomwrite(G,'test.dcm',metadata,'CreateMode','Copy', 'MultiframeSingleFile',true);
If I comment out the filter the script runs smoothly (input file = output file) but as soon as the filter is introduced is returns a distorted image (which I suspect is a wrap-around artefact as described here (pg5): ). My attempts to pad (before the FFT) and subsequently crop (after extracting the real components of the array post-IFFT) a 3D matrix have been unsuccessful. Any advice?

Answers (3)

Michal on 20 Jan 2020
Hi Alex,
not sure it's relevant any more... after this long :) but still:
as far as I can see in filter function documentation, this function works along the first non-singleton dimension of the data.
This means that if your fft in F is 128x128x128, the filter is applied only along the rows. I'm not sure how exacly this transfers from frequency to spatial domain, but it's probably cause fro trouble, because you filter only partially, and the high frequency content remains in the other 2 dimensions.... At least thats how I see it.
Maybe you can try filtering once along each dimension, shifting the dimensions each run 1:3... but I'm not sure BW filter is separable this way.
I'm looking for a 3D BW filter myself, or I'd point you to one.
Good luck,

Abhishek Ballaney
Abhishek Ballaney on 31 Jan 2018
  1 Comment
Alex Doruyter
Alex Doruyter on 31 Jan 2018
Thanks Abhishek. I already use the butter command and I've seen the documentation. My issue relates to the distorted output I get from applying that filter. I think I need to zero pad my original 3D matrix (and then crop it back again later) but have been unsuccessful in this.

Sign in to comment.

Alex Doruyter
Alex Doruyter on 30 Jan 2020
Hi Michal
Thanks for your reply. I eventually addressed my application with the assistance of Patrick Dupont from KUL. The solution is not truly a 3D BW filter, since it only filters the 2D slices of a 3D volume, but it turned out that this was actually what I needed in any case. In the event that the script is useful to you, I post it below.
Best wishes,
% TB_Butterworth
% This script will read a 3D DICOM image, apply a butterworth filter and
% saves the filtered image again as DICOM
% author: Alex Doruyter and Patrick Dupont
% date: March 2018
% history:
% @(#)TB_Butterworth.m v0.01 last modified: 2018/03/19
close all
%------- Settings ---------------------------------------------------------
cutoff = 0.44; % cut off frequency (/cm) of butterworth filter
order = 5; % order of the butterworth filter
B = 1;
figures_on = 0; % if set to 1, the middle slice will be shown for QC purposes
%-------- END of Settings -------------------------------------------------
[name,path] = uigetfile('*.dcm', 'Select image to filter','Multiselect','on');
filename = fullfile(path,name);
for f = 1:length(filename)
% read dicom data
location = char(filename(f))
Y = dicomread(location); % now we have the DICOM image
metadata = dicominfo(location); %now the metadata of the original file is known
img3D = squeeze(Y);
[dimx,dimy,dimz] = size(img3D);
newdim = 2*[dimx dimy dimz] + 1; % to include zeropadding
% read the voxel size
tmp = metadata.PixelSpacing;
size_x = tmp(1);
size_y = tmp(2);
% calculate cut-off taking into account the voxel size (assuming isotropic
% voxels)
cutoff_per_pixel = size_x*cutoff/10;
d = cutoff_per_pixel*newdim(1);
fftimg3D = fftshift(fftn(double(img3D),newdim));
[x,y,z] = meshgrid(-floor(newdim(1)/2):floor(newdim(1)/2),-floor(newdim(2)/2):floor(newdim(2)/2),-floor(newdim(3)/2):floor(newdim(3)/2));
D = sqrt(x.^2 + y.^2 + z.^2); % Distance to the centre
butterworth = 1./(1+B*((D./d).^(2*order))); % calculate the butterworth filter
fftimg3D_filtered = fftimg3D .* butterworth; % filter in the frequency domain
img3D_filterd = real(ifftn(ifftshift(fftimg3D_filtered))); % transform back to the spatial domain
if figures_on == 1
middle_slice = round(dimz/2);
img = img3D(:,:,middle_slice);
fftimg2D = fftshift(fft2(double(img),newdim(1),newdim(2)));
[x2,y2] = meshgrid(-floor(newdim(1)/2):floor(newdim(1)/2),-floor(newdim(2)/2):floor(newdim(2)/2));
D2 = sqrt(x2.^2 + y2.^2); % Distance to the centre
butterworth2D = 1./(1+B*((D2./d).^(2*order)));
fftimg2D_filtered = fftimg2D .* butterworth2D;
title('original image')
title('original 2D spectrum')
title('filtered 2D spectrum')
title('filtered image')
% reintroduce the fourth dimension of size "1" to indicate that the DICOM should be read as greyscale
outputimg = reshape(img3D_filterd(1:dimx,1:dimy,1:dimz),[dimx dimy 1 dimz]);
% change the outputimg to fit again with the original dataformat
format_class = class(Y);
eval(['outputimg = ' format_class '(outputimg);']);
%Provide a new descriptor indicating that the data has been filtered with a Butterworth filter.
metadata.StudyDescription = strcat('bw',metadata.StudyDescription);
%Provide a new prefix to the file to indicate it has been filtered.
newfilename = strcat(location(1:end-4),'_bw.dcm');
%save new version of (filtered) dicom
dicomwrite(outputimg,newfilename,metadata,'CreateMode','Copy', 'MultiframeSingleFile',true);
  1 Comment
Michal on 18 May 2022
Thank you Alex, it may still be useful to me :) Good luck with your application!

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!