restoration.c 72.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
257
258
259
260
261
262
263
264
// restoration unit. rsb is the stored stripe boundaries (the saved output from
// the deblocker). stripe_height is the height of each stripe. 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)
static int setup_processing_stripe_boundary(
    const RestorationTileLimits *limits, const RestorationStripeBoundaries *rsb,
    int stripe_height, int ss_y, int use_highbd, uint8_t *data8, int stride,
    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
270
271
272
273
274
  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.
275
  const int buf_stride = rsb->stripe_boundary_stride;
276
277
278
279
280
  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;
281

282
283
  assert(CONFIG_HIGHBITDEPTH || !use_highbd);

284
285
  // Replace the pixels above the top of the stripe, unless this is the top of
  // the image.
286
287
  // 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.
288
289
  if (stripe_index > 0) {
    const int above_buf_y = 2 * (stripe_index - 1);
290
    uint8_t *data8_tl = data8 + (limits->v_start - 3) * stride + data_x0_off;
291

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

  // Replace the pixels below the bottom of the stripe if necessary. This might
  // not be needed if the stripe is less than stripe_height high (which might
  // happen on the bottom of a loop restoration unit), in which case
  // rows_needed_below might be negative.
307
308
  // Similarly to above, we expand 2 lines from rb->stripe_boundary_below into
  // 3 lines of below pixels. This time we duplicate the bottommost row.
309
  const int stripe_bottom = stripe_height * (1 + stripe_index) - tile_offset;
310
  const int rows_needed_below = AOMMIN(limits->v_end + 3 - stripe_bottom, 3);
311
312
313
314
315

  const int below_buf_y = 2 * stripe_index;
  uint8_t *data8_bl = data8 + stripe_bottom * stride + data_x0_off;

  for (int i = 0; i < rows_needed_below; ++i) {
316
317
    const int src_row = AOMMIN(1, i);
    const int buf_off = buf_x0_off + (below_buf_y + src_row) * buf_stride;
318
    const uint8_t *src = rsb->stripe_boundary_below + (buf_off << use_highbd);
319
320
    uint8_t *dst8 = data8_bl + i * stride;
    // Save old pixels, then replace with data from boundary_below_buf
321
    memcpy(rlbs->tmp_save_below[i], REAL_PTR(use_highbd, dst8), line_size);
322
    memcpy(REAL_PTR(use_highbd, dst8), src, line_size);
323
  }
324

325
326
  // Finally, return the actual height of this stripe.
  return AOMMIN(limits->v_end, stripe_bottom) - limits->v_start;
327
328
329
330
}

// This function restores the boundary lines modified by
// setup_processing_stripe_boundary.
331
static void restore_processing_stripe_boundary(
332
333
334
    const RestorationTileLimits *limits, const RestorationLineBuffers *rlbs,
    int stripe_height, int ss_y, int use_highbd, uint8_t *data8, int stride) {
  const int tile_offset = RESTORATION_TILE_OFFSET >> ss_y;
335
336
337
338
339
340
  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;
341

342
343
  assert(CONFIG_HIGHBITDEPTH || !use_highbd);

344
  if (stripe_index > 0) {
345
346
    uint8_t *data8_tl = data8 + (limits->v_start - 3) * stride + data_x0_off;
    for (int i = 0; i < 3; ++i) {
347
348
      uint8_t *dst8 = data8_tl + i * stride;
      // Save old pixels, then replace with data from boundary_above_buf
349
      memcpy(REAL_PTR(use_highbd, dst8), rlbs->tmp_save_above[i], line_size);
350
351
    }
  }
352
353

  const int stripe_bottom = stripe_height * (1 + stripe_index) - tile_offset;
354
  const int rows_needed_below = AOMMIN(limits->v_end + 3 - stripe_bottom, 3);
355
356
357
358
359
360

  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
361
    memcpy(REAL_PTR(use_highbd, dst8), rlbs->tmp_save_below[i], line_size);
362
363
  }
}
364
#undef REAL_PTR
365
366
#endif

367
368
static void stepdown_wiener_kernel(const InterpKernel orig, InterpKernel vert,
                                   int boundary_dist, int istop) {
369
  memcpy(vert, orig, sizeof(InterpKernel));
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
  switch (boundary_dist) {
    case 0:
      vert[WIENER_HALFWIN] += vert[2] + vert[1] + vert[0];
      vert[2] = vert[1] = vert[0] = 0;
      break;
    case 1:
      vert[2] += vert[1] + vert[0];
      vert[1] = vert[0] = 0;
      break;
    case 2:
      vert[1] += vert[0];
      vert[0] = 0;
      break;
    default: break;
  }
  if (!istop) {
    int tmp;
    tmp = vert[0];
    vert[0] = vert[WIENER_WIN - 1];
    vert[WIENER_WIN - 1] = tmp;
    tmp = vert[1];
    vert[1] = vert[WIENER_WIN - 2];
    vert[WIENER_WIN - 2] = tmp;
    tmp = vert[2];
    vert[2] = vert[WIENER_WIN - 3];
    vert[WIENER_WIN - 3] = tmp;
  }
397
398
}

