Add GPU resident radius GICP covariance preparation

CHANGELOG.md

@@ -21,6 +21,15 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`).
 
 ### Changed
 
+- Add explicit resident radius covariance preparation for GICP under
+  `register-gicp-gpu`. Inclusive bounded GPU neighborhoods, centered/scaled
+  compensated central sums and plane eigenvectors feed resident equations
+  without covariance readback/re-upload. Caller capacity includes self and
+  duplicates; overflow rejects all output without truncation. Sparse rows
+  preserve epsilon*I. Radius must be positive finite with normal f32 square,
+  epsilon=1e-4..=1. CPU trees/6x6 solves remain; legacy routes retain behavior.
+  Device f32 accumulation and Jacobi iteration can differ from grid results.
+
 - Add explicit resident k-NN covariance preparation for GICP under
   `register-gicp-gpu`. GPU bounded neighborhoods, centered/scaled sample
   covariance and plane eigenvectors feed the existing resident equation loop

README.md

@@ -672,6 +672,15 @@ readback or re-upload occurs. This path requires covariance_radius=None,
 k=3..=64 and epsilon=1e-4..=1; CPU builds trees and solves 6x6 equations.
 See [resident covariance measurements](notes/2026-10-13_gpu_gicp_resident_covariances.md).
 
+`GicpRegistration::align_with_gpu_resident_radius_covariances(&runtime, &source,
+&target, neighbor_capacity)` prepares radius neighborhoods and covariances
+without CPU readback/re-upload. Set `covariance_radius=Some(radius)` and choose
+a positive capacity including self/duplicates; overflow rejects the result.
+Radius must have a finite normal f32 square, epsilon must be in 1e-4..=1.
+Sparse neighborhoods retain the grid route's epsilon*I fallback. CPU tree
+construction and 6x6 solves remain; only 256-byte equations are downloaded.
+See [resident radius evidence](notes/2026-10-13_gpu_gicp_resident_radius_covariances.md).
+
 Point-to-plane ICP also exposes
 `PointToPlaneIcp::align_with_gpu_reduction(&runtime, &source, &target)` under
 `register-icp-point-to-plane-gpu`. Transforms, nearest matching and plane

crates/spatialrust-gpu/src/kernels/gicp_covariance_tests.rs

@@ -89,7 +89,7 @@ fn matrix(c: [f32; 6]) -> Mat3<f64> {
 
 #[test]
 fn resident_covariance_plans_bound_counts_arity_epsilon_and_dispatch() {
-    assert_eq!(std::mem::size_of::<Params>(), 16);
+    assert_eq!(std::mem::size_of::<Params>(), 32);
     let limits = wgpu::Limits::default();
     assert_eq!(
         covariance_plan(65537, 65554, 20, 0.001, &limits).unwrap(),
@@ -296,76 +296,106 @@ fn malformed_resident_neighbor_rows_invalidate_both_covariance_sets() {
     let source = [[0.0; 3], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]];
     let st = tree(&source);
     let tt = tree(&source[..2]);
-    for second in [false, true] {
-        for case in 0..3 {
-            let positions = GpuIcpInput::upload(&runtime, &source, &tt).unwrap();
-            let source_tree = GpuKdTree::upload(&runtime, &st).unwrap();
-            let sn = source_tree.nearest_k_resident(&runtime, &positions.source, 3).unwrap();
-            let tn = positions.tree.nearest_k_resident(&runtime, &positions.target, 3).unwrap();
-            let statement = match case {
-                0 => "data[0]=0u;",
-                1 => "data[2]=0xffffffffu;",
-                _ => "data[3]=0x7f800000u;",
-            };
-            let shader=runtime.device().create_shader_module(wgpu::ShaderModuleDescriptor{label:None,source:wgpu::ShaderSource::Wgsl(format!("@group(0) @binding(0) var<storage,read_write> data:array<u32>; @compute @workgroup_size(1) fn main() {{{statement}}}").into())});
-            let pipeline =
-                runtime.device().create_compute_pipeline(&wgpu::ComputePipelineDescriptor {
+    for radius_mode in [0, 1] {
+        for second in [false, true] {
+            for case in 0..if radius_mode == 0 { 3 } else { 6 } {
+                let positions = GpuIcpInput::upload(&runtime, &source, &tt).unwrap();
+                let source_tree = GpuKdTree::upload(&runtime, &st).unwrap();
+                let sn = if radius_mode == 0 {
+                    source_tree.nearest_k_resident(&runtime, &positions.source, 3).unwrap()
+                } else {
+                    source_tree
+                        .radius_neighbors_resident(&runtime, &positions.source, 2.0, 3)
+                        .unwrap()
+                };
+                let tn = if radius_mode == 0 {
+                    positions.tree.nearest_k_resident(&runtime, &positions.target, 3).unwrap()
+                } else {
+                    positions
+                        .tree
+                        .radius_neighbors_resident(&runtime, &positions.target, 2.0, 3)
+                        .unwrap()
+                };
+                let statement = match case {
+                    0 => "data[0]=0u;",
+                    1 => "data[2]=0xffffffffu;",
+                    2 => "data[3]=0x7f800000u;",
+                    3 => "data[3]=0x40a00000u;",
+                    4 => "data[1]=0u;",
+                    _ => "data[1]=3u;data[0]=3u;data[2]=0xffffffffu;",
+                };
+                let shader=runtime.device().create_shader_module(wgpu::ShaderModuleDescriptor{label:None,source:wgpu::ShaderSource::Wgsl(format!("@group(0) @binding(0) var<storage,read_write> data:array<u32>; @compute @workgroup_size(1) fn main() {{{statement}}}").into())});
+                let pipeline =
+                    runtime.device().create_compute_pipeline(&wgpu::ComputePipelineDescriptor {
+                        label: None,
+                        layout: None,
+                        module: &shader,
+                        entry_point: Some("main"),
+                        compilation_options: Default::default(),
+                        cache: None,
+                    });
+                let buffer = if second { tn.buffer() } else { sn.buffer() };
+                let group = runtime.device().create_bind_group(&wgpu::BindGroupDescriptor {
                     label: None,
-                    layout: None,
-                    module: &shader,
-                    entry_point: Some("main"),
-                    compilation_options: Default::default(),
-                    cache: None,
+                    layout: &pipeline.get_bind_group_layout(0),
+                    entries: &[wgpu::BindGroupEntry {
+                        binding: 0,
+                        resource: buffer.as_entire_binding(),
+                    }],
                 });
-            let buffer = if second { tn.buffer() } else { sn.buffer() };
-            let group = runtime.device().create_bind_group(&wgpu::BindGroupDescriptor {
-                label: None,
-                layout: &pipeline.get_bind_group_layout(0),
-                entries: &[wgpu::BindGroupEntry {
-                    binding: 0,
-                    resource: buffer.as_entire_binding(),
-                }],
-            });
-            let mut encoder = runtime.device().create_command_encoder(&Default::default());
-            {
-                let mut pass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor {
+                let mut encoder = runtime.device().create_command_encoder(&Default::default());
+                {
+                    let mut pass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor {
+                        label: None,
+                        timestamp_writes: None,
+                    });
+                    pass.set_pipeline(&pipeline);
+                    pass.set_bind_group(0, &group, &[]);
+                    pass.dispatch_workgroups(1, 1, 1);
+                }
+                runtime.queue().submit(Some(encoder.finish()));
+                let covariances = runtime.device().create_buffer(&wgpu::BufferDescriptor {
                     label: None,
-                    timestamp_writes: None,
+                    size: 6 * 32,
+                    usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
+                    mapped_at_creation: false,
                 });
-                pass.set_pipeline(&pipeline);
-                pass.set_bind_group(0, &group, &[]);
-                pass.dispatch_workgroups(1, 1, 1);
+                encode_covariances(
+                    &runtime,
+                    &positions,
+                    [&sn, &tn],
+                    &covariances,
+                    Params {
+                        sources: 3,
+                        targets: 2,
+                        neighbors: 3,
+                        epsilon: 0.001,
+                        radius_squared: 4.0,
+                        radius_mode,
+                        padding: [0; 2],
+                    },
+                );
+                let cov = words(&runtime, &covariances);
+                assert_eq!(cov[0], 1);
+                assert!(cov[1..].iter().all(|&v| v == 0));
+                let transfers = positions.upload_transfers();
+                let input = GpuGicpInput {
+                    positions,
+                    covariances,
+                    weight_scale: 6000.0,
+                    covariance_offset: 1,
+                    transfers,
+                };
+                let output =
+                    input.equations_resident(&runtime, Isometry3::identity(), 2.0).unwrap();
+                let w = words(&runtime, output.buffer());
+                assert_eq!(w[0], 16);
+                assert!(w[1..].iter().all(|&v| v == 0));
+                assert!(output.download(&runtime).is_err());
             }
-            runtime.queue().submit(Some(encoder.finish()));
-            let covariances = runtime.device().create_buffer(&wgpu::BufferDescriptor {
-                label: None,
-                size: 6 * 32,
-                usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
-                mapped_at_creation: false,
-            });
-            encode_covariances(
-                &runtime,
-                &positions,
-                [&sn, &tn],
-                &covariances,
-                Params { sources: 3, targets: 2, neighbors: 3, epsilon: 0.001 },
-            );
-            let cov = words(&runtime, &covariances);
-            assert_eq!(cov[0], 1);
-            assert!(cov[1..].iter().all(|&v| v == 0));
-            let transfers = positions.upload_transfers();
-            let input = GpuGicpInput {
-                positions,
-                covariances,
-                weight_scale: 6000.0,
-                covariance_offset: 1,
-                transfers,
-            };
-            let output = input.equations_resident(&runtime, Isometry3::identity(), 2.0).unwrap();
-            let w = words(&runtime, output.buffer());
-            assert_eq!(w[0], 16);
-            assert!(w[1..].iter().all(|&v| v == 0));
-            assert!(output.download(&runtime).is_err());
         }
     }
 }
+
+#[path = "gicp_radius_covariance_tests.rs"]
+mod radius;

crates/spatialrust-gpu/src/kernels/gicp_covariances.rs

@@ -12,6 +12,15 @@ struct Params {
     targets: u32,
     neighbors: u32,
     epsilon: f32,
+    radius_squared: f32,
+    radius_mode: u32,
+    padding: [u32; 2],
+}
+
+#[derive(Clone, Copy)]
+enum Preparation {
+    Knn(u32),
+    Radius { radius: f32, capacity: u32 },
 }
 
 fn covariance_plan(
@@ -21,17 +30,26 @@ fn covariance_plan(
     epsilon: f32,
     limits: &wgpu::Limits,
 ) -> SpatialResult<u64> {
-    if source == 0
-        || target == 0
-        || !(3..=MAX_GPU_K_NEIGHBORS).contains(&k)
-        || !(1e-4..=1.0).contains(&epsilon)
-    {
+    if !(3..=MAX_GPU_K_NEIGHBORS).contains(&k) {
+        return Err(invalid());
+    }
+    bounded_covariance_plan(source, target, k, epsilon, limits)
+}
+
+fn bounded_covariance_plan(
+    source: u64,
+    target: u64,
+    capacity: u32,
+    epsilon: f32,
+    limits: &wgpu::Limits,
+) -> SpatialResult<u64> {
+    if source == 0 || target == 0 || capacity == 0 || !(1e-4..=1.0).contains(&epsilon) {
         return Err(invalid());
     }
     let total = source.checked_add(target).and_then(|v| v.checked_add(1)).ok_or_else(invalid)?;
     plan(source, total, limits)?;
     let bound = limits.max_buffer_size.min(u64::from(limits.max_storage_buffer_binding_size));
-    let stride = 8 + 8 * u64::from(k);
+    let stride = 8 + 8 * u64::from(capacity);
     for count in [source, target] {
         if count > u64::from(u32::MAX)
             || count.div_ceil(64) > u64::from(limits.max_compute_workgroups_per_dimension)
@@ -50,6 +68,10 @@ fn covariance_plan(
     Ok(total * 32)
 }
 
+fn valid_radius(radius: f32) -> bool {
+    radius.is_finite() && radius > 0.0 && (radius * radius).is_normal()
+}
+
 impl GpuGicpInput {
     /// Explicitly uploads finite source/target trees and original-order positions,
     /// then computes both k-NN plane-regularized covariances entirely on GPU.
@@ -70,16 +92,58 @@ impl GpuGicpInput {
         target: &KdTree,
         k: u32,
         epsilon: f32,
+    ) -> SpatialResult<Self> {
+        Self::upload_prepared(runtime, source, target, Preparation::Knn(k), epsilon)
+    }
+
+    /// Uploads original positions/trees and prepares complete inclusive radius
+    /// neighborhoods and plane covariances on GPU, without payload readback.
+    /// `capacity` is a positive per-point bound, including self and duplicates.
+    /// Overflow invalidates both covariance sets and every equation download;
+    /// rows are never silently truncated. Radius must be finite/positive with
+    /// a normal f32 square; epsilon must be in 1e-4..=1.
+    /// Fewer than three neighbors use epsilon*I, as in the existing grid route.
+    /// Other rows use centered/scaled compensated sums and 32 max-pivot Jacobi
+    /// steps with physical 1e-20 stopping on the unnormalized central sum.
+    /// f32 arithmetic, traversal order and repeated eigenvalues can differ
+    /// from the legacy grid's raw means and 16 three-pivot sweeps.
+    /// Transfers, ownership and CPU solve boundaries match the k-NN method.
+    pub fn upload_with_gpu_radius_covariances(
+        runtime: &WgpuRuntime,
+        source: &KdTree,
+        target: &KdTree,
+        radius: f32,
+        capacity: u32,
+        epsilon: f32,
+    ) -> SpatialResult<Self> {
+        if !valid_radius(radius) {
+            return Err(invalid());
+        }
+        Self::upload_prepared(
+            runtime,
+            source,
+            target,
+            Preparation::Radius { radius, capacity },
+            epsilon,
+        )
+    }
+
+    fn upload_prepared(
+        runtime: &WgpuRuntime,
+        source: &KdTree,
+        target: &KdTree,
+        preparation: Preparation,
+        epsilon: f32,
     ) -> SpatialResult<Self> {
         let (x, y, z) = source.layout().positions();
         let targets = target.layout().positions().0.len();
-        let bytes = covariance_plan(
-            x.len() as u64,
-            targets as u64,
-            k,
-            epsilon,
-            &runtime.device().limits(),
-        )?;
+        let (capacity, radius) = match preparation {
+            Preparation::Knn(k) => (k, None),
+            Preparation::Radius { radius, capacity } => (capacity, Some(radius)),
+        };
+        let planner = if radius.is_some() { bounded_covariance_plan } else { covariance_plan };
+        let bytes =
+            planner(x.len() as u64, targets as u64, capacity, epsilon, &runtime.device().limits())?;
         let bound = runtime
             .device()
             .limits()
@@ -93,16 +157,42 @@ impl GpuGicpInput {
         let points: Vec<_> = (0..x.len()).map(|i| [x[i], y[i], z[i]]).collect();
         let positions = GpuIcpInput::upload(runtime, &points, target)?;
         let source_tree = GpuKdTree::upload(runtime, source)?;
-        let source_neighbors = source_tree.nearest_k_resident(runtime, &positions.source, k)?;
-        let target_neighbors = positions.tree.nearest_k_resident(runtime, &positions.target, k)?;
+        let (source_neighbors, target_neighbors) = if let Some(radius) = radius {
+            (
+                source_tree.radius_neighbors_resident(
+                    runtime,
+                    &positions.source,
+                    radius,
+                    capacity,
+                )?,
+                positions.tree.radius_neighbors_resident(
+                    runtime,
+                    &positions.target,
+                    radius,
+                    capacity,
+                )?,
+            )
+        } else {
+            (
+                source_tree.nearest_k_resident(runtime, &positions.source, capacity)?,
+                positions.tree.nearest_k_resident(runtime, &positions.target, capacity)?,
+            )
+        };
         let covariances = runtime.device().create_buffer(&wgpu::BufferDescriptor {
             label: Some("gicp-resident-covariances"),
             size: bytes,
             usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC,
             mapped_at_creation: false,
         });
-        let params =
-            Params { sources: x.len() as u32, targets: targets as u32, neighbors: k, epsilon };
+        let params = Params {
+            sources: x.len() as u32,
+            targets: targets as u32,
+            neighbors: capacity,
+            epsilon,
+            radius_squared: radius.map_or(0.0, |r| r * r),
+            radius_mode: u32::from(radius.is_some()),
+            padding: [0; 2],
+        };
         encode_covariances(
             runtime,
             &positions,