pickrst.c 43.7 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

24
25
#include "av1/common/onyxc_int.h"
#include "av1/common/quant_common.h"
26

27
28
29
30
#include "av1/encoder/encoder.h"
#include "av1/encoder/picklpf.h"
#include "av1/encoder/pickrst.h"
#include "av1/encoder/quantize.h"
31

32
33
34
typedef double (*search_restore_type)(const YV12_BUFFER_CONFIG *src,
                                      AV1_COMP *cpi, int filter_level,
                                      int partial_frame, RestorationInfo *info,
35
36
                                      double *best_tile_cost,
                                      YV12_BUFFER_CONFIG *dst_frame);
37

38
const int frame_level_restore_bits[RESTORE_TYPES] = { 2, 2, 3, 3, 2 };
39
40

static int64_t sse_restoration_tile(const YV12_BUFFER_CONFIG *src,
41
42
                                    const YV12_BUFFER_CONFIG *dst,
                                    const AV1_COMMON *cm, int h_start,
43
44
                                    int width, int v_start, int height,
                                    int y_only) {
45
46
47
  int64_t filt_err;
#if CONFIG_AOM_HIGHBITDEPTH
  if (cm->use_highbitdepth) {
48
49
    filt_err =
        aom_highbd_get_y_sse_part(src, dst, h_start, width, v_start, height);
50
51
52
53
54
55
56
57
58
    if (!y_only) {
      filt_err += aom_highbd_get_u_sse_part(
          src, dst, h_start >> cm->subsampling_x, width >> cm->subsampling_x,
          v_start >> cm->subsampling_y, height >> cm->subsampling_y);
      filt_err += aom_highbd_get_v_sse_part(
          src, dst, h_start >> cm->subsampling_x, width >> cm->subsampling_x,
          v_start >> cm->subsampling_y, height >> cm->subsampling_y);
    }
    return filt_err;
59
60
  }
#endif  // CONFIG_AOM_HIGHBITDEPTH
61
62
63
64
65
66
67
68
69
  filt_err = aom_get_y_sse_part(src, dst, h_start, width, v_start, height);
  if (!y_only) {
    filt_err += aom_get_u_sse_part(
        src, dst, h_start >> cm->subsampling_x, width >> cm->subsampling_x,
        v_start >> cm->subsampling_y, height >> cm->subsampling_y);
    filt_err += aom_get_u_sse_part(
        src, dst, h_start >> cm->subsampling_x, width >> cm->subsampling_x,
        v_start >> cm->subsampling_y, height >> cm->subsampling_y);
  }
70
71
72
73
74
  return filt_err;
}

static int64_t try_restoration_tile(const YV12_BUFFER_CONFIG *src,
                                    AV1_COMP *const cpi, RestorationInfo *rsi,
75
                                    int y_only, int partial_frame, int tile_idx,
76
77
                                    int subtile_idx, int subtile_bits,
                                    YV12_BUFFER_CONFIG *dst_frame) {
78
79
80
81
82
83
84
85
  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;
  const int ntiles = av1_get_rest_ntiles(cm->width, cm->height, &tile_width,
                                         &tile_height, &nhtiles, &nvtiles);
  (void)ntiles;

86
  av1_loop_restoration_frame(cm->frame_to_show, cm, rsi, y_only, partial_frame,
87
                             dst_frame);
88
89
90
91
  av1_get_rest_tile_limits(tile_idx, subtile_idx, subtile_bits, nhtiles,
                           nvtiles, tile_width, tile_height, cm->width,
                           cm->height, 0, 0, &h_start, &h_end, &v_start,
                           &v_end);
92
  filt_err = sse_restoration_tile(src, dst_frame, cm, h_start, h_end - h_start,
93
                                  v_start, v_end - v_start, y_only);
94
95
96
97
98

  return filt_err;
}

static int64_t try_restoration_frame(const YV12_BUFFER_CONFIG *src,
Yaowu Xu's avatar
Yaowu Xu committed
99
                                     AV1_COMP *const cpi, RestorationInfo *rsi,
100
                                     int y_only, int partial_frame,
101
                                     YV12_BUFFER_CONFIG *dst_frame) {
Yaowu Xu's avatar
Yaowu Xu committed
102
  AV1_COMMON *const cm = &cpi->common;
103
  int64_t filt_err;
104
  av1_loop_restoration_frame(cm->frame_to_show, cm, rsi, y_only, partial_frame,
105
                             dst_frame);
Yaowu Xu's avatar
Yaowu Xu committed
106
#if CONFIG_AOM_HIGHBITDEPTH
107
  if (cm->use_highbitdepth) {
108
    filt_err = aom_highbd_get_y_sse(src, dst_frame);
109
110
111
112
113
    if (!y_only) {
      filt_err += aom_highbd_get_u_sse(src, dst_frame);
      filt_err += aom_highbd_get_v_sse(src, dst_frame);
    }
    return filt_err;
114
  }
Yaowu Xu's avatar
Yaowu Xu committed
115
#endif  // CONFIG_AOM_HIGHBITDEPTH
116
117
118
119
120
  filt_err = aom_get_y_sse(src, dst_frame);
  if (!y_only) {
    filt_err += aom_get_u_sse(src, dst_frame);
    filt_err += aom_get_v_sse(src, dst_frame);
  }
121
122
123
  return filt_err;
}

