0
$\begingroup$

I'm quite new to MATLAB and I'm trying to obtain an implementation of the computed tomography without using any built-in functions such as radon() and iradon(). I guess I understood what Radon Transform does: it simply adds up projections taken at some angles (in matrix columns in MATLAB) and produces a sinogram as a result. But I am at a loss at this point. My attempt never produces anything that is similar to a sinogram. What am I doing wrong? I cannot get ceil(128*sqrt2) as the maximum value for the detector size as well (of one unit each on the projection axis).

Side note: Well, I admit that the code is lousy (my image, too), but I tried my best in the last week to decipher MATLAB. Thank you very much in advance for any kind of help and suggestions.

A=zeros(128,128);
for k=1:128
    for n=1:128
        if k>70&&k<100&&n>70&&n<100
            A(k,n)=1;
            else if k>20&&k<35&&n>30&&n<55
            A(k,n)=0.7;
                else if k*k+n*n>100*40&&k*k+n*n<100*50
                    A(k,n)=0.8;
                    else if (k-50)^2+(n-50)^2<=400
                        A(k,n)=0.5;
                        else if (k-100)^2+(n-20)^2>=100&&(k-100)^2+(n-20)^2<=400
                                A(k,n)=0.8;
                            end
                    end
                end        
            end   
        end    
    end
end
imshow(A)
M=size(A);
[ilength, iwidth]=size(A);
idiag=sqrt(128^2+128^2);
lengthpad=ceil(idiag-ilength)+2;
widthpad=ceil(idiag-iwidth)+2;
pad=zeros(ilength+lengthpad,iwidth+widthpad);
pad(ceil(lengthpad/2):(ceil(lengthpad/2)+ilength-1),ceil(widthpad/2):(ceil(widthpad/2)+iwidth-1)) = A;
figure, imshow(pad)
N=100;

% angles through which we will perform the projection
for j=0:N-1
            thetaj=j*(pi/N);
            costhetaj=cos(thetaj);
            sinthetaj=sin(thetaj);


    %finding the midpoints of the pixels on the image 
    for iy=1:size(pad)
        for i=1:size(pad)

            y_img=iy-pad/2-0.5;
            x_img=i-pad/2-0.5;
            img=sum(costhetaj*x_img*pad(i,iy)+sinthetaj*y_img*pad(i,iy),1);
            last_img(i,iy)=sum(img);

        end
        end
end
    figure, imshow(last_img)
    figure, plot(last_img)

Edit: I added zero padding to the code, but I still don't know why I get half of the image as output.padded original image output

$\endgroup$
  • $\begingroup$ Can I please ask if this was resolved? $\endgroup$ – A_A Apr 30 '18 at 8:16
  • $\begingroup$ Hi, thank you very much for your help again! I managed to do the padding and rotate the image without any problem, and I obtained the sinogram with the help of the code you provided. I guess thank you messages are not appreciated here, this is why I never replied, sorry. $\endgroup$ – takashi z. May 7 '18 at 19:16
  • $\begingroup$ Thanks for letting me know and I am glad to hear you found the answer useful. You can upvote and / or accept it from the controls on the left of the answer post. This will also stop it from circulating in the board as "unanswered". All the best. $\endgroup$ – A_A May 7 '18 at 21:05
0
$\begingroup$

There are two points that stand out from the provided code:

  1. The determination of detector_size

    • You seem to be trying to track the "widest line of sums" (?). But we know that this is the diagonal and it occurs at 45 degrees in square images. Therefore, we can "pre-allocate" our detector to have that maximum size anyway.
  2. The rotation of the image

    • The code seems to be along the right track, taking the "image grid coordinates" (i,iy) and rotating them to (x_img,y_img) but:
      1. There is no centre of rotation in the calculation (?). This usually is taken to be at the centre of the image.
      2. While some rotation seems to be taking place, the data is still read from the original locations of A (involving some scaling too) (?).

There are two approaches that you can take here:

  1. Rotate the image
  2. Rotate the "detector" (but do it properly)

The first approach is dead simple but a bit wasteful in terms of memory and computation. It goes a little something like this:

  1. Given a $M\times N$ image ($I$), construct a $D \times D$ image ($A$) where $D = \sqrt{M^2 + N^2}$

  2. Now, centre $I$ into $A$. This basically translates $\left(0,0\right)$ of $I$ to $\left(\frac{D-M}{2},\frac{D-N}{2}\right)$ in $A$.

  3. For each angle $\theta \in \Theta$, rotate the image (using imrotate) to produce $A_r$. Now, note here that depending on your frame of reference, you might have to transform $\theta$. That is because, if you sum the columns of the rotated image, then your detector is positioned at $+ \frac{\pi}{2}$. But that is just a simple linear transformation to align your results with convention. Also important to note here that it is worth rotating $A$ afresh each time rather than keep rotating $A_r$ by increments, which would progressively distort the image due to successive applications of interpolation.

  4. Once you have $A_r$, produce the sum of the "detector" by adding up its columns.

