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
aom_noise_psd_get_default_value(int block_size,float factor)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
aom_noise_tx_malloc(int block_size)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
aom_noise_tx_forward(struct aom_noise_tx_t * noise_tx,const float * data)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
aom_noise_tx_filter(struct aom_noise_tx_t * noise_tx,const float * psd)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
aom_noise_tx_inverse(struct aom_noise_tx_t * noise_tx,float * data)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
aom_noise_tx_add_energy(const struct aom_noise_tx_t * noise_tx,float * psd)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
aom_noise_tx_free(struct aom_noise_tx_t * noise_tx)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
aom_normalized_cross_correlation(const double * a,const double * b,int n)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
aom_noise_data_validate(const double * data,int w,int h)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_malloc(sizeof(*mean_x) * w);
164 var_x = (double *)aom_malloc(sizeof(*var_x) * w);
165 mean_y = (double *)aom_malloc(sizeof(*mean_x) * h);
166 var_y = (double *)aom_malloc(sizeof(*var_y) * h);
167
168 memset(mean_x, 0, sizeof(*mean_x) * w);
169 memset(var_x, 0, sizeof(*var_x) * w);
170 memset(mean_y, 0, sizeof(*mean_y) * h);
171 memset(var_y, 0, sizeof(*var_y) * h);
172
173 for (y = 0; y < h; ++y) {
174 for (x = 0; x < w; ++x) {
175 const double d = data[y * w + x];
176 var_x[x] += d * d;
177 var_y[y] += d * d;
178 mean_x[x] += d;
179 mean_y[y] += d;
180 var += d * d;
181 mean += d;
182 }
183 }
184 mean /= (w * h);
185 var = var / (w * h) - mean * mean;
186
187 for (y = 0; y < h; ++y) {
188 mean_y[y] /= h;
189 var_y[y] = var_y[y] / h - mean_y[y] * mean_y[y];
190 if (fabs(var_y[y] - var) >= kVarianceThreshold) {
191 fprintf(stderr, "Variance distance too large %f %f\n", var_y[y], var);
192 ret_value = 0;
193 break;
194 }
195 if (fabs(mean_y[y] - mean) >= kMeanThreshold) {
196 fprintf(stderr, "Mean distance too large %f %f\n", mean_y[y], mean);
197 ret_value = 0;
198 break;
199 }
200 }
201
202 for (x = 0; x < w; ++x) {
203 mean_x[x] /= w;
204 var_x[x] = var_x[x] / w - mean_x[x] * mean_x[x];
205 if (fabs(var_x[x] - var) >= kVarianceThreshold) {
206 fprintf(stderr, "Variance distance too large %f %f\n", var_x[x], var);
207 ret_value = 0;
208 break;
209 }
210 if (fabs(mean_x[x] - mean) >= kMeanThreshold) {
211 fprintf(stderr, "Mean distance too large %f %f\n", mean_x[x], mean);
212 ret_value = 0;
213 break;
214 }
215 }
216
217 aom_free(mean_x);
218 aom_free(mean_y);
219 aom_free(var_x);
220 aom_free(var_y);
221
222 return ret_value;
223 }
224