124
125
126
127
static int64_t get_pixel_proj_error(int32_t *src, int width, int height,
                                    int src_stride, int32_t *dgd,
                                    int dgd_stride, int32_t *flt1,
                                    int flt1_stride, int32_t *flt2,
128
129
130
131
132
133
134
                                    int flt2_stride, int *xqd) {
  int i, j;
  int64_t err = 0;
  int xq[2];
  decode_xq(xqd, xq);
  for (i = 0; i < height; ++i) {
    for (j = 0; j < width; ++j) {
135
136
137
138
      const int32_t s = (int32_t)src[i * src_stride + j];
      const int32_t u = (int32_t)dgd[i * dgd_stride + j];
      const int32_t f1 = (int32_t)flt1[i * flt1_stride + j] - u;
      const int32_t f2 = (int32_t)flt2[i * flt2_stride + j] - u;
139
      const int64_t v = xq[0] * f1 + xq[1] * f2 + (u << SGRPROJ_PRJ_BITS);
140
      const int32_t e =
141
142
143
144
145
146
147
148
          ROUND_POWER_OF_TWO(v, SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS) -
          ROUND_POWER_OF_TWO(s, SGRPROJ_RST_BITS);
      err += e * e;
    }
  }
  return err;
}

149
150
151
static void get_proj_subspace(int32_t *src, int width, int height,
                              int src_stride, int32_t *dgd, int dgd_stride,
                              int32_t *flt1, int flt1_stride, int32_t *flt2,
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
                              int flt2_stride, int *xq) {
  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;

  xq[0] = -(1 << SGRPROJ_PRJ_BITS) / 4;
  xq[1] = (1 << SGRPROJ_PRJ_BITS) - xq[0];
  for (i = 0; i < height; ++i) {
    for (j = 0; j < width; ++j) {
      const double u = (double)dgd[i * dgd_stride + j];
      const double s = (double)src[i * src_stride + j] - 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;
    }
  }
  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) {
  xqd[0] = -xq[0];
  xqd[0] = clamp(xqd[0], SGRPROJ_PRJ_MIN0, SGRPROJ_PRJ_MAX0);
  xqd[1] = (1 << SGRPROJ_PRJ_BITS) + xqd[0] - xq[1];
  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,
199
                                          int *eps, int *xqd, void *srcbuf,
200
                                          void *rstbuf) {
201
202
203
204
  int32_t *srd = (int32_t *)srcbuf;
  int32_t *dgd = (int32_t *)rstbuf;
  int32_t *flt1 = dgd + RESTORATION_TILEPELS_MAX;
  int32_t *flt2 = flt1 + RESTORATION_TILEPELS_MAX;
205
206
207
208
209
210
211
212
213
214
215
  uint8_t *tmpbuf2 = (uint8_t *)(flt2 + RESTORATION_TILEPELS_MAX);
  int i, j, ep, bestep = 0;
  int64_t err, besterr = -1;
  int exqd[2], bestxqd[2] = { 0, 0 };
  for (ep = 0; ep < SGRPROJ_PARAMS; ep++) {
    int exq[2];
    if (bit_depth > 8) {
      uint16_t *src = CONVERT_TO_SHORTPTR(src8);
      uint16_t *dat = CONVERT_TO_SHORTPTR(dat8);
      for (i = 0; i < height; ++i) {
        for (j = 0; j < width; ++j) {
216
217
218
          flt1[i * width + j] = (int32_t)dat[i * dat_stride + j];
          flt2[i * width + j] = (int32_t)dat[i * dat_stride + j];
          dgd[i * width + j] = (int32_t)dat[i * dat_stride + j]
219
                               << SGRPROJ_RST_BITS;
220
          srd[i * width + j] = (int32_t)src[i * src_stride + j]
221
222
223
224
225
226
227
228
229
230
                               << SGRPROJ_RST_BITS;
        }
      }
    } else {
      uint8_t *src = src8;
      uint8_t *dat = dat8;
      for (i = 0; i < height; ++i) {
        for (j = 0; j < width; ++j) {
          const int k = i * width + j;
          const int l = i * dat_stride + j;
231
232
233
234
          flt1[k] = (int32_t)dat[l];
          flt2[k] = (int32_t)dat[l];
          dgd[k] = (int32_t)dat[l] << SGRPROJ_RST_BITS;
          srd[k] = (int32_t)src[i * src_stride + j] << SGRPROJ_RST_BITS;
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
        }
      }
    }
    av1_selfguided_restoration(flt1, width, height, width, bit_depth,
                               sgr_params[ep].r1, sgr_params[ep].e1, tmpbuf2);
    av1_selfguided_restoration(flt2, width, height, width, bit_depth,
                               sgr_params[ep].r2, sgr_params[ep].e2, tmpbuf2);
    get_proj_subspace(srd, width, height, width, dgd, width, flt1, width, flt2,
                      width, exq);
    encode_xq(exq, exqd);
    err = get_pixel_proj_error(srd, width, height, width, dgd, width, flt1,
                               width, flt2, width, exqd);
    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,
                             int filter_level, int partial_frame,
261
262
                             RestorationInfo *info, double *best_tile_cost,
                             YV12_BUFFER_CONFIG *dst_frame) {
263
264
265
266
267
268
  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;
269
  RestorationInfo *rsi = &cpi->rst_search;
270
271
  int tile_idx, tile_width, tile_height, nhtiles, nvtiles;
  int h_start, h_end, v_start, v_end;
272
  // Allocate for the src buffer at high precision
273
274
275
276
277
278
279
280
  const int ntiles = av1_get_rest_ntiles(cm->width, cm->height, &tile_width,
                                         &tile_height, &nhtiles, &nvtiles);
  //  Make a copy of the unfiltered / processed recon buffer
  aom_yv12_copy_y(cm->frame_to_show, &cpi->last_frame_uf);
  av1_loop_filter_frame(cm->frame_to_show, cm, &cpi->td.mb.e_mbd, filter_level,
                        1, partial_frame);
  aom_yv12_copy_y(cm->frame_to_show, &cpi->last_frame_db);

281
  rsi->frame_restoration_type = RESTORE_SGRPROJ;
282
283

  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx)
284
    rsi->sgrproj_info[tile_idx].level = 0;
285
286
287
288
289
  // 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);
290
    err = sse_restoration_tile(src, cm->frame_to_show, cm, h_start,
291
                               h_end - h_start, v_start, v_end - v_start, 1);
292
293
294
295
296
297
298
299
300
301
302
303
304
    // #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,
#if CONFIG_AOM_HIGHBITDEPTH
        cm->bit_depth,
#else
        8,
#endif  // CONFIG_AOM_HIGHBITDEPTH
305
306
307
308
        &rsi->sgrproj_info[tile_idx].ep, rsi->sgrproj_info[tile_idx].xqd,
        cpi->extra_rstbuf, cm->rst_internal.tmpbuf);
    rsi->sgrproj_info[tile_idx].level = 1;
    err = try_restoration_tile(src, cpi, rsi, 1, partial_frame, tile_idx, 0, 0,
309
                               dst_frame);
