RI-imaging / radontea

python library for tomographic image reconstruction
https://radontea.readthedocs.io
BSD 3-Clause "New" or "Revised" License
14 stars 7 forks source link

Bug in `utils.compute_angle_weights_1d` #7

Closed paloha closed 2 years ago

paloha commented 2 years ago

I believe this function should return a flat line in case all the angles in theta are equidistant.

This is true in this example:

theta = np.arange(0., 360., 1.)
>> [0.  1.  2.  ...  357.  358.  359.]
weights = utils.compute_angle_weights_1d(np.deg2rad())
>> [1.  1.  1. ... 1.  1.  1.]

But there seems to be a problem in the following cases:

1) Not a full circle

But it is not the case on the edge angles (first and last in the vector) if we do not do a full circle:

theta = np.arange(0., 136., 1.)  
>> [0.  1.  2. ...  133.  134.  135.]
weights = utils.compute_angle_weights_1d(np.deg2rad())
>> [17.37777778  0.75555556  0.75555556  ...  0.75555556  0.75555556 17.37777778]

# theta = np.arange(-90., 46., 1.)  # This theta vector produces the same result

image

2) Edges lying on the same line

And then there is another problem, if the edge angles lie on the same line through the origin of the unit circle:

theta = np.arange(-90., 91., 1.)
>> [-90.  -89.  -88.  ...  88.  89.  90.]
weights = utils.compute_angle_weights_1d(np.deg2rad())
>> [0.50277778 1.00555556 1.00555556  ...  1.00555556 1.00555556  0.50277778]

image

# This theta vector produces the same result as above
theta = np.arange(-80., 101., 1.) 

image

3) Edges are the same angle

Also, when the edge angles are the same, we get a 0 weight on the 90°, which does not make sense IMO:

theta = np.arange(-90., 271., 1.)
>>
weights = utils.compute_angle_weights_1d(np.deg2rad())
>> [-90.  -89.  -88.  ...  89.  90.  91.  ...  268.  269.  270.]
>> [1.00277778  1.00277778  1.00277778  ...  1.00277778  0.  1.00277778  ...  1.00277778  1.00277778  1.00277778]

image

Summary

Either this is not a correct implementation, or I do not understand the intentions with this function. When I try to reconstruct a tomogram from the projections recorded at the usual tilt angles (-60°, 60°), the two edge tilts get extremely big weights and the reconstruction looks terrible. Thanks for looking at this.

paulmueller commented 2 years ago

I will try to explain what is happening:

  1. It is assumed that a full 180° must be covered for a proper reconstruction. The idea is to give each angle a weight that is proportional to the distance to its neighboring angles. Obviously, this idea does not work if the angles are strongly unevenly distributed (e.g. dense coverage from 0° to 90°, nothing between 90° and 180°). But it works for the example cases in the docs. This basically answers the behavior in case (1).

  2. For case (2), this also makes perfect sense. If there are two angles at the same location, then the corresponding recordings only get half of the weight.

  3. Case (3) might actually be a bug. The coverage of the weights is larger than 180° and this case is probably not treated well.

For your -60°/+60° case, you should probably disable weighting the angles (if they are evenly distributed). Otherwise, I would be open to introducing another keyword argument that defines the desired spread of the angles...

Cheers!

paloha commented 2 years ago

@paulmueller thank you for looking into this.

Yes, I disabled weighting. I just wanted to report the behavior I found odd when I was testing this function out.

In the point 1., I understand that 180° must be covered for a proper reconstruction, unfortunately, this is not what you can measure in a real microscope. Case 2. makes sense, OK. Case 3. looks like a bug indeed. Might be just solved by raising an error, because I also can not imagine a case where somebody would actually have coverage >180°. Cheers

paulmueller commented 2 years ago

Maybe the docs could also be clearer about point 1.