Mercurial > repos > pointsource_algs / file comparison
comparison: src/kernels/gaussian.rs
src/kernels/gaussian.rs
- changeset 31
- 6105b5cd8d89
- parent 0
- eb3c7813b67a
- child 33
- aec67cdd6b14
equal
deleted
inserted
replaced
| 26:acf57c458740 | 31:6105b5cd8d89 |
|---|---|
| 186 SupportProductFirst(ref cut, | 186 SupportProductFirst(ref cut, |
| 187 ref gaussian)) = self; | 187 ref gaussian)) = self; |
| 188 let a = cut.r.value(); | 188 let a = cut.r.value(); |
| 189 let b = ind.r.value(); | 189 let b = ind.r.value(); |
| 190 let σ = gaussian.variance.value().sqrt(); | 190 let σ = gaussian.variance.value().sqrt(); |
| 191 let π = F::PI; | |
| 192 let t = F::SQRT_2 * σ; | 191 let t = F::SQRT_2 * σ; |
| 193 let c = σ * (8.0/π).sqrt(); | 192 let c = 0.5; // 1/(σ√(2π) * σ√(π/2) = 1/2 |
| 194 | 193 |
| 195 // This is just a product of one-dimensional versions | 194 // This is just a product of one-dimensional versions |
| 196 let unscaled = y.product_map(|x| { | 195 y.product_map(|x| { |
| 197 let c1 = -(a.min(b + x)); //(-a).max(-x-b); | 196 let c1 = -(a.min(b + x)); //(-a).max(-x-b); |
| 198 let c2 = a.min(b - x); | 197 let c2 = a.min(b - x); |
| 199 if c1 >= c2 { | 198 if c1 >= c2 { |
| 200 0.0 | 199 0.0 |
| 201 } else { | 200 } else { |
| 202 let e1 = F::cast_from(erf((c1 / t).as_())); | 201 let e1 = F::cast_from(erf((c1 / t).as_())); |
| 203 let e2 = F::cast_from(erf((c2 / t).as_())); | 202 let e2 = F::cast_from(erf((c2 / t).as_())); |
| 204 debug_assert!(e2 >= e1); | 203 debug_assert!(e2 >= e1); |
| 205 c * (e2 - e1) | 204 c * (e2 - e1) |
| 206 } | 205 } |
| 207 }); | 206 }) |
| 208 | |
| 209 unscaled / gaussian.scale() | |
| 210 } | 207 } |
| 211 } | 208 } |
| 212 | 209 |
| 213 impl<F : Float, R, C, S, const N : usize> Apply<Loc<F, N>> | 210 impl<F : Float, R, C, S, const N : usize> Apply<Loc<F, N>> |
| 214 for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> | 211 for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> |
