pickrst.c 43.9 KB
Newer Older
1
/*
2
 * Copyright (c) 2016, Alliance for Open Media. All rights reserved
3
 *
4
5
6
7
8
9
 * 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
 */

#include <assert.h>
13
#include <float.h>
14
15
16
#include <limits.h>
#include <math.h>

Yaowu Xu's avatar
Yaowu Xu committed
17
#include "./aom_scale_rtcd.h"
18

19
#include "aom_dsp/psnr.h"
Yaowu Xu's avatar
Yaowu Xu committed
20
21
#include "aom_dsp/aom_dsp_common.h"
#include "aom_mem/aom_mem.h"
22
#include "aom_ports/mem.h"
23
#include "aom_ports/system_state.h"
24

25
26
#include "av1/common/onyxc_int.h"
#include "av1/common/quant_common.h"
27
#include "av1/common/restoration.h"
28

29
#include "av1/encoder/av1_quantize.h"
30
31
32
#include "av1/encoder/encoder.h"
#include "av1/encoder/picklpf.h"
#include "av1/encoder/pickrst.h"
33

34
35
36
// When set to RESTORE_WIENER or RESTORE_SGRPROJ only those are allowed.
// When set to RESTORE_NONE (0) we allow switchable.
const RestorationType force_restore_type = RESTORE_NONE;
37
38
39
40

// Number of Wiener iterations
#define NUM_WIENER_ITERS 10

41
typedef double (*search_restore_type)(const YV12_BUFFER_CONFIG *src,
42
43
                                      AV1_COMP *cpi, int partial_frame,
                                      RestorationInfo *info,
44
                                      RestorationType *rest_level,
45
46
                                      double *best_tile_cost,
                                      YV12_BUFFER_CONFIG *dst_frame);
47

48
const int frame_level_restore_bits[RESTORE_TYPES] = { 2, 2, 2, 2 };
49
50

static int64_t sse_restoration_tile(const YV12_BUFFER_CONFIG *src,
51
52
                                    const YV12_BUFFER_CONFIG *dst,
                                    const AV1_COMMON *cm, int h_start,
53
                                    int width, int v_start, int height,
54
55
                                    int components_pattern) {
  int64_t filt_err = 0;
56
57
58
59
  (void)cm;
  // Y and UV components cannot be mixed
  assert(components_pattern == 1 || components_pattern == 2 ||
         components_pattern == 4 || components_pattern == 6);
60
#if CONFIG_HIGHBITDEPTH
61
  if (cm->use_highbitdepth) {
62
63
64
65
66
    if ((components_pattern >> AOM_PLANE_Y) & 1) {
      filt_err +=
          aom_highbd_get_y_sse_part(src, dst, h_start, width, v_start, height);
    }
    if ((components_pattern >> AOM_PLANE_U) & 1) {
67
68
      filt_err +=
          aom_highbd_get_u_sse_part(src, dst, h_start, width, v_start, height);
69
70
    }
    if ((components_pattern >> AOM_PLANE_V) & 1) {
71
72
      filt_err +=
          aom_highbd_get_v_sse_part(src, dst, h_start, width, v_start, height);
73
74
    }
    return filt_err;
75
  }
76
#endif  // CONFIG_HIGHBITDEPTH
77
78
79
80
  if ((components_pattern >> AOM_PLANE_Y) & 1) {
    filt_err += aom_get_y_sse_part(src, dst, h_start, width, v_start, height);
  }
  if ((components_pattern >> AOM_PLANE_U) & 1) {
81
    filt_err += aom_get_u_sse_part(src, dst, h_start, width, v_start, height);
82
83
  }
  if ((components_pattern >> AOM_PLANE_V) & 1) {
84
    filt_err += aom_get_v_sse_part(src, dst, h_start, width, v_start, height);
85
  }
86
87
88
  return filt_err;
}

89
90
static int64_t sse_restoration_frame(AV1_COMMON *const cm,
                                     const YV12_BUFFER_CONFIG *src,
91
92
93
                                     const YV12_BUFFER_CONFIG *dst,
                                     int components_pattern) {
  int64_t filt_err = 0;
94
#if CONFIG_HIGHBITDEPTH
95
96
97
98
99
100
101
102
103
104
105
106
  if (cm->use_highbitdepth) {
    if ((components_pattern >> AOM_PLANE_Y) & 1) {
      filt_err += aom_highbd_get_y_sse(src, dst);
    }
    if ((components_pattern >> AOM_PLANE_U) & 1) {
      filt_err += aom_highbd_get_u_sse(src, dst);
    }
    if ((components_pattern >> AOM_PLANE_V) & 1) {
      filt_err += aom_highbd_get_v_sse(src, dst);
    }
    return filt_err;
  }
107
108
#else
  (void)cm;
109
#endif  // CONFIG_HIGHBITDEPTH
110
111
112
113
114
115
116
117
118
119
120
121
  if ((components_pattern >> AOM_PLANE_Y) & 1) {
    filt_err = aom_get_y_sse(src, dst);
  }
  if ((components_pattern >> AOM_PLANE_U) & 1) {
    filt_err += aom_get_u_sse(src, dst);
  }
  if ((components_pattern >> AOM_PLANE_V) & 1) {
    filt_err += aom_get_v_sse(src, dst);
  }
  return filt_err;
}

122
123
static int64_t try_restoration_tile(const YV12_BUFFER_CONFIG *src,
                                    AV1_COMP *const cpi, RestorationInfo *rsi,
124
125
126
                                    int components_pattern, int partial_frame,
                                    int tile_idx, int subtile_idx,
                                    int subtile_bits,
127
                                    YV12_BUFFER_CONFIG *dst_frame) {
128
129
130
131
  AV1_COMMON *const cm = &cpi->common;
  int64_t filt_err;
  int tile_width, tile_height, nhtiles, nvtiles;
  int h_start, h_end, v_start, v_end;
132
133
134
135
136
137
138
139
140
141
142
143
144
  int ntiles, width, height;

  // Y and UV components cannot be mixed
  assert(components_pattern == 1 || components_pattern == 2 ||
         components_pattern == 4 || components_pattern == 6);

  if (components_pattern == 1) {  // Y only
    width = src->y_crop_width;
    height = src->y_crop_height;
  } else {  // Color
    width = src->uv_crop_width;
    height = src->uv_crop_height;
  }
145
146
147
  ntiles = av1_get_rest_ntiles(
      width, height, cm->rst_info[components_pattern > 1].restoration_tilesize,
      &tile_width, &tile_height, &nhtiles, &nvtiles);
148
149
  (void)ntiles;

150
151
  av1_loop_restoration_frame(cm->frame_to_show, cm, rsi, components_pattern,
                             partial_frame, dst_frame);
152
  av1_get_rest_tile_limits(tile_idx, subtile_idx, subtile_bits, nhtiles,
153
154
                           nvtiles, tile_width, tile_height, width, height, 0,
                           0, &h_start, &h_end, &v_start, &v_end);
155
  filt_err = sse_restoration_tile(src, dst_frame, cm, h_start, h_end - h_start,
156
                                  v_start, v_end - v_start, components_pattern);
157
158
159
160
161

  return filt_err;
}

