Hi All,

I am new to Julia, coming from Matlab, recently started working on OPT (Optical Projection Tomography) project where one of requirements is open source.

I tried to implement something similar to Matlab radon/iradon pair, using as a base the C++ code developed previously by T. Knopp (https://github.com/tknopp/numcpp).

For debugging, i mirrored Julia implementation of fbp in Matlab.

Quality-wise, sinograms and reconstructions of phantom are in good agreement with each other and with original.

(there is one minor scale issue on iradon, as I didnâ€™t find suitable â€śimresizeâ€ť, but that can be fixed later)

The main drawback is that the â€śnaĂŻveâ€ť Julia functions are much slower than optimized Matlabâ€™s.

For 1000x1000 phantom image, the ratios of processing times were

radon - 16.7

iradon - 20.5

iradon translated to Matlab from Julia - 23.5

That is not practical for the project.

The main issue looks like simply slow brackets.

I paste the code and test program below, hope to hear any suggestion on speed-up.

(Itâ€™s amateurish but sort of does the job).

I think sooner or later the question of radon/iradon will appear in Julia anyway.

```
# tomography.jl
########################################################
function radon(A,angles)
angles = angles + pi/2;
w,h = size(A);
Nr_div2 = floor(sqrt(w*w+h*h)/2);
Nr = Int64(Nr_div2*2+1);
L = length(angles);
sinogram = zeros(Nr,L);
SIN = sin(angles);
COS = cos(angles);
R = linspace(-0.5,0.5,Nr)*Nr;
RSIN = zeros(Nr,L);
RCOS = zeros(Nr,L);
w2 = w/2;
h2 = h/2;
for a=1:L
for k=1:Nr
RSIN[k,a]=R[k]*SIN[a];
RCOS[k,a]=R[k]*COS[a];
end
end
for a=1:L
for k=1:Nr
for l=1:Nr
x = RCOS[k,a] - RSIN[l,a] + w2;
y = RSIN[k,a] + RCOS[l,a] + h2;
x1 = floor(x);
x2 = ceil(x);
y1 = floor(y);
y2 = ceil(y);
if !(x1 <= 0 || x2 > w || y1 <= 0 || y2 > h)
# BILINEAR INTERPOLATION - STARTS
# x1 y1 - btm left
# x1 y2 - top left
# x2 y1 - btm rght
# x2 y2 - top rght
f11 = A[Int64(x1),Int64(y1)];
f21 = A[Int64(x1),Int64(y2)]; # that looks like, correct way
f12 = A[Int64(x2),Int64(y1)];
f22 = A[Int64(x2),Int64(y2)];
value = 0;
if x2==x1 && y2==y1
value = f11;
elseif x2==x1
value = 1/(y2-y1)*( f11*(y2-y) + f22*(y-y1) );
elseif y2==y1
value = 1/(x2-x1)*( f11*(x2-x) + f22*(x-x1) );
else
value = 1/(x2-x1)/(y2-y1)*( f11*(x2-x)*(y2-y) +
f21*(x-x1)*(y2-y) +
f12*(x2-x)*(y-y1) +
f22*(x-x1)*(y-y1) );
end
# BILINEAR INTERPOLATION - ENDS
sinogram[k,a] += value;
end
end
end
end
return sinogram;
end
#########################################
function iradon(A,angles) # input angles [rad]
angles = angles + pi/2;
N, nAngles = size(A);
I = zeros(N,N); # reconstruction
x = linspace(-0.5,0.5,N);
filter = abs(linspace(-1, 1, N));
# FT domain filtering
for t=1:length(angles)
fhat = fftshift(fft(slice(A,:,t)));
A[:,t] = real(ifft(ifftshift(fhat.*filter)));
end
XCOS = zeros(N,length(angles));
XSIN = zeros(N,length(angles));
for k=1:N
for a=1:length(angles)
XCOS[k,a]=x[k]*cos(angles[a]);
XSIN[k,a]=x[k]*sin(angles[a]);
end
end
for m=1:N
for n=1:N
for t=1:nAngles
r = XCOS[m,t] + XSIN[n,t];
index = Int64(round((r/sqrt(2.0) + 0.5)*N)); # sqrt(2) magnified
# index = Int64(round((r + 0.5)*N)); # to keep pixel resolution
if index>0 && index<=N
I[m,n] += A[index,t];
end
end
end
end
return I;
end
```

test program:

```
include("tomography.jl")
using Images, TestImages, ImageView
angles = (0:1:359)/360*2*pi; # [rad]
z = shepp_logan(1000; highContrast=true);
ImageView.imshow(z)
z = float(z);
tic();
sinogram = radon(z,angles);
display(toc());
v1 = minimum(sinogram);
v2 = maximum(sinogram);
todisplay = (sinogram-v1)/(v2-v1);
ImageView.imshow(todisplay);
tic();
reconstruction = iradon(sinogram,angles);
display(toc());
v1 = minimum(reconstruction);
v2 = maximum(reconstruction);
todisplay = (reconstruction-v1)/(v2-v1);
ImageView.imshow(todisplay);
```