restoration.c 69.1 KB
Newer Older
1
/*
Yaowu Xu's avatar
Yaowu Xu committed
2
3
4
5
6
7
8
9
 * Copyright (c) 2016, Alliance for Open Media. All rights reserved
 *
 * This source code is subject to the terms of the BSD 2 Clause License and
 * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
 * was not distributed with this source code in the LICENSE file, you can
 * obtain it at www.aomedia.org/license/software. If the Alliance for Open
 * Media Patent License 1.0 was not distributed with this source code in the
 * PATENTS file, you can obtain it at www.aomedia.org/license/patent.
10
11
12
13
14
 *
 */

#include <math.h>

Yaowu Xu's avatar
Yaowu Xu committed
15
16
#include "./aom_config.h"
#include "./aom_dsp_rtcd.h"
17
#include "./aom_scale_rtcd.h"
18
#include "av1/common/onyxc_int.h"
19
20
21
#if CONFIG_FRAME_SUPERRES
#include "av1/common/resize.h"
#endif
22
#include "av1/common/restoration.h"
Yaowu Xu's avatar
Yaowu Xu committed
23
24
#include "aom_dsp/aom_dsp_common.h"
#include "aom_mem/aom_mem.h"
25

26
#include "aom_ports/mem.h"
27

28
const sgr_params_type sgr_params[SGRPROJ_PARAMS] = {
29
30
31
32
33
34
35
#if USE_HIGHPASS_IN_SGRPROJ
  // corner, edge, r2, eps2
  { -1, 2, 1, 1 }, { -1, 2, 1, 2 }, { -1, 2, 1, 3 }, { -1, 2, 1, 4 },
  { -1, 2, 1, 5 }, { -2, 3, 1, 2 }, { -2, 3, 1, 3 }, { -2, 3, 1, 4 },
  { -2, 3, 1, 5 }, { -2, 3, 1, 6 }, { -3, 4, 1, 3 }, { -3, 4, 1, 4 },
  { -3, 4, 1, 5 }, { -3, 4, 1, 6 }, { -3, 4, 1, 7 }, { -3, 4, 1, 8 }
#else
36
// r1, eps1, r2, eps2
37
#if MAX_RADIUS == 2
38
39
  { 2, 12, 1, 4 },  { 2, 15, 1, 6 },  { 2, 18, 1, 8 },  { 2, 20, 1, 9 },
  { 2, 22, 1, 10 }, { 2, 25, 1, 11 }, { 2, 35, 1, 12 }, { 2, 45, 1, 13 },
40
  { 2, 55, 1, 14 }, { 2, 65, 1, 15 }, { 2, 75, 1, 16 }, { 2, 30, 1, 6 },
41
42
  { 2, 50, 1, 12 }, { 2, 60, 1, 13 }, { 2, 70, 1, 14 }, { 2, 80, 1, 15 },
#else
43
44
45
46
  { 2, 12, 1, 4 },  { 2, 15, 1, 6 },  { 2, 18, 1, 8 },  { 2, 20, 1, 9 },
  { 2, 22, 1, 10 }, { 2, 25, 1, 11 }, { 2, 35, 1, 12 }, { 2, 45, 1, 13 },
  { 2, 55, 1, 14 }, { 2, 65, 1, 15 }, { 2, 75, 1, 16 }, { 3, 30, 1, 10 },
  { 3, 50, 1, 12 }, { 3, 50, 2, 25 }, { 3, 60, 2, 35 }, { 3, 70, 2, 45 },
47
#endif  // MAX_RADIUS == 2
48
#endif
49
50
};

51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#if CONFIG_MAX_TILE
static void tile_width_and_height(const AV1_COMMON *cm, int is_uv, int sb_w,
                                  int sb_h, int *px_w, int *px_h) {
  const int scaled_sb_w = sb_w << MAX_MIB_SIZE_LOG2;
  const int scaled_sb_h = sb_h << MAX_MIB_SIZE_LOG2;

  const int ss_x = is_uv && cm->subsampling_x;
  const int ss_y = is_uv && cm->subsampling_y;

  *px_w = (scaled_sb_w + ss_x) >> ss_x;
  *px_h = (scaled_sb_h + ss_y) >> ss_y;
#if CONFIG_FRAME_SUPERRES
  if (!av1_superres_unscaled(cm)) {
    av1_calculate_unscaled_superres_size(px_w, px_h,
                                         cm->superres_scale_denominator);
  }
#endif  // CONFIG_FRAME_SUPERRES
}
#endif  // CONFIG_MAX_TILE

// Count horizontal or vertical units per tile (use a width or height for
// tile_size, respectively). We basically want to divide the tile size by the
// size of a restoration unit. Rather than rounding up unconditionally as you
// might expect, we round to nearest, which models the way a right or bottom
// restoration unit can extend to up to 150% its normal width or height. The
// max with 1 is to deal with tiles that are smaller than half of a restoration
// unit.
static int count_units_in_tile(int unit_size, int tile_size) {
  return AOMMAX((tile_size + (unit_size >> 1)) / unit_size, 1);
}

void av1_alloc_restoration_struct(AV1_COMMON *cm, RestorationInfo *rsi,
                                  int is_uv) {
#if CONFIG_MAX_TILE
  // We need to allocate enough space for restoration units to cover the
  // largest tile. Without CONFIG_MAX_TILE, this is always the tile at the
  // top-left and we can use av1_get_tile_rect. With CONFIG_MAX_TILE, we have
  // to do the computation ourselves, iterating over the tiles and keeping
  // track of the largest width and height, then upscaling.
  int max_sb_w = 0;
  int max_sb_h = 0;
  for (int i = 0; i < cm->tile_cols; ++i) {
    const int sb_w = cm->tile_col_start_sb[i + 1] - cm->tile_col_start_sb[i];
    max_sb_w = AOMMAX(max_sb_w, sb_w);
  }
  for (int i = 0; i < cm->tile_rows; ++i) {
    const int sb_h = cm->tile_row_start_sb[i + 1] - cm->tile_row_start_sb[i];
    max_sb_h = AOMMAX(max_sb_h, sb_h);
  }

  int max_tile_w, max_tile_h;
  tile_width_and_height(cm, is_uv, max_sb_w, max_sb_h, &max_tile_w,
                        &max_tile_h);
#else
  TileInfo tile_info;
  av1_tile_init(&tile_info, cm, 0, 0);

  const AV1PixelRect tile_rect = av1_get_tile_rect(&tile_info, cm, is_uv);
  assert(tile_rect.left == 0 && tile_rect.top == 0);

  const int max_tile_w = tile_rect.right;
  const int max_tile_h = tile_rect.bottom;
#endif  // CONFIG_MAX_TILE

  // To calculate hpertile and vpertile (horizontal and vertical units per
  // tile), we basically want to divide the largest tile width or height by the
  // size of a restoration unit. Rather than rounding up unconditionally as you
  // might expect, we round to nearest, which models the way a right or bottom
  // restoration unit can extend to up to 150% its normal width or height. The
  // max with 1 is to deal with tiles that are smaller than half of a
  // restoration unit.
  const int unit_size = rsi->restoration_unit_size;
  const int hpertile = count_units_in_tile(unit_size, max_tile_w);
  const int vpertile = count_units_in_tile(unit_size, max_tile_h);

  rsi->units_per_tile = hpertile * vpertile;
  rsi->horz_units_per_tile = hpertile;
  rsi->vert_units_per_tile = vpertile;

  const int ntiles = cm->tile_rows * cm->tile_cols;
  const int nunits = ntiles * rsi->units_per_tile;

  aom_free(rsi->unit_info);
  CHECK_MEM_ERROR(cm, rsi->unit_info, (RestorationUnitInfo *)aom_malloc(
                                          sizeof(*rsi->unit_info) * nunits));
136
137
138
}

