%---TeX-start-of-header--- \documentclass{report} \usepackage{epsfig,mathematica,figi} \mmanobreak \begin{comment} \begin{mathematica} zero = 0.00001; \ ohide = DisplayFunction->Identity; \ oshow = DisplayFunction -> $DisplayFunction; \ ojoin = PlotJoined -> True; \ oall = PlotRange->All; \ osurf = HiddenSurface->False; \ deg = N[Degree]; \ pi = N[Pi]; \end{mathematica} \begin{mathematica} <{{-1,1},{-1,1}}, AspectRatio->1 ]; \end{mathematica} Function \texttt{InEllipse} returns the density value for given ellipse and a pair of coordinates; it they are out of this ellipse function returns 0. \begin{mathematica} InEllipse = Compile[{x0, y0, a, b, alpha, rho, x, y}, Module[{x1, y1}, x1 = x - x0; y1 = y - y0; If[((x1 Cos[alpha] + y1 Sin[alpha]) / a)^2 + ((-x1 Sin[alpha] + y1 Cos[alpha]) / b)^2 <= 1., rho, 0. ] ]]; \end{mathematica} The following function calculates the total sum of ellipses densities at a given point. \begin{mathematica} Model[modeldata_List, x_, y_] := Plus @@ ((InEllipse @@ Join[#, {x, y}])& /@ modeldata); \end{mathematica} Function \texttt{SampledModel} returns an array of values of the model sampled on grid $n \times n$ over a range $[-w, w]$. \begin{mathematica} SampledModel[modeldata_List, n_Integer, w_] := Table[ Model[modeldata, x, y], {y, -w, w, 2w/(n-1)}, {x, -w, w, 2w/(n-1)}]; \end{mathematica} %------------------------------------------------------------------------ \section{Projections} Function \texttt{EllipseProjection} finds the value of the projection at angle $th$ at point $t$ for an ellipse defined by the set of parameters. \begin{mathematica} EllipseProjection = Compile[{x0, y0, a, b, alpha, rho, theta, t}, Module[{angle, s, t1}, s = Sqrt[x0^2 + y0^2]; t1 = If[s != 0, (t - s Cos[N[ArcTan[x0, y0]] - theta])^2, t^2]; angle = theta - alpha; at = (a Cos[angle])^2 + (b Sin[angle])^2; If[t1 <= at, 2. rho a b / at Sqrt[at - t1] , 0.] ]]; \end{mathematica} Function \texttt{ModelProjection} finds the resultant projection values after using all ellipses in the model described by $modeldata$ parameter. \begin{mathematica} ModelProjection[modeldata_List, th_, t_] := Plus @@ ((EllipseProjection @@ Join[#, {th, t}])& /@ modeldata); \end{mathematica} Function \texttt{Slice} calculates a complete set of $K$ projections using $n$ rays for projection width $w$. \begin{mathematica} Slice[modeldata_List, K_Integer, n_Integer, w_] := Table[ ModelProjection[modeldata, th, t], {th, 0., (K - 1)/K pi, pi / K}, {t, -w, w, 2 w/(n - 1)} ]; \end{mathematica} %------------------------------------------------------------------------ \section{Filtering} The following function calculates values of the $h_1$ sequence. \begin{mathematica} h1[n_Integer, tau_] := If[n == 0, 1/(4. tau^2), If[OddQ[n], -1/(n pi tau)^2, 0.] ]; \end{mathematica} Just to check if it is justified to skip IFFT/FFT steps proposed for calculation of $H_1$ in Kak let us plot $|w|$ versus FFT transform of the $h_1$ sequence: \begin{mathematica} tau = .1; \ h1seq = Pad[Table[h1[n, tau], {n, -64, 64}], 256]; \ PSTeX[ sh[Plot[Abs[w], {w, -1/(2 tau), 1/(2 tau)}], ListPlot[rot[tau FFT[h1seq], 1/tau], ojoin, oall]], "h1compare"]; \end{mathematica} \figi{h1compare}{Comparison of frequency sequences obtained by direct calculation and by FFT.} They are nearly identical, as a matter of fact it is the sampled version which is slightly inaccurate. It follows that one could use the sequence $H_1$ generated by the following function, but I will not do it, because then it would be hard to estimate aliasing influence from projection signal. \begin{mathematica} H1[n_Integer] := Table[0.5 - Abs[w - 0.5], {w, 0., 1. - 1./n, 1./n}]; \end{mathematica} Function \texttt{FilterSlice} calculates filtered values of all projections. \begin{mathematica} FilterSlice[slice_List, w_] := Module[{n, H, tau}, n = Length[slice[[1]]]; tau = 2 w / (n - 1); H = tau FFT @ ZeroPad[ Table[h1[i, tau], {i, -n/2, n/2 - 1}], 2 n]; (IFFT[ FFT[ ZeroPad[#, 2 n]] H])& /@ slice ]; \end{mathematica} %------------------------------------------------------------------------ \section{First attempt at reconstruction} Function \texttt{InterpolateSlice} for a given filtered slice data $fsl$ returns an array of $K$ interpolation functions of order $order$ defined over projection width $w$. \begin{mathematica} InterpolateSlice[fsl_List, order_Integer, w_] := Module[{l, idx}, l = Length[fsl[[1]]]; idx = Table[2. w (2. i - l - 1)/(l - 1), {i, 1, l}]; Interpolation[ Transpose[{idx, #}], InterpolationOrder->order]& /@ fsl ]; \end{mathematica} The following function is the visible trace of my futile attempts to speed the calculations up. \begin{mathematica} cossin = Compile[{x, y, th}, x Cos[th] + y Sin[th]]; \end{mathematica} Function \texttt{Reconstruction} calculates values of reconstructed model at point $(x,y)$ for the given interpolated slice data. \begin{mathematica} Reconstruction[isl_List, x_, y_] := Module[{K = Length[isl]}, pi / K Sum[ isl[[k+1]] [cossin[x, y, pi k / K]], {k, 0, K-1}] ]; \end{mathematica} The three above functions work just great, but really slooow, so below there is my next attempt. %------------------------------------------------------------------------ \section{Second and faster attempt at reconstruction} Function \texttt{ind} finds an index of the element of the sequence which is closest to the desired position. \begin{mathematica} ind = Compile[{x, y, w, th, {n, _Integer}}, Round[ ((n - 1) (x Cos[th] + y Sin[th]) / w + 1 + n) / 2] ]; \end{mathematica} Function \texttt{Reconstruction} calculates the value of the reconstructed image at the point $(x,y)$. \begin{mathematica} Reconstruction[fsl_List, w_, x_, y_] := Module[ {K = Length[fsl], n = Length[fsl[[1]]]}, pi / K Sum[ fsl[[k+1, ind[x, y, w, pi k / K, n] ]], {k, 0, K-1}] ]; \end{mathematica} Function \texttt{SampledReconstruction} given an interpolated slice returns an array of $n\times n$ values approximating image. \begin{mathematica} SampledReconstruction[isl_List, n_Integer, w_] := Table[ Reconstruction[isl, w, x, y], {x, -1, 1, 2/(n-1)}, {y, -1, 1, 2/(n-1)} ]; \end{mathematica} \chapter{Results} \section{Shepp-Logan model} The variable \texttt{modeldata} keeps data about the Shepp-Logan head model. \begin{mathematica} (modeldata = {{ 0, 0, 0.69, 0.92, 0 deg, 1}, {0, -0.0184, 0.6624, 0.874, 0 deg, -0.98}, {0.22, 0, 0.11, 0.31, -18 deg, -0.02}, {-0.22, 0, 0.16, 0.41, 18 deg, -0.02}, {0, 0.35, 0.21, 0.25, 0, 0.01}, {0, 0.1, 0.046, 0.046, 0, 0.01}, {0, -0.1, 0.046, 0.046, 0, 0.01}, {-0.08, -0.605, 0.046, 0.023, 0, 0.01}, {0, -0.605, 0.023, 0.023, 0, 0.01}, {0.06, -0.605, 0.023, 0.046, 0 deg, 0.01}}) // MatrixForm . Out[22]//MatrixForm= > 0 0 0.69 0.92 0 1 0 -0.0184 0.6624 0.874 0 -0.98 0.22 0 0.11 0.31 -0.314159 -0.02 -0.22 0 0.16 0.41 0.314159 -0.02 0 0.35 0.21 0.25 0 0.01 0 0.1 0.046 0.046 0 0.01 0 -0.1 0.046 0.046 0 0.01 -0.08 -0.605 0.046 0.023 0 0.01 0 -0.605 0.023 0.023 0 0.01 0.06 -0.605 0.023 0.046 0 0.01 \end{mathematica} \begin{comment} \begin{mathematica} PSTeX[showmodel[modeldata], "model"]; \end{mathematica} \end{comment} \figi{model}{The Shepp-Logan head phantom.} As one can see we can assume that this model is defined over $[-1, 1]\times [-1, 1]$ square. \begin{mathematica} W = 1.; \end{mathematica} The following plot visualizes projections obtained for ellipse ``e'' of the model: \begin{mathematica} PSTeX[Plot3D[EllipseProjection @@ Join[ modeldata[[5]], {th, t}], {t, -W, W}, {th, 0., pi}, PlotPoints->{50,32}, AxesLabel->{"t", "theta", "p"}, Compiled->False], "eproj"]; \end{mathematica} \figi{eproj}{Plot of projections of the ellipse ``e'' in function of the angle of the projection.} The function \texttt{SampledModel} is used for creation of a $128\times 128$ bitmap image of the model. \begin{mathematica} WritePGM[Reverse @ Round[255 normalize [SampledModel[modeldata, 128, W]]], "smodel"]; \end{mathematica} \figi{smodel}{Bitmap image ($128\times 128$) of the model (after histogram equalization).} The $128 \times 128$ reconstruction is shown in Figure~\ref{sreco}. \figi{sreco}{Sampled reconstruction of the Shepp-Logan model (after histogram equalization).} The comparison of the values on the $y = -0.605$ line for original and reconstructed model is shown in Figure~\ref{comp}. \begin{comment} \begin{mathematica} PSTeX[Plot[{ Reconstruction[f32, 2., x, -0.605], Model[modeldata, x, -0.605]}, {x, -W, W}, PlotPoints->100(*, PlotRange->{0, .04}*) ], "comp"]; \end{mathematica} \end{comment} \figi{comp}{Comparison of true and reconstructed values on the $y=-0.605$ line.} It is easy to notice that reconstructed values are slightly higher than these obtained by Kak. But he simply didn't write by what means he obtained his image -- his mysterious remarks about windowing and interpolation suggest that he used both techniques for improving image. This model is very demanding because the narrow rim of high values representing skull definitely increases the small values inside the head what presumably is the main reason of obtained distortions. \section{Demonstration of aliasing effects} Variable \texttt{ellipsedata} keeps the definition of the model consisting of single ellipse with density 1. \begin{mathematica} ellipsedata = {{0, 0, 0.15, 0.25, 0, 1}}; \end{mathematica} \begin{mathematica} PSTeX[showmodel[ellipsedata], "emodel"]; \end{mathematica} \figi{emodel}{Model used for demonstration of aliasing effects.} The following function calculates the sampled reconstruction of the above model using $K$ projections and $n$ rays. \begin{mathematica} erec[K_, n_] := SampledReconstruction[ FilterSlice[ Slice[ellipsedata, K, n, W], W], 128, 2.]; \end{mathematica} Next statement writes whole set of images in sizes from $32\times 32$ to $512\times 512$. It took approximately 50 hours of CPU time to calculate it (I distributed it among several machines). \begin{mathematica} WritePGM[Transpose @ Reverse @ Round[255 normalize[ erec[#[[1]], #[[2]]] ]], "e" <> ToString[#[[1]]] <> "_" <> ToString[#[[2]]]] & /@ Distribute[{#,#}, List]&[{32, 64, 128, 256, 512}]; \end{mathematica} The results are shown in Figure~\ref{alias}. \newcommand{\efig}[2]{\psfig{figure=e#1_#2,width=1.5in}\,} \newcommand{\eline}[1]{\efig{64}{#1}\efig{128}{#1}\efig{256}{#1}\efig{512}{#1}\smallskip} \begin{figure}[p!] \centering \eline{32} \eline{64} \eline{128} \eline{256} \eline{512} \caption{Twenty $128\times 128$ reconstructions of the single ellipse model for values of $n$ changing from 64 to 512 (columns) and values of $K$ changing from 32 to 512 (rows).} \label{alias} \end{figure} The results are completely consistent with Figure~21 of the reference. Quite good reconstructions are obtained already for the case $K=256$ and $n=256$. The pictures look so bad beacuse I used histogram equalization (not simple clipping) to make the streaks better visible, and the printing was performed with 300dpi resolution, what is simply not enough in this case. \end{document} \begin{mathematica} Timing[s32 = Slice[modeldata, 100, 128, W];] \end{mathematica} \begin{mathematica} ListPlot3D[s32,osurf]; \end{mathematica} \begin{mathematica} f32 = FilterSlice[s32, W]; \end{mathematica} \begin{mathematica} ListPlot3D[f32,osurf]; \end{mathematica} \begin{mathematica} WritePGM[Transpose @ Reverse @ Round[255 normalize [SampledReconstruction[f32, 128, 2.]]], "sreco"]; \end{mathematica}