And that is it.

It is dead simple because you don't have to handle the rotation of the image which, by itself, is not difficult but it involves interpolation which can get messy. The rotated image grid will never land on integer pixel locations and so you are going to have to apply some form of interpolation that goes beyond nearest neighbour. For linear interpolation (for example), this would involve determining the patch of the 4 nearest pixels to the target location, then obtaining the distances to those pixels and finally creating the weighted sum.

It is "wasteful" in CPU because $A$ is bigger than $I$ and all pixels are rotated even if their contribution to a sum is in fact zero.

It is "wasteful" in memory because $A$ is bigger than $I$ and for large differences of $M,N$, large square images are produced.

The immediate improvement you can apply to this is to rotate the "scanning grid" (as you are trying to do now) but when you obtain a slice off of the image to centre it in the "detector slice" which would have to be $D$ pixels wide anyway. In this case of course, you still have to handle interpolation. But in this approach, any coordinates landing outside of the original image are simply skipped.

Hope this helps

EDIT:

The rotation is still not performed correctly:

% Load the image
Q = imread("testImage.jpg");
% Resize the image, in case it is too big
Q = imresize(Q,0.05);
Q_size = size(Q);
% Find its diagonal
D = ceil(sqrt(sum(Q_size.^2)));
% Create a square image of size equal to the diagonal
I = zeros(D,D);
I_size = size(I);
% Center the image in the square
offset = ceil((D-Q_size)./2);
% Copy the original image across
I(offset(1):offset(1)+Q_size(1)-1,offset(2):offset(2)+Q_size(2)-1) = Q;

I contains the padded image.

Using imrotate:

% Theta angles to obtain sums at
theta = 0:11.25:180;
% Initialise the output. One row for each angle, one column for every summing "sensor".
R1 = zeros(length(theta), D);
for k=1:length(theta)
    R1(k,:)=sum(imrotate(I,theta(k),'nearest','crop'));    
end;

figure;imshow(R1./max(max(R1)))

The above finishes very quickly. Notice nearest as interpolation in imrotate and crop so that imrotate does not add yet more padding.

Not using imrotate:

% Rotating the detector
theta_rad = (theta./180).*pi; %This simply converts theta (from above) in radians
% Initialise the output. One row for each angle, one column for every summing "sensor".
R2 = zeros(length(theta_rad), D);
% THE CENTRE OF THE IMAGE
I_centre = round(I_size./2);
% FOR EVERY ANGLE
for k=1:length(theta_rad)
    % m,n scan the image left to right top to bottom...BUT...
    for m=0:(I_size(2)-1)
        for n=0:(I_size(1)-1)
            % THEY GET REMAPPED TO U,V
            % Remap m,n to their rotated versions.            
            % I do the rotation in polar coordinates here
            rho = sqrt(sum((I_centre-[n,m]).^2));
            phi = atan2(n - I_centre(1), m - I_centre(2));
            % This point now maps to
            u = round(I_centre(2) + rho .* cos(phi-theta_rad(k)));
            v = round(I_centre(1) - rho .* sin(phi-theta_rad(k)));
            % Check that the point is inside the image
            if (u>0 && u<I_size(1) && v>0 && v<I_size(2))
                % NOTICE THAT I IS READ FROM THE REMAPPED LOCATIONS
                R2(k,m+1) = R2(k,m+1) + I(v,u);
            end;
        end;
    end;
end;

figure;imshow(R2./max(max(R2)))
$\endgroup$
  • $\begingroup$ Thank you for the answer, A_A, it really helped. I prefer not to use imrotate, either, but with your explanation, it looks easy. Now then, is what I need to do to pad the image with zeros to prevent data loss and determine a centre of rotation (I thought it was taken as the centre automatically)? And I'm afraid I still don't get why I get half of the image as output. Is it because the rotation uses the original image each time? As for the detector_size, I just tried to calculate the number of detectors needed for each angle, but it doesn't work. It does not even show the correct value for max. $\endgroup$ – takashi z. Mar 29 '18 at 10:51
  • $\begingroup$ @takashiz. please see update. $\endgroup$ – A_A Mar 29 '18 at 14:12

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.