void av1_free_restoration_struct(RestorationInfo *rst_info) {
139
140
  aom_free(rst_info->unit_info);
  rst_info->unit_info = NULL;
141
}
142
143
144

// TODO(debargha): This table can be substantially reduced since only a few
// values are actually used.
David Barker's avatar
David Barker committed
145
int sgrproj_mtable[MAX_EPS][MAX_NELEM];
146
147
148
149
150
151
152
153
154
155

static void GenSgrprojVtable() {
  int e, n;
  for (e = 1; e <= MAX_EPS; ++e)
    for (n = 1; n <= MAX_NELEM; ++n) {
      const int n2e = n * n * e;
      sgrproj_mtable[e - 1][n - 1] =
          (((1 << SGRPROJ_MTABLE_BITS) + n2e / 2) / n2e);
    }
}
156
157

void av1_loop_restoration_precal() { GenSgrprojVtable(); }
158

159
160
static void extend_frame_lowbd(uint8_t *data, int width, int height, int stride,
                               int border_horz, int border_vert) {
161
162
163
164
  uint8_t *data_p;
  int i;
  for (i = 0; i < height; ++i) {
    data_p = data + i * stride;
165
166
    memset(data_p - border_horz, data_p[0], border_horz);
    memset(data_p + width, data_p[width - 1], border_horz);
167
  }
168
169
170
  data_p = data - border_horz;
  for (i = -border_vert; i < 0; ++i) {
    memcpy(data_p + i * stride, data_p, width + 2 * border_horz);
171
  }
172
  for (i = height; i < height + border_vert; ++i) {
173
    memcpy(data_p + i * stride, data_p + (height - 1) * stride,
174
           width + 2 * border_horz);
175
176
177
  }
}

178
179
180
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
211
212
213
#if CONFIG_HIGHBITDEPTH
static void extend_frame_highbd(uint16_t *data, int width, int height,
                                int stride, int border_horz, int border_vert) {
  uint16_t *data_p;
  int i, j;
  for (i = 0; i < height; ++i) {
    data_p = data + i * stride;
    for (j = -border_horz; j < 0; ++j) data_p[j] = data_p[0];
    for (j = width; j < width + border_horz; ++j) data_p[j] = data_p[width - 1];
  }
  data_p = data - border_horz;
  for (i = -border_vert; i < 0; ++i) {
    memcpy(data_p + i * stride, data_p,
           (width + 2 * border_horz) * sizeof(uint16_t));
  }
  for (i = height; i < height + border_vert; ++i) {
    memcpy(data_p + i * stride, data_p + (height - 1) * stride,
           (width + 2 * border_horz) * sizeof(uint16_t));
  }
}
#endif

void extend_frame(uint8_t *data, int width, int height, int stride,
                  int border_horz, int border_vert, int highbd) {
#if !CONFIG_HIGHBITDEPTH
  assert(highbd == 0);
  (void)highbd;
#else
  if (highbd)
    extend_frame_highbd(CONVERT_TO_SHORTPTR(data), width, height, stride,
                        border_horz, border_vert);
  else
#endif
  extend_frame_lowbd(data, width, height, stride, border_horz, border_vert);
}

214
215
216
217
static void copy_tile_lowbd(int width, int height, const uint8_t *src,
                            int src_stride, uint8_t *dst, int dst_stride) {
  for (int i = 0; i < height; ++i)
    memcpy(dst + i * dst_stride, src + i * src_stride, width);
218
219
220
}

#if CONFIG_HIGHBITDEPTH
221
222
223
224
static void copy_tile_highbd(int width, int height, const uint16_t *src,
                             int src_stride, uint16_t *dst, int dst_stride) {
  for (int i = 0; i < height; ++i)
    memcpy(dst + i * dst_stride, src + i * src_stride, width * sizeof(*dst));
225
226
227
}
#endif

228
229
static void copy_tile(int width, int height, const uint8_t *src, int src_stride,
                      uint8_t *dst, int dst_stride, int highbd) {
230
231
232
233
234
#if !CONFIG_HIGHBITDEPTH
  assert(highbd == 0);
  (void)highbd;
#else
  if (highbd)
235
    copy_tile_highbd(width, height, CONVERT_TO_SHORTPTR(src), src_stride,
236
237
238
                     CONVERT_TO_SHORTPTR(dst), dst_stride);
  else
#endif
239
  copy_tile_lowbd(width, height, src, src_stride, dst, dst_stride);
240
}
241

242
243
244
#if CONFIG_STRIPED_LOOP_RESTORATION
#define REAL_PTR(hbd, d) ((hbd) ? (uint8_t *)CONVERT_TO_SHORTPTR(d) : (d))