310
311
312
313
314
315
    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) {
      sgrproj_info[tile_idx].level = 0;
    } else {
316
      memcpy(&sgrproj_info[tile_idx], &rsi->sgrproj_info[tile_idx],
317
318
319
320
321
322
             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);
    }
323
    rsi->sgrproj_info[tile_idx].level = 0;
324
325
  }
  // Cost for Sgrproj filtering
326
  bits = frame_level_restore_bits[rsi->frame_restoration_type]
327
328
329
330
         << AV1_PROB_COST_SHIFT;
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
    bits +=
        av1_cost_bit(RESTORE_NONE_SGRPROJ_PROB, sgrproj_info[tile_idx].level);
331
    memcpy(&rsi->sgrproj_info[tile_idx], &sgrproj_info[tile_idx],
332
333
334
335
336
           sizeof(sgrproj_info[tile_idx]));
    if (sgrproj_info[tile_idx].level) {
      bits += (SGRPROJ_BITS << AV1_PROB_COST_SHIFT);
    }
  }
337
  err = try_restoration_frame(src, cpi, rsi, 1, partial_frame, dst_frame);
338
339
340
341
342
343
  cost_sgrproj = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);

  aom_yv12_copy_y(&cpi->last_frame_uf, cm->frame_to_show);
  return cost_sgrproj;
}

344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
static int64_t compute_sse(uint8_t *dgd, int width, int height, int dgd_stride,
                           uint8_t *src, int src_stride) {
  int64_t sse = 0;
  int i, j;
  for (i = 0; i < height; ++i) {
    for (j = 0; j < width; ++j) {
      const int diff =
          (int)dgd[i * dgd_stride + j] - (int)src[i * src_stride + j];
      sse += diff * diff;
    }
  }
  return sse;
}

#if CONFIG_AOM_HIGHBITDEPTH
static int64_t compute_sse_highbd(uint16_t *dgd, int width, int height,
                                  int dgd_stride, uint16_t *src,
                                  int src_stride) {
  int64_t sse = 0;
  int i, j;
  for (i = 0; i < height; ++i) {
    for (j = 0; j < width; ++j) {
      const int diff =
          (int)dgd[i * dgd_stride + j] - (int)src[i * src_stride + j];
      sse += diff * diff;
    }
  }
  return sse;
}
#endif  // CONFIG_AOM_HIGHBITDEPTH

