FBP for 3D Euler geometry?

[robtovey]

Hi,

As the heading suggests, I am trying to perform a FBP on data that is computed on an Euler geometry (2d grid of angles, no translations) but there are a couple of points where the machinery breaks straight away.

There is also no example that I can see here which suggests that someone has probably known about this hole and not gotten around to it yet. If someone could give me a couple of pointers then I'd be happy to do the leg work on this, I can't imagine that it is a huge amount of work but if I'm mistaken then please correct me!

As far as I can tell, there are two conceptual changes. An Euler geometry can have a 2D grid of angles and does not have a rotation axis. Other than this I think the implementation is exactly the same just with a couple of minor tweaks which I'm hoping someone else will see clearer than me. Any advice on the 'proper' way to deal this case and which lines are causing the problems I've listed below would be great.

My understanding is that the important function is fbp_filter_op and the offending lines are:

• 367: an Euler geometry can have a 2D angular grid which doesn't have a length, can this be substituted by volume?
• 392: an Euler geometry doesn't have a rotation axis, can we just jump to axes=[1,2]?
• 420: if we don't have a rot_dir then abs_freq should be something like sqrt(x[1]**2+x[2]**2)?
• 450: generates an error I don't understand. 'operands could not be broadcast together with shapes (3,) (4,).' I guess this is also something to do with some Euler geometries having 2D angular domains? I think ran_shp needs to be length 4 too?
• 474: multiplying operators together gives an error that 'cell sides of domain and range differ significantly'. I don't really understand this one either but it is an issue with the zero-padded Fourier transform, the image and domain are no longer equal.

I have attached my partial solution in case that is of use, the lines I've changed are marked with # HERE flags. filtered_back_projection.py.txt Thanks again for your help.

[adler-j]

This is a very interesting issue and I really liked the research you did. I'll skip the questions (for now), not because they are bad but because I need more background to be able to answer properly.

To start doing this, we need to first figure out what we're supposed to do mathematically. The regular FBP algorithm is really only valid for parallel beam geometries and its extension to fan beam and especially cone beam is a bit complicated. However, extending it to e.g. a Euler geometry is even more complicated and would probably require a bit of effort. Do you have some reference to an article or book where this is done so we have something to follow? Perhaps an implementation in some other library?

[ozanoktem]

I think Euler angles Is just another way to parametrise a 3D parallell beam geometry commonly used in both electron tomography and single particle analysis.

[ozanoktem]

The note FBP_in_3D.pdf derives explicit expressions for the filter in FBP for inverting the ray transform restricted to a general parallel beam line complex in R^n. One may also consider consulting this paper that treats a more general setting.

[robtovey]

I don't have a reference for this but the way I see it, with 2 angles, you have one angle which corresponds to the tilt axis and another which corresponds to a tilt angle (alternatively the first angle selects a great circle and the second is the coordinate within it). We know the FBP for a simple single tilt axis so, in a noiseless case, if we then integrate these over all of our tilt axes we should just be summing up lots of copies of the exact reconstruction which will just be off by a constant factor, I think \pi.

There is a question of what happens when you have a sub-sampled sphere but I guess the linearity just lets you do the obvious.

Extending to 3 angles should similarly just be off by a constant, probably 2\pi^2.

Thanks for the references Ozan, I can't read the IP paper at the moment (despite the wonders of shiboleth...) but the attachment seems to agree that not much changes as we go up to 2 angles other than sorting out some constants depending on the spacing of the uniform grid (I think this is the role of \omega'(\tau)/\omega(\tau)).

Does that sound reasonable?

[robtovey]

Sorry, I'm going to disagree with myself now, I really shouldn't try to do maths past 10pm...

Thinking about it as a stack of 1D FBPs is wrong. If you re-parametrised the sphere then you would get different 1D filters for each projection which might still give you the same answer in an analytical noise-free example but not in general.

The important thing is that in 3D we naturally have a 2D sample on the sphere and the filter is simply the 2D ramp function. Hopefully the following is correct:

Also, if you have restricted domain:

Now we just follow the standard methods, as in Ozan's first reference, and realise the FBP can invert this filter in the data domain as well where the angular indicator function vanishes and this is the method used by the current code.

Broadly speaking, this is the method followed here on page 26/27 although the final exercise suggests c = 0.5 for all n which I guess is a typo...