245
246
247
248
// With striped loop restoration, the filtering for each 64-pixel stripe gets
// most of its input from the output of CDEF (stored in data8), but pixels just
// above and below the stripe come straight from the deblocker. These have been
// stored away in separate buffers.
249
//
250
251
252
253
254
255
// This function modifies data8 (which was the output from CDEF) by copying in
// the boundary pixels. Before doing so, it saves the pixels that get
// overwritten into a temporary buffer. They will be restored again by
// restore_processing_stripe_boundary.
//
// limits gives the rectangular limits of the remaining stripes for the current
256
// restoration unit. rsb is the stored stripe boundaries (the saved output from
257
258
259
260
// the deblocker). ss_y is true if we're on a chroma plane with vertical
// subsampling. use_highbd is true if the data has 2 bytes per pixel. rlbs
// contain scratch buffers to hold the CDEF data (written back to the frame by
// restore_processing_stripe_boundary)
261
262
static int setup_processing_stripe_boundary(
    const RestorationTileLimits *limits, const RestorationStripeBoundaries *rsb,
263
    int ss_y, int use_highbd, uint8_t *data8, int stride,
264
    RestorationLineBuffers *rlbs) {
265
266
267
  // Which stripe is this? limits->v_start is the top of the stripe in pixel
  // units, but we add tile_offset to get the number of pixels from the top of
  // the first stripe, which lies off the image.
268
  const int tile_offset = RESTORATION_TILE_OFFSET >> ss_y;
269
  const int stripe_height = RESTORATION_PROC_UNIT_SIZE >> ss_y;
270
271
272
273
274
275
  const int stripe_index = (limits->v_start + tile_offset) / stripe_height;

  // Horizontal offsets within the line buffers. The buffer logically starts at
  // column -RESTORATION_EXTRA_HORZ. We'll start our copy from the column
  // limits->h_start - RESTORATION_EXTRA_HORZ and copy up to the column
  // limits->h_end + RESTORATION_EXTRA_HORZ.
276
  const int buf_stride = rsb->stripe_boundary_stride;
277
278
279
280
281
  const int buf_x0_off = limits->h_start;
  const int line_width =
      (limits->h_end - limits->h_start) + 2 * RESTORATION_EXTRA_HORZ;
  const int line_size = line_width << use_highbd;
  const int data_x0_off = limits->h_start - RESTORATION_EXTRA_HORZ;
282

283
284
  assert(CONFIG_HIGHBITDEPTH || !use_highbd);

285
286
287
288
289
  // Replace RESTORATION_BORDER pixels above the top of the stripe, unless this
  // is the top of the image. We expand 2 lines from rsb->stripe_boundary_above
  // to fill 3 lines of above pixels. This is done by duplicating the topmost
  // of the 2 lines (see the AOMMAX call when calculating src_row, which gets
  // the values 0, 0, 1 for i = -3, -2, -1).
290
291
  if (stripe_index > 0) {
    const int above_buf_y = 2 * (stripe_index - 1);
292
    uint8_t *data8_tl = data8 + limits->v_start * stride + data_x0_off;
293

294
295
    for (int i = -RESTORATION_BORDER; i < 0; ++i) {
      const int src_row = AOMMAX(i + RESTORATION_CTX_VERT, 0);
296
      const int buf_off = buf_x0_off + (above_buf_y + src_row) * buf_stride;
297
      const uint8_t *src = rsb->stripe_boundary_above + (buf_off << use_highbd);
298
      uint8_t *dst8 = data8_tl + i * stride;
299
300
301
      // Save old pixels, then replace with data from stripe_boundary_above
      memcpy(rlbs->tmp_save_above[i + RESTORATION_BORDER],
             REAL_PTR(use_highbd, dst8), line_size);
302
      memcpy(REAL_PTR(use_highbd, dst8), src, line_size);
303
304
    }
  }
305

306
307
308
309
310
311
312
313
314
315
316
  // Replace up to RESTORATION_BORDER pixels below the bottom of the
  // stripe. When replacing the maximum amount, the second buffer row is
  // repeated, so src_row gets the values 0, 1, 1 for i = 0, 1, 2.
  //
  // We might not write that many rows if the stripe isn't of full height
  // (which might happen at the bottom of a restoration unit).
  const int full_stripe_bottom =
      stripe_height * (1 + stripe_index) - tile_offset;
  const int stripe_bottom = AOMMIN(limits->v_end, full_stripe_bottom);
  const int rows_needed_below =
      RESTORATION_BORDER - (full_stripe_bottom - stripe_bottom);
317

318
319
  const int below_buf_y = RESTORATION_CTX_VERT * stripe_index;
  uint8_t *data8_bl = data8 + full_stripe_bottom * stride + data_x0_off;
320
321

  for (int i = 0; i < rows_needed_below; ++i) {
322
    const int src_row = AOMMIN(i, RESTORATION_CTX_VERT - 1);
323
    const int buf_off = buf_x0_off + (below_buf_y + src_row) * buf_stride;
324
    const uint8_t *src = rsb->stripe_boundary_below + (buf_off << use_highbd);
325
    uint8_t *dst8 = data8_bl + i * stride;
326
    // Save old pixels, then replace with data from stripe_boundary_below
327
    memcpy(rlbs->tmp_save_below[i], REAL_PTR(use_highbd, dst8), line_size);
328
    memcpy(REAL_PTR(use_highbd, dst8), src, line_size);
329
  }
330

331
  // Finally, return the actual height of this stripe.
332
  return AOMMIN(limits->v_end, full_stripe_bottom) - limits->v_start;
333
334
335
336
}

// This function restores the boundary lines modified by
// setup_processing_stripe_boundary.
337
static void restore_processing_stripe_boundary(
338
    const RestorationTileLimits *limits, const RestorationLineBuffers *rlbs,
339
    int ss_y, int use_highbd, uint8_t *data8, int stride) {
340
  const int tile_offset = RESTORATION_TILE_OFFSET >> ss_y;
341
  const int stripe_height = RESTORATION_PROC_UNIT_SIZE >> ss_y;
342
343
344
345
346
347
  const int stripe_index = (limits->v_start + tile_offset) / stripe_height;

  const int line_width =
      (limits->h_end - limits->h_start) + 2 * RESTORATION_EXTRA_HORZ;
  const int line_size = line_width << use_highbd;
  const int data_x0_off = limits->h_start - RESTORATION_EXTRA_HORZ;
348

349
350
  assert(CONFIG_HIGHBITDEPTH || !use_highbd);

351
  if (stripe_index > 0) {
352
353
    uint8_t *data8_tl = data8 + limits->v_start * stride + data_x0_off;
    for (int i = -RESTORATION_BORDER; i < 0; ++i) {
354
355
      uint8_t *dst8 = data8_tl + i * stride;
      // Save old pixels, then replace with data from boundary_above_buf
356
357
      memcpy(REAL_PTR(use_highbd, dst8),
             rlbs->tmp_save_above[i + RESTORATION_BORDER], line_size);
358
359
    }
  }
360

361
362
363
364
365
  const int full_stripe_bottom =
      stripe_height * (1 + stripe_index) - tile_offset;
  const int stripe_bottom = AOMMIN(limits->v_end, full_stripe_bottom);
  const int rows_needed_below =
      RESTORATION_BORDER - (full_stripe_bottom - stripe_bottom);
366
367
368
369
370

  uint8_t *data8_bl = data8 + stripe_bottom * stride + data_x0_off;
  for (int i = 0; i < rows_needed_below; ++i) {
    uint8_t *dst8 = data8_bl + i * stride;
    // Save old pixels, then replace with data from boundary_below_buf
371
    memcpy(REAL_PTR(use_highbd, dst8), rlbs->tmp_save_below[i], line_size);
372
373
374
375
  }
}
#endif

376
377
#if USE_WIENER_HIGH_INTERMEDIATE_PRECISION
#define wiener_convolve8_add_src aom_convolve8_add_src_hip
378
#else
379
#define wiener_convolve8_add_src aom_convolve8_add_src
380
381
#endif

382
383
384
385
386
387
388
389
390
391
392
static void wiener_filter_stripe(const RestorationUnitInfo *rui,
                                 int stripe_width, int stripe_height,
                                 int procunit_width, const uint8_t *src,
                                 int src_stride, uint8_t *dst, int dst_stride,
                                 int32_t *tmpbuf, int bit_depth) {
  (void)tmpbuf;
  (void)bit_depth;
  assert(bit_depth == 8);

  for (int j = 0; j < stripe_width; j += procunit_width) {
    int w = AOMMIN(procunit_width, (stripe_width - j + 15) & ~15);
393
394
395
    const uint8_t *src_p = src + j;
    uint8_t *dst_p = dst + j;
    wiener_convolve8_add_src(src_p, src_stride, dst_p, dst_stride,
396
                             rui->wiener_info.hfilter, 16,
397
                             rui->wiener_info.vfilter, 16, w, stripe_height);
398
  }
399
}
400

