in libvmaf/src/feature/integer_adm.c [660:778]
static void adm_decouple(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 adm_dwt_band_t *ref = &buf->ref_dwt2;
const adm_dwt_band_t *dis = &buf->dis_dwt2;
const adm_dwt_band_t *r = &buf->decouple_r;
const adm_dwt_band_t *a = &buf->decouple_a;
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) {
int16_t oh = ref->band_h[i * stride + j];
int16_t ov = ref->band_v[i * stride + j];
int16_t od = ref->band_d[i * stride + j];
int16_t th = dis->band_h[i * stride + j];
int16_t tv = dis->band_v[i * stride + j];
int16_t td = dis->band_d[i * stride + j];
int16_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;
/**
* angle_flag is calculated in floating-point by converting fixed-point variables back to
* floating-point
*/
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
*/
int32_t tmp_kh = (oh == 0) ?
32768 : (((int64_t)div_lookup[oh + 32768] * th) + 16384) >> 15;
int32_t tmp_kv = (ov == 0) ?
32768 : (((int64_t)div_lookup[ov + 32768] * tv) + 16384) >> 15;
int32_t tmp_kd = (od == 0) ?
32768 : (((int64_t)div_lookup[od + 32768] * td) + 16384) >> 15;
int32_t kh = tmp_kh < 0 ? 0 : (tmp_kh > 32768 ? 32768 : tmp_kh);
int32_t kv = tmp_kv < 0 ? 0 : (tmp_kv > 32768 ? 32768 : tmp_kv);
int32_t kd = tmp_kd < 0 ? 0 : (tmp_kd > 32768 ? 32768 : tmp_kd);
/**
* kh,kv,kd are in Q15 type and oh,ov,od are in Q16 type hence shifted by
* 15 to make result Q16
*/
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;
}
}
}