static void search_domaintxfmrf_restoration(uint8_t *dgd8, int width,
                                            int height, int dgd_stride,
                                            uint8_t *src8, int src_stride,
378
379
                                            int bit_depth, int *sigma_r,
                                            void *fltbuf, void *rstbuf) {
380
381
382
383
384
  const int first_p_step = 8;
  const int second_p_range = first_p_step >> 1;
  const int second_p_step = 2;
  const int third_p_range = second_p_step >> 1;
  const int third_p_step = 1;
385
  int p, best_p0, best_p = -1;
386
387
  int64_t best_sse = INT64_MAX, sse;
  if (bit_depth == 8) {
388
389
    uint8_t *flt = (uint8_t *)fltbuf;
    int32_t *tmpbuf = (int32_t *)rstbuf;
390
391
392
393
    uint8_t *dgd = dgd8;
    uint8_t *src = src8;
    // First phase
    for (p = first_p_step / 2; p < DOMAINTXFMRF_PARAMS; p += first_p_step) {
394
      av1_domaintxfmrf_restoration(dgd, width, height, dgd_stride, p, flt,
395
                                   width, tmpbuf);
396
      sse = compute_sse(flt, width, height, width, src, src_stride);
397
398
399
400
401
402
403
404
405
406
      if (sse < best_sse || best_p == -1) {
        best_p = p;
        best_sse = sse;
      }
    }
    // Second Phase
    best_p0 = best_p;
    for (p = best_p0 - second_p_range; p <= best_p0 + second_p_range;
         p += second_p_step) {
      if (p < 0 || p == best_p || p >= DOMAINTXFMRF_PARAMS) continue;
407
      av1_domaintxfmrf_restoration(dgd, width, height, dgd_stride, p, flt,
408
                                   width, tmpbuf);
409
      sse = compute_sse(flt, width, height, width, src, src_stride);
410
411
412
413
414
415
416
417
418
419
      if (sse < best_sse) {
        best_p = p;
        best_sse = sse;
      }
    }
    // Third Phase
    best_p0 = best_p;
    for (p = best_p0 - third_p_range; p <= best_p0 + third_p_range;
         p += third_p_step) {
      if (p < 0 || p == best_p || p >= DOMAINTXFMRF_PARAMS) continue;
420
      av1_domaintxfmrf_restoration(dgd, width, height, dgd_stride, p, flt,
421
                                   width, tmpbuf);
422
      sse = compute_sse(flt, width, height, width, src, src_stride);
423
424
425
426
427
428
429
      if (sse < best_sse) {
        best_p = p;
        best_sse = sse;
      }
    }
  } else {
#if CONFIG_AOM_HIGHBITDEPTH
430
431
    uint16_t *flt = (uint16_t *)fltbuf;
    int32_t *tmpbuf = (int32_t *)rstbuf;
432
433
434
435
    uint16_t *dgd = CONVERT_TO_SHORTPTR(dgd8);
    uint16_t *src = CONVERT_TO_SHORTPTR(src8);
    // First phase
    for (p = first_p_step / 2; p < DOMAINTXFMRF_PARAMS; p += first_p_step) {
436
      av1_domaintxfmrf_restoration_highbd(dgd, width, height, dgd_stride, p,
437
438
                                          bit_depth, flt, width, tmpbuf);
      sse = compute_sse_highbd(flt, width, height, width, src, src_stride);
439
440
441
442
443
444
445
446
447
448
      if (sse < best_sse || best_p == -1) {
        best_p = p;
        best_sse = sse;
      }
    }
    // Second Phase
    best_p0 = best_p;
    for (p = best_p0 - second_p_range; p <= best_p0 + second_p_range;
         p += second_p_step) {
      if (p < 0 || p == best_p || p >= DOMAINTXFMRF_PARAMS) continue;
449
      av1_domaintxfmrf_restoration_highbd(dgd, width, height, dgd_stride, p,
450
451
                                          bit_depth, flt, width, tmpbuf);
      sse = compute_sse_highbd(flt, width, height, width, src, src_stride);
452
453
454
455
456
457
458
459
460
461
      if (sse < best_sse) {
        best_p = p;
        best_sse = sse;
      }
    }
    // Third Phase
    best_p0 = best_p;
    for (p = best_p0 - third_p_range; p <= best_p0 + third_p_range;
         p += third_p_step) {
      if (p < 0 || p == best_p || p >= DOMAINTXFMRF_PARAMS) continue;
462
      av1_domaintxfmrf_restoration_highbd(dgd, width, height, dgd_stride, p,
463
464
                                          bit_depth, flt, width, tmpbuf);
      sse = compute_sse_highbd(flt, width, height, width, src, src_stride);
465
466
467
468
469
470
471
472
473
474
475
476
477
478
      if (sse < best_sse) {
        best_p = p;
        best_sse = sse;
      }
    }
#else
    assert(0);
#endif  // CONFIG_AOM_HIGHBITDEPTH
  }
  *sigma_r = best_p;
}

static double search_domaintxfmrf(const YV12_BUFFER_CONFIG *src, AV1_COMP *cpi,
                                  int filter_level, int partial_frame,
479
480
                                  RestorationInfo *info, double *best_tile_cost,
                                  YV12_BUFFER_CONFIG *dst_frame) {
481
  DomaintxfmrfInfo *domaintxfmrf_info = info->domaintxfmrf_info;
482
483
  double cost_norestore, cost_domaintxfmrf;
  int64_t err;
484
485
486
487
  int bits;
  MACROBLOCK *x = &cpi->td.mb;
  AV1_COMMON *const cm = &cpi->common;
  const YV12_BUFFER_CONFIG *dgd = cm->frame_to_show;
488
  RestorationInfo *rsi = &cpi->rst_search;
489
490
491
492
493
494
495
496
497
498
  int tile_idx, tile_width, tile_height, nhtiles, nvtiles;
  int h_start, h_end, v_start, v_end;
  const int ntiles = av1_get_rest_ntiles(cm->width, cm->height, &tile_width,
                                         &tile_height, &nhtiles, &nvtiles);
  //  Make a copy of the unfiltered / processed recon buffer
  aom_yv12_copy_y(cm->frame_to_show, &cpi->last_frame_uf);
  av1_loop_filter_frame(cm->frame_to_show, cm, &cpi->td.mb.e_mbd, filter_level,
                        1, partial_frame);
  aom_yv12_copy_y(cm->frame_to_show, &cpi->last_frame_db);

499
500
  rsi->frame_restoration_type = RESTORE_DOMAINTXFMRF;
  rsi->domaintxfmrf_info = cpi->rst_search.domaintxfmrf_info;
501
502

  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx)
503
    rsi->domaintxfmrf_info[tile_idx].level = 0;
504
505
506
507
508
  // Compute best Domaintxfm 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);
509
    err = sse_restoration_tile(src, cm->frame_to_show, cm, h_start,
510
                               h_end - h_start, v_start, v_end - v_start, 1);
511
512
513
514
515
516
517
518
519
520
521
522
523
524
    // #bits when a tile is not restored
    bits = av1_cost_bit(RESTORE_NONE_DOMAINTXFMRF_PROB, 0);
    cost_norestore = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);
    best_tile_cost[tile_idx] = DBL_MAX;

    search_domaintxfmrf_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,
#if CONFIG_AOM_HIGHBITDEPTH
        cm->bit_depth,
#else
        8,
#endif  // CONFIG_AOM_HIGHBITDEPTH
525
526
        &rsi->domaintxfmrf_info[tile_idx].sigma_r, cpi->extra_rstbuf,
        cm->rst_internal.tmpbuf);
527

528
529
    rsi->domaintxfmrf_info[tile_idx].level = 1;
    err = try_restoration_tile(src, cpi, rsi, 1, partial_frame, tile_idx, 0, 0,
530
                               dst_frame);
