looprestoration_tmpl.c 21.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
/*
 * Copyright © 2018, VideoLAN and dav1d authors
 * Copyright © 2018, Two Orioles, LLC
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include "config.h"

#include <stdlib.h>

#include "common/intops.h"

#include "src/looprestoration.h"
#include "src/tables.h"

37 38
// 256 * 1.5 + 3 + 3 = 390
#define REST_UNIT_STRIDE (390)
39 40 41

// TODO Reuse p when no padding is needed (add and remove lpf pixels in p)
// TODO Chroma only requires 2 rows of padding.
42
static void padding(pixel *dst, const pixel *p, const ptrdiff_t p_stride,
43
                    const pixel (*left)[4],
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
                    const pixel *lpf, const ptrdiff_t lpf_stride,
                    int unit_w, const int stripe_h, const enum LrEdgeFlags edges)
{
    const int have_left = !!(edges & LR_HAVE_LEFT);
    const int have_right = !!(edges & LR_HAVE_RIGHT);

    // Copy more pixels if we don't have to pad them
    unit_w += 3 * have_left + 3 * have_right;
    pixel *dst_l = dst + 3 * !have_left;
    p -= 3 * have_left;
    lpf -= 3 * have_left;

    if (edges & LR_HAVE_TOP) {
        // Copy previous loop filtered rows
        const pixel *const above_1 = lpf;
        const pixel *const above_2 = above_1 + PXSTRIDE(lpf_stride);
        pixel_copy(dst_l, above_1, unit_w);
61 62
        pixel_copy(dst_l + REST_UNIT_STRIDE, above_1, unit_w);
        pixel_copy(dst_l + 2 * REST_UNIT_STRIDE, above_2, unit_w);
63 64 65
    } else {
        // Pad with first row
        pixel_copy(dst_l, p, unit_w);
66 67
        pixel_copy(dst_l + REST_UNIT_STRIDE, p, unit_w);
        pixel_copy(dst_l + 2 * REST_UNIT_STRIDE, p, unit_w);
68 69 70 71 72
        if (have_left) {
            pixel_copy(dst_l, &left[0][1], 3);
            pixel_copy(dst_l + REST_UNIT_STRIDE, &left[0][1], 3);
            pixel_copy(dst_l + 2 * REST_UNIT_STRIDE, &left[0][1], 3);
        }
73 74
    }

75
    pixel *dst_tl = dst_l + 3 * REST_UNIT_STRIDE;
76 77 78 79
    if (edges & LR_HAVE_BOTTOM) {
        // Copy next loop filtered rows
        const pixel *const below_1 = lpf + 6 * PXSTRIDE(lpf_stride);
        const pixel *const below_2 = below_1 + PXSTRIDE(lpf_stride);
80 81 82
        pixel_copy(dst_tl + stripe_h * REST_UNIT_STRIDE, below_1, unit_w);
        pixel_copy(dst_tl + (stripe_h + 1) * REST_UNIT_STRIDE, below_2, unit_w);
        pixel_copy(dst_tl + (stripe_h + 2) * REST_UNIT_STRIDE, below_2, unit_w);
83 84 85
    } else {
        // Pad with last row
        const pixel *const src = p + (stripe_h - 1) * PXSTRIDE(p_stride);
86 87 88
        pixel_copy(dst_tl + stripe_h * REST_UNIT_STRIDE, src, unit_w);
        pixel_copy(dst_tl + (stripe_h + 1) * REST_UNIT_STRIDE, src, unit_w);
        pixel_copy(dst_tl + (stripe_h + 2) * REST_UNIT_STRIDE, src, unit_w);
89 90 91 92 93
        if (have_left) {
            pixel_copy(dst_tl + stripe_h * REST_UNIT_STRIDE, &left[stripe_h - 1][1], 3);
            pixel_copy(dst_tl + (stripe_h + 1) * REST_UNIT_STRIDE, &left[stripe_h - 1][1], 3);
            pixel_copy(dst_tl + (stripe_h + 2) * REST_UNIT_STRIDE, &left[stripe_h - 1][1], 3);
        }
94 95 96 97
    }

    // Inner UNIT_WxSTRIPE_H
    for (int j = 0; j < stripe_h; j++) {
98
        pixel_copy(dst_tl + 3 * have_left, p + 3 * have_left, unit_w - 3 * have_left);
99
        dst_tl += REST_UNIT_STRIDE;
100 101 102 103 104 105 106 107 108
        p += PXSTRIDE(p_stride);
    }

    if (!have_right) {
        pixel *pad = dst_l + unit_w;
        pixel *row_last = &dst_l[unit_w - 1];
        // Pad 3x(STRIPE_H+6) with last column
        for (int j = 0; j < stripe_h + 6; j++) {
            pixel_set(pad, *row_last, 3);
109 110
            pad += REST_UNIT_STRIDE;
            row_last += REST_UNIT_STRIDE;
111 112 113 114 115 116 117
        }
    }

    if (!have_left) {
        // Pad 3x(STRIPE_H+6) with first column
        for (int j = 0; j < stripe_h + 6; j++) {
            pixel_set(dst, *dst_l, 3);
118 119
            dst += REST_UNIT_STRIDE;
            dst_l += REST_UNIT_STRIDE;
120
        }
121 122 123 124 125 126
    } else {
        dst += 3 * REST_UNIT_STRIDE;
        for (int j = 0; j < stripe_h; j++) {
            pixel_copy(dst, &left[j][1], 3);
            dst += REST_UNIT_STRIDE;
        }
127 128 129 130 131 132 133 134
    }
}

// FIXME Could split into luma and chroma specific functions,
// (since first and last tops are always 0 for chroma)
// FIXME Could implement a version that requires less temporary memory
// (should be possible to implement with only 6 rows of temp storage)
static void wiener_c(pixel *p, const ptrdiff_t p_stride,
135
                     const pixel (*const left)[4],
136 137 138
                     const pixel *lpf, const ptrdiff_t lpf_stride,
                     const int w, const int h,
                     const int16_t filterh[7], const int16_t filterv[7],
Ronald S. Bultje's avatar
Ronald S. Bultje committed
139
                     const enum LrEdgeFlags edges HIGHBD_DECL_SUFFIX)
140
{
141 142 143
    // Wiener filtering is applied to a maximum stripe height of 64 + 3 pixels
    // of padding above and below
    pixel tmp[70 /*(64 + 3 + 3)*/ * REST_UNIT_STRIDE];