static int64_t try_restoration_frame(const YV12_BUFFER_CONFIG *src,
Yaowu Xu's avatar
Yaowu Xu committed
162
                                     AV1_COMP *const cpi, RestorationInfo *rsi,
163
                                     int components_pattern, int partial_frame,
164
                                     YV12_BUFFER_CONFIG *dst_frame) {
Yaowu Xu's avatar
Yaowu Xu committed
165
  AV1_COMMON *const cm = &cpi->common;
166
  int64_t filt_err;
167
168
  av1_loop_restoration_frame(cm->frame_to_show, cm, rsi, components_pattern,
                             partial_frame, dst_frame);
169
  filt_err = sse_restoration_frame(cm, src, dst_frame, components_pattern);
170
171
172
  return filt_err;
}

173
174
175
176
177
static int64_t get_pixel_proj_error(uint8_t *src8, int width, int height,
                                    int src_stride, uint8_t *dat8,
                                    int dat_stride, int bit_depth,
                                    int32_t *flt1, int flt1_stride,
                                    int32_t *flt2, int flt2_stride, int *xqd) {
178
179
180
181
  int i, j;
  int64_t err = 0;
  int xq[2];
  decode_xq(xqd, xq);
182
183
184
185
186
187
188
189
190
  if (bit_depth == 8) {
    const uint8_t *src = src8;
    const uint8_t *dat = dat8;
    for (i = 0; i < height; ++i) {
      for (j = 0; j < width; ++j) {
        const int32_t u =
            (int32_t)(dat[i * dat_stride + j] << SGRPROJ_RST_BITS);
        const int32_t f1 = (int32_t)flt1[i * flt1_stride + j] - u;
        const int32_t f2 = (int32_t)flt2[i * flt2_stride + j] - u;
David Barker's avatar
David Barker committed
191
        const int32_t v = xq[0] * f1 + xq[1] * f2 + (u << SGRPROJ_PRJ_BITS);
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
        const int32_t e =
            ROUND_POWER_OF_TWO(v, SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS) -
            src[i * src_stride + j];
        err += e * e;
      }
    }
  } else {
    const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
    const uint16_t *dat = CONVERT_TO_SHORTPTR(dat8);
    for (i = 0; i < height; ++i) {
      for (j = 0; j < width; ++j) {
        const int32_t u =
            (int32_t)(dat[i * dat_stride + j] << SGRPROJ_RST_BITS);
        const int32_t f1 = (int32_t)flt1[i * flt1_stride + j] - u;
        const int32_t f2 = (int32_t)flt2[i * flt2_stride + j] - u;
David Barker's avatar
David Barker committed
207
        const int32_t v = xq[0] * f1 + xq[1] * f2 + (u << SGRPROJ_PRJ_BITS);
208
209
210
211
212
        const int32_t e =
            ROUND_POWER_OF_TWO(v, SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS) -
            src[i * src_stride + j];
        err += e * e;
      }
213
214
215
216
217
    }
  }
  return err;
}

218
219
220
221
static void get_proj_subspace(uint8_t *src8, int width, int height,
                              int src_stride, uint8_t *dat8, int dat_stride,
                              int bit_depth, int32_t *flt1, int flt1_stride,
                              int32_t *flt2, int flt2_stride, int *xq) {
222
223
224
225
226
227
228
  int i, j;
  double H[2][2] = { { 0, 0 }, { 0, 0 } };
  double C[2] = { 0, 0 };
  double Det;
  double x[2];
  const int size = width * height;

229
230
  aom_clear_system_state();

231
232
233
  // Default
  xq[0] = 0;
  xq[1] = 0;
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
  if (bit_depth == 8) {
    const uint8_t *src = src8;
    const uint8_t *dat = dat8;
    for (i = 0; i < height; ++i) {
      for (j = 0; j < width; ++j) {
        const double u = (double)(dat[i * dat_stride + j] << SGRPROJ_RST_BITS);
        const double s =
            (double)(src[i * src_stride + j] << SGRPROJ_RST_BITS) - u;
        const double f1 = (double)flt1[i * flt1_stride + j] - u;
        const double f2 = (double)flt2[i * flt2_stride + j] - u;
        H[0][0] += f1 * f1;
        H[1][1] += f2 * f2;
        H[0][1] += f1 * f2;
        C[0] += f1 * s;
        C[1] += f2 * s;
      }
    }
  } else {
    const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
    const uint16_t *dat = CONVERT_TO_SHORTPTR(dat8);
    for (i = 0; i < height; ++i) {
      for (j = 0; j < width; ++j) {
        const double u = (double)(dat[i * dat_stride + j] << SGRPROJ_RST_BITS);
        const double s =
            (double)(src[i * src_stride + j] << SGRPROJ_RST_BITS) - u;
        const double f1 = (double)flt1[i * flt1_stride + j] - u;
        const double f2 = (double)flt2[i * flt2_stride + j] - u;
        H[0][0] += f1 * f1;
        H[1][1] += f2 * f2;
        H[0][1] += f1 * f2;
        C[0] += f1 * s;
        C[1] += f2 * s;
      }
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
    }
  }
  H[0][0] /= size;
  H[0][1] /= size;
  H[1][1] /= size;
  H[1][0] = H[0][1];
  C[0] /= size;
  C[1] /= size;
  Det = (H[0][0] * H[1][1] - H[0][1] * H[1][0]);
  if (Det < 1e-8) return;  // ill-posed, return default values
  x[0] = (H[1][1] * C[0] - H[0][1] * C[1]) / Det;
  x[1] = (H[0][0] * C[1] - H[1][0] * C[0]) / Det;
  xq[0] = (int)rint(x[0] * (1 << SGRPROJ_PRJ_BITS));
  xq[1] = (int)rint(x[1] * (1 << SGRPROJ_PRJ_BITS));
}