531
532
533
534
535
536
    bits = DOMAINTXFMRF_PARAMS_BITS << AV1_PROB_COST_SHIFT;
    bits += av1_cost_bit(RESTORE_NONE_DOMAINTXFMRF_PROB, 1);
    cost_domaintxfmrf = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);
    if (cost_domaintxfmrf >= cost_norestore) {
      domaintxfmrf_info[tile_idx].level = 0;
    } else {
537
      memcpy(&domaintxfmrf_info[tile_idx], &rsi->domaintxfmrf_info[tile_idx],
538
539
540
541
542
543
544
             sizeof(domaintxfmrf_info[tile_idx]));
      bits = DOMAINTXFMRF_PARAMS_BITS << AV1_PROB_COST_SHIFT;
      best_tile_cost[tile_idx] = RDCOST_DBL(
          x->rdmult, x->rddiv,
          (bits + cpi->switchable_restore_cost[RESTORE_DOMAINTXFMRF]) >> 4,
          err);
    }
545
    rsi->domaintxfmrf_info[tile_idx].level = 0;
546
547
  }
  // Cost for Domaintxfmrf filtering
548
  bits = frame_level_restore_bits[rsi->frame_restoration_type]
549
550
551
552
         << AV1_PROB_COST_SHIFT;
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
    bits += av1_cost_bit(RESTORE_NONE_DOMAINTXFMRF_PROB,
                         domaintxfmrf_info[tile_idx].level);
553
    memcpy(&rsi->domaintxfmrf_info[tile_idx], &domaintxfmrf_info[tile_idx],
554
555
556
557
558
           sizeof(domaintxfmrf_info[tile_idx]));
    if (domaintxfmrf_info[tile_idx].level) {
      bits += (DOMAINTXFMRF_PARAMS_BITS << AV1_PROB_COST_SHIFT);
    }
  }
559
  err = try_restoration_frame(src, cpi, rsi, 1, partial_frame, dst_frame);
560
561
562
563
564
565
  cost_domaintxfmrf = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);

  aom_yv12_copy_y(&cpi->last_frame_uf, cm->frame_to_show);
  return cost_domaintxfmrf;
}

566
567
static double find_average(uint8_t *src, int h_start, int h_end, int v_start,
                           int v_end, int stride) {
568
569
570
  uint64_t sum = 0;
  double avg = 0;
  int i, j;
571
572
573
  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));
574
575
576
  return avg;
}

577
578
579
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) {
580
  int i, j, k, l;
581
  double Y[WIENER_WIN2];
582
583
  const double avg =
      find_average(dgd, h_start, h_end, v_start, v_end, dgd_stride);
584

585
586
  memset(M, 0, sizeof(*M) * WIENER_WIN2);
  memset(H, 0, sizeof(*H) * WIENER_WIN2 * WIENER_WIN2);
587
588
  for (i = v_start; i < v_end; i++) {
    for (j = h_start; j < h_end; j++) {
589
590
      const double X = (double)src[i * src_stride + j] - avg;
      int idx = 0;
591
592
      for (k = -WIENER_HALFWIN; k <= WIENER_HALFWIN; k++) {
        for (l = -WIENER_HALFWIN; l <= WIENER_HALFWIN; l++) {
593
594
595
596
          Y[idx] = (double)dgd[(i + l) * dgd_stride + (j + k)] - avg;
          idx++;
        }
      }
597
      for (k = 0; k < WIENER_WIN2; ++k) {
598
        M[k] += Y[k] * X;
599
600
        H[k * WIENER_WIN2 + k] += Y[k] * Y[k];
        for (l = k + 1; l < WIENER_WIN2; ++l) {
601
          double value = Y[k] * Y[l];
602
603
          H[k * WIENER_WIN2 + l] += value;
          H[l * WIENER_WIN2 + k] += value;
604
605
606
607
608
609
        }
      }
    }
  }
}

Yaowu Xu's avatar
Yaowu Xu committed
610
#if CONFIG_AOM_HIGHBITDEPTH
611
612
static double find_average_highbd(uint16_t *src, int h_start, int h_end,
                                  int v_start, int v_end, int stride) {
613
614
615
  uint64_t sum = 0;
  double avg = 0;
  int i, j;
616
617
618
  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));
619
620
621
  return avg;
}

622
623
624
625
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) {
626
  int i, j, k, l;
627
  double Y[WIENER_WIN2];
628
629
  uint16_t *src = CONVERT_TO_SHORTPTR(src8);
  uint16_t *dgd = CONVERT_TO_SHORTPTR(dgd8);
630
631
  const double avg =
      find_average_highbd(dgd, h_start, h_end, v_start, v_end, dgd_stride);
632

633
634
  memset(M, 0, sizeof(*M) * WIENER_WIN2);
  memset(H, 0, sizeof(*H) * WIENER_WIN2 * WIENER_WIN2);
635
636
  for (i = v_start; i < v_end; i++) {
    for (j = h_start; j < h_end; j++) {
637
638
      const double X = (double)src[i * src_stride + j] - avg;
      int idx = 0;
639
640
      for (k = -WIENER_HALFWIN; k <= WIENER_HALFWIN; k++) {
        for (l = -WIENER_HALFWIN; l <= WIENER_HALFWIN; l++) {
641
642
643
644
          Y[idx] = (double)dgd[(i + l) * dgd_stride + (j + k)] - avg;
          idx++;
        }
      }
645
      for (k = 0; k < WIENER_WIN2; ++k) {
646
        M[k] += Y[k] * X;
647
648
        H[k * WIENER_WIN2 + k] += Y[k] * Y[k];
        for (l = k + 1; l < WIENER_WIN2; ++l) {
649
          double value = Y[k] * Y[l];
650
651
          H[k * WIENER_WIN2 + l] += value;
          H[l * WIENER_WIN2 + k] += value;
652
653
654
655
656
        }
      }
    }
  }
}
Yaowu Xu's avatar
Yaowu Xu committed
657
#endif  // CONFIG_AOM_HIGHBITDEPTH
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679

