-
Notifications
You must be signed in to change notification settings - Fork 97
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Implement
medfilt into AzimuthalIntegrator
- Loading branch information
Showing
2 changed files
with
305 additions
and
5 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -30,7 +30,7 @@ | |
| __contact__ = "[email protected]" | ||
| __license__ = "MIT" | ||
| __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" | ||
| __date__ = "10/10/2024" | ||
| __date__ = "05/12/2024" | ||
| __status__ = "stable" | ||
| __docformat__ = 'restructuredtext' | ||
|
|
||
|
|
@@ -49,6 +49,7 @@ | |
| from ..io.ponifile import PoniFile | ||
| error = None | ||
| from ..method_registry import IntegrationMethod | ||
| from ..utils.decorators import deprecated | ||
| # | ||
| from .load_engines import ocl_sort | ||
| #ocl_azim_csr, ocl_azim_lut, , histogram, splitBBox, \ | ||
|
|
@@ -1147,8 +1148,8 @@ def integrate2d_ng(self, data, npt_rad, npt_azim=360, | |
|
|
||
| integrate2d = _integrate2d_ng = integrate2d_ng | ||
|
|
||
|
|
||
| def medfilt1d(self, data, npt_rad=1024, npt_azim=512, | ||
| @deprecated(since_version="2024.12", only_once=True, replacement="medfilt1d_ng", deprecated_since="2024.12") | ||
| def medfilt1d_legacy(self, data, npt_rad=1024, npt_azim=512, | ||
| correctSolidAngle=True, | ||
| radial_range=None, azimuth_range=None, | ||
| polarization_factor=None, dark=None, flat=None, | ||
|
|
@@ -1285,6 +1286,305 @@ def medfilt1d(self, data, npt_rad=1024, npt_azim=512, | |
| result._set_normalization_factor(normalization_factor) | ||
| return result | ||
|
|
||
| medfilt1d = medfilt1d_legacy | ||
|
|
||
| def medfilt1d_ng(self, data, | ||
| npt=1024, | ||
| correctSolidAngle=True, | ||
| polarization_factor=None, | ||
| variance=None, | ||
| error_model=ErrorModel.NO, | ||
| radial_range=None, | ||
| azimuth_range=None, | ||
| dark=None, | ||
| flat=None, | ||
| absorption=None, | ||
| method=("full", "csr", "cython"), | ||
| unit=units.Q, | ||
| percentile=50, | ||
| dummy=None, | ||
| delta_dummy=None, | ||
| mask=None, | ||
| normalization_factor=1.0, | ||
| metadata=None, | ||
| safe=True, | ||
| **kwargs): | ||
| """Performs iteratively the 1D integration with variance propagation | ||
| and performs a sigm-clipping at each iteration, i.e. | ||
| all pixel which intensity differs more than thres*std is | ||
| discarded for next iteration. | ||
| Keep only pixels with intensty: | ||
| ``|I - <I>| < thres * σ(I)`` | ||
| This enforces a symmetric, bell-shaped distibution (i.e. gaussian-like) | ||
| and is very good at extracting background or amorphous isotropic scattering | ||
| out of Bragg peaks. | ||
| :param data: input image as numpy array | ||
| :param npt_rad: number of radial points | ||
| :param bool correctSolidAngle: correct for solid angle of each pixel if True | ||
| :param float polarization_factor: polarization factor between: | ||
| -1 (vertical) | ||
| +1 (horizontal). | ||
| - 0 for circular polarization or random, | ||
| - None for no correction, | ||
| - True for using the former correction | ||
| :param radial_range: The lower and upper range of the radial unit. If not provided, range is simply (data.min(), data.max()). Values outside the range are ignored. | ||
| :type radial_range: (float, float), optional | ||
| :param azimuth_range: The lower and upper range of the azimuthal angle in degree. If not provided, range is simply (data.min(), data.max()). Values outside the range are ignored. | ||
| :type azimuth_range: (float, float), optional | ||
| :param ndarray dark: dark noise image | ||
| :param ndarray flat: flat field image | ||
| :param ndarray absorption: Detector absorption (image) | ||
| :param ndarray variance: the variance of the signal | ||
| :param str error_model: can be "poisson" to assume a poissonian detector (variance=I) or "azimuthal" to take the std² in each ring (better, more expenive) | ||
| :param unit: unit to be used for integration | ||
| :param IntegrationMethod method: IntegrationMethod instance or 3-tuple with (splitting, algorithm, implementation) | ||
| :param percentile: which percentile use for cutting out. | ||
| percentil can be a 2-tuple to specify a region to average out, | ||
| like: (25,75) to average the second and third quartile. | ||
| :param mask: masked out pixels array | ||
| :param float normalization_factor: Value of a normalization monitor | ||
| :param metadata: any other metadata, | ||
| :type metadata: JSON serializable dict | ||
| :param safe: set to False to skip some tests | ||
| :return: Integrate1D like result like | ||
| The difference with the previous `medfilt_legacy` implementation is that there is no 2D regrouping. | ||
| """ | ||
| for k in kwargs: | ||
| if k == "npt_azim": | ||
| logger.warning("'npt_azim' argument is not used in sigma_clip_ng as not 2D intergration is performed anymore") | ||
| else: | ||
| logger.warning("Got unknown argument %s %s", k, kwargs[k]) | ||
|
|
||
| error_model = ErrorModel.parse(error_model) | ||
| if variance is not None: | ||
| assert variance.size == data.size | ||
| error_model = ErrorModel.VARIANCE | ||
|
|
||
| unit = units.to_unit(unit) | ||
| if radial_range: | ||
| radial_range = tuple(radial_range[i] / unit.scale for i in (0, -1)) | ||
| if azimuth_range is not None: | ||
| azimuth_range = self.normalize_azimuth_range(azimuth_range) | ||
| try: | ||
| quant_min = min(percentile)/100 | ||
| quant_max = max(percentile)/100 | ||
| except: | ||
| quant_min = quant_max = percentile/100.0 | ||
|
|
||
| method = self._normalize_method(method, dim=1, default=self.DEFAULT_METHOD_1D) | ||
|
|
||
| if mask is None: | ||
| has_mask = "from detector" | ||
| mask = self.mask | ||
| mask_crc = self.detector.get_mask_crc() | ||
| if mask is None: | ||
| has_mask = False | ||
| mask_crc = None | ||
| else: | ||
| has_mask = "user provided" | ||
| mask = numpy.ascontiguousarray(mask) | ||
| mask_crc = crc32(mask) | ||
|
|
||
| if dark is None: | ||
| dark = self.detector.darkcurrent | ||
| if dark is None: | ||
| has_dark = False | ||
| else: | ||
| has_dark = "from detector" | ||
| else: | ||
| has_dark = "provided" | ||
|
|
||
| if flat is None: | ||
| flat = self.detector.flatfield | ||
| if dark is None: | ||
| has_flat = False | ||
| else: | ||
| has_flat = "from detector" | ||
| else: | ||
| has_flat = "provided" | ||
|
|
||
| if correctSolidAngle: | ||
| solidangle = self.solidAngleArray(data.shape, correctSolidAngle) | ||
| else: | ||
| solidangle = None | ||
|
|
||
| if polarization_factor is None: | ||
| polarization = polarization_crc = None | ||
| else: | ||
| polarization, polarization_crc = self.polarization(data.shape, polarization_factor, with_checksum=True) | ||
|
|
||
| if (method.algo_lower == "csr"): | ||
| "This is the only method implemented for now ..." | ||
| # Prepare LUT if needed! | ||
| # initialize the CSR integrator in Cython as it may be needed later on. | ||
| cython_method = IntegrationMethod.select_method(method.dimension, method.split_lower, method.algo_lower, "cython")[0] | ||
| if cython_method not in self.engines: | ||
| cython_engine = self.engines[cython_method] = Engine() | ||
| else: | ||
| cython_engine = self.engines[cython_method] | ||
| with cython_engine.lock: | ||
| # Validate that the engine used is the proper one | ||
| cython_integr = cython_engine.engine | ||
| cython_reset = None | ||
| if cython_integr is None: | ||
| cython_reset = "of first initialization" | ||
| if (not cython_reset) and safe: | ||
| if cython_integr.unit != unit: | ||
| cython_reset = "unit was changed" | ||
| if cython_integr.bins != npt: | ||
| cython_reset = "number of points changed" | ||
| if cython_integr.size != data.size: | ||
| cython_reset = "input image size changed" | ||
| if cython_integr.empty != self._empty: | ||
| cython_reset = "empty value changed " | ||
| if (mask is not None) and (not cython_integr.check_mask): | ||
| cython_reset = "mask but CSR was without mask" | ||
| elif (mask is None) and (cython_integr.check_mask): | ||
| cython_reset = "no mask but CSR has mask" | ||
| elif (mask is not None) and (cython_integr.mask_checksum != mask_crc): | ||
| cython_reset = "mask changed" | ||
| if (radial_range is None) and (cython_integr.pos0_range is not None): | ||
| cython_reset = "radial_range was defined in CSR" | ||
| elif (radial_range is not None) and cython_integr.pos0_range != (min(radial_range), max(radial_range)): | ||
| cython_reset = "radial_range is defined but not the same as in CSR" | ||
| if (azimuth_range is None) and (cython_integr.pos1_range is not None): | ||
| cython_reset = "azimuth_range not defined and CSR had azimuth_range defined" | ||
| elif (azimuth_range is not None) and cython_integr.pos1_range != (min(azimuth_range), max(azimuth_range)): | ||
| cython_reset = "azimuth_range requested and CSR's azimuth_range don't match" | ||
| if cython_reset: | ||
| logger.info("AI.sigma_clip_ng: Resetting Cython integrator because %s", cython_reset) | ||
| split = method.split_lower | ||
| if split == "pseudo": | ||
| split = "full" | ||
| try: | ||
| cython_integr = self.setup_sparse_integrator(data.shape, npt, mask=mask, | ||
| mask_checksum=mask_crc, | ||
| unit=unit, split=split, algo="CSR", | ||
| pos0_range=radial_range, | ||
| pos1_range=azimuth_range, | ||
| empty=self._empty, | ||
| scale=False) | ||
| except MemoryError: # CSR method is hungry... | ||
| logger.warning("MemoryError: falling back on forward implementation") | ||
| cython_integr = None | ||
| self.reset_engines() | ||
| method = self.DEFAULT_METHOD_1D | ||
| else: | ||
| cython_engine.set_engine(cython_integr) | ||
| if method not in self.engines: | ||
| # instanciated the engine | ||
| engine = self.engines[method] = Engine() | ||
| else: | ||
| engine = self.engines[method] | ||
| with engine.lock: | ||
| # Validate that the engine used is the proper one | ||
| integr = engine.engine | ||
| reset = None | ||
| # This whole block uses CSR, Now we should treat all the various implementation: Cython, OpenCL and finally Python. | ||
|
|
||
| # Validate that the engine used is the proper one | ||
| if integr is None: | ||
| reset = "of first initialization" | ||
| if (not reset) and safe: | ||
| if integr.unit != unit: | ||
| reset = "unit was changed" | ||
| if integr.bins != npt: | ||
| reset = "number of points changed" | ||
| if integr.size != data.size: | ||
| reset = "input image size changed" | ||
| if integr.empty != self._empty: | ||
| reset = "empty value changed " | ||
| if (mask is not None) and (not integr.check_mask): | ||
| reset = "mask but CSR was without mask" | ||
| elif (mask is None) and (integr.check_mask): | ||
| reset = "no mask but CSR has mask" | ||
| elif (mask is not None) and (integr.mask_checksum != mask_crc): | ||
| reset = "mask changed" | ||
| # TODO | ||
| if (radial_range is None) and (integr.pos0_range is not None): | ||
| reset = "radial_range was defined in CSR" | ||
| elif (radial_range is not None) and integr.pos0_range != (min(radial_range), max(radial_range)): | ||
| reset = "radial_range is defined but not the same as in CSR" | ||
| if (azimuth_range is None) and (integr.pos1_range is not None): | ||
| reset = "azimuth_range not defined and CSR had azimuth_range defined" | ||
| elif (azimuth_range is not None) and integr.pos1_range != (min(azimuth_range), max(azimuth_range)): | ||
| reset = "azimuth_range requested and CSR's azimuth_range don't match" | ||
|
|
||
| if reset: | ||
| logger.info("ai.sigma_clip_ng: Resetting ocl_csr integrator because %s", reset) | ||
| csr_integr = self.engines[cython_method].engine | ||
| if method.impl_lower == "opencl": | ||
| try: | ||
| integr = method.class_funct_ng.klass(csr_integr.lut, | ||
| image_size=data.size, | ||
| checksum=csr_integr.lut_checksum, | ||
| empty=self._empty, | ||
| unit=unit, | ||
| mask_checksum=csr_integr.mask_checksum, | ||
| bin_centers=csr_integr.bin_centers, | ||
| platformid=method.target[0], | ||
| deviceid=method.target[1]) | ||
| except MemoryError: | ||
| logger.warning("MemoryError: falling back on default forward implementation") | ||
| self.reset_engines() | ||
| method = self.DEFAULT_METHOD_1D | ||
| else: | ||
| # Copy some properties from the cython integrator | ||
| integr.pos0_range = csr_integr.pos0_range | ||
| integr.pos1_range = csr_integr.pos1_range | ||
| engine.set_engine(integr) | ||
| elif method.impl_lower in ("python", "cython"): | ||
| integr = method.class_funct_ng.klass(lut=csr_integr.lut, | ||
| image_size=data.size, | ||
| empty=self._empty, | ||
| unit=unit, | ||
| mask_checksum=csr_integr.mask_checksum, | ||
| bin_centers=csr_integr.bin_centers) | ||
| # Copy some properties from the cython integrator | ||
| integr.pos0_range = csr_integr.pos0_range | ||
| integr.pos1_range = csr_integr.pos1_range | ||
| engine.set_engine(integr) | ||
| else: | ||
| logger.error(f"Implementation {method.impl_lower} not supported") | ||
| else: | ||
| integr = self.engines[method].engine | ||
| kwargs = {"dark":dark, "dummy":dummy, "delta_dummy":delta_dummy, | ||
| "variance":variance, "dark_variance":None, | ||
| "flat":flat, "solidangle":solidangle, "polarization":polarization, "absorption":absorption, | ||
| "error_model":error_model, "normalization_factor":normalization_factor, | ||
| "quant_min":quant_min, "quant_max":quant_max} | ||
|
|
||
| intpl = integr.medfilt(data, **kwargs) | ||
| else: | ||
| raise RuntimeError("Not yet implemented. Sorry") | ||
| result = Integrate1dResult(intpl.position * unit.scale, intpl.intensity, intpl.sem) | ||
| result._set_method_called("sigma_clip_ng") | ||
| result._set_method(method) | ||
| result._set_compute_engine(str(method)) | ||
| result._set_percentile(percentile) | ||
| result._set_unit(unit) | ||
| result._set_has_mask_applied(has_mask) | ||
| result._set_has_dark_correction(has_dark) | ||
| result._set_has_flat_correction(has_flat) | ||
| result._set_metadata(metadata) | ||
| result._set_sum_signal(intpl.signal) | ||
| result._set_sum_normalization(intpl.normalization) | ||
| result._set_sum_normalization2(intpl.norm_sq) | ||
| result._set_std(intpl.std) | ||
| result._set_sem(intpl.sem) | ||
| result._set_sum_variance(intpl.variance) | ||
| result._set_count(intpl.count) | ||
| result._set_polarization_factor(polarization_factor) | ||
| result._set_normalization_factor(normalization_factor) | ||
| result._set_error_model(error_model) | ||
| return result | ||
|
|
||
| def sigma_clip_legacy(self, data, npt_rad=1024, npt_azim=512, | ||
| correctSolidAngle=True, polarization_factor=None, | ||
| radial_range=None, azimuth_range=None, | ||
|
|
||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters