Thread Subject: uniaxial rotation of 3D array

Hi Group,

I encountered a problem that I could not figure out a way to solve. It
appears to be simple (it might just be as being smart has nothing to
do with me ... ). But the detailed execution keeps evading my fingers
or brain... I'd love if you could shed me some lights. Thanks.

Suppose that I have a 3D object, which can be written as I = I (x, y,
z). We can envision this as a real space object. I need to obtain its
2D projection along the z axis,

I_prime (x, y) = Int_z (I (x, y, z)).

I also would like to have the projection of the object after it is
rotated along the z axis with a specified angle \phi, and obtain a
similar projection image after the rotation. I am, however, confused
about how to do this rotation in the original Euclid's coordinates.
Namely, how to execute R_z (\phi) I (x, y, z), where R_z(\phi) is an
operation on 3D array I(x, y, z).

I hope I have made myself clear. Please do not hesitate if you have
any question. Thanks!

Regards,
Frank

On Sep 23, 4:38 pm, Frank Chang <etaght...@gmail.com> wrote:
> Hi Group,
>
> I encountered a problem that I could not figure out a way to solve. It
> appears to be simple (it might just be as being smart has nothing to
> do with me ... ). But the detailed execution keeps evading my fingers
> or brain... I'd love if you could shed me some lights. Thanks.
>
> Suppose that I have a 3D object, which can be written as  I = I (x, y,
> z).  We can envision this as a real space object. I need to obtain its
> 2D projection along the z axis,
>
> I_prime (x, y) = Int_z (I (x, y, z)).
>
> I also would like to have the projection of the object after it is
> rotated along the z axis with a specified angle \phi, and obtain a
> similar projection image after the rotation. I am, however, confused
> about how to do this rotation in the original Euclid's coordinates.
> Namely, how to execute R_z (\phi) I (x, y, z), where R_z(\phi) is an
> operation on 3D array I(x, y, z).
>
> I hope I have made myself clear. Please do not hesitate if you have
> any question. Thanks!
>
> Regards,
> Frank

------------------------------------------------------------------------------------------------------
Frank:
Have you looked at the Radon Transform (the basis of computed
tomography)? You're in luck because this is exactly what it does. It
appears to be built in to the base MATLAB product.
Good luck,
ImageAnalyst

Hi ImageAnalyst,

Thanks for your reply. I have never heard of Radon Transform. So this
is completely new to me.

I checked the help file of radon(), it appears that it does a 1D
projection of 2D image. A google for "3D radon matlab" does not give
obvious hints that the original radon function is applicable to 3D
projection. It would be great if you could lead me a bit further.
Anyway, thank you for pointing this out!

Thanks,
Frank

On Sep 23, 5:38 pm, ImageAnalyst <imageanal...@mailinator.com> wrote:
> On Sep 23, 4:38 pm, Frank Chang <etaght...@gmail.com> wrote:
>
>
>
> > Hi Group,
>
> > I encountered a problem that I could not figure out a way to solve. It
> > appears to be simple (it might just be as being smart has nothing to
> > do with me ... ). But the detailed execution keeps evading my fingers
> > or brain... I'd love if you could shed me some lights. Thanks.
>
> > Suppose that I have a 3D object, which can be written as  I = I (x, y,
> > z).  We can envision this as a real space object. I need to obtain its
> > 2D projection along the z axis,
>
> > I_prime (x, y) = Int_z (I (x, y, z)).
>
> > I also would like to have the projection of the object after it is
> > rotated along the z axis with a specified angle \phi, and obtain a
> > similar projection image after the rotation. I am, however, confused
> > about how to do this rotation in the original Euclid's coordinates.
> > Namely, how to execute R_z (\phi) I (x, y, z), where R_z(\phi) is an
> > operation on 3D array I(x, y, z).
>
> > I hope I have made myself clear. Please do not hesitate if you have
> > any question. Thanks!
>
> > Regards,
> > Frank
>
> ------------------------------------------------------------------------------------------------------
> Frank:
> Have you looked at the Radon Transform (the basis of computed
> tomography)?  You're in luck because this is exactly what it does.  It
> appears to be built in to the base MATLAB product.
> Good luck,
> ImageAnalyst

