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

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;
clc;
[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); 
end

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): http://apps.usd.edu/coglab/schieber/psyc707/pdf/2D-FFT.pdf ). 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?

Answer 1

[Michal]

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,

Michal

Answer 2

[Abhishek Ballaney]

https://in.mathworks.com/help/signal/ref/butter.html

[Alex Doruyter]

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.

Answer 3

[Alex Doruyter]

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,

Alex.

% 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
% 
clear
close all
clc;

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

       figure
       subplot(2,2,1)
       imagesc(img3D(:,:,middle_slice))
       title('original image')
       subplot(2,2,2)
       imagesc(log(1+abs(fftimg2D)))
       title('original 2D spectrum')
       subplot(2,2,3)
       imagesc(log(1+abs(fftimg2D_filtered)))
       title('filtered 2D spectrum')
       subplot(2,2,4)
       imagesc(img3D_filterd(1:dimx,1:dimy,middle_slice))
       title('filtered image')
    end


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

[Michal]

Thank you Alex, it may still be useful to me :) Good luck with your application!