adm_tools.h

Go to the documentation of this file.

 
#include <math.h>
#include "common/macros.h"
#include "barten_csf_tools.h"
#include "adm_csf_tools.h"
 
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327
#endif
 
#pragma once
 
#ifndef ADM_TOOLS_H_
#define ADM_TOOLS_H_
 
// i = 0, j = 0: indices y: 1,0,1, x: 1,0,1
#define ADM_CM_THRESH_S_0_0(angles,flt_angles,src_px_stride,accum,w,h,i,j) \
{ \
    *accum = 0; \
    for (int theta = 0; theta < 3; ++theta) { \
        const float *src_ptr = angles[theta]; \
        const float *flt_ptr = flt_angles[theta]; \
        float sum = 0; \
        sum += flt_ptr[src_px_stride + 1]; \
        sum += flt_ptr[src_px_stride]; \
        sum += flt_ptr[src_px_stride + 1]; \
        sum += flt_ptr[1]; \
        sum += FLOAT_ONE_BY_15 * fabsf(src_ptr[0]); \
        sum += flt_ptr[1]; \
        sum += flt_ptr[src_px_stride + 1]; \
        sum += flt_ptr[src_px_stride]; \
        sum += flt_ptr[src_px_stride + 1]; \
        *accum += sum; \
    } \
}
 
// i = 0, j = w-1: indices y: 1,0,1, x: w-2, w-1, w-1 
#define ADM_CM_THRESH_S_0_W_M_1(angles,flt_angles,src_px_stride,accum,w,h,i,j) \
{ \
    *accum = 0; \
    for (int theta = 0; theta < 3; ++theta) { \
        const float *src_ptr = angles[theta]; \
        const float *flt_ptr = flt_angles[theta]; \
        float sum = 0; \
        sum += flt_ptr[src_px_stride + w - 2]; \
        sum += flt_ptr[src_px_stride + w - 1]; \
        sum += flt_ptr[src_px_stride + w - 1]; \
        sum += flt_ptr[w - 2]; \
        sum += FLOAT_ONE_BY_15 * fabsf(src_ptr[w - 1]); \
        sum += flt_ptr[w - 1]; \
        sum += flt_ptr[src_px_stride + w - 2]; \
        sum += flt_ptr[src_px_stride + w - 1]; \
        sum += flt_ptr[src_px_stride + w - 1]; \
        *accum += sum; \
    } \
}
 
// i = 0, j = 1, ..., w-2: indices y: 1,0,1, x: j-1,j,j+1 
#define ADM_CM_THRESH_S_0_J(angles,flt_angles,src_px_stride,accum,w,h,i,j) \
{ \
    *accum = 0; \
    for (int theta = 0; theta < 3; ++theta) { \
        const float *src_ptr = angles[theta]; \
        const float *flt_ptr = flt_angles[theta]; \
        float sum = 0; \
        sum += flt_ptr[src_px_stride + j - 1]; \
        sum += flt_ptr[src_px_stride + j]; \
        sum += flt_ptr[src_px_stride + j + 1]; \
        sum += flt_ptr[j - 1]; \
        sum += FLOAT_ONE_BY_15 * fabsf(src_ptr[j]); \
        sum += flt_ptr[j + 1]; \
        sum += flt_ptr[src_px_stride + j - 1]; \
        sum += flt_ptr[src_px_stride + j]; \
        sum += flt_ptr[src_px_stride + j + 1];  \
        *accum += sum; \
    } \
}
 
// i = h-1, j = 0: indices y: h-2,h-1,h-1, x: 1,0,1 
#define ADM_CM_THRESH_S_H_M_1_0(angles,flt_angles,src_px_stride,accum,w,h,i,j) \
{ \
    *accum = 0; \
    for (int theta = 0; theta < 3; ++theta) { \
        const float *src_ptr = angles[theta]; \
        const float *flt_ptr = flt_angles[theta]; \
        float sum = 0; \
        src_ptr += (src_px_stride * (h - 2)); \
        flt_ptr += (src_px_stride * (h - 2)); \
        sum += flt_ptr[1]; \
        sum += flt_ptr[0]; \
        sum += flt_ptr[1]; \
        src_ptr += src_px_stride; \
        flt_ptr += src_px_stride; \
        sum += flt_ptr[1]; \
        sum += FLOAT_ONE_BY_15 * fabsf(src_ptr[0]); \
        sum += flt_ptr[1]; \
        sum += flt_ptr[1]; \
        sum += flt_ptr[0]; \
        sum += flt_ptr[1]; \
        *accum += sum; \
    } \
}
 
