Mercurial > repos > pointsource_algs / file diff

--- a/src/kernels/gaussian.rs	Tue Dec 31 09:34:24 2024 -0500
+++ b/src/kernels/gaussian.rs	Tue Aug 01 10:32:12 2023 +0300
@@ -220,12 +220,11 @@
         let a = cut.r.value();
         let b = ind.r.value();
         let σ = gaussian.variance.value().sqrt();
-        let π = F::PI;
         let t = F::SQRT_2 * σ;
-        let c = σ * (8.0/π).sqrt();
+        let c = 0.5; // 1/(σ√(2π) * σ√(π/2) = 1/2
         
         // This is just a product of one-dimensional versions
-        let unscaled = y.product_map(|x| {
+        y.product_map(|x| {
             let c1 = -(a.min(b + x)); //(-a).max(-x-b);
             let c2 = a.min(b - x);
             if c1 >= c2 {
@@ -236,9 +235,7 @@
                 debug_assert!(e2 >= e1);
                 c * (e2 - e1)
             }
-        });
-        
-        unscaled / gaussian.scale()
+        })
     }
 }
 
@@ -276,10 +273,9 @@
         let a = cut.r.value();
         let b = ind.r.value();
         let σ = gaussian.variance.value().sqrt();
-        let π = F::PI;
         let t = F::SQRT_2 * σ;
-        let c = σ * (8.0/π).sqrt();
-        let cd = (8.0).sqrt(); // σ * (8.0/π).sqrt() / t * (√2/π)
+        let c = 0.5; // 1/(σ√(2π) * σ√(π/2) = 1/2
+        let c_div_t = c / t;
         
         // Calculate the values for all component functions of the
         // product. This is just the loop from apply above.
@@ -306,9 +302,9 @@
                 // from the chain rule
                 let de1 = (-(c1/t).powi(2)).exp();
                 let de2 = (-(c2/t).powi(2)).exp();
-                cd * (de1 - de2)
+                c_div_t * (de1 - de2)
             }
-        }) / gaussian.scale()
+        })
     }
 }