401
402
/* Calculate windowed sums (if sqr=0) or sums of squares (if sqr=1)
   over the input. The window is of size (2r + 1)x(2r + 1), and we
403
   specialize to r = 1, 2, 3. A default function is used for r > 3.
404
405
406
407
408
409
410
411
412
413
414
415
416
417

   Each loop follows the same format: We keep a window's worth of input
   in individual variables and select data out of that as appropriate.
*/
static void boxsum1(int32_t *src, int width, int height, int src_stride,
                    int sqr, int32_t *dst, int dst_stride) {
  int i, j, a, b, c;

  // Vertical sum over 3-pixel regions, from src into dst.
  if (!sqr) {
    for (j = 0; j < width; ++j) {
      a = src[j];
      b = src[src_stride + j];
      c = src[2 * src_stride + j];
418

419
420
421
422
423
424
425
426
427
428
429
430
431
432
      dst[j] = a + b;
      for (i = 1; i < height - 2; ++i) {
        // Loop invariant: At the start of each iteration,
        // a = src[(i - 1) * src_stride + j]
        // b = src[(i    ) * src_stride + j]
        // c = src[(i + 1) * src_stride + j]
        dst[i * dst_stride + j] = a + b + c;
        a = b;
        b = c;
        c = src[(i + 2) * src_stride + j];
      }
      dst[i * dst_stride + j] = a + b + c;
      dst[(i + 1) * dst_stride + j] = b + c;
    }
433
  } else {
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
    for (j = 0; j < width; ++j) {
      a = src[j] * src[j];
      b = src[src_stride + j] * src[src_stride + j];
      c = src[2 * src_stride + j] * src[2 * src_stride + j];

      dst[j] = a + b;
      for (i = 1; i < height - 2; ++i) {
        dst[i * dst_stride + j] = a + b + c;
        a = b;
        b = c;
        c = src[(i + 2) * src_stride + j] * src[(i + 2) * src_stride + j];
      }
      dst[i * dst_stride + j] = a + b + c;
      dst[(i + 1) * dst_stride + j] = b + c;
    }
  }

  // Horizontal sum over 3-pixel regions of dst
  for (i = 0; i < height; ++i) {
    a = dst[i * dst_stride];
    b = dst[i * dst_stride + 1];
    c = dst[i * dst_stride + 2];

    dst[i * dst_stride] = a + b;
    for (j = 1; j < width - 2; ++j) {
      // Loop invariant: At the start of each iteration,
      // a = src[i * src_stride + (j - 1)]
      // b = src[i * src_stride + (j    )]
      // c = src[i * src_stride + (j + 1)]
      dst[i * dst_stride + j] = a + b + c;
      a = b;
      b = c;
      c = dst[i * dst_stride + (j + 2)];
    }
    dst[i * dst_stride + j] = a + b + c;
    dst[i * dst_stride + (j + 1)] = b + c;
  }
}

static void boxsum2(int32_t *src, int width, int height, int src_stride,
                    int sqr, int32_t *dst, int dst_stride) {
  int i, j, a, b, c, d, e;

  // Vertical sum over 5-pixel regions, from src into dst.
  if (!sqr) {
    for (j = 0; j < width; ++j) {
      a = src[j];
      b = src[src_stride + j];
      c = src[2 * src_stride + j];
      d = src[3 * src_stride + j];
      e = src[4 * src_stride + j];

      dst[j] = a + b + c;
      dst[dst_stride + j] = a + b + c + d;
      for (i = 2; i < height - 3; ++i) {
        // Loop invariant: At the start of each iteration,
        // a = src[(i - 2) * src_stride + j]
        // b = src[(i - 1) * src_stride + j]
        // c = src[(i    ) * src_stride + j]
        // d = src[(i + 1) * src_stride + j]
        // e = src[(i + 2) * src_stride + j]
        dst[i * dst_stride + j] = a + b + c + d + e;
        a = b;
        b = c;
        c = d;
        d = e;
        e = src[(i + 3) * src_stride + j];
      }
      dst[i * dst_stride + j] = a + b + c + d + e;
      dst[(i + 1) * dst_stride + j] = b + c + d + e;
      dst[(i + 2) * dst_stride + j] = c + d + e;
    }
  } else {
    for (j = 0; j < width; ++j) {
      a = src[j] * src[j];
      b = src[src_stride + j] * src[src_stride + j];
      c = src[2 * src_stride + j] * src[2 * src_stride + j];
      d = src[3 * src_stride + j] * src[3 * src_stride + j];
      e = src[4 * src_stride + j] * src[4 * src_stride + j];

      dst[j] = a + b + c;
      dst[dst_stride + j] = a + b + c + d;
      for (i = 2; i < height - 3; ++i) {
        dst[i * dst_stride + j] = a + b + c + d + e;
        a = b;
        b = c;
        c = d;
        d = e;
        e = src[(i + 3) * src_stride + j] * src[(i + 3) * src_stride + j];
      }
      dst[i * dst_stride + j] = a + b + c + d + e;
      dst[(i + 1) * dst_stride + j] = b + c + d + e;
      dst[(i + 2) * dst_stride + j] = c + d + e;
    }
  }

  // Horizontal sum over 5-pixel regions of dst
  for (i = 0; i < height; ++i) {
    a = dst[i * dst_stride];
    b = dst[i * dst_stride + 1];
    c = dst[i * dst_stride + 2];
    d = dst[i * dst_stride + 3];
    e = dst[i * dst_stride + 4];

    dst[i * dst_stride] = a + b + c;
    dst[i * dst_stride + 1] = a + b + c + d;
    for (j = 2; j < width - 3; ++j) {
      // Loop invariant: At the start of each iteration,
      // a = src[i * src_stride + (j - 2)]
      // b = src[i * src_stride + (j - 1)]
      // c = src[i * src_stride + (j    )]
      // d = src[i * src_stride + (j + 1)]
      // e = src[i * src_stride + (j + 2)]
      dst[i * dst_stride + j] = a + b + c + d + e;
      a = b;
      b = c;
      c = d;
      d = e;
      e = dst[i * dst_stride + (j + 3)];
    }
    dst[i * dst_stride + j] = a + b + c + d + e;
    dst[i * dst_stride + (j + 1)] = b + c + d + e;
    dst[i * dst_stride + (j + 2)] = c + d + e;
  }
}