144 145
    pixel *tmp_ptr = tmp;

146
    padding(tmp, p, p_stride, left, lpf, lpf_stride, w, h, edges);
147 148 149

    // Values stored between horizontal and vertical filtering don't
    // fit in a uint8_t.
150
    uint16_t hor[70 /*(64 + 3 + 3)*/ * REST_UNIT_STRIDE];
151 152
    uint16_t *hor_ptr = hor;

Ronald S. Bultje's avatar
Ronald S. Bultje committed
153 154
    const int bitdepth = bitdepth_from_max(bitdepth_max);
    const int round_bits_h = 3 + (bitdepth == 12) * 2;
155
    const int rounding_off_h = 1 << (round_bits_h - 1);
Ronald S. Bultje's avatar
Ronald S. Bultje committed
156
    const int clip_limit = 1 << (bitdepth + 1 + 7 - round_bits_h);
157 158
    for (int j = 0; j < h + 6; j++) {
        for (int i = 0; i < w; i++) {
Ronald S. Bultje's avatar
Ronald S. Bultje committed
159
            int sum = (tmp_ptr[i + 3] << 7) + (1 << (bitdepth + 6));
160 161 162 163 164 165

            for (int k = 0; k < 7; k++) {
                sum += tmp_ptr[i + k] * filterh[k];
            }

            hor_ptr[i] =
166
                iclip((sum + rounding_off_h) >> round_bits_h, 0, clip_limit - 1);
167
        }
168 169
        tmp_ptr += REST_UNIT_STRIDE;
        hor_ptr += REST_UNIT_STRIDE;
170 171
    }

Ronald S. Bultje's avatar
Ronald S. Bultje committed
172
    const int round_bits_v = 11 - (bitdepth == 12) * 2;
173
    const int rounding_off_v = 1 << (round_bits_v - 1);
