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

9 次查看（过去 30 天）
Alex Doruyter
Alex Doruyter 2018-1-30
评论： Michal 2022-5-18
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?

请先登录，再回答此问题。

回答（3 个）

Michal
Michal 2020-1-20
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
Abhishek Ballaney
Abhishek Ballaney 2018-1-31
https://in.mathworks.com/help/signal/ref/butter.html
  1 个评论
Alex Doruyter
Alex Doruyter 2018-1-31
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.

请先登录，再进行评论。

Alex Doruyter
Alex Doruyter 2020-1-30
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

类别

在 Help Center 和 File Exchange 中查找有关 DICOM Format 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by