Asked by Bibek
on 27 Dec 2011

I have a CT image data(transverse slices),CT of 512*512*30. I know how to find CT # for any slice k corresponding to ith row and jth column, CT#=CT(i,j,k).I now trace a ray,starting from source point(0,0,100) then goes through CT image data and finally strikes to detector point(say 10,25,-50). Spatial coordinates for source and detector points are determined with isocenter taken at (0,0,0).Using linspace I can find spatial coordinates(i.e.,x,y,z)of points between source point and detector point. So, I have,say 105 intermediate points and I know their x,y,z coordinates.Say, 15 points lie outside of CT image data(i.e.in air) before passing through it and another 15 points also lie outside of CT image data(i.e.in air) after exiting through it. For remaining 75 points(assuming a ray goes through each slice), how can I assign CT # to each points(with their known spatial coordinates). Basically it means how can I get CT#=CT(x,y,z,15)for a point in 15th slice with its spatial coordinates x,y,z. I tried with, c = improfile(I,xi,yi,n)but got error message when tried to extend it, c = improfile(I,xi,yi,zi,n). Any help will be appreciated. Thanks!!

*No products are associated with this question.*

Answer by Image Analyst
on 27 Dec 2011

Why not just use trilinear interpolation? I don't think there's a built-in function for it, but it would be easy enough to program up yourself.

Opportunities for recent engineering grads.

## 6 Comments

## Naz (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/24805#comment_54793

I understood that you have some number of rays and you want to assign to each coordinate a number of the ray that passes through it. Is that right?

It depends on how you design your system. If your x-ray is a point source, than unlimited number of rays pass through each coordinate. If you have an array source of parallel rays, then you can assume some width of the ray and consequently each ray will pass through a limited set of coordinates around the line of propagation (no scattering).

## Naz (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/24805#comment_54795

In other words you can imagine each ray as a bar rather than a line.

## Bibek (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/24805#comment_54828

I think I am trying to do something different. I traced a ray through CT slice data. I also know 3 dimensional spatial coordinates of points through which the ray traverses. I now need to find out CT number corresponding to these points.For example if the ray passes through a point P with spatial coordinates (-10,15,-50)then I need to figure out what is CT # corresponding to P.

## Naz (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/24805#comment_54844

And what is CT#?

## Image Analyst (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/24805#comment_54863

Those are the digital numbers the Computed Tomography unit gives, called Hounsfield Units: http://en.wikipedia.org/wiki/Hounsfield_scale.

Bibek, how did my suggestion work out for you?

## Bibek (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/24805#comment_55019

In matlab help I found the code, VI = interp3(X,Y,Z,V,XI,YI,ZI). So, in my case V=CT, it is of size 512*512*30,X=x_coordinates, Y=y_coordinates and Z=z_coordinates. All these X,Y,Z are matrices I created myself and are also of size 512*512*30. Since I am trying to find CT # corresponding to a point with x,y,z coordinates(1.5,25,15), I kept XI=1.5,YI=25,ZI=15. But I got a error message:

vi=interp3(x_Coordinates,y_Coordinates,z_Coordinates,CT,1.5,25,15);

??? Error using ==> interp3 at 138

X, Y and Z must be matrices produced by MESHGRID. Use

TriScatteredInterp instead

of INTERP3 for scattered data.

In my case I create x,y,z coordinates as

X9=1:1:512;

x_Coordinates=repmat(X9,[512,1,75]);% first row corresponds to 712 z-coordinates(all zeros) for 1st image and so on

Y9=1:1:512;

y_Coordinates=repmat(Y9,[512,1,75]);

Z9=0:2.5:185;

for i=1:75;

z_Coordinates(:,:,i)=repmat(Z9(1,i),[512,512,1]);

end

Since I don't know anything about meshgrid, I am not sure I will be able to create x,y,z coordinates like above using meshgrid.

When I tried with TriScatteredInterp I got this error message:

vi=TriScatteredInterp(x_Coordinates,y_Coordinates,z_Coordinates,CT,1.5,25,15);

??? Error using ==> TriScatteredInterp

Incorrect number of arguments for TriScatteredInterp.

I don't know where I did a mistake. Probably I was wrong in defining X,Y,Z and XI,YI,ZI. Any help will be appreciated.