fx/
warp.rs

1// warp -- warp framework of the distort effects: inverse mapping per destination pixel, n x n
2// supersampling, bilinear sampling in premultiplied space, edge behaviour; the f64 maths the
3// fx kernels share                                                          [src/fx/fx_warp.s]
4//
5//   prepare(a) -> aux                bounds for centre / radius (the selection's bounding box)
6//   Warp::new(a)                     the WARP record of a render call (fxw_init)
7//   wr.center(x, y)                  centre from an FP_PAN pair (fxw_center)
8//   wr.run(a, dst, map) -> 0 | ERR_CANCEL      render the ROI through map (fxw_run)
9//   sincos(x) -> (sin, cos)          f64, |x| <= 1e9 (fdlibm kernels)        [fxm_sincos]
10//   log(x)                           f64 ln, x clamped to [DBL_MIN, DBL_MAX]   [fxm_log]
11//
12// Coordinates: pixel (x, y) is the POINT (x, y) (pixel centres at integers); the image covers
13// [-0.5, w - 0.5) x [-0.5, h - 0.5).  For every destination pixel of the ROI run takes n x n
14// sample points p = (x - cx + ox_i, y - cy + oy_i) relative to the centre (paint.net's rotated
15// grid: oy_i = (i + 1) / (n*n + 1) - 0.5, ox_i = frac(n (i + 1) / (n*n + 1)) - 0.5; n = 1:
16// one sample at the centre) and calls map(wr, cur, pt) -> identity?  pt = the point relative to
17// the centre on entry, the SOURCE point on return (identity: the sample is the pixel's own
18// source pixel, exactly).  A source point is sampled bilinearly from the 2 x 2 pixels around it
19// with f32 weights, alpha-weighted (each tap contributes w (c a, a)); taps outside the image
20// follow the edge mode (CLAMP / WRAP / REFLECT per tap; PRIMARY / SECONDARY / TRANSPARENT /
21// ORIGINAL replace a sample outside the image area by that colour).  The pixel = sum of its
22// samples: a = round(S_a / (n n)), c = round(S_c / S_a); a = 0 -> 0.  Pixels farther than
23// rin + 0.75 from the centre (rin > 0) are copied unchanged.
24
25use alloc::boxed::Box;
26
27use rt::fx::{cvtd64, cvttd64, maxsd, minsd};
28use simd::{Simd, Sse2};
29
30use crate::*;
31
32pub const FXW_CLAMP: i32 = 0; // clamp taps to the image
33pub const FXW_WRAP: i32 = 1; // tile the image
34pub const FXW_REFLECT: i32 = 2; // mirror (period 2n, edge pixels repeated)
35pub const FXW_PRIMARY: i32 = 3; // samples outside the image = col1
36pub const FXW_SECONDARY: i32 = 4; // samples outside the image = col2
37pub const FXW_TRANSPARENT: i32 = 5; // samples outside the image = transparent
38pub const FXW_ORIGINAL: i32 = 6; // samples outside the image = the destination's own pixel
39pub const FXW_NEDGE: i32 = 7;
40
41const LO: f64 = -1073741824.0;
42const HI: f64 = 1073741824.0;
43const SIGN: u64 = 0x8000_0000_0000_0000;
44
45#[inline(always)]
46fn neg(v: f64) -> f64 {
47    f64::from_bits(v.to_bits() ^ SIGN)
48}
49
50/// the aux of prepare: the bounds x0, y0, x1, y1
51pub struct Bounds(pub [i32; 4]);
52
53/// the bounding box of the selection (clipped to the image) or the whole image  [fxw_prepare]
54pub fn prepare(a: &FxArgs) -> Result<Option<Aux>, i32> {
55    let (w, h) = (a.src.w, a.src.h);
56    let mut b = [0, 0, w, h];
57    if let Some(m) = a.sel {
58        if let Some(bb) = m.bbox() {
59            let c = [bb[0].max(0), bb[1].max(0), bb[2].min(w), bb[3].min(h)];
60            if c[0] < c[2] && c[1] < c[3] {
61                b = c;
62            }
63        }
64    }
65    Ok(Some(Box::new(Bounds(b))))
66}
67
68/// floor of an f64 as i64 (|x| <= 2^30 after clamping)
69#[inline(always)]
70fn floor64(x: f64) -> i64 {
71    let t = cvttd64(x);
72    if t as f64 > x { t - 1 } else { t }
73}
74
75type V = <Sse2 as Simd>::V;
76
77/// straight pixel -> 4 x f32 (b a, g a, r a, a)  [PX2PM]
78#[inline(always)]
79fn px2pm(p: u32) -> V {
80    let s = Sse2;
81    let z = s.zero();
82    lanes2pm(s.unpacklo16(s.unpacklo8(s.from_i32(p as i32), z), z))
83}
84
85/// a pixel as 4 x i32 lanes (b, g, r, a) -> 4 x f32 (b a, g a, r a, a)  [PX2PM]
86#[inline(always)]
87fn lanes2pm(p: V) -> V {
88    let s = Sse2;
89    let v = s.cvt_i32_f32(p);
90    let a = s.shuf32::<0xFF>(v);
91    let m012 = s.load_i32(&[-1, -1, -1, 0]);
92    let one3 = s.load_i32(&[0, 0, 0, 0x3F80_0000]);
93    s.mulf(s.or(s.and(v, m012), one3), a)
94}
95
96/// the bilinear mix of the 4 premultiplied taps (00, 10, 01, 11) with weights wt  [.ws_mix]
97#[inline(always)]
98fn mix(t0: V, t1: V, t2: V, t3: V, wt: &[f32; 4]) -> V {
99    let s = Sse2;
100    let acc = s.addf(s.zero(), s.mulf(t0, s.splatf(wt[0])));
101    let acc = s.addf(acc, s.mulf(t1, s.splatf(wt[1])));
102    let acc = s.addf(acc, s.mulf(t2, s.splatf(wt[2])));
103    s.addf(acc, s.mulf(t3, s.splatf(wt[3])))
104}
105
106/// the WARP record of one render call
107pub struct Warp<'a> {
108    pub src: &'a Surface,
109    pub w: i32,
110    pub h: i32,
111    pub x0: i32,
112    pub y0: i32,
113    pub x1: i32,
114    pub y1: i32,
115    pub job: Option<&'a (dyn Cancel + Sync)>,
116    /// quality 1..5: n x n samples per pixel
117    pub q: i32,
118    pub edge: i32,
119    pub col1: u32,
120    pub col2: u32,
121    pub cx: f64,
122    pub cy: f64,
123    /// radius of influence (0 = map every pixel)
124    pub rin: f64,
125    /// default radius = min(bounds w, h) / 2
126    pub rdef: f64,
127    pub bx0: i32,
128    pub by0: i32,
129    pub bx1: i32,
130    pub by1: i32,
131    /// effect parameters (WR_P)
132    pub p: [f64; 16],
133    pub seed: u32,
134    cv1: [u32; 2],
135}
136
137/// the pixel and sample a map function is called for (WR_CURX / CURY / CURS)
138#[derive(Clone, Copy)]
139pub struct Cur {
140    pub x: i32,
141    pub y: i32,
142    pub s: i32,
143}
144
145#[inline(always)]
146fn clip0(v: i32, lim: i32) -> i32 {
147    if v < 0 {
148        0
149    } else if v > lim {
150        lim
151    } else {
152        v
153    }
154}
155
156/// edge index of a tap coordinate (CLAMP / WRAP / REFLECT)  [EDGE_IDX]
157#[inline(always)]
158fn edge_idx(v: i64, n: i64, mode: i32) -> i64 {
159    match mode {
160        FXW_WRAP => {
161            let r = rem(v, n);
162            if r < 0 { r + n } else { r }
163        }
164        FXW_REFLECT => {
165            let n2 = n + n;
166            let mut r = rem(v, n2);
167            if r < 0 {
168                r += n2;
169            }
170            if r >= n { n2 - 1 - r } else { r }
171        }
172        _ => {
173            if v < 0 {
174                0
175            } else if v >= n {
176                n - 1
177            } else {
178                v
179            }
180        }
181    }
182}
183
184/// the edge indices of the taps v and v + 1 (EDGE_IDX twice; perf-fx: one remainder, none
185/// for v within one period below or above)
186#[inline(always)]
187fn edge_pair(v: i64, n: i64, mode: i32) -> (i64, i64) {
188    match mode {
189        FXW_WRAP => {
190            // r = v mod n in [0, n) (= EDGE_IDX), (v + 1) mod n follows from it
191            let r = if (v as u64) < n as u64 {
192                v
193            } else if v < 0 && v >= -n {
194                v + n
195            } else {
196                let r = rem(v, n);
197                if r < 0 { r + n } else { r }
198            };
199            (r, if r + 1 == n { 0 } else { r + 1 })
200        }
201        FXW_REFLECT => {
202            let n2 = n + n;
203            let m = if (v as u64) < n2 as u64 {
204                v
205            } else if v < 0 && v >= -n2 {
206                v + n2
207            } else {
208                let m = rem(v, n2);
209                if m < 0 { m + n2 } else { m }
210            };
211            let m1 = if m + 1 == n2 { 0 } else { m + 1 };
212            (if m >= n { n2 - 1 - m } else { m }, if m1 >= n { n2 - 1 - m1 } else { m1 })
213        }
214        _ => (edge_idx(v, n, mode), edge_idx(v + 1, n, mode)),
215    }
216}
217
218/// v % n (truncating, n >= 1): perf-fx: a 32-bit idiv when both fit (the taps are clamped to
219/// +-(2^30 + 1), so always for images below 2^30 pixels; the same remainder as the asm's
220/// 64-bit idiv)
221#[inline(always)]
222fn rem(v: i64, n: i64) -> i64 {
223    if v as i32 as i64 == v && n as i32 as i64 == n { ((v as i32) % (n as i32)) as i64 } else { v % n }
224}
225
226impl<'a> Warp<'a> {
227    /// the record from the arguments: ROI clipped to the image, bounds (aux of prepare or the
228    /// image), rdef, centre of the bounds, q = 1, edge CLAMP, rin = 0  [fxw_init]
229    pub fn new(a: &FxArgs<'a>) -> Warp<'a> {
230        let (w, h) = (a.src.w, a.src.h);
231        let mut b = [0, 0, w, h];
232        if let Some(Bounds(c)) = a.aux::<Bounds>() {
233            if c[0] >= 0 && c[1] >= 0 && c[0] < c[2] && c[1] < c[3] && c[2] <= w && c[3] <= h {
234                b = *c;
235            }
236        }
237        let rd = (b[2] - b[0]).min(b[3] - b[1]);
238        Warp {
239            src: a.src,
240            w,
241            h,
242            x0: clip0(a.x0, w),
243            y0: clip0(a.y0, h),
244            x1: clip0(a.x1, w),
245            y1: clip0(a.y1, h),
246            job: a.job,
247            q: 1,
248            edge: FXW_CLAMP,
249            col1: 0,
250            col2: 0,
251            cx: (b[0] + b[2] - 1) as f64 * 0.5,
252            cy: (b[1] + b[3] - 1) as f64 * 0.5,
253            rin: 0.0,
254            rdef: rd as f64 * 0.5,
255            bx0: b[0],
256            by0: b[1],
257            bx1: b[2],
258            by1: b[3],
259            p: [0.0; 16],
260            seed: 0,
261            cv1: [0; 2],
262        }
263    }
264
265    /// centre from an FP_PAN value (-1 = left / top edge of the bounds, 1 = right / bottom;
266    /// clamped to -2 .. 2): cx = bx0 - 0.5 + bw (1 + x) / 2  [fxw_center]
267    pub fn center(&mut self, x: f64, y: f64) {
268        let x = minsd(maxsd(x, -2.0), 2.0);
269        let y = minsd(maxsd(y, -2.0), 2.0);
270        self.cx = (x + 1.0) * (self.bx1 - self.bx0) as f64 * 0.5 + (self.bx0 as f64 - 0.5);
271        self.cy = (y + 1.0) * (self.by1 - self.by0) as f64 * 0.5 + (self.by0 as f64 - 0.5);
272    }
273
274    #[inline(always)]
275    fn spx(&self, x: i64, y: i64) -> u32 {
276        px(&self.src.px[y as usize * self.src.stride as usize..], x as usize)
277    }
278
279    /// the 4 taps (x0 y0, x1 y0, x0 y1, x1 y1) of a sample whose 2 x 2 pixels are not all
280    /// inside, through the edge mode  [wr_sample .ws_slow]
281    #[inline(always)]
282    fn edge_taps(&self, ix: i64, iy: i64, mode: i32) -> [u32; 4] {
283        let (w, h) = (self.w as i64, self.h as i64);
284        let (x0, x1) = edge_pair(ix, w, mode);
285        let (y0, y1) = edge_pair(iy, h, mode);
286        [self.spx(x0, y0), self.spx(x1, y0), self.spx(x0, y1), self.spx(x1, y1)]
287    }
288
289    /// premultiplied bilinear sample at the absolute point (sx, sy) (weights sum 1)  [wr_sample]
290    /// (perf-fx: inline in the sample loop; a call cost the Win64 build 8 pushes + 3 xmm saves
291    /// and the vector result through memory per sample)
292    #[inline(always)]
293    fn sample(&self, cur: Cur, sx: f64, sy: f64) -> V {
294        let sx = minsd(maxsd(sx, LO), HI);
295        let sy = minsd(maxsd(sy, LO), HI);
296        let mut mode = self.edge;
297        if mode >= FXW_PRIMARY {
298            let out = sx < -0.5 || sy < -0.5 || !(sx < self.w as f64 - 0.5) || !(sy < self.h as f64 - 0.5);
299            if out {
300                return match mode {
301                    FXW_PRIMARY => px2pm(self.cv1[0]),
302                    FXW_SECONDARY => px2pm(self.cv1[1]),
303                    FXW_ORIGINAL => px2pm(self.spx(cur.x as i64, cur.y as i64)),
304                    _ => Sse2.zero(),
305                };
306            }
307            mode = FXW_CLAMP;
308        }
309        let ix = floor64(sx);
310        let fx = (sx - ix as f64) as f32;
311        let iy = floor64(sy);
312        let fy = (sy - iy as f64) as f32;
313        let wt = [(1.0 - fx) * (1.0 - fy), fx * (1.0 - fy), (1.0 - fx) * fy, fx * fy];
314        let (w, h) = (self.w as i64, self.h as i64);
315        if (ix as u64) < (w - 1) as u64 && (iy as u64) < (h - 1) as u64 {
316            // perf-fx: the 2 pixels of each row with one 8-byte load (movq, as the asm), kept
317            // in registers
318            let s = Sse2;
319            let z = s.zero();
320            let st = self.src.stride as usize;
321            let o = iy as usize * st + ix as usize * 4;
322            let r0 = s.widen_px(s.load_raw(&self.src.px[o..]));
323            let r1 = s.widen_px(s.load_raw(&self.src.px[o + st..]));
324            mix(
325                lanes2pm(s.unpacklo16(r0, z)),
326                lanes2pm(s.unpackhi16(r0, z)),
327                lanes2pm(s.unpacklo16(r1, z)),
328                lanes2pm(s.unpackhi16(r1, z)),
329                &wt,
330            )
331        } else {
332            let t = self.edge_taps(ix, iy, mode);
333            mix(px2pm(t[0]), px2pm(t[1]), px2pm(t[2]), px2pm(t[3]), &wt)
334        }
335    }
336
337    /// render the ROI through map -> 0 | ERR_CANCEL  [fxw_run]
338    pub fn run(&mut self, dst: &mut Dst, mut map: impl FnMut(&Warp, Cur, &mut [f64; 2]) -> bool) -> i32 {
339        let n = self.q.clamp(1, 5);
340        self.q = n;
341        let ns = (n * n) as usize;
342        let mut ofs = [[0f64; 2]; 25];
343        if ns != 1 {
344            let d = ns as i32 as f64 + 1.0;
345            for (i, o) in ofs.iter_mut().enumerate().take(ns) {
346                let y = (i as i32 + 1) as f64 / d;
347                let mut x = y * n as f64;
348                x -= rt::fx::cvttd(x) as f64;
349                *o = [x - 0.5, y - 0.5];
350            }
351        }
352        self.cv1 = [self.col1, self.col2];
353        if self.edge as u32 >= FXW_NEDGE as u32 {
354            self.edge = FXW_CLAMP;
355        }
356        let rlim = if self.rin > 0.0 {
357            let r = self.rin + 0.75;
358            r * r
359        } else {
360            -1.0
361        };
362        let me: &Warp = self;
363        for y in me.y0..me.y1 {
364            if me.job.is_some_and(|j| j.cancelled()) {
365                return ERR_CANCEL;
366            }
367            let dy = y as f64 - me.cy;
368            let dx0 = dst.x0;
369            let drow = dst.row(y);
370            for x in me.x0..me.x1 {
371                let dx = x as f64 - me.cx;
372                let o = (x - dx0) as usize;
373                if !(rlim < 0.0) && !(dx * dx + dy * dy <= rlim) {
374                    put(drow, o, me.spx(x as i64, y as i64));
375                    continue;
376                }
377                // perf-fx: first the source points of every sample (map is a pure function of
378                // its arguments), then the samples: the taps of a pixel's samples are fetched
379                // close together, so their cache misses overlap; same sums in the same order
380                let mut sp = [[0f64; 2]; 25];
381                let mut own = [false; 25];
382                for (s, of) in ofs.iter().enumerate().take(ns) {
383                    let mut pt = [of[0] + dx, of[1] + dy];
384                    own[s] = map(me, Cur { x, y, s: s as i32 }, &mut pt);
385                    sp[s] = [pt[0] + me.cx, pt[1] + me.cy];
386                }
387                let mut acc = Sse2.zero();
388                for s in 0..ns {
389                    let v = if own[s] {
390                        px2pm(me.spx(x as i64, y as i64))
391                    } else {
392                        me.sample(Cur { x, y, s: s as i32 }, sp[s][0], sp[s][1])
393                    };
394                    acc = Sse2.addf(acc, v);
395                }
396                put(drow, o, finish(acc, ns as i32));
397            }
398        }
399        0
400    }
401}
402
403/// sum (S_b, S_g, S_r, S_a) of n samples -> straight pixel: rint(S_c / S_a), rint(S_a / n),
404/// saturated to i16, then to u8  [wr_finish]
405#[inline]
406fn finish(acc: V, n: i32) -> u32 {
407    let s = Sse2;
408    let mut t = [0f32; 4];
409    s.store_f32(&mut t, acc);
410    let sa = t[3];
411    if !(sa > 0.0) {
412        return 0;
413    }
414    let d = s.load_f32(&[sa, sa, sa, n as f32]);
415    let v = s.cvtr_f32_i32(s.divf(acc, d));
416    let q = s.packs32(v, v);
417    let p = s.to_i32(s.packus16(q, q)) as u32;
418    if p < 0x0100_0000 { 0 } else { p }
419}
420
421// ---------------------------------------------------------------- f64 maths
422
423const PIO2_1: f64 = 1.57079632673412561417e+00;
424const PIO2_1T: f64 = 6.07710050650619224932e-11;
425const S1: f64 = -1.66666666666666324348e-01;
426const S2: f64 = 8.33333333332248946124e-03;
427const S3: f64 = -1.98412698298579493134e-04;
428const S4: f64 = 2.75573137070700676789e-06;
429const S5: f64 = -2.50507602534068634195e-08;
430const S6: f64 = 1.58969099521155010221e-10;
431const C1: f64 = 4.16666666666666019037e-02;
432const C2: f64 = -1.38888888888741095749e-03;
433const C3: f64 = 2.48015872894767294178e-05;
434const C4: f64 = -2.75573143513906633035e-07;
435const C5: f64 = 2.08757232129817482790e-09;
436const C6: f64 = -1.13596475577881948265e-11;
437
438/// (sin x, cos x): k = rint(x 2/pi), Cody-Waite reduction, fdlibm kernels; x clamped to
439/// +-1e9 (NaN -> 0)  [fxm_sincos]
440pub fn sincos(x: f64) -> (f64, f64) {
441    let x = if x.is_nan() { 0.0 } else { minsd(maxsd(x, -1.0e9), 1.0e9) };
442    let k = cvtd64(x * 0.63661977236758134308);
443    let kd = k as f64;
444    let mut r = x - kd * PIO2_1;
445    r -= kd * PIO2_1T;
446    let z = r * r;
447    let ps = ((((S6 * z + S5) * z + S4) * z + S3) * z + S2) * z + S1;
448    let s = ps * z * r + r;
449    let pc = ((((C6 * z + C5) * z + C4) * z + C3) * z + C2) * z + C1;
450    let c = (1.0 - z * 0.5) + pc * z * z;
451    match k & 3 {
452        1 => (c, neg(s)),
453        2 => (neg(s), neg(c)),
454        3 => (neg(c), s),
455        _ => (s, c),
456    }
457}
458
459const LN2HI: f64 = 6.93147180369123816490e-01;
460const LN2LO: f64 = 1.90821492927058770002e-10;
461const LG1: f64 = 6.666666666666735130e-01;
462const LG2: f64 = 3.999999999940941908e-01;
463const LG3: f64 = 2.857142874366239149e-01;
464const LG4: f64 = 2.222219843214978396e-01;
465const LG5: f64 = 1.818357216161805012e-01;
466const LG6: f64 = 1.531383769920937332e-01;
467const LG7: f64 = 1.479819860511658591e-01;
468
469/// ln x (fdlibm e_log.c); x clamped to [DBL_MIN, DBL_MAX] (NaN, 0, negatives -> ln DBL_MIN)
470/// [fxm_log]
471pub fn log(x: f64) -> f64 {
472    let x = minsd(maxsd(x, 2.2250738585072014e-308), 1.7976931348623157e+308);
473    let b = x.to_bits();
474    let mut k = (b >> 52) as i32 - 1023;
475    let mut m = b & 0x000F_FFFF_FFFF_FFFF;
476    if m > 0x0006_A09E_667F_3BCD {
477        m |= 0x3FE0_0000_0000_0000;
478        k += 1;
479    } else {
480        m |= 0x3FF0_0000_0000_0000;
481    }
482    let f = f64::from_bits(m) - 1.0;
483    let s = f / (f + 2.0);
484    let z = s * s;
485    let w = z * z;
486    let t1 = ((LG6 * w + LG4) * w + LG2) * w;
487    let t2 = (((LG7 * w + LG5) * w + LG3) * w + LG1) * z;
488    let r = t2 + t1;
489    let hfsq = f * f * 0.5;
490    let kd = k as f64;
491    let a = (hfsq + r) * s + kd * LN2LO;
492    let b2 = hfsq - a - f;
493    kd * LN2HI - b2
494}
fx/warp.rs

fx/
warp.rs