Ronald S. Bultje's avatar
Ronald S. Bultje committed
174
    const int round_offset = 1 << (bitdepth + (round_bits_v - 1));
175 176
    for (int i = 0; i < w; i++) {
        for (int j = 0; j < h; j++) {
177
            int sum = (hor[(j + 3) * REST_UNIT_STRIDE + i] << 7) - round_offset;
178 179

            for (int k = 0; k < 7; k++) {
180
                sum += hor[(j + k) * REST_UNIT_STRIDE + i] * filterv[k];
181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
            }

            p[j * PXSTRIDE(p_stride) + i] =
                iclip_pixel((sum + rounding_off_v) >> round_bits_v);
        }
    }
}

// Sum over a 3x3 area
// The dst and src pointers are positioned 3 pixels above and 3 pixels to the
// left of the top left corner. However, the self guided filter only needs 1
// pixel above and one pixel to the left. As for the pixels below and to the
// right they must be computed in the sums, but don't need to be stored.
//
// Example for a 4x4 block:
//      x x x x x x x x x x
//      x c c c c c c c c x
//      x i s s s s s s i x
//      x i s s s s s s i x
//      x i s s s s s s i x
//      x i s s s s s s i x
//      x i s s s s s s i x
//      x i s s s s s s i x
//      x c c c c c c c c x
//      x x x x x x x x x x
//
// s: Pixel summed and stored
// i: Pixel summed and stored (between loops)
// c: Pixel summed not stored
// x: Pixel not summed not stored
211
static void boxsum3(coef *dst, const pixel *src, const int w, const int h) {
212
    // We skip the first row, as it is never used
213 214
    src += REST_UNIT_STRIDE;
    dst += REST_UNIT_STRIDE;
215 216 217 218 219

    // We skip the first and last columns, as they are never used
    for (int x = 1; x < w - 1; x++) {
        coef *ds = dst + x;
        const pixel *s = src + x;
220
        int a = s[0], b = s[REST_UNIT_STRIDE];
221 222 223 224

        // We skip the first 2 rows, as they are skipped in the next loop and
        // we don't need the last 2 row as it is skipped in the next loop
        for (int y = 2; y < h - 2; y++) {
225 226 227
            s += REST_UNIT_STRIDE;
            const int c = s[REST_UNIT_STRIDE];
            ds += REST_UNIT_STRIDE;
228 229 230 231 232 233 234
            *ds = a + b + c;
            a = b;
            b = c;
        }
     }

    // We skip the first 2 rows as they are never read
235
    dst += REST_UNIT_STRIDE;
236 237 238 239 240 241 242 243 244 245 246 247
    // We skip the last 2 rows as it is never read
    for (int y = 2; y < h - 2; y++) {
        int a = dst[1], b = dst[2];

        // We don't store the first column as it is never read and
        // we don't store the last 2 columns as they are never read
        for (int x = 2; x < w - 2; x++) {
            const int c = dst[x + 1];
            dst[x] = a + b + c;
            a = b;
            b = c;
        }
248
        dst += REST_UNIT_STRIDE;
249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
    }
}

