Home | History | Annotate | Download | only in aom_dsp 
      1 /*
      2  * Copyright (c) 2017, Alliance for Open Media. All rights reserved
      3  *
      4  * This source code is subject to the terms of the BSD 2 Clause License and
      5  * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
      6  * was not distributed with this source code in the LICENSE file, you can
      7  * obtain it at www.aomedia.org/license/software. If the Alliance for Open
      8  * Media Patent License 1.0 was not distributed with this source code in the
      9  * PATENTS file, you can obtain it at www.aomedia.org/license/patent.
     10  */
     11 
     12 #include <math.h>
     13 
     14 #include <stdio.h>
     15 #include <stdlib.h>
     16 #include <string.h>
     17 
     18 #include "aom_dsp/noise_util.h"
     19 #include "aom_dsp/fft_common.h"
     20 #include "aom_mem/aom_mem.h"
     21 #include "config/aom_dsp_rtcd.h"
     22 
     23 float aom_noise_psd_get_default_value(int block_size, float factor) {
     24   return (factor * factor / 10000) * block_size * block_size / 8;
     25 }
     26 
     27 // Internal representation of noise transform. It keeps track of the
     28 // transformed data and a temporary working buffer to use during the
     29 // transform.
     30 struct aom_noise_tx_t {
     31   float *tx_block;
     32   float *temp;
     33   int block_size;
     34   void (*fft)(const float *, float *, float *);
     35   void (*ifft)(const float *, float *, float *);
     36 };
     37 
     38 struct aom_noise_tx_t *aom_noise_tx_malloc(int block_size) {
     39   struct aom_noise_tx_t *noise_tx =
     40       (struct aom_noise_tx_t *)aom_malloc(sizeof(struct aom_noise_tx_t));
     41   if (!noise_tx) return NULL;
     42   memset(noise_tx, 0, sizeof(*noise_tx));
     43   switch (block_size) {
     44     case 2:
     45       noise_tx->fft = aom_fft2x2_float;
     46       noise_tx->ifft = aom_ifft2x2_float;
     47       break;
     48     case 4:
     49       noise_tx->fft = aom_fft4x4_float;
     50       noise_tx->ifft = aom_ifft4x4_float;
     51       break;
     52     case 8:
     53       noise_tx->fft = aom_fft8x8_float;
     54       noise_tx->ifft = aom_ifft8x8_float;
     55       break;
     56     case 16:
     57       noise_tx->fft = aom_fft16x16_float;
     58       noise_tx->ifft = aom_ifft16x16_float;
     59       break;
     60     case 32:
     61       noise_tx->fft = aom_fft32x32_float;
     62       noise_tx->ifft = aom_ifft32x32_float;
     63       break;
     64     default:
     65       aom_free(noise_tx);
     66       fprintf(stderr, "Unsupported block size %d\n", block_size);
     67       return NULL;
     68   }
     69   noise_tx->block_size = block_size;
     70   noise_tx->tx_block = (float *)aom_memalign(
     71       32, 2 * sizeof(*noise_tx->tx_block) * block_size * block_size);
     72   noise_tx->temp = (float *)aom_memalign(
     73       32, 2 * sizeof(*noise_tx->temp) * block_size * block_size);
     74   if (!noise_tx->tx_block || !noise_tx->temp) {
     75     aom_noise_tx_free(noise_tx);
     76     return NULL;
     77   }
     78   // Clear the buffers up front. Some outputs of the forward transform are
     79   // real only (the imaginary component will never be touched)
     80   memset(noise_tx->tx_block, 0,
     81          2 * sizeof(*noise_tx->tx_block) * block_size * block_size);
     82   memset(noise_tx->temp, 0,
     83          2 * sizeof(*noise_tx->temp) * block_size * block_size);
     84   return noise_tx;
     85 }
     86 
     87 void aom_noise_tx_forward(struct aom_noise_tx_t *noise_tx, const float *data) {
     88   noise_tx->fft(data, noise_tx->temp, noise_tx->tx_block);
     89 }
     90 
     91 void aom_noise_tx_filter(struct aom_noise_tx_t *noise_tx, const float *psd) {
     92   const int block_size = noise_tx->block_size;
     93   const float kBeta = 1.1f;
     94   const float kEps = 1e-6f;
     95   for (int y = 0; y < block_size; ++y) {
     96     for (int x = 0; x < block_size; ++x) {
     97       int i = y * block_size + x;
     98       float *c = noise_tx->tx_block + 2 * i;
     99       const float p = c[0] * c[0] + c[1] * c[1];
    100       if (p > kBeta * psd[i] && p > 1e-6) {
    101         noise_tx->tx_block[2 * i + 0] *= (p - psd[i]) / AOMMAX(p, kEps);
    102         noise_tx->tx_block[2 * i + 1] *= (p - psd[i]) / AOMMAX(p, kEps);
    103       } else {
    104         noise_tx->tx_block[2 * i + 0] *= (kBeta - 1.0f) / kBeta;
    105         noise_tx->tx_block[2 * i + 1] *= (kBeta - 1.0f) / kBeta;
    106       }
    107     }
    108   }
    109 }
    110 
    111 void aom_noise_tx_inverse(struct aom_noise_tx_t *noise_tx, float *data) {
    112   const int n = noise_tx->block_size * noise_tx->block_size;
    113   noise_tx->ifft(noise_tx->tx_block, noise_tx->temp, data);
    114   for (int i = 0; i < n; ++i) {
    115     data[i] /= n;
    116   }
    117 }
    118 
    119 void aom_noise_tx_add_energy(const struct aom_noise_tx_t *noise_tx,
    120                              float *psd) {
    121   const int block_size = noise_tx->block_size;
    122   for (int yb = 0; yb < block_size; ++yb) {
    123     for (int xb = 0; xb <= block_size / 2; ++xb) {
    124       float *c = noise_tx->tx_block + 2 * (yb * block_size + xb);
    125       psd[yb * block_size + xb] += c[0] * c[0] + c[1] * c[1];
    126     }
    127   }
    128 }
    129 
    130 void aom_noise_tx_free(struct aom_noise_tx_t *noise_tx) {
    131   if (!noise_tx) return;
    132   aom_free(noise_tx->tx_block);
    133   aom_free(noise_tx->temp);
    134   aom_free(noise_tx);
    135 }
    136 
    137 double aom_normalized_cross_correlation(const double *a, const double *b,
    138                                         int n) {
    139   double c = 0;
    140   double a_len = 0;
    141   double b_len = 0;
    142   for (int i = 0; i < n; ++i) {
    143     a_len += a[i] * a[i];
    144     b_len += b[i] * b[i];
    145     c += a[i] * b[i];
    146   }
    147   return c / (sqrt(a_len) * sqrt(b_len));
    148 }
    149 
    150 int aom_noise_data_validate(const double *data, int w, int h) {
    151   const double kVarianceThreshold = 2;
    152   const double kMeanThreshold = 2;
    153 
    154   int x = 0, y = 0;
    155   int ret_value = 1;
    156   double var = 0, mean = 0;
    157   double *mean_x, *mean_y, *var_x, *var_y;
    158 
    159   // Check that noise variance is not increasing in x or y
    160   // and that the data is zero mean.
    161   mean_x = (double *)aom_malloc(sizeof(*mean_x) * w);
    162   var_x = (double *)aom_malloc(sizeof(*var_x) * w);
    163   mean_y = (double *)aom_malloc(sizeof(*mean_x) * h);
    164   var_y = (double *)aom_malloc(sizeof(*var_y) * h);
    165 
    166   memset(mean_x, 0, sizeof(*mean_x) * w);
    167   memset(var_x, 0, sizeof(*var_x) * w);
    168   memset(mean_y, 0, sizeof(*mean_y) * h);
    169   memset(var_y, 0, sizeof(*var_y) * h);
    170 
    171   for (y = 0; y < h; ++y) {
    172     for (x = 0; x < w; ++x) {
    173       const double d = data[y * w + x];
    174       var_x[x] += d * d;
    175       var_y[y] += d * d;
    176       mean_x[x] += d;
    177       mean_y[y] += d;
    178       var += d * d;
    179       mean += d;
    180     }
    181   }
    182   mean /= (w * h);
    183   var = var / (w * h) - mean * mean;
    184 
    185   for (y = 0; y < h; ++y) {
    186     mean_y[y] /= h;
    187     var_y[y] = var_y[y] / h - mean_y[y] * mean_y[y];
    188     if (fabs(var_y[y] - var) >= kVarianceThreshold) {
    189       fprintf(stderr, "Variance distance too large %f %f\n", var_y[y], var);
    190       ret_value = 0;
    191       break;
    192     }
    193     if (fabs(mean_y[y] - mean) >= kMeanThreshold) {
    194       fprintf(stderr, "Mean distance too large %f %f\n", mean_y[y], mean);
    195       ret_value = 0;
    196       break;
    197     }
    198   }
    199 
    200   for (x = 0; x < w; ++x) {
    201     mean_x[x] /= w;
    202     var_x[x] = var_x[x] / w - mean_x[x] * mean_x[x];
    203     if (fabs(var_x[x] - var) >= kVarianceThreshold) {
    204       fprintf(stderr, "Variance distance too large %f %f\n", var_x[x], var);
    205       ret_value = 0;
    206       break;
    207     }
    208     if (fabs(mean_x[x] - mean) >= kMeanThreshold) {
    209       fprintf(stderr, "Mean distance too large %f %f\n", mean_x[x], mean);
    210       ret_value = 0;
    211       break;
    212     }
    213   }
    214 
    215   aom_free(mean_x);
    216   aom_free(mean_y);
    217   aom_free(var_x);
    218   aom_free(var_y);
    219 
    220   return ret_value;
    221 }
    222