restoration.c 76.6 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
19
#include "av1/common/onyxc_int.h"
#include "av1/common/restoration.h"
Yaowu Xu's avatar
Yaowu Xu committed
20
21
#include "aom_dsp/aom_dsp_common.h"
#include "aom_mem/aom_mem.h"
22

23
#include "aom_ports/mem.h"
24

25
const sgr_params_type sgr_params[SGRPROJ_PARAMS] = {
26
27
28
29
30
31
32
#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
33
// r1, eps1, r2, eps2
34
#if MAX_RADIUS == 2
35
36
37
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 },
  { 2, 55, 1, 14 }, { 2, 65, 1, 15 }, { 2, 75, 1, 16 }, { 2, 30, 1, 2 },
  { 2, 50, 1, 12 }, { 2, 60, 1, 13 }, { 2, 70, 1, 14 }, { 2, 80, 1, 15 },
#else
40
41
42
43
  { 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 },
44
#endif  // MAX_RADIUS == 2
45
#endif
46
47
};

clang-format's avatar
clang-format committed
48
49
typedef void (*restore_func_type)(uint8_t *data8, int width, int height,
                                  int stride, RestorationInternal *rst,
50
                                  uint8_t *dst8, int dst_stride);
51
#if CONFIG_HIGHBITDEPTH
clang-format's avatar
clang-format committed
52
53
typedef void (*restore_func_highbd_type)(uint8_t *data8, int width, int height,
                                         int stride, RestorationInternal *rst,
54
55
                                         int bit_depth, uint8_t *dst8,
                                         int dst_stride);
56
#endif  // CONFIG_HIGHBITDEPTH
57

58
59
int av1_alloc_restoration_struct(AV1_COMMON *cm, RestorationInfo *rst_info,
                                 int width, int height) {
60
61
  const int ntiles = av1_get_rest_ntiles(
      width, height, rst_info->restoration_tilesize, NULL, NULL, NULL, NULL);
Alex Converse's avatar
Alex Converse committed
62
  aom_free(rst_info->restoration_type);
63
  CHECK_MEM_ERROR(cm, rst_info->restoration_type,
Alex Converse's avatar
Alex Converse committed
64
                  (RestorationType *)aom_malloc(
65
                      sizeof(*rst_info->restoration_type) * ntiles));
66
67
68
69
  aom_free(rst_info->wiener_info);
  CHECK_MEM_ERROR(
      cm, rst_info->wiener_info,
      (WienerInfo *)aom_memalign(16, sizeof(*rst_info->wiener_info) * ntiles));
70
  memset(rst_info->wiener_info, 0, sizeof(*rst_info->wiener_info) * ntiles);
Alex Converse's avatar
Alex Converse committed
71
  aom_free(rst_info->sgrproj_info);
72
73
  CHECK_MEM_ERROR(
      cm, rst_info->sgrproj_info,
Alex Converse's avatar
Alex Converse committed
74
      (SgrprojInfo *)aom_malloc(sizeof(*rst_info->sgrproj_info) * ntiles));
75
76
77
78
79
80
81
82
83
84
  return ntiles;
}

void av1_free_restoration_struct(RestorationInfo *rst_info) {
  aom_free(rst_info->restoration_type);
  rst_info->restoration_type = NULL;
  aom_free(rst_info->wiener_info);
  rst_info->wiener_info = NULL;
  aom_free(rst_info->sgrproj_info);
  rst_info->sgrproj_info = NULL;
85
}
86
87
88

// TODO(debargha): This table can be substantially reduced since only a few
// values are actually used.
David Barker's avatar
David Barker committed
89
int sgrproj_mtable[MAX_EPS][MAX_NELEM];
90
91
92
93
94
95
96
97
98
99

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);
    }
}
100
101

void av1_loop_restoration_precal() { GenSgrprojVtable(); }
102

103
static void loop_restoration_init(RestorationInternal *rst, int kf) {
104
  rst->keyframe = kf;
105
106
}

107
108
void extend_frame(uint8_t *data, int width, int height, int stride,
                  int border_horz, int border_vert) {
109
110
111
112
  uint8_t *data_p;
  int i;
  for (i = 0; i < height; ++i) {
    data_p = data + i * stride;
113
114
    memset(data_p - border_horz, data_p[0], border_horz);
    memset(data_p + width, data_p[width - 1], border_horz);
115
  }
116
117
118
  data_p = data - border_horz;
  for (i = -border_vert; i < 0; ++i) {
    memcpy(data_p + i * stride, data_p, width + 2 * border_horz);
119
  }
120
  for (i = height; i < height + border_vert; ++i) {
121
    memcpy(data_p + i * stride, data_p + (height - 1) * stride,
122
           width + 2 * border_horz);
123
124
125
  }
}