399
400
#if USE_WIENER_HIGH_INTERMEDIATE_PRECISION
#define wiener_convolve8_add_src aom_convolve8_add_src_hip
401
#else
402
#define wiener_convolve8_add_src aom_convolve8_add_src
403
404
#endif

405
406
407
408
409
410
411
412
413
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);

414
415
416
  const int mid_height =
      stripe_height - (WIENER_HALFWIN - WIENER_BORDER_VERT) * 2;
  assert(mid_height > 0);
417
418
  for (int j = 0; j < stripe_width; j += procunit_width) {
    int w = AOMMIN(procunit_width, (stripe_width - j + 15) & ~15);
419
420
421
422
    const uint8_t *src_p = src + j;
    uint8_t *dst_p = dst + j;
    for (int b = 0; b < WIENER_HALFWIN - WIENER_BORDER_VERT; ++b) {
      InterpKernel vertical_top;
423
      stepdown_wiener_kernel(rui->wiener_info.vfilter, vertical_top,
424
425
                             WIENER_BORDER_VERT + b, 1);
      wiener_convolve8_add_src(src_p, src_stride, dst_p, dst_stride,
426
427
                               rui->wiener_info.hfilter, 16, vertical_top, 16,
                               w, 1);
428
429
      src_p += src_stride;
      dst_p += dst_stride;
430
431
    }

432
    wiener_convolve8_add_src(src_p, src_stride, dst_p, dst_stride,
433
434
                             rui->wiener_info.hfilter, 16,
                             rui->wiener_info.vfilter, 16, w, mid_height);
435
436
437
438
439
    src_p += src_stride * mid_height;
    dst_p += dst_stride * mid_height;

    for (int b = WIENER_HALFWIN - WIENER_BORDER_VERT - 1; b >= 0; --b) {
      InterpKernel vertical_bot;
440
      stepdown_wiener_kernel(rui->wiener_info.vfilter, vertical_bot,
441
442
                             WIENER_BORDER_VERT + b, 0);
      wiener_convolve8_add_src(src_p, src_stride, dst_p, dst_stride,
443
444
                               rui->wiener_info.hfilter, 16, vertical_bot, 16,
                               w, 1);
445
446
447
      src_p += src_stride;
      dst_p += dst_stride;
    }
448
  }
449
}
450

