Bug in a tutorial with a CT scan

[fedor-goncharov]

Hi guys,

I tried to use your lib for classical Radon transform in 2D along straight lines, so I reproduced the code given in the tutorial at PyLops CT-Tutorial

I believe that there is a problem with this tutorial - the image of the resulting sinogram is not how sinograms usually look like.

Your image

Usual sinogram

A signal in a sinogram (you tutorial calls it "Data") decays to zero on one axis (vertical axis in images above) but does not decay to zero along other axis (horizontal axis in images above) (PyLOps sinogram has clear zeros at left and right ends along horizontral axis) since it it corresponds to an angular variable (I don't want to go to math here - one who implemented Radon transform should understand this quickly).

P. S. Maybe inside integration a Jacobian normalization was forgotten or smth like that since decay to zero is still quite smooth.

[fedor-goncharov]

I checked - it is indeed a bug with a jacobian integration - library is ok :)
To fix it - you have to multiply every horizontal line by 1/sin(theta), i.e.,

correct_sinogram = np.divide(data, np.sin(theta), axis = 0)

I will do a pull-request on this correction if you don't mind. Until the pull is done, I don't close the issue myself.

[mrava87]

Hello @fedor-goncharov, thanks for pointing this out.

The Radon implementation used here is slightly different than what you would usually do for Radons in CT imaging (basically uses Cartesian coordinates instead of polar). In the example we use it in a bit of a special way… the point you make about the Jacobian makes sense to me, and indeed perhaps this is what is missing here.

Happy for you to go ahead and make a PR adding this scaling to the data. Note however 2 things:

• the data would need to scaled back before inversion of the operator should be multiplies by the scaling as otherwise we won’t get back the original model
• Given that you divide by sin(theta) this would lead to infinite scaling at 0 and 180 degrees, how you plan to handle that?

[fedor-goncharov]

Hello @mrava87,

Thanks for your reply. I agree, that the correct point of the tutorial is to use normal equations / split-bregman for matrices that are given by your library which are close to Radon Transforms in CT but not exactly the same. In this scope everything is fine with the tutorial - but it makes it further away from CT scenario - that's why I would like to make this small fix.

So I think 2 fixes are possible:

• short - just to put a phrase 'Note that in classical CT a jacobian normalization must be applied and maybe display an image with corrected sinogram - leave equations as they are.
• longer - define completely new operator (if possible just finding another parametrization of straight lines e.g. - with parametrization simply along the line but not using the angle variable) - in this case no division by zero problem at all

P. S. I actually wonder if in the definition of the curve function one really needs to reference 'ny' variable since it is not in the variable list.... (a bit not clean, wonder if it is possible to clean)

• Division by zero can be fixed just removing the end-points if we keep old parameters. Removing this little data will not affect reconstructions.

[mrava87]

Right, for the division by zero I was more worried about the scaling getting very large when approaching the end-points, of course the end-points themselves could be just removed :)

Regarding the short/long options, I would suggest going for the short and I explain why. The initial idea of this tutorial was to show how we could adapt the already available Radon operator to operate almost like the one commonly used for CT. However, we are not really interested to build a new operator for CT as it is much easier there to leverage existing, high-performance implementation in libraries like Astra-toolbox. I actually have had a branch for quite some time (see https://github.com/mrava87/pylops/tree/feature-astra) to bring in this and adapt the tutorial (no more strange decaying at edges... unfortunately I am currently struggling to make CI and Readthedocs work with Astra, this is the reason why this is not in yet.

[mrava87]

Hello @fedor-goncharov, I havent heard back from you. Are you still interested to contribute your example as discussed above?

[fedor-goncharov]

Hi, I just finished the corrections but struggle to follow the contributions guide - e.g., to run tests - 'setup.py' file is missing (whether I should do it at all since I implement no new class).

[mrava87]

Hello, good good.

Mmh the setup.py was removed two days ago as now PEP enforces a new way to build packages and so we needed to adapt to it for new versions of setuptools to work. I have just fixed the Makefile (https://github.com/PyLops/pylops/pull/534), I guess this is what you refer to as not being able to run tests?

However, as you say, if you do not implement new classed but just change a tutorial/example that step isn't really needed :) You can just make a PR and I can help finalizing it if anything else is needed