126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
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
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
#if CONFIG_STRIPED_LOOP_RESTORATION

// This function setup a processing stripe by replacing the vertical
// stripe boundary (2 lines above and 2 lines below) by data coming
// from the above/below buffers. Before doing so the original
// frame data is saved into a temporary buffer, such that it
// can be restored by the restore_processing_stripe_boundary
// function after the filtering of the processing stripe.
// Returns the height of the processing stripe
static int setup_processing_stripe_boundary(int y0, int v_end, int h_start,
                                            int h_end, uint8_t *data,
                                            int stride,
                                            RestorationInternal *rst,
                                            int use_highbd) {
  int y, y_stripe_topmost, stripe_index, i;
  int tile_offset = RESTORATION_TILE_OFFSET >> rst->subsampling_y;
  int stripe_height = rst->rsi->procunit_height;
  int comp = rst->component;
  uint8_t *boundary_above_buf = rst->stripe_boundary_above[comp];
  uint8_t *boundary_below_buf = rst->stripe_boundary_below[comp];
  int boundary_stride = rst->stripe_boundary_stride[comp];
  int x0 = h_start - RESTORATION_EXTRA_HORZ;
  int x1 = h_end + RESTORATION_EXTRA_HORZ;

  stripe_index = (y0 + tile_offset) / stripe_height;
  y_stripe_topmost = stripe_index * stripe_height - tile_offset;
  boundary_above_buf +=
      ((stripe_index - 1) * 2 * boundary_stride + RESTORATION_EXTRA_HORZ)
      << use_highbd;
  boundary_below_buf +=
      (stripe_index * 2 * boundary_stride + RESTORATION_EXTRA_HORZ)
      << use_highbd;

  // setup the 2 lines above the stripe
  for (i = 0; i < 2; i++) {
    y = y_stripe_topmost - 2 + i;
    if (y >= 0 && y < y0 && y >= y0 - 2) {
      uint8_t *p = data + ((y * stride + x0) << use_highbd);
      uint8_t *new_data =
          boundary_above_buf + ((i * boundary_stride + x0) << use_highbd);
      // printf("above %3d %3d: %08x %08x : %08x %08x\n", y, x0,
      // ((uint32_t*)p)[0], ((uint32_t*)p)[1], ((uint32_t*)new_data)[0],
      // ((uint32_t*)new_data)[1]);
      // Save old pixels
      memcpy(rst->tmp_save_above[i], p, (x1 - x0) << use_highbd);
      // Replace width pixels from boundary_above_buf
      memcpy(p, new_data, (x1 - x0) << use_highbd);
    }
  }
  // setup the 2 lines below the stripe
  for (i = 0; i < 2; i++) {
    y = y_stripe_topmost + stripe_height + i;
    if (y < v_end + 2) {
      uint8_t *p = data + ((y * stride + x0) << use_highbd);
      uint8_t *new_data =
          boundary_below_buf + ((i * boundary_stride + x0) << use_highbd);
      // printf("below %3d %3d: %08x %08x : %08x %08x\n", y, x0,
      // ((uint32_t*)p)[0], ((uint32_t*)p)[1], ((uint32_t*)new_data)[0],
      // ((uint32_t*)new_data)[1]);
      // Save old pixels
      memcpy(rst->tmp_save_below[i], p, (x1 - x0) << use_highbd);
      // Replace width pixels from boundary_below_buf
      memcpy(p, new_data, (x1 - x0) << use_highbd);
    }
  }
  // Return actual stripe height
  return AOMMIN(v_end, y_stripe_topmost + stripe_height) - y0;
}