// i = h-1, j = w-1: indices y: h-2,h-1,h-1, x: w-2, w-1, w-1 
#define ADM_CM_THRESH_S_H_M_1_W_M_1(angles,flt_angles,src_px_stride,accum,w,h,i,j) \
{ \
    *accum = 0; \
    for (int theta = 0; theta < 3; ++theta) { \
        const float *src_ptr = angles[theta]; \
        const float *flt_ptr = flt_angles[theta]; \
        float sum = 0; \
        src_ptr += (src_px_stride * (h - 2)); \
        flt_ptr += (src_px_stride * (h - 2)); \
        sum += flt_ptr[w - 2]; \
        sum += flt_ptr[w - 1]; \
        sum += flt_ptr[w - 1]; \
        src_ptr += src_px_stride; \
        flt_ptr += src_px_stride; \
        sum += flt_ptr[w - 2]; \
        sum += FLOAT_ONE_BY_15 * fabsf(src_ptr[w - 1]); \
        sum += flt_ptr[w - 1]; \
        sum += flt_ptr[w - 2]; \
        sum += flt_ptr[w - 1]; \
        sum += flt_ptr[w - 1]; \
        *accum += sum; \
    } \
}
 
// i = h-1, j = 1, ..., w-2: indices y: h-2,h-1,h-1, x: j-1,j,j+1 
#define ADM_CM_THRESH_S_H_M_1_J(angles,flt_angles,src_px_stride,accum,w,h,i,j) \
{ \
    *accum = 0; \
    for (int theta = 0; theta < 3; ++theta) { \
        const float *src_ptr = angles[theta]; \
        const float *flt_ptr = flt_angles[theta]; \
        float sum = 0; \
        src_ptr += (src_px_stride * (h - 2)); \
        flt_ptr += (src_px_stride * (h - 2)); \
        sum += flt_ptr[j - 1];\
        sum += flt_ptr[j]; \
        sum += flt_ptr[j + 1]; \
        src_ptr += src_px_stride; \
        flt_ptr += src_px_stride; \
        sum += flt_ptr[j - 1]; \
        sum += FLOAT_ONE_BY_15 * fabsf(src_ptr[j]); \
        sum += flt_ptr[j + 1]; \
        sum += flt_ptr[j - 1]; \
        sum += flt_ptr[j]; \
        sum += flt_ptr[j + 1]; \
        *accum += sum; \
    } \
}
 
// i = 1,..,h-2, j = 1,..,w-2: indices y: i-1,i,i+1, x: j-1,j,j+1 
#define ADM_CM_THRESH_S_I_J(angles,flt_angles,src_px_stride,accum,w,h,i,j) \
{ \
    *accum = 0; \
    for (int theta = 0; theta < 3; ++theta) { \
        const float *src_ptr = angles[theta]; \
        const float *flt_ptr = flt_angles[theta]; \
        float sum = 0; \
        src_ptr += (src_px_stride * (i - 1)); \
        flt_ptr += (src_px_stride * (i - 1)); \
        sum += flt_ptr[j - 1]; \
        sum += flt_ptr[j]; \
        sum += flt_ptr[j + 1]; \
        src_ptr += src_px_stride; \
        flt_ptr += src_px_stride; \
        sum += flt_ptr[j - 1]; \
        sum += FLOAT_ONE_BY_15 * fabsf(src_ptr[j]); \
        sum += flt_ptr[j + 1]; \
        src_ptr += src_px_stride; \
        flt_ptr += src_px_stride; \
        sum += flt_ptr[j - 1]; \
        sum += flt_ptr[j]; \
        sum += flt_ptr[j + 1]; \
        *accum += sum; \
    } \
}
 
// i = 1,..,h-2, j = 0: indices y: i-1,i,i+1, x: 1,0,1 
#define ADM_CM_THRESH_S_I_0(angles,flt_angles,src_px_stride,accum,w,h,i,j) \
{ \
    *accum = 0; \
    for (int theta = 0; theta < 3; ++theta) { \
        const float *src_ptr = angles[theta]; \
        const float *flt_ptr = flt_angles[theta]; \
        float sum = 0; \
        src_ptr += (src_px_stride * (i - 1)); \
        flt_ptr += (src_px_stride * (i - 1)); \
        sum += flt_ptr[1]; \
        sum += flt_ptr[0]; \
        sum += flt_ptr[1]; \
        src_ptr += src_px_stride; \
        flt_ptr += src_px_stride; \
        sum += flt_ptr[1]; \
        sum += FLOAT_ONE_BY_15 * fabsf(src_ptr[0]); \
        sum += flt_ptr[1]; \
        src_ptr += src_px_stride; \
        flt_ptr += src_px_stride; \
        sum += flt_ptr[1]; \
        sum += flt_ptr[0]; \
        sum += flt_ptr[1]; \
        *accum += sum; \
    } \
}
 