void encode_xq(int *xq, int *xqd) {
284
  xqd[0] = xq[0];
285
  xqd[0] = clamp(xqd[0], SGRPROJ_PRJ_MIN0, SGRPROJ_PRJ_MAX0);
286
  xqd[1] = (1 << SGRPROJ_PRJ_BITS) - xqd[0] - xq[1];
287
288
289
290
291
292
  xqd[1] = clamp(xqd[1], SGRPROJ_PRJ_MIN1, SGRPROJ_PRJ_MAX1);
}

static void search_selfguided_restoration(uint8_t *dat8, int width, int height,
                                          int dat_stride, uint8_t *src8,
                                          int src_stride, int bit_depth,
293
294
                                          int *eps, int *xqd, int32_t *rstbuf) {
  int32_t *flt1 = rstbuf;
295
  int32_t *flt2 = flt1 + RESTORATION_TILEPELS_MAX;
296
  int32_t *tmpbuf2 = flt2 + RESTORATION_TILEPELS_MAX;
297
  int ep, bestep = 0;
298
299
  int64_t err, besterr = -1;
  int exqd[2], bestxqd[2] = { 0, 0 };
300

301
302
  for (ep = 0; ep < SGRPROJ_PARAMS; ep++) {
    int exq[2];
303
#if CONFIG_HIGHBITDEPTH
304
305
    if (bit_depth > 8) {
      uint16_t *dat = CONVERT_TO_SHORTPTR(dat8);
306
307
#if USE_HIGHPASS_IN_SGRPROJ
      av1_highpass_filter_highbd(dat, width, height, dat_stride, flt1, width,
308
                                 sgr_params[ep].corner, sgr_params[ep].edge);
309
#else
310
311
312
      av1_selfguided_restoration_highbd(dat, width, height, dat_stride, flt1,
                                        width, bit_depth, sgr_params[ep].r1,
                                        sgr_params[ep].e1, tmpbuf2);
313
#endif  // USE_HIGHPASS_IN_SGRPROJ
314
315
316
      av1_selfguided_restoration_highbd(dat, width, height, dat_stride, flt2,
                                        width, bit_depth, sgr_params[ep].r2,
                                        sgr_params[ep].e2, tmpbuf2);
317
    } else {
318
#endif
319
320
#if USE_HIGHPASS_IN_SGRPROJ
      av1_highpass_filter(dat8, width, height, dat_stride, flt1, width,
321
                          sgr_params[ep].corner, sgr_params[ep].edge);
322
323
324
325
#else
    av1_selfguided_restoration(dat8, width, height, dat_stride, flt1, width,
                               sgr_params[ep].r1, sgr_params[ep].e1, tmpbuf2);
#endif  // USE_HIGHPASS_IN_SGRPROJ
326
      av1_selfguided_restoration(dat8, width, height, dat_stride, flt2, width,
327
                                 sgr_params[ep].r2, sgr_params[ep].e2, tmpbuf2);
328
#if CONFIG_HIGHBITDEPTH
329
    }
330
#endif
331
332
    get_proj_subspace(src8, width, height, src_stride, dat8, dat_stride,
                      bit_depth, flt1, width, flt2, width, exq);
333
    encode_xq(exq, exqd);
334
335
336
    err =
        get_pixel_proj_error(src8, width, height, src_stride, dat8, dat_stride,
                             bit_depth, flt1, width, flt2, width, exqd);
337
338
339
340
341
342
343
344
345
346
347
348
349
    if (besterr == -1 || err < besterr) {
      bestep = ep;
      besterr = err;
      bestxqd[0] = exqd[0];
      bestxqd[1] = exqd[1];
    }
  }
  *eps = bestep;
  xqd[0] = bestxqd[0];
  xqd[1] = bestxqd[1];
}

static double search_sgrproj(const YV12_BUFFER_CONFIG *src, AV1_COMP *cpi,
350
351
                             int partial_frame, RestorationInfo *info,
                             RestorationType *type, double *best_tile_cost,
352
                             YV12_BUFFER_CONFIG *dst_frame) {
353
354
355
356
357
358
  SgrprojInfo *sgrproj_info = info->sgrproj_info;
  double err, cost_norestore, cost_sgrproj;
  int bits;
  MACROBLOCK *x = &cpi->td.mb;
  AV1_COMMON *const cm = &cpi->common;
  const YV12_BUFFER_CONFIG *dgd = cm->frame_to_show;
359
  RestorationInfo *rsi = &cpi->rst_search[0];
360
361
  int tile_idx, tile_width, tile_height, nhtiles, nvtiles;
  int h_start, h_end, v_start, v_end;
362
  // Allocate for the src buffer at high precision
363
364
365
  const int ntiles = av1_get_rest_ntiles(
      cm->width, cm->height, cm->rst_info[0].restoration_tilesize, &tile_width,
      &tile_height, &nhtiles, &nvtiles);
366
  rsi->frame_restoration_type = RESTORE_SGRPROJ;
367

368
369
370
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
    rsi->restoration_type[tile_idx] = RESTORE_NONE;
  }
371
372
373
374
375
  // Compute best Sgrproj filters for each tile
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
    av1_get_rest_tile_limits(tile_idx, 0, 0, nhtiles, nvtiles, tile_width,
                             tile_height, cm->width, cm->height, 0, 0, &h_start,
                             &h_end, &v_start, &v_end);
376
    err = sse_restoration_tile(src, cm->frame_to_show, cm, h_start,
377
                               h_end - h_start, v_start, v_end - v_start, 1);
378
379
380
381
382
383
384
385
    // #bits when a tile is not restored
    bits = av1_cost_bit(RESTORE_NONE_SGRPROJ_PROB, 0);
    cost_norestore = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);
    best_tile_cost[tile_idx] = DBL_MAX;
    search_selfguided_restoration(
        dgd->y_buffer + v_start * dgd->y_stride + h_start, h_end - h_start,
        v_end - v_start, dgd->y_stride,
        src->y_buffer + v_start * src->y_stride + h_start, src->y_stride,
386
#if CONFIG_HIGHBITDEPTH
387
388
389
        cm->bit_depth,
#else
        8,
390
#endif  // CONFIG_HIGHBITDEPTH
391
        &rsi->sgrproj_info[tile_idx].ep, rsi->sgrproj_info[tile_idx].xqd,
392
        cm->rst_internal.tmpbuf);
