noise_util.c (7000B)
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 c0 = AOMMAX((float)fabs(c[0]), 1e-8f); 100 const float c1 = AOMMAX((float)fabs(c[1]), 1e-8f); 101 const float p = c0 * c0 + c1 * c1; 102 if (p > kBeta * psd[i] && p > 1e-6) { 103 noise_tx->tx_block[2 * i + 0] *= (p - psd[i]) / AOMMAX(p, kEps); 104 noise_tx->tx_block[2 * i + 1] *= (p - psd[i]) / AOMMAX(p, kEps); 105 } else { 106 noise_tx->tx_block[2 * i + 0] *= (kBeta - 1.0f) / kBeta; 107 noise_tx->tx_block[2 * i + 1] *= (kBeta - 1.0f) / kBeta; 108 } 109 } 110 } 111 } 112 113 void aom_noise_tx_inverse(struct aom_noise_tx_t *noise_tx, float *data) { 114 const int n = noise_tx->block_size * noise_tx->block_size; 115 noise_tx->ifft(noise_tx->tx_block, noise_tx->temp, data); 116 for (int i = 0; i < n; ++i) { 117 data[i] /= n; 118 } 119 } 120 121 void aom_noise_tx_add_energy(const struct aom_noise_tx_t *noise_tx, 122 float *psd) { 123 const int block_size = noise_tx->block_size; 124 for (int yb = 0; yb < block_size; ++yb) { 125 for (int xb = 0; xb <= block_size / 2; ++xb) { 126 float *c = noise_tx->tx_block + 2 * (yb * block_size + xb); 127 psd[yb * block_size + xb] += c[0] * c[0] + c[1] * c[1]; 128 } 129 } 130 } 131 132 void aom_noise_tx_free(struct aom_noise_tx_t *noise_tx) { 133 if (!noise_tx) return; 134 aom_free(noise_tx->tx_block); 135 aom_free(noise_tx->temp); 136 aom_free(noise_tx); 137 } 138 139 double aom_normalized_cross_correlation(const double *a, const double *b, 140 int n) { 141 double c = 0; 142 double a_len = 0; 143 double b_len = 0; 144 for (int i = 0; i < n; ++i) { 145 a_len += a[i] * a[i]; 146 b_len += b[i] * b[i]; 147 c += a[i] * b[i]; 148 } 149 return c / (sqrt(a_len) * sqrt(b_len)); 150 } 151 152 int aom_noise_data_validate(const double *data, int w, int h) { 153 const double kVarianceThreshold = 2; 154 const double kMeanThreshold = 2; 155 156 int x = 0, y = 0; 157 int ret_value = 1; 158 double var = 0, mean = 0; 159 double *mean_x, *mean_y, *var_x, *var_y; 160 161 // Check that noise variance is not increasing in x or y 162 // and that the data is zero mean. 163 mean_x = (double *)aom_calloc(w, sizeof(*mean_x)); 164 var_x = (double *)aom_calloc(w, sizeof(*var_x)); 165 mean_y = (double *)aom_calloc(h, sizeof(*mean_x)); 166 var_y = (double *)aom_calloc(h, sizeof(*var_y)); 167 if (!(mean_x && var_x && mean_y && var_y)) { 168 aom_free(mean_x); 169 aom_free(mean_y); 170 aom_free(var_x); 171 aom_free(var_y); 172 return 0; 173 } 174 175 for (y = 0; y < h; ++y) { 176 for (x = 0; x < w; ++x) { 177 const double d = data[y * w + x]; 178 var_x[x] += d * d; 179 var_y[y] += d * d; 180 mean_x[x] += d; 181 mean_y[y] += d; 182 var += d * d; 183 mean += d; 184 } 185 } 186 mean /= (w * h); 187 var = var / (w * h) - mean * mean; 188 189 for (y = 0; y < h; ++y) { 190 mean_y[y] /= h; 191 var_y[y] = var_y[y] / h - mean_y[y] * mean_y[y]; 192 if (fabs(var_y[y] - var) >= kVarianceThreshold) { 193 fprintf(stderr, "Variance distance too large %f %f\n", var_y[y], var); 194 ret_value = 0; 195 break; 196 } 197 if (fabs(mean_y[y] - mean) >= kMeanThreshold) { 198 fprintf(stderr, "Mean distance too large %f %f\n", mean_y[y], mean); 199 ret_value = 0; 200 break; 201 } 202 } 203 204 for (x = 0; x < w; ++x) { 205 mean_x[x] /= w; 206 var_x[x] = var_x[x] / w - mean_x[x] * mean_x[x]; 207 if (fabs(var_x[x] - var) >= kVarianceThreshold) { 208 fprintf(stderr, "Variance distance too large %f %f\n", var_x[x], var); 209 ret_value = 0; 210 break; 211 } 212 if (fabs(mean_x[x] - mean) >= kMeanThreshold) { 213 fprintf(stderr, "Mean distance too large %f %f\n", mean_x[x], mean); 214 ret_value = 0; 215 break; 216 } 217 } 218 219 aom_free(mean_x); 220 aom_free(mean_y); 221 aom_free(var_x); 222 aom_free(var_y); 223 224 return ret_value; 225 }