third_party/dav1d/src/warpmv.c

   1 /*
   2  * Copyright © 2018, VideoLAN and dav1d authors
   3  * Copyright © 2018, Two Orioles, LLC
   4  * All rights reserved.
   5  *
   6  * Redistribution and use in source and binary forms, with or without
   7  * modification, are permitted provided that the following conditions are met:
   8  *
   9  * 1. Redistributions of source code must retain the above copyright notice, this
  10  *    list of conditions and the following disclaimer.
  11  *
  12  * 2. Redistributions in binary form must reproduce the above copyright notice,
  13  *    this list of conditions and the following disclaimer in the documentation
  14  *    and/or other materials provided with the distribution.
  15  *
  16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  19  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
  20  * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
  22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
  23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  26  */
  27
  28 #include "config.h"
  29
  30 #include <stdlib.h>
  31
  32 #include "common/intops.h"
  33
  34 #include "src/warpmv.h"
  35
  36 static const uint16_t div_lut[257] = {
  37     16384, 16320, 16257, 16194, 16132, 16070, 16009, 15948, 15888, 15828, 15768,
  38     15709, 15650, 15592, 15534, 15477, 15420, 15364, 15308, 15252, 15197, 15142,
  39     15087, 15033, 14980, 14926, 14873, 14821, 14769, 14717, 14665, 14614, 14564,
  40     14513, 14463, 14413, 14364, 14315, 14266, 14218, 14170, 14122, 14075, 14028,
  41     13981, 13935, 13888, 13843, 13797, 13752, 13707, 13662, 13618, 13574, 13530,
  42     13487, 13443, 13400, 13358, 13315, 13273, 13231, 13190, 13148, 13107, 13066,
  43     13026, 12985, 12945, 12906, 12866, 12827, 12788, 12749, 12710, 12672, 12633,
  44     12596, 12558, 12520, 12483, 12446, 12409, 12373, 12336, 12300, 12264, 12228,
  45     12193, 12157, 12122, 12087, 12053, 12018, 11984, 11950, 11916, 11882, 11848,
  46     11815, 11782, 11749, 11716, 11683, 11651, 11619, 11586, 11555, 11523, 11491,
  47     11460, 11429, 11398, 11367, 11336, 11305, 11275, 11245, 11215, 11185, 11155,
  48     11125, 11096, 11067, 11038, 11009, 10980, 10951, 10923, 10894, 10866, 10838,
  49     10810, 10782, 10755, 10727, 10700, 10673, 10645, 10618, 10592, 10565, 10538,
  50     10512, 10486, 10460, 10434, 10408, 10382, 10356, 10331, 10305, 10280, 10255,
  51     10230, 10205, 10180, 10156, 10131, 10107, 10082, 10058, 10034, 10010,  9986,
  52      9963,  9939,  9916,  9892,  9869,  9846,  9823,  9800,  9777,  9754,  9732,
  53      9709,  9687,  9664,  9642,  9620,  9598,  9576,  9554,  9533,  9511,  9489,
  54      9468,  9447,  9425,  9404,  9383,  9362,  9341,  9321,  9300,  9279,  9259,
  55      9239,  9218,  9198,  9178,  9158,  9138,  9118,  9098,  9079,  9059,  9039,
  56      9020,  9001,  8981,  8962,  8943,  8924,  8905,  8886,  8867,  8849,  8830,
  57      8812,  8793,  8775,  8756,  8738,  8720,  8702,  8684,  8666,  8648,  8630,
  58      8613,  8595,  8577,  8560,  8542,  8525,  8508,  8490,  8473,  8456,  8439,
  59      8422,  8405,  8389,  8372,  8355,  8339,  8322,  8306,  8289,  8273,  8257,
  60      8240,  8224,  8208,  8192,
  61 };
  62
  63 static inline int iclip_wmp(const int v) {
  64     const int cv = iclip(v, INT16_MIN, INT16_MAX);
  65
  66     return apply_sign((abs(cv) + 32) >> 6, cv) * (1 << 6);
  67 }
  68
  69 static inline int resolve_divisor_32(const unsigned d, int *const shift) {
  70     *shift = ulog2(d);
  71     const int e = d - (1 << *shift);
  72     const int f = *shift > 8 ? (e + (1 << (*shift - 9))) >> (*shift - 8) :
  73                                e << (8 - *shift);
  74     assert(f <= 256);
  75     *shift += 14;
  76     // Use f as lookup into the precomputed table of multipliers
  77     return div_lut[f];
  78 }
  79
  80 int dav1d_get_shear_params(Dav1dWarpedMotionParams *const wm) {
  81     const int32_t *const mat = wm->matrix;
  82
  83     if (mat[2] <= 0) return 1;
  84
  85     wm->u.p.alpha = iclip_wmp(mat[2] - 0x10000);
  86     wm->u.p.beta = iclip_wmp(mat[3]);
  87
  88     int shift;
  89     const int y = apply_sign(resolve_divisor_32(abs(mat[2]), &shift), mat[2]);
  90     const int64_t v1 = ((int64_t) mat[4] * 0x10000) * y;
  91     const int rnd = (1 << shift) >> 1;
  92     wm->u.p.gamma = iclip_wmp(apply_sign64((int) ((llabs(v1) + rnd) >> shift), v1));
  93     const int64_t v2 = ((int64_t) mat[3] * mat[4]) * y;
  94     wm->u.p.delta = iclip_wmp(mat[5] -
  95                           apply_sign64((int) ((llabs(v2) + rnd) >> shift), v2) -
  96                           0x10000);
  97
  98     return (4 * abs(wm->u.p.alpha) + 7 * abs(wm->u.p.beta) >= 0x10000) ||
  99            (4 * abs(wm->u.p.gamma) + 4 * abs(wm->u.p.delta) >= 0x10000);
 100 }
 101
 102 static int resolve_divisor_64(const uint64_t d, int *const shift) {
 103     *shift = u64log2(d);
 104     const int64_t e = d - (1LL << *shift);
 105     const int64_t f = *shift > 8 ? (e + (1LL << (*shift - 9))) >> (*shift - 8) :
 106                                    e << (8 - *shift);
 107     assert(f <= 256);
 108     *shift += 14;
 109     // Use f as lookup into the precomputed table of multipliers
 110     return div_lut[f];
 111 }
 112
 113 static int get_mult_shift_ndiag(const int64_t px,
 114                                 const int idet, const int shift)
 115 {
 116     const int64_t v1 = px * idet;
 117     const int v2 = apply_sign64((int) ((llabs(v1) +
 118                                         ((1LL << shift) >> 1)) >> shift),
 119                                 v1);
 120     return iclip(v2, -0x1fff, 0x1fff);
 121 }
 122
 123 static int get_mult_shift_diag(const int64_t px,
 124                                const int idet, const int shift)
 125 {
 126     const int64_t v1 = px * idet;
 127     const int v2 = apply_sign64((int) ((llabs(v1) +
 128                                         ((1LL << shift) >> 1)) >> shift),
 129                                 v1);
 130     return iclip(v2, 0xe001, 0x11fff);
 131 }
 132
 133 void dav1d_set_affine_mv2d(const int bw4, const int bh4,
 134                            const mv mv, Dav1dWarpedMotionParams *const wm,
 135                            const int bx4, const int by4)
 136 {
 137     int32_t *const mat = wm->matrix;
 138     const int rsuy = 2 * bh4 - 1;
 139     const int rsux = 2 * bw4 - 1;
 140     const int isuy = by4 * 4 + rsuy;
 141     const int isux = bx4 * 4 + rsux;
 142
 143     mat[0] = iclip(mv.x * 0x2000 - (isux * (mat[2] - 0x10000) + isuy * mat[3]),
 144                    -0x800000, 0x7fffff);
 145     mat[1] = iclip(mv.y * 0x2000 - (isux * mat[4] + isuy * (mat[5] - 0x10000)),
 146                    -0x800000, 0x7fffff);
 147 }
 148
 149 int dav1d_find_affine_int(const int (*pts)[2][2], const int np,
 150                           const int bw4, const int bh4,
 151                           const mv mv, Dav1dWarpedMotionParams *const wm,
 152                           const int bx4, const int by4)
 153 {
 154     int32_t *const mat = wm->matrix;
 155     int a[2][2] = { { 0, 0 }, { 0, 0 } };
 156     int bx[2] = { 0, 0 };
 157     int by[2] = { 0, 0 };
 158     const int rsuy = 2 * bh4 - 1;
 159     const int rsux = 2 * bw4 - 1;
 160     const int suy = rsuy * 8;
 161     const int sux = rsux * 8;
 162     const int duy = suy + mv.y;
 163     const int dux = sux + mv.x;
 164     const int isuy = by4 * 4 + rsuy;
 165     const int isux = bx4 * 4 + rsux;
 166
 167     for (int i = 0; i < np; i++) {
 168         const int dx = pts[i][1][0] - dux;
 169         const int dy = pts[i][1][1] - duy;
 170         const int sx = pts[i][0][0] - sux;
 171         const int sy = pts[i][0][1] - suy;
 172         if (abs(sx - dx) < 256 && abs(sy - dy) < 256) {
 173             a[0][0] += ((sx * sx) >> 2) + sx * 2 + 8;
 174             a[0][1] += ((sx * sy) >> 2) + sx + sy + 4;
 175             a[1][1] += ((sy * sy) >> 2) + sy * 2 + 8;
 176             bx[0] += ((sx * dx) >> 2) + sx + dx + 8;
 177             bx[1] += ((sy * dx) >> 2) + sy + dx + 4;
 178             by[0] += ((sx * dy) >> 2) + sx + dy + 4;
 179             by[1] += ((sy * dy) >> 2) + sy + dy + 8;
 180         }
 181     }
 182
 183     // compute determinant of a
 184     const int64_t det = (int64_t) a[0][0] * a[1][1] - (int64_t) a[0][1] * a[0][1];
 185     if (det == 0) return 1;
 186     int shift, idet = apply_sign64(resolve_divisor_64(llabs(det), &shift), det);
 187     shift -= 16;
 188     if (shift < 0) {
 189         idet <<= -shift;
 190         shift = 0;
 191     }
 192
 193     // solve the least-squares
 194     mat[2] = get_mult_shift_diag((int64_t) a[1][1] * bx[0] -
 195                                  (int64_t) a[0][1] * bx[1], idet, shift);
 196     mat[3] = get_mult_shift_ndiag((int64_t) a[0][0] * bx[1] -
 197                                   (int64_t) a[0][1] * bx[0], idet, shift);
 198     mat[4] = get_mult_shift_ndiag((int64_t) a[1][1] * by[0] -
 199                                   (int64_t) a[0][1] * by[1], idet, shift);
 200     mat[5] = get_mult_shift_diag((int64_t) a[0][0] * by[1] -
 201                                  (int64_t) a[0][1] * by[0], idet, shift);
 202
 203     mat[0] = iclip(mv.x * 0x2000 - (isux * (mat[2] - 0x10000) + isuy * mat[3]),
 204                    -0x800000, 0x7fffff);
 205     mat[1] = iclip(mv.y * 0x2000 - (isux * mat[4] + isuy * (mat[5] - 0x10000)),
 206                    -0x800000, 0x7fffff);
 207
 208     return 0;
 209 }