393
    rsi->restoration_type[tile_idx] = RESTORE_SGRPROJ;
394
    err = try_restoration_tile(src, cpi, rsi, 1, partial_frame, tile_idx, 0, 0,
395
                               dst_frame);
396
397
398
399
    bits = SGRPROJ_BITS << AV1_PROB_COST_SHIFT;
    bits += av1_cost_bit(RESTORE_NONE_SGRPROJ_PROB, 1);
    cost_sgrproj = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);
    if (cost_sgrproj >= cost_norestore) {
400
      type[tile_idx] = RESTORE_NONE;
401
    } else {
402
      type[tile_idx] = RESTORE_SGRPROJ;
403
      memcpy(&sgrproj_info[tile_idx], &rsi->sgrproj_info[tile_idx],
404
405
406
407
408
409
             sizeof(sgrproj_info[tile_idx]));
      bits = SGRPROJ_BITS << AV1_PROB_COST_SHIFT;
      best_tile_cost[tile_idx] = RDCOST_DBL(
          x->rdmult, x->rddiv,
          (bits + cpi->switchable_restore_cost[RESTORE_SGRPROJ]) >> 4, err);
    }
410
    rsi->restoration_type[tile_idx] = RESTORE_NONE;
411
412
  }
  // Cost for Sgrproj filtering
413
  bits = frame_level_restore_bits[rsi->frame_restoration_type]
414
415
416
         << AV1_PROB_COST_SHIFT;
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
    bits +=
417
        av1_cost_bit(RESTORE_NONE_SGRPROJ_PROB, type[tile_idx] != RESTORE_NONE);
418
    memcpy(&rsi->sgrproj_info[tile_idx], &sgrproj_info[tile_idx],
419
           sizeof(sgrproj_info[tile_idx]));
420
    if (type[tile_idx] == RESTORE_SGRPROJ) {
421
422
      bits += (SGRPROJ_BITS << AV1_PROB_COST_SHIFT);
    }
423
    rsi->restoration_type[tile_idx] = type[tile_idx];
424
  }
425
  err = try_restoration_frame(src, cpi, rsi, 1, partial_frame, dst_frame);
426
427
428
429
430
  cost_sgrproj = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);

  return cost_sgrproj;
}

431
432
static double find_average(uint8_t *src, int h_start, int h_end, int v_start,
                           int v_end, int stride) {
433
434
435
  uint64_t sum = 0;
  double avg = 0;
  int i, j;
436
  aom_clear_system_state();
437
438
439
  for (i = v_start; i < v_end; i++)
    for (j = h_start; j < h_end; j++) sum += src[i * stride + j];
  avg = (double)sum / ((v_end - v_start) * (h_end - h_start));
440
441
442
  return avg;
}

443
444
445
static void compute_stats(uint8_t *dgd, uint8_t *src, int h_start, int h_end,
                          int v_start, int v_end, int dgd_stride,
                          int src_stride, double *M, double *H) {
446
  int i, j, k, l;
447
  double Y[WIENER_WIN2];
448
449
  const double avg =
      find_average(dgd, h_start, h_end, v_start, v_end, dgd_stride);
450

451
452
  memset(M, 0, sizeof(*M) * WIENER_WIN2);
  memset(H, 0, sizeof(*H) * WIENER_WIN2 * WIENER_WIN2);
453
454
  for (i = v_start; i < v_end; i++) {
    for (j = h_start; j < h_end; j++) {
455
456
      const double X = (double)src[i * src_stride + j] - avg;
      int idx = 0;
457
458
      for (k = -WIENER_HALFWIN; k <= WIENER_HALFWIN; k++) {
        for (l = -WIENER_HALFWIN; l <= WIENER_HALFWIN; l++) {
459
460
461
462
          Y[idx] = (double)dgd[(i + l) * dgd_stride + (j + k)] - avg;
          idx++;
        }
      }
463
      for (k = 0; k < WIENER_WIN2; ++k) {
464
        M[k] += Y[k] * X;
465
466
        H[k * WIENER_WIN2 + k] += Y[k] * Y[k];
        for (l = k + 1; l < WIENER_WIN2; ++l) {
467
468
469
470
          // H is a symmetric matrix, so we only need to fill out the upper
          // triangle here. We can copy it down to the lower triangle outside
          // the (i, j) loops.
          H[k * WIENER_WIN2 + l] += Y[k] * Y[l];
471
472
473
474
        }
      }
    }
  }
475
476
477
478
479
  for (k = 0; k < WIENER_WIN2; ++k) {
    for (l = k + 1; l < WIENER_WIN2; ++l) {
      H[l * WIENER_WIN2 + k] = H[k * WIENER_WIN2 + l];
    }
  }
480
481
}

482
#if CONFIG_HIGHBITDEPTH
483
484
static double find_average_highbd(uint16_t *src, int h_start, int h_end,
                                  int v_start, int v_end, int stride) {
485
486
487
  uint64_t sum = 0;
  double avg = 0;
  int i, j;
488
  aom_clear_system_state();
489
490
491
  for (i = v_start; i < v_end; i++)
    for (j = h_start; j < h_end; j++) sum += src[i * stride + j];
  avg = (double)sum / ((v_end - v_start) * (h_end - h_start));
492
493
494
  return avg;
}

