Info on the CT-data from an older generation Siemens Somatom Scanner A. Faridani, (faridani@math.orst.edu), May 10, 2007 The data were collected in 1986. The data are in the binary files phantom.ctd (resolution phantom) pelvis.ctd (human pelvis cross-section) Some explanations: Fan-beam scanning geometry, as described in Natterer's 1986 book 'The Mathematics of Computerized Tomography'. Radius of source circle RS = 2.868. 512 rays per source position. The data are 2 byte integers, stored in records of 1024 bytes, i.e. each record consists of the 512 measurements corresponding to one source position. phantom.ctd: 720 source positions, equispaced over 360 degrees pelvis.ctd: 360 source positions, equispaced over 360 degrees Numbering of rays is from left to right when standing at the source and looking towards the origin. First ray is tangent to unit circle, 257th ray passes through origin, rays are angularly equispaced. Spacing is d_alpha=arcsin(1./RS)/256. The whole fan is most likely shifted by d_alpha/4. The following MATLAB routine illustrates how the data can be read and a reconstruction be computed. It does not take a shift by d_alpha/4 into account. The program is given for illustrative purposes and should not be expected to yield the best possible reconstruction. Below are three M-files: fbpreadat.m is the main routine; slkernel.m implements the Shepp-Logan convolution kernel; and window3.m displays the reconstructed image. These three M-files should be stored separately in the same directory as the data files phantom.ctd or pelvis.ctd. With the current settings, typing the command >> fbpreadat at the Matlab prompt will then produce a reconstruction from the projection data in phantom.ctd. %----- Cut here for fbpreadat.m ------------------------ %fbpreadat.m %Fan-beam filtered backprojection algorithm % for the standard lattice % Data are read from a file % Last revision: May 10, 2007 % % Reconstructed image is stored in matrix P. % % Examples for input parameters: % Parameter for data from old generation Siemens scanner: %Pelvis data: datfile = 'pelvis.ctd' % p= 360, q=256, R=2.868, % roi=[-.9 .9 -.9 .9]; % %Resolution phantom: datfile = 'phantom.ctd' % p=720, q=256, R=2.868 % roi=[-.5 .5 -.5 .5]; %SPECIFY INPUT PARAMETERS HERE: datfile = 'phantom.ctd'; p=720 %number of view angles between 0 and 2*pi q=256 %d=arcsin(1/R)/q, d = angular detector spacing R = 2.868 % Radius of source circle. % Support of object is contained in unit circle MX=256; MY = 256; %matrix dimensions for image roi=[-.5 .5 -.5 .5]; %roi=[xmin xmax ymin ymax] %region of interest where %reconstruction is computed circle = 1; % If circle = 1 image computed only inside % circle inscribed in roi. %END OF INPUT SECTION b=0.5*pi*q/asin(1/R) fid = fopen(datfile,'r'); rps=1/b; if MX > 1 hx = (roi(2)-roi(1))/(MX-1); xrange = roi(1) + hx*[0:MX-1]; else hx = 0; xrange = roi(1); end if MY > 1 hy = (roi(4)-roi(3))/(MY-1); yrange = flipud((roi(3) + hy*[0:MY-1])'); else hx = 0; yrange = roi(3); end center = [(roi(1)+roi(2)), (roi(3)+roi(4))]/2; x1 = ones(MY,1)*xrange; %x-coordinate matrix x2 = yrange*ones(1,MX); %y-coordinate matrix if circle == 1 re = min([roi(2)-roi(1),roi(4)-roi(3)])/2; chi = ((x1-center(1)).^2 + (x2-center(2)).^2 <= re^2); %char. fct of roi; else chi = isfinite(x1); end x1 = x1(chi); x2 = x2(chi); P = zeros(MY,MX);Pchi = P(chi); h = asin(1/R)/q; beta = h*[-q:q-1]; s = R*sin(beta); cosbeta = cos(beta); bs = b*sin(h*[-2*q:2*q-1]); wb = slkernel(bs)/(rps^2); %compute discrete convolution kernel. for j = 1:p if (mod(j,10) == 0) j end alphaj = (2*pi*(j-1)/p); aj = R*[cos(alphaj);sin(alphaj)]; phi = alphaj + pi/2 - beta; % RF = fread(fid,2*q,'integer*2'); RF = fread(fid,2*q,'integer*2','ieee-be'); RF = flipud(RF)'; RF = RF.*cosbeta; % convolution C = conv(RF,wb); Q = h*C(2*q+1:4*q); Q(2*q+1)=0; % interpolation and backprojection Q = [real(Q)'; 0]; amx1 = aj(1)-x1; amx2 = aj(2)-x2; amxabsq = amx1.^2 + amx2.^2; %||a-x||^2 t = real(acos( (amx1*aj(1) + amx2*aj(2))./(R*sqrt(amxabsq)) )); ik = find(-aj(2)*x1 + aj(1)*x2 < 0); t(ik)=-t(ik); k1 = floor(t/h); u = (t/h-k1); k = max(1,k1+q+1); k = min(k,2*q); Pupdate = ((1-u).*Q(k)+u.*Q(k+1)); Pchi=Pchi+ Pupdate./amxabsq; end % j-loop P(chi) = Pchi*(R*2*pi/p); % Display image pmin = min(min(P)); pmax = max(max(P)); window3(pmin,pmax,roi,P); % view the computed image %----- Cut here for slkernel.m ------------------------ function u = slkernel(t) u = zeros(size(t)); i1 = abs(abs(t)-pi/2)<=1.e-6; u(i1) = ones(size(u(i1)))/pi; t1 = t(abs(abs(t)-pi/2)>1.e-6); v = (pi/2 - t1.*sin(t1))./((pi/2)^2 - t1.^2); u(abs(abs(t)-pi/2)>1.e-6) = v; u = u/(2*pi^3); % --- Cut here for window3.m function pic1 = window3(mi,ma,roi,pic); %function pic1 = window3(mi,ma,roi,pic); % displays image pic with coordinates given by roi % roi = [xmin xmax ymin ymax] x = [roi(1), roi(2)]; y = [roi(3), roi(4)]; colors = 128; co = colors-1; pic1 = pic - mi*ones(size(pic)); pic1 = (co/(ma-mi))*pic1; P = (pic1 >= 0); pic1 = pic1.*P; P = (pic1 <= co); pic1 = pic1.*P + co*(ones(size(pic1)) - P); colormap(gray(colors)); image(x,fliplr(y),flipud(pic1)); axis('square');