Add utilities to aom_dsp for modeling correlated noise.
The auto-regressive model allows for different window shapes
and different lag sizes.
Although most likely to be used as a reference for modeling
noise in AV1, the model is currently parameterized more generally
than AV1 needs.
I will add an example (hopefully with a denoiser) in future
commits.
Change-Id: I1ba1067543601c2c01db4970d42766bb35da77f0
diff --git a/aom_dsp/noise_util.c b/aom_dsp/noise_util.c
new file mode 100644
index 0000000..d1a0ccf
--- /dev/null
+++ b/aom_dsp/noise_util.c
@@ -0,0 +1,149 @@
+/*
+ * Copyright (c) 2017, Alliance for Open Media. All rights reserved
+ *
+ * This source code is subject to the terms of the BSD 2 Clause License and
+ * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
+ * was not distributed with this source code in the LICENSE file, you can
+ * obtain it at www.aomedia.org/license/software. If the Alliance for Open
+ * Media Patent License 1.0 was not distributed with this source code in the
+ * PATENTS file, you can obtain it at www.aomedia.org/license/patent.
+ */
+
+#include <math.h>
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "aom_dsp/noise_util.h"
+#include "aom_mem/aom_mem.h"
+
+// Return normally distrbuted values with standard deviation of sigma.
+double aom_randn(double sigma) {
+ while (1) {
+ const double u = 2.0 * ((double)rand()) / RAND_MAX - 1.0;
+ const double v = 2.0 * ((double)rand()) / RAND_MAX - 1.0;
+ const double s = u * u + v * v;
+ if (s > 0 && s < 1) {
+ return sigma * (u * sqrt(-2.0 * log(s) / s));
+ }
+ }
+ return 0;
+}
+
+double aom_normalized_cross_correlation(const double *a, const double *b,
+ int n) {
+ double c = 0;
+ double a_len = 0;
+ double b_len = 0;
+ for (int i = 0; i < n; ++i) {
+ a_len += a[i] * a[i];
+ b_len += b[i] * b[i];
+ c += a[i] * b[i];
+ }
+ return c / (sqrt(a_len) * sqrt(b_len));
+}
+
+void aom_noise_synth(int lag, int n, const int (*coords)[2],
+ const double *coeffs, double *data, int w, int h) {
+ const int pad_size = 3 * lag;
+ const int padded_w = w + pad_size;
+ const int padded_h = h + pad_size;
+ int x = 0, y = 0;
+ double *padded = (double *)aom_malloc(padded_w * padded_h * sizeof(*padded));
+
+ for (y = 0; y < padded_h; ++y) {
+ for (x = 0; x < padded_w; ++x) {
+ padded[y * padded_w + x] = aom_randn(1.0);
+ }
+ }
+ for (y = lag; y < padded_h; ++y) {
+ for (x = lag; x < padded_w; ++x) {
+ double sum = 0;
+ int i = 0;
+ for (i = 0; i < n; ++i) {
+ const int dx = coords[i][0];
+ const int dy = coords[i][1];
+ sum += padded[(y + dy) * padded_w + (x + dx)] * coeffs[i];
+ }
+ padded[y * padded_w + x] += sum;
+ }
+ }
+ // Copy over the padded rows to the output
+ for (y = 0; y < h; ++y) {
+ memcpy(data + y * w, padded + y * padded_w, sizeof(*data) * w);
+ }
+ aom_free(padded);
+}
+
+int aom_noise_data_validate(const double *data, int w, int h) {
+ const double kVarianceThreshold = 2;
+ const double kMeanThreshold = 2;
+
+ int x = 0, y = 0;
+ int ret_value = 1;
+ double var = 0, mean = 0;
+ double *mean_x, *mean_y, *var_x, *var_y;
+
+ // Check that noise variance is not increasing in x or y
+ // and that the data is zero mean.
+ mean_x = (double *)aom_malloc(sizeof(*mean_x) * w);
+ var_x = (double *)aom_malloc(sizeof(*var_x) * w);
+ mean_y = (double *)aom_malloc(sizeof(*mean_x) * h);
+ var_y = (double *)aom_malloc(sizeof(*var_y) * h);
+
+ memset(mean_x, 0, sizeof(*mean_x) * w);
+ memset(var_x, 0, sizeof(*var_x) * w);
+ memset(mean_y, 0, sizeof(*mean_y) * h);
+ memset(var_y, 0, sizeof(*var_y) * h);
+
+ for (y = 0; y < h; ++y) {
+ for (x = 0; x < w; ++x) {
+ const double d = data[y * w + x];
+ var_x[x] += d * d;
+ var_y[y] += d * d;
+ mean_x[x] += d;
+ mean_y[y] += d;
+ var += d * d;
+ mean += d;
+ }
+ }
+ mean /= (w * h);
+ var = var / (w * h) - mean * mean;
+
+ for (y = 0; y < h; ++y) {
+ mean_y[y] /= h;
+ var_y[y] = var_y[y] / h - mean_y[y] * mean_y[y];
+ if (fabs(var_y[y] - var) >= kVarianceThreshold) {
+ fprintf(stderr, "Variance distance too large %f %f\n", var_y[y], var);
+ ret_value = 0;
+ break;
+ }
+ if (fabs(mean_y[y] - mean) >= kMeanThreshold) {
+ fprintf(stderr, "Mean distance too large %f %f\n", mean_y[y], mean);
+ ret_value = 0;
+ break;
+ }
+ }
+
+ for (x = 0; x < w; ++x) {
+ mean_x[x] /= w;
+ var_x[x] = var_x[x] / w - mean_x[x] * mean_x[x];
+ if (fabs(var_x[x] - var) >= kVarianceThreshold) {
+ fprintf(stderr, "Variance distance too large %f %f\n", var_x[x], var);
+ ret_value = 0;
+ break;
+ }
+ if (fabs(mean_x[x] - mean) >= kMeanThreshold) {
+ fprintf(stderr, "Mean distance too large %f %f\n", mean_x[x], mean);
+ ret_value = 0;
+ break;
+ }
+ }
+
+ aom_free(mean_x);
+ aom_free(mean_y);
+ aom_free(var_x);
+ aom_free(var_y);
+
+ return ret_value;
+}