495
496
497
498
static void compute_stats_highbd(uint8_t *dgd8, uint8_t *src8, int h_start,
                                 int h_end, int v_start, int v_end,
                                 int dgd_stride, int src_stride, double *M,
                                 double *H) {
499
  int i, j, k, l;
500
  double Y[WIENER_WIN2];
501
502
  uint16_t *src = CONVERT_TO_SHORTPTR(src8);
  uint16_t *dgd = CONVERT_TO_SHORTPTR(dgd8);
503
504
  const double avg =
      find_average_highbd(dgd, h_start, h_end, v_start, v_end, dgd_stride);
505

506
507
  memset(M, 0, sizeof(*M) * WIENER_WIN2);
  memset(H, 0, sizeof(*H) * WIENER_WIN2 * WIENER_WIN2);
508
509
  for (i = v_start; i < v_end; i++) {
    for (j = h_start; j < h_end; j++) {
510
511
      const double X = (double)src[i * src_stride + j] - avg;
      int idx = 0;
512
513
      for (k = -WIENER_HALFWIN; k <= WIENER_HALFWIN; k++) {
        for (l = -WIENER_HALFWIN; l <= WIENER_HALFWIN; l++) {
514
515
516
517
          Y[idx] = (double)dgd[(i + l) * dgd_stride + (j + k)] - avg;
          idx++;
        }
      }
518
      for (k = 0; k < WIENER_WIN2; ++k) {
519
        M[k] += Y[k] * X;
520
521
        H[k * WIENER_WIN2 + k] += Y[k] * Y[k];
        for (l = k + 1; l < WIENER_WIN2; ++l) {
522
523
524
525
          // H is a symmetric matrix, so we only need to fill out the upper
          // triangle here. We can copy it down to the lower triangle outside
          // the (i, j) loops.
          H[k * WIENER_WIN2 + l] += Y[k] * Y[l];
526
527
528
529
        }
      }
    }
  }
530
531
532
533
534
  for (k = 0; k < WIENER_WIN2; ++k) {
    for (l = k + 1; l < WIENER_WIN2; ++l) {
      H[l * WIENER_WIN2 + k] = H[k * WIENER_WIN2 + l];
    }
  }
535
}
536
#endif  // CONFIG_HIGHBITDEPTH
537
538
539
540
541

// Solves Ax = b, where x and b are column vectors
static int linsolve(int n, double *A, int stride, double *b, double *x) {
  int i, j, k;
  double c;
542
543
544

  aom_clear_system_state();

545
546
  // Forward elimination
  for (k = 0; k < n - 1; k++) {
547
548
549
550
551
552
553
554
555
556
557
558
559
    // Bring the largest magitude to the diagonal position
    for (i = n - 1; i > k; i--) {
      if (fabs(A[(i - 1) * stride + k]) < fabs(A[i * stride + k])) {
        for (j = 0; j < n; j++) {
          c = A[i * stride + j];
          A[i * stride + j] = A[(i - 1) * stride + j];
          A[(i - 1) * stride + j] = c;
        }
        c = b[i];
        b[i] = b[i - 1];
        b[i - 1] = c;
      }
    }
560
    for (i = k; i < n - 1; i++) {
561
      if (fabs(A[k * stride + k]) < 1e-10) return 0;
562
      c = A[(i + 1) * stride + k] / A[k * stride + k];
563
      for (j = 0; j < n; j++) A[(i + 1) * stride + j] -= c * A[k * stride + j];
564
565
566
567
568
      b[i + 1] -= c * b[k];
    }
  }
  // Backward substitution
  for (i = n - 1; i >= 0; i--) {
569
    if (fabs(A[i * stride + i]) < 1e-10) return 0;
570
    c = 0;
571
    for (j = i + 1; j <= n - 1; j++) c += A[i * stride + j] * x[j];
572
573
    x[i] = (b[i] - c) / A[i * stride + i];
  }
574

575
576
577
578
  return 1;
}

static INLINE int wrap_index(int i) {
579
  return (i >= WIENER_HALFWIN1 ? WIENER_WIN - 1 - i : i);
580
581
582
583
584
}

// Fix vector b, update vector a
static void update_a_sep_sym(double **Mc, double **Hc, double *a, double *b) {
  int i, j;
585
  double S[WIENER_WIN];
586
  double A[WIENER_HALFWIN1], B[WIENER_HALFWIN1 * WIENER_HALFWIN1];
Aamir Anis's avatar
Aamir Anis committed
587
  int w, w2;
588
589
  memset(A, 0, sizeof(A));
  memset(B, 0, sizeof(B));
590
591
  for (i = 0; i < WIENER_WIN; i++) {
    for (j = 0; j < WIENER_WIN; ++j) {
592
593
594
595
      const int jj = wrap_index(j);
      A[jj] += Mc[i][j] * b[i];
    }
  }
596
597
  for (i = 0; i < WIENER_WIN; i++) {
    for (j = 0; j < WIENER_WIN; j++) {
598
      int k, l;
599
600
      for (k = 0; k < WIENER_WIN; ++k)
        for (l = 0; l < WIENER_WIN; ++l) {
601
602
          const int kk = wrap_index(k);
          const int ll = wrap_index(l);
603
604
          B[ll * WIENER_HALFWIN1 + kk] +=
              Hc[j * WIENER_WIN + i][k * WIENER_WIN2 + l] * b[i] * b[j];
605
606
607
        }
    }
  }
Aamir Anis's avatar
Aamir Anis committed
608
  // Normalization enforcement in the system of equations itself
609
  w = WIENER_WIN;
Aamir Anis's avatar
Aamir Anis committed
610
611
  w2 = (w >> 1) + 1;
  for (i = 0; i < w2 - 1; ++i)
612
613
    A[i] -=
        A[w2 - 1] * 2 + B[i * w2 + w2 - 1] - 2 * B[(w2 - 1) * w2 + (w2 - 1)];
Aamir Anis's avatar
Aamir Anis committed
614
615
616
617
618
619
620
621
622
  for (i = 0; i < w2 - 1; ++i)
    for (j = 0; j < w2 - 1; ++j)
      B[i * w2 + j] -= 2 * (B[i * w2 + (w2 - 1)] + B[(w2 - 1) * w2 + j] -
                            2 * B[(w2 - 1) * w2 + (w2 - 1)]);
  if (linsolve(w2 - 1, B, w2, A, S)) {
    S[w2 - 1] = 1.0;
    for (i = w2; i < w; ++i) {
      S[i] = S[w - 1 - i];
      S[w2 - 1] -= 2 * S[i];
623
    }
Aamir Anis's avatar
Aamir Anis committed
624
    memcpy(a, S, w * sizeof(*a));
625
626
627
628
629
630
  }
}

