usnistgov / PyHyperScattering

Tools for hyperspectral x-ray and neutron scattering data loading, reduction, slicing, and visualization.
Other
6 stars 8 forks source link

Incorrect q axis in WPIntegrator in PyHyperScattering version 0.2.0 #84

Closed PriyankaKetkar closed 5 months ago

PriyankaKetkar commented 1 year ago

Bringing up a bug found in PyHyper version 0.2.0. In 0.2.0, the q values of scattering curves simulated using NRSS are shifted lower than they should be, and switching to version 0.1.8 or 0.1.10 solved the issue. Bug is likely in the Warp Polar Integrator, which shifts q values by a factor of sqrt(2). More context below if needed.

I initially was running into a problem that my NRSS simulation Bragg peaks (red) were at lower q values than my experimental RSoXS primary Bragg peak (black). It was confirmed that experimental RSoXS primary Bragg peak is matching with the AFM spacing after fourier transforming (gray line). It also was confirmed that I am entering the voxel size (PhysSize = 1.95312 nm) correctly into the NRSS simulation and also the PyHyper reader. After switching to PyHyper version 0.1.8 or 0.1.10, the red curve shifted to the right and now lines up with the black curve.

image

phonghnguyen commented 1 year ago

My guess is that it has to do with these lines of code in WPIntegrator:


        if self.MACHINE_HAS_CUDA:
            TwoD = self.warp_polar_gpu(img_to_integ,center=(center_x,center_y), radius = np.nanmax(img_to_integ.shape))
        else:
            TwoD = skimage.transform.warp_polar(img_to_integ,center=(center_x,center_y), radius = np.nanmax(img_to_integ.shape))

The q-axis seem to be implicitly generated up to the given radius and this takes the radius to be the maximum dimension of the resulting simulated qx, qy image. The actual radius for a square image should be proportionally sqrt(2)/2 as per this drawing:

IMG_4844

Multiplying by sqrt(2)/2 would only work for square qx, qy images. A suitable fix that maintains functionality with non-square detector images might be something along these lines:

        if self.MACHINE_HAS_CUDA:
            TwoD = self.warp_polar_gpu(img_to_integ,center=(center_x,center_y), radius = np.sqrt((img_to_integ.shape[0] - center_x)**2 + (img_to_integ.shape[1] - center_y)**2))
        else:
            TwoD = skimage.transform.warp_polar(img_to_integ,center=(center_x,center_y), radius = np.sqrt((img_to_integ.shape[0] - center_x)**2 + (img_to_integ.shape[1] - center_y)**2))
pbeaucage commented 1 year ago

This looks like it should work, thanks for figuring it out so quickly! @phonghnguyen can you put this into a PR so that @PriyankaKetkar can test?

pdudenas commented 1 year ago

The current fix looks to have the default behavior when radius is left as None. Should we implement this as a kwarg instead of hardcoded?

phonghnguyen commented 1 year ago

I agree that having radius supplied as a kwarg would be best. There might be other cases where a different specification of the radius might make sense. We could include examples of different cases in the documentation e.g., non-square images or when the beam center is not at the center of the image).

pbeaucage commented 5 months ago

I believe this was fixed in #85. If not, please feel free to reopen.