feat: spatial resolution measurement #107
Conversation
Example using fake dataOriginal image "Measured image" (fake obtained by convolution the ideal image by a point-spread-function) Computed modulation transfer function we can see the cutoff frequency in this example is at about 0.2 lines per pixel. |
|
This still needs some more work - mainly applying it to real test data to check that it work in practice. It also needs tests. |
|
As we discussed earlier, we should just assume the input is:
or
Then what the interface/notebook should do becomes quite simple right...? |
|
Also the resolution should be dependent on the wavelength in principal but this function doesn't care if its 2D or 3D right...? or should be broadcast the perfect image first...? |
I think this makes most sense, the cropping can be done elsewhere. If the resolution is different for different wavelengths then the resolution function can be called multiple times, once for each wavelength. |
nvaytet
left a comment
nvaytet
left a comment
There was a problem hiding this comment.
Changes look fine, but we should try and not have any 'magic' numbers without adding explanations about where they came from.
src/ess/imaging/tools/resolution.py
Outdated
| ) | ||
|
|
||
|
|
||
| def _radial_profile(data: NDArray): |
|
|
||
|
|
||
| def _radial_profile(data: NDArray): | ||
| '''Integrate ellipses around center of image.''' |
There was a problem hiding this comment.
Does this mean we need to center the image before giving it to the function?
Can this be tricky in some cases?
Would passing the (x,y) coordinates of the center as arguments (we can set them to 0 by default) make things easier in some cases?
There was a problem hiding this comment.
It's not necessary, it's applied in the frequency domain and then we get a centered image automatically.
| _measured = measured_image.values | ||
| # Can't do inplace because dtype of sum might be different from dtype of input | ||
| _measured = _measured / _measured.sum() | ||
| _reference = (open_beam_image * target).to(unit=measured_image.unit).values |
There was a problem hiding this comment.
Can you explain why we are multiplying target by open beam image?
Alternatively, can you add a link in the docstring where you got the implementation from (paper or online article).
There was a problem hiding this comment.
There is a motivation for the procedure on this page: https://git.esss.dk/dram/demo-drawer/odin-demo/-/blob/main/STAP-2025-fall/siemens-star-resolution-estimation.ipynb?ref_type=heads
Unfortunately because of gitlab rendering limits the page is not displayed correctly.
Here is a screenshot of the relevant section:
There was a problem hiding this comment.
Can you add all this in the docstring? It would be useful to have all the math in the API reference docs.
It can be under a Notes section in the docstring if you like, so that we first have all the parameters and return types?
| # a kind of counts I think in our unit system that is best | ||
| # represented as 'dimensionless'. | ||
| coords={'frequency': sc.linspace('frequency', 0, (1 / 2) ** 0.5, len(_mtf))}, | ||
| # We're only interested in frequencies below 0.5 oscillations per pixel |
There was a problem hiding this comment.
As suggested above, a link to an article would help understand where the number 0.5 comes from.
There was a problem hiding this comment.
(1/2)**0.5 is the largest magnitude of the spatial frequency returned by the 2D FFT. It's one period per one diagonal step in the grid.
0.5 is the largest signal frequency magnitude that can be represented in every direction of the 2D grid.
We want the MTF to be independent of direction, so the higher frequencies are not really "physical", that's why it's limited to 0.5.
There was a problem hiding this comment.
Adding this explanation as a comment would be useful 👍
|
|
||
| def estimate_cut_off_frequency(mtf: sc.DataArray) -> sc.Variable: | ||
| '''Estimates the cut off frequency of | ||
| the modulation transfer function (mtf). |
There was a problem hiding this comment.
Can you add a link to a reference with the method to estimate the cut-off frequency in the docstring?
There was a problem hiding this comment.
I don't have a reference for it. It's just something I came up with. Didn't find anything after making a quick search for it, and figured this is good enough. It's not like we can exactly determine the cutoff frequency in a noisy experiment anyway.
| maxiters = 100 | ||
| for _ in range(maxiters): | ||
| p = np.polyfit(_freq[m], _mtf[m], 1, w=w[m]) | ||
| if abs(-p[1] / p[0] - fc) < 1e-4: |
There was a problem hiding this comment.
Should the threshold parameter be an input argument that can be tuned by the user?
If not, why was the number 1e-4 chosen?
There was a problem hiding this comment.
I don't think it makes sense for the user to tune that parameter. The parameter represents the difference between the previous iteration and the next, so reducing it means we don't stop the iteration before the difference is even smaller than 1e-4. But the method does not produce a result with 4 significant digits anyway, so reducing it does not reduce the error.
You could increase the threshold, to reduce the number of iterations and reduce runtime, but it's not like the algorithm takes any time as it is, didn't measure the runtime but it's in ms. So that doesn't make much sense from the users perspective.
If we get requests for this in the future we can always add it later. I think it's better to start with a simpler interface and add complexity as requested, than it is to start with a complex interface.
There was a problem hiding this comment.
I agree that a simple interface is better. Can you just add a small comment as to why this number works well? (basically the end of the first paragraph in your response?)
| _freq = np.concat([[0.0], mtf.coords['frequency'].values]) | ||
| _mtf = np.concat([[1.0], mtf.values]) | ||
| # The line should go through (0, 1), so give it a big weight. | ||
| w = np.concat([[10 * len(mtf)], np.ones(len(mtf))]) |
There was a problem hiding this comment.
The weight should be something that is significantly larger than the weight of all other points together.
But not so large that we run into numerical issues.
I picked 10 x total_weight and it seems to work well, I didn't try other values.
There was a problem hiding this comment.
Add this explanation as a comment would still be useful. It would at least make it clear to someone reading this that the number 10 may not be a hard requirement, but a value someone can experiment with.
Saying that a value was obtained from trial and error is completely legitimate.
| "siemens-star-measured.h5": "md5:8e333d36c7c102f474b2b66cb785f5e8", | ||
| "siemens-star-openbeam.h5": "md5:ee429b2c247aeaafb0ef3ca4171f2e6a", |
There was a problem hiding this comment.
I think we have a similar issue with some of the other files in here, so it's not specific to these changes: we should probably add descriptions on how these images were created (link to a script we used if we have that, or link to where we got the image from).
Co-authored-by: Neil Vaytet <39047984+nvaytet@users.noreply.github.com>
Fixes #17
We don't have perfect requirements here yet, but I have some idea about what is needed, so started implementing that.
The "one number to estimate resolution" implemented is
f_cthe cutoff frequency of the MTF (modulation transfer function) and MTF50 the frequency where MTF decreases below0.5.