560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
static void boxsum3(int32_t *src, int width, int height, int src_stride,
                    int sqr, int32_t *dst, int dst_stride) {
  int i, j, a, b, c, d, e, f, g;

  // Vertical sum over 7-pixel regions, from src into dst.
  if (!sqr) {
    for (j = 0; j < width; ++j) {
      a = src[j];
      b = src[1 * src_stride + j];
      c = src[2 * src_stride + j];
      d = src[3 * src_stride + j];
      e = src[4 * src_stride + j];
      f = src[5 * src_stride + j];
      g = src[6 * src_stride + j];

      dst[j] = a + b + c + d;
      dst[dst_stride + j] = a + b + c + d + e;
      dst[2 * dst_stride + j] = a + b + c + d + e + f;
      for (i = 3; i < height - 4; ++i) {
        dst[i * dst_stride + j] = a + b + c + d + e + f + g;
        a = b;
        b = c;
        c = d;
        d = e;
        e = f;
        f = g;
        g = src[(i + 4) * src_stride + j];
      }
      dst[i * dst_stride + j] = a + b + c + d + e + f + g;
      dst[(i + 1) * dst_stride + j] = b + c + d + e + f + g;
      dst[(i + 2) * dst_stride + j] = c + d + e + f + g;
      dst[(i + 3) * dst_stride + j] = d + e + f + g;
    }
  } else {
    for (j = 0; j < width; ++j) {
      a = src[j] * src[j];
      b = src[1 * src_stride + j] * src[1 * src_stride + j];
      c = src[2 * src_stride + j] * src[2 * src_stride + j];
      d = src[3 * src_stride + j] * src[3 * src_stride + j];
      e = src[4 * src_stride + j] * src[4 * src_stride + j];
      f = src[5 * src_stride + j] * src[5 * src_stride + j];
      g = src[6 * src_stride + j] * src[6 * src_stride + j];

      dst[j] = a + b + c + d;
      dst[dst_stride + j] = a + b + c + d + e;
      dst[2 * dst_stride + j] = a + b + c + d + e + f;
      for (i = 3; i < height - 4; ++i) {
        dst[i * dst_stride + j] = a + b + c + d + e + f + g;
        a = b;
        b = c;
        c = d;
        d = e;
        e = f;
        f = g;
        g = src[(i + 4) * src_stride + j] * src[(i + 4) * src_stride + j];
      }
      dst[i * dst_stride + j] = a + b + c + d + e + f + g;
      dst[(i + 1) * dst_stride + j] = b + c + d + e + f + g;
      dst[(i + 2) * dst_stride + j] = c + d + e + f + g;
      dst[(i + 3) * dst_stride + j] = d + e + f + g;
    }
  }

  // Horizontal sum over 7-pixel regions of dst
  for (i = 0; i < height; ++i) {
    a = dst[i * dst_stride];
    b = dst[i * dst_stride + 1];
    c = dst[i * dst_stride + 2];
    d = dst[i * dst_stride + 3];
    e = dst[i * dst_stride + 4];
    f = dst[i * dst_stride + 5];
    g = dst[i * dst_stride + 6];

    dst[i * dst_stride] = a + b + c + d;
    dst[i * dst_stride + 1] = a + b + c + d + e;
    dst[i * dst_stride + 2] = a + b + c + d + e + f;
    for (j = 3; j < width - 4; ++j) {
      dst[i * dst_stride + j] = a + b + c + d + e + f + g;
      a = b;
      b = c;
      c = d;
      d = e;
      e = f;
      f = g;
      g = dst[i * dst_stride + (j + 4)];
    }
    dst[i * dst_stride + j] = a + b + c + d + e + f + g;
    dst[i * dst_stride + (j + 1)] = b + c + d + e + f + g;
    dst[i * dst_stride + (j + 2)] = c + d + e + f + g;
    dst[i * dst_stride + (j + 3)] = d + e + f + g;
  }
}

// Generic version for any r. To be removed after experiments are done.
static void boxsumr(int32_t *src, int width, int height, int src_stride, int r,
                    int sqr, int32_t *dst, int dst_stride) {
  int32_t *tmp = aom_malloc(width * height * sizeof(*tmp));
  int tmp_stride = width;
  int i, j;
  if (sqr) {
    for (j = 0; j < width; ++j) tmp[j] = src[j] * src[j];
    for (j = 0; j < width; ++j)
      for (i = 1; i < height; ++i)
        tmp[i * tmp_stride + j] =
            tmp[(i - 1) * tmp_stride + j] +
            src[i * src_stride + j] * src[i * src_stride + j];
  } else {
    memcpy(tmp, src, sizeof(*tmp) * width);
    for (j = 0; j < width; ++j)
      for (i = 1; i < height; ++i)
        tmp[i * tmp_stride + j] =
            tmp[(i - 1) * tmp_stride + j] + src[i * src_stride + j];
  }
  for (i = 0; i <= r; ++i)
    memcpy(&dst[i * dst_stride], &tmp[(i + r) * tmp_stride],
           sizeof(*tmp) * width);
  for (i = r + 1; i < height - r; ++i)
    for (j = 0; j < width; ++j)
      dst[i * dst_stride + j] =
          tmp[(i + r) * tmp_stride + j] - tmp[(i - r - 1) * tmp_stride + j];
  for (i = height - r; i < height; ++i)
    for (j = 0; j < width; ++j)
      dst[i * dst_stride + j] = tmp[(height - 1) * tmp_stride + j] -
                                tmp[(i - r - 1) * tmp_stride + j];

  for (i = 0; i < height; ++i) tmp[i * tmp_stride] = dst[i * dst_stride];
  for (i = 0; i < height; ++i)
    for (j = 1; j < width; ++j)
      tmp[i * tmp_stride + j] =
          tmp[i * tmp_stride + j - 1] + dst[i * src_stride + j];

  for (j = 0; j <= r; ++j)
    for (i = 0; i < height; ++i)
      dst[i * dst_stride + j] = tmp[i * tmp_stride + j + r];
  for (j = r + 1; j < width - r; ++j)
    for (i = 0; i < height; ++i)
      dst[i * dst_stride + j] =
          tmp[i * tmp_stride + j + r] - tmp[i * tmp_stride + j - r - 1];
  for (j = width - r; j < width; ++j)
    for (i = 0; i < height; ++i)
      dst[i * dst_stride + j] =
          tmp[i * tmp_stride + width - 1] - tmp[i * tmp_stride + j - r - 1];
  aom_free(tmp);
}

705
706
707
708
709
710
static void boxsum(int32_t *src, int width, int height, int src_stride, int r,
                   int sqr, int32_t *dst, int dst_stride) {
  if (r == 1)
    boxsum1(src, width, height, src_stride, sqr, dst, dst_stride);
  else if (r == 2)
    boxsum2(src, width, height, src_stride, sqr, dst, dst_stride);
711
712
713
714
  else if (r == 3)
    boxsum3(src, width, height, src_stride, sqr, dst, dst_stride);
  else
    boxsumr(src, width, height, src_stride, r, sqr, dst, dst_stride);
715
716
717
718
}

static void boxnum(int width, int height, int r, int8_t *num, int num_stride) {
  int i, j;
719
720
721
  for (i = 0; i <= r; ++i) {
    for (j = 0; j <= r; ++j) {
      num[i * num_stride + j] = (r + 1 + i) * (r + 1 + j);
722
723
724
725
726
727
      num[i * num_stride + (width - 1 - j)] = num[i * num_stride + j];
      num[(height - 1 - i) * num_stride + j] = num[i * num_stride + j];
      num[(height - 1 - i) * num_stride + (width - 1 - j)] =
          num[i * num_stride + j];
    }
  }
728
729
  for (j = 0; j <= r; ++j) {
    const int val = (2 * r + 1) * (r + 1 + j);
730
731
732
733
734
    for (i = r + 1; i < height - r; ++i) {
      num[i * num_stride + j] = val;
      num[i * num_stride + (width - 1 - j)] = val;
    }
  }
735
736
  for (i = 0; i <= r; ++i) {
    const int val = (2 * r + 1) * (r + 1 + i);
737
738
739
740
741
742
743
    for (j = r + 1; j < width - r; ++j) {
      num[i * num_stride + j] = val;
      num[(height - 1 - i) * num_stride + j] = val;
    }
  }
  for (i = r + 1; i < height - r; ++i) {
    for (j = r + 1; j < width - r; ++j) {
744
      num[i * num_stride + j] = (2 * r + 1) * (2 * r + 1);
745
746
747
748
    }
  }
}