// This function restores the boundary lines modified by
// setup_processing_stripe_boundary.
static void restore_processing_stripe_boundary(int y0, int v_end, int h_start,
                                               int h_end, uint8_t *data,
                                               int stride,
                                               RestorationInternal *rst,
                                               int use_highbd) {
  int y, y_stripe_topmost, i, stripe_index;
  int tile_offset = 8 >> rst->subsampling_y;
  int stripe_height = rst->rsi->procunit_height;
  int x0 = h_start - RESTORATION_EXTRA_HORZ;
  int x1 = h_end + RESTORATION_EXTRA_HORZ;

  stripe_index = (y0 + tile_offset) / stripe_height;
  y_stripe_topmost = stripe_index * stripe_height - tile_offset;

  // restore the 2 lines above the stripe
  for (i = 0; i < 2; i++) {
    y = y_stripe_topmost - 2 + i;
    if (y >= 0 && y < y0 && y >= y0 - 2) {
      uint8_t *p = data + ((y * stride + x0) << use_highbd);
      memcpy(p, rst->tmp_save_above[i], (x1 - x0) << use_highbd);
    }
  }
  // restore the 2 lines below the stripe
  for (i = 0; i < 2; i++) {
    y = y_stripe_topmost + stripe_height + i;
    if (y < v_end + 2) {
      uint8_t *p = data + ((y * stride + x0) << use_highbd);
      memcpy(p, rst->tmp_save_below[i], (x1 - x0) << use_highbd);
    }
  }
}

#endif

231
232
static void loop_copy_tile(uint8_t *data, int tile_idx, int width, int height,
                           int stride, RestorationInternal *rst, uint8_t *dst,
233
                           int dst_stride) {
234
235
  const int tile_width = rst->tile_width;
  const int tile_height = rst->tile_height;
236
237
  RestorationTileLimits limits =
      av1_get_rest_tile_limits(tile_idx, rst->nhtiles, rst->nvtiles, tile_width,
238
239
240
#if CONFIG_STRIPED_LOOP_RESTORATION
                               tile_height, width, height, rst->subsampling_y);
#else
241
                               tile_height, width, height);
242
#endif
243
244
245
  for (int i = limits.v_start; i < limits.v_end; ++i)
    memcpy(dst + i * dst_stride + limits.h_start,
           data + i * stride + limits.h_start, limits.h_end - limits.h_start);
246
247
}

248
249
static void stepdown_wiener_kernel(const InterpKernel orig, InterpKernel vert,
                                   int boundary_dist, int istop) {
250
  memcpy(vert, orig, sizeof(InterpKernel));
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
  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;
  }
278
279
}

280
281
static void loop_wiener_filter_tile(uint8_t *data, int tile_idx, int width,
                                    int height, int stride,
282
                                    RestorationInternal *rst, uint8_t *dst,
283
                                    int dst_stride) {
284
  const int procunit_width = rst->rsi->procunit_width;
285
286
287
#if CONFIG_STRIPED_LOOP_RESTORATION
  int procunit_height;
#else
288
  const int procunit_height = rst->rsi->procunit_height;
289
#endif
290
291
  const int tile_width = rst->tile_width;
  const int tile_height = rst->tile_height;
292
  if (rst->rsi->restoration_type[tile_idx] == RESTORE_NONE) {
293
    loop_copy_tile(data, tile_idx, width, height, stride, rst, dst, dst_stride);
294
295
    return;
  }
296
  InterpKernel vertical_topbot;
297
298
  RestorationTileLimits limits =
      av1_get_rest_tile_limits(tile_idx, rst->nhtiles, rst->nvtiles, tile_width,
299
300
301
#if CONFIG_STRIPED_LOOP_RESTORATION
                               tile_height, width, height, rst->subsampling_y);
#else
302
                               tile_height, width, height);
303
304
#endif

305
306
  // Convolve the whole tile (done in blocks here to match the requirements
  // of the vectorized convolve functions, but the result is equivalent)