// Fix vector a, update vector b
static void update_b_sep_sym(double **Mc, double **Hc, double *a, double *b) {
  int i, j;
631
  double S[WIENER_WIN];
632
  double A[WIENER_HALFWIN1], B[WIENER_HALFWIN1 * WIENER_HALFWIN1];
Aamir Anis's avatar
Aamir Anis committed
633
  int w, w2;
634
635
  memset(A, 0, sizeof(A));
  memset(B, 0, sizeof(B));
636
  for (i = 0; i < WIENER_WIN; i++) {
637
    const int ii = wrap_index(i);
638
    for (j = 0; j < WIENER_WIN; j++) A[ii] += Mc[i][j] * a[j];
639
640
  }

641
642
  for (i = 0; i < WIENER_WIN; i++) {
    for (j = 0; j < WIENER_WIN; j++) {
643
644
645
      const int ii = wrap_index(i);
      const int jj = wrap_index(j);
      int k, l;
646
647
648
649
      for (k = 0; k < WIENER_WIN; ++k)
        for (l = 0; l < WIENER_WIN; ++l)
          B[jj * WIENER_HALFWIN1 + ii] +=
              Hc[i * WIENER_WIN + j][k * WIENER_WIN2 + l] * a[k] * a[l];
650
651
    }
  }
Aamir Anis's avatar
Aamir Anis committed
652
  // Normalization enforcement in the system of equations itself
653
654
  w = WIENER_WIN;
  w2 = WIENER_HALFWIN1;
Aamir Anis's avatar
Aamir Anis committed
655
  for (i = 0; i < w2 - 1; ++i)
656
657
    A[i] -=
        A[w2 - 1] * 2 + B[i * w2 + w2 - 1] - 2 * B[(w2 - 1) * w2 + (w2 - 1)];
Aamir Anis's avatar
Aamir Anis committed
658
659
660
661
662
663
664
665
666
  for (i = 0; i < w2 - 1; ++i)
    for (j = 0; j < w2 - 1; ++j)
      B[i * w2 + j] -= 2 * (B[i * w2 + (w2 - 1)] + B[(w2 - 1) * w2 + j] -
                            2 * B[(w2 - 1) * w2 + (w2 - 1)]);
  if (linsolve(w2 - 1, B, w2, A, S)) {
    S[w2 - 1] = 1.0;
    for (i = w2; i < w; ++i) {
      S[i] = S[w - 1 - i];
      S[w2 - 1] -= 2 * S[i];
667
    }
Aamir Anis's avatar
Aamir Anis committed
668
    memcpy(b, S, w * sizeof(*b));
669
670
671
  }
}

672
673
static int wiener_decompose_sep_sym(double *M, double *H, double *a,
                                    double *b) {
674
  static const double init_filt[WIENER_WIN] = {
675
    0.035623, -0.127154, 0.211436, 0.760190, 0.211436, -0.127154, 0.035623,
676
677
  };
  int i, j, iter;
678
679
680
681
682
683
684
  double *Hc[WIENER_WIN2];
  double *Mc[WIENER_WIN];
  for (i = 0; i < WIENER_WIN; i++) {
    Mc[i] = M + i * WIENER_WIN;
    for (j = 0; j < WIENER_WIN; j++) {
      Hc[i * WIENER_WIN + j] =
          H + i * WIENER_WIN * WIENER_WIN2 + j * WIENER_WIN;
685
686
    }
  }
687
688
  memcpy(a, init_filt, sizeof(*a) * WIENER_WIN);
  memcpy(b, init_filt, sizeof(*b) * WIENER_WIN);
689
690

  iter = 1;
691
  while (iter < NUM_WIENER_ITERS) {
692
693
694
695
    update_a_sep_sym(Mc, Hc, a, b);
    update_b_sep_sym(Mc, Hc, a, b);
    iter++;
  }
696
  return 1;
697
698
}

699
// Computes the function x'*H*x - x'*M for the learned 2D filter x, and compares
Aamir Anis's avatar
Aamir Anis committed
700
701
// against identity filters; Final score is defined as the difference between
// the function values
702
703
static double compute_score(double *M, double *H, InterpKernel vfilt,
                            InterpKernel hfilt) {
704
  double ab[WIENER_WIN * WIENER_WIN];
Aamir Anis's avatar
Aamir Anis committed
705
706
707
708
  int i, k, l;
  double P = 0, Q = 0;
  double iP = 0, iQ = 0;
  double Score, iScore;
709
  double a[WIENER_WIN], b[WIENER_WIN];
710
711
712

  aom_clear_system_state();

713
714
715
716
717
718
  a[WIENER_HALFWIN] = b[WIENER_HALFWIN] = 1.0;
  for (i = 0; i < WIENER_HALFWIN; ++i) {
    a[i] = a[WIENER_WIN - i - 1] = (double)vfilt[i] / WIENER_FILT_STEP;
    b[i] = b[WIENER_WIN - i - 1] = (double)hfilt[i] / WIENER_FILT_STEP;
    a[WIENER_HALFWIN] -= 2 * a[i];
    b[WIENER_HALFWIN] -= 2 * b[i];
Aamir Anis's avatar
Aamir Anis committed
719
  }
720
721
  for (k = 0; k < WIENER_WIN; ++k) {
    for (l = 0; l < WIENER_WIN; ++l) ab[k * WIENER_WIN + l] = a[l] * b[k];
Aamir Anis's avatar
Aamir Anis committed
722
  }
723
  for (k = 0; k < WIENER_WIN2; ++k) {
Aamir Anis's avatar
Aamir Anis committed
724
    P += ab[k] * M[k];
725
726
    for (l = 0; l < WIENER_WIN2; ++l)
      Q += ab[k] * H[k * WIENER_WIN2 + l] * ab[l];
Aamir Anis's avatar
Aamir Anis committed
727
728
729
  }
  Score = Q - 2 * P;

730
731
  iP = M[WIENER_WIN2 >> 1];
  iQ = H[(WIENER_WIN2 >> 1) * WIENER_WIN2 + (WIENER_WIN2 >> 1)];
Aamir Anis's avatar
Aamir Anis committed
732
733
734
735
736
  iScore = iQ - 2 * iP;

  return Score - iScore;
}

737
static void quantize_sym_filter(double *f, InterpKernel fi) {
738
  int i;
739
740
  for (i = 0; i < WIENER_HALFWIN; ++i) {
    fi[i] = RINT(f[i] * WIENER_FILT_STEP);
741
742
743
744
745
  }
  // Specialize for 7-tap filter
  fi[0] = CLIP(fi[0], WIENER_FILT_TAP0_MINV, WIENER_FILT_TAP0_MAXV);
  fi[1] = CLIP(fi[1], WIENER_FILT_TAP1_MINV, WIENER_FILT_TAP1_MAXV);
  fi[2] = CLIP(fi[2], WIENER_FILT_TAP2_MINV, WIENER_FILT_TAP2_MAXV);
746
747
748
749
  // Satisfy filter constraints
  fi[WIENER_WIN - 1] = fi[0];
  fi[WIENER_WIN - 2] = fi[1];
  fi[WIENER_WIN - 3] = fi[2];
750
751
  // The central element has an implicit +WIENER_FILT_STEP
  fi[3] = -2 * (fi[0] + fi[1] + fi[2]);
752
753
754
}

