8000 iradon doesn't quantitatively agree with MATLAB · Issue #3742 · scikit-image/scikit-image · GitHub
[go: up one dir, main page]
Skip to content

iradon doesn't quantitatively agree with MATLAB #3742

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Open
philipstarkey opened this issue Feb 11, 2019 · 8 comments
Open

iradon doesn't quantitatively agree with MATLAB #3742

philipstarkey opened this issue Feb 11, 2019 · 8 comments
Labels
⚖️ Needs decision 🐛 Bug

Comments

@philipstarkey
Copy link

Description

I've looked at PR #3067 and PR #3099 (addressing #2853) and PR #683 and while the first two do fix some issues with the iradon function and the last improves feature parity between the skimage iradon and the MATLAB iradon, the output between the two (at least for my dataset) does not quantitatively agree.

The error appears to be cause by 3 off-by-one errors that I believe probably crept in when the code was originally ported from MATLAB code. I think whoever did this probably missed the fact that the A:B style syntax in MATLAB is used both as slicing (in which case you need to reduce A by 1 in Python to handle the 0 indexing of Python) and for defining a range (in which case you want to increment B by 1 in Python to account for the fact that ranges in MATLAB include the last integer whereas in Python they do not).

The fix (also requires PR #3067) is to the following two lines:
[X, Y] = np.mgrid[0:output_size, 0:output_size]
and
x = np.arange(radon_filtered.shape[0]) - mid_index
which must become
[X, Y] = np.mgrid[1:output_size+1, 1:output_size+1]
and
x = np.arange(1,radon_filtered.shape[0]+1) - mid_index

This increments the start of the ranges by 1 to undo the incorrect translation (so now the indices would match the original MATLAB code) and then increments the end point of the ranges by 1 to apply the correct translation between MATLAB and Python code. This change (merged with PR #3067) then gives me quantitatively the same output as MATLAB.

Caveat: I can't provide the data I'm working on as it's not data I'm allowed to make public. However I can tell you that radon_image is a 26x30 and I'm calling iradon with theta = np.arange(0,180,6), output_size=20, filter='ramp', interpolation='linear', circle=False

Version information

# your output here
>>> import sys; print(sys.version)
3.6.4 |Anaconda, Inc.| (default, Jan 16 2018, 10:22:32) [MSC v.1900 64 bit (AMD64)]
>>> import platform; print(platform.platform())
Windows-10-10.0.17134-SP0
>>> import skimage; print("scikit-image version: {}".format(skimage.__version__))
scikit-image version: 0.14.2
>>> import numpy; print("numpy version: {}".format(numpy.__version__))
numpy version: 1.15.0
@hmaarrfk
Copy link
Member

Hey. Thanks for the detailed report.

I’m trying to get one of the remaining PRs through that you listed with regards to iradon.

If you can answer that last question I had in the review that would help it move forward.

It seems that it still wouldn’t resolve the issues you are pointing out. to ensure the issue is resolved to your satisfaction, we will need some kind of test. It doesn’t have to be with real data. Feel free to create a phantom datasets or use numpy random functionality with a predetermined seed.

If checking against the full matlab array is too complicated, maybe check against a few pixels?

@philipstarkey
Copy link
Author

Hi @hmaarrfk, thanks for getting back to me.

Before I create a test, it might be worth getting some feedback from people on whether the aim of skimage is to replicate the behaviour of MATLAB or not. I spoke to someone at my University yesterday who said that it's possible the off-by-one errors are not actually bugs, just a change in convention (although he wasn't certain without looking into it more).

If the consensus is that it should agree with MATLAB I'm happy to make the changes and create a test, but I'd like clarification on that before I put in the work :)