307
308
309
310
311
312
313
314
315
  for (int i = limits.v_start; i < limits.v_end; i += procunit_height) {
#if CONFIG_STRIPED_LOOP_RESTORATION
    int h = setup_processing_stripe_boundary(
        i, limits.v_end, limits.h_start, limits.h_end, data, stride, rst, 0);
    h = ALIGN_POWER_OF_TWO(h, 1);
    procunit_height = h;
#else
    int h = AOMMIN(procunit_height, (limits.v_end - i + 15) & ~15);
#endif
316
317
    for (int j = limits.h_start; j < limits.h_end; j += procunit_width) {
      int w = AOMMIN(procunit_width, (limits.h_end - j + 15) & ~15);
318
319
      const uint8_t *data_p = data + i * stride + j;
      uint8_t *dst_p = dst + i * dst_stride + j;
320
321
322
323
      // Note h is at least 16
      for (int b = 0; b < WIENER_HALFWIN - WIENER_BORDER_VERT; ++b) {
        stepdown_wiener_kernel(rst->rsi->wiener_info[tile_idx].vfilter,
                               vertical_topbot, WIENER_BORDER_VERT + b, 1);
324
#if USE_WIENER_HIGH_INTERMEDIATE_PRECISION
325
326
327
        aom_convolve8_add_src_hip(data_p, stride, dst_p, dst_stride,
                                  rst->rsi->wiener_info[tile_idx].hfilter, 16,
                                  vertical_topbot, 16, w, 1);
328
#else
329
330
331
        aom_convolve8_add_src(data_p, stride, dst_p, dst_stride,
                              rst->rsi->wiener_info[tile_idx].hfilter, 16,
                              vertical_topbot, 16, w, 1);
332
#endif  // USE_WIENER_HIGH_INTERMEDIATE_PRECISION
333
334
335
        data_p += stride;
        dst_p += dst_stride;
      }
336
337
338
339
#if USE_WIENER_HIGH_INTERMEDIATE_PRECISION
      aom_convolve8_add_src_hip(data_p, stride, dst_p, dst_stride,
                                rst->rsi->wiener_info[tile_idx].hfilter, 16,
                                rst->rsi->wiener_info[tile_idx].vfilter, 16, w,
340
                                h - (WIENER_HALFWIN - WIENER_BORDER_VERT) * 2);
341
342
343
344
#else
      aom_convolve8_add_src(data_p, stride, dst_p, dst_stride,
                            rst->rsi->wiener_info[tile_idx].hfilter, 16,
                            rst->rsi->wiener_info[tile_idx].vfilter, 16, w,
345
                            h - (WIENER_HALFWIN - WIENER_BORDER_VERT) * 2);
346
#endif  // USE_WIENER_HIGH_INTERMEDIATE_PRECISION
347
348
349
350
351
      data_p += stride * (h - (WIENER_HALFWIN - WIENER_BORDER_VERT) * 2);
      dst_p += dst_stride * (h - (WIENER_HALFWIN - WIENER_BORDER_VERT) * 2);
      for (int b = WIENER_HALFWIN - WIENER_BORDER_VERT - 1; b >= 0; --b) {
        stepdown_wiener_kernel(rst->rsi->wiener_info[tile_idx].vfilter,
                               vertical_topbot, WIENER_BORDER_VERT + b, 0);
352
#if USE_WIENER_HIGH_INTERMEDIATE_PRECISION
353
354
355
        aom_convolve8_add_src_hip(data_p, stride, dst_p, dst_stride,
                                  rst->rsi->wiener_info[tile_idx].hfilter, 16,
                                  vertical_topbot, 16, w, 1);
356
#else
357
358
359
        aom_convolve8_add_src(data_p, stride, dst_p, dst_stride,
                              rst->rsi->wiener_info[tile_idx].hfilter, 16,
                              vertical_topbot, 16, w, 1);
360
#endif  // USE_WIENER_HIGH_INTERMEDIATE_PRECISION
361
362
363
        data_p += stride;
        dst_p += dst_stride;
      }
364
    }
365
366
367
368
369
#if CONFIG_STRIPED_LOOP_RESTORATION
    restore_processing_stripe_boundary(i, limits.v_end, limits.h_start,
                                       limits.h_end, data, stride, rst, 0);
#endif
  }
370
371
}

clang-format's avatar
clang-format committed
372
static void loop_wiener_filter(uint8_t *data, int width, int height, int stride,
373
374
375
                               RestorationInternal *rst, uint8_t *dst,
                               int dst_stride) {
  int tile_idx;
376
377
  extend_frame(data, width, height, stride, WIENER_BORDER_HORZ,
               WIENER_BORDER_VERT);
378
  for (tile_idx = 0; tile_idx < rst->ntiles; ++tile_idx) {
379
380
    loop_wiener_filter_tile(data, tile_idx, width, height, stride, rst, dst,
                            dst_stride);
381
  }
382
}
383