static double search_wiener_uv(const YV12_BUFFER_CONFIG *src, AV1_COMP *cpi,
755
                               int partial_frame, int plane,
756
                               RestorationInfo *info, RestorationType *type,
757
758
759
760
761
762
                               YV12_BUFFER_CONFIG *dst_frame) {
  WienerInfo *wiener_info = info->wiener_info;
  AV1_COMMON *const cm = &cpi->common;
  RestorationInfo *rsi = cpi->rst_search;
  int64_t err;
  int bits;
763
  double cost_wiener, cost_norestore, cost_wiener_frame, cost_norestore_frame;
764
765
766
767
768
769
770
771
772
773
774
  MACROBLOCK *x = &cpi->td.mb;
  double M[WIENER_WIN2];
  double H[WIENER_WIN2 * WIENER_WIN2];
  double vfilterd[WIENER_WIN], hfilterd[WIENER_WIN];
  const YV12_BUFFER_CONFIG *dgd = cm->frame_to_show;
  const int width = src->uv_crop_width;
  const int height = src->uv_crop_height;
  const int src_stride = src->uv_stride;
  const int dgd_stride = dgd->uv_stride;
  double score;
  int tile_idx, tile_width, tile_height, nhtiles, nvtiles;
775
  int h_start, h_end, v_start, v_end;
776
777
778
  const int ntiles =
      av1_get_rest_ntiles(width, height, cm->rst_info[1].restoration_tilesize,
                          &tile_width, &tile_height, &nhtiles, &nvtiles);
779
780
781
782
  assert(width == dgd->uv_crop_width);
  assert(height == dgd->uv_crop_height);

  rsi[plane].frame_restoration_type = RESTORE_NONE;
783
  err = sse_restoration_frame(cm, src, cm->frame_to_show, (1 << plane));
784
  bits = 0;
785
  cost_norestore_frame = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);
786
787

  rsi[plane].frame_restoration_type = RESTORE_WIENER;
788

789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
    rsi[plane].restoration_type[tile_idx] = RESTORE_NONE;
  }

  // Compute best Wiener filters for each tile
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
    av1_get_rest_tile_limits(tile_idx, 0, 0, nhtiles, nvtiles, tile_width,
                             tile_height, width, height, 0, 0, &h_start, &h_end,
                             &v_start, &v_end);
    err = sse_restoration_tile(src, cm->frame_to_show, cm, h_start,
                               h_end - h_start, v_start, v_end - v_start,
                               1 << plane);
    // #bits when a tile is not restored
    bits = av1_cost_bit(RESTORE_NONE_WIENER_PROB, 0);
    cost_norestore = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);
    // best_tile_cost[tile_idx] = DBL_MAX;

    av1_get_rest_tile_limits(tile_idx, 0, 0, nhtiles, nvtiles, tile_width,
                             tile_height, width, height, WIENER_HALFWIN,
                             WIENER_HALFWIN, &h_start, &h_end, &v_start,
                             &v_end);
810
    if (plane == AOM_PLANE_U) {
811
#if CONFIG_HIGHBITDEPTH
812
813
814
815
      if (cm->use_highbitdepth)
        compute_stats_highbd(dgd->u_buffer, src->u_buffer, h_start, h_end,
                             v_start, v_end, dgd_stride, src_stride, M, H);
      else
816
#endif  // CONFIG_HIGHBITDEPTH
817
818
        compute_stats(dgd->u_buffer, src->u_buffer, h_start, h_end, v_start,
                      v_end, dgd_stride, src_stride, M, H);
819
    } else if (plane == AOM_PLANE_V) {
820
#if CONFIG_HIGHBITDEPTH
821
822
823
824
      if (cm->use_highbitdepth)
        compute_stats_highbd(dgd->v_buffer, src->v_buffer, h_start, h_end,
                             v_start, v_end, dgd_stride, src_stride, M, H);
      else
825
#endif  // CONFIG_HIGHBITDEPTH
826
827
        compute_stats(dgd->v_buffer, src->v_buffer, h_start, h_end, v_start,
                      v_end, dgd_stride, src_stride, M, H);
828
829
830
    } else {
      assert(0);
    }
831
832
833
834
835
836

    type[tile_idx] = RESTORE_WIENER;

    if (!wiener_decompose_sep_sym(M, H, vfilterd, hfilterd)) {
      type[tile_idx] = RESTORE_NONE;
      continue;
837
    }
838
839
    quantize_sym_filter(vfilterd, rsi[plane].wiener_info[tile_idx].vfilter);
    quantize_sym_filter(hfilterd, rsi[plane].wiener_info[tile_idx].hfilter);
840

841
842
843
844
845
846
847
848
849
    // Filter score computes the value of the function x'*A*x - x'*b for the
    // learned filter and compares it against identity filer. If there is no
    // reduction in the function, the filter is reverted back to identity
    score = compute_score(M, H, rsi[plane].wiener_info[tile_idx].vfilter,
                          rsi[plane].wiener_info[tile_idx].hfilter);
    if (score > 0.0) {
      type[tile_idx] = RESTORE_NONE;
      continue;
    }
850

851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
    rsi[plane].restoration_type[tile_idx] = RESTORE_WIENER;
    err = try_restoration_tile(src, cpi, rsi, 1 << plane, partial_frame,
                               tile_idx, 0, 0, dst_frame);
    bits = WIENER_FILT_BITS << AV1_PROB_COST_SHIFT;
    bits += av1_cost_bit(RESTORE_NONE_WIENER_PROB, 1);
    cost_wiener = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);
    if (cost_wiener >= cost_norestore) {
      type[tile_idx] = RESTORE_NONE;
    } else {
      type[tile_idx] = RESTORE_WIENER;
      memcpy(&wiener_info[tile_idx], &rsi[plane].wiener_info[tile_idx],
             sizeof(wiener_info[tile_idx]));
    }
    rsi[plane].restoration_type[tile_idx] = RESTORE_NONE;
  }
  // Cost for Wiener filtering
  bits = 0;
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
    bits +=
        av1_cost_bit(RESTORE_NONE_WIENER_PROB, type[tile_idx] != RESTORE_NONE);
    memcpy(&rsi[plane].wiener_info[tile_idx], &wiener_info[tile_idx],
           sizeof(wiener_info[tile_idx]));
    if (type[tile_idx] == RESTORE_WIENER) {
      bits += (WIENER_FILT_BITS << AV1_PROB_COST_SHIFT);
    }
    rsi[plane].restoration_type[tile_idx] = type[tile_idx];
