How can I apply a Butterworth filter to a 3D DICOM image?
1 次查看(过去 30 天)
显示 更早的评论
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?
0 个评论
回答(3 个)
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
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!