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;