On Sep 23, 9:43 pm, Frank Chang <etaght...@gmail.com> wrote:
> Hi ImageAnalyst,
>
> Thanks for your reply. I have never heard of Radon Transform. So this
> is completely new to me.
>
> I checked the help file of radon(), it appears that it does a 1D
> projection of 2D image. A google for "3D radon matlab" does not give
> obvious hints that the original radon function is applicable to 3D
> projection. It would be great if you could lead me a bit further.
> Anyway, thank you for pointing this out!
>
> Thanks,
> Frank
>
> On Sep 23, 5:38 pm, ImageAnalyst <imageanal...@mailinator.com> wrote:
>
>
>
> > On Sep 23, 4:38 pm, Frank Chang <etaght...@gmail.com> wrote:
>
> > > Hi Group,
>
> > > I encountered a problem that I could not figure out a way to solve. It
> > > appears to be simple (it might just be as being smart has nothing to
> > > do with me ... ). But the detailed execution keeps evading my fingers
> > > or brain... I'd love if you could shed me some lights. Thanks.
>
> > > Suppose that I have a 3D object, which can be written as  I = I (x, y,
> > > z).  We can envision this as a real space object. I need to obtain its
> > > 2D projection along the z axis,
>
> > > I_prime (x, y) = Int_z (I (x, y, z)).
>
> > > I also would like to have the projection of the object after it is
> > > rotated along the z axis with a specified angle \phi, and obtain a
> > > similar projection image after the rotation. I am, however, confused
> > > about how to do this rotation in the original Euclid's coordinates.
> > > Namely, how to execute R_z (\phi) I (x, y, z), where R_z(\phi) is an
> > > operation on 3D array I(x, y, z).
>
> > > I hope I have made myself clear. Please do not hesitate if you have
> > > any question. Thanks!
>
> > > Regards,
> > > Frank
>
> > ---------------------------------------------------------------------------­---------------------------
> > Frank:
> > Have you looked at the Radon Transform (the basis of computed
> > tomography)?  You're in luck because this is exactly what it does.  It
> > appears to be built in to the base MATLAB product.
> > Good luck,
> > ImageAnalyst- Hide quoted text -
>
> - Show quoted text -

------------------------------------------------------------------------------------
Frank:
Yes, it appears that the MATLAB version only handles projections along
2D images, not 3D images like you want. So you'll have to look
elsewhere for 3D Radon transform code. You might search the CT
literature or perhaps Harrison Barret's or James Green's books on
medical image processing math.
Good luck,
ImageAnalyst

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <ae253e3e-9127-4490-a08c-90196c4a2515@l34g2000vba.googlegroups.com>...
> On Sep 23, 9:43?pm, Frank Chang <etaght...@gmail.com> wrote:
> > Hi ImageAnalyst,
> >
> > Thanks for your reply. I have never heard of Radon Transform. So this
> > is completely new to me.
> >
> > I checked the help file of radon(), it appears that it does a 1D
> > projection of 2D image. A google for "3D radon matlab" does not give
> > obvious hints that the original radon function is applicable to 3D
> > projection. It would be great if you could lead me a bit further.
> > Anyway, thank you for pointing this out!
> >
> > Thanks,
> > Frank
> >
> > On Sep 23, 5:38?pm, ImageAnalyst <imageanal...@mailinator.com> wrote:
> >
> >
> >
> > > On Sep 23, 4:38?pm, Frank Chang <etaght...@gmail.com> wrote:
> >
> > > > Hi Group,
> >
> > > > I encountered a problem that I could not figure out a way to solve. It
> > > > appears to be simple (it might just be as being smart has nothing to
> > > > do with me ... ). But the detailed execution keeps evading my fingers
> > > > or brain... I'd love if you could shed me some lights. Thanks.
> >
> > > > Suppose that I have a 3D object, which can be written as ?I = I (x, y,
> > > > z). ?We can envision this as a real space object. I need to obtain its
> > > > 2D projection along the z axis,
> >
> > > > I_prime (x, y) = Int_z (I (x, y, z)).
> >
> > > > I also would like to have the projection of the object after it is
> > > > rotated along the z axis with a specified angle \phi, and obtain a
> > > > similar projection image after the rotation. I am, however, confused
> > > > about how to do this rotation in the original Euclid's coordinates.
> > > > Namely, how to execute R_z (\phi) I (x, y, z), where R_z(\phi) is an
> > > > operation on 3D array I(x, y, z).
> >
> > > > I hope I have made myself clear. Please do not hesitate if you have
> > > > any question. Thanks!
> >
> > > > Regards,
> > > > Frank
> >
> > > ---------------------------------------------------------------------------?---------------------------
> > > Frank:
> > > Have you looked at the Radon Transform (the basis of computed
> > > tomography)? ?You're in luck because this is exactly what it does. ?It
> > > appears to be built in to the base MATLAB product.
> > > Good luck,
> > > ImageAnalyst- Hide quoted text -
> >
> > - Show quoted text -
>
> ------------------------------------------------------------------------------------
> Frank:
> Yes, it appears that the MATLAB version only handles projections along
> 2D images, not 3D images like you want. So you'll have to look
> elsewhere for 3D Radon transform code. You might search the CT
> literature or perhaps Harrison Barret's or James Green's books on
> medical image processing math.
> Good luck,
> ImageAnalyst

