Mercurial > repos > pointsource_algs / file comparison
comparison: src/kernels/gaussian.rs
src/kernels/gaussian.rs
- changeset 3
- 0778a71cbb6a
- parent 0
- eb3c7813b67a
equal
deleted
inserted
replaced
| 0:eb3c7813b67a | 3:0778a71cbb6a |
|---|---|
| 14 Bounds, | 14 Bounds, |
| 15 LocalAnalysis, | 15 LocalAnalysis, |
| 16 GlobalAnalysis, | 16 GlobalAnalysis, |
| 17 Weighted, | 17 Weighted, |
| 18 Bounded, | 18 Bounded, |
| 19 Taylor2Model, | |
| 20 Taylor2ModelParams, | |
| 19 }; | 21 }; |
| 20 use alg_tools::mapping::Apply; | 22 use alg_tools::mapping::Apply; |
| 21 use alg_tools::maputil::array_init; | 23 use alg_tools::maputil::array_init; |
| 22 | 24 |
| 23 use crate::fourier::Fourier; | 25 use crate::fourier::Fourier; |
| 222 fn apply(&self, y : Loc<F, N>) -> F { | 224 fn apply(&self, y : Loc<F, N>) -> F { |
| 223 self.apply(&y) | 225 self.apply(&y) |
| 224 } | 226 } |
| 225 } | 227 } |
| 226 | 228 |
| 229 #[replace_float_literals(F::cast_from(literal))] | |
| 230 impl<'a, F : Float, R, C, S, const N : usize> Taylor2Model<F, N> | |
| 231 for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> | |
| 232 where R : Constant<Type=F>, | |
| 233 C : Constant<Type=F>, | |
| 234 S : Constant<Type=F>, | |
| 235 Loc<Loc<F, 3>, N> : Product2Taylor<F, N> { | |
| 236 | |
| 237 fn taylor2model(&self, y : Loc<F, N>) -> Taylor2ModelParams<F, N> { | |
| 238 let Convolution(ref ind, | |
| 239 SupportProductFirst(ref cut, | |
| 240 ref gaussian)) = self; | |
| 241 let a = cut.r.value(); | |
| 242 let b = ind.r.value(); | |
| 243 let σ = gaussian.variance.value().sqrt(); | |
| 244 let π = F::PI; | |
| 245 let t = F::SQRT_2 * σ; | |
| 246 let c = σ * (8.0/π).sqrt(); | |
| 247 let sc = gaussian.scale(); | |
| 248 | |
| 249 // This is just a product of one-dimensional versions | |
| 250 y.map(|x| { | |
| 251 let c1 = -(a.min(b + x)); //(-a).max(-x-b); | |
| 252 let c2 = a.min(b - x); | |
| 253 if c1 >= c2 { | |
| 254 Loc([0.0, 0.0, 0.0]) | |
| 255 } else { | |
| 256 let s1 = c1 / t; | |
| 257 let s2 = c2 / t; | |
| 258 let e1 = F::cast_from(erf(s1.as_())); | |
| 259 let e2 = F::cast_from(erf(s2.as_())); | |
| 260 let d1 = F::FRAC_2_SQRT_PI * (-s1*s1).exp(); | |
| 261 let d2 = F::FRAC_2_SQRT_PI * (-s2*s2).exp(); | |
| 262 let dd1 = -2.0 * s1 * d1; | |
| 263 let dd2 = -2.0 * s2 * s2; | |
| 264 debug_assert!(e2 >= e1); | |
| 265 Loc([c * (e2 - e1), c * (d2 - d1), c * (dd2 - dd1)]) | |
| 266 } | |
| 267 }).product2taylor_scaled(sc) | |
| 268 } | |
| 269 | |
| 270 } | |
| 271 | |
| 227 impl<F : Float, R, C, S, const N : usize> | 272 impl<F : Float, R, C, S, const N : usize> |
| 228 Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> | 273 Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> |
| 229 where R : Constant<Type=F>, | 274 where R : Constant<Type=F>, |
| 230 C : Constant<Type=F>, | 275 C : Constant<Type=F>, |
| 231 S : Constant<Type=F> { | 276 S : Constant<Type=F> { |
