D parameter in fanbeam, ifanbeam. similar to radon and iradon!

Hello,

I am hoping someone can help me understand the "D" parameter that is used in fanbeam and ifanbeam. The help files say:

"D is the distance in pixels from the fan-beamvertex to the center of rotation. The center of rotation is the centerpixel of the image, defined as floor((size(I)+1)/2). D must be large enough to ensure that the fan-beam vertex is outside of the image at all rotation angles. As a guideline, try making D a few pixels larger than half the image diagonal dimension, calculated as follows sqrt(size(I,1)^2 + size(I,2)^2)"

So, I'm trying to understand where this value of D comes from, so I'm starting with the more simplified case which is a radon transform (parallel projections). Using : http://www.mathworks.com/help/images/examples/reconstructing-an-image-from-projection-data.html I am able to follow the example until we get to the value of N = 367. When I calculate the diagonal of the image (256x256) I get 363 (rounded up). Can someone help me understand where this 367 comes from? I would REALLY appreciate it!

You might be wondering why this matters to me. Basically, the image I'm using fanbeam on is full of values that are [per cm]. When I pass over the image with each beam, I need to know how long the actual beam is so that I can multiply the value by the integration length.

Thanks so much!! Kristin

[Sean de Wolski]

If Matt J doesn't answer this by EOD, I will do my best. But he's our SME for this the various CT projections here.

I'm really only versed in parallel beam..

[Kristin Busa]

Thanks, I appreciate it Sean!

[Kristin Busa]

I found some more info in the radon.m code for grandfathered syntax:

"R = RADON(I,THETA,N) returns a Radon transform with the projection computed at N points. R has N rows. If you do not specify N, the number of points the projection is computed at is: 2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3

This number is sufficient to compute the projection at unit intervals, even along the diagonal."

When using I = P from the phantom(256), the value 376 is returned. So, I suppose this is where it comes from, but I guess I'm still not sure why the +3 is required, perhaps just as a "buffer".

Sean/Matt, if you have any more info about why this process is the way it is, I'd appreciate. Otherwise, I've at least found where the N = 367 comes from. Thanks!

Answer 1

[Matt J]

What value of D do I use when calling the function ifanbeam? L1 or L2?

L1. By the way, if you run "doc fanbeam", you would see essentially the same diagram as the one you've drawn, but with L1 labelled as D. Unfortunately, the same diagram is not available in "doc ifanbeam"

What's also interesting is that because we're not measuring a solid object, the area that's shaded green is also captured by all of the fans and is not necessarily = 0. However, the output of ifanbeam is inherently a rectangular matrix (image).

Yes, but you could use a larger reconstruction grid (one that covers the entire green shaded region) if you wanted to reconstruct all of the non-zero material. The only reason, by the way, that you are allowed to reconstruct on a grid that is smaller than your green shaded region is because the particular algorithm used by ifanbeam is capable of reconstructing each pixel independently of all the others, and that flexibility comes at a price. In some reconstruction algorithms, however, all non-zero pixel values are obtained interdependently and so you would have to use a reconstruction grid large enough to cover all non-zero points.

BUT, the data needs to reflect that the length of the beam going through the rectangle. I.e., Lb ~= Ly. After I perform fanbeam on the data, I need to be able to multiply the integrated values by the length of the beam which passed through the matrix of the data. Does that make sense? Or does fanbeam.m already account for this difference in distance?

FANBEAM does take intersection length into account, but remember, it's doing a line integral along each ray through the object. So, it's not normally the total length of intersection with the image grid that is relevant, but rather the length of the intersection of each ray with each material. For example, if you look at your ray Ly, you can see that that ray intersects at least 3 different non-zero materials, and the intersection length with each is different. Approximately speaking, FANBEAM multiplies the length of each of those sub-segments of Ly with the image value of the material they each intersect. Those weighted lengths are summed to obtain the approximate line integral along Ly, and that's what FANBEAM returns as output. Of course, FANBEAM measures those lengths in pixels, as opposed to cm or mm, so there is a global scaling that you would have to perform if you want to work with physical units of length.

Answer 2

[Matt J]

So, I'm trying to understand where this value of D comes from

I may not be interpreting the question correctly. "Where does it come from" is a little vague to me. But basically, D is supposed to measure a physical dimension in the imaging system.

When you do fanbeam projection, it's usually because you're dealing with a system in which an x-ray source moves in a circular orbit around an object firing a fan-shaped beam of x-rays at it. The parameter D is supposed to represent the radius of this circular orbit measured in pixel lengths.

