diff -r afb3ce40cff6 -r 66b8ac0ebe00 src/sliding_fb.rs
--- a/src/sliding_fb.rs	Thu Jul 30 09:27:44 2026 -0500
+++ b/src/sliding_fb.rs	Fri Jul 31 10:41:47 2026 -0500
@@ -28,6 +28,7 @@
 use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping, Mapping, RealMapping};
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use alg_tools::norms::Norm;
+use alg_tools::sets::{Set, Cube};
 use anyhow::ensure;
 use std::ops::ControlFlow;
 
@@ -230,6 +231,7 @@
     Domain: Space + Clone,
 {
     type TransportStepLength;
+    type Dom : Set<Domain>;
 
     /// Constrution of initial transport `γ1` from initial measure `μ` and `v=F'(μ)`
     /// with step lengh τ and transport step length `θ_or_adaptive`.
@@ -242,6 +244,7 @@
         τθ_or_adaptive: &mut Self::TransportStepLength,
         v: &PreadjointCodomain,
         tconfig: &TransportConfig<F>,
+        domain: Option<&Self::Dom>,
     );
 
     /// A posteriori transport adaptation.
@@ -281,6 +284,7 @@
     RNDM<N, F>: SpikeMerging<F>,
 {
     type TransportStepLength = TransportStepLength<F>;
+    type Dom = Cube<N, F>;
 
     fn get_transport_steplength(
         &self,
@@ -306,8 +310,9 @@
         τθ_or_adaptive: &mut TransportStepLength<F>,
         v: &M,
         tconfig: &TransportConfig<F>,
+        domain: Option<&Cube<N, F>>,
     ) {
-        γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig);
+        γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig, domain);
     }
 
     fn aposteriori_transport(
@@ -387,6 +392,7 @@
     RNDM<N, F>: SpikeMerging<F>,
 {
     type TransportStepLength = TransportStepLength<F>;
+    type Dom = Cube<N, F>;
 
     fn get_transport_steplength(
         &self,
@@ -412,8 +418,9 @@
         τθ_or_adaptive: &mut TransportStepLength<F>,
         v: &M,
         tconfig: &TransportConfig<F>,
+        domain: Option<&Cube<N, F>>,
     ) {
-        γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig);
+        γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig, domain);
     }
 
     fn aposteriori_transport(
@@ -744,6 +751,7 @@
         μ: &RNDM<N, F>,
         τθ_or_adaptive: &mut TransportStepLength<F>,
         tconfig: &TransportConfig<F>,
+        domain: Option<&Cube<N, F>>,
     ) where
         D: DifferentiableRealMapping<N, F>,
     {
@@ -810,7 +818,8 @@
                         if ρ.fail_count < tconfig.max_fail {
                             let dv_x = v.differential(&ρ.x);
                             let g = &dv_x * (ρ.α_γ.signum() * τθ);
-                            ρ.y = ρ.x - g;
+                            let tmp = ρ.x - g;
+                            ρ.y = if domain.is_none_or(|d| d.contains(&tmp)) { tmp } else { ρ.x };
                             let n = g.norm2();
                             if n >= F::EPSILON {
                                 // Estimate Lipschitz factor of ∇v
@@ -835,7 +844,8 @@
         if let Some(τθ) = simple_τθ {
             for ρ in self.iter_mut() {
                 if ρ.fail_count <= tconfig.max_fail {
-                    ρ.y = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * τθ);
+                    let tmp = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * τθ);
+                    ρ.y = if domain.is_none_or(|d| d.contains(&tmp)) { tmp } else { ρ.x };
                 }
             }
         }
@@ -964,10 +974,11 @@
 /// Iteratively solve the pointsource localisation problem using sliding forward-backward
 /// splitting
 ///
-/// The parametrisation is as for [`pointsource_fb_reg`].
+/// The parametrisation is as for [`pointsource_fb_reg`], except an additional domain restriction
+/// may be given.
 /// Inertia is currently not supported.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_sliding_fb_reg<F, I, Dat, Reg, Plot, P, const N: usize>(
+pub fn pointsource_sliding_fb_dom<F, I, Dat, Reg, Plot, P, const N: usize>(
     f: &Dat,
     reg: &Reg,
     prox_penalty: &P,
@@ -975,6 +986,7 @@
     iterator: I,
     mut plotter: Plot,
     μ0: Option<RNDM<N, F>>,
+    domain: Option<&P::Dom>,
 ) -> DynResult<RNDM<N, F>>
 where
     F: Float + ToNalgebraRealField,
@@ -1037,6 +1049,7 @@
             &mut τθ_or_adaptive,
             &v,
             &config.transport,
+            domain
         );
 
         let mut attempts = 0;
@@ -1148,3 +1161,29 @@
     //postprocess(μ, &config.insertion, f)
     postprocess(μ, &config.insertion, |μ̃| f.apply(μ̃))
 }
+
+#[deprecated]
+/// Version of `pointsource_sliding_fb` without domain restriction
+pub fn pointsource_sliding_fb_reg<F, I, Dat, Reg, Plot, P, const N: usize>(
+    f: &Dat,
+    reg: &Reg,
+    prox_penalty: &P,
+    config: &SlidingFBConfig<F>,
+    iterator: I,
+    plotter: Plot,
+    μ0: Option<RNDM<N, F>>,
+) -> DynResult<RNDM<N, F>>
+where
+    F: Float + ToNalgebraRealField,
+    I: AlgIteratorFactory<IterInfo<F>>,
+    Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F> + BoundedCurvature<F>,
+    Dat::DerivativeDomain:
+        DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>,
+    //for<'a> Dat::Differential<'a>: Lipschitz<&'a P, FloatType = F>,
+    RNDM<N, F>: SpikeMerging<F>,
+    Reg: SlidingRegTerm<Loc<N, F>, F>,
+    P: TransportProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>,
+    Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>,
+{
+    pointsource_sliding_fb_dom(f, reg, prox_penalty, config, iterator, plotter, μ0, None)
+}