// 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;
  // Partial pivoting
  for (i = n - 1; i > 0; i--) {
    if (A[(i - 1) * stride] < A[i * stride]) {
      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;
    }
  }
  // Forward elimination
  for (k = 0; k < n - 1; k++) {
    for (i = k; i < n - 1; i++) {
      c = A[(i + 1) * stride + k] / A[k * stride + k];
680
      for (j = 0; j < n; j++) A[(i + 1) * stride + j] -= c * A[k * stride + j];
681
682
683
684
685
      b[i + 1] -= c * b[k];
    }
  }
  // Backward substitution
  for (i = n - 1; i >= 0; i--) {
686
    if (fabs(A[i * stride + i]) < 1e-10) return 0;
687
    c = 0;
688
    for (j = i + 1; j <= n - 1; j++) c += A[i * stride + j] * x[j];
689
690
691
692
693
694
    x[i] = (b[i] - c) / A[i * stride + i];
  }
  return 1;
}

static INLINE int wrap_index(int i) {
695
  return (i >= WIENER_HALFWIN1 ? WIENER_WIN - 1 - i : i);
696
697
698
699
700
}

// Fix vector b, update vector a
static void update_a_sep_sym(double **Mc, double **Hc, double *a, double *b) {
  int i, j;
701
702
  double S[WIENER_WIN];
  double A[WIENER_WIN], B[WIENER_WIN2];
Aamir Anis's avatar
Aamir Anis committed
703
  int w, w2;
704
705
  memset(A, 0, sizeof(A));
  memset(B, 0, sizeof(B));
706
707
  for (i = 0; i < WIENER_WIN; i++) {
    for (j = 0; j < WIENER_WIN; ++j) {
708
709
710
711
      const int jj = wrap_index(j);
      A[jj] += Mc[i][j] * b[i];
    }
  }
712
713
  for (i = 0; i < WIENER_WIN; i++) {
    for (j = 0; j < WIENER_WIN; j++) {
714
      int k, l;
715
716
      for (k = 0; k < WIENER_WIN; ++k)
        for (l = 0; l < WIENER_WIN; ++l) {
717
718
          const int kk = wrap_index(k);
          const int ll = wrap_index(l);
719
720
          B[ll * WIENER_HALFWIN1 + kk] +=
              Hc[j * WIENER_WIN + i][k * WIENER_WIN2 + l] * b[i] * b[j];
721
722
723
        }
    }
  }
Aamir Anis's avatar
Aamir Anis committed
724
  // Normalization enforcement in the system of equations itself
725
  w = WIENER_WIN;
Aamir Anis's avatar
Aamir Anis committed
726
727
  w2 = (w >> 1) + 1;
  for (i = 0; i < w2 - 1; ++i)
728
729
    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
730
731
732
733
734
735
736
737
738
  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];
739
    }
Aamir Anis's avatar
Aamir Anis committed
740
    memcpy(a, S, w * sizeof(*a));
741
742
743
744
745
746
  }
}

// Fix vector a, update vector b
static void update_b_sep_sym(double **Mc, double **Hc, double *a, double *b) {
  int i, j;
747
748
  double S[WIENER_WIN];
  double A[WIENER_WIN], B[WIENER_WIN2];
Aamir Anis's avatar
Aamir Anis committed
749
  int w, w2;
750
751
  memset(A, 0, sizeof(A));
  memset(B, 0, sizeof(B));
752
  for (i = 0; i < WIENER_WIN; i++) {
753
    const int ii = wrap_index(i);
754
    for (j = 0; j < WIENER_WIN; j++) A[ii] += Mc[i][j] * a[j];
755
756
  }

757
758
  for (i = 0; i < WIENER_WIN; i++) {
    for (j = 0; j < WIENER_WIN; j++) {
759
760
761
      const int ii = wrap_index(i);
      const int jj = wrap_index(j);
      int k, l;
762
763
764
765
      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];
766
767
    }
  }
Aamir Anis's avatar
Aamir Anis committed
768
  // Normalization enforcement in the system of equations itself
769
770
  w = WIENER_WIN;
  w2 = WIENER_HALFWIN1;
Aamir Anis's avatar
Aamir Anis committed
771
  for (i = 0; i < w2 - 1; ++i)
772
773
    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
774
775
776
777
778
779
780
781
782
  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];
783
    }
Aamir Anis's avatar
Aamir Anis committed
784
    memcpy(b, S, w * sizeof(*b));
785
786
787
  }
}

788
789
static int wiener_decompose_sep_sym(double *M, double *H, double *a,
                                    double *b) {
790
  static const double init_filt[WIENER_WIN] = {
791
    0.035623, -0.127154, 0.211436, 0.760190, 0.211436, -0.127154, 0.035623,
792
793
  };
  int i, j, iter;
794
795
796
797
798
799
800
  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;
801
802
    }
  }
803
804
  memcpy(a, init_filt, sizeof(*a) * WIENER_WIN);
  memcpy(b, init_filt, sizeof(*b) * WIENER_WIN);
805
806
807
808
809
810
811

  iter = 1;
  while (iter < 10) {
    update_a_sep_sym(Mc, Hc, a, b);
    update_b_sep_sym(Mc, Hc, a, b);
    iter++;
  }
812
  return 1;
813
814
}