451
452
/* 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
453
   specialize to r = 1, 2, 3. A default function is used for r > 3.
454
455
456
457
458
459
460
461
462
463
464
465
466
467

   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];
468

469
470
471
472
473
474
475
476
477
478
479
480
481
482
      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;
    }
483
  } else {
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
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
    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;
  }
}

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
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
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);
}

755
756
757
758
759
760
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);
761
762
763
764
  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);
765
766
767
768
}

static void boxnum(int width, int height, int r, int8_t *num, int num_stride) {
  int i, j;
769
770
771
  for (i = 0; i <= r; ++i) {
    for (j = 0; j <= r; ++j) {
      num[i * num_stride + j] = (r + 1 + i) * (r + 1 + j);
772
773
774
775
776
777
      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];
    }
  }
778
779
  for (j = 0; j <= r; ++j) {
    const int val = (2 * r + 1) * (r + 1 + j);
780
781
782
783
784
    for (i = r + 1; i < height - r; ++i) {
      num[i * num_stride + j] = val;
      num[i * num_stride + (width - 1 - j)] = val;
    }
  }
785
786
  for (i = 0; i <= r; ++i) {
    const int val = (2 * r + 1) * (r + 1 + i);
787
788
789
790
791
792
793
    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) {
794
      num[i * num_stride + j] = (2 * r + 1) * (2 * r + 1);
795
796
797
798
    }
  }
}

799
void decode_xq(const int *xqd, int *xq) {
800
  xq[0] = xqd[0];
801
802
803
  xq[1] = (1 << SGRPROJ_PRJ_BITS) - xq[0] - xqd[1];
}

David Barker's avatar
David Barker committed
804
const int32_t x_by_xplus1[256] = {
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
  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
825
const int32_t one_by_x[MAX_NELEM] = {
826
  4096, 2048, 1365, 1024, 819, 683, 585, 512, 455, 410, 372, 341, 315,
827
828
829
830
831
  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
832
833
};

834
static void av1_selfguided_restoration_internal(int32_t *dgd, int width,
835
836
                                                int height, int dgd_stride,
                                                int32_t *dst, int dst_stride,
837
                                                int bit_depth, int r, int eps) {
838
839
840
  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
841
842
843
844
  // 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.
845
  int buf_stride = ((width_ext + 3) & ~3) + 16;
846
847
848
849
850
  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];
851
852
  int8_t *num = num_ + SGRPROJ_BORDER_VERT * num_stride + SGRPROJ_BORDER_HORZ;
  int i, j;
853

854
855
856
  // Don't filter tiles with dimensions < 5 on any axis
  if ((width < 5) || (height < 5)) return;

857
858
859
860
861
  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);
862
  assert(r <= 3);
863
864
  A += SGRPROJ_BORDER_VERT * buf_stride + SGRPROJ_BORDER_HORZ;
  B += SGRPROJ_BORDER_VERT * buf_stride + SGRPROJ_BORDER_HORZ;
865
866
  for (i = 0; i < height; ++i) {
    for (j = 0; j < width; ++j) {
David Barker's avatar
David Barker committed
867
      const int k = i * buf_stride + j;
868
      const int n = num[i * num_stride + j];
869

870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
      // 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);
898
899
900
901
902
    }
  }
  i = 0;
  j = 0;
  {
David Barker's avatar
David Barker committed
903
    const int k = i * buf_stride + j;
904
905
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
906
    const int nb = 3;
907
    const int32_t a =
David Barker's avatar
David Barker committed
908
        3 * A[k] + 2 * A[k + 1] + 2 * A[k + buf_stride] + A[k + buf_stride + 1];
909
    const int32_t b =
David Barker's avatar
David Barker committed
910
        3 * B[k] + 2 * B[k + 1] + 2 * B[k + buf_stride] + B[k + buf_stride + 1];
911
    const int32_t v = a * dgd[l] + b;
912
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
913
914
915
916
  }
  i = 0;
  j = width - 1;
  {
David Barker's avatar
David Barker committed
917
    const int k = i * buf_stride + j;
918
919
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
920
    const int nb = 3;
921
    const int32_t a =
David Barker's avatar
David Barker committed
922
        3 * A[k] + 2 * A[k - 1] + 2 * A[k + buf_stride] + A[k + buf_stride - 1];
923
    const int32_t b =
David Barker's avatar
David Barker committed
924
        3 * B[k] + 2 * B[k - 1] + 2 * B[k + buf_stride] + B[k + buf_stride - 1];
925
    const int32_t v = a * dgd[l] + b;
926
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
927
928
929
930
  }
  i = height - 1;
  j = 0;
  {
David Barker's avatar
David Barker committed
931
    const int k = i * buf_stride + j;
932
933
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
934
    const int nb = 3;
935
    const int32_t a =
David Barker's avatar
David Barker committed
936
        3 * A[k] + 2 * A[k + 1] + 2 * A[k - buf_stride] + A[k - buf_stride + 1];
937
    const int32_t b =
David Barker's avatar
David Barker committed
938
        3 * B[k] + 2 * B[k + 1] + 2 * B[k - buf_stride] + B[k - buf_stride + 1];
939
    const int32_t v = a * dgd[l] + b;
940
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
941
942
943
944
  }
  i = height - 1;
  j = width - 1;
  {
David Barker's avatar
David Barker committed
945
    const int k = i * buf_stride + j;
946
947
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
948
    const int nb = 3;
949
    const int32_t a =
David Barker's avatar
David Barker committed
950
        3 * A[k] + 2 * A[k - 1] + 2 * A[k - buf_stride] + A[k - buf_stride - 1];
951
    const int32_t b =
David Barker's avatar
David Barker committed
952
        3 * B[k] + 2 * B[k - 1] + 2 * B[k - buf_stride] + B[k - buf_stride - 1];
953
    const int32_t v = a * dgd[l] + b;
954
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
955
956
957
  }
  i = 0;
  for (j = 1; j < width - 1; ++j) {
David Barker's avatar
David Barker committed
958
    const int k = i * buf_stride + j;
959
960
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
961
    const int nb = 3;
David Barker's avatar
David Barker committed
962
963
964
965
    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];
966
    const int32_t v = a * dgd[l] + b;
967
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
968
969
970
  }
  i = height - 1;
  for (j = 1; j < width - 1; ++j) {
David Barker's avatar
David Barker committed
971
    const int k = i * buf_stride + j;
972
973
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
974
    const int nb = 3;
David Barker's avatar
David Barker committed
975
976
977
978
    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];
979
    const int32_t v = a * dgd[l] + b;
980
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
981
982
983
  }
  j = 0;
  for (i = 1; i < height - 1; ++i) {
David Barker's avatar
David Barker committed
984
    const int k = i * buf_stride + j;
985
986
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
987
    const int nb = 3;
David Barker's avatar
David Barker committed
988
989
990
991
    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];
992
    const int32_t v = a * dgd[l] + b;
993
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
994
995
996
  }
  j = width - 1;
  for (i = 1; i < height - 1; ++i) {
David Barker's avatar
David Barker committed
997
    const int k = i * buf_stride + j;
998
999
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
1000
    const int nb = 3;
David Barker's avatar
David Barker committed
1001
1002
1003
1004
    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];
1005
    const int32_t v = a * dgd[l] + b;
1006
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
1007
1008
1009
  }
  for (i = 1; i < height - 1; ++i) {
    for (j = 1; j < width - 1; ++j) {
David Barker's avatar
David Barker committed
1010
      const int k = i * buf_stride + j;
1011
1012
      const int l = i * dgd_stride + j;
      const int m = i * dst_stride + j;
1013
      const int nb = 5;
1014
      const int32_t a =
David Barker's avatar
David Barker committed
1015
1016
1017
1018
          (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]) *
1019
              3;
1020
      const int32_t b =
David Barker's avatar
David Barker committed
1021
1022
1023
1024
          (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]) *
1025
              3;
1026
      const int32_t v = a * dgd[l] + b;
1027
      dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
1028
1029
1030
1031
    }
  }
}

1032
void av1_selfguided_restoration_c(const uint8_t *dgd, int width, int height,
1033
                                  int stride, int32_t *dst, int dst_stride,
1034
1035
                                  int r, int eps) {
  int32_t dgd32_[RESTORATION_PROC_UNIT_PELS];
1036
1037
  const int dgd32_stride = width + 2 * SGRPROJ_BORDER_HORZ;
  int32_t *dgd32 =
1038
      dgd32_ + dgd32_stride * SGRPROJ_BORDER_VERT + SGRPROJ_BORDER_HORZ;
1039
  int i, j;
1040
1041
1042
  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];
1043
1044
    }
  }
1045
  av1_selfguided_restoration_internal(dgd32, width, height, dgd32_stride, dst,
1046
                                      dst_stride, 8, r, eps);
1047
1048
}

1049
1050
1051
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) {
1052
  int i, j;
1053
  const int center = (1 << SGRPROJ_RST_BITS) - 4 * (corner + edge);
1054
1055
1056
1057
1058

  i = 0;
  j = 0;
  {
    const int k = i * stride + j;
1059
1060
1061
1062
    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]);
1063
1064
1065
1066
1067
  }
  i = 0;
  j = width - 1;
  {
    const int k = i * stride + j;
1068
1069
1070
1071
    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]);
1072
1073
1074
1075
1076
  }
  i = height - 1;
  j = 0;
  {
    const int k = i * stride + j;
1077
1078
1079
1080
    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]);
1081
1082
1083
1084
1085
  }
  i = height - 1;
  j = width - 1;
  {
    const int k = i * stride + j;
1086
1087
1088
1089
    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]);
1090
1091
1092
1093
  }
  i = 0;
  for (j = 1; j < width - 1; ++j) {
    const int k = i * stride + j;
1094
1095
1096
1097
1098
    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]);
1099
1100
1101
1102
  }
  i = height - 1;
  for (j = 1; j < width - 1; ++j) {
    const int k = i * stride + j;
1103
1104
1105
1106
1107
    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]);
1108
1109
1110
1111
  }
  j = 0;
  for (i = 1; i < height - 1; ++i) {
    const int k = i * stride + j;
1112
1113
1114
1115
1116
    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]);
1117
1118
1119
1120
  }
  j = width - 1;
  for (i = 1; i < height - 1; ++i) {
    const int k = i * stride + j;
1121
1122
1123
1124
1125
    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]);
1126
1127
1128
1129
  }
  for (i = 1; i < height - 1; ++i) {
    for (j = 1; j < width - 1; ++j) {
      const int k = i * stride + j;
1130
1131
1132
1133
1134
1135
      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]);
1136
1137
1138
1139
    }
  }
}

1140
void apply_selfguided_restoration_c(const uint8_t *dat, int width, int height,
1141
1142
1143
                                    int stride, int eps, const int *xqd,
                                    uint8_t *dst, int dst_stride,
                                    int32_t *tmpbuf) {
1144
  int xq[2];
1145
  int32_t *flt1 = tmpbuf;
1146
  int32_t *flt2 = flt1 + RESTORATION_TILEPELS_MAX;
1147
  int i, j;