// i = 1,..,h-2, j = w-1: indices y: i-1,i,i+1, x: w-2,w-1,w-1 
#define ADM_CM_THRESH_S_I_W_M_1(angles,flt_angles,src_px_stride,accum,w,h,i,j) \
{ \
    *accum = 0; \
    for (int theta = 0; theta < 3; ++theta) { \
        const float *src_ptr = angles[theta]; \
        const float *flt_ptr = flt_angles[theta]; \
        float sum = 0; \
        src_ptr += (src_px_stride * (i-1)); \
        flt_ptr += (src_px_stride * (i - 1)); \
        sum += flt_ptr[w - 2]; \
        sum += flt_ptr[w - 1]; \
        sum += flt_ptr[w - 1]; \
        src_ptr += src_px_stride; \
        flt_ptr += src_px_stride; \
        sum += flt_ptr[w - 2]; \
        sum += FLOAT_ONE_BY_15 * fabsf(src_ptr[w - 1]); \
        sum += flt_ptr[w - 1]; \
        src_ptr += src_px_stride; \
        flt_ptr += src_px_stride; \
        sum += flt_ptr[w - 2]; \
        sum += flt_ptr[w - 1]; \
        sum += flt_ptr[w - 1]; \
        *accum += sum; \
    } \
}
 
typedef struct adm_dwt_band_t_s {
    float *band_a; /* Low-pass V + low-pass H. */
    float *band_v; /* Low-pass V + high-pass H. */
    float *band_h; /* High-pass V + low-pass H. */
    float *band_d; /* High-pass V + high-pass H. */
} adm_dwt_band_t_s;
 
typedef struct adm_dwt_band_t_d {
    double *band_a; /* Low-pass V + low-pass H. */
    double *band_v; /* Low-pass V + high-pass H. */
    double *band_h; /* High-pass V + low-pass H. */
    double *band_d; /* High-pass V + high-pass H. */
} adm_dwt_band_t_d;
 
float adm_sum_cube_s(const float *x, int w, int h, int stride, double border_factor, double adm_p_norm);
 
void adm_decouple_s(const adm_dwt_band_t_s *ref, const adm_dwt_band_t_s *dis, const adm_dwt_band_t_s *r, const adm_dwt_band_t_s *a, int w, int h, int ref_stride, int dis_stride, int r_stride, int a_stride, double border_factor, double adm_enhn_gain_limit);
 
void adm_csf_s(const adm_dwt_band_t_s *src, const adm_dwt_band_t_s *dst, const adm_dwt_band_t_s *flt, int orig_h, int scale, int w, int h, int src_stride, int dst_stride, double border_factor, double adm_norm_view_dist, int adm_ref_display_height, int adm_csf_mode, double luminance_level, double adm_csf_scale, double adm_csf_diag_scale, double adm_f1s0, double adm_f1s1, double adm_f1s2, double adm_f1s3, double adm_f2s0, double adm_f2s1, double adm_f2s2, double adm_f2s3);
 
void adm_cm_thresh_s(const adm_dwt_band_t_s *src, float *dst, int w, int h, int src_stride, int dst_stride);
 
float adm_csf_den_scale_s(const adm_dwt_band_t_s *src, int orig_h, int scale, int w, int h, int src_stride, double border_factor, double adm_norm_view_dist, int adm_ref_display_height, int adm_csf_mode, double luminance_level, double adm_csf_scale, double adm_csf_diag_scale, double adm_noise_weight, double adm_p_norm, double adm_f1s0, double adm_f1s1, double adm_f1s2, double adm_f1s3, double adm_f2s0, double adm_f2s1, double adm_f2s2, double adm_f2s3);
 