Aamir Anis's avatar
Aamir Anis committed
815
816
817
// Computes the function x'*A*x - x'*b for the learned filters, and compares
// against identity filters; Final score is defined as the difference between
// the function values
818
static double compute_score(double *M, double *H, int *vfilt, int *hfilt) {
819
  double ab[WIENER_WIN * WIENER_WIN];
Aamir Anis's avatar
Aamir Anis committed
820
821
822
823
824
  int i, k, l;
  double P = 0, Q = 0;
  double iP = 0, iQ = 0;
  double Score, iScore;
  int w;
825
826
827
828
829
830
831
832
  double a[WIENER_WIN], b[WIENER_WIN];
  w = WIENER_WIN;
  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
833
834
  }
  for (k = 0; k < w; ++k) {
835
    for (l = 0; l < w; ++l) ab[k * w + l] = a[l] * b[k];
Aamir Anis's avatar
Aamir Anis committed
836
837
838
  }
  for (k = 0; k < w * w; ++k) {
    P += ab[k] * M[k];
839
    for (l = 0; l < w * w; ++l) Q += ab[k] * H[k * w * w + l] * ab[l];
Aamir Anis's avatar
Aamir Anis committed
840
841
842
843
844
845
846
847
848
849
  }
  Score = Q - 2 * P;

  iP = M[(w * w) >> 1];
  iQ = H[((w * w) >> 1) * w * w + ((w * w) >> 1)];
  iScore = iQ - 2 * iP;

  return Score - iScore;
}

850
851
static void quantize_sym_filter(double *f, int *fi) {
  int i;
852
853
  for (i = 0; i < WIENER_HALFWIN; ++i) {
    fi[i] = RINT(f[i] * WIENER_FILT_STEP);
854
855
856
857
858
859
860
  }
  // 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);
}

861
862
static double search_wiener(const YV12_BUFFER_CONFIG *src, AV1_COMP *cpi,
                            int filter_level, int partial_frame,
863
864
                            RestorationInfo *info, double *best_tile_cost,
                            YV12_BUFFER_CONFIG *dst_frame) {
865
  WienerInfo *wiener_info = info->wiener_info;
Yaowu Xu's avatar
Yaowu Xu committed
866
  AV1_COMMON *const cm = &cpi->common;
867
  RestorationInfo *rsi = &cpi->rst_search;
868
869
  int64_t err;
  int bits;
870
  double cost_wiener, cost_norestore;
871
  MACROBLOCK *x = &cpi->td.mb;
872
873
874
  double M[WIENER_WIN2];
  double H[WIENER_WIN2 * WIENER_WIN2];
  double vfilterd[WIENER_WIN], hfilterd[WIENER_WIN];
875
876
877
878
879
  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
880
  double score;
881
  int tile_idx, tile_width, tile_height, nhtiles, nvtiles;
882
  int h_start, h_end, v_start, v_end;
883
  int i;
884

885
886
  const int ntiles = av1_get_rest_ntiles(width, height, &tile_width,
                                         &tile_height, &nhtiles, &nvtiles);
887
888
889
890
891
892
  assert(width == dgd->y_crop_width);
  assert(height == dgd->y_crop_height);
  assert(width == src->y_crop_width);
  assert(height == src->y_crop_height);

  //  Make a copy of the unfiltered / processed recon buffer
Yaowu Xu's avatar
Yaowu Xu committed
893
894
895
896
  aom_yv12_copy_y(cm->frame_to_show, &cpi->last_frame_uf);
  av1_loop_filter_frame(cm->frame_to_show, cm, &cpi->td.mb.e_mbd, filter_level,
                        1, partial_frame);
  aom_yv12_copy_y(cm->frame_to_show, &cpi->last_frame_db);
897

898
  rsi->frame_restoration_type = RESTORE_WIENER;
899

900
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx)
901
    rsi->wiener_info[tile_idx].level = 0;
902
903
904

  // Compute best Wiener filters for each tile
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
905
906
907
    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);
908
    err = sse_restoration_tile(src, cm->frame_to_show, cm, h_start,
909
                               h_end - h_start, v_start, v_end - v_start, 1);
910
911
    // #bits when a tile is not restored
    bits = av1_cost_bit(RESTORE_NONE_WIENER_PROB, 0);
912
    cost_norestore = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);
913
    best_tile_cost[tile_idx] = DBL_MAX;
914
915

    av1_get_rest_tile_limits(tile_idx, 0, 0, nhtiles, nvtiles, tile_width,
916
917
918
                             tile_height, width, height, WIENER_HALFWIN,
                             WIENER_HALFWIN, &h_start, &h_end, &v_start,
                             &v_end);
Yaowu Xu's avatar
Yaowu Xu committed
919
#if CONFIG_AOM_HIGHBITDEPTH
920
921
922
923
    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
Yaowu Xu's avatar
Yaowu Xu committed
924
#endif  // CONFIG_AOM_HIGHBITDEPTH
925
926
927
      compute_stats(dgd->y_buffer, src->y_buffer, h_start, h_end, v_start,
                    v_end, dgd_stride, src_stride, M, H);

928
    wiener_info[tile_idx].level = 1;
929
    if (!wiener_decompose_sep_sym(M, H, vfilterd, hfilterd)) {
930
      wiener_info[tile_idx].level = 0;
931
932
      continue;
    }
933
934
    quantize_sym_filter(vfilterd, rsi->wiener_info[tile_idx].vfilter);
    quantize_sym_filter(hfilterd, rsi->wiener_info[tile_idx].hfilter);
935
936
937
938

    // 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
939
940
    score = compute_score(M, H, rsi->wiener_info[tile_idx].vfilter,
                          rsi->wiener_info[tile_idx].hfilter);
941
    if (score > 0.0) {
942
      wiener_info[tile_idx].level = 0;
943
944
      continue;
    }
945

946
947
    rsi->wiener_info[tile_idx].level = 1;
    err = try_restoration_tile(src, cpi, rsi, 1, partial_frame, tile_idx, 0, 0,
948
                               dst_frame);
949
950
951
    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);
