How can I apply a Butterworth filter to a 3D DICOM image?
4 views (last 30 days)
I want to apply a Butterworth filter to a 3d image (single file, multiframe DICOM), in this case a perfusion SPECT of the lungs.
%Apply a Butterworth filter to a volumetric (3D) DICOM image.
[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);
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?
Michal on 20 Jan 2020
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.