warpmv.c (9052B)
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 }