// Sum over a 5x5 area
// The dst and src pointers are positioned 3 pixels above and 3 pixels to the
// left of the top left corner. However, the self guided filter only needs 1
// pixel above and one pixel to the left. As for the pixels below and to the
// right they must be computed in the sums, but don't need to be stored.
//
// Example for a 4x4 block:
//      c c c c c c c c c c
//      c c c c c c c c c c
//      i i s s s s s s i i
//      i i s s s s s s i i
//      i i s s s s s s i i
//      i i s s s s s s i i
//      i i s s s s s s i i
//      i i s s s s s s i i
//      c c c c c c c c c c
//      c c c c c c c c c c
//
// s: Pixel summed and stored
// i: Pixel summed and stored (between loops)
// c: Pixel summed not stored
// x: Pixel not summed not stored
274
static void boxsum5(coef *dst, const pixel *const src, const int w, const int h) {
275
    // We skip the first row, as it is never used
276
    dst += REST_UNIT_STRIDE;
277 278 279

    for (int x = 0; x < w; x++) {
        coef *ds = dst + x;
280 281 282 283
        const pixel *s = src + 3 * REST_UNIT_STRIDE + x;
        int a = s[-3 * REST_UNIT_STRIDE];
        int b = s[-2 * REST_UNIT_STRIDE];
        int c = s[-1 * REST_UNIT_STRIDE];
284 285 286 287 288
        int d = s[0];

        // We skip the first 2 rows, as they are skipped in the next loop and
        // we don't need the last 2 row as it is skipped in the next loop
        for (int y = 2; y < h - 2; y++) {
289
            s += REST_UNIT_STRIDE;
290
            const int e = *s;
291
            ds += REST_UNIT_STRIDE;
292 293 294 295 296 297 298 299 300
            *ds = a + b + c + d + e;
            a = b;
            b = c;
            c = d;
            d = e;
        }
    }

    // We skip the first 2 rows as they are never read
301
    dst += REST_UNIT_STRIDE;
302 303 304 305 306 307 308 309 310 311 312 313 314 315
    for (int y = 2; y < h - 2; y++) {
        int a = dst[0];
        int b = dst[1];
        int c = dst[2];
        int d = dst[3];

        for (int x = 2; x < w - 2; x++) {
            const int e = dst[x + 2];
            dst[x] = a + b + c + d + e;
            a = b;
            b = c;
            c = d;
            d = e;
        }
316
        dst += REST_UNIT_STRIDE;
317 318 319 320
    }
}

// See boxsum3 function comments for details on row and column skipping
321
static void boxsum3sqr(int32_t *dst, const pixel *src, const int w, const int h) {
322
    // We skip the first row, as it is never used
323 324
    src += REST_UNIT_STRIDE;
    dst += REST_UNIT_STRIDE;
325 326 327

    // We skip the first and last columns, as they are never used
    for (int x = 1; x < w - 1; x++) {
328
        int32_t *ds = dst + x;
329 330
        const pixel *s = src + x;
        int a = s[0] * s[0];
331
        int b = s[REST_UNIT_STRIDE] * s[REST_UNIT_STRIDE];
332 333 334 335

        // We skip the first row, as it is skipped in the next loop and
        // we don't need the last row as it is skipped in the next loop
        for (int y = 2; y < h - 2; y++) {
336 337 338
            s += REST_UNIT_STRIDE;
            const int c = s[REST_UNIT_STRIDE] * s[REST_UNIT_STRIDE];
            ds += REST_UNIT_STRIDE;
339 340 341 342 343 344 345
            *ds = a + b + c;
            a = b;
            b = c;
        }
     }

    // We skip the first row as it is never read
346
    dst += REST_UNIT_STRIDE;
347 348 349 350 351 352 353 354 355 356 357 358
    // We skip the last row as it is never read
    for (int y = 2; y < h - 2; y++) {
        int a = dst[1], b = dst[2];

        // We don't store the first column as it is never read and
        // we don't store the last 2 columns as they are never read
        for (int x = 2; x < w - 2; x++) {
            const int c = dst[x + 1];
            dst[x] = a + b + c;
            a = b;
            b = c;
        }
359
        dst += REST_UNIT_STRIDE;
360 361 362 363
    }
}

// See boxsum5 function comments for details on row and column skipping
364 365
static void boxsum5sqr(int32_t *dst, const pixel *const src, const int w,
                       const int h)