749
void decode_xq(const int *xqd, int *xq) {
750
  xq[0] = xqd[0];
751
752
753
  xq[1] = (1 << SGRPROJ_PRJ_BITS) - xq[0] - xqd[1];
}

David Barker's avatar
David Barker committed
754
const int32_t x_by_xplus1[256] = {
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
  0,   128, 171, 192, 205, 213, 219, 224, 228, 230, 233, 235, 236, 238, 239,
  240, 241, 242, 243, 243, 244, 244, 245, 245, 246, 246, 247, 247, 247, 247,
  248, 248, 248, 248, 249, 249, 249, 249, 249, 250, 250, 250, 250, 250, 250,
  250, 251, 251, 251, 251, 251, 251, 251, 251, 251, 251, 252, 252, 252, 252,
  252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 252, 253, 253,
  253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253,
  253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 253, 254, 254, 254,
  254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254,
  254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254,
  254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254,
  254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254, 254,
  254, 254, 254, 254, 254, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
  255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
  255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
  255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
  255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
  255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
  256,
};

David Barker's avatar
David Barker committed
775
const int32_t one_by_x[MAX_NELEM] = {
776
  4096, 2048, 1365, 1024, 819, 683, 585, 512, 455, 410, 372, 341, 315,
777
778
779
780
781
  293,  273,  256,  241,  228, 216, 205, 195, 186, 178, 171, 164,
#if MAX_RADIUS > 2
  158,  152,  146,  141,  137, 132, 128, 124, 120, 117, 114, 111, 108,
  105,  102,  100,  98,   95,  93,  91,  89,  87,  85,  84
#endif  // MAX_RADIUS > 2
782
783
};

784
static void av1_selfguided_restoration_internal(int32_t *dgd, int width,
785
786
                                                int height, int dgd_stride,
                                                int32_t *dst, int dst_stride,
787
                                                int bit_depth, int r, int eps) {
788
789
790
  const int width_ext = width + 2 * SGRPROJ_BORDER_HORZ;
  const int height_ext = height + 2 * SGRPROJ_BORDER_VERT;
  const int num_stride = width_ext;
David Barker's avatar
David Barker committed
791
792
793
794
  // Adjusting the stride of A and B here appears to avoid bad cache effects,
  // leading to a significant speed improvement.
  // We also align the stride to a multiple of 16 bytes, for consistency
  // with the SIMD version of this function.
795
  int buf_stride = ((width_ext + 3) & ~3) + 16;
796
797
798
799
800
  int32_t A_[RESTORATION_PROC_UNIT_PELS];
  int32_t B_[RESTORATION_PROC_UNIT_PELS];
  int32_t *A = A_;
  int32_t *B = B_;
  int8_t num_[RESTORATION_PROC_UNIT_PELS];
801
802
  int8_t *num = num_ + SGRPROJ_BORDER_VERT * num_stride + SGRPROJ_BORDER_HORZ;
  int i, j;
803

804
805
806
  // Don't filter tiles with dimensions < 5 on any axis
  if ((width < 5) || (height < 5)) return;

807
808
809
810
811
  boxsum(dgd - dgd_stride * SGRPROJ_BORDER_VERT - SGRPROJ_BORDER_HORZ,
         width_ext, height_ext, dgd_stride, r, 0, B, buf_stride);
  boxsum(dgd - dgd_stride * SGRPROJ_BORDER_VERT - SGRPROJ_BORDER_HORZ,
         width_ext, height_ext, dgd_stride, r, 1, A, buf_stride);
  boxnum(width_ext, height_ext, r, num_, num_stride);
812
  assert(r <= 3);
813
814
  A += SGRPROJ_BORDER_VERT * buf_stride + SGRPROJ_BORDER_HORZ;
  B += SGRPROJ_BORDER_VERT * buf_stride + SGRPROJ_BORDER_HORZ;
815
816
  for (i = 0; i < height; ++i) {
    for (j = 0; j < width; ++j) {
David Barker's avatar
David Barker committed
817
      const int k = i * buf_stride + j;
818
      const int n = num[i * num_stride + j];
819

820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
      // a < 2^16 * n < 2^22 regardless of bit depth
      uint32_t a = ROUND_POWER_OF_TWO(A[k], 2 * (bit_depth - 8));
      // b < 2^8 * n < 2^14 regardless of bit depth
      uint32_t b = ROUND_POWER_OF_TWO(B[k], bit_depth - 8);

      // Each term in calculating p = a * n - b * b is < 2^16 * n^2 < 2^28,
      // and p itself satisfies p < 2^14 * n^2 < 2^26.
      // Note: Sometimes, in high bit depth, we can end up with a*n < b*b.
      // This is an artefact of rounding, and can only happen if all pixels
      // are (almost) identical, so in this case we saturate to p=0.
      uint32_t p = (a * n < b * b) ? 0 : a * n - b * b;
      uint32_t s = sgrproj_mtable[eps - 1][n - 1];

      // p * s < (2^14 * n^2) * round(2^20 / n^2 eps) < 2^34 / eps < 2^32
      // as long as eps >= 4. So p * s fits into a uint32_t, and z < 2^12
      // (this holds even after accounting for the rounding in s)
      const uint32_t z = ROUND_POWER_OF_TWO(p * s, SGRPROJ_MTABLE_BITS);

      A[k] = x_by_xplus1[AOMMIN(z, 255)];  // < 2^8

      // SGRPROJ_SGR - A[k] < 2^8, B[k] < 2^(bit_depth) * n,
      // one_by_x[n - 1] = round(2^12 / n)
      // => the product here is < 2^(20 + bit_depth) <= 2^32,
      // and B[k] is set to a value < 2^(8 + bit depth)
      B[k] = (int32_t)ROUND_POWER_OF_TWO((uint32_t)(SGRPROJ_SGR - A[k]) *
                                             (uint32_t)B[k] *
                                             (uint32_t)one_by_x[n - 1],
                                         SGRPROJ_RECIP_BITS);
848
849
850
851
852
    }
  }
  i = 0;
  j = 0;
  {
David Barker's avatar
David Barker committed
853
    const int k = i * buf_stride + j;
854
855
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
856
    const int nb = 3;
857
    const int32_t a =
David Barker's avatar
David Barker committed
858
        3 * A[k] + 2 * A[k + 1] + 2 * A[k + buf_stride] + A[k + buf_stride + 1];
859
    const int32_t b =
David Barker's avatar
David Barker committed
860
        3 * B[k] + 2 * B[k + 1] + 2 * B[k + buf_stride] + B[k + buf_stride + 1];
861
    const int32_t v = a * dgd[l] + b;
862
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
863
864
865
866
  }
  i = 0;
  j = width - 1;
  {
David Barker's avatar
David Barker committed
867
    const int k = i * buf_stride + j;
868
869
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
870
    const int nb = 3;
871
    const int32_t a =
David Barker's avatar
David Barker committed
872
        3 * A[k] + 2 * A[k - 1] + 2 * A[k + buf_stride] + A[k + buf_stride - 1];
873
    const int32_t b =
David Barker's avatar
David Barker committed
874
        3 * B[k] + 2 * B[k - 1] + 2 * B[k + buf_stride] + B[k + buf_stride - 1];
875
    const int32_t v = a * dgd[l] + b;
876
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
877
878
879
880
  }
  i = height - 1;
  j = 0;
  {
David Barker's avatar
David Barker committed
881
    const int k = i * buf_stride + j;
882
883
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
884
    const int nb = 3;
885
    const int32_t a =
David Barker's avatar
David Barker committed
886
        3 * A[k] + 2 * A[k + 1] + 2 * A[k - buf_stride] + A[k - buf_stride + 1];
887
    const int32_t b =
David Barker's avatar
David Barker committed
888
        3 * B[k] + 2 * B[k + 1] + 2 * B[k - buf_stride] + B[k - buf_stride + 1];
889
    const int32_t v = a * dgd[l] + b;
890
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
891
892
893
894
  }
  i = height - 1;
  j = width - 1;
  {
David Barker's avatar
David Barker committed
895
    const int k = i * buf_stride + j;
896
897
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
898
    const int nb = 3;
899
    const int32_t a =
David Barker's avatar
David Barker committed
900
        3 * A[k] + 2 * A[k - 1] + 2 * A[k - buf_stride] + A[k - buf_stride - 1];
901
    const int32_t b =
David Barker's avatar
David Barker committed
902
        3 * B[k] + 2 * B[k - 1] + 2 * B[k - buf_stride] + B[k - buf_stride - 1];
903
    const int32_t v = a * dgd[l] + b;
904
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
905
906
907
  }
  i = 0;
  for (j = 1; j < width - 1; ++j) {
David Barker's avatar
David Barker committed
908
    const int k = i * buf_stride + j;
909
910
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
911
    const int nb = 3;
David Barker's avatar
David Barker committed
912
913
914
915
    const int32_t a = A[k] + 2 * (A[k - 1] + A[k + 1]) + A[k + buf_stride] +
                      A[k + buf_stride - 1] + A[k + buf_stride + 1];
    const int32_t b = B[k] + 2 * (B[k - 1] + B[k + 1]) + B[k + buf_stride] +
                      B[k + buf_stride - 1] + B[k + buf_stride + 1];
916
    const int32_t v = a * dgd[l] + b;
917
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
918
919
920
  }
  i = height - 1;
  for (j = 1; j < width - 1; ++j) {
David Barker's avatar
David Barker committed
921
    const int k = i * buf_stride + j;
922
923
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
924
    const int nb = 3;
David Barker's avatar
David Barker committed
925
926
927
928
    const int32_t a = A[k] + 2 * (A[k - 1] + A[k + 1]) + A[k - buf_stride] +
                      A[k - buf_stride - 1] + A[k - buf_stride + 1];
    const int32_t b = B[k] + 2 * (B[k - 1] + B[k + 1]) + B[k - buf_stride] +
                      B[k - buf_stride - 1] + B[k - buf_stride + 1];
929
    const int32_t v = a * dgd[l] + b;
930
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
931
932
933
  }
  j = 0;
  for (i = 1; i < height - 1; ++i) {
David Barker's avatar
David Barker committed
934
    const int k = i * buf_stride + j;
935
936
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
937
    const int nb = 3;
David Barker's avatar
David Barker committed
938
939
940
941
    const int32_t a = A[k] + 2 * (A[k - buf_stride] + A[k + buf_stride]) +
                      A[k + 1] + A[k - buf_stride + 1] + A[k + buf_stride + 1];
    const int32_t b = B[k] + 2 * (B[k - buf_stride] + B[k + buf_stride]) +
                      B[k + 1] + B[k - buf_stride + 1] + B[k + buf_stride + 1];
942
    const int32_t v = a * dgd[l] + b;
943
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
944
945
946
  }
  j = width - 1;
  for (i = 1; i < height - 1; ++i) {
David Barker's avatar
David Barker committed
947
    const int k = i * buf_stride + j;
948
949
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
950
    const int nb = 3;
David Barker's avatar
David Barker committed
951
952
953
954
    const int32_t a = A[k] + 2 * (A[k - buf_stride] + A[k + buf_stride]) +
                      A[k - 1] + A[k - buf_stride - 1] + A[k + buf_stride - 1];
    const int32_t b = B[k] + 2 * (B[k - buf_stride] + B[k + buf_stride]) +
                      B[k - 1] + B[k - buf_stride - 1] + B[k + buf_stride - 1];
955
    const int32_t v = a * dgd[l] + b;
956
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
957
958
959
  }
  for (i = 1; i < height - 1; ++i) {
    for (j = 1; j < width - 1; ++j) {
David Barker's avatar
David Barker committed
960
      const int k = i * buf_stride + j;
961
962
      const int l = i * dgd_stride + j;
      const int m = i * dst_stride + j;
963
      const int nb = 5;
964
      const int32_t a =
David Barker's avatar
David Barker committed
965
966
967
968
          (A[k] + A[k - 1] + A[k + 1] + A[k - buf_stride] + A[k + buf_stride]) *
              4 +
          (A[k - 1 - buf_stride] + A[k - 1 + buf_stride] +
           A[k + 1 - buf_stride] + A[k + 1 + buf_stride]) *
969
              3;
970
      const int32_t b =
David Barker's avatar
David Barker committed
971
972
973
974
          (B[k] + B[k - 1] + B[k + 1] + B[k - buf_stride] + B[k + buf_stride]) *
              4 +
          (B[k - 1 - buf_stride] + B[k - 1 + buf_stride] +
           B[k + 1 - buf_stride] + B[k + 1 + buf_stride]) *
975
              3;
976
      const int32_t v = a * dgd[l] + b;
977
      dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
978
979
980
981
    }
  }
}

