1
$\begingroup$

I'm working on finding the center of rotation of a set of test tomographic projections in order to perform 2D reconstruction. I'm trying to implement the algorithm to find center of rotation as described in Vo et al. 2014, which involves Fourier transforming the sinogram and searching frequency space outside the "bow tie" double-wedge region for minimum artifacts that result from wrong center of rotation. However, when I Fourier transformed a sample sinogram generated from Shepp-Logan model, the resulting frequency space has lots of artifacts outside the double-wedge region even when center of rotation is correct. Through visual inspection I've also noticed that the amount of artifacts are sensitive to number of projections in the sinogram.

Phantom Sinogram Below: FT of sinogram with 200 projections, from 0 degree to 180 degree

Fourier transform of the sinogram Below: FT of sinogram with 180 projections, from 0 degree to 180 degree

FT of sinogram

The FFT'ed sinogram in the paper shows distinct zero regions outside the double wedge area, which is mathematically true as it is proven here. However, I'm not sure why my FFT'ed sinograms do not show similar zero regions as in the paper.

This is my first time posting here, so please let me know what other information is needed for others to help me. Thank you.

Below is my code


P=phantom(256);

angs=linspace(0,179,200);

pjs=radon(P,angs);

pjs=rot90(pjs,3);

pjs_inv=fliplr(pjs);

pjshift_inv=circshift(pjs_inv,[0 0]);

sinogram=[pjs; pjshift_inv];

figure(1),imagesc(sinogram);axis image,colormap gray

ft=real(fftshift(fft2(sinogram)));

figure(2), imshow(ft),colormap gray,axis image,axis on

$\endgroup$
1
$\begingroup$

Try imtool() and adjust the contrast to show ft similar as in the paper. I set eliminate outliers to 5%, then (contrast) window minimum to 0.0. May I suggest to ask the author about the settings for the display range in his paper.

display range adapted ft

Now to the signal processing.

The mentioned paper by Vo is public available. Moreover Vo published the algorithm in mathematica code. Note that the code by Vo does a gaussian filtering, i.e. a lowpass of the sinogram. IMHO to achieve similar results as in a real scan.

Generally a real scan shows lowpass behaviour, i.e. the high frequencies in the projections in the sinogram are lost due to the size of the focal point inside the x-ray tube, the detector system and radiation scatter inside the object. In contrast the forward projection algoritm: Those may keep high frequencies in the projections, depending on their base algorithm.

sinogram1 = GaussianFilter[ sinogram1, {{3,2}}];
(* Apply smooth filter: recommend blur more in the vertical direction *)

The Mathematica help tells, that those 3 and 2 steer the xy-sigma of the GaussianFilter(). One can go theroretically through that to find the values of the matrix, but I prefer to run the above filter over a dirac impulse to see the resulting filtermatrix and establish it similar in matlab for use in imfilter(). But I have no access to mathematica. May be s.o can translate the Mathematica GaussianFilter() into a matlab fspecial().

So let's filter the sinogram with a lowpass to ged rid of the high frequencies introduced by radon() :

imtool(ft); %% eliminate 5% outliers and set min to 0.0
h=fspecial('gaussian',[2 3]); %% wider in projection direction
sino_lp = imfilter(sinogram,h,'replicate');
ft_lp=real(fftshift(fft2(sino_lp)));
imtool(ft_lp);

Such a fourier transformed sinogram has lower values in the high frequencies.

$\endgroup$
1
  • $\begingroup$ @lennon310: May I ask to add tags like "sinogram" and "interpolation" $\endgroup$
    – Stephan
    May 22 '15 at 16:46
1
$\begingroup$

The chapter "Alaising Artifacts And Noise In CT Images" gives in formula #10 for an equal radial and channel resolution that the number of parallel projections is pi/2 times the number of projection channels for 360°, or pi/4 for 180°. In order to reconstruct an image, a fan beam system needs data from 180° + the fan angle, e.g. 50°. Now (230°/180°) * pi/4 is about 1.0, i.e. a fan beam system needs about the same amount of fan projections as it has channels.

P has dimension 256x256 and radon() forward projects this into 367 ~ 256xsqrt(2) projection channels, i.e. the number of channels becomes the diameter of the square given by P. Please checkout the radon() help.

Formula #10 derives from those 367 then 550 projections onwards for 360°. I round this down to the next power of 2 to achieve better fft2() performance to 256 angles on 180°, but radon() will take longer for those 256 than your given 200 angels.

Extending the sinogram to 512 channels and 512 'repeating' projections over pi reduces alaising in the fourier space.

Note that your angs are not evenly distributed in 360°, but have a gap close to 180°.

...
angs=linspace(0,180,257);
angs(end)=[]; %% now 256 even distributed angles in range [0..180[
...
ft=real(fftshift(fft2(sinogram,512,512)));
...
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.