It seems to me that, if you are dealing with parallel (as opposed to fan or cone beam) projections, you should be able to use the "radon" command on each 2-D slice ("image") along the z axis, i.e. the process is repeated at all discrete heights (z) along the object and the 2-D projections are obtained by stacking the 1-D projections one on top of another.

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <ae253e3e-9127-4490-a08c-90196c4a2515@l34g2000vba.googlegroups.com>...
> On Sep 23, 9:43?pm, Frank Chang <etaght...@gmail.com> wrote:
> > Hi ImageAnalyst,
> >
> > Thanks for your reply. I have never heard of Radon Transform. So this
> > is completely new to me.
> >
> > I checked the help file of radon(), it appears that it does a 1D
> > projection of 2D image. A google for "3D radon matlab" does not give
> > obvious hints that the original radon function is applicable to 3D
> > projection. It would be great if you could lead me a bit further.
> > Anyway, thank you for pointing this out!
> >
> > Thanks,
> > Frank
> >
> > On Sep 23, 5:38?pm, ImageAnalyst <imageanal...@mailinator.com> wrote:
> >
> >
> >
> > > On Sep 23, 4:38?pm, Frank Chang <etaght...@gmail.com> wrote:
> >
> > > > Hi Group,
> >
> > > > I encountered a problem that I could not figure out a way to solve. It
> > > > appears to be simple (it might just be as being smart has nothing to
> > > > do with me ... ). But the detailed execution keeps evading my fingers
> > > > or brain... I'd love if you could shed me some lights. Thanks.
> >
> > > > Suppose that I have a 3D object, which can be written as ?I = I (x, y,
> > > > z). ?We can envision this as a real space object. I need to obtain its
> > > > 2D projection along the z axis,
> >
> > > > I_prime (x, y) = Int_z (I (x, y, z)).
> >
> > > > I also would like to have the projection of the object after it is
> > > > rotated along the z axis with a specified angle \phi, and obtain a
> > > > similar projection image after the rotation. I am, however, confused
> > > > about how to do this rotation in the original Euclid's coordinates.
> > > > Namely, how to execute R_z (\phi) I (x, y, z), where R_z(\phi) is an
> > > > operation on 3D array I(x, y, z).
> >
> > > > I hope I have made myself clear. Please do not hesitate if you have
> > > > any question. Thanks!
> >
> > > > Regards,
> > > > Frank
> >
> > > ---------------------------------------------------------------------------?---------------------------
> > > Frank:
> > > Have you looked at the Radon Transform (the basis of computed
> > > tomography)? ?You're in luck because this is exactly what it does. ?It
> > > appears to be built in to the base MATLAB product.
> > > Good luck,
> > > ImageAnalyst- Hide quoted text -
> >
> > - Show quoted text -
>
> ------------------------------------------------------------------------------------
> Frank:
> Yes, it appears that the MATLAB version only handles projections along
> 2D images, not 3D images like you want. So you'll have to look
> elsewhere for 3D Radon transform code. You might search the CT
> literature or perhaps Harrison Barret's or James Green's books on
> medical image processing math.
> Good luck,
> ImageAnalyst