982
void av1_selfguided_restoration_c(const uint8_t *dgd, int width, int height,
983
                                  int stride, int32_t *dst, int dst_stride,
984
985
                                  int r, int eps) {
  int32_t dgd32_[RESTORATION_PROC_UNIT_PELS];
986
987
  const int dgd32_stride = width + 2 * SGRPROJ_BORDER_HORZ;
  int32_t *dgd32 =
988
      dgd32_ + dgd32_stride * SGRPROJ_BORDER_VERT + SGRPROJ_BORDER_HORZ;
989
  int i, j;
990
991
992
  for (i = -SGRPROJ_BORDER_VERT; i < height + SGRPROJ_BORDER_VERT; ++i) {
    for (j = -SGRPROJ_BORDER_HORZ; j < width + SGRPROJ_BORDER_HORZ; ++j) {
      dgd32[i * dgd32_stride + j] = dgd[i * stride + j];
993
994
    }
  }
995
  av1_selfguided_restoration_internal(dgd32, width, height, dgd32_stride, dst,
996
                                      dst_stride, 8, r, eps);
997
998
}

999
1000
1001
void av1_highpass_filter_c(const uint8_t *dgd, int width, int height,
                           int stride, int32_t *dst, int dst_stride, int corner,
                           int edge) {
1002
  int i, j;
1003
  const int center = (1 << SGRPROJ_RST_BITS) - 4 * (corner + edge);
1004
1005
1006
1007
1008

  i = 0;
  j = 0;
  {
    const int k = i * stride + j;
1009
1010
1011
1012
    const int l = i * dst_stride + j;
    dst[l] =
        center * dgd[k] + edge * (dgd[k + 1] + dgd[k + stride] + dgd[k] * 2) +
        corner * (dgd[k + stride + 1] + dgd[k + 1] + dgd[k + stride] + dgd[k]);
1013
1014
1015
1016
1017
  }
  i = 0;
  j = width - 1;
  {
    const int k = i * stride + j;
1018
1019
1020
1021
    const int l = i * dst_stride + j;
    dst[l] =
        center * dgd[k] + edge * (dgd[k - 1] + dgd[k + stride] + dgd[k] * 2) +
        corner * (dgd[k + stride - 1] + dgd[k - 1] + dgd[k + stride] + dgd[k]);
1022
1023
1024
1025
1026
  }
  i = height - 1;
  j = 0;
  {
    const int k = i * stride + j;
1027
1028
1029
1030
    const int l = i * dst_stride + j;
    dst[l] =
        center * dgd[k] + edge * (dgd[k + 1] + dgd[k - stride] + dgd[k] * 2) +
        corner * (dgd[k - stride + 1] + dgd[k + 1] + dgd[k - stride] + dgd[k]);
1031
1032
1033
1034
1035
  }
  i = height - 1;
  j = width - 1;
  {
    const int k = i * stride + j;
1036
1037
1038
1039
    const int l = i * dst_stride + j;
    dst[l] =
        center * dgd[k] + edge * (dgd[k - 1] + dgd[k - stride] + dgd[k] * 2) +
        corner * (dgd[k - stride - 1] + dgd[k - 1] + dgd[k - stride] + dgd[k]);
1040
1041
1042
1043
  }
  i = 0;
  for (j = 1; j < width - 1; ++j) {
    const int k = i * stride + j;
1044
1045
1046
1047
1048
    const int l = i * dst_stride + j;
    dst[l] = center * dgd[k] +
             edge * (dgd[k - 1] + dgd[k + stride] + dgd[k + 1] + dgd[k]) +
             corner * (dgd[k + stride - 1] + dgd[k + stride + 1] + dgd[k - 1] +
                       dgd[k + 1]);
1049
1050
1051
1052
  }
  i = height - 1;
  for (j = 1; j < width - 1; ++j) {
    const int k = i * stride + j;
1053
1054
1055
1056
1057
    const int l = i * dst_stride + j;
    dst[l] = center * dgd[k] +
             edge * (dgd[k - 1] + dgd[k - stride] + dgd[k + 1] + dgd[k]) +
             corner * (dgd[k - stride - 1] + dgd[k - stride + 1] + dgd[k - 1] +
                       dgd[k + 1]);
1058
1059
1060
1061
  }
  j = 0;
  for (i = 1; i < height - 1; ++i) {
    const int k = i * stride + j;
1062
1063
1064
1065
1066
    const int l = i * dst_stride + j;
    dst[l] = center * dgd[k] +
             edge * (dgd[k - stride] + dgd[k + 1] + dgd[k + stride] + dgd[k]) +
             corner * (dgd[k + stride + 1] + dgd[k - stride + 1] +
                       dgd[k - stride] + dgd[k + stride]);
1067
1068
1069
1070
  }
  j = width - 1;
  for (i = 1; i < height - 1; ++i) {
    const int k = i * stride + j;
1071
1072
1073
1074
1075
    const int l = i * dst_stride + j;
    dst[l] = center * dgd[k] +
             edge * (dgd[k - stride] + dgd[k - 1] + dgd[k + stride] + dgd[k]) +
             corner * (dgd[k + stride - 1] + dgd[k - stride - 1] +
                       dgd[k - stride] + dgd[k + stride]);
1076
1077
1078
1079
  }
  for (i = 1; i < height - 1; ++i) {
    for (j = 1; j < width - 1; ++j) {
      const int k = i * stride + j;
1080
1081
1082
1083
1084
1085
      const int l = i * dst_stride + j;
      dst[l] =
          center * dgd[k] +
          edge * (dgd[k - stride] + dgd[k - 1] + dgd[k + stride] + dgd[k + 1]) +
          corner * (dgd[k + stride - 1] + dgd[k - stride - 1] +
                    dgd[k - stride + 1] + dgd[k + stride + 1]);
1086
1087
1088
1089
    }
  }
}