float adm_cm_s(const adm_dwt_band_t_s *src, const adm_dwt_band_t_s *dst, const adm_dwt_band_t_s *csf_a, int w, int h, int src_stride, int dst_stride, int csf_a_stride, double border_factor, int scale, double adm_norm_view_dist, int adm_ref_display_height, int adm_csf_mode, double luminance_level, double adm_csf_scale, double adm_csf_diag_scale, double adm_noise_weight, int adm_bypass_cm, double adm_p_norm, double adm_f1s0, double adm_f1s1, double adm_f1s2, double adm_f1s3, double adm_f2s0, double adm_f2s1, double adm_f2s2, double adm_f2s3);
 
void dwt2_src_indices_filt_s(int **src_ind_y, int **src_ind_x, int w, int h);
 
void adm_dwt2_s(const float *src, const adm_dwt_band_t_s *dst, int **ind_y, int **ind_x, int w, int h, int src_stride, int dst_stride);
 
void adm_dwt2_lo_s(const float *src, const adm_dwt_band_t_s *dst, int **ind_y, int **ind_x, int w, int h, int src_stride, int dst_stride);
 
void adm_dwt2_d(const double *src, const adm_dwt_band_t_d *dst, int **ind_y, int **ind_x, int w, int h, int src_stride, int dst_stride);
 
/* ================= */
/* Noise floor model */
/* ================= */
 
/*
 * The following dwt visibility threshold parameters are taken from
 * "Visibility of Wavelet Quantization Noise"
 * by A. B. Watson, G. Y. Yang, J. A. Solomon and J. Villasenor
 * IEEE Trans. on Image Processing, Vol. 6, No 8, Aug. 1997
 * Page 1170, formula (7) and corresponding Table IV
 * Table IV has 2 entries for Cb and Cr thresholds
 * Chose those corresponding to subject "sfl" since they are lower
 * These thresholds were obtained and modeled for the 7-9 biorthogonal wavelet basis
 */
struct dwt_model_params {
    float a;
    float k;
    float f0;
    float g[4];
};
 
// 0 -> Y, 1 -> Cb, 2 -> Cr
static const struct dwt_model_params dwt_7_9_YCbCr_threshold[3] = {
    { .a = 0.495, .k = 0.466, .f0 = 0.401, .g = { 1.501, 1.0, 0.534, 1.0} },
    { .a = 1.633, .k = 0.353, .f0 = 0.209, .g = { 1.520, 1.0, 0.502, 1.0} },
    { .a = 0.944, .k = 0.521, .f0 = 0.404, .g = { 1.868, 1.0, 0.516, 1.0} }
};
 
/*
 * The following dwt basis function amplitudes, A(lambda,theta), are taken from
 * "Visibility of Wavelet Quantization Noise"
 * by A. B. Watson, G. Y. Yang, J. A. Solomon and J. Villasenor
 * IEEE Trans. on Image Processing, Vol. 6, No 8, Aug. 1997
 * Page 1172, Table V
 * The table has been transposed, i.e. it can be used directly to obtain A[lambda][theta]
 * These amplitudes were calculated for the 7-9 biorthogonal wavelet basis
 */
static const float dwt_7_9_basis_function_amplitudes[6][4] = {
    { 0.62171,  0.67234,  0.72709,  0.67234  },
    { 0.34537,  0.41317,  0.49428,  0.41317  },
    { 0.18004,  0.22727,  0.28688,  0.22727  },
    { 0.091401, 0.11792,  0.15214,  0.11792  },
    { 0.045943, 0.059758, 0.077727, 0.059758 },
    { 0.023013, 0.030018, 0.039156, 0.030018 }
};
 
/*
 * lambda = 0 (finest scale), 1, 2, 3 (coarsest scale);
 * theta = 0 (ll), 1 (lh - vertical), 2 (hh - diagonal), 3(hl - horizontal).
 */
static FORCE_INLINE float dwt_quant_step(const struct dwt_model_params *params,
        int lambda, int theta, double adm_norm_view_dist, int adm_ref_display_height)
{
    // Formula (1), page 1165 - display visual resolution (DVR), in pixels/degree of visual angle. This should be 56.55
    float r = adm_norm_view_dist * adm_ref_display_height * M_PI / 180.0;
 
    // Formula (9), page 1171
    float temp = log10(pow(2.0,lambda+1)*params->f0*params->g[theta]/r);
    float Q = 2.0*params->a*pow(10.0,params->k*temp*temp)/dwt_7_9_basis_function_amplitudes[lambda][theta];
 
    return Q;
}
 
#endif /* ADM_TOOLS_H_ */