366 367
{
    // We skip the first row, as it is never used
368
    dst += REST_UNIT_STRIDE;
369 370

    for (int x = 0; x < w; x++) {
371
        int32_t *ds = dst + x;
372 373 374 375
        const pixel *s = src + 3 * REST_UNIT_STRIDE + x;
        int a = s[-3 * REST_UNIT_STRIDE] * s[-3 * REST_UNIT_STRIDE];
        int b = s[-2 * REST_UNIT_STRIDE] * s[-2 * REST_UNIT_STRIDE];
        int c = s[-1 * REST_UNIT_STRIDE] * s[-1 * REST_UNIT_STRIDE];
376 377 378 379 380
        int d = s[0] * s[0];

        // We skip the first 2 rows, as they are skipped in the next loop and
        // we don't need the last 2 row as it is skipped in the next loop
        for (int y = 2; y < h - 2; y++) {
381
            s += REST_UNIT_STRIDE;
382
            const int e = s[0] * s[0];
383
            ds += REST_UNIT_STRIDE;
384 385 386 387 388 389 390 391 392
            *ds = a + b + c + d + e;
            a = b;
            b = c;
            c = d;
            d = e;
        }
    }

    // We skip the first 2 rows as they are never read
393
    dst += REST_UNIT_STRIDE;
394 395 396 397 398 399 400 401 402 403 404 405 406 407
    for (int y = 2; y < h - 2; y++) {
        int a = dst[0];
        int b = dst[1];
        int c = dst[2];
        int d = dst[3];

        for (int x = 2; x < w - 2; x++) {
            const int e = dst[x + 2];
            dst[x] = a + b + c + d + e;
            a = b;
            b = c;
            c = d;
            d = e;
        }
408
        dst += REST_UNIT_STRIDE;
409 410 411
    }
}

Ronald S. Bultje's avatar
Ronald S. Bultje committed
412
static void selfguided_filter(coef *dst, const pixel *src,
413
                              const ptrdiff_t src_stride, const int w,
Ronald S. Bultje's avatar
Ronald S. Bultje committed
414 415
                              const int h, const int n, const int s
                              HIGHBD_DECL_SUFFIX)
