From 6265ae75064645bf87a66270ac39c4378ea47480 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Thu, 11 Apr 2024 10:08:39 +0100 Subject: [PATCH 01/10] Add variance over mean ratio to FindSXPEaksConvolve For a constant background the Poisson distributed counts should have Variance = mean (i.e. ratio = 1) - in the region of a Bragg peak the variance/mean > 1. This idea is taken from DIALS software package for x-ray diffraction as described Winter (2018) https://doi.org/10.1107/S2059798317017235 - generalised to 3D TOF Laue --- .../plugins/algorithms/FindSXPeaksConvolve.py | 78 +++++++++++++------ 1 file changed, 55 insertions(+), 23 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py index c16fa3898182..3cb147ebed28 100644 --- a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py +++ b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py @@ -19,6 +19,7 @@ EnabledWhenProperty, PropertyCriterion, logger, + StringListValidator, ) import numpy as np from scipy.ndimage import label, maximum_position, binary_closing, sum_labels, uniform_filter1d @@ -109,16 +110,30 @@ def PyInit(self): validator=FloatBoundedValidator(lower=0.0), doc="Minimum peak size as a fraction of the kernel size.", ) + self.declareProperty( + name="PeakFindingStrategy", + defaultValue="IOverSigma", + direction=Direction.Input, + validator=StringListValidator(["VarianceOverMean", "IOverSigma"]), + doc="To-Do.", + ) + self.declareProperty( + name="ThresholdVarianceOverMean", + defaultValue=3.0, + direction=Direction.Input, + validator=FloatBoundedValidator(lower=1.0), + doc="To-Do", + ) def PyExec(self): # get input ws = self.getProperty("InputWorkspace").value - threshold_i_over_sig = self.getProperty("ThresholdIoverSigma").value nrows = self.getProperty("NRows").value ncols = self.getProperty("NCols").value nfwhm = self.getProperty("NFWHM").value get_nbins_from_b2bexp_params = self.getProperty("GetNBinsFromBackToBackParams").value min_frac_size = self.getProperty("MinFracSize").value + peak_finding_strategy = self.getProperty("PeakFindingStrategy").value # create output table workspace peaks = self.exec_child_alg("CreatePeaksWorkspace", InstrumentWorkspace=ws, NumberOfPeaks=0, OutputWorkspace="_peaks") @@ -148,20 +163,31 @@ def PyExec(self): # get data in detector coords peak_data = array_converter.get_peak_data(dummy_pk, detid, bank.getName(), bank.xpixels(), bank.ypixels(), 1, 1) _, y, esq, _ = peak_data.get_data_arrays() # 3d arrays [rows x cols x tof] - # perform convolutions to integrate kernel/shoebox - # pad with nearest so don't get peaks at edge when -ve values go outside data extent - kernel = make_kernel(nrows, ncols, nbins) - - yconv = convolve(y, kernel, mode="valid") - econv = np.sqrt(convolve(esq, kernel**2, mode="valid")) - - with np.errstate(divide="ignore", invalid="ignore"): - intens_over_sig = yconv / econv # ignore 0/0 which produces NaN (recall NaN > x = False) - intens_over_sig[~np.isfinite(intens_over_sig)] = 0 - intens_over_sig = uniform_filter1d(intens_over_sig, size=3, axis=2, mode="nearest") - - # find peaks above threshold I/sigma - pk_mask = intens_over_sig > threshold_i_over_sig + if peak_finding_strategy == "IOverSigma": + kernel = make_kernel(nrows, ncols, nbins) # integration shoebox kernel + yconv = convolve(y, kernel, mode="valid") + econv = np.sqrt(convolve(esq, kernel**2, mode="valid")) + with np.errstate(divide="ignore", invalid="ignore"): + ratio = yconv / econv # ignore 0/0 which produces NaN (recall NaN > x = False) + threshold = self.getProperty("ThresholdIoverSigma").value + else: + # Variance = (E[X^2] - E[X]^2) + # Evaluate ratio of Variance/mean (should be 1 for Poisson distribution if constant bg) + kernel = np.ones((nrows, ncols, nbins)) / (nrows * ncols * nbins) + # scale to raw counts + with np.errstate(divide="ignore", invalid="ignore"): + scale = y / esq + scale[~np.isfinite(scale)] = 0 + ycnts = y * scale + avg_y = convolve(ycnts, kernel, mode="valid") + avg_ysq = convolve(ycnts**2, kernel, mode="valid") + ratio = (avg_ysq - avg_y**2) / avg_y + threshold = self.getProperty("ThresholdVarianceOverMean").value + # perform final smoothing + ratio[~np.isfinite(ratio)] = 0 + ratio = uniform_filter1d(ratio, size=3, axis=2, mode="nearest") + # identify peaks above threshold + pk_mask = ratio > threshold pk_mask = binary_closing(pk_mask) # removes holes - helps merge close peaks labels, nlabels = label(pk_mask) # identify contiguous nearest-neighbour connected regions # identify labels of peaks above min size @@ -170,15 +196,12 @@ def PyExec(self): ilabels = np.flatnonzero(npixels > min_size) + 1 # find index of maximum in I/sigma for each valid peak (label index in ilabels) - imaxs = maximum_position(intens_over_sig, labels, ilabels) + imaxs = maximum_position(ratio, labels, ilabels) # add peaks to table for ipk in range(len(imaxs)): irow, icol, itof = imaxs[ipk] - peak_conv_intens = yconv[irow, icol, itof] - sig_conv_intens = econv[irow, icol, itof] - # map convolution output to the input irow += (kernel.shape[0] - 1) // 2 icol += (kernel.shape[1] - 1) // 2 @@ -193,6 +216,12 @@ def PyExec(self): itof_lo = np.clip(itof - kernel.shape[2] // 2, a_min=0, a_max=y.shape[2]) itof_hi = np.clip(itof + kernel.shape[2] // 2, a_min=0, a_max=y.shape[2]) ypk = y[irow_lo:irow_hi, icol_lo:icol_hi, itof_lo:itof_hi] + if peak_finding_strategy == "VarianceOverMean": + # perform additional check as in Winter (2018) https://doi.org/10.1107/S2059798317017235 + background_cnts = avg_y[tuple(imaxs[ipk])] + ypk_cnts = ycnts[irow_lo:irow_hi, icol_lo:icol_hi, itof_lo:itof_hi] + if not np.any(ypk_cnts > background_cnts + np.sqrt(background_cnts)): + continue # skip peak # integrate over TOF and select detector with max intensity imax_det = np.argmax(ypk.sum(axis=2)) irow_max, icol_max = np.unravel_index(imax_det, ypk.shape[0:-1]) @@ -209,10 +238,13 @@ def PyExec(self): DetectorID=int(peak_data.detids[irow_max, icol_max]), ) # set intensity of peak (rough estimate) - pk = peaks.getPeak(peaks.getNumberPeaks() - 1) - bin_width = xspec[itof + 1] - xspec[itof] - pk.setIntensity(peak_conv_intens * bin_width) - pk.setSigmaIntensity(sig_conv_intens * bin_width) + if peak_finding_strategy == "IOverSigma": + peak_conv_intens = yconv[tuple(imaxs[ipk])] + sig_conv_intens = econv[tuple(imaxs[ipk])] + pk = peaks.getPeak(peaks.getNumberPeaks() - 1) + bin_width = xspec[itof + 1] - xspec[itof] + pk.setIntensity(peak_conv_intens * bin_width) + pk.setSigmaIntensity(sig_conv_intens * bin_width) # remove dummy peak self.exec_child_alg("DeleteTableRows", TableWorkspace=peaks, Rows=[irow_to_del]) From e26041daee0f899621ef9f8b9104ee9b628acab0 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 15 Apr 2024 13:05:57 +0100 Subject: [PATCH 02/10] Implement usign 1D convolutions (with uniform_filter) for performance --- .../plugins/algorithms/FindSXPeaksConvolve.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py index 3cb147ebed28..8e8c274ea89c 100644 --- a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py +++ b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py @@ -22,7 +22,7 @@ StringListValidator, ) import numpy as np -from scipy.ndimage import label, maximum_position, binary_closing, sum_labels, uniform_filter1d +from scipy.ndimage import label, maximum_position, binary_closing, sum_labels, uniform_filter1d, uniform_filter from scipy.signal import convolve from IntegratePeaksSkew import InstrumentArrayConverter, get_fwhm_from_back_to_back_params @@ -179,9 +179,13 @@ def PyExec(self): scale = y / esq scale[~np.isfinite(scale)] = 0 ycnts = y * scale - avg_y = convolve(ycnts, kernel, mode="valid") - avg_ysq = convolve(ycnts**2, kernel, mode="valid") + avg_y = uniform_filter(ycnts, size=(nrows, ncols, nbins), mode="nearest") + avg_ysq = uniform_filter(ycnts**2, size=(nrows, ncols, nbins), mode="nearest") ratio = (avg_ysq - avg_y**2) / avg_y + # crop to valid region + ratio = ratio[ + (nrows - 1) // 2 : -(nrows - 1) // 2, (ncols - 1) // 2 : -(ncols - 1) // 2, (nbins - 1) // 2 : -(nbins - 1) // 2 + ] threshold = self.getProperty("ThresholdVarianceOverMean").value # perform final smoothing ratio[~np.isfinite(ratio)] = 0 From eb641f76cb66198a87edb5b98c586726d08447f7 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 15 Apr 2024 14:16:18 +0100 Subject: [PATCH 03/10] Add argument doc strings and refactor to reduce size of main exec --- .../plugins/algorithms/FindSXPeaksConvolve.py | 80 +++++++++++-------- 1 file changed, 48 insertions(+), 32 deletions(-) diff --git a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py index 8e8c274ea89c..afa64dcb97d9 100644 --- a/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py +++ b/Framework/PythonInterface/plugins/algorithms/FindSXPeaksConvolve.py @@ -115,15 +115,21 @@ def PyInit(self): defaultValue="IOverSigma", direction=Direction.Input, validator=StringListValidator(["VarianceOverMean", "IOverSigma"]), - doc="To-Do.", + doc="PeakFindingStrategy=IOverSigma will find peaks by integrating data using a shoebox kernel and looking" + " for peaks with I/sigma > ThresholdIoverSigma. PeakFindingStrategy=VarianceOverMean will look for " + "peaks with local variance/mean > ThresholdVarianceOverMean.", ) self.declareProperty( name="ThresholdVarianceOverMean", defaultValue=3.0, direction=Direction.Input, validator=FloatBoundedValidator(lower=1.0), - doc="To-Do", + doc="Threshold value for variance/mean used to identify peaks.", ) + use_variance_over_mean = EnabledWhenProperty("PeakFindingStrategy", PropertyCriterion.IsNotDefault) + self.setPropertySettings("ThresholdVarianceOverMean", use_variance_over_mean) + use_intens_over_sigma = EnabledWhenProperty("PeakFindingStrategy", PropertyCriterion.IsDefault) + self.setPropertySettings("ThresholdIoverSigma", use_intens_over_sigma) def PyExec(self): # get input @@ -164,29 +170,11 @@ def PyExec(self): peak_data = array_converter.get_peak_data(dummy_pk, detid, bank.getName(), bank.xpixels(), bank.ypixels(), 1, 1) _, y, esq, _ = peak_data.get_data_arrays() # 3d arrays [rows x cols x tof] if peak_finding_strategy == "IOverSigma": - kernel = make_kernel(nrows, ncols, nbins) # integration shoebox kernel - yconv = convolve(y, kernel, mode="valid") - econv = np.sqrt(convolve(esq, kernel**2, mode="valid")) - with np.errstate(divide="ignore", invalid="ignore"): - ratio = yconv / econv # ignore 0/0 which produces NaN (recall NaN > x = False) threshold = self.getProperty("ThresholdIoverSigma").value + ratio, yconv, econv, kernel_shape = calculate_intens_over_sigma(y, esq, nrows, ncols, nbins) else: - # Variance = (E[X^2] - E[X]^2) - # Evaluate ratio of Variance/mean (should be 1 for Poisson distribution if constant bg) - kernel = np.ones((nrows, ncols, nbins)) / (nrows * ncols * nbins) - # scale to raw counts - with np.errstate(divide="ignore", invalid="ignore"): - scale = y / esq - scale[~np.isfinite(scale)] = 0 - ycnts = y * scale - avg_y = uniform_filter(ycnts, size=(nrows, ncols, nbins), mode="nearest") - avg_ysq = uniform_filter(ycnts**2, size=(nrows, ncols, nbins), mode="nearest") - ratio = (avg_ysq - avg_y**2) / avg_y - # crop to valid region - ratio = ratio[ - (nrows - 1) // 2 : -(nrows - 1) // 2, (ncols - 1) // 2 : -(ncols - 1) // 2, (nbins - 1) // 2 : -(nbins - 1) // 2 - ] threshold = self.getProperty("ThresholdVarianceOverMean").value + ratio, ycnts, avg_y, kernel_shape = calculate_variance_over_mean(y, esq, nrows, ncols, nbins) # perform final smoothing ratio[~np.isfinite(ratio)] = 0 ratio = uniform_filter1d(ratio, size=3, axis=2, mode="nearest") @@ -195,7 +183,7 @@ def PyExec(self): pk_mask = binary_closing(pk_mask) # removes holes - helps merge close peaks labels, nlabels = label(pk_mask) # identify contiguous nearest-neighbour connected regions # identify labels of peaks above min size - min_size = int(min_frac_size * kernel.size) + min_size = int(min_frac_size * np.prod(kernel_shape)) npixels = sum_labels(pk_mask, labels, range(1, nlabels + 1)) ilabels = np.flatnonzero(npixels > min_size) + 1 @@ -207,18 +195,18 @@ def PyExec(self): irow, icol, itof = imaxs[ipk] # map convolution output to the input - irow += (kernel.shape[0] - 1) // 2 - icol += (kernel.shape[1] - 1) // 2 - itof += (kernel.shape[2] - 1) // 2 + irow += (kernel_shape[0] - 1) // 2 + icol += (kernel_shape[1] - 1) // 2 + itof += (kernel_shape[2] - 1) // 2 # find peak position # get data in kernel window around index with max I/sigma - irow_lo = np.clip(irow - kernel.shape[0] // 2, a_min=0, a_max=y.shape[0]) - irow_hi = np.clip(irow + kernel.shape[0] // 2, a_min=0, a_max=y.shape[0]) - icol_lo = np.clip(icol - kernel.shape[1] // 2, a_min=0, a_max=y.shape[1]) - icol_hi = np.clip(icol + kernel.shape[1] // 2, a_min=0, a_max=y.shape[1]) - itof_lo = np.clip(itof - kernel.shape[2] // 2, a_min=0, a_max=y.shape[2]) - itof_hi = np.clip(itof + kernel.shape[2] // 2, a_min=0, a_max=y.shape[2]) + irow_lo = np.clip(irow - kernel_shape[0] // 2, a_min=0, a_max=y.shape[0]) + irow_hi = np.clip(irow + kernel_shape[0] // 2, a_min=0, a_max=y.shape[0]) + icol_lo = np.clip(icol - kernel_shape[1] // 2, a_min=0, a_max=y.shape[1]) + icol_hi = np.clip(icol + kernel_shape[1] // 2, a_min=0, a_max=y.shape[1]) + itof_lo = np.clip(itof - kernel_shape[2] // 2, a_min=0, a_max=y.shape[2]) + itof_hi = np.clip(itof + kernel_shape[2] // 2, a_min=0, a_max=y.shape[2]) ypk = y[irow_lo:irow_hi, icol_lo:icol_hi, itof_lo:itof_hi] if peak_finding_strategy == "VarianceOverMean": # perform additional check as in Winter (2018) https://doi.org/10.1107/S2059798317017235 @@ -271,6 +259,34 @@ def exec_child_alg(self, alg_name, **kwargs): return None +def calculate_intens_over_sigma(y, esq, nrows, ncols, nbins): + kernel = make_kernel(nrows, ncols, nbins) # integration shoebox kernel + yconv = convolve(y, kernel, mode="valid") + econv = np.sqrt(convolve(esq, kernel**2, mode="valid")) + with np.errstate(divide="ignore", invalid="ignore"): + intens_over_sigma = yconv / econv # ignore 0/0 which produces NaN (recall NaN > x = False) + return intens_over_sigma, yconv, econv, kernel.shape + + +def calculate_variance_over_mean(y, esq, nrows, ncols, nbins): + # Evaluate ratio of Variance/mean (should be 1 for Poisson distribution if constant bg) + # Peak finding criterion used in DIALS, Winter (2018) https://doi.org/10.1107/S2059798317017235 + with np.errstate(divide="ignore", invalid="ignore"): + # scale to raw counts + scale = y / esq + scale[~np.isfinite(scale)] = 0 + ycnts = y * scale + # calculate variance = (E[X^2] - E[X]^2) + avg_y = uniform_filter(ycnts, size=(nrows, ncols, nbins), mode="nearest") + avg_ysq = uniform_filter(ycnts**2, size=(nrows, ncols, nbins), mode="nearest") + var_over_mean = (avg_ysq - avg_y**2) / avg_y + # crop to valid region + var_over_mean = var_over_mean[ + (nrows - 1) // 2 : -(nrows - 1) // 2, (ncols - 1) // 2 : -(ncols - 1) // 2, (nbins - 1) // 2 : -(nbins - 1) // 2 + ] + return var_over_mean, ycnts, avg_y, (nrows, ncols, nbins) + + def make_kernel(nrows, ncols, nbins): # make a kernel with background subtraction shell of approx. same number of elements as peak region kernel_shape, (nrows_bg, ncols_bg, nbins_bg) = get_kernel_shape(nrows, ncols, nbins) From 917cf2ca512238e8dc5d3c8d0ee1d3afa9f79e79 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 15 Apr 2024 14:17:58 +0100 Subject: [PATCH 04/10] Add unit test --- .../algorithms/FindSXPeaksConvolveTest.py | 25 +++++++++++++------ 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/FindSXPeaksConvolveTest.py b/Framework/PythonInterface/test/python/plugins/algorithms/FindSXPeaksConvolveTest.py index 84eab925724d..62caeb9b1e8a 100644 --- a/Framework/PythonInterface/test/python/plugins/algorithms/FindSXPeaksConvolveTest.py +++ b/Framework/PythonInterface/test/python/plugins/algorithms/FindSXPeaksConvolveTest.py @@ -59,15 +59,13 @@ def setUpClass(cls): def tearDownClass(cls): AnalysisDataService.clear() - def _assert_found_correct_peaks(self, peak_ws): + def _assert_found_correct_peaks(self, peak_ws, integrated=True): self.assertEqual(peak_ws.getNumberPeaks(), 1) - peak_ws = SortPeaksWorkspace( - InputWorkspace=peak_ws, OutputWorkspace=peak_ws.name(), ColumnNameToSortBy="DetID", SortAscending=False - ) pk = peak_ws.getPeak(0) self.assertEqual(pk.getDetectorID(), 74) - self.assertAlmostEqual(pk.getTOF(), 5.0, delta=1e-8) - self.assertAlmostEqual(pk.getIntensityOverSigma(), 7.4826, delta=1e-4) + if integrated: + self.assertAlmostEqual(pk.getTOF(), 5.0, delta=1e-8) + self.assertAlmostEqual(pk.getIntensityOverSigma(), 7.4826, delta=1e-4) def test_exec_specify_nbins(self): out = FindSXPeaksConvolve( @@ -93,10 +91,23 @@ def test_exec_IoverSigma_threshold(self): def test_exec_min_frac_size(self): out = FindSXPeaksConvolve( - InputWorkspace=self.ws, PeaksWorkspace="peaks1", NRows=3, NCols=3, NBins=5, ThresholdIoverSigma=3.0, MinFracSize=0.5 + InputWorkspace=self.ws, PeaksWorkspace="peaks4", NRows=3, NCols=3, NBins=5, ThresholdIoverSigma=3.0, MinFracSize=0.5 ) self.assertEqual(out.getNumberPeaks(), 0) + def test_exec_VarOverMean(self): + out = FindSXPeaksConvolve( + InputWorkspace=self.ws, + PeaksWorkspace="peaks5", + NRows=5, + NCols=5, + Nbins=3, + PeakFindingStrategy="VarianceOverMean", + ThresholdVarianceOverMean=3.0, + MinFracSize=0.02, + ) + self._assert_found_correct_peaks(out, integrated=False) + if __name__ == "__main__": unittest.main() From a842637d64a3581420ca654a67e8b27169c2fd57 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 15 Apr 2024 14:41:07 +0100 Subject: [PATCH 05/10] Use FindSXPeaksConvolve as default in sxd and update system test --- .../tests/framework/SXDAnalysis.py | 6 +++--- .../SXD23767_found_peaks_convolve.nxs.md5 | 1 + scripts/Diffraction/single_crystal/sxd.py | 20 +++++++++++-------- 3 files changed, 16 insertions(+), 11 deletions(-) create mode 100644 Testing/SystemTests/tests/framework/reference/SXD23767_found_peaks_convolve.nxs.md5 diff --git a/Testing/SystemTests/tests/framework/SXDAnalysis.py b/Testing/SystemTests/tests/framework/SXDAnalysis.py index 7c9afd239c91..1fecc3e232c6 100644 --- a/Testing/SystemTests/tests/framework/SXDAnalysis.py +++ b/Testing/SystemTests/tests/framework/SXDAnalysis.py @@ -21,16 +21,16 @@ def cleanup(self): def runTest(self): ws = Load(Filename="SXD23767.raw", LoadMonitors="Exclude") - self.peaks = SXD.find_sx_peaks(ws, nstd=6) + self.peaks = SXD.find_sx_peaks(ws, ThresholdVarianceOverMean=1.5) FindUBUsingFFT(PeaksWorkspace=self.peaks, MinD=1, MaxD=10, Tolerance=0.15) SelectCellOfType(PeaksWorkspace=self.peaks, CellType="Cubic", Centering="F", Apply=True) OptimizeLatticeForCellType(PeaksWorkspace=self.peaks, CellType="Cubic", Apply=True) self.nindexed, *_ = IndexPeaks(PeaksWorkspace=self.peaks, Tolerance=0.1, CommonUBForAll=True) def validate(self): - self.assertEqual(214, self.nindexed) + self.assertEqual(171, self.nindexed) latt = SXD.retrieve(self.peaks).sample().getOrientedLattice() - a, alpha = 5.6541, 90 # published value for NaCl is a=6.6402 but the detector positions haven't been calibrated + a, alpha = 5.6533, 90 # published value for NaCl is a=6.6402 but the detector positions haven't been calibrated self.assertAlmostEqual(a, latt.a(), delta=1e-5) self.assertAlmostEqual(a, latt.b(), delta=1e-5) self.assertAlmostEqual(a, latt.c(), delta=1e-5) diff --git a/Testing/SystemTests/tests/framework/reference/SXD23767_found_peaks_convolve.nxs.md5 b/Testing/SystemTests/tests/framework/reference/SXD23767_found_peaks_convolve.nxs.md5 new file mode 100644 index 000000000000..443842a84cac --- /dev/null +++ b/Testing/SystemTests/tests/framework/reference/SXD23767_found_peaks_convolve.nxs.md5 @@ -0,0 +1 @@ +23bbde7e64f3d3fcb00984fd91e89e62 diff --git a/scripts/Diffraction/single_crystal/sxd.py b/scripts/Diffraction/single_crystal/sxd.py index 78a23ccc089a..28c7db83cb3c 100644 --- a/scripts/Diffraction/single_crystal/sxd.py +++ b/scripts/Diffraction/single_crystal/sxd.py @@ -336,16 +336,20 @@ def remove_peaks_on_detector_edge(peaks, nedge): mantid.DeleteTableRows(TableWorkspace=peaks, Rows=iremove, EnableLogging=False) @staticmethod - def find_sx_peaks(ws, bg=None, nstd=None, lambda_min=0.45, lambda_max=3, nbunch=3, out_peaks=None, **kwargs): + def find_sx_peaks(ws, out_peaks=None, **kwargs): wsname = SXD.retrieve(ws).name() - ws_rb = wsname + "_rb" - mantid.Rebunch(InputWorkspace=ws, OutputWorkspace=ws_rb, NBunch=nbunch, EnableLogging=False) - SXD.mask_detector_edges(ws_rb) - SXD.crop_ws(ws_rb, lambda_min, lambda_max, xunit="Wavelength") if out_peaks is None: - out_peaks = wsname + "_peaks" # need to do this so not "_rb_peaks" - BaseSX.find_sx_peaks(ws_rb, bg, nstd, out_peaks, **kwargs) - mantid.DeleteWorkspace(ws_rb, EnableLogging=False) + out_peaks = wsname + "_peaks" + default_kwargs = { + "NRows": 7, + "NCols": 7, + "GetNBinsFromBackToBackParams": True, + "NFWHM": 6, + "PeakFindingStrategy": "VarianceOverMean", + "ThresholdVarianceOverMean": 2, + } + kwargs = {**default_kwargs, **kwargs} # will overwrite default with provided if duplicate keys + mantid.FindSXPeaksConvolve(InputWorkspace=wsname, PeaksWorkspace=out_peaks, **kwargs) return out_peaks @staticmethod From 0e6c1862e9aa54a542e30a3d6efeb874d9ef68bc Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 15 Apr 2024 14:54:03 +0100 Subject: [PATCH 06/10] Update documentation --- .../algorithms/FindSXPeaksConvolve-v1.rst | 27 ++++++++++++++----- ...ISIS_SingleCrystalDiffraction_Workflow.rst | 2 +- 2 files changed, 21 insertions(+), 8 deletions(-) diff --git a/docs/source/algorithms/FindSXPeaksConvolve-v1.rst b/docs/source/algorithms/FindSXPeaksConvolve-v1.rst index 77485b607754..5d67d2889c02 100644 --- a/docs/source/algorithms/FindSXPeaksConvolve-v1.rst +++ b/docs/source/algorithms/FindSXPeaksConvolve-v1.rst @@ -10,8 +10,18 @@ Description ----------- This is an algorithm to find single-crystal Bragg peaks in a :ref:`MatrixWorkspace ` with detector -banks of type :ref:`RectangularDetector ` (e.g. SXD, TOPAZ) by integrating the data using a -convolution with a shoebox kernel. +banks of type :ref:`RectangularDetector ` (e.g. SXD, TOPAZ). + +There are two peak finding stratgeyies set by ``PeakFindingStrategy``: + +a. ``PeakFindingStrategy="IOverSigma"`` - by integrating the data by convolution with a shoebox kernel and looking for + regions with statistically significant I/sigma (larger than ``ThresholdIoverSigma``). Note a valid peak would be + expected to have intensity/sigma > 3 and stronger peaks will have a larger intensity/sigma. + +b. ``PeakFindingStrategy="VarianceOverMean"`` - by looking for regions with ratio of variance/mean larger than + ``ThresholdVarianceOverMean`` - note for a poisson distributed counts with a constant count-rate the ratio is + expected to be 1. This peak finding criterion taken from DIALS [1]_. + The size of the kernel is defined in the input to the algorithm and should match the approximate extent of a typical peak. The size on the detector is governed by ``NRows`` and ``NCols`` which are in units of pixels. @@ -26,15 +36,14 @@ The size of the kernel along the TOF dimension can be specified in one of two wa Note to use method 2, back-to-back exponential coefficients must be defined in the Parameters.xml file for the instrument. -The integration requires a background shell with negative weights, such that the total kernel size is increased by a -factor 1.25 along each dimension (such that there are approximately the same number of elements in the kernel and the -background shell). The integral of the kernel and background shell together is zero. +The shoebox integration (for ``PeakFindingStrategy="IOverSigma"``) requires a background shell with negative weights, +such that the total kernel size is increased by a factor 1.25 along each dimension (such that there are approximately +the same number of elements in the kernel and the background shell). The integral of the kernel and background shell +together is zero. The integrated intensity is evaluated by convolution of the signal array with the kernel and the square of the error on the integrated intensity is determined by convolution of the squared error array with the squared value of the kernel. -The threshold for peak detection is given by ``ThresholdIoverSigma`` which is the cutoff ratio of intensity/sigma (i.e. -a valid peak would be expected to have intensity/sigma > 3). Stronger peaks will have a larger intensity/sigma. Usage ----- @@ -56,6 +65,10 @@ Usage Found 261 peaks +References +---------- + +.. [1] Winter, G., et al. Acta Crystallographica Section D: Structural Biology 74.2 (2018): 85-97. .. categories:: diff --git a/docs/source/techniques/ISIS_SingleCrystalDiffraction_Workflow.rst b/docs/source/techniques/ISIS_SingleCrystalDiffraction_Workflow.rst index 95b272ffad8e..2f250d08ba3c 100644 --- a/docs/source/techniques/ISIS_SingleCrystalDiffraction_Workflow.rst +++ b/docs/source/techniques/ISIS_SingleCrystalDiffraction_Workflow.rst @@ -37,7 +37,7 @@ Here is an example using the ``SXD`` class that finds peaks and then removes dup # find peaks using SXD static method - determines peak threshold from # the standard deviation of the intensity distribution - peaks_ws = SXD.find_sx_peaks(ws, nstd=6) + peaks_ws = SXD.find_sx_peaks(ws, ThresholdVarianceOverMean=2.0) SXD.remove_duplicate_peaks_by_qlab(peaks_ws, q_tol=0.05) # find a UB From 396ca195e5a2b9adf18466da6ac78f78048b2629 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 15 Apr 2024 15:22:28 +0100 Subject: [PATCH 07/10] Add release notes --- .../v6.10.0/Diffraction/Single_Crystal/New_features/37171.rst | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 docs/source/release/v6.10.0/Diffraction/Single_Crystal/New_features/37171.rst diff --git a/docs/source/release/v6.10.0/Diffraction/Single_Crystal/New_features/37171.rst b/docs/source/release/v6.10.0/Diffraction/Single_Crystal/New_features/37171.rst new file mode 100644 index 000000000000..257fa714f384 --- /dev/null +++ b/docs/source/release/v6.10.0/Diffraction/Single_Crystal/New_features/37171.rst @@ -0,0 +1,2 @@ +- Add option to find peaks using the ratio of variance/ mean in :ref:`FindSXPeaksConvolve ` - this is a peak finding criterion used in DIALS software Winter, G., et al. Acta Crystallographica Section D: Structural Biology 74.2 (2018): 85-97. +- :ref:`FindSXPeaksConvolve ` is the default peak finding algorithm in the SXD reduction class. \ No newline at end of file From dfd02fb9c305b1cc3141b0ec347dc527bfeb6dd0 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Mon, 15 Apr 2024 15:26:00 +0100 Subject: [PATCH 08/10] Fix typo in docs --- docs/source/algorithms/FindSXPeaksConvolve-v1.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/algorithms/FindSXPeaksConvolve-v1.rst b/docs/source/algorithms/FindSXPeaksConvolve-v1.rst index 5d67d2889c02..c12064c61bbc 100644 --- a/docs/source/algorithms/FindSXPeaksConvolve-v1.rst +++ b/docs/source/algorithms/FindSXPeaksConvolve-v1.rst @@ -12,7 +12,7 @@ Description This is an algorithm to find single-crystal Bragg peaks in a :ref:`MatrixWorkspace ` with detector banks of type :ref:`RectangularDetector ` (e.g. SXD, TOPAZ). -There are two peak finding stratgeyies set by ``PeakFindingStrategy``: +There are two peak finding strategies set by ``PeakFindingStrategy``: a. ``PeakFindingStrategy="IOverSigma"`` - by integrating the data by convolution with a shoebox kernel and looking for regions with statistically significant I/sigma (larger than ``ThresholdIoverSigma``). Note a valid peak would be From 95608247336fee4e75f5f5f7ce0c8309d73fc4d9 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 1 May 2024 09:47:34 +0100 Subject: [PATCH 09/10] Remove unnecessary reference file to avoid flaky tests Test passes on windows not on jenkins linux, asserts test the functionality sufficiently --- Testing/SystemTests/tests/framework/SXDAnalysis.py | 1 - 1 file changed, 1 deletion(-) diff --git a/Testing/SystemTests/tests/framework/SXDAnalysis.py b/Testing/SystemTests/tests/framework/SXDAnalysis.py index 1fecc3e232c6..f165f93a9668 100644 --- a/Testing/SystemTests/tests/framework/SXDAnalysis.py +++ b/Testing/SystemTests/tests/framework/SXDAnalysis.py @@ -37,7 +37,6 @@ def validate(self): self.assertAlmostEqual(alpha, latt.alpha(), delta=1e-10) self.assertAlmostEqual(alpha, latt.beta(), delta=1e-10) self.assertAlmostEqual(alpha, latt.gamma(), delta=1e-10) - return self.peaks, "SXD23767_found_peaks.nxs" class SXDDetectorCalibration(systemtesting.MantidSystemTest): From f653bed9e142f01fb311c9b9bba34f0524283c95 Mon Sep 17 00:00:00 2001 From: Richard Waite Date: Wed, 1 May 2024 13:21:09 +0100 Subject: [PATCH 10/10] Update technique docs as per PR review --- .../techniques/ISIS_SingleCrystalDiffraction_Workflow.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/techniques/ISIS_SingleCrystalDiffraction_Workflow.rst b/docs/source/techniques/ISIS_SingleCrystalDiffraction_Workflow.rst index 2f250d08ba3c..ac6e60bb57f8 100644 --- a/docs/source/techniques/ISIS_SingleCrystalDiffraction_Workflow.rst +++ b/docs/source/techniques/ISIS_SingleCrystalDiffraction_Workflow.rst @@ -36,7 +36,7 @@ Here is an example using the ``SXD`` class that finds peaks and then removes dup ws = Load(Filename='SXD33335.nxs', OutputWorkspace='SXD33335') # find peaks using SXD static method - determines peak threshold from - # the standard deviation of the intensity distribution + # the ratio of local variance over mean peaks_ws = SXD.find_sx_peaks(ws, ThresholdVarianceOverMean=2.0) SXD.remove_duplicate_peaks_by_qlab(peaks_ws, q_tol=0.05)