text: /*
   * 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;
  }
  