diff --git a/Sources/CuboBragg/BraggDisks.swift b/Sources/CuboBragg/BraggDisks.swift index b2f257e..de18ff3 100644 --- a/Sources/CuboBragg/BraggDisks.swift +++ b/Sources/CuboBragg/BraggDisks.swift @@ -38,11 +38,11 @@ public struct BraggDetector: Sendable { singleBox = SingleBox(rows: detectorRows, columns: detectorColumns, kernel: kernel.values) } - /// Whether `BraggPrecision.single` can run on a detector of this size: vDSP's complex DFT must support both lengths - /// (f·2ⁿ with f in 1, 3, 5, 15, ...: 64, 128, 192, 240, 256, 320, 384, 480, 512, ...). Other sizes are refused by the - /// `.single` entry points with `CuboBraggError.invalidOption`; `.double` runs on any size. + /// Whether `BraggPrecision.single` can run on a detector of this size: every positive size. Lengths vDSP's complex DFT + /// takes (f·2ⁿ with f in 1, 3, 5, 15, ...: 64, 128, 192, 240, 256, 320, 384, 480, 512, ...) run on it directly; any + /// other runs through Bluestein's chirp-z on a longer vDSP DFT, at the same float32 tolerances and a higher cost. public static func supportsSinglePrecision(detectorRows: Int, detectorColumns: Int) -> Bool { - FloatFFT2D.supports(rows: detectorRows, columns: detectorColumns) + detectorRows > 0 && detectorColumns > 0 } /// The correlation of one pattern: the real map `CC` and the Fourier-domain product `cc` (see the file header). @@ -116,7 +116,7 @@ public struct BraggDetector: Sendable { return result } - /// The Float32 path's state for the whole-cube pass (refuses an unsupported size here, before any tile is read). + /// The Float32 path's state for the whole-cube pass (built here, before any tile is read). func singleState() throws -> SingleState { try singleBox.state() } func correlate(values input: [Double], options o: BraggDiskOptions) throws -> (Correlation, Int) { @@ -266,7 +266,7 @@ public enum BraggDisks { throw CuboBraggError.shapeMismatch("the cube's detector is \(shape.detectorRows) × \(shape.detectorColumns), the kernel \(kernel.detectorRows) × \(kernel.detectorColumns)") } let detector = try BraggDetector(kernel: kernel) - let single = options.precision == .single ? try detector.singleState() : nil // refuses an unsupported size before any read + let single = options.precision == .single ? try detector.singleState() : nil // the Float plan, built once before any read let n = shape.detectorCount let total = shape.scanCount var all: [[BraggPeak]] = [] diff --git a/Sources/CuboBragg/BraggTypes.swift b/Sources/CuboBragg/BraggTypes.swift index 69bf4b4..670c2a5 100644 --- a/Sources/CuboBragg/BraggTypes.swift +++ b/Sources/CuboBragg/BraggTypes.swift @@ -66,10 +66,10 @@ public enum BraggPrecision: Sendable, Hashable, CaseIterable { case double /// Float32 (Accelerate: vDSP convolution and complex DFT) for the pattern copy, the non-finite fill, both Gaussian /// filters, the 2-D transform, the correlation product and the maxima search; the few accepted peaks are refined in - /// Double. Available only for detector sizes vDSP's DFT supports on both axes (f·2ⁿ with f in 1, 3, 5, 15, ...: 64, 128, - /// 192, 240, 256, 320, 384, 480, 512, ...; `BraggDetector.supportsSinglePrecision(detectorRows:detectorColumns:)` asks); - /// any other size is refused with `CuboBraggError.invalidOption`, never computed another way. Agrees with `.double` - /// to about 1e-5 of the correlation maximum, hence peak positions to about 1e-4 px (docs/bragg.md has the measured gaps). + /// Double. Runs on every detector size: lengths vDSP's DFT takes (f·2ⁿ with f in 1, 3, 5, 15, ...: 64, 128, 192, 240, + /// 256, 320, 384, 480, 512, ...) directly, any other through Bluestein's chirp-z on a longer vDSP DFT (slower, same + /// tolerances). Agrees with `.double` to about 1e-5 of the correlation maximum, hence peak positions to about 1e-4 px + /// (docs/bragg.md has the measured gaps). case single } diff --git a/Sources/CuboBragg/CuboBragg.docc/CuboBragg.md b/Sources/CuboBragg/CuboBragg.docc/CuboBragg.md index 24f9dd1..9383a66 100644 --- a/Sources/CuboBragg/CuboBragg.docc/CuboBragg.md +++ b/Sources/CuboBragg/CuboBragg.docc/CuboBragg.md @@ -10,9 +10,8 @@ correlation, the filter cascade, and parabolic or matrix-DFT (multicorr) subpixe Opt in to the faster Float32 path with `options.precision = .single` (``BraggPrecision``): vDSP convolution and complex DFT, a reusable scratch per concurrent worker, about ten times faster than Double at 256 × 256 and equal to it to about 1e-5 of the -correlation maximum. It needs a detector size vDSP supports (64, 128, 192, 240, 256, ...; ask with -``BraggDetector/supportsSinglePrecision(detectorRows:detectorColumns:)``) and is refused with ``CuboBraggError/invalidOption(_:)`` -otherwise. ``BraggDetector/detectWithFunnel(_:options:)`` returns the peaks of `detect` with a ``BraggFunnel``, the count of +correlation maximum. It runs on every detector size: sizes vDSP's DFT takes (64, 128, 192, 240, 256, ...) directly, any other +through Bluestein's chirp-z, slower and at the same tolerances. ``BraggDetector/detectWithFunnel(_:options:)`` returns the peaks of `detect` with a ``BraggFunnel``, the count of maxima after each stage of the selection, so a short or empty result can be explained; ``BraggDetector/fillNonFinite(_:rows:columns:)`` is the median fill of non-finite pixels, public. diff --git a/Sources/CuboBragg/CuboBragg.docc/Finding-Bragg-disks.md b/Sources/CuboBragg/CuboBragg.docc/Finding-Bragg-disks.md index c3ba824..54cf686 100644 --- a/Sources/CuboBragg/CuboBragg.docc/Finding-Bragg-disks.md +++ b/Sources/CuboBragg/CuboBragg.docc/Finding-Bragg-disks.md @@ -28,8 +28,8 @@ let peaks = vectors[0, 0] // the ``` ``BraggDiskOptions`` defaults are py4DSTEM's. Choose the subpixel method with ``BraggSubpixel`` and the Float32 path with -`options.precision = .single` (``BraggPrecision``). Float32 needs a detector size vDSP supports; ask first with -``BraggDetector/supportsSinglePrecision(detectorRows:detectorColumns:)``. +`options.precision = .single` (``BraggPrecision``). Float32 runs on every detector size; sizes vDSP's DFT takes +(64, 128, 192, 240, 256, ...) are the fastest. ## When a result is short or empty diff --git a/Sources/CuboBragg/SingleDetection.swift b/Sources/CuboBragg/SingleDetection.swift index 23b019c..6ba9f34 100644 --- a/Sources/CuboBragg/SingleDetection.swift +++ b/Sources/CuboBragg/SingleDetection.swift @@ -80,9 +80,9 @@ final class SingleState: @unchecked Sendable { private var free: [SingleScratch] = [] init(rows: Int, columns: Int, kernel: [Double]) throws { - guard let fft = FloatFFT2D(rows: rows, columns: columns) else { - throw CuboBraggError.invalidOption("precision .single needs a detector size vDSP's DFT supports on both axes (f·2ⁿ, f in 1, 3, 5, 15, ...); \(rows) × \(columns) is not one (use .double)") - } + guard rows > 0, columns > 0 else { throw CuboBraggError.shapeMismatch("the kernel is empty") } + // vDSP's DFT where it takes the length (the same calls and bits as `FloatFFT2D(rows:columns:)`), Bluestein otherwise. + let fft = FloatFFT2D(anyRows: rows, columns: columns) self.rows = rows; self.columns = columns; self.fft = fft let n = rows * columns let scratch = SingleScratch(count: n) diff --git a/Tests/CuboBraggTests/SinglePrecisionTests.swift b/Tests/CuboBraggTests/SinglePrecisionTests.swift index f559858..bccea51 100644 --- a/Tests/CuboBraggTests/SinglePrecisionTests.swift +++ b/Tests/CuboBraggTests/SinglePrecisionTests.swift @@ -12,7 +12,8 @@ import Testing func single(_ o: BraggDiskOptions) -> BraggDiskOptions { var s = o; s.precision = .single; return s } - /// Scenes vDSP's DFT supports (64 × 64); S2 (40 × 56) is refused. + /// Every scene runs in `.single`: the 64 × 64 ones on vDSP's DFT directly, S2 (40 × 56: vDSP has no DFT of length 56) + /// through Bluestein (CuboKit#60). func supported(_ c: [String: Any]) -> Bool { let s = fixture.scene(c["scene"] as! String) return BraggDetector.supportsSinglePrecision(detectorRows: s.rows, detectorColumns: s.columns) @@ -34,7 +35,8 @@ import Testing strongPosition = max(strongPosition, g.strongPosition); strongIntensity = max(strongIntensity, g.strongIntensity); checked += 1 } print("CuboBragg .single vs py4DSTEM: \(checked) cases, maximum position gap \(worstPosition) px (peaks above 1e-3 of the brightest: \(strongPosition)), maximum relative intensity gap \(worstIntensity) (those: \(strongIntensity))") - #expect(checked >= 24) + #expect(checked >= 29) // S2 (40 × 56, Bluestein) included since CuboKit#60 + } /// The same inputs through `.double` and `.single`: peaks agree (py4DSTEM's cases and the app's variants). @@ -54,7 +56,7 @@ import Testing strongPosition = max(strongPosition, g.strongPosition); strongIntensity = max(strongIntensity, g.strongIntensity); checked += 1 } print("CuboBragg .single vs .double: \(checked) cases, maximum position gap \(worstPosition) px (peaks above 1e-3 of the brightest: \(strongPosition)), maximum relative intensity gap \(worstIntensity) (those: \(strongIntensity))") - #expect(checked >= 30) + #expect(checked >= 39) // S2 (40 × 56, Bluestein) included since CuboKit#60 } /// The correlation itself: the real map and the Fourier product, as a fraction of their largest magnitude. @@ -80,35 +82,86 @@ import Testing #expect(worstFourier <= SingleTolerance.fourier, "Fourier gap \(worstFourier)") } - @Test func refusesUnsupportedSizesAndKeepsDouble() async throws { - let s2 = fixture.scene("S2") // 40 × 56: vDSP has no DFT of length 56 - #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 40, detectorColumns: 56)) - #expect(BraggDetector.supportsSinglePrecision(detectorRows: 64, detectorColumns: 64)) - #expect(BraggDetector.supportsSinglePrecision(detectorRows: 256, detectorColumns: 256)) - #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 56, detectorColumns: 56)) - #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 127, detectorColumns: 128)) - #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 0, detectorColumns: 64)) - var o = BraggDiskOptions(); o.precision = .single - #expect(throws: CuboBraggError.self) { try BraggDisks.detect(s2.patterns[0], kernel: s2.kernel, options: o) } - #expect(throws: CuboBraggError.self) { try BraggDetector(kernel: s2.kernel).detectWithFunnel(s2.patterns[0], options: o) } - #expect(throws: CuboBraggError.self) { try BraggDetector(kernel: s2.kernel).correlate(s2.patterns[0], options: o) } - do { _ = try BraggDisks.detect(s2.patterns[0], kernel: s2.kernel, options: o) } catch let e as CuboBraggError { - guard case .invalidOption(let message) = e else { Issue.record("wrong error \(e)"); return } - #expect(message.contains("40 × 56")) + /// CuboKit#60: `.single` runs on every positive detector size (vDSP's DFT where it takes the length, Bluestein otherwise); + /// `.double` stays the default and options are validated first, as in `.double`. + @Test func everySizeIsSupportedAndDoubleStaysDefault() throws { + for (rows, columns) in [(64, 64), (256, 256), (40, 56), (56, 56), (127, 128), (22, 21), (165, 40), (1, 1)] { + #expect(BraggDetector.supportsSinglePrecision(detectorRows: rows, detectorColumns: columns), "\(rows) × \(columns)") } - // The whole-cube pass refuses before reading a tile. - let shape = try CubeShape(scanRows: 1, scanColumns: 1, detectorRows: 40, detectorColumns: 56) - let cube = try MemoryCube(shape: shape, values: [Float](repeating: 0, count: 40 * 56)) - await #expect(throws: CuboBraggError.self) { _ = try await BraggDisks.detectAll(in: cube, kernel: s2.kernel, options: o) } - // .double is untouched on that size, and is the default. + #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 0, detectorColumns: 64)) + #expect(!BraggDetector.supportsSinglePrecision(detectorRows: 64, detectorColumns: -1)) #expect(BraggDiskOptions().precision == .double) - _ = try BraggDisks.detect(s2.patterns[0], kernel: s2.kernel) - // Options are validated first, as in .double. - o.upsampleFactor = 2 + var o = BraggDiskOptions(); o.precision = .single; o.upsampleFactor = 2 #expect(throws: CuboBraggError.self) { try BraggDisks.detect(fixture.scene("S1").patterns[0], kernel: fixture.scene("S1").kernel, options: o) } } - /// Non-square supported grids (the transposed spectrum must be undone correctly) and a Gaussian wider than the image + /// CuboKit#60: detectors vDSP's DFT cannot take on one or both axes (22 × 21, 165 × 40, S2's 40 × 56): the `.single` + /// correlation map and Fourier product agree with `.double` at the written float32 tolerances, the peaks of the S1 data + /// tiled onto the grid agree within a stated bound (below), and the whole-cube pass equals the per-pattern pass. + @Test func nonVDSPSizesAgreeWithDouble() async throws { + let s = fixture.scene("S1") + var worstMap = 0.0, worstFourier = 0.0, worstPosition = 0.0, worstIntensity = 0.0, peaksChecked = 0 + for (rows, columns, sigma) in [(22, 21, 0.0), (22, 21, 1.5), (165, 40, 0.0), (165, 40, 2.0), (40, 56, 0.0)] { + #expect(!FloatFFT2D.supports(rows: rows, columns: columns), "\(rows) × \(columns) should need Bluestein") + let kernel = DetectorImage(detectorRows: rows, detectorColumns: columns, + values: (0..<(rows * columns)).map { s.kernel.values[$0 % 4096] }, path: .cpu) + let detector = try BraggDetector(kernel: kernel) + var o = BraggDiskOptions(); o.patternSigma = sigma + var os = o; os.precision = .single + var floats: [Float] = [] + for pattern in s.patterns { + let values = (0..<(rows * columns)).map { Double(Float(pattern.values[$0 % 4096])) } + floats += values.map { Float($0) } + let image = DetectorImage(detectorRows: rows, detectorColumns: columns, values: values, path: .cpu) + let d = try detector.correlate(image, options: o) + let f = try detector.correlate(image, options: os) + let mapScale = d.real.max()! + let fourierScale = zip(d.fourierRe, d.fourierIm).map { hypot($0, $1) }.max()! + for i in d.real.indices { + worstMap = max(worstMap, abs(f.real[i] - d.real[i]) / mapScale) + worstFourier = max(worstFourier, hypot(f.fourierRe[i] - d.fourierRe[i], f.fourierIm[i] - d.fourierIm[i]) / fourierScale) + } + for subpixel in [BraggSubpixel.parabolic, .multicorr] { + var p = BraggDiskOptions.mac4DSTEM; p.patternSigma = sigma; p.subpixel = subpixel + p.edgeBoundary = 2; p.minimumPeakSpacing = 4 + var ps = p; ps.precision = .single + let dp = try detector.detect(image, options: p) + let fp = try detector.detect(image, options: ps) + let id = "\(rows) × \(columns) σ \(sigma) \(subpixel)" + #expect(dp.count == fp.count, "\(id): double \(dp.count) peaks, single \(fp.count)") + guard dp.count == fp.count else { continue } + // The tiled scene has broad, low-contrast maxima; there even vDSP's own sizes (160 × 40) exceed the strict + // per-peak tier (measured 8.4e-6 relative intensity), so these peaks get a stated bound of 2e-4 px (the + // app's own tolerance against py4DSTEM) and 2e-5 (measured here 9.9e-5 px, 8.0e-6). The real non-vDSP + // scene S2 (40 × 56) is held to the strict tier in `singleAgreesWithDouble`. + for (a, e) in zip(fp, dp) { + let position = max(abs(a.row - e.row), abs(a.column - e.column)) + let intensity = abs(a.intensity - e.intensity) / max(abs(e.intensity), 1e-300) + #expect(position <= 2e-4 && intensity <= 2e-5, + "\(id): peak at \(e.row), \(e.column): position gap \(position), intensity gap \(intensity)") + worstPosition = max(worstPosition, position); worstIntensity = max(worstIntensity, intensity) + } + peaksChecked += dp.count + } + } + // The whole-cube pass in `.single` on this size is the per-pattern pass. + let count = s.patterns.count + let shape = try CubeShape(scanRows: count, scanColumns: 1, detectorRows: rows, detectorColumns: columns) + let cube = try MemoryCube(shape: shape, values: floats) + var p = BraggDiskOptions.mac4DSTEM; p.patternSigma = sigma; p.precision = .single + let all = try await BraggDisks.detectAll(in: cube, kernel: kernel, options: p, tiling: .rows(2)) + let n = rows * columns + for k in 0.. 0) + } + + /// Non-square vDSP grids (the transposed spectrum must be undone correctly) and a Gaussian wider than the image /// (the reflection folds repeatedly): the correlation map and its Fourier product agree with `.double`. @Test func nonSquareAndWideFilter() throws { let s = fixture.scene("S1") diff --git a/docs/bragg.md b/docs/bragg.md index 76970a6..b47ae97 100644 --- a/docs/bragg.md +++ b/docs/bragg.md @@ -53,10 +53,12 @@ hybrid power in `powf`), `max(Re, 0)` and the maxima search run in Float. The fe from the Float neighbourhood, the matrix-DFT upsampling on the Float complex correlation widened to Double. The reported intensity is the Float correlation, widened. -- **Sizes.** The transform is the exact N-point DFT at the pattern's own size, so only detector sizes vDSP supports on both axes are - available: f·2ⁿ with f in 1, 3, 5, 15, ... (64, 128, 192, 240, 256, 320, 384, 480, 512, ...; not 56, 100, 127, 130). - `BraggDetector.supportsSinglePrecision(detectorRows:detectorColumns:)` asks vDSP itself. Any other size is refused with - `CuboBraggError.invalidOption` (naming the size), never computed another way; `.double` runs on every size. +- **Sizes.** The transform is the exact N-point DFT at the pattern's own size, on every detector size (CuboCore's + `FloatFFT2D(anyRows:columns:)`). A length vDSP's DFT takes (f·2ⁿ with f in 1, 3, 5, 15, ...: 64, 128, 192, 240, 256, 320, + 384, 480, 512, ...) runs on it directly, with the same calls and bits as before; any other (56, 100, 127, 130, 165, ...) + runs as Bluestein's chirp-z through a vDSP DFT of the shortest length M ≥ 2N − 1 vDSP takes, slower but held to the same + tolerances (`SinglePrecisionTests.nonVDSPSizesAgreeWithDouble`: 22 × 21, 165 × 40, 40 × 56). + `BraggDetector.supportsSinglePrecision(detectorRows:detectorColumns:)` is true for every positive size. - **Scratch.** A detector keeps a pool of scratch buffers (six `rows × columns` Float planes); a call checks one out and returns it. `BraggDisks.detectAll` therefore allocates one per concurrent worker, not one per pattern, and with `.single` hands each Float pattern straight from the tile to the detector (no widening). The Float plan and kernel spectrum are built on the first diff --git a/tools/consumer-check/Sources/consumer-check/main.swift b/tools/consumer-check/Sources/consumer-check/main.swift index 830b460..22a135f 100644 --- a/tools/consumer-check/Sources/consumer-check/main.swift +++ b/tools/consumer-check/Sources/consumer-check/main.swift @@ -121,7 +121,8 @@ do { && BraggPrecision.allCases == [.double, .single] && singlePeaks.count == 2 && abs(singlePeaks[0].row - doublePeaks[0].row) < 1e-3 && abs(singlePeaks[1].column - doublePeaks[1].column) < 1e-3) check("BraggDetector.supportsSinglePrecision", BraggDetector.supportsSinglePrecision(detectorRows: 48, detectorColumns: 48) - && !BraggDetector.supportsSinglePrecision(detectorRows: 56, detectorColumns: 40)) + && BraggDetector.supportsSinglePrecision(detectorRows: 56, detectorColumns: 40) + && !BraggDetector.supportsSinglePrecision(detectorRows: 0, detectorColumns: 40)) let plainPeaks = try detector.detect(image, options: braggOptions) check("BraggDetector.detectWithFunnel(DetectorImage)", doublePeaks == plainPeaks && doubleFunnel.acceptedCount == 2 && doubleFunnel.localMaximumCount >= 2 && singleFunnel.acceptedCount == 2)