silx-kit / pyFAI

Fast Azimuthal Integration in Python
Other
104 stars 94 forks source link

medfilt1d: issue in OpenCL implementation #2187

Closed pierrepaleo closed 2 months ago

pierrepaleo commented 2 months ago

Hi,

There seems to be a problem with medfilt1d when using the OpenCL implementation.
The last values can be arbitrary (eg. 1e38) - probably a symptom of memory allocated but never written. This issue does not occur with Cython. Setting dummy=0 in medfilt kwargs mitigates the issue.

Attached is an example - sorry I could not easily make it simpler on dummy data. I can attach the necessary files (mask, poni and detector) if needed.

Example code (click to expand) ```python import matplotlib.pyplot as plt from silx.io.url import DataUrl from silx.io import get_data import pyFAI.azimuthalIntegrator from pyFAI.method_registry import IntegrationMethod def parse_poni_file(poni_file): conf = {} with open(poni_file) as opened_file: for line in opened_file: line = line.strip() if line.startswith("#") or (":" not in line): continue words = line.split(":", 1) key = words[0].strip().lower() try: value = words[1].strip() except Exception as error: print("Error %s with line: %s", error, line) conf[key] = value return conf def create_azimuthal_integrator(poni_file, detector_file, mask_file): """ Create a AzimuthalIntegrator object from a configuration. Parameters ---------- poni_file: str Path to the PONI file detector_file: str Path to the detector definition (HDF5 file) mask_file: str Path to the mask file Returns ------- azimuthal_integrator: pyFAI.azimuthalIntegrator.AzimuthalIntegrator Azimuthal integrator instance """ # Create azimuthal integrator. # pyFAI.load(ai_config.poni_file) cannot be used directly if the detector file therein # can't be accessed. So we have to use this trick instead poni = parse_poni_file(poni_file) ai_kwargs = {name: float(poni[name]) for name in ["poni1", "poni2", "rot1", "rot2", "rot3", "wavelength"]} ai_kwargs["dist"] = poni["distance"] azimuthal_integrator = pyFAI.azimuthalIntegrator.AzimuthalIntegrator() for param_name, value in ai_kwargs.items(): setattr(azimuthal_integrator, param_name, value) # Set detector detector = pyFAI.detector_factory(detector_file) azimuthal_integrator.detector = detector # Set mask mask_file = mask_file if isinstance(mask_file, str): mask_file = DataUrl(file_path=mask_file, scheme="fabio") mask_array = get_data(mask_file) detector.mask = mask_array return azimuthal_integrator def select_integration_method(dim=1, split="full", algo="csr", impl="opencl", target=0): ai_methods = IntegrationMethod.select_method( dim=dim, split=split, algo=algo, impl=impl, target=target ) ai_method = ai_methods[0] return ai_method if __name__ == "__main__": poni_file = "/data/visitor/ma5629/id15a/20230425/processed/pdf.poni" detector_file = "/data/visitor/ma5629/id15a/20230425/processed/id15_pilatus.h5" mask_file = "/data/visitor/ma5629/id15a/20230425/processed/mask_for_pdf_integration.edf" azim = create_azimuthal_integrator(poni_file, detector_file, mask_file) image = get_data( DataUrl( file_path="/data/visitor/ma5629/id15a/20230425/NSI2_2_pdfct_-0.005/NSI2_2_pdfct_-0.005_0001/NSI2_2_pdfct_-0.005_0001.h5", data_path="1.1/measurement/pilatus", data_slice=slice(500, 501), scheme="silx" ) )[0] ai_method = select_integration_method(impl="opencl") # works when "cython" is selected result = azim.medfilt1d( image, unit="q_A^-1", polarization_factor=1.0, correctSolidAngle=True, method=ai_method, npt_rad=2000, npt_azim=2000, percentile=(5, 95), ) plt.figure() plt.plot(result.radial, result.intensity) plt.show() ```
kif commented 2 months ago

I manage to reproduce the bug: Figure 3 (2)

kif commented 2 months ago

The bug is linked to to the range in percentile

kif commented 2 months ago

When percentile is a 2-tuple, sorter.trimmed_mean_horizontal is called; while when percentile is a single value , sorter.filter_horizontal is called. So I guess the error is in sorter.trimmed_mean_horizontal