384
385
/* 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
386
   specialize to r = 1, 2, 3. A default function is used for r > 3.
387
388
389
390
391
392
393
394
395
396
397
398
399
400

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

402
403
404
405
406
407
408
409
410
411
412
413
414
415
      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;
    }
416
  } else {
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
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
    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;
  }
}

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
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
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);
}

688
689
690
691
692
693
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);
694
695
696
697
  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);
698
699
700
701
}

static void boxnum(int width, int height, int r, int8_t *num, int num_stride) {
  int i, j;
702
703
704
  for (i = 0; i <= r; ++i) {
    for (j = 0; j <= r; ++j) {
      num[i * num_stride + j] = (r + 1 + i) * (r + 1 + j);
705
706
707
708
709
710
      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];
    }
  }
711
712
  for (j = 0; j <= r; ++j) {
    const int val = (2 * r + 1) * (r + 1 + j);
713
714
715
716
717
    for (i = r + 1; i < height - r; ++i) {
      num[i * num_stride + j] = val;
      num[i * num_stride + (width - 1 - j)] = val;
    }
  }
718
719
  for (i = 0; i <= r; ++i) {
    const int val = (2 * r + 1) * (r + 1 + i);
720
721
722
723
724
725
726
    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) {
727
      num[i * num_stride + j] = (2 * r + 1) * (2 * r + 1);
728
729
730
731
732
    }
  }
}

void decode_xq(int *xqd, int *xq) {
733
  xq[0] = xqd[0];
734
735
736
  xq[1] = (1 << SGRPROJ_PRJ_BITS) - xq[0] - xqd[1];
}

David Barker's avatar
David Barker committed
737
const int32_t x_by_xplus1[256] = {
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
  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
758
const int32_t one_by_x[MAX_NELEM] = {
759
  4096, 2048, 1365, 1024, 819, 683, 585, 512, 455, 410, 372, 341, 315,
760
761
762
763
764
  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
765
766
};

767
static void av1_selfguided_restoration_internal(int32_t *dgd, int width,
768
769
                                                int height, int dgd_stride,
                                                int32_t *dst, int dst_stride,
770
                                                int bit_depth, int r, int eps) {
771
772
773
  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
774
775
776
777
  // 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.
778
  int buf_stride = ((width_ext + 3) & ~3) + 16;
779
780
781
782
783
  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];
784
785
  int8_t *num = num_ + SGRPROJ_BORDER_VERT * num_stride + SGRPROJ_BORDER_HORZ;
  int i, j;
786

787
788
789
  // Don't filter tiles with dimensions < 5 on any axis
  if ((width < 5) || (height < 5)) return;

790
791
792
793
794
  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);
795
  assert(r <= 3);
796
797
  A += SGRPROJ_BORDER_VERT * buf_stride + SGRPROJ_BORDER_HORZ;
  B += SGRPROJ_BORDER_VERT * buf_stride + SGRPROJ_BORDER_HORZ;
798
799
  for (i = 0; i < height; ++i) {
    for (j = 0; j < width; ++j) {
David Barker's avatar
David Barker committed
800
      const int k = i * buf_stride + j;
801
      const int n = num[i * num_stride + j];
802

803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
      // 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);
831
832
833
834
835
    }
  }
  i = 0;
  j = 0;
  {
David Barker's avatar
David Barker committed
836
    const int k = i * buf_stride + j;
837
838
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
839
    const int nb = 3;
840
    const int32_t a =
David Barker's avatar
David Barker committed
841
        3 * A[k] + 2 * A[k + 1] + 2 * A[k + buf_stride] + A[k + buf_stride + 1];
842
    const int32_t b =
David Barker's avatar
David Barker committed
843
        3 * B[k] + 2 * B[k + 1] + 2 * B[k + buf_stride] + B[k + buf_stride + 1];
844
    const int32_t v = a * dgd[l] + b;
845
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
846
847
848
849
  }
  i = 0;
  j = width - 1;
  {
David Barker's avatar
David Barker committed
850
    const int k = i * buf_stride + j;
851
852
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
853
    const int nb = 3;
854
    const int32_t a =
David Barker's avatar
David Barker committed
855
        3 * A[k] + 2 * A[k - 1] + 2 * A[k + buf_stride] + A[k + buf_stride - 1];
856
    const int32_t b =
David Barker's avatar
David Barker committed
857
        3 * B[k] + 2 * B[k - 1] + 2 * B[k + buf_stride] + B[k + buf_stride - 1];
858
    const int32_t v = a * dgd[l] + b;
859
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
860
861
862
863
  }
  i = height - 1;
  j = 0;
  {
David Barker's avatar
David Barker committed
864
    const int k = i * buf_stride + j;
865
866
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
867
    const int nb = 3;
868
    const int32_t a =
David Barker's avatar
David Barker committed
869
        3 * A[k] + 2 * A[k + 1] + 2 * A[k - buf_stride] + A[k - buf_stride + 1];
870
    const int32_t b =
David Barker's avatar
David Barker committed
871
        3 * B[k] + 2 * B[k + 1] + 2 * B[k - buf_stride] + B[k - buf_stride + 1];
872
    const int32_t v = a * dgd[l] + b;
873
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
874
875
876
877
  }
  i = height - 1;
  j = width - 1;
  {
David Barker's avatar
David Barker committed
878
    const int k = i * buf_stride + j;
879
880
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
881
    const int nb = 3;
882
    const int32_t a =
David Barker's avatar
David Barker committed
883
        3 * A[k] + 2 * A[k - 1] + 2 * A[k - buf_stride] + A[k - buf_stride - 1];
884
    const int32_t b =
David Barker's avatar
David Barker committed
885
        3 * B[k] + 2 * B[k - 1] + 2 * B[k - buf_stride] + B[k - buf_stride - 1];
886
    const int32_t v = a * dgd[l] + b;
887
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
888
889
890
  }
  i = 0;
  for (j = 1; j < width - 1; ++j) {
David Barker's avatar
David Barker committed
891
    const int k = i * buf_stride + j;
892
893
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
894
    const int nb = 3;
David Barker's avatar
David Barker committed
895
896
897
898
    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];
899
    const int32_t v = a * dgd[l] + b;
900
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
901
902
903
  }
  i = height - 1;
  for (j = 1; j < width - 1; ++j) {
David Barker's avatar
David Barker committed
904
    const int k = i * buf_stride + j;
905
906
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
907
    const int nb = 3;
David Barker's avatar
David Barker committed
908
909
910
911
    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];
912
    const int32_t v = a * dgd[l] + b;
913
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
914
915
916
  }
  j = 0;
  for (i = 1; i < height - 1; ++i) {
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;
David Barker's avatar
David Barker committed
921
922
923
924
    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];
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
  }
  j = width - 1;
  for (i = 1; i < height - 1; ++i) {
David Barker's avatar
David Barker committed
930
    const int k = i * buf_stride + j;
931
932
    const int l = i * dgd_stride + j;
    const int m = i * dst_stride + j;
933
    const int nb = 3;
David Barker's avatar
David Barker committed
934
935
936
937
    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];
938
    const int32_t v = a * dgd[l] + b;
939
    dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
940
941
942
  }
  for (i = 1; i < height - 1; ++i) {
    for (j = 1; j < width - 1; ++j) {
David Barker's avatar
David Barker committed
943
      const int k = i * buf_stride + j;
944
945
      const int l = i * dgd_stride + j;
      const int m = i * dst_stride + j;
946
      const int nb = 5;
947
      const int32_t a =
David Barker's avatar
David Barker committed
948
949
950
951
          (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]) *
952
              3;
953
      const int32_t b =
David Barker's avatar
David Barker committed
954
955
956
957
          (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]) *
958
              3;
959
      const int32_t v = a * dgd[l] + b;
960
      dst[m] = ROUND_POWER_OF_TWO(v, SGRPROJ_SGR_BITS + nb - SGRPROJ_RST_BITS);
961
962
963
964
    }
  }
}

965
966
void av1_selfguided_restoration_c(uint8_t *dgd, int width, int height,
                                  int stride, int32_t *dst, int dst_stride,
967
968
                                  int r, int eps) {
  int32_t dgd32_[RESTORATION_PROC_UNIT_PELS];
969
970
  const int dgd32_stride = width + 2 * SGRPROJ_BORDER_HORZ;
  int32_t *dgd32 =
971
      dgd32_ + dgd32_stride * SGRPROJ_BORDER_VERT + SGRPROJ_BORDER_HORZ;
972
  int i, j;
973
974
975
  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];
976
977
    }
  }
978
  av1_selfguided_restoration_internal(dgd32, width, height, dgd32_stride, dst,
979
                                      dst_stride, 8, r, eps);
980
981
}

982
983
void av1_highpass_filter_c(uint8_t *dgd, int width, int height, int stride,
                           int32_t *dst, int dst_stride, int corner, int edge) {
984
  int i, j;
985
  const int center = (1 << SGRPROJ_RST_BITS) - 4 * (corner + edge);
986
987
988
989
990

  i = 0;
  j = 0;
  {
    const int k = i * stride + j;
991
992
993
994
    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]);
995
996
997
998
999
  }
  i = 0;
  j = width - 1;
  {
    const int k = i * stride + j;
1000
1001
1002
1003
    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]);
1004
1005
1006
1007
1008
  }
  i = height - 1;
  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 = height - 1;
  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
  }
  i = 0;
  for (j = 1; j < width - 1; ++j) {
    const int k = i * stride + j;
1026
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 + 1] + dgd[k]) +
             corner * (dgd[k + stride - 1] + dgd[k + stride + 1] + dgd[k - 1] +
                       dgd[k + 1]);
1031
1032
1033
1034
  }
  i = height - 1;
  for (j = 1; j < width - 1; ++j) {
    const int k = i * stride + j;
1035
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 + 1] + dgd[k]) +
             corner * (dgd[k - stride - 1] + dgd[k - stride + 1] + dgd[k - 1] +
                       dgd[k + 1]);
1040
1041
1042
1043
  }
  j = 0;
  for (i = 1; i < height - 1; ++i) {
    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 - stride] + dgd[k + 1] + dgd[k + stride] + dgd[k]) +
             corner * (dgd[k + stride + 1] + dgd[k - stride + 1] +
                       dgd[k - stride] + dgd[k + stride]);
1049
1050
1051
1052
  }
  j = width - 1;
  for (i = 1; i < height - 1; ++i) {
    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 - stride] + dgd[k - 1] + dgd[k + stride] + dgd[k]) +
             corner * (dgd[k + stride - 1] + dgd[k - stride - 1] +
                       dgd[k - stride] + dgd[k + stride]);
1058
1059
1060
1061
  }
  for (i = 1; i < height - 1; ++i) {
    for (j = 1; j < width - 1; ++j) {
      const int k = i * stride + j;
1062
1063
1064
1065
1066
1067
      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]);
1068
1069
1070
1071
    }
  }
}

David Barker's avatar
David Barker committed
1072
void apply_selfguided_restoration_c(uint8_t *dat, int width, int height,
1073
1074
                                    int stride, int eps, int *xqd, uint8_t *dst,
                                    int dst_stride, int32_t *tmpbuf) {
1075
  int xq[2];
1076
  int32_t *flt1 = tmpbuf;
1077
  int32_t *flt2 = flt1 + RESTORATION_TILEPELS_MAX;
1078
  int i, j;
1079
  assert(width * height <= RESTORATION_TILEPELS_MAX);
1080
1081
#if USE_HIGHPASS_IN_SGRPROJ
  av1_highpass_filter_c(dat, width, height, stride, flt1, width,
1082
                        sgr_params[eps].corner, sgr_params[eps].edge);
1083
#else
1084
  av1_selfguided_restoration_c(dat, width, height, stride, flt1, width,
1085
                               sgr_params[eps].r1, sgr_params[eps].e1);
1086
#endif  // USE_HIGHPASS_IN_SGRPROJ
1087
  av1_selfguided_restoration_c(dat, width, height, stride, flt2, width,
1088
                               sgr_params[eps].r2, sgr_params[eps].e2);
1089
1090
1091
1092
1093
  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;
1094
1095
1096
1097
      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
1098
      const int32_t v = xq[0] * f1 + xq[1] * f2 + (u << SGRPROJ_PRJ_BITS);
1099
1100
      const int16_t w =
          (int16_t)ROUND_POWER_OF_TWO(v, SGRPROJ_PRJ_BITS + SGRPROJ_RST_BITS);
1101
      dst[m] = clip_pixel(w);
1102
1103
1104
1105
1106
1107
    }
  }
}

static void loop_sgrproj_filter_tile(uint8_t *data, int tile_idx, int width,
                                     int height, int stride,
1108
1109
                                     RestorationInternal *rst, uint8_t *dst,
                                     int dst_stride) {
1110
  const int procunit_width = rst->rsi->procunit_width;
1111
1112
1113
#if CONFIG_STRIPED_LOOP_RESTORATION
  int procunit_height;
#else
1114
  const int procunit_height = rst->rsi->procunit_height;
1115
#endif
1116
1117
  const int tile_width = rst->tile_width;
  const int tile_height = rst->tile_height;