There is also the possibility that if this is a convention change and not a bug, that it will mess-up a lot of peoples code (although given the filters are also currently wrong, maybe that doesn't matter so much?) and that the forward radon transform might also need updating. Not being an expert in radon transforms though I can't be certain.

@stefanv
Copy link
Member
stefanv commented Feb 13, 2019

The aim is to behave correctly, regardless of what MATLAB does! So, we have to figure out what desired behavior is outside of any software system, and then ensure we implement that.

This module definitely needs more review, so your help is greatly valued.

@hmaarrfk
Copy link
Member

Right, my goal is to get the mathematics correct.

I agree that there are some issues with the module. I would just like to get that PR pushed in since I think it definitely addresses some of the issues.

Without a concrete little bit of code (should be about 10 lines long), it is really difficult to tell if anything would continue being wrong once PR #3067 gets merged in.

@philipstarkey
Copy link
Author

@hmaarrfk I've attached some test code to demonstrate the issue. This code:

  • Uses the MATLAB Python interface to directly call iradon in MATLAB, so that results can be compared. In case you don't have access to MATLAB, I've also included a csv of the results from MATLAB.
  • Runs a comparison between MATLAB and the current skimage.transform.iradon function, the one from PR transform: fix filter bias #3067 (which I copy/pasted into my code) and a version I've modified to make agree with MATLAB results.
  • Uses a (fake) sinogram with both an even and odd number of rows
  • Compares with a tolerance of 1e-3 (I set this large to be safe, but my modified iradon function agrees with MATLAB to better than 1e-12)

The changes I made in addition to PR #3067 were:

  1. Updated the initial calculation of mid_index
  2. Updated the initial calculation of the X,Y grids
  3. Added padding to radon_filtered before calculating the backprojections
  4. Updated the calculation for the interpolation x-coordinates (this, combined with 1 and 3 ensure the x-coordinates are now equally spaced around 0 as the padding in 3. ensures there are always an odd number of rows)

I think the key here is probably the additional padding (3.). I hadn't done this in my very first tests when reporting this issue, as it didn't seem to be necessary but it turns out it changes the results of the last row and column significantly (which I hadn't been looking at closely). Once the padding is applied, I think that the changes 1. and 4. flow naturally from that. The change in 2. might just be a convention, I'm honestly not sure.

Of course none of this is to say that MATLAB is correct. But hopefully this is enough information for you to determine if the mathematics of the skimage radon transform are correct or not.

If my proposed changes are just a change of convention in padding and/or determining the center point, then it might be worth exposing this through new keyword arguments so users can choose what to do. For example, it looks like there is possibly an alternative convention where the center is chosen as the boundary between pixels (if I'm reading the code of radontea correctly)

iradon_test.zip

@hmaarrfk
Copy link
Member
hmaarrfk commented Mar 3, 2019

Matlab creates a matlab engine for python, only to ignore numpy....

matlab.double([[t] for t in theta])

crazyness.

These potential off by one errors are a little crazy.

Ok, I kinda see your point, I would like to better understand if:

  1. Matlab round trips well.
  2. Scikit-image round trips well.

I think I need to get that PR through first. It seems like the right thing to do for now, I Just don't understand one line of it.

@hmaarrfk
Copy link
Member
hmaarrfk commented Mar 3, 2019

Seems like padding is good clue. Might have to do with the diagonal of the image. not too sure.
I'm trying to find where they talk about padding in
https://engineering.purdue.edu/~malcolm/pct/CTI_Ch03.pdf
https://engineering.purdue.edu/~malcolm/pct/CTI_Ch04.pdf

Also, seems like somebody else is looking at these kinds of errors:
https://stackoverflow.com/questions/50527580/inverse-radon-transformation-in-python

@AlbertZhangHIT
Copy link
AlbertZhangHIT commented Apr 12, 2019 •
edited
Loading

Hi, I also encountered a problem that may be caused by the iradon function. I tried to reconstruct an images in dimension of 512x512 from its radon measure in 50 views, i.e., theta=np.linspace(0., 180., 50, endpoint=False); rec = iradon(radon(image, theta, circle=False), theta, output_size=512, circle=False). The reconstruction is quite different from its matlab version, using theta=linspace(0, 180*(1-1/50), 50); rec = iradon(radon(image, theta), theta, 512). I checked the thetas, they are the same.

Note: There are some difference between the radons. 92C6 The dimension of those outputs are different: 725x50 for skimage radon, 729x50 for matlab radon so I cannot directly check the difference of radon measure.

Version Information
3.6.4 |Anaconda, Inc.| (default, Jan 16 2018, 10:22:32) [MSC v.1900 64 bit (AMD64)]
Windows-10-10.0.17134-SP0
scikit-image version: 0.14.1
numpy version: 1.15.4

@scikit-image scikit-image locked and limited conversation to collaborators Oct 18, 2021
@grlee77 grlee77 reopened this Mar 21, 2022
@lagru lagru added the 🐛 Bug label Sep 16, 2023
Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Assignees
No one assigned
Labels
⚖️ Needs decision 🐛 Bug
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
8 participants
@stefanv @hmaarrfk @jarrodmillman @philipstarkey @rfezzani @grlee77 @lagru @AlbertZhangHIT
0