416
{
Luc Trudeau's avatar
Luc Trudeau committed
417 418
    const int sgr_one_by_x = n == 25 ? 164 : 455;

419 420 421 422
    // Selfguided filter is applied to a maximum stripe height of 64 + 3 pixels
    // of padding above and below
    int32_t A_[70 /*(64 + 3 + 3)*/ * REST_UNIT_STRIDE];
    int32_t *A = A_ + 3 * REST_UNIT_STRIDE + 3;
423 424
    // By inverting A and B after the boxsums, B can be of size coef instead
    // of int32_t
425 426
    coef B_[70 /*(64 + 3 + 3)*/ * REST_UNIT_STRIDE];
    coef *B = B_ + 3 * REST_UNIT_STRIDE + 3;
427 428 429

    const int step = (n == 25) + 1;
    if (n == 25) {
430 431
        boxsum5(B_, src, w + 6, h + 6);
        boxsum5sqr(A_, src, w + 6, h + 6);
432
    } else {
433 434
        boxsum3(B_, src, w + 6, h + 6);
        boxsum3sqr(A_, src, w + 6, h + 6);
435
    }
Ronald S. Bultje's avatar
Ronald S. Bultje committed
436
    const int bitdepth_min_8 = bitdepth_from_max(bitdepth_max) - 8;
437

438 439
    int32_t *AA = A - REST_UNIT_STRIDE;
    coef *BB = B - REST_UNIT_STRIDE;
440 441 442
    for (int j = -1; j < h + 1; j+= step) {
        for (int i = -1; i < w + 1; i++) {
            const int a =
Ronald S. Bultje's avatar
Ronald S. Bultje committed
443
                (AA[i] + ((1 << (2 * bitdepth_min_8)) >> 1)) >> (2 * bitdepth_min_8);
444
            const int b =
Ronald S. Bultje's avatar
Ronald S. Bultje committed
445
                (BB[i] + ((1 << bitdepth_min_8) >> 1)) >> bitdepth_min_8;
446

447
            const unsigned p = imax(a * n - b * b, 0);
Ronald S. Bultje's avatar
Ronald S. Bultje committed
448
            const unsigned z = (p * s + (1 << 19)) >> 20;
Henrik Gramner's avatar
Henrik Gramner committed
449
            const unsigned x = dav1d_sgr_x_by_x[imin(z, 255)];
450 451

            // This is where we invert A and B, so that B is of size coef.
Henrik Gramner's avatar
Henrik Gramner committed
452 453
            AA[i] = (x * BB[i] * sgr_one_by_x + (1 << 11)) >> 12;
            BB[i] = 256 - x;
454
        }
455 456
        AA += step * REST_UNIT_STRIDE;
        BB += step * REST_UNIT_STRIDE;
457 458
    }

459
    src += 3 * REST_UNIT_STRIDE + 3;
460
    if (n == 25) {
461
        int j = 0;
462 463 464 465
#define SIX_NEIGHBORS(P, i)\
    ((P[i - REST_UNIT_STRIDE]     + P[i + REST_UNIT_STRIDE]) * 6 +   \
     (P[i - 1 - REST_UNIT_STRIDE] + P[i - 1 + REST_UNIT_STRIDE] +    \
      P[i + 1 - REST_UNIT_STRIDE] + P[i + 1 + REST_UNIT_STRIDE]) * 5)
466
        for (; j < h - 1; j+=2) {
467
            for (int i = 0; i < w; i++) {
Ronald S. Bultje's avatar
Ronald S. Bultje committed
468 469
                const int a = SIX_NEIGHBORS(B, i);
                const int b = SIX_NEIGHBORS(A, i);
470 471
                dst[i] = (a * src[i] + b + (1 << 8)) >> 9;
            }
472 473 474 475
            dst += 384 /* Maximum restoration width is 384 (256 * 1.5) */;
            src += REST_UNIT_STRIDE;
            B += REST_UNIT_STRIDE;
            A += REST_UNIT_STRIDE;
476
            for (int i = 0; i < w; i++) {
Ronald S. Bultje's avatar
Ronald S. Bultje committed
477 478
                const int a = B[i] * 6 + (B[i - 1] + B[i + 1]) * 5;
                const int b = A[i] * 6 + (A[i - 1] + A[i + 1]) * 5;
479 480
                dst[i] = (a * src[i] + b + (1 << 7)) >> 8;
            }
481 482 483 484
            dst += 384 /* Maximum restoration width is 384 (256 * 1.5) */;
            src += REST_UNIT_STRIDE;
            B += REST_UNIT_STRIDE;
            A += REST_UNIT_STRIDE;
485
        }
486 487
        if (j + 1 == h) { // Last row, when number of rows is odd
            for (int i = 0; i < w; i++) {
Ronald S. Bultje's avatar
Ronald S. Bultje committed
488 489
                const int a = SIX_NEIGHBORS(B, i);
                const int b = SIX_NEIGHBORS(A, i);
490 491 492 493
                dst[i] = (a * src[i] + b + (1 << 8)) >> 9;
            }
        }
#undef SIX_NEIGHBORS
494
    } else {
495 496 497 498
#define EIGHT_NEIGHBORS(P, i)\
    ((P[i] + P[i - 1] + P[i + 1] + P[i - REST_UNIT_STRIDE] + P[i + REST_UNIT_STRIDE]) * 4 + \
     (P[i - 1 - REST_UNIT_STRIDE] + P[i - 1 + REST_UNIT_STRIDE] +                           \
      P[i + 1 - REST_UNIT_STRIDE] + P[i + 1 + REST_UNIT_STRIDE]) * 3)
499 500
        for (int j = 0; j < h; j++) {
            for (int i = 0; i < w; i++) {
Ronald S. Bultje's avatar
Ronald S. Bultje committed
501 502
                const int a = EIGHT_NEIGHBORS(B, i);
                const int b = EIGHT_NEIGHBORS(A, i);
503 504
                dst[i] = (a * src[i] + b + (1 << 8)) >> 9;
            }
505 506 507 508
            dst += 384;
            src += REST_UNIT_STRIDE;
            B += REST_UNIT_STRIDE;
            A += REST_UNIT_STRIDE;
509 510
        }
    }
511
#undef NINE_NEIGHBORS
512 513 514
}

static void selfguided_c(pixel *p, const ptrdiff_t p_stride,
515
                         const pixel (*const left)[4],
516 517
                         const pixel *lpf, const ptrdiff_t lpf_stride,
                         const int w, const int h, const int sgr_idx,
Ronald S. Bultje's avatar
Ronald S. Bultje committed
518 519
                         const int16_t sgr_w[2], const enum LrEdgeFlags edges
                         HIGHBD_DECL_SUFFIX)