877
  }
878
  err = try_restoration_frame(src, cpi, rsi, 1 << plane, partial_frame,
879
                              dst_frame);
880
881
882
883
884
  cost_wiener_frame = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);

  if (cost_wiener_frame < cost_norestore_frame) {
    info->frame_restoration_type = RESTORE_WIENER;
  } else {
885
886
887
    info->frame_restoration_type = RESTORE_NONE;
  }

888
889
  return info->frame_restoration_type == RESTORE_WIENER ? cost_wiener_frame
                                                        : cost_norestore_frame;
890
891
}

892
static double search_wiener(const YV12_BUFFER_CONFIG *src, AV1_COMP *cpi,
893
894
                            int partial_frame, RestorationInfo *info,
                            RestorationType *type, double *best_tile_cost,
895
                            YV12_BUFFER_CONFIG *dst_frame) {
896
  WienerInfo *wiener_info = info->wiener_info;
Yaowu Xu's avatar
Yaowu Xu committed
897
  AV1_COMMON *const cm = &cpi->common;
898
  RestorationInfo *rsi = cpi->rst_search;
899
900
  int64_t err;
  int bits;
901
  double cost_wiener, cost_norestore;
902
  MACROBLOCK *x = &cpi->td.mb;
903
904
905
  double M[WIENER_WIN2];
  double H[WIENER_WIN2 * WIENER_WIN2];
  double vfilterd[WIENER_WIN], hfilterd[WIENER_WIN];
906
907
908
909
910
  const YV12_BUFFER_CONFIG *dgd = cm->frame_to_show;
  const int width = cm->width;
  const int height = cm->height;
  const int src_stride = src->y_stride;
  const int dgd_stride = dgd->y_stride;
Aamir Anis's avatar
Aamir Anis committed
911
  double score;
912
  int tile_idx, tile_width, tile_height, nhtiles, nvtiles;
913
  int h_start, h_end, v_start, v_end;
914
915
916
  const int ntiles =
      av1_get_rest_ntiles(width, height, cm->rst_info[0].restoration_tilesize,
                          &tile_width, &tile_height, &nhtiles, &nvtiles);
917
918
919
920
921
  assert(width == dgd->y_crop_width);
  assert(height == dgd->y_crop_height);
  assert(width == src->y_crop_width);
  assert(height == src->y_crop_height);

922
  rsi->frame_restoration_type = RESTORE_WIENER;
923

924
925
926
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
    rsi->restoration_type[tile_idx] = RESTORE_NONE;
  }
927

928
// Construct a (WIENER_HALFWIN)-pixel border around the frame
929
#if CONFIG_HIGHBITDEPTH
930
931
932
933
934
935
936
  if (cm->use_highbitdepth)
    extend_frame_highbd(CONVERT_TO_SHORTPTR(dgd->y_buffer), width, height,
                        dgd_stride);
  else
#endif
    extend_frame(dgd->y_buffer, width, height, dgd_stride);

937
938
  // Compute best Wiener filters for each tile
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
939
940
941
    av1_get_rest_tile_limits(tile_idx, 0, 0, nhtiles, nvtiles, tile_width,
                             tile_height, width, height, 0, 0, &h_start, &h_end,
                             &v_start, &v_end);
942
    err = sse_restoration_tile(src, cm->frame_to_show, cm, h_start,
943
                               h_end - h_start, v_start, v_end - v_start, 1);
944
945
    // #bits when a tile is not restored
    bits = av1_cost_bit(RESTORE_NONE_WIENER_PROB, 0);
946
    cost_norestore = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);
947
    best_tile_cost[tile_idx] = DBL_MAX;
948
949

    av1_get_rest_tile_limits(tile_idx, 0, 0, nhtiles, nvtiles, tile_width,
950
951
                             tile_height, width, height, 0, 0, &h_start, &h_end,
                             &v_start, &v_end);
952
#if CONFIG_HIGHBITDEPTH
953
954
955
956
    if (cm->use_highbitdepth)
      compute_stats_highbd(dgd->y_buffer, src->y_buffer, h_start, h_end,
                           v_start, v_end, dgd_stride, src_stride, M, H);
    else
957
#endif  // CONFIG_HIGHBITDEPTH
958
959
960
      compute_stats(dgd->y_buffer, src->y_buffer, h_start, h_end, v_start,
                    v_end, dgd_stride, src_stride, M, H);

961
962
    type[tile_idx] = RESTORE_WIENER;

963
    if (!wiener_decompose_sep_sym(M, H, vfilterd, hfilterd)) {
964
      type[tile_idx] = RESTORE_NONE;
965
966
      continue;
    }
967
968
    quantize_sym_filter(vfilterd, rsi->wiener_info[tile_idx].vfilter);
    quantize_sym_filter(hfilterd, rsi->wiener_info[tile_idx].hfilter);
969
970
971
972

    // Filter score computes the value of the function x'*A*x - x'*b for the
    // learned filter and compares it against identity filer. If there is no
    // reduction in the function, the filter is reverted back to identity
973
974
    score = compute_score(M, H, rsi->wiener_info[tile_idx].vfilter,
                          rsi->wiener_info[tile_idx].hfilter);
975
    if (score > 0.0) {
976
      type[tile_idx] = RESTORE_NONE;
977
978
      continue;
    }
979

980
    rsi->restoration_type[tile_idx] = RESTORE_WIENER;
981
    err = try_restoration_tile(src, cpi, rsi, 1, partial_frame, tile_idx, 0, 0,
982
                               dst_frame);
983
984
985
    bits = WIENER_FILT_BITS << AV1_PROB_COST_SHIFT;
    bits += av1_cost_bit(RESTORE_NONE_WIENER_PROB, 1);
    cost_wiener = RDCOST_DBL(x->rdmult