in libvmaf/src/feature/integer_adm.c [789:930]
static void adm_decouple_s123(AdmBuffer *buf, int w, int h, int stride,
double adm_enhn_gain_limit)
{
const float cos_1deg_sq = cos(1.0 * M_PI / 180.0) * cos(1.0 * M_PI / 180.0);
const i4_adm_dwt_band_t *ref = &buf->i4_ref_dwt2;
const i4_adm_dwt_band_t *dis = &buf->i4_dis_dwt2;
const i4_adm_dwt_band_t *r = &buf->i4_decouple_r;
const i4_adm_dwt_band_t *a = &buf->i4_decouple_a;
/* The computation of the score is not required for the regions
which lie outside the frame borders */
int left = w * ADM_BORDER_FACTOR - 0.5 - 1; // -1 for filter tap
int top = h * ADM_BORDER_FACTOR - 0.5 - 1;
int right = w - left + 2; // +2 for filter tap
int bottom = h - top + 2;
if (left < 0) {
left = 0;
}
if (right > w) {
right = w;
}
if (top < 0) {
top = 0;
}
if (bottom > h) {
bottom = h;
}
int64_t ot_dp, o_mag_sq, t_mag_sq;
for (int i = top; i < bottom; ++i)
{
for (int j = left; j < right; ++j)
{
int32_t oh = ref->band_h[i * stride + j];
int32_t ov = ref->band_v[i * stride + j];
int32_t od = ref->band_d[i * stride + j];
int32_t th = dis->band_h[i * stride + j];
int32_t tv = dis->band_v[i * stride + j];
int32_t td = dis->band_d[i * stride + j];
int32_t rst_h, rst_v, rst_d;
/* Determine if angle between (oh,ov) and (th,tv) is less than 1 degree.
* Given that u is the angle (oh,ov) and v is the angle (th,tv), this can
* be done by testing the inequvality.
*
* { (u.v.) >= 0 } AND { (u.v)^2 >= cos(1deg)^2 * ||u||^2 * ||v||^2 }
*
* Proof:
*
* cos(theta) = (u.v) / (||u|| * ||v||)
*
* IF u.v >= 0 THEN
* cos(theta)^2 = (u.v)^2 / (||u||^2 * ||v||^2)
* (u.v)^2 = cos(theta)^2 * ||u||^2 * ||v||^2
*
* IF |theta| < 1deg THEN
* (u.v)^2 >= cos(1deg)^2 * ||u||^2 * ||v||^2
* END
* ELSE
* |theta| > 90deg
* END
*/
ot_dp = (int64_t)oh * th + (int64_t)ov * tv;
o_mag_sq = (int64_t)oh * oh + (int64_t)ov * ov;
t_mag_sq = (int64_t)th * th + (int64_t)tv * tv;
int angle_flag = (((float)ot_dp / 4096.0) >= 0.0f) &&
(((float)ot_dp / 4096.0) * ((float)ot_dp / 4096.0) >=
cos_1deg_sq * ((float)o_mag_sq / 4096.0) * ((float)t_mag_sq / 4096.0));
/**
* Division th/oh is carried using lookup table and converted to multiplication
* int64 / int32 is converted to multiplication using following method
* num /den :
* DenAbs = Abs(den)
* MSBDen = MSB(DenAbs) (gives position of first 1 bit form msb side)
* If (DenAbs < (1 << 15))
* Round = (1<<14)
* Score = (num * div_lookup[den] + Round ) >> 15
* else
* RoundD = (1<< (16 - MSBDen))
* Round = (1<< (14 + (17 - MSBDen))
* Score = (num * div_lookup[(DenAbs + RoundD )>>(17 - MSBDen)]*sign(Denominator) + Round)
* >> ((15 + (17 - MSBDen))
*/
int32_t kh_shift = 0;
int32_t kv_shift = 0;
int32_t kd_shift = 0;
uint32_t abs_oh = abs(oh);
uint32_t abs_ov = abs(ov);
uint32_t abs_od = abs(od);
int8_t kh_sign = (oh < 0 ? -1 : 1);
int8_t kv_sign = (ov < 0 ? -1 : 1);
int8_t kd_sign = (od < 0 ? -1 : 1);
uint16_t kh_msb = (abs_oh < (32768) ? abs_oh : get_best15_from32(abs_oh, &kh_shift));
uint16_t kv_msb = (abs_ov < (32768) ? abs_ov : get_best15_from32(abs_ov, &kv_shift));
uint16_t kd_msb = (abs_od < (32768) ? abs_od : get_best15_from32(abs_od, &kd_shift));
int64_t tmp_kh = (oh == 0) ? 32768 : (((int64_t)div_lookup[kh_msb + 32768] * th)*(kh_sign) +
(1 << (14 + kh_shift))) >> (15 + kh_shift);
int64_t tmp_kv = (ov == 0) ? 32768 : (((int64_t)div_lookup[kv_msb + 32768] * tv)*(kv_sign) +
(1 << (14 + kv_shift))) >> (15 + kv_shift);
int64_t tmp_kd = (od == 0) ? 32768 : (((int64_t)div_lookup[kd_msb + 32768] * td)*(kd_sign) +
(1 << (14 + kd_shift))) >> (15 + kd_shift);
int64_t kh = tmp_kh < 0 ? 0 : (tmp_kh > 32768 ? 32768 : tmp_kh);
int64_t kv = tmp_kv < 0 ? 0 : (tmp_kv > 32768 ? 32768 : tmp_kv);
int64_t kd = tmp_kd < 0 ? 0 : (tmp_kd > 32768 ? 32768 : tmp_kd);
rst_h = ((kh * oh) + 16384) >> 15;
rst_v = ((kv * ov) + 16384) >> 15;
rst_d = ((kd * od) + 16384) >> 15;
const float rst_h_f = ((float)kh / 32768) * ((float)oh / 64);
const float rst_v_f = ((float)kv / 32768) * ((float)ov / 64);
const float rst_d_f = ((float)kd / 32768) * ((float)od / 64);
if (angle_flag && (rst_h_f > 0.)) rst_h = MIN((rst_h * adm_enhn_gain_limit), th);
if (angle_flag && (rst_h_f < 0.)) rst_h = MAX((rst_h * adm_enhn_gain_limit), th);
if (angle_flag && (rst_v_f > 0.)) rst_v = MIN(rst_v * adm_enhn_gain_limit, tv);
if (angle_flag && (rst_v_f < 0.)) rst_v = MAX(rst_v * adm_enhn_gain_limit, tv);
if (angle_flag && (rst_d_f > 0.)) rst_d = MIN(rst_d * adm_enhn_gain_limit, td);
if (angle_flag && (rst_d_f < 0.)) rst_d = MAX(rst_d * adm_enhn_gain_limit, td);
r->band_h[i * stride + j] = rst_h;
r->band_v[i * stride + j] = rst_v;
r->band_d[i * stride + j] = rst_d;
a->band_h[i * stride + j] = th - rst_h;
a->band_v[i * stride + j] = tv - rst_v;
a->band_d[i * stride + j] = td - rst_d;
}
}
}