Coverage Report

Created: 2022-08-24 06:15

/src/aom/av1/encoder/pickrst.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright (c) 2016, Alliance for Open Media. All rights reserved
3
 *
4
 * This source code is subject to the terms of the BSD 2 Clause License and
5
 * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
6
 * was not distributed with this source code in the LICENSE file, you can
7
 * obtain it at www.aomedia.org/license/software. If the Alliance for Open
8
 * Media Patent License 1.0 was not distributed with this source code in the
9
 * PATENTS file, you can obtain it at www.aomedia.org/license/patent.
10
 */
11
12
#include <assert.h>
13
#include <float.h>
14
#include <limits.h>
15
#include <math.h>
16
17
#include "config/aom_scale_rtcd.h"
18
#include "config/av1_rtcd.h"
19
20
#include "aom_dsp/aom_dsp_common.h"
21
#include "aom_dsp/binary_codes_writer.h"
22
#include "aom_dsp/mathutils.h"
23
#include "aom_dsp/psnr.h"
24
#include "aom_mem/aom_mem.h"
25
#include "aom_ports/mem.h"
26
#include "av1/common/av1_common_int.h"
27
#include "av1/common/quant_common.h"
28
#include "av1/common/restoration.h"
29
30
#include "av1/encoder/av1_quantize.h"
31
#include "av1/encoder/encoder.h"
32
#include "av1/encoder/picklpf.h"
33
#include "av1/encoder/pickrst.h"
34
35
// When set to RESTORE_WIENER or RESTORE_SGRPROJ only those are allowed.
36
// When set to RESTORE_TYPES we allow switchable.
37
static const RestorationType force_restore_type = RESTORE_TYPES;
38
39
// Number of Wiener iterations
40
0
#define NUM_WIENER_ITERS 5
41
42
// Penalty factor for use of dual sgr
43
0
#define DUAL_SGR_PENALTY_MULT 0.01
44
45
// Working precision for Wiener filter coefficients
46
0
#define WIENER_TAP_SCALE_FACTOR ((int64_t)1 << 16)
47
48
0
#define SGRPROJ_EP_GRP1_START_IDX 0
49
0
#define SGRPROJ_EP_GRP1_END_IDX 9
50
0
#define SGRPROJ_EP_GRP1_SEARCH_COUNT 4
51
0
#define SGRPROJ_EP_GRP2_3_SEARCH_COUNT 2
52
static const int sgproj_ep_grp1_seed[SGRPROJ_EP_GRP1_SEARCH_COUNT] = { 0, 3, 6,
53
                                                                       9 };
54
static const int sgproj_ep_grp2_3[SGRPROJ_EP_GRP2_3_SEARCH_COUNT][14] = {
55
  { 10, 10, 11, 11, 12, 12, 13, 13, 13, 13, -1, -1, -1, -1 },
56
  { 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 15 }
57
};
58
59
typedef int64_t (*sse_extractor_type)(const YV12_BUFFER_CONFIG *a,
60
                                      const YV12_BUFFER_CONFIG *b);
61
typedef int64_t (*sse_part_extractor_type)(const YV12_BUFFER_CONFIG *a,
62
                                           const YV12_BUFFER_CONFIG *b,
63
                                           int hstart, int width, int vstart,
64
                                           int height);
65
typedef uint64_t (*var_part_extractor_type)(const YV12_BUFFER_CONFIG *a,
66
                                            int hstart, int width, int vstart,
67
                                            int height);