The requirement that "D must be large enough to ensure that the fan-beam vertex is outside of the image at all rotation angles" sounds a bit ill-phrased, because D is normally a fixed physical dimension of the system, not a parameter you adjust. It's the image grid size that you should be adjusting relative to D. Basically,the doc is saying that you have to use an image grid that is cropped small enough to fit completely within the circular orbit of the source. In practical CT systems, your object would normally fit inside this region so as not to get hit by the rotating gantry holding the source, so it's a very light restriction.

Basically, the image I'm using fanbeam on is full of values that are [per cm]. When I pass over the image with each beam, I need to know how long the actual beam is so that I can multiply the value by the integration length.

I don't think so. If you're projecting an attenuation image expressed in per cm, the length of a pixel in cm is what you would normally use as the integration length, not the length of the beam.

[Kristin Busa]

Matt, thanks for your response. I apologize that my question is unclear -- I suppose that reflects on my general confusion when using these functions!

My use of these functions is very similar to that of a medical x-ray/CT scanner. Instead of using x-ray projections to reconstruct a 2D cross-sectional image, I'm using laser absorption projections to reconstruct a 2D cross-sectional image of a gas' absorption.

I have two situations in which I work with the fanbeam/ifanbeam functions.

1. I have physically collected data using this type of fanbeam setup (illustrated in the image below). L1 was fixed by our equipment. The total fan angle theta_fan was chosen to encompass length, L2, when the sensor is perpendicular to the diagonal. This is the largest angle the fan will see as it rotates around the measurement space. What value of D do I use when calling the function ifanbeam? L1 or L2? Or something else all together? What's also interesting is that because we're not measuring a solid object, the area that's shaded green is also captured by all of the fans and is not necessarily = 0. However, the output of ifanbeam is inherently a rectangular matrix (image).

2. I have a matrix of a fixed number of pixels. This matrix is a computer-simulated solution (from a separate source) of what my solution "may" look like. The pixels in this matrix can be directly scaled to physical length. I want to perform a theoretical fanbeam measurement on this matrix, so I use fanbeam. This is nice because now I don't have any data where that green-shaded area is, BUT, the data needs to reflect that the length of the beam going through the rectangle. I.e., Lb ~= Ly. After I perform fanbeam on the data, I need to be able to multiply the integrated values by the length of the beam which passed through the matrix of the data. Does that make sense? Or does fanbeam.m already account for this difference in distance?

Thanks for your patience with me, I really do appreciate it.

Kristin

Answer 3

[Kristin Busa]

Matt,

Thanks so much for all of your responses. They have been really helpful as I learn to use these functions. I have one more question for you, then I'll leave you alone! :)

My question is regarding situation #2 I described above where I have my matrix of values that I wish to apply fanbeam to. So my scale from the matrix to real life is 1 pixel = .1 cm. This allows me to perform that global scaling you mentioned above. This situation is mimicking a real-life experiment, in which the distance from sensor to the center of the measurement was 44.2 cm (442 pixel lengths). We measured 21 rays with 1.3 degree separation, giving a total fan angle of 26 degrees. We rotated the fan in increments of 9 degrees. When I try to call fanbeam using these parameters:

fanbeam(data,442,'FanRotationIncrement',9,'FanSensorGeometry','arc', 'FanSensorSpacing',1.3)

I get the error "Subscripted assignment dimension mismatch." This also happens for larger values of D, for example D = 500. If I use a smaller D, for example, D = 420, it works fine. Any idea why this is the case? Also, specifying D leads to the program determining the # of rays and fans.. since I've already physically done a certain setup, is there a way to specify the number of rays and fans? (I've learned to guess/check w/ values of D to get my desired sinogram dimensions, I was just wondering if you had a more analytical approach).

Thanks again for your help. Kudos! Kristin

[Matt J]

Hi Kristin,

Before moving on to new questions, we should close this question first. Please accept-click my Answer above (since you seem to be satisfied with it) and then post your new question in a separate thread.

It would also be helpful to this board if you closed some of your other open questions, e.g., this one

http://www.mathworks.com/matlabcentral/answers/51546-using-ifanbeam-for-a-reconstruction-getting-false-negative-values

unless you plan to continue them.

Finally, when you post your new question, you should probably mention the dimensions of the "data" array so that it might be possible for us to run equivalent code. You should also clarify where your "sensor" is. Does it form an arc-shaped array of cells lying along the arc of the fanbeam or is it at the fan beam's vertex. Note that L1=D is supposed to be the distance to the vertex.