1090
void apply_selfguided_restoration_c(const uint8_t *dat, int width, int height,
1091
1092
1093
                                    int stride, int eps, const int *xqd,
                                    uint8_t *dst, int dst_stride,
                                    int32_t *tmpbuf) {
1094
  int xq[2];
1095
  int32_t *flt1 = tmpbuf;
1096
  int32_t *flt2 = flt1 + RESTORATION_TILEPELS_MAX;
1097
  int i, j;
1098
  assert(width * height <= RESTORATION_TILEPELS_MAX);
1099
1100
#if USE_HIGHPASS_IN_SGRPROJ
  av1_highpass_filter_c(dat, width, height, stride, flt1, width,
1101
                        sgr_params[eps].corner, sgr_params[eps].edge);
1102
#else
1103
  av1_selfguided_restoration_c(dat, width, height, stride, flt1, width,
1104
                               sgr_params[eps].r1, sgr_params[eps].e1);
1105
#endif  // USE_HIGHPASS_IN_SGRPROJ
1106
  av1_selfguided_restoration_c(dat, width, height, stride, flt2, width,
1107
                               sgr_params[eps].r2, sgr_params[eps].e2);
1108
1109
1110
1111
1112
  decode_xq(xqd, xq);
  for (i = 0; i < height; ++i) {
    for (j = 0; j < width; ++j) {
      const int k = i * width + j;
      const int l = i * stride + j;
1113
1114
1115
1116
      const int m = i * dst_stride + j;
      const int32_t u = ((int32_t)dat[l] << SGRPROJ_RST_BITS);
      const int32_t f1 = (int32_t)flt1[k] - u;
      const int32_t f2 = (int32_t)flt2[k] - u;
David Barker's avatar
David Barker committed
1117
      const int32_t v = xq[0] * f1 + xq[1] * f2 + (u << SGRPROJ_PRJ_BITS);
1118
1119
      const int16_t w =
          (int16_t)ROUND_POWER_OF_TWO(v, SGRPROJ_PRJ_BITS + SGRPROJ_RST_BITS);
1120
      dst[m] = clip_pixel(w);
1121
1122
1123
1124
    }
  }
}

1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
static void sgrproj_filter_stripe(const RestorationUnitInfo *rui,
                                  int stripe_width, int stripe_height,
                                  int procunit_width, const uint8_t *src,
                                  int src_stride, uint8_t *dst, int dst_stride,
                                  int32_t *tmpbuf, int bit_depth) {
  (void)bit_depth;
  assert(bit_depth == 8);

  for (int j = 0; j < stripe_width; j += procunit_width) {
    int w = AOMMIN(procunit_width, stripe_width - j);
    apply_selfguided_restoration(src + j, w, stripe_height, src_stride,
                                 rui->sgrproj_info.ep, rui->sgrproj_info.xqd,
                                 dst + j, dst_stride, tmpbuf);
1138
1139
1140
  }
}

1141
#if CONFIG_HIGHBITDEPTH
1142
#if USE_WIENER_HIGH_INTERMEDIATE_PRECISION
1143
#define wiener_highbd_convolve8_add_src aom_highbd_convolve8_add_src_hip
1144
#else