68
69
#if CONFIG_AV1_HIGHBITDEPTH
70
#define NUM_EXTRACTORS (3 * (1 + 1))
71
#else
72
#define NUM_EXTRACTORS 3
73
#endif
74
static const sse_part_extractor_type sse_part_extractors[NUM_EXTRACTORS] = {
75
  aom_get_y_sse_part,        aom_get_u_sse_part,
76
  aom_get_v_sse_part,
77
#if CONFIG_AV1_HIGHBITDEPTH
78
  aom_highbd_get_y_sse_part, aom_highbd_get_u_sse_part,
79
  aom_highbd_get_v_sse_part,
80
#endif
81
};
82
static const var_part_extractor_type var_part_extractors[NUM_EXTRACTORS] = {
83
  aom_get_y_var,        aom_get_u_var,        aom_get_v_var,
84
#if CONFIG_AV1_HIGHBITDEPTH
85
  aom_highbd_get_y_var, aom_highbd_get_u_var, aom_highbd_get_v_var,
86
#endif
87
};
88
89
static int64_t sse_restoration_unit(const RestorationTileLimits *limits,
90
                                    const YV12_BUFFER_CONFIG *src,
91
                                    const YV12_BUFFER_CONFIG *dst, int plane,
92
0
                                    int highbd) {
93
0
  return sse_part_extractors[3 * highbd + plane](
94
0
      src, dst, limits->h_start, limits->h_end - limits->h_start,
95
0
      limits->v_start, limits->v_end - limits->v_start);
96
0
}
97
98
static uint64_t var_restoration_unit(const RestorationTileLimits *limits,
99
                                     const YV12_BUFFER_CONFIG *src, int plane,
100
0
                                     int highbd) {
101
0
  return var_part_extractors[3 * highbd + plane](
102
0
      src, limits->h_start, limits->h_end - limits->h_start, limits->v_start,
103
0
      limits->v_end - limits->v_start);
104
0
}
105
106
typedef struct {
107
  // The best coefficients for Wiener or Sgrproj restoration
108
  WienerInfo wiener;
109
  SgrprojInfo sgrproj;
110
111
  // The sum of squared errors for this rtype.
112
  int64_t sse[RESTORE_SWITCHABLE_TYPES];
113
114
  // The rtype to use for this unit given a frame rtype as
115
  // index. Indices: WIENER, SGRPROJ, SWITCHABLE.
116
  RestorationType best_rtype[RESTORE_TYPES - 1];
117
118
  // This flag will be set based on the speed feature
119
  // 'prune_sgr_based_on_wiener'. 0 implies no pruning and 1 implies pruning.
120
  uint8_t skip_sgr_eval;
121
} RestUnitSearchInfo;
122
123
typedef struct {
124
  const YV12_BUFFER_CONFIG *src;
125
  YV12_BUFFER_CONFIG *dst;
126
127
  const AV1_COMMON *cm;
128
  const MACROBLOCK *x;
129
  int plane;
130
  int plane_width;
131
  int plane_height;
132
  RestUnitSearchInfo *rusi;
133
134
  // Speed features
135
  const LOOP_FILTER_SPEED_FEATURES *lpf_sf;
136
137
  uint8_t *dgd_buffer;
138
  int dgd_stride;
139
  const uint8_t *src_buffer;
140
  int src_stride;
141
142
  // sse and bits are initialised by reset_rsc in search_rest_type
143
  int64_t sse;
144
  int64_t bits;
145
  int tile_y0, tile_stripe0;
146
147
  // sgrproj and wiener are initialised by rsc_on_tile when starting the first
148
  // tile in the frame.
149
  SgrprojInfo sgrproj;
150
  WienerInfo wiener;
151
  AV1PixelRect tile_rect;
152
} RestSearchCtxt;
153
154
0
static AOM_INLINE void rsc_on_tile(void *priv) {
155
0
  RestSearchCtxt *rsc = (RestSearchCtxt *)priv;
156
0
  set_default_sgrproj(&rsc->sgrproj);
157
0
  set_default_wiener(&rsc->wiener);
158
0
  rsc->tile_stripe0 = 0;
159
0
}
160
161
0
static AOM_INLINE void reset_rsc(RestSearchCtxt *rsc) {
162
0
  rsc->sse = 0;
163
0
  rsc->bits = 0;
164
0
}
165
166
static AOM_INLINE void init_rsc(const YV12_BUFFER_CONFIG *src,
167
                                const AV1_COMMON *cm, const MACROBLOCK *x,
168
                                const LOOP_FILTER_SPEED_FEATURES *lpf_sf,
169
                                int plane, RestUnitSearchInfo *rusi,
170
0
                                YV12_BUFFER_CONFIG *dst, RestSearchCtxt *rsc) {
171
0
  rsc->src = src;
172
0
  rsc->dst = dst;
173
0
  rsc->cm = cm;
174
0
  rsc->x = x;
175
0
  rsc->plane = plane;
176
0
  rsc->rusi = rusi;
177
0
  rsc->lpf_sf = lpf_sf;
178
179
0
  const YV12_BUFFER_CONFIG *dgd = &cm->cur_frame->buf;
180
0
  const int is_uv = plane != AOM_PLANE_Y;
181
0
  rsc->plane_width = src->crop_widths[is_uv];
182
0
  rsc->plane_height = src->crop_heights[is_uv];
183
0
  rsc->src_buffer = src->buffers[plane];
184
0
  rsc->src_stride = src->strides[is_uv];
185
0
  rsc->dgd_buffer = dgd->buffers[plane];
186
0
  rsc->dgd_stride = dgd->strides[is_uv];
187
0
  rsc->tile_rect = av1_whole_frame_rect(cm, is_uv);
188
0
  assert(src->crop_widths[is_uv] == dgd->crop_widths[is_uv]);
189
0
  assert(src->crop_heights[is_uv] == dgd->crop_heights[is_uv]);
190
0
}
191
192
static int64_t try_restoration_unit(const RestSearchCtxt *rsc,
193
                                    const RestorationTileLimits *limits,
194
                                    const AV1PixelRect *tile_rect,
195
0
                                    const RestorationUnitInfo *rui) {
196
0
  const AV1_COMMON *const cm = rsc->cm;
197
0
  const int plane = rsc->plane;
198
0
  const int is_uv = plane > 0;
199
0
  const RestorationInfo *rsi = &cm->rst_info[plane];
200
0
  RestorationLineBuffers rlbs;
201
0
  const int bit_depth = cm->seq_params->bit_depth;
202
0
  const int highbd = cm->seq_params->use_highbitdepth;
203
204
0
  const YV12_BUFFER_CONFIG *fts = &cm->cur_frame->buf;
205
  // TODO(yunqing): For now, only use optimized LR filter in decoder. Can be
206
  // also used in encoder.
207
0
  const int optimized_lr = 0;
208
209
0
  av1_loop_restoration_filter_unit(
210
0
      limits, rui, &rsi->boundaries, &rlbs, tile_rect, rsc->tile_stripe0,
211
0
      is_uv && cm->seq_params->subsampling_x,
212
0
      is_uv && cm->seq_params->subsampling_y, highbd, bit_depth,
213
0
      fts->buffers[plane], fts->strides[is_uv], rsc->dst->buffers[plane],
214
0
      rsc->dst->strides[is_uv], cm->rst_tmpbuf, optimized_lr);
215
216
0
  return sse_restoration_unit(limits, rsc->src, rsc->dst, plane, highbd);
217
0
}
218
219
int64_t av1_lowbd_pixel_proj_error_c(const uint8_t *src8, int width, int height,
220
                                     int src_stride, const uint8_t *dat8,
221
                                     int dat_stride, int32_t *flt0,
222
                                     int flt0_stride, int32_t *flt1,
223
                                     int flt1_stride, int xq[2],
224
0
                                     const sgr_params_type *params) {
225
0
  int i, j;
226
0
  const uint8_t *src = src8;
227
0
  const uint8_t *dat = dat8;
228
0
  int64_t err = 0;
229
0
  if (params->r[0] > 0 && params->r[1] > 0) {
230
0
    for (i = 0; i < height; ++i) {
231
0
      for (j = 0; j < width; ++j) {
232
0
        assert(flt1[j] < (1 << 15) && flt1[j] > -(1 << 15));
233
0
        assert(flt0[j] < (1 << 15) && flt0[j] > -(1 << 15));
234
0
        const int32_t u = (int32_t)(dat[j] << SGRPROJ_RST_BITS);
235
0
        int32_t v = u << SGRPROJ_PRJ_BITS;
236
0
        v += xq[0] * (flt0[j] - u) + xq[1] * (flt1[j] - u);
237
0
        const int32_t e =
238
0
            ROUND_POWER_OF_TWO(v, SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS) - src[j];
239
0
        err += ((int64_t)e * e);
240
0
      }
241
0
      dat += dat_stride;
242
0
      src += src_stride;
243
0
      flt0 += flt0_stride;
244
0
      flt1 += flt1_stride;
245
0
    }
246
0
  } else if (params->r[0] > 0) {
247
0
    for (i = 0; i < height; ++i) {
248
0
      for (j = 0; j < width; ++j) {
249
0
        assert(flt0[j] < (1 << 15) && flt0[j] > -(1 << 15));
250
0
        const int32_t u = (int32_t)(dat[j] << SGRPROJ_RST_BITS);
251
0
        int32_t v = u << SGRPROJ_PRJ_BITS;
252
0
        v += xq[0] * (flt0[j] - u);
253
0
        const int32_t e =
254
0
            ROUND_POWER_OF_TWO(v, SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS) - src[j];
255
0
        err += ((int64_t)e * e);
256
0
      }
257
0
      dat += dat_stride;
258
0
      src += src_stride;
259
0
      flt0 += flt0_stride;
260
0
    }
261
0
  } else if (params->r[1] > 0) {
262
0
    for (i = 0; i < height; ++i) {
263
0
      for (j = 0; j < width; ++j) {
264
0
        assert(flt1[j] < (1 << 15) && flt1[j] > -(1 << 15));
265
0
        const int32_t u = (int32_t)(dat[j] << SGRPROJ_RST_BITS);
266
0
        int32_t v = u << SGRPROJ_PRJ_BITS;
267
0
        v += xq[1] * (flt1[j] - u);
268
0
        const int32_t e =
269
0
            ROUND_POWER_OF_TWO(v, SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS) - src[j];
270
0
        err += ((int64_t)e * e);
271
0
      }
272
0
      dat += dat_stride;
273
0
      src += src_stride;
274
0
      flt1 += flt1_stride;
275
0
    }
276
0
  } else {
277
0
    for (i = 0; i < height; ++i) {
278
0
      for (j = 0; j < width; ++j) {
279
0
        const int32_t e = (int32_t)(dat[j]) - src[j];
280
0
        err += ((int64_t)e * e);
281
0
      }
282
0
      dat += dat_stride;
283
0
      src += src_stride;
284
0
    }
285
0
  }
286
287
0
  return err;
288
0
}
289
290
#if CONFIG_AV1_HIGHBITDEPTH
291
int64_t av1_highbd_pixel_proj_error_c(const uint8_t *src8, int width,
292
                                      int height, int src_stride,
293
                                      const uint8_t *dat8, int dat_stride,
294
                                      int32_t *flt0, int flt0_stride,
295
                                      int32_t *flt1, int flt1_stride, int xq[2],
296
0
                                      const sgr_params_type *params) {
297
0
  const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
298
0
  const uint16_t *dat = CONVERT_TO_SHORTPTR(dat8);
299
0
  int i, j;
300
0
  int64_t err = 0;
301
0
  const int32_t half = 1 << (SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS - 1);
302
0
  if (params->r[0] > 0 && params->r[1] > 0) {
303
0
    int xq0 = xq[0];
304
0
    int xq1 = xq[1];
305
0
    for (i = 0; i < height; ++i) {
306
0
      for (j = 0; j < width; ++j) {
307
0
        const int32_t d = dat[j];
308
0
        const int32_t s = src[j];
309
0
        const int32_t u = (int32_t)(d << SGRPROJ_RST_BITS);
310
0
        int32_t v0 = flt0[j] - u;
311
0
        int32_t v1 = flt1[j] - u;
312
0
        int32_t v = half;
313
0
        v += xq0 * v0;
314
0
        v += xq1 * v1;
315
0
        const int32_t e = (v >> (SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS)) + d - s;
316
0
        err += ((int64_t)e * e);
317
0
      }
318
0
      dat += dat_stride;
319
0
      flt0 += flt0_stride;
320
0
      flt1 += flt1_stride;
321
0
      src += src_stride;
322
0
    }
323
0
  } else if (params->r[0] > 0 || params->r[1] > 0) {
324
0
    int exq;
325
0
    int32_t *flt;
326
0
    int flt_stride;
327
0
    if (params->r[0] > 0) {
328
0
      exq = xq[0];
329
0
      flt = flt0;
330
0
      flt_stride = flt0_stride;
331
0
    } else {
332
0
      exq = xq[1];
333
0
      flt = flt1;
334
0
      flt_stride = flt1_stride;
335
0
    }
336
0
    for (i = 0; i < height; ++i) {
337
0
      for (j = 0; j < width; ++j) {
338
0
        const int32_t d = dat[j];
339
0
        const int32_t s = src[j];
340
0
        const int32_t u = (int32_t)(d << SGRPROJ_RST_BITS);
341
0
        int32_t v = half;
342
0
        v += exq * (flt[j] - u);
343
0
        const int32_t e = (v >> (SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS)) + d - s;
344
0
        err += ((int64_t)e * e);
345
0
      }
346
0
      dat += dat_stride;
347
0
      flt += flt_stride;
348
0
      src += src_stride;
349
0
    }
350
0
  } else {
351
0
    for (i = 0; i < height; ++i) {
352
0
      for (j = 0; j < width; ++j) {
353
0
        const int32_t d = dat[j];
354
0
        const int32_t s = src[j];
355
0
        const int32_t e = d - s;
356
0
        err += ((int64_t)e * e);
357
0
      }
358
0
      dat += dat_stride;
359
0
      src += src_stride;
360
0
    }
361
0
  }
362
0
  return err;
363
0
}
364
#endif  // CONFIG_AV1_HIGHBITDEPTH
365
366
static int64_t get_pixel_proj_error(const uint8_t *src8, int width, int height,
367
                                    int src_stride, const uint8_t *dat8,
368
                                    int dat_stride, int use_highbitdepth,
369
                                    int32_t *flt0, int flt0_stride,
370
                                    int32_t *flt1, int flt1_stride, int *xqd,
371
0
                                    const sgr_params_type *params) {
372
0
  int xq[2];
373
0
  av1_decode_xq(xqd, xq, params);
374
375
0
#if CONFIG_AV1_HIGHBITDEPTH
376
0
  if (use_highbitdepth) {
377
0
    return av1_highbd_pixel_proj_error(src8, width, height, src_stride, dat8,
378
0
                                       dat_stride, flt0, flt0_stride, flt1,
379
0
                                       flt1_stride, xq, params);
380
381
0
  } else {
382
0
    return av1_lowbd_pixel_proj_error(src8, width, height, src_stride, dat8,
383
0
                                      dat_stride, flt0, flt0_stride, flt1,
384
0
                                      flt1_stride, xq, params);
385
0
  }
386
#else
387
  (void)use_highbitdepth;
388
  return av1_lowbd_pixel_proj_error(src8, width, height, src_stride, dat8,
389
                                    dat_stride, flt0, flt0_stride, flt1,
390
                                    flt1_stride, xq, params);
391
#endif
392
0
}
393
394
#define USE_SGRPROJ_REFINEMENT_SEARCH 1
395
static int64_t finer_search_pixel_proj_error(
396
    const uint8_t *src8, int width, int height, int src_stride,
397
    const uint8_t *dat8, int dat_stride, int use_highbitdepth, int32_t *flt0,
398
    int flt0_stride, int32_t *flt1, int flt1_stride, int start_step, int *xqd,
399
0
    const sgr_params_type *params) {
400
0
  int64_t err = get_pixel_proj_error(
401
0
      src8, width, height, src_stride, dat8, dat_stride, use_highbitdepth, flt0,
402
0
      flt0_stride, flt1, flt1_stride, xqd, params);
403
0
  (void)start_step;
404
0
#if USE_SGRPROJ_REFINEMENT_SEARCH
405
0
  int64_t err2;
406
0
  int tap_min[] = { SGRPROJ_PRJ_MIN0, SGRPROJ_PRJ_MIN1 };
407
0
  int tap_max[] = { SGRPROJ_PRJ_MAX0, SGRPROJ_PRJ_MAX1 };
408
0
  for (int s = start_step; s >= 1; s >>= 1) {
409
0
    for (int p = 0; p < 2; ++p) {
410
0
      if ((params->r[0] == 0 && p == 0) || (params->r[1] == 0 && p == 1)) {
411
0
        continue;
412
0
      }
413
0
      int skip = 0;
414
0
      do {
415
0
        if (xqd[p] - s >= tap_min[p]) {
416
0
          xqd[p] -= s;
417
0
          err2 =
418
0
              get_pixel_proj_error(src8, width, height, src_stride, dat8,
419
0
                                   dat_stride, use_highbitdepth, flt0,
420
0
                                   flt0_stride, flt1, flt1_stride, xqd, params);
421
0
          if (err2 > err) {
422
0
            xqd[p] += s;
423
0
          } else {
424
0
            err = err2;
425
0
            skip = 1;
426
            // At the highest step size continue moving in the same direction
427
0
            if (s == start_step) continue;
428
0
          }
429
0
        }
430
0
        break;
431
0
      } while (1);
432
0
      if (skip) break;
433
0
      do {
434
0
        if (xqd[p] + s <= tap_max[p]) {
435
0
          xqd[p] += s;
436
0
          err2 =
437
0
              get_pixel_proj_error(src8, width, height, src_stride, dat8,
438
0
                                   dat_stride, use_highbitdepth, flt0,
439
0
                                   flt0_stride, flt1, flt1_stride, xqd, params);
440
0
          if (err2 > err) {
441
0
            xqd[p] -= s;
442
0
          } else {
443
0
            err = err2;
444
            // At the highest step size continue moving in the same direction
445
0
            if (s == start_step) continue;
446
0
          }
447
0
        }
448
0
        break;
449
0
      } while (1);
450
0
    }
451
0
  }
452
0
#endif  // USE_SGRPROJ_REFINEMENT_SEARCH
453
0
  return err;
454
0
}
455
456
0
static int64_t signed_rounded_divide(int64_t dividend, int64_t divisor) {
457
0
  if (dividend < 0)
458
0
    return (dividend - divisor / 2) / divisor;
459
0
  else
460
0
    return (dividend + divisor / 2) / divisor;
461
0
}
462
463
static AOM_INLINE void calc_proj_params_r0_r1_c(
464
    const uint8_t *src8, int width, int height, int src_stride,
465
    const uint8_t *dat8, int dat_stride, int32_t *flt0, int flt0_stride,
466
0
    int32_t *flt1, int flt1_stride, int64_t H[2][2], int64_t C[2]) {
467
0
  const int size = width * height;
468
0
  const uint8_t *src = src8;
469
0
  const uint8_t *dat = dat8;
470
0
  for (int i = 0; i < height; ++i) {
471
0
    for (int j = 0; j < width; ++j) {
472
0
      const int32_t u = (int32_t)(dat[i * dat_stride + j] << SGRPROJ_RST_BITS);
473
0
      const int32_t s =
474
0
          (int32_t)(src[i * src_stride + j] << SGRPROJ_RST_BITS) - u;
475
0
      const int32_t f1 = (int32_t)flt0[i * flt0_stride + j] - u;
476
0
      const int32_t f2 = (int32_t)flt1[i * flt1_stride + j] - u;
477
0
      H[0][0] += (int64_t)f1 * f1;
478
0
      H[1][1] += (int64_t)f2 * f2;
479
0
      H[0][1] += (int64_t)f1 * f2;
480
0
      C[0] += (int64_t)f1 * s;
481
0
      C[1] += (int64_t)f2 * s;
482
0
    }
483
0
  }
484
0
  H[0][0] /= size;
485
0
  H[0][1] /= size;
486
0
  H[1][1] /= size;
487
0
  H[1][0] = H[0][1];
488
0
  C[0] /= size;
489
0
  C[1] /= size;
490
0
}
491
492
static AOM_INLINE void calc_proj_params_r0_r1_high_bd_c(
493
    const uint8_t *src8, int width, int height, int src_stride,
494
    const uint8_t *dat8, int dat_stride, int32_t *flt0, int flt0_stride,
495
0
    int32_t *flt1, int flt1_stride, int64_t H[2][2], int64_t C[2]) {
496
0
  const int size = width * height;
497
0
  const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
498
0
  const uint16_t *dat = CONVERT_TO_SHORTPTR(dat8);
499
0
  for (int i = 0; i < height; ++i) {
500
0
    for (int j = 0; j < width; ++j) {
501
0
      const int32_t u = (int32_t)(dat[i * dat_stride + j] << SGRPROJ_RST_BITS);
502
0
      const int32_t s =
503
0
          (int32_t)(src[i * src_stride + j] << SGRPROJ_RST_BITS) - u;
504
0
      const int32_t f1 = (int32_t)flt0[i * flt0_stride + j] - u;
505
0
      const int32_t f2 = (int32_t)flt1[i * flt1_stride + j] - u;
506
0
      H[0][0] += (int64_t)f1 * f1;
507
0
      H[1][1] += (int64_t)f2 * f2;
508
0
      H[0][1] += (int64_t)f1 * f2;
509
0
      C[0] += (int64_t)f1 * s;
510
0
      C[1] += (int64_t)f2 * s;
511
0
    }
512
0
  }
513
0
  H[0][0] /= size;
514
0
  H[0][1] /= size;
515
0
  H[1][1] /= size;
516
0
  H[1][0] = H[0][1];
517
0
  C[0] /= size;
518
0
  C[1] /= size;
519
0
}
520
521
static AOM_INLINE void calc_proj_params_r0_c(const uint8_t *src8, int width,
522
                                             int height, int src_stride,
523
                                             const uint8_t *dat8,
524
                                             int dat_stride, int32_t *flt0,
525
                                             int flt0_stride, int64_t H[2][2],
526
0
                                             int64_t C[2]) {
527
0
  const int size = width * height;
528
0
  const uint8_t *src = src8;
529
0
  const uint8_t *dat = dat8;
530
0
  for (int i = 0; i < height; ++i) {
531
0
    for (int j = 0; j < width; ++j) {
532
0
      const int32_t u = (int32_t)(dat[i * dat_stride + j] << SGRPROJ_RST_BITS);
533
0
      const int32_t s =
534
0
          (int32_t)(src[i * src_stride + j] << SGRPROJ_RST_BITS) - u;
535
0
      const int32_t f1 = (int32_t)flt0[i * flt0_stride + j] - u;
536
0
      H[0][0] += (int64_t)f1 * f1;
537
0
      C[0] += (int64_t)f1 * s;
538
0
    }
539
0
  }
540
0
  H[0][0] /= size;
541
0
  C[0] /= size;
542
0
}
543
544
static AOM_INLINE void calc_proj_params_r0_high_bd_c(
545
    const uint8_t *src8, int width, int height, int src_stride,
546
    const uint8_t *dat8, int dat_stride, int32_t *flt0, int flt0_stride,
547
0
    int64_t H[2][2], int64_t C[2]) {
548
0
  const int size = width * height;
549
0
  const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
550
0
  const uint16_t *dat = CONVERT_TO_SHORTPTR(dat8);
551
0
  for (int i = 0; i < height; ++i) {
552
0
    for (int j = 0; j < width; ++j) {
553
0
      const int32_t u = (int32_t)(dat[i * dat_stride + j] << SGRPROJ_RST_BITS);
554
0
      const int32_t s =
555
0
          (int32_t)(src[i * src_stride + j] << SGRPROJ_RST_BITS) - u;
556
0
      const int32_t f1 = (int32_t)flt0[i * flt0_stride + j] - u;
557
0
      H[0][0] += (int64_t)f1 * f1;
558
0
      C[0] += (int64_t)f1 * s;
559
0
    }
560
0
  }
561
0
  H[0][0] /= size;
562
0
  C[0] /= size;
563
0
}
564
565
static AOM_INLINE void calc_proj_params_r1_c(const uint8_t *src8, int width,
566
                                             int height, int src_stride,
567
                                             const uint8_t *dat8,
568
                                             int dat_stride, int32_t *flt1,
569
                                             int flt1_stride, int64_t H[2][2],
570
0
                                             int64_t C[2]) {
571
0
  const int size = width * height;
572
0
  const uint8_t *src = src8;
573
0
  const uint8_t *dat = dat8;
574
0
  for (int i = 0; i < height; ++i) {
575
0
    for (int j = 0; j < width; ++j) {
576
0
      const int32_t u = (int32_t)(dat[i * dat_stride + j] << SGRPROJ_RST_BITS);
577
0
      const int32_t s =
578
0
          (int32_t)(src[i * src_stride + j] << SGRPROJ_RST_BITS) - u;
579
0
      const int32_t f2 = (int32_t)flt1[i * flt1_stride + j] - u;
580
0
      H[1][1] += (int64_t)f2 * f2;
581
0
      C[1] += (int64_t)f2 * s;
582
0
    }
583
0
  }
584
0
  H[1][1] /= size;
585
0
  C[1] /= size;
586
0
}
587
588
static AOM_INLINE void calc_proj_params_r1_high_bd_c(
589
    const uint8_t *src8, int width, int height, int src_stride,
590
    const uint8_t *dat8, int dat_stride, int32_t *flt1, int flt1_stride,
591
0
    int64_t H[2][2], int64_t C[2]) {
592
0
  const int size = width * height;
593
0
  const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
594
0
  const uint16_t *dat = CONVERT_TO_SHORTPTR(dat8);
595
0
  for (int i = 0; i < height; ++i) {
596
0
    for (int j = 0; j < width; ++j) {
597
0
      const int32_t u = (int32_t)(dat[i * dat_stride + j] << SGRPROJ_RST_BITS);
598
0
      const int32_t s =
599
0
          (int32_t)(src[i * src_stride + j] << SGRPROJ_RST_BITS) - u;
600
0
      const int32_t f2 = (int32_t)flt1[i * flt1_stride + j] - u;
601
0
      H[1][1] += (int64_t)f2 * f2;
602
0
      C[1] += (int64_t)f2 * s;
603
0
    }
604
0
  }
605
0
  H[1][1] /= size;
606
0
  C[1] /= size;
607
0
}
608
609
// The function calls 3 subfunctions for the following cases :
610
// 1) When params->r[0] > 0 and params->r[1] > 0. In this case all elements
611
// of C and H need to be computed.
612
// 2) When only params->r[0] > 0. In this case only H[0][0] and C[0] are
613
// non-zero and need to be computed.
614
// 3) When only params->r[1] > 0. In this case only H[1][1] and C[1] are
615
// non-zero and need to be computed.
616
void av1_calc_proj_params_c(const uint8_t *src8, int width, int height,
617
                            int src_stride, const uint8_t *dat8, int dat_stride,
618
                            int32_t *flt0, int flt0_stride, int32_t *flt1,
619
                            int flt1_stride, int64_t H[2][2], int64_t C[2],
620
0
                            const sgr_params_type *params) {
621
0
  if ((params->r[0] > 0) && (params->r[1] > 0)) {
622
0
    calc_proj_params_r0_r1_c(src8, width, height, src_stride, dat8, dat_stride,
623
0
                             flt0, flt0_stride, flt1, flt1_stride, H, C);
624
0
  } else if (params->r[0] > 0) {
625
0
    calc_proj_params_r0_c(src8, width, height, src_stride, dat8, dat_stride,
626
0
                          flt0, flt0_stride, H, C);
627
0
  } else if (params->r[1] > 0) {
628
0
    calc_proj_params_r1_c(src8, width, height, src_stride, dat8, dat_stride,
629
0
                          flt1, flt1_stride, H, C);
630
0
  }
631
0
}
632
633
void av1_calc_proj_params_high_bd_c(const uint8_t *src8, int width, int height,
634
                                    int src_stride, const uint8_t *dat8,
635
                                    int dat_stride, int32_t *flt0,
636
                                    int flt0_stride, int32_t *flt1,
637
                                    int flt1_stride, int64_t H[2][2],
638
                                    int64_t C[2],
639
0
                                    const sgr_params_type *params) {
640
0
  if ((params->r[0] > 0) && (params->r[1] > 0)) {
641
0
    calc_proj_params_r0_r1_high_bd_c(src8, width, height, src_stride, dat8,
642
0
                                     dat_stride, flt0, flt0_stride, flt1,
643
0
                                     flt1_stride, H, C);
644
0
  } else if (params->r[0] > 0) {
645
0
    calc_proj_params_r0_high_bd_c(src8, width, height, src_stride, dat8,
646
0
                                  dat_stride, flt0, flt0_stride, H, C);
647
0
  } else if (params->r[1] > 0) {
648
0
    calc_proj_params_r1_high_bd_c(src8, width, height, src_stride, dat8,
649
0
                                  dat_stride, flt1, flt1_stride, H, C);
650
0
  }
651
0
}
652
653
static AOM_INLINE void get_proj_subspace(const uint8_t *src8, int width,
654
                                         int height, int src_stride,
655
                                         const uint8_t *dat8, int dat_stride,
656
                                         int use_highbitdepth, int32_t *flt0,
657
                                         int flt0_stride, int32_t *flt1,
658
                                         int flt1_stride, int *xq,
659
0
                                         const sgr_params_type *params) {
660
0
  int64_t H[2][2] = { { 0, 0 }, { 0, 0 } };
661
0
  int64_t C[2] = { 0, 0 };
662
663
  // Default values to be returned if the problem becomes ill-posed
664
0
  xq[0] = 0;
665
0
  xq[1] = 0;
666
667
0
  if (!use_highbitdepth) {
668
0
    if ((width & 0x7) == 0) {
669
0
      av1_calc_proj_params(src8, width, height, src_stride, dat8, dat_stride,
670
0
                           flt0, flt0_stride, flt1, flt1_stride, H, C, params);
671
0
    } else {
672
0
      av1_calc_proj_params_c(src8, width, height, src_stride, dat8, dat_stride,
673
0
                             flt0, flt0_stride, flt1, flt1_stride, H, C,
674
0
                             params);
675
0
    }
676
0
  }
677
0
#if CONFIG_AV1_HIGHBITDEPTH
678
0
  else {  // NOLINT
679
0
    if ((width & 0x7) == 0) {
680
0
      av1_calc_proj_params_high_bd(src8, width, height, src_stride, dat8,
681
0
                                   dat_stride, flt0, flt0_stride, flt1,
682
0
                                   flt1_stride, H, C, params);
683
0
    } else {
684
0
      av1_calc_proj_params_high_bd_c(src8, width, height, src_stride, dat8,
685
0
                                     dat_stride, flt0, flt0_stride, flt1,
686
0
                                     flt1_stride, H, C, params);
687
0
    }
688
0
  }
689
0
#endif
690
691
0
  if (params->r[0] == 0) {
692
    // H matrix is now only the scalar H[1][1]
693
    // C vector is now only the scalar C[1]
694
0
    const int64_t Det = H[1][1];
695
0
    if (Det == 0) return;  // ill-posed, return default values
696
0
    xq[0] = 0;
697
0
    xq[1] = (int)signed_rounded_divide(C[1] * (1 << SGRPROJ_PRJ_BITS), Det);
698
0
  } else if (params->r[1] == 0) {
699
    // H matrix is now only the scalar H[0][0]
700
    // C vector is now only the scalar C[0]
701
0
    const int64_t Det = H[0][0];
702
0
    if (Det == 0) return;  // ill-posed, return default values
703
0
    xq[0] = (int)signed_rounded_divide(C[0] * (1 << SGRPROJ_PRJ_BITS), Det);
704
0
    xq[1] = 0;
705
0
  } else {
706
0
    const int64_t Det = H[0][0] * H[1][1] - H[0][1] * H[1][0];
707
0
    if (Det == 0) return;  // ill-posed, return default values
708
709
    // If scaling up dividend would overflow, instead scale down the divisor
710
0
    const int64_t div1 = H[1][1] * C[0] - H[0][1] * C[1];
711
0
    if ((div1 > 0 && INT64_MAX / (1 << SGRPROJ_PRJ_BITS) < div1) ||
712
0
        (div1 < 0 && INT64_MIN / (1 << SGRPROJ_PRJ_BITS) > div1))
713
0
      xq[0] = (int)signed_rounded_divide(div1, Det / (1 << SGRPROJ_PRJ_BITS));
714
0
    else
715
0
      xq[0] = (int)signed_rounded_divide(div1 * (1 << SGRPROJ_PRJ_BITS), Det);
716
717
0
    const int64_t div2 = H[0][0] * C[1] - H[1][0] * C[0];
718
0
    if ((div2 > 0 && INT64_MAX / (1 << SGRPROJ_PRJ_BITS) < div2) ||
719
0
        (div2 < 0 && INT64_MIN / (1 << SGRPROJ_PRJ_BITS) > div2))
720
0
      xq[1] = (int)signed_rounded_divide(div2, Det / (1 << SGRPROJ_PRJ_BITS));
721
0
    else
722
0
      xq[1] = (int)signed_rounded_divide(div2 * (1 << SGRPROJ_PRJ_BITS), Det);
723
0
  }
724
0
}
725
726
static AOM_INLINE void encode_xq(int *xq, int *xqd,
727
0
                                 const sgr_params_type *params) {
728
0
  if (params->r[0] == 0) {
729
0
    xqd[0] = 0;
730
0
    xqd[1] = clamp((1 << SGRPROJ_PRJ_BITS) - xq[1], SGRPROJ_PRJ_MIN1,
731
0
                   SGRPROJ_PRJ_MAX1);
732
0
  } else if (params->r[1] == 0) {
733
0
    xqd[0] = clamp(xq[0], SGRPROJ_PRJ_MIN0, SGRPROJ_PRJ_MAX0);
734
0
    xqd[1] = clamp((1 << SGRPROJ_PRJ_BITS) - xqd[0], SGRPROJ_PRJ_MIN1,
735
0
                   SGRPROJ_PRJ_MAX1);
736
0
  } else {
737
0
    xqd[0] = clamp(xq[0], SGRPROJ_PRJ_MIN0, SGRPROJ_PRJ_MAX0);
738
0
    xqd[1] = clamp((1 << SGRPROJ_PRJ_BITS) - xqd[0] - xq[1], SGRPROJ_PRJ_MIN1,
739
0
                   SGRPROJ_PRJ_MAX1);
740
0
  }
741
0
}
742
743
// Apply the self-guided filter across an entire restoration unit.
744
static AOM_INLINE void apply_sgr(int sgr_params_idx, const uint8_t *dat8,
745
                                 int width, int height, int dat_stride,
746
                                 int use_highbd, int bit_depth, int pu_width,
747
                                 int pu_height, int32_t *flt0, int32_t *flt1,
748
0
                                 int flt_stride) {
749
0
  for (int i = 0; i < height; i += pu_height) {
750
0
    const int h = AOMMIN(pu_height, height - i);
751
0
    int32_t *flt0_row = flt0 + i * flt_stride;
752
0
    int32_t *flt1_row = flt1 + i * flt_stride;
753
0
    const uint8_t *dat8_row = dat8 + i * dat_stride;
754
755
    // Iterate over the stripe in blocks of width pu_width
756
0
    for (int j = 0; j < width; j += pu_width) {
757
0
      const int w = AOMMIN(pu_width, width - j);
758
0
      const int ret = av1_selfguided_restoration(
759
0
          dat8_row + j, w, h, dat_stride, flt0_row + j, flt1_row + j,
760
0
          flt_stride, sgr_params_idx, bit_depth, use_highbd);
761
0
      (void)ret;
762
0
      assert(!ret);
763
0
    }
764
0
  }
765
0
}
766
767
static AOM_INLINE void compute_sgrproj_err(
768
    const uint8_t *dat8, const int width, const int height,
769
    const int dat_stride, const uint8_t *src8, const int src_stride,
770
    const int use_highbitdepth, const int bit_depth, const int pu_width,
771
    const int pu_height, const int ep, int32_t *flt0, int32_t *flt1,
772
0
    const int flt_stride, int *exqd, int64_t *err) {
773
0
  int exq[2];
774
0
  apply_sgr(ep, dat8, width, height, dat_stride, use_highbitdepth, bit_depth,
775
0
            pu_width, pu_height, flt0, flt1, flt_stride);
776
0
  const sgr_params_type *const params = &av1_sgr_params[ep];
777
0
  get_proj_subspace(src8, width, height, src_stride, dat8, dat_stride,
778
0
                    use_highbitdepth, flt0, flt_stride, flt1, flt_stride, exq,
779
0
                    params);
780
0
  encode_xq(exq, exqd, params);
781
0
  *err = finer_search_pixel_proj_error(
782
0
      src8, width, height, src_stride, dat8, dat_stride, use_highbitdepth, flt0,
783
0
      flt_stride, flt1, flt_stride, 2, exqd, params);
784
0
}
785
786
static AOM_INLINE void get_best_error(int64_t *besterr, const int64_t err,
787
                                      const int *exqd, int *bestxqd,
788
0
                                      int *bestep, const int ep) {
789
0
  if (*besterr == -1 || err < *besterr) {
790
0
    *bestep = ep;
791
0
    *besterr = err;
792
0
    bestxqd[0] = exqd[0];
793
0
    bestxqd[1] = exqd[1];
794
0
  }
795
0
}
796
797
static SgrprojInfo search_selfguided_restoration(
798
    const uint8_t *dat8, int width, int height, int dat_stride,
799
    const uint8_t *src8, int src_stride, int use_highbitdepth, int bit_depth,
800
0
    int pu_width, int pu_height, int32_t *rstbuf, int enable_sgr_ep_pruning) {
801
0
  int32_t *flt0 = rstbuf;
802
0
  int32_t *flt1 = flt0 + RESTORATION_UNITPELS_MAX;
803
0
  int ep, idx, bestep = 0;
804
0
  int64_t besterr = -1;
805
0
  int exqd[2], bestxqd[2] = { 0, 0 };
806
0
  int flt_stride = ((width + 7) & ~7) + 8;
807
0
  assert(pu_width == (RESTORATION_PROC_UNIT_SIZE >> 1) ||
808
0
         pu_width == RESTORATION_PROC_UNIT_SIZE);
809
0
  assert(pu_height == (RESTORATION_PROC_UNIT_SIZE >> 1) ||
810
0
         pu_height == RESTORATION_PROC_UNIT_SIZE);
811
0
  if (!enable_sgr_ep_pruning) {
812
0
    for (ep = 0; ep < SGRPROJ_PARAMS; ep++) {
813
0
      int64_t err;
814
0
      compute_sgrproj_err(dat8, width, height, dat_stride, src8, src_stride,
815
0
                          use_highbitdepth, bit_depth, pu_width, pu_height, ep,
816
0
                          flt0, flt1, flt_stride, exqd, &err);
817
0
      get_best_error(&besterr, err, exqd, bestxqd, &bestep, ep);
818
0
    }
819
0
  } else {
820
    // evaluate first four seed ep in first group
821
0
    for (idx = 0; idx < SGRPROJ_EP_GRP1_SEARCH_COUNT; idx++) {
822
0
      ep = sgproj_ep_grp1_seed[idx];
823
0
      int64_t err;
824
0
      compute_sgrproj_err(dat8, width, height, dat_stride, src8, src_stride,
825
0
                          use_highbitdepth, bit_depth, pu_width, pu_height, ep,
826
0
                          flt0, flt1, flt_stride, exqd, &err);
827
0
      get_best_error(&besterr, err, exqd, bestxqd, &bestep, ep);
828
0
    }
829
    // evaluate left and right ep of winner in seed ep
830
0
    int bestep_ref = bestep;
831
0
    for (ep = bestep_ref - 1; ep < bestep_ref + 2; ep += 2) {
832
0
      if (ep < SGRPROJ_EP_GRP1_START_IDX || ep > SGRPROJ_EP_GRP1_END_IDX)
833
0
        continue;
834
0
      int64_t err;
835
0
      compute_sgrproj_err(dat8, width, height, dat_stride, src8, src_stride,
836
0
                          use_highbitdepth, bit_depth, pu_width, pu_height, ep,
837
0
                          flt0, flt1, flt_stride, exqd, &err);
838
0
      get_best_error(&besterr, err, exqd, bestxqd, &bestep, ep);
839
0
    }
840
    // evaluate last two group
841
0
    for (idx = 0; idx < SGRPROJ_EP_GRP2_3_SEARCH_COUNT; idx++) {
842
0
      ep = sgproj_ep_grp2_3[idx][bestep];
843
0
      int64_t err;
844
0
      compute_sgrproj_err(dat8, width, height, dat_stride, src8, src_stride,
845
0
                          use_highbitdepth, bit_depth, pu_width, pu_height, ep,
846
0
                          flt0, flt1, flt_stride, exqd, &err);
847
0
      get_best_error(&besterr, err, exqd, bestxqd, &bestep, ep);
848
0
    }
849
0
  }
850
851
0
  SgrprojInfo ret;
852
0
  ret.ep = bestep;
853
0
  ret.xqd[0] = bestxqd[0];
854
0
  ret.xqd[1] = bestxqd[1];
855
0
  return ret;
856
0
}
857
858
static int count_sgrproj_bits(SgrprojInfo *sgrproj_info,
859
0
                              SgrprojInfo *ref_sgrproj_info) {
860
0
  int bits = SGRPROJ_PARAMS_BITS;
861
0
  const sgr_params_type *params = &av1_sgr_params[sgrproj_info->ep];
862
0
  if (params->r[0] > 0)
863
0
    bits += aom_count_primitive_refsubexpfin(
864
0
        SGRPROJ_PRJ_MAX0 - SGRPROJ_PRJ_MIN0 + 1, SGRPROJ_PRJ_SUBEXP_K,
865
0
        ref_sgrproj_info->xqd[0] - SGRPROJ_PRJ_MIN0,
866
0
        sgrproj_info->xqd[0] - SGRPROJ_PRJ_MIN0);
867
0
  if (params->r[1] > 0)
868
0
    bits += aom_count_primitive_refsubexpfin(
869
0
        SGRPROJ_PRJ_MAX1 - SGRPROJ_PRJ_MIN1 + 1, SGRPROJ_PRJ_SUBEXP_K,
870
0
        ref_sgrproj_info->xqd[1] - SGRPROJ_PRJ_MIN1,
871
0
        sgrproj_info->xqd[1] - SGRPROJ_PRJ_MIN1);
872
0
  return bits;
873
0
}
874
875
static AOM_INLINE void search_sgrproj(const RestorationTileLimits *limits,
876
                                      const AV1PixelRect *tile,
877
                                      int rest_unit_idx, void *priv,
878
                                      int32_t *tmpbuf,
879
0
                                      RestorationLineBuffers *rlbs) {
880
0
  (void)rlbs;
881
0
  RestSearchCtxt *rsc = (RestSearchCtxt *)priv;
882
0
  RestUnitSearchInfo *rusi = &rsc->rusi[rest_unit_idx];
883
884
0
  const MACROBLOCK *const x = rsc->x;
885
0
  const AV1_COMMON *const cm = rsc->cm;
886
0
  const int highbd = cm->seq_params->use_highbitdepth;
887
0
  const int bit_depth = cm->seq_params->bit_depth;
888
889
0
  const int64_t bits_none = x->mode_costs.sgrproj_restore_cost[0];
890
  // Prune evaluation of RESTORE_SGRPROJ if 'skip_sgr_eval' is set
891
0
  if (rusi->skip_sgr_eval) {
892
0
    rsc->bits += bits_none;
893
0
    rsc->sse += rusi->sse[RESTORE_NONE];
894
0
    rusi->best_rtype[RESTORE_SGRPROJ - 1] = RESTORE_NONE;
895
0
    rusi->sse[RESTORE_SGRPROJ] = INT64_MAX;
896
0
    return;
897
0
  }
898
899
0
  uint8_t *dgd_start =
900
0
      rsc->dgd_buffer + limits->v_start * rsc->dgd_stride + limits->h_start;
901
0
  const uint8_t *src_start =
902
0
      rsc->src_buffer + limits->v_start * rsc->src_stride + limits->h_start;
903
904
0
  const int is_uv = rsc->plane > 0;
905
0
  const int ss_x = is_uv && cm->seq_params->subsampling_x;
906
0
  const int ss_y = is_uv && cm->seq_params->subsampling_y;
907
0
  const int procunit_width = RESTORATION_PROC_UNIT_SIZE >> ss_x;
908
0
  const int procunit_height = RESTORATION_PROC_UNIT_SIZE >> ss_y;
909
910
0
  rusi->sgrproj = search_selfguided_restoration(
911
0
      dgd_start, limits->h_end - limits->h_start,
912
0
      limits->v_end - limits->v_start, rsc->dgd_stride, src_start,
913
0
      rsc->src_stride, highbd, bit_depth, procunit_width, procunit_height,
914
0
      tmpbuf, rsc->lpf_sf->enable_sgr_ep_pruning);
915
916
0
  RestorationUnitInfo rui;
917
0
  rui.restoration_type = RESTORE_SGRPROJ;
918
0
  rui.sgrproj_info = rusi->sgrproj;
919
920
0
  rusi->sse[RESTORE_SGRPROJ] = try_restoration_unit(rsc, limits, tile, &rui);
921
922
0
  const int64_t bits_sgr = x->mode_costs.sgrproj_restore_cost[1] +
923
0
                           (count_sgrproj_bits(&rusi->sgrproj, &rsc->sgrproj)
924
0
                            << AV1_PROB_COST_SHIFT);
925
0
  double cost_none = RDCOST_DBL_WITH_NATIVE_BD_DIST(
926
0
      x->rdmult, bits_none >> 4, rusi->sse[RESTORE_NONE], bit_depth);
927
0
  double cost_sgr = RDCOST_DBL_WITH_NATIVE_BD_DIST(
928
0
      x->rdmult, bits_sgr >> 4, rusi->sse[RESTORE_SGRPROJ], bit_depth);
929
0
  if (rusi->sgrproj.ep < 10)
930
0
    cost_sgr *=
931
0
        (1 + DUAL_SGR_PENALTY_MULT * rsc->lpf_sf->dual_sgr_penalty_level);
932
933
0
  RestorationType rtype =
934
0
      (cost_sgr < cost_none) ? RESTORE_SGRPROJ : RESTORE_NONE;
935
0
  rusi->best_rtype[RESTORE_SGRPROJ - 1] = rtype;
936
937
0
  rsc->sse += rusi->sse[rtype];
938
0
  rsc->bits += (cost_sgr < cost_none) ? bits_sgr : bits_none;
939
0
  if (cost_sgr < cost_none) rsc->sgrproj = rusi->sgrproj;
940
0
}
941
942
void acc_stat_one_line(const uint8_t *dgd, const uint8_t *src, int dgd_stride,
943
                       int h_start, int h_end, uint8_t avg,
944
                       const int wiener_halfwin, const int wiener_win2,
945
0
                       int32_t *M_int32, int32_t *H_int32, int count) {
946
0
  int j, k, l;
947
0
  int16_t Y[WIENER_WIN2];
948
949
0
  for (j = h_start; j < h_end; j++) {
950
0
    const int16_t X = (int16_t)src[j] - (int16_t)avg;
951
0
    int idx = 0;
952
0
    for (k = -wiener_halfwin; k <= wiener_halfwin; k++) {
953
0
      for (l = -wiener_halfwin; l <= wiener_halfwin; l++) {
954
0
        Y[idx] =
955
0
            (int16_t)dgd[(count + l) * dgd_stride + (j + k)] - (int16_t)avg;
956
0
        idx++;
957
0
      }
958
0
    }
959
0
    assert(idx == wiener_win2);
960
0
    for (k = 0; k < wiener_win2; ++k) {
961
0
      M_int32[k] += (int32_t)Y[k] * X;
962
0
      for (l = k; l < wiener_win2; ++l) {
963
        // H is a symmetric matrix, so we only need to fill out the upper
964
        // triangle here. We can copy it down to the lower triangle outside
965
        // the (i, j) loops.
966
0
        H_int32[k * wiener_win2 + l] += (int32_t)Y[k] * Y[l];
967
0
      }
968
0
    }
969
0
  }
970
0
}
971
972
void av1_compute_stats_c(int wiener_win, const uint8_t *dgd, const uint8_t *src,
973
                         int h_start, int h_end, int v_start, int v_end,
974
                         int dgd_stride, int src_stride, int64_t *M, int64_t *H,
975
0
                         int use_downsampled_wiener_stats) {
976
0
  int i, k, l;
977
0
  const int wiener_win2 = wiener_win * wiener_win;
978
0
  const int wiener_halfwin = (wiener_win >> 1);
979
0
  uint8_t avg = find_average(dgd, h_start, h_end, v_start, v_end, dgd_stride);
980
0
  int32_t M_row[WIENER_WIN2] = { 0 };
981
0
  int32_t H_row[WIENER_WIN2 * WIENER_WIN2] = { 0 };
982
0
  int downsample_factor =
983
0
      use_downsampled_wiener_stats ? WIENER_STATS_DOWNSAMPLE_FACTOR : 1;
984
985
0
  memset(M, 0, sizeof(*M) * wiener_win2);
986
0
  memset(H, 0, sizeof(*H) * wiener_win2 * wiener_win2);
987
988
0
  for (i = v_start; i < v_end; i = i + downsample_factor) {
989
0
    if (use_downsampled_wiener_stats &&
990
0
        (v_end - i < WIENER_STATS_DOWNSAMPLE_FACTOR)) {
991
0
      downsample_factor = v_end - i;
992
0
    }
993
994
0
    memset(M_row, 0, sizeof(int32_t) * WIENER_WIN2);
995
0
    memset(H_row, 0, sizeof(int32_t) * WIENER_WIN2 * WIENER_WIN2);
996
0
    acc_stat_one_line(dgd, src + i * src_stride, dgd_stride, h_start, h_end,
997
0
                      avg, wiener_halfwin, wiener_win2, M_row, H_row, i);
998
999
0
    for (k = 0; k < wiener_win2; ++k) {
1000
      // Scale M matrix based on the downsampling factor
1001
0
      M[k] += ((int64_t)M_row[k] * downsample_factor);
1002
0
      for (l = k; l < wiener_win2; ++l) {
1003
        // H is a symmetric matrix, so we only need to fill out the upper
1004
        // triangle here. We can copy it down to the lower triangle outside
1005
        // the (i, j) loops.
1006
        // Scale H Matrix based on the downsampling factor
1007
0
        H[k * wiener_win2 + l] +=
1008
0
            ((int64_t)H_row[k * wiener_win2 + l] * downsample_factor);
1009
0
      }
1010
0
    }
1011
0
  }
1012
1013
0
  for (k = 0; k < wiener_win2; ++k) {
1014
0
    for (l = k + 1; l < wiener_win2; ++l) {
1015
0
      H[l * wiener_win2 + k] = H[k * wiener_win2 + l];
1016
0
    }
1017
0
  }
1018
0
}
1019
1020
#if CONFIG_AV1_HIGHBITDEPTH
1021
void av1_compute_stats_highbd_c(int wiener_win, const uint8_t *dgd8,
1022
                                const uint8_t *src8, int h_start, int h_end,
1023
                                int v_start, int v_end, int dgd_stride,
1024
                                int src_stride, int64_t *M, int64_t *H,
1025
0
                                aom_bit_depth_t bit_depth) {
1026
0
  int i, j, k, l;
1027
0
  int32_t Y[WIENER_WIN2];
1028
0
  const int wiener_win2 = wiener_win * wiener_win;
1029
0
  const int wiener_halfwin = (wiener_win >> 1);
1030
0
  const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
1031
0
  const uint16_t *dgd = CONVERT_TO_SHORTPTR(dgd8);
1032
0
  uint16_t avg =
1033
0
      find_average_highbd(dgd, h_start, h_end, v_start, v_end, dgd_stride);
1034
1035
0
  uint8_t bit_depth_divider = 1;
1036
0
  if (bit_depth == AOM_BITS_12)
1037
0
    bit_depth_divider = 16;
1038
0
  else if (bit_depth == AOM_BITS_10)
1039
0
    bit_depth_divider = 4;
1040
1041
0
  memset(M, 0, sizeof(*M) * wiener_win2);
1042
0
  memset(H, 0, sizeof(*H) * wiener_win2 * wiener_win2);
1043
0
  for (i = v_start; i < v_end; i++) {
1044
0
    for (j = h_start; j < h_end; j++) {
1045
0
      const int32_t X = (int32_t)src[i * src_stride + j] - (int32_t)avg;
1046
0
      int idx = 0;
1047
0
      for (k = -wiener_halfwin; k <= wiener_halfwin; k++) {
1048
0
        for (l = -wiener_halfwin; l <= wiener_halfwin; l++) {
1049
0
          Y[idx] = (int32_t)dgd[(i + l) * dgd_stride + (j + k)] - (int32_t)avg;
1050
0
          idx++;
1051
0
        }
1052
0
      }
1053
0
      assert(idx == wiener_win2);
1054
0
      for (k = 0; k < wiener_win2; ++k) {
1055
0
        M[k] += (int64_t)Y[k] * X;
1056
0
        for (l = k; l < wiener_win2; ++l) {
1057
          // H is a symmetric matrix, so we only need to fill out the upper
1058
          // triangle here. We can copy it down to the lower triangle outside
1059
          // the (i, j) loops.
1060
0
          H[k * wiener_win2 + l] += (int64_t)Y[k] * Y[l];
1061
0
        }
1062
0
      }
1063
0
    }
1064
0
  }
1065
0
  for (k = 0; k < wiener_win2; ++k) {
1066
0
    M[k] /= bit_depth_divider;
1067
0
    H[k * wiener_win2 + k] /= bit_depth_divider;
1068
0
    for (l = k + 1; l < wiener_win2; ++l) {
1069
0
      H[k * wiener_win2 + l] /= bit_depth_divider;
1070
0
      H[l * wiener_win2 + k] = H[k * wiener_win2 + l];
1071
0
    }
1072
0
  }
1073
0
}
1074
#endif  // CONFIG_AV1_HIGHBITDEPTH
1075
1076
0
static INLINE int wrap_index(int i, int wiener_win) {
1077
0
  const int wiener_halfwin1 = (wiener_win >> 1) + 1;
1078
0
  return (i >= wiener_halfwin1 ? wiener_win - 1 - i : i);
1079
0
}
1080
1081
// Solve linear equations to find Wiener filter tap values
1082
// Taps are output scaled by WIENER_FILT_STEP
1083
static int linsolve_wiener(int n, int64_t *A, int stride, int64_t *b,
1084
0
                           int32_t *x) {
1085
0
  for (int k = 0; k < n - 1; k++) {
1086
    // Partial pivoting: bring the row with the largest pivot to the top
1087
0
    for (int i = n - 1; i > k; i--) {
1088
      // If row i has a better (bigger) pivot than row (i-1), swap them
1089
0
      if (llabs(A[(i - 1) * stride + k]) < llabs(A[i * stride + k])) {
1090
0
        for (int j = 0; j < n; j++) {
1091
0
          const int64_t c = A[i * stride + j];
1092
0
          A[i * stride + j] = A[(i - 1) * stride + j];
1093
0
          A[(i - 1) * stride + j] = c;
1094
0
        }
1095
0
        const int64_t c = b[i];
1096
0
        b[i] = b[i - 1];
1097
0
        b[i - 1] = c;
1098
0
      }
1099
0
    }
1100
    // Forward elimination (convert A to row-echelon form)
1101
0
    for (int i = k; i < n - 1; i++) {
1102
0
      if (A[k * stride + k] == 0) return 0;
1103
0
      const int64_t c = A[(i + 1) * stride + k];
1104
0
      const int64_t cd = A[k * stride + k];
1105
0
      for (int j = 0; j < n; j++) {
1106
0
        A[(i + 1) * stride + j] -= c / 256 * A[k * stride + j] / cd * 256;
1107
0
      }
1108
0
      b[i + 1] -= c * b[k] / cd;
1109
0
    }
1110
0
  }
1111
  // Back-substitution
1112
0
  for (int i = n - 1; i >= 0; i--) {
1113
0
    if (A[i * stride + i] == 0) return 0;
1114
0
    int64_t c = 0;
1115
0
    for (int j = i + 1; j <= n - 1; j++) {
1116
0
      c += A[i * stride + j] * x[j] / WIENER_TAP_SCALE_FACTOR;
1117
0
    }
1118
    // Store filter taps x in scaled form.
1119
0
    x[i] = (int32_t)(WIENER_TAP_SCALE_FACTOR * (b[i] - c) / A[i * stride + i]);
1120
0
  }
1121
1122
0
  return 1;
1123
0
}
1124
1125
// Fix vector b, update vector a
1126
static AOM_INLINE void update_a_sep_sym(int wiener_win, int64_t **Mc,
1127
0
                                        int64_t **Hc, int32_t *a, int32_t *b) {
1128
0
  int i, j;
1129
0
  int32_t S[WIENER_WIN];
1130
0
  int64_t A[WIENER_HALFWIN1], B[WIENER_HALFWIN1 * WIENER_HALFWIN1];
1131
0
  const int wiener_win2 = wiener_win * wiener_win;
1132
0
  const int wiener_halfwin1 = (wiener_win >> 1) + 1;
1133
0
  memset(A, 0, sizeof(A));
1134
0
  memset(B, 0, sizeof(B));
1135
0
  for (i = 0; i < wiener_win; i++) {
1136
0
    for (j = 0; j < wiener_win; ++j) {
1137
0
      const int jj = wrap_index(j, wiener_win);
1138
0
      A[jj] += Mc[i][j] * b[i] / WIENER_TAP_SCALE_FACTOR;
1139
0
    }
1140
0
  }
1141
0
  for (i = 0; i < wiener_win; i++) {
1142
0
    for (j = 0; j < wiener_win; j++) {
1143
0
      int k, l;
1144
0
      for (k = 0; k < wiener_win; ++k) {
1145
0
        for (l = 0; l < wiener_win; ++l) {
1146
0
          const int kk = wrap_index(k, wiener_win);
1147
0
          const int ll = wrap_index(l, wiener_win);
1148
0
          B[ll * wiener_halfwin1 + kk] +=
1149
0
              Hc[j * wiener_win + i][k * wiener_win2 + l] * b[i] /
1150
0
              WIENER_TAP_SCALE_FACTOR * b[j] / WIENER_TAP_SCALE_FACTOR;
1151
0
        }
1152
0
      }
1153
0
    }
1154
0
  }
1155
  // Normalization enforcement in the system of equations itself
1156
0
  for (i = 0; i < wiener_halfwin1 - 1; ++i) {
1157
0
    A[i] -=
1158
0
        A[wiener_halfwin1 - 1] * 2 +
1159
0
        B[i * wiener_halfwin1 + wiener_halfwin1 - 1] -
1160
0
        2 * B[(wiener_halfwin1 - 1) * wiener_halfwin1 + (wiener_halfwin1 - 1)];
1161
0
  }
1162
0
  for (i = 0; i < wiener_halfwin1 - 1; ++i) {
1163
0
    for (j = 0; j < wiener_halfwin1 - 1; ++j) {
1164
0
      B[i * wiener_halfwin1 + j] -=
1165
0
          2 * (B[i * wiener_halfwin1 + (wiener_halfwin1 - 1)] +
1166
0
               B[(wiener_halfwin1 - 1) * wiener_halfwin1 + j] -
1167
0
               2 * B[(wiener_halfwin1 - 1) * wiener_halfwin1 +
1168
0
                     (wiener_halfwin1 - 1)]);
1169
0
    }
1170
0
  }
1171
0
  if (linsolve_wiener(wiener_halfwin1 - 1, B, wiener_halfwin1, A, S)) {
1172
0
    S[wiener_halfwin1 - 1] = WIENER_TAP_SCALE_FACTOR;
1173
0
    for (i = wiener_halfwin1; i < wiener_win; ++i) {
1174
0
      S[i] = S[wiener_win - 1 - i];
1175
0
      S[wiener_halfwin1 - 1] -= 2 * S[i];
1176
0
    }
1177
0
    memcpy(a, S, wiener_win * sizeof(*a));
1178
0
  }
1179
0
}
1180
1181
// Fix vector a, update vector b
1182
static AOM_INLINE void update_b_sep_sym(int wiener_win, int64_t **Mc,
1183
0
                                        int64_t **Hc, int32_t *a, int32_t *b) {
1184
0
  int i, j;
1185
0
  int32_t S[WIENER_WIN];
1186
0
  int64_t A[WIENER_HALFWIN1], B[WIENER_HALFWIN1 * WIENER_HALFWIN1];
1187
0
  const int wiener_win2 = wiener_win * wiener_win;
1188
0
  const int wiener_halfwin1 = (wiener_win >> 1) + 1;
1189
0
  memset(A, 0, sizeof(A));
1190
0
  memset(B, 0, sizeof(B));
1191
0
  for (i = 0; i < wiener_win; i++) {
1192
0
    const int ii = wrap_index(i, wiener_win);
1193
0
    for (j = 0; j < wiener_win; j++) {
1194
0
      A[ii] += Mc[i][j] * a[j] / WIENER_TAP_SCALE_FACTOR;
1195
0
    }
1196
0
  }
1197
1198
0
  for (i = 0; i < wiener_win; i++) {
1199
0
    for (j = 0; j < wiener_win; j++) {
1200
0
      const int ii = wrap_index(i, wiener_win);
1201
0
      const int jj = wrap_index(j, wiener_win);
1202
0
      int k, l;
1203
0
      for (k = 0; k < wiener_win; ++k) {
1204
0
        for (l = 0; l < wiener_win; ++l) {
1205
0
          B[jj * wiener_halfwin1 + ii] +=
1206
0
              Hc[i * wiener_win + j][k * wiener_win2 + l] * a[k] /
1207
0
              WIENER_TAP_SCALE_FACTOR * a[l] / WIENER_TAP_SCALE_FACTOR;
1208
0
        }
1209
0
      }
1210
0
    }
1211
0
  }
1212
  // Normalization enforcement in the system of equations itself
1213
0
  for (i = 0; i < wiener_halfwin1 - 1; ++i) {
1214
0
    A[i] -=
1215
0
        A[wiener_halfwin1 - 1] * 2 +
1216
0
        B[i * wiener_halfwin1 + wiener_halfwin1 - 1] -
1217
0
        2 * B[(wiener_halfwin1 - 1) * wiener_halfwin1 + (wiener_halfwin1 - 1)];
1218
0
  }
1219
0
  for (i = 0; i < wiener_halfwin1 - 1; ++i) {
1220
0
    for (j = 0; j < wiener_halfwin1 - 1; ++j) {
1221
0
      B[i * wiener_halfwin1 + j] -=
1222
0
          2 * (B[i * wiener_halfwin1 + (wiener_halfwin1 - 1)] +
1223
0
               B[(wiener_halfwin1 - 1) * wiener_halfwin1 + j] -
1224
0
               2 * B[(wiener_halfwin1 - 1) * wiener_halfwin1 +
1225
0
                     (wiener_halfwin1 - 1)]);
1226
0
    }
1227
0
  }
1228
0
  if (linsolve_wiener(wiener_halfwin1 - 1, B, wiener_halfwin1, A, S)) {
1229
0
    S[wiener_halfwin1 - 1] = WIENER_TAP_SCALE_FACTOR;
1230
0
    for (i = wiener_halfwin1; i < wiener_win; ++i) {
1231
0
      S[i] = S[wiener_win - 1 - i];
1232
0
      S[wiener_halfwin1 - 1] -= 2 * S[i];
1233
0
    }
1234
0
    memcpy(b, S, wiener_win * sizeof(*b));
1235
0
  }
1236
0
}
1237
1238
static void wiener_decompose_sep_sym(int wiener_win, int64_t *M, int64_t *H,
1239
0
                                     int32_t *a, int32_t *b) {
1240
0
  static const int32_t init_filt[WIENER_WIN] = {
1241
0
    WIENER_FILT_TAP0_MIDV, WIENER_FILT_TAP1_MIDV, WIENER_FILT_TAP2_MIDV,
1242
0
    WIENER_FILT_TAP3_MIDV, WIENER_FILT_TAP2_MIDV, WIENER_FILT_TAP1_MIDV,
1243
0
    WIENER_FILT_TAP0_MIDV,
1244
0
  };
1245
0
  int64_t *Hc[WIENER_WIN2];
1246
0
  int64_t *Mc[WIENER_WIN];
1247
0
  int i, j, iter;
1248
0
  const int plane_off = (WIENER_WIN - wiener_win) >> 1;
1249
0
  const int wiener_win2 = wiener_win * wiener_win;
1250
0
  for (i = 0; i < wiener_win; i++) {
1251
0
    a[i] = b[i] =
1252
0
        WIENER_TAP_SCALE_FACTOR / WIENER_FILT_STEP * init_filt[i + plane_off];
1253
0
  }
1254
0
  for (i = 0; i < wiener_win; i++) {
1255
0
    Mc[i] = M + i * wiener_win;
1256
0
    for (j = 0; j < wiener_win; j++) {
1257
0
      Hc[i * wiener_win + j] =
1258
0
          H + i * wiener_win * wiener_win2 + j * wiener_win;
1259
0
    }
1260
0
  }
1261
1262
0
  iter = 1;
1263
0
  while (iter < NUM_WIENER_ITERS) {
1264
0
    update_a_sep_sym(wiener_win, Mc, Hc, a, b);
1265
0
    update_b_sep_sym(wiener_win, Mc, Hc, a, b);
1266
0
    iter++;
1267
0
  }
1268
0
}
1269
1270
// Computes the function x'*H*x - x'*M for the learned 2D filter x, and compares
1271
// against identity filters; Final score is defined as the difference between
1272
// the function values
1273
static int64_t compute_score(int wiener_win, int64_t *M, int64_t *H,
1274
0
                             InterpKernel vfilt, InterpKernel hfilt) {
1275
0
  int32_t ab[WIENER_WIN * WIENER_WIN];
1276
0
  int16_t a[WIENER_WIN], b[WIENER_WIN];
1277
0
  int64_t P = 0, Q = 0;
1278
0
  int64_t iP = 0, iQ = 0;
1279
0
  int64_t Score, iScore;
1280
0
  int i, k, l;
1281
0
  const int plane_off = (WIENER_WIN - wiener_win) >> 1;
1282
0
  const int wiener_win2 = wiener_win * wiener_win;
1283
1284
0
  a[WIENER_HALFWIN] = b[WIENER_HALFWIN] = WIENER_FILT_STEP;
1285
0
  for (i = 0; i < WIENER_HALFWIN; ++i) {
1286
0
    a[i] = a[WIENER_WIN - i - 1] = vfilt[i];
1287
0
    b[i] = b[WIENER_WIN - i - 1] = hfilt[i];
1288
0
    a[WIENER_HALFWIN] -= 2 * a[i];
1289
0
    b[WIENER_HALFWIN] -= 2 * b[i];
1290
0
  }
1291
0
  memset(ab, 0, sizeof(ab));
1292
0
  for (k = 0; k < wiener_win; ++k) {
1293
0
    for (l = 0; l < wiener_win; ++l)
1294
0
      ab[k * wiener_win + l] = a[l + plane_off] * b[k + plane_off];
1295
0
  }
1296
0
  for (k = 0; k < wiener_win2; ++k) {
1297
0
    P += ab[k] * M[k] / WIENER_FILT_STEP / WIENER_FILT_STEP;
1298
0
    for (l = 0; l < wiener_win2; ++l) {
1299
0
      Q += ab[k] * H[k * wiener_win2 + l] * ab[l] / WIENER_FILT_STEP /
1300
0
           WIENER_FILT_STEP / WIENER_FILT_STEP / WIENER_FILT_STEP;
1301
0
    }
1302
0
  }
1303
0
  Score = Q - 2 * P;
1304
1305
0
  iP = M[wiener_win2 >> 1];
1306
0
  iQ = H[(wiener_win2 >> 1) * wiener_win2 + (wiener_win2 >> 1)];
1307
0
  iScore = iQ - 2 * iP;
1308
1309
0
  return Score - iScore;
1310
0
}
1311
1312
static AOM_INLINE void finalize_sym_filter(int wiener_win, int32_t *f,
1313
0
                                           InterpKernel fi) {
1314
0
  int i;
1315
0
  const int wiener_halfwin = (wiener_win >> 1);
1316
1317
0
  for (i = 0; i < wiener_halfwin; ++i) {
1318
0
    const int64_t dividend = (int64_t)f[i] * WIENER_FILT_STEP;
1319
0
    const int64_t divisor = WIENER_TAP_SCALE_FACTOR;
1320
    // Perform this division with proper rounding rather than truncation
1321
0
    if (dividend < 0) {
1322
0
      fi[i] = (int16_t)((dividend - (divisor / 2)) / divisor);
1323
0
    } else {
1324
0
      fi[i] = (int16_t)((dividend + (divisor / 2)) / divisor);
1325
0
    }
1326
0
  }
1327
  // Specialize for 7-tap filter
1328
0
  if (wiener_win == WIENER_WIN) {
1329
0
    fi[0] = CLIP(fi[0], WIENER_FILT_TAP0_MINV, WIENER_FILT_TAP0_MAXV);
1330
0
    fi[1] = CLIP(fi[1], WIENER_FILT_TAP1_MINV, WIENER_FILT_TAP1_MAXV);
1331
0
    fi[2] = CLIP(fi[2], WIENER_FILT_TAP2_MINV, WIENER_FILT_TAP2_MAXV);
1332
0
  } else {
1333
0
    fi[2] = CLIP(fi[1], WIENER_FILT_TAP2_MINV, WIENER_FILT_TAP2_MAXV);
1334
0
    fi[1] = CLIP(fi[0], WIENER_FILT_TAP1_MINV, WIENER_FILT_TAP1_MAXV);
1335
0
    fi[0] = 0;
1336
0
  }
1337
  // Satisfy filter constraints
1338
0
  fi[WIENER_WIN - 1] = fi[0];
1339
0
  fi[WIENER_WIN - 2] = fi[1];
1340
0
  fi[WIENER_WIN - 3] = fi[2];
1341
  // The central element has an implicit +WIENER_FILT_STEP
1342
0
  fi[3] = -2 * (fi[0] + fi[1] + fi[2]);
1343
0
}
1344
1345
static int count_wiener_bits(int wiener_win, WienerInfo *wiener_info,
1346
0
                             WienerInfo *ref_wiener_info) {
1347
0
  int bits = 0;
1348
0
  if (wiener_win == WIENER_WIN)
1349
0
    bits += aom_count_primitive_refsubexpfin(
1350
0
        WIENER_FILT_TAP0_MAXV - WIENER_FILT_TAP0_MINV + 1,
1351
0
        WIENER_FILT_TAP0_SUBEXP_K,
1352
0
        ref_wiener_info->vfilter[0] - WIENER_FILT_TAP0_MINV,
1353
0
        wiener_info->vfilter[0] - WIENER_FILT_TAP0_MINV);
1354
0
  bits += aom_count_primitive_refsubexpfin(
1355
0
      WIENER_FILT_TAP1_MAXV - WIENER_FILT_TAP1_MINV + 1,
1356
0
      WIENER_FILT_TAP1_SUBEXP_K,
1357
0
      ref_wiener_info->vfilter[1] - WIENER_FILT_TAP1_MINV,
1358
0
      wiener_info->vfilter[1] - WIENER_FILT_TAP1_MINV);
1359
0
  bits += aom_count_primitive_refsubexpfin(
1360
0
      WIENER_FILT_TAP2_MAXV - WIENER_FILT_TAP2_MINV + 1,
1361
0
      WIENER_FILT_TAP2_SUBEXP_K,
1362
0
      ref_wiener_info->vfilter[2] - WIENER_FILT_TAP2_MINV,
1363
0
      wiener_info->vfilter[2] - WIENER_FILT_TAP2_MINV);
1364
0
  if (wiener_win == WIENER_WIN)
1365
0
    bits += aom_count_primitive_refsubexpfin(
1366
0
        WIENER_FILT_TAP0_MAXV - WIENER_FILT_TAP0_MINV + 1,
1367
0
        WIENER_FILT_TAP0_SUBEXP_K,
1368
0
        ref_wiener_info->hfilter[0] - WIENER_FILT_TAP0_MINV,
1369
0
        wiener_info->hfilter[0] - WIENER_FILT_TAP0_MINV);
1370
0
  bits += aom_count_primitive_refsubexpfin(
1371
0
      WIENER_FILT_TAP1_MAXV - WIENER_FILT_TAP1_MINV + 1,
1372
0
      WIENER_FILT_TAP1_SUBEXP_K,
1373
0
      ref_wiener_info->hfilter[1] - WIENER_FILT_TAP1_MINV,
1374
0
      wiener_info->hfilter[1] - WIENER_FILT_TAP1_MINV);
1375
0
  bits += aom_count_primitive_refsubexpfin(
1376
0
      WIENER_FILT_TAP2_MAXV - WIENER_FILT_TAP2_MINV + 1,
1377
0
      WIENER_FILT_TAP2_SUBEXP_K,
1378
0
      ref_wiener_info->hfilter[2] - WIENER_FILT_TAP2_MINV,
1379
0
      wiener_info->hfilter[2] - WIENER_FILT_TAP2_MINV);
1380
0
  return bits;
1381
0
}
1382
1383
#define USE_WIENER_REFINEMENT_SEARCH 1
1384
static int64_t finer_tile_search_wiener(const RestSearchCtxt *rsc,
1385
                                        const RestorationTileLimits *limits,
1386
                                        const AV1PixelRect *tile,
1387
                                        RestorationUnitInfo *rui,
1388
0
                                        int wiener_win) {
1389
0
  const int plane_off = (WIENER_WIN - wiener_win) >> 1;
1390
0
  int64_t err = try_restoration_unit(rsc, limits, tile, rui);
1391
0
#if USE_WIENER_REFINEMENT_SEARCH
1392
0
  int64_t err2;
1393
0
  int tap_min[] = { WIENER_FILT_TAP0_MINV, WIENER_FILT_TAP1_MINV,
1394
0
                    WIENER_FILT_TAP2_MINV };
1395
0
  int tap_max[] = { WIENER_FILT_TAP0_MAXV, WIENER_FILT_TAP1_MAXV,
1396
0
                    WIENER_FILT_TAP2_MAXV };
1397
1398
0
  WienerInfo *plane_wiener = &rui->wiener_info;
1399
1400
  // printf("err  pre = %"PRId64"\n", err);
1401
0
  const int start_step = 4;
1402
0
  for (int s = start_step; s >= 1; s >>= 1) {
1403
0
    for (int p = plane_off; p < WIENER_HALFWIN; ++p) {
1404
0
      int skip = 0;
1405
0
      do {
1406
0
        if (plane_wiener->hfilter[p] - s >= tap_min[p]) {
1407
0
          plane_wiener->hfilter[p] -= s;
1408
0
          plane_wiener->hfilter[WIENER_WIN - p - 1] -= s;
1409
0
          plane_wiener->hfilter[WIENER_HALFWIN] += 2 * s;
1410
0
          err2 = try_restoration_unit(rsc, limits, tile, rui);
1411
0
          if (err2 > err) {
1412
0
            plane_wiener->hfilter[p] += s;
1413
0
            plane_wiener->hfilter[WIENER_WIN - p - 1] += s;
1414
0
            plane_wiener->hfilter[WIENER_HALFWIN] -= 2 * s;
1415
0
          } else {
1416
0
            err = err2;
1417
0
            skip = 1;
1418
            // At the highest step size continue moving in the same direction
1419
0
            if (s == start_step) continue;
1420
0
          }
1421
0
        }
1422
0
        break;
1423
0
      } while (1);
1424
0
      if (skip) break;
1425
0
      do {
1426
0
        if (plane_wiener->hfilter[p] + s <= tap_max[p]) {
1427
0
          plane_wiener->hfilter[p] += s;
1428
0
          plane_wiener->hfilter[WIENER_WIN - p - 1] += s;
1429
0
          plane_wiener->hfilter[WIENER_HALFWIN] -= 2 * s;
1430
0
          err2 = try_restoration_unit(rsc, limits, tile, rui);
1431
0
          if (err2 > err) {
1432
0
            plane_wiener->hfilter[p] -= s;
1433
0
            plane_wiener->hfilter[WIENER_WIN - p - 1] -= s;
1434
0
            plane_wiener->hfilter[WIENER_HALFWIN] += 2 * s;
1435
0
          } else {
1436
0
            err = err2;
1437
            // At the highest step size continue moving in the same direction
1438
0
            if (s == start_step) continue;
1439
0
          }
1440
0
        }
1441
0
        break;
1442
0
      } while (1);
1443
0
    }
1444
0
    for (int p = plane_off; p < WIENER_HALFWIN; ++p) {
1445
0
      int skip = 0;
1446
0
      do {
1447
0
        if (plane_wiener->vfilter[p] - s >= tap_min[p]) {
1448
0
          plane_wiener->vfilter[p] -= s;
1449
0
          plane_wiener->vfilter[WIENER_WIN - p - 1] -= s;
1450
0
          plane_wiener->vfilter[WIENER_HALFWIN] += 2 * s;
1451
0
          err2 = try_restoration_unit(rsc, limits, tile, rui);
1452
0
          if (err2 > err) {
1453
0
            plane_wiener->vfilter[p] += s;
1454
0
            plane_wiener->vfilter[WIENER_WIN - p - 1] += s;
1455
0
            plane_wiener->vfilter[WIENER_HALFWIN] -= 2 * s;
1456
0
          } else {
1457
0
            err = err2;
1458
0
            skip = 1;
1459
            // At the highest step size continue moving in the same direction
1460
0
            if (s == start_step) continue;
1461
0
          }
1462
0
        }
1463
0
        break;
1464
0
      } while (1);
1465
0
      if (skip) break;
1466
0
      do {
1467
0
        if (plane_wiener->vfilter[p] + s <= tap_max[p]) {
1468
0
          plane_wiener->vfilter[p] += s;
1469
0
          plane_wiener->vfilter[WIENER_WIN - p - 1] += s;
1470
0
          plane_wiener->vfilter[WIENER_HALFWIN] -= 2 * s;
1471
0
          err2 = try_restoration_unit(rsc, limits, tile, rui);
1472
0
          if (err2 > err) {
1473
0
            plane_wiener->vfilter[p] -= s;
1474
0
            plane_wiener->vfilter[WIENER_WIN - p - 1] -= s;
1475
0
            plane_wiener->vfilter[WIENER_HALFWIN] += 2 * s;
1476
0
          } else {
1477
0
            err = err2;
1478
            // At the highest step size continue moving in the same direction
1479
0
            if (s == start_step) continue;
1480
0
          }
1481
0
        }
1482
0
        break;
1483
0
      } while (1);
1484
0
    }
1485
0
  }
1486
  // printf("err post = %"PRId64"\n", err);
1487
0
#endif  // USE_WIENER_REFINEMENT_SEARCH
1488
0
  return err;
1489
0
}
1490
1491
static AOM_INLINE void search_wiener(const RestorationTileLimits *limits,
1492
                                     const AV1PixelRect *tile_rect,
1493
                                     int rest_unit_idx, void *priv,
1494
                                     int32_t *tmpbuf,
1495
0
                                     RestorationLineBuffers *rlbs) {
1496
0
  (void)tmpbuf;
1497
0
  (void)rlbs;
1498
0
  RestSearchCtxt *rsc = (RestSearchCtxt *)priv;
1499
0
  RestUnitSearchInfo *rusi = &rsc->rusi[rest_unit_idx];
1500
1501
0
  const MACROBLOCK *const x = rsc->x;
1502
0
  const int64_t bits_none = x->mode_costs.wiener_restore_cost[0];
1503
1504
  // Skip Wiener search for low variance contents
1505
0
  if (rsc->lpf_sf->prune_wiener_based_on_src_var) {
1506
0
    const int scale[3] = { 0, 1, 2 };
1507
    // Obtain the normalized Qscale
1508
0
    const int qs = av1_dc_quant_QTX(rsc->cm->quant_params.base_qindex, 0,
1509
0
                                    rsc->cm->seq_params->bit_depth) >>
1510
0
                   3;
1511
    // Derive threshold as sqr(normalized Qscale) * scale / 16,
1512
0
    const uint64_t thresh =
1513
0
        (qs * qs * scale[rsc->lpf_sf->prune_wiener_based_on_src_var]) >> 4;
1514
0
    const int highbd = rsc->cm->seq_params->use_highbitdepth;
1515
0
    const uint64_t src_var =
1516
0
        var_restoration_unit(limits, rsc->src, rsc->plane, highbd);
1517
    // Do not perform Wiener search if source variance is lower than threshold
1518
    // or if the reconstruction error is zero
1519
0
    int prune_wiener = (src_var < thresh) || (rusi->sse[RESTORE_NONE] == 0);
1520
0
    if (prune_wiener) {
1521
0
      rsc->bits += bits_none;
1522
0
      rsc->sse += rusi->sse[RESTORE_NONE];
1523
0
      rusi->best_rtype[RESTORE_WIENER - 1] = RESTORE_NONE;
1524
0
      rusi->sse[RESTORE_WIENER] = INT64_MAX;
1525
0
      if (rsc->lpf_sf->prune_sgr_based_on_wiener == 2) rusi->skip_sgr_eval = 1;
1526
0
      return;
1527
0
    }
1528
0
  }
1529
1530
0
  const int wiener_win =
1531
0
      (rsc->plane == AOM_PLANE_Y) ? WIENER_WIN : WIENER_WIN_CHROMA;
1532
1533
0
  int reduced_wiener_win = wiener_win;
1534
0
  if (rsc->lpf_sf->reduce_wiener_window_size) {
1535
0
    reduced_wiener_win =
1536
0
        (rsc->plane == AOM_PLANE_Y) ? WIENER_WIN_REDUCED : WIENER_WIN_CHROMA;
1537
0
  }
1538
1539
0
  int64_t M[WIENER_WIN2];
1540
0
  int64_t H[WIENER_WIN2 * WIENER_WIN2];
1541
0
  int32_t vfilter[WIENER_WIN], hfilter[WIENER_WIN];
1542
1543
0
#if CONFIG_AV1_HIGHBITDEPTH
1544
0
  const AV1_COMMON *const cm = rsc->cm;
1545
0
  if (cm->seq_params->use_highbitdepth) {
1546
    // TODO(any) : Add support for use_downsampled_wiener_stats SF in HBD
1547
    // functions
1548
0
    av1_compute_stats_highbd(reduced_wiener_win, rsc->dgd_buffer,
1549
0
                             rsc->src_buffer, limits->h_start, limits->h_end,
1550
0
                             limits->v_start, limits->v_end, rsc->dgd_stride,
1551
0
                             rsc->src_stride, M, H, cm->seq_params->bit_depth);
1552
0
  } else {
1553
0
    av1_compute_stats(reduced_wiener_win, rsc->dgd_buffer, rsc->src_buffer,
1554
0
                      limits->h_start, limits->h_end, limits->v_start,
1555
0
                      limits->v_end, rsc->dgd_stride, rsc->src_stride, M, H,
1556
0
                      rsc->lpf_sf->use_downsampled_wiener_stats);
1557
0
  }
1558
#else
1559
  av1_compute_stats(reduced_wiener_win, rsc->dgd_buffer, rsc->src_buffer,
1560
                    limits->h_start, limits->h_end, limits->v_start,
1561
                    limits->v_end, rsc->dgd_stride, rsc->src_stride, M, H,
1562
                    rsc->lpf_sf->use_downsampled_wiener_stats);
1563
#endif
1564
1565
0
  wiener_decompose_sep_sym(reduced_wiener_win, M, H, vfilter, hfilter);
1566
1567
0
  RestorationUnitInfo rui;
1568
0
  memset(&rui, 0, sizeof(rui));
1569
0
  rui.restoration_type = RESTORE_WIENER;
1570
0
  finalize_sym_filter(reduced_wiener_win, vfilter, rui.wiener_info.vfilter);
1571
0
  finalize_sym_filter(reduced_wiener_win, hfilter, rui.wiener_info.hfilter);
1572
1573
  // Filter score computes the value of the function x'*A*x - x'*b for the
1574
  // learned filter and compares it against identity filer. If there is no
1575
  // reduction in the function, the filter is reverted back to identity
1576
0
  if (compute_score(reduced_wiener_win, M, H, rui.wiener_info.vfilter,
1577
0
                    rui.wiener_info.hfilter) > 0) {
1578
0
    rsc->bits += bits_none;
1579
0
    rsc->sse += rusi->sse[RESTORE_NONE];
1580
0
    rusi->best_rtype[RESTORE_WIENER - 1] = RESTORE_NONE;
1581
0
    rusi->sse[RESTORE_WIENER] = INT64_MAX;
1582
0
    if (rsc->lpf_sf->prune_sgr_based_on_wiener == 2) rusi->skip_sgr_eval = 1;
1583
0
    return;
1584
0
  }
1585
1586
0
  rusi->sse[RESTORE_WIENER] = finer_tile_search_wiener(
1587
0
      rsc, limits, tile_rect, &rui, reduced_wiener_win);
1588
0
  rusi->wiener = rui.wiener_info;
1589
1590
0
  if (reduced_wiener_win != WIENER_WIN) {
1591
0
    assert(rui.wiener_info.vfilter[0] == 0 &&
1592
0
           rui.wiener_info.vfilter[WIENER_WIN - 1] == 0);
1593
0
    assert(rui.wiener_info.hfilter[0] == 0 &&
1594
0
           rui.wiener_info.hfilter[WIENER_WIN - 1] == 0);
1595
0
  }
1596
1597
0
  const int64_t bits_wiener =
1598
0
      x->mode_costs.wiener_restore_cost[1] +
1599
0
      (count_wiener_bits(wiener_win, &rusi->wiener, &rsc->wiener)
1600
0
       << AV1_PROB_COST_SHIFT);
1601
1602
0
  double cost_none = RDCOST_DBL_WITH_NATIVE_BD_DIST(
1603
0
      x->rdmult, bits_none >> 4, rusi->sse[RESTORE_NONE],
1604
0
      rsc->cm->seq_params->bit_depth);
1605
0
  double cost_wiener = RDCOST_DBL_WITH_NATIVE_BD_DIST(
1606
0
      x->rdmult, bits_wiener >> 4, rusi->sse[RESTORE_WIENER],
1607
0
      rsc->cm->seq_params->bit_depth);
1608
1609
0
  RestorationType rtype =
1610
0
      (cost_wiener < cost_none) ? RESTORE_WIENER : RESTORE_NONE;
1611
0
  rusi->best_rtype[RESTORE_WIENER - 1] = rtype;
1612
1613
  // Set 'skip_sgr_eval' based on rdcost ratio of RESTORE_WIENER and
1614
  // RESTORE_NONE or based on best_rtype
1615
0
  if (rsc->lpf_sf->prune_sgr_based_on_wiener == 1) {
1616
0
    rusi->skip_sgr_eval = cost_wiener > (1.01 * cost_none);
1617
0
  } else if (rsc->lpf_sf->prune_sgr_based_on_wiener == 2) {
1618
0
    rusi->skip_sgr_eval = rusi->best_rtype[RESTORE_WIENER - 1] == RESTORE_NONE;
1619
0
  }
1620
1621
0
  rsc->sse += rusi->sse[rtype];
1622
0
  rsc->bits += (cost_wiener < cost_none) ? bits_wiener : bits_none;
1623
0
  if (cost_wiener < cost_none) rsc->wiener = rusi->wiener;
1624
0
}
1625
1626
static AOM_INLINE void search_norestore(const RestorationTileLimits *limits,
1627
                                        const AV1PixelRect *tile_rect,
1628
                                        int rest_unit_idx, void *priv,
1629
                                        int32_t *tmpbuf,
1630
0
                                        RestorationLineBuffers *rlbs) {
1631
0
  (void)tile_rect;
1632
0
  (void)tmpbuf;
1633
0
  (void)rlbs;
1634
1635
0
  RestSearchCtxt *rsc = (RestSearchCtxt *)priv;
1636
0
  RestUnitSearchInfo *rusi = &rsc->rusi[rest_unit_idx];
1637
1638
0
  const int highbd = rsc->cm->seq_params->use_highbitdepth;
1639
0
  rusi->sse[RESTORE_NONE] = sse_restoration_unit(
1640
0
      limits, rsc->src, &rsc->cm->cur_frame->buf, rsc->plane, highbd);
1641
1642
0
  rsc->sse += rusi->sse[RESTORE_NONE];
1643
0
}
1644
1645
static AOM_INLINE void search_switchable(const RestorationTileLimits *limits,
1646
                                         const AV1PixelRect *tile_rect,
1647
                                         int rest_unit_idx, void *priv,
1648
                                         int32_t *tmpbuf,
1649
0
                                         RestorationLineBuffers *rlbs) {
1650
0
  (void)limits;
1651
0
  (void)tile_rect;
1652
0
  (void)tmpbuf;
1653
0
  (void)rlbs;
1654
0
  RestSearchCtxt *rsc = (RestSearchCtxt *)priv;
1655
0
  RestUnitSearchInfo *rusi = &rsc->rusi[rest_unit_idx];
1656
1657
0
  const MACROBLOCK *const x = rsc->x;
1658
1659
0
  const int wiener_win =
1660
0
      (rsc->plane == AOM_PLANE_Y) ? WIENER_WIN : WIENER_WIN_CHROMA;
1661
1662
0
  double best_cost = 0;
1663
0
  int64_t best_bits = 0;
1664
0
  RestorationType best_rtype = RESTORE_NONE;
1665
1666
0
  for (RestorationType r = 0; r < RESTORE_SWITCHABLE_TYPES; ++r) {
1667
    // Check for the condition that wiener or sgrproj search could not
1668
    // find a solution or the solution was worse than RESTORE_NONE.
1669
    // In either case the best_rtype will be set as RESTORE_NONE. These
1670
    // should be skipped from the test below.
1671
0
    if (r > RESTORE_NONE) {
1672
0
      if (rusi->best_rtype[r - 1] == RESTORE_NONE) continue;
1673
0
    }
1674
1675
0
    const int64_t sse = rusi->sse[r];
1676
0
    int64_t coeff_pcost = 0;
1677
0
    switch (r) {
1678
0
      case RESTORE_NONE: coeff_pcost = 0; break;
1679
0
      case RESTORE_WIENER:
1680
0
        coeff_pcost =
1681
0
            count_wiener_bits(wiener_win, &rusi->wiener, &rsc->wiener);
1682
0
        break;
1683
0
      case RESTORE_SGRPROJ:
1684
0
        coeff_pcost = count_sgrproj_bits(&rusi->sgrproj, &rsc->sgrproj);
1685
0
        break;
1686
0
      default: assert(0); break;
1687
0
    }
1688
0
    const int64_t coeff_bits = coeff_pcost << AV1_PROB_COST_SHIFT;
1689
0
    const int64_t bits = x->mode_costs.switchable_restore_cost[r] + coeff_bits;
1690
0
    double cost = RDCOST_DBL_WITH_NATIVE_BD_DIST(
1691
0
        x->rdmult, bits >> 4, sse, rsc->cm->seq_params->bit_depth);
1692
0
    if (r == RESTORE_SGRPROJ && rusi->sgrproj.ep < 10)
1693
0
      cost *= (1 + DUAL_SGR_PENALTY_MULT * rsc->lpf_sf->dual_sgr_penalty_level);
1694
0
    if (r == 0 || cost < best_cost) {
1695
0
      best_cost = cost;
1696
0
      best_bits = bits;
1697
0
      best_rtype = r;
1698
0
    }
1699
0
  }
1700
1701
0
  rusi->best_rtype[RESTORE_SWITCHABLE - 1] = best_rtype;
1702
1703
0
  rsc->sse += rusi->sse[best_rtype];
1704
0
  rsc->bits += best_bits;
1705
0
  if (best_rtype == RESTORE_WIENER) rsc->wiener = rusi->wiener;
1706
0
  if (best_rtype == RESTORE_SGRPROJ) rsc->sgrproj = rusi->sgrproj;
1707
0
}
1708
1709
static AOM_INLINE void copy_unit_info(RestorationType frame_rtype,
1710
                                      const RestUnitSearchInfo *rusi,
1711
0
                                      RestorationUnitInfo *rui) {
1712
0
  assert(frame_rtype > 0);
1713
0
  rui->restoration_type = rusi->best_rtype[frame_rtype - 1];
1714
0
  if (rui->restoration_type == RESTORE_WIENER)
1715
0
    rui->wiener_info = rusi->wiener;
1716
0
  else
1717
0
    rui->sgrproj_info = rusi->sgrproj;
1718
0
}
1719
1720
0
static double search_rest_type(RestSearchCtxt *rsc, RestorationType rtype) {
1721
0
  static const rest_unit_visitor_t funs[RESTORE_TYPES] = {
1722
0
    search_norestore, search_wiener, search_sgrproj, search_switchable
1723
0
  };
1724
1725
0
  reset_rsc(rsc);
1726
0
  rsc_on_tile(rsc);
1727
1728
0
  av1_foreach_rest_unit_in_plane(rsc->cm, rsc->plane, funs[rtype], rsc,
1729
0
                                 &rsc->tile_rect, rsc->cm->rst_tmpbuf, NULL);
1730
0
  return RDCOST_DBL_WITH_NATIVE_BD_DIST(
1731
0
      rsc->x->rdmult, rsc->bits >> 4, rsc->sse, rsc->cm->seq_params->bit_depth);
1732
0
}
1733
1734
0
static int rest_tiles_in_plane(const AV1_COMMON *cm, int plane) {
1735
0
  const RestorationInfo *rsi = &cm->rst_info[plane];
1736
0
  return rsi->units_per_tile;
1737
0
}
1738
1739
0
void av1_pick_filter_restoration(const YV12_BUFFER_CONFIG *src, AV1_COMP *cpi) {
1740
0
  AV1_COMMON *const cm = &cpi->common;
1741
0
  MACROBLOCK *const x = &cpi->td.mb;
1742
0
  const SequenceHeader *const seq_params = cm->seq_params;
1743
0
  const int num_planes = av1_num_planes(cm);
1744
0
  assert(!cm->features.all_lossless);
1745
1746
0
  av1_fill_lr_rates(&x->mode_costs, x->e_mbd.tile_ctx);
1747
1748
0
  int ntiles[2];
1749
0
  for (int is_uv = 0; is_uv < 2; ++is_uv)
1750
0
    ntiles[is_uv] = rest_tiles_in_plane(cm, is_uv);
1751
1752
0
  assert(ntiles[1] <= ntiles[0]);
1753
0
  RestUnitSearchInfo *rusi =
1754
0
      (RestUnitSearchInfo *)aom_memalign(16, sizeof(*rusi) * ntiles[0]);
1755
1756
  // If the restoration unit dimensions are not multiples of
1757
  // rsi->restoration_unit_size then some elements of the rusi array may be
1758
  // left uninitialised when we reach copy_unit_info(...). This is not a
1759
  // problem, as these elements are ignored later, but in order to quiet
1760
  // Valgrind's warnings we initialise the array below.
1761
0
  memset(rusi, 0, sizeof(*rusi) * ntiles[0]);
1762
0
  x->rdmult = cpi->rd.RDMULT;
1763
1764
  // Allocate the frame buffer trial_frame_rst, which is used to temporarily
1765
  // store the loop restored frame.
1766
0
  if (aom_realloc_frame_buffer(
1767
0
          &cpi->trial_frame_rst, cm->superres_upscaled_width,
1768
0
          cm->superres_upscaled_height, seq_params->subsampling_x,
1769
0
          seq_params->subsampling_y, seq_params->use_highbitdepth,
1770
0
          AOM_RESTORATION_FRAME_BORDER, cm->features.byte_alignment, NULL, NULL,
1771
0
          NULL, 0))
1772
0
    aom_internal_error(cm->error, AOM_CODEC_MEM_ERROR,
1773
0
                       "Failed to allocate trial restored frame buffer");
1774
1775
0
  RestSearchCtxt rsc;
1776
0
  const int plane_start = AOM_PLANE_Y;
1777
0
  const int plane_end = num_planes > 1 ? AOM_PLANE_V : AOM_PLANE_Y;
1778
0
  for (int plane = plane_start; plane <= plane_end; ++plane) {
1779
0
    init_rsc(src, &cpi->common, x, &cpi->sf.lpf_sf, plane, rusi,
1780
0
             &cpi->trial_frame_rst, &rsc);
1781
1782
0
    const int plane_ntiles = ntiles[plane > 0];
1783
0
    const RestorationType num_rtypes =
1784
0
        (plane_ntiles > 1) ? RESTORE_TYPES : RESTORE_SWITCHABLE_TYPES;
1785
1786
0
    double best_cost = 0;
1787
0
    RestorationType best_rtype = RESTORE_NONE;
1788
1789
0
    const int highbd = rsc.cm->seq_params->use_highbitdepth;
1790
0
    if ((plane && !cpi->sf.lpf_sf.disable_loop_restoration_chroma) ||
1791
0
        (!plane && !cpi->sf.lpf_sf.disable_loop_restoration_luma)) {
1792
0
      av1_extend_frame(rsc.dgd_buffer, rsc.plane_width, rsc.plane_height,
1793
0
                       rsc.dgd_stride, RESTORATION_BORDER, RESTORATION_BORDER,
1794
0
                       highbd);
1795
1796
0
      for (RestorationType r = 0; r < num_rtypes; ++r) {
1797
0
        if ((force_restore_type != RESTORE_TYPES) && (r != RESTORE_NONE) &&
1798
0
            (r != force_restore_type))
1799
0
          continue;
1800
1801
0
        double cost = search_rest_type(&rsc, r);
1802
1803
0
        if (r == 0 || cost < best_cost) {
1804
0
          best_cost = cost;
1805
0
          best_rtype = r;
1806
0
        }
1807
0
      }
1808
0
    }
1809
1810
0
    cm->rst_info[plane].frame_restoration_type = best_rtype;
1811
0
    if (force_restore_type != RESTORE_TYPES)
1812
0
      assert(best_rtype == force_restore_type || best_rtype == RESTORE_NONE);
1813
1814
0
    if (best_rtype != RESTORE_NONE) {
1815
0
      for (int u = 0; u < plane_ntiles; ++u) {
1816
0
        copy_unit_info(best_rtype, &rusi[u], &cm->rst_info[plane].unit_info[u]);
1817
0
      }
1818
0
    }
1819
0
  }
1820
1821
0
  aom_free(rusi);
1822
0
}