520
{
521 522 523 524
    // Selfguided filter is applied to a maximum stripe height of 64 + 3 pixels
    // of padding above and below
    pixel tmp[70 /*(64 + 3 + 3)*/ * REST_UNIT_STRIDE];

525
    padding(tmp, p, p_stride, left, lpf, lpf_stride, w, h, edges);
526

527 528
    // Selfguided filter outputs to a maximum stripe height of 64 and a
    // maximum restoration width of 384 (256 * 1.5)
Ronald S. Bultje's avatar
Ronald S. Bultje committed
529
    coef dst[64 * 384];
530 531

    // both r1 and r0 can't be zero
532 533
    if (!dav1d_sgr_params[sgr_idx][0]) {
        const int s1 = dav1d_sgr_params[sgr_idx][3];
Ronald S. Bultje's avatar
Ronald S. Bultje committed
534
        selfguided_filter(dst, tmp, REST_UNIT_STRIDE, w, h, 9, s1 HIGHBD_TAIL_SUFFIX);
535 536 537
        const int w1 = (1 << 7) - sgr_w[1];
        for (int j = 0; j < h; j++) {
            for (int i = 0; i < w; i++) {
Ronald S. Bultje's avatar
Ronald S. Bultje committed
538 539
                const int u = (p[i] << 4);
                const int v = (u << 7) + w1 * (dst[j * 384 + i] - u);
540 541 542 543
                p[i] = iclip_pixel((v + (1 << 10)) >> 11);
            }
            p += PXSTRIDE(p_stride);
        }
544 545
    } else if (!dav1d_sgr_params[sgr_idx][1]) {
        const int s0 = dav1d_sgr_params[sgr_idx][2];
Ronald S. Bultje's avatar
Ronald S. Bultje committed
546
        selfguided_filter(dst, tmp, REST_UNIT_STRIDE, w, h, 25, s0 HIGHBD_TAIL_SUFFIX);
547 548 549
        const int w0 = sgr_w[0];
        for (int j = 0; j < h; j++) {
            for (int i = 0; i < w; i++) {
Ronald S. Bultje's avatar
Ronald S. Bultje committed
550 551
                const int u = (p[i] << 4);
                const int v = (u << 7) + w0 * (dst[j * 384 + i] - u);
552 553 554 555 556
                p[i] = iclip_pixel((v + (1 << 10)) >> 11);
            }
            p += PXSTRIDE(p_stride);
        }
    } else {
Ronald S. Bultje's avatar
Ronald S. Bultje committed
557
        coef dst1[64 * 384];
558 559
        const int s0 = dav1d_sgr_params[sgr_idx][2];
        const int s1 = dav1d_sgr_params[sgr_idx][3];
560 561
        const int w0 = sgr_w[0];
        const int w1 = (1 << 7) - w0 - sgr_w[1];
Ronald S. Bultje's avatar
Ronald S. Bultje committed
562 563
        selfguided_filter(dst, tmp, REST_UNIT_STRIDE, w, h, 25, s0 HIGHBD_TAIL_SUFFIX);
        selfguided_filter(dst1, tmp, REST_UNIT_STRIDE, w, h, 9, s1 HIGHBD_TAIL_SUFFIX);
564 565
        for (int j = 0; j < h; j++) {
            for (int i = 0; i < w; i++) {
Ronald S. Bultje's avatar
Ronald S. Bultje committed
566 567 568
                const int u = (p[i] << 4);
                const int v = (u << 7) + w0 * (dst[j * 384 + i] - u) +
                              w1 * (dst1[j * 384 + i] - u);
569 570 571 572 573 574 575 576 577 578
                p[i] = iclip_pixel((v + (1 << 10)) >> 11);
            }
            p += PXSTRIDE(p_stride);
        }
    }
}

void bitfn(dav1d_loop_restoration_dsp_init)(Dav1dLoopRestorationDSPContext *const c) {
    c->wiener = wiener_c;
    c->selfguided = selfguided_c;
579

580 581 582 583
#if HAVE_ASM
#if ARCH_AARCH64 || ARCH_ARM
    bitfn(dav1d_loop_restoration_dsp_init_arm)(c);
#elif ARCH_X86
584 585
    bitfn(dav1d_loop_restoration_dsp_init_x86)(c);
#endif
586
#endif
587
}