It seems to me that, if you are dealing with parallel (as opposed to fan or cone beam) projections, you should be able to use the "radon" command on each 2-D slice ("image") along the z axis, i.e. the process is repeated at all discrete heights (z) along the object and the 2-D projections are obtained by stacking the 1-D projections one on top of another.

Frank Chang <etaghtron@gmail.com> wrote in message <e32f3b1f-551c-427a-8f7d-172ce4be1f50@2g2000prl.googlegroups.com>...
> Hi Group,
>
> I encountered a problem that I could not figure out a way to solve. It
> appears to be simple (it might just be as being smart has nothing to
> do with me ... ). But the detailed execution keeps evading my fingers
> or brain... I'd love if you could shed me some lights. Thanks.
>
> Suppose that I have a 3D object, which can be written as I = I (x, y,
> z). We can envision this as a real space object. I need to obtain its
> 2D projection along the z axis,
>
> I_prime (x, y) = Int_z (I (x, y, z)).
>
> I also would like to have the projection of the object after it is
> rotated along the z axis with a specified angle \phi, and obtain a
> similar projection image after the rotation. I am, however, confused
> about how to do this rotation in the original Euclid's coordinates.
> Namely, how to execute R_z (\phi) I (x, y, z), where R_z(\phi) is an
> operation on 3D array I(x, y, z).
>

Rotate an object involve interpolation of the intensity on a rotated grid. So the first step is to compute a rotated grid. A nice trick is to use function VIEW to compute the transformation matrix from the two provided rotation angles, then multiply this matrix on the original grid coordinates to get the rotated one. You might want to place the origin at the rotation center before. Next use INTERP3 to interpolate the intensity on the new grid. Finally use what ever discrete integration such as TRAPZ or SUM to integrate the intensity along an axis.

Bruno

>
> Rotate an object involve interpolation of the intensity on a rotated grid. So the first step is to compute a rotated grid. A nice trick is to use function VIEW to compute the transformation matrix from the two provided rotation angles, then multiply this matrix on the original grid coordinates to get the rotated one. You might want to place the origin at the rotation center before. Next use INTERP3 to interpolate the intensity on the new grid. Finally use what ever discrete integration such as TRAPZ or SUM to integrate the intensity along an axis.
>
> Bruno

You are making the problem more complicated than it actually is. The "radon" command does the interpolation for you.
Here is some quick code I wrote to test this on a sphere:
[x, y, z] = meshgrid (-50:49, -50:49, -50:49);
C = sqrt(x.^2+y.^2+z.^2);
I = (C<40); % Sphere of radius 40 in matrix of size 100x100x100, center at x=y=z=50

theta = 0:3:180; % Projections every 3?

P = radon(I(:, :, 1),theta); % Just to get the dimensions of P

% P_all will hold all the slice projections
% e.g. P_all(:, :, 1) holds all the 1-D projections for z=1, column by
% column, i.e. P_all(:, 1, 1) e.g. holds the 1-D projection for theta=0?
% at z=1
P_all = rand(size(P, 1), size(P, 2), size(I, 3)); % create P_all w/ rand or zeros e.g.

for z = 1:100
    P_all(:, :, z) = radon(I(:, :, z),theta); % Fill P_all with all projections of slice z
end

% Perform backprojections slice by slice:
for z = 1:100
    I_recon(:, :, z) = iradon(P_all(:, :, z),theta,'linear','Hamming');
end

% Plot a few slices at different heights (z):
set(gcf, 'Position', [10 10 1660 960]);set(gcf, 'Position', [10 10 1660 960]);
subplot(2,3,1), imshow(I(:, :, 15),[]), title('Original at z=15')
subplot(2,3,2), imshow(I(:, :, 50),[]), title('Original at z=50')
subplot(2,3,3), imshow(I(:, :, 70),[]), title('Original at z=70')
subplot(2,3,4), imshow(I_recon(:, :, 15),[]), title('Back-proj at z=15')
subplot(2,3,5), imshow(I_recon(:, :, 50),[]), title('Back-proj at z=50')
subplot(2,3,6), imshow(I_recon(:, :, 70),[]), title('Back-proj at z=70')