952
    if (cost_wiener >= cost_norestore) {
953
954
955
      wiener_info[tile_idx].level = 0;
    } else {
      wiener_info[tile_idx].level = 1;
956
957
958
959
960
      for (i = 0; i < WIENER_HALFWIN; ++i) {
        wiener_info[tile_idx].vfilter[i] =
            rsi->wiener_info[tile_idx].vfilter[i];
        wiener_info[tile_idx].hfilter[i] =
            rsi->wiener_info[tile_idx].hfilter[i];
961
      }
962
963
964
965
966
      bits = WIENER_FILT_BITS << AV1_PROB_COST_SHIFT;
      best_tile_cost[tile_idx] = RDCOST_DBL(
          x->rdmult, x->rddiv,
          (bits + cpi->switchable_restore_cost[RESTORE_WIENER]) >> 4, err);
    }
967
    rsi->wiener_info[tile_idx].level = 0;
968
  }
969
  // Cost for Wiener filtering
970
  bits = frame_level_restore_bits[rsi->frame_restoration_type]
971
         << AV1_PROB_COST_SHIFT;
972
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
973
    bits += av1_cost_bit(RESTORE_NONE_WIENER_PROB, wiener_info[tile_idx].level);
974
    rsi->wiener_info[tile_idx].level = wiener_info[tile_idx].level;
975
    if (wiener_info[tile_idx].level) {
976
      bits += (WIENER_FILT_BITS << AV1_PROB_COST_SHIFT);
977
978
979
980
981
      for (i = 0; i < WIENER_HALFWIN; ++i) {
        rsi->wiener_info[tile_idx].vfilter[i] =
            wiener_info[tile_idx].vfilter[i];
        rsi->wiener_info[tile_idx].hfilter[i] =
            wiener_info[tile_idx].hfilter[i];
982
983
      }
    }
Aamir Anis's avatar
Aamir Anis committed
984
  }
985
  err = try_restoration_frame(src, cpi, rsi, 1, partial_frame, dst_frame);
986
  cost_wiener = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);
987

Yaowu Xu's avatar
Yaowu Xu committed
988
  aom_yv12_copy_y(&cpi->last_frame_uf, cm->frame_to_show);
989
  return cost_wiener;
990
991
}

992
993
static double search_norestore(const YV12_BUFFER_CONFIG *src, AV1_COMP *cpi,
                               int filter_level, int partial_frame,
994
995
                               RestorationInfo *info, double *best_tile_cost,
                               YV12_BUFFER_CONFIG *dst_frame) {
996
997
  double err, cost_norestore;
  int bits;
998
  MACROBLOCK *x = &cpi->td.mb;
999
1000
  AV1_COMMON *const cm = &cpi->common;
  int tile_idx, tile_width, tile_height, nhtiles, nvtiles;
1001
1002
1003
  int h_start, h_end, v_start, v_end;
  const int ntiles = av1_get_rest_ntiles(cm->width, cm->height, &tile_width,
                                         &tile_height, &nhtiles, &nvtiles);
1004
  (void)info;
1005
  (void)dst_frame;
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016

  //  Make a copy of the unfiltered / processed recon buffer
  aom_yv12_copy_y(cm->frame_to_show, &cpi->last_frame_uf);
  av1_loop_filter_frame(cm->frame_to_show, cm, &cpi->td.mb.e_mbd, filter_level,
                        1, partial_frame);
  aom_yv12_copy_y(cm->frame_to_show, &cpi->last_frame_db);

  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);
1017
    err = sse_restoration_tile(src, cm->frame_to_show, cm, h_start,
1018
                               h_end - h_start, v_start, v_end - v_start, 1);
1019
    best_tile_cost[tile_idx] =
1020
1021
        RDCOST_DBL(x->rdmult, x->rddiv,
                   (cpi->switchable_restore_cost[RESTORE_NONE] >> 4), err);
1022
1023
  }
  // RD cost associated with no restoration
1024
  err = sse_restoration_tile(src, cm->frame_to_show, cm, 0, cm->width, 0,
1025
                             cm->height, 1);
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
  bits = frame_level_restore_bits[RESTORE_NONE] << AV1_PROB_COST_SHIFT;
  cost_norestore = RDCOST_DBL(x->rdmult, x->rddiv, (bits >> 4), err);
  aom_yv12_copy_y(&cpi->last_frame_uf, cm->frame_to_show);
  return cost_norestore;
}

static double search_switchable_restoration(
    AV1_COMP *cpi, int filter_level, int partial_frame, RestorationInfo *rsi,
    double *tile_cost[RESTORE_SWITCHABLE_TYPES]) {
  AV1_COMMON *const cm = &cpi->common;
  MACROBLOCK *x = &cpi->td.mb;
  double cost_switchable = 0;
  int r, bits, tile_idx;
  const int ntiles =
      av1_get_rest_ntiles(cm->width, cm->height, NULL, NULL, NULL, NULL);

  //  Make a copy of the unfiltered / processed recon buffer
  aom_yv12_copy_y(cm->frame_to_show, &cpi->last_frame_uf);
  av1_loop_filter_frame(cm->frame_to_show, cm, &cpi->td.mb.e_mbd, filter_level,
                        1, partial_frame);
  aom_yv12_copy_y(cm->frame_to_show, &cpi->last_frame_db);

  rsi->frame_restoration_type = RESTORE_SWITCHABLE;
  bits = frame_level_restore_bits[rsi->frame_restoration_type]
         << AV1_PROB_COST_SHIFT;
  cost_switchable = RDCOST_DBL(x->rdmult, x->rddiv, bits >> 4, 0);
  for (tile_idx = 0; tile_idx < ntiles; ++tile_idx) {
    double best_cost = tile_cost[RESTORE_NONE][tile_idx];
    rsi->restoration_type[tile_idx] = RESTORE_NONE;
    for (r = 1; r < RESTORE_SWITCHABLE_TYPES; r++) {
      if (tile_cost[r][tile_idx] < best_cost) {
        rsi->restoration_type[tile_idx] = r</