xref: /aosp_15_r20/external/libaom/test/noise_model_test.cc (revision 77c1e3ccc04c968bd2bc212e87364f250e820521)
1 /*
2  * Copyright (c) 2018, 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 <limits.h>
13 #include <math.h>
14 #include <algorithm>
15 #include <vector>
16 
17 #include "aom_dsp/noise_model.h"
18 #include "aom_dsp/noise_util.h"
19 #include "config/aom_dsp_rtcd.h"
20 #include "gtest/gtest.h"
21 #include "test/acm_random.h"
22 
23 namespace {
24 
25 // Return normally distrbuted values with standard deviation of sigma.
randn(libaom_test::ACMRandom * random,double sigma)26 double randn(libaom_test::ACMRandom *random, double sigma) {
27   while (true) {
28     const double u = 2.0 * ((double)random->Rand31() /
29                             testing::internal::Random::kMaxRange) -
30                      1.0;
31     const double v = 2.0 * ((double)random->Rand31() /
32                             testing::internal::Random::kMaxRange) -
33                      1.0;
34     const double s = u * u + v * v;
35     if (s > 0 && s < 1) {
36       return sigma * (u * sqrt(-2.0 * log(s) / s));
37     }
38   }
39 }
40 
41 // Synthesizes noise using the auto-regressive filter of the given lag,
42 // with the provided n coefficients sampled at the given coords.
noise_synth(libaom_test::ACMRandom * random,int lag,int n,const int (* coords)[2],const double * coeffs,double * data,int w,int h)43 void noise_synth(libaom_test::ACMRandom *random, int lag, int n,
44                  const int (*coords)[2], const double *coeffs, double *data,
45                  int w, int h) {
46   const int pad_size = 3 * lag;
47   const int padded_w = w + pad_size;
48   const int padded_h = h + pad_size;
49   int x = 0, y = 0;
50   std::vector<double> padded(padded_w * padded_h);
51 
52   for (y = 0; y < padded_h; ++y) {
53     for (x = 0; x < padded_w; ++x) {
54       padded[y * padded_w + x] = randn(random, 1.0);
55     }
56   }
57   for (y = lag; y < padded_h; ++y) {
58     for (x = lag; x < padded_w; ++x) {
59       double sum = 0;
60       int i = 0;
61       for (i = 0; i < n; ++i) {
62         const int dx = coords[i][0];
63         const int dy = coords[i][1];
64         sum += padded[(y + dy) * padded_w + (x + dx)] * coeffs[i];
65       }
66       padded[y * padded_w + x] += sum;
67     }
68   }
69   // Copy over the padded rows to the output
70   for (y = 0; y < h; ++y) {
71     memcpy(data + y * w, &padded[0] + y * padded_w, sizeof(*data) * w);
72   }
73 }
74 
get_noise_psd(double * noise,int width,int height,int block_size)75 std::vector<float> get_noise_psd(double *noise, int width, int height,
76                                  int block_size) {
77   float *block =
78       (float *)aom_memalign(32, block_size * block_size * sizeof(block));
79   std::vector<float> psd(block_size * block_size);
80   if (block == nullptr) {
81     EXPECT_NE(block, nullptr);
82     return psd;
83   }
84   int num_blocks = 0;
85   struct aom_noise_tx_t *tx = aom_noise_tx_malloc(block_size);
86   if (tx == nullptr) {
87     EXPECT_NE(tx, nullptr);
88     return psd;
89   }
90   for (int y = 0; y <= height - block_size; y += block_size / 2) {
91     for (int x = 0; x <= width - block_size; x += block_size / 2) {
92       for (int yy = 0; yy < block_size; ++yy) {
93         for (int xx = 0; xx < block_size; ++xx) {
94           block[yy * block_size + xx] = (float)noise[(y + yy) * width + x + xx];
95         }
96       }
97       aom_noise_tx_forward(tx, &block[0]);
98       aom_noise_tx_add_energy(tx, &psd[0]);
99       num_blocks++;
100     }
101   }
102   for (int yy = 0; yy < block_size; ++yy) {
103     for (int xx = 0; xx <= block_size / 2; ++xx) {
104       psd[yy * block_size + xx] /= num_blocks;
105     }
106   }
107   // Fill in the data that is missing due to symmetries
108   for (int xx = 1; xx < block_size / 2; ++xx) {
109     psd[(block_size - xx)] = psd[xx];
110   }
111   for (int yy = 1; yy < block_size; ++yy) {
112     for (int xx = 1; xx < block_size / 2; ++xx) {
113       psd[(block_size - yy) * block_size + (block_size - xx)] =
114           psd[yy * block_size + xx];
115     }
116   }
117   aom_noise_tx_free(tx);
118   aom_free(block);
119   return psd;
120 }
121 
122 }  // namespace
123 
TEST(NoiseStrengthSolver,GetCentersTwoBins)124 TEST(NoiseStrengthSolver, GetCentersTwoBins) {
125   aom_noise_strength_solver_t solver;
126   aom_noise_strength_solver_init(&solver, 2, 8);
127   EXPECT_NEAR(0, aom_noise_strength_solver_get_center(&solver, 0), 1e-5);
128   EXPECT_NEAR(255, aom_noise_strength_solver_get_center(&solver, 1), 1e-5);
129   aom_noise_strength_solver_free(&solver);
130 }
131 
TEST(NoiseStrengthSolver,GetCentersTwoBins10bit)132 TEST(NoiseStrengthSolver, GetCentersTwoBins10bit) {
133   aom_noise_strength_solver_t solver;
134   aom_noise_strength_solver_init(&solver, 2, 10);
135   EXPECT_NEAR(0, aom_noise_strength_solver_get_center(&solver, 0), 1e-5);
136   EXPECT_NEAR(1023, aom_noise_strength_solver_get_center(&solver, 1), 1e-5);
137   aom_noise_strength_solver_free(&solver);
138 }
139 
TEST(NoiseStrengthSolver,GetCenters256Bins)140 TEST(NoiseStrengthSolver, GetCenters256Bins) {
141   const int num_bins = 256;
142   aom_noise_strength_solver_t solver;
143   aom_noise_strength_solver_init(&solver, num_bins, 8);
144 
145   for (int i = 0; i < 256; ++i) {
146     EXPECT_NEAR(i, aom_noise_strength_solver_get_center(&solver, i), 1e-5);
147   }
148   aom_noise_strength_solver_free(&solver);
149 }
150 
151 // Tests that the noise strength solver returns the identity transform when
152 // given identity-like constraints.
TEST(NoiseStrengthSolver,ObserveIdentity)153 TEST(NoiseStrengthSolver, ObserveIdentity) {
154   const int num_bins = 256;
155   aom_noise_strength_solver_t solver;
156   ASSERT_EQ(1, aom_noise_strength_solver_init(&solver, num_bins, 8));
157 
158   // We have to add a big more strength to constraints at the boundary to
159   // overcome any regularization.
160   for (int j = 0; j < 5; ++j) {
161     aom_noise_strength_solver_add_measurement(&solver, 0, 0);
162     aom_noise_strength_solver_add_measurement(&solver, 255, 255);
163   }
164   for (int i = 0; i < 256; ++i) {
165     aom_noise_strength_solver_add_measurement(&solver, i, i);
166   }
167   EXPECT_EQ(1, aom_noise_strength_solver_solve(&solver));
168   for (int i = 2; i < num_bins - 2; ++i) {
169     EXPECT_NEAR(i, solver.eqns.x[i], 0.1);
170   }
171 
172   aom_noise_strength_lut_t lut;
173   EXPECT_EQ(1, aom_noise_strength_solver_fit_piecewise(&solver, 2, &lut));
174 
175   ASSERT_EQ(2, lut.num_points);
176   EXPECT_NEAR(0.0, lut.points[0][0], 1e-5);
177   EXPECT_NEAR(0.0, lut.points[0][1], 0.5);
178   EXPECT_NEAR(255.0, lut.points[1][0], 1e-5);
179   EXPECT_NEAR(255.0, lut.points[1][1], 0.5);
180 
181   aom_noise_strength_lut_free(&lut);
182   aom_noise_strength_solver_free(&solver);
183 }
184 
TEST(NoiseStrengthSolver,SimplifiesCurve)185 TEST(NoiseStrengthSolver, SimplifiesCurve) {
186   const int num_bins = 256;
187   aom_noise_strength_solver_t solver;
188   EXPECT_EQ(1, aom_noise_strength_solver_init(&solver, num_bins, 8));
189 
190   // Create a parabolic input
191   for (int i = 0; i < 256; ++i) {
192     const double x = (i - 127.5) / 63.5;
193     aom_noise_strength_solver_add_measurement(&solver, i, x * x);
194   }
195   EXPECT_EQ(1, aom_noise_strength_solver_solve(&solver));
196 
197   // First try to fit an unconstrained lut
198   aom_noise_strength_lut_t lut;
199   EXPECT_EQ(1, aom_noise_strength_solver_fit_piecewise(&solver, -1, &lut));
200   ASSERT_LE(20, lut.num_points);
201   aom_noise_strength_lut_free(&lut);
202 
203   // Now constrain the maximum number of points
204   const int kMaxPoints = 9;
205   EXPECT_EQ(1,
206             aom_noise_strength_solver_fit_piecewise(&solver, kMaxPoints, &lut));
207   ASSERT_EQ(kMaxPoints, lut.num_points);
208 
209   // Check that the input parabola is still well represented
210   EXPECT_NEAR(0.0, lut.points[0][0], 1e-5);
211   EXPECT_NEAR(4.0, lut.points[0][1], 0.1);
212   for (int i = 1; i < lut.num_points - 1; ++i) {
213     const double x = (lut.points[i][0] - 128.) / 64.;
214     EXPECT_NEAR(x * x, lut.points[i][1], 0.1);
215   }
216   EXPECT_NEAR(255.0, lut.points[kMaxPoints - 1][0], 1e-5);
217 
218   EXPECT_NEAR(4.0, lut.points[kMaxPoints - 1][1], 0.1);
219   aom_noise_strength_lut_free(&lut);
220   aom_noise_strength_solver_free(&solver);
221 }
222 
TEST(NoiseStrengthLut,LutInitNegativeOrZeroSize)223 TEST(NoiseStrengthLut, LutInitNegativeOrZeroSize) {
224   aom_noise_strength_lut_t lut;
225   ASSERT_FALSE(aom_noise_strength_lut_init(&lut, -1));
226   ASSERT_FALSE(aom_noise_strength_lut_init(&lut, 0));
227 }
228 
TEST(NoiseStrengthLut,LutEvalSinglePoint)229 TEST(NoiseStrengthLut, LutEvalSinglePoint) {
230   aom_noise_strength_lut_t lut;
231   ASSERT_TRUE(aom_noise_strength_lut_init(&lut, 1));
232   ASSERT_EQ(1, lut.num_points);
233   lut.points[0][0] = 0;
234   lut.points[0][1] = 1;
235   EXPECT_EQ(1, aom_noise_strength_lut_eval(&lut, -1));
236   EXPECT_EQ(1, aom_noise_strength_lut_eval(&lut, 0));
237   EXPECT_EQ(1, aom_noise_strength_lut_eval(&lut, 1));
238   aom_noise_strength_lut_free(&lut);
239 }
240 
TEST(NoiseStrengthLut,LutEvalMultiPointInterp)241 TEST(NoiseStrengthLut, LutEvalMultiPointInterp) {
242   const double kEps = 1e-5;
243   aom_noise_strength_lut_t lut;
244   ASSERT_TRUE(aom_noise_strength_lut_init(&lut, 4));
245   ASSERT_EQ(4, lut.num_points);
246 
247   lut.points[0][0] = 0;
248   lut.points[0][1] = 0;
249 
250   lut.points[1][0] = 1;
251   lut.points[1][1] = 1;
252 
253   lut.points[2][0] = 2;
254   lut.points[2][1] = 1;
255 
256   lut.points[3][0] = 100;
257   lut.points[3][1] = 1001;
258 
259   // Test lower boundary
260   EXPECT_EQ(0, aom_noise_strength_lut_eval(&lut, -1));
261   EXPECT_EQ(0, aom_noise_strength_lut_eval(&lut, 0));
262 
263   // Test first part that should be identity
264   EXPECT_NEAR(0.25, aom_noise_strength_lut_eval(&lut, 0.25), kEps);
265   EXPECT_NEAR(0.75, aom_noise_strength_lut_eval(&lut, 0.75), kEps);
266 
267   // This is a constant section (should evaluate to 1)
268   EXPECT_NEAR(1.0, aom_noise_strength_lut_eval(&lut, 1.25), kEps);
269   EXPECT_NEAR(1.0, aom_noise_strength_lut_eval(&lut, 1.75), kEps);
270 
271   // Test interpolation between to non-zero y coords.
272   EXPECT_NEAR(1, aom_noise_strength_lut_eval(&lut, 2), kEps);
273   EXPECT_NEAR(251, aom_noise_strength_lut_eval(&lut, 26.5), kEps);
274   EXPECT_NEAR(751, aom_noise_strength_lut_eval(&lut, 75.5), kEps);
275 
276   // Test upper boundary
277   EXPECT_EQ(1001, aom_noise_strength_lut_eval(&lut, 100));
278   EXPECT_EQ(1001, aom_noise_strength_lut_eval(&lut, 101));
279 
280   aom_noise_strength_lut_free(&lut);
281 }
282 
TEST(NoiseModel,InitSuccessWithValidSquareShape)283 TEST(NoiseModel, InitSuccessWithValidSquareShape) {
284   aom_noise_model_params_t params = { AOM_NOISE_SHAPE_SQUARE, 2, 8, 0 };
285   aom_noise_model_t model;
286 
287   EXPECT_TRUE(aom_noise_model_init(&model, params));
288 
289   const int kNumCoords = 12;
290   const int kCoords[][2] = { { -2, -2 }, { -1, -2 }, { 0, -2 },  { 1, -2 },
291                              { 2, -2 },  { -2, -1 }, { -1, -1 }, { 0, -1 },
292                              { 1, -1 },  { 2, -1 },  { -2, 0 },  { -1, 0 } };
293   EXPECT_EQ(kNumCoords, model.n);
294   for (int i = 0; i < kNumCoords; ++i) {
295     const int *coord = kCoords[i];
296     EXPECT_EQ(coord[0], model.coords[i][0]);
297     EXPECT_EQ(coord[1], model.coords[i][1]);
298   }
299   aom_noise_model_free(&model);
300 }
301 
TEST(NoiseModel,InitSuccessWithValidDiamondShape)302 TEST(NoiseModel, InitSuccessWithValidDiamondShape) {
303   aom_noise_model_t model;
304   aom_noise_model_params_t params = { AOM_NOISE_SHAPE_DIAMOND, 2, 8, 0 };
305   EXPECT_TRUE(aom_noise_model_init(&model, params));
306   EXPECT_EQ(6, model.n);
307   const int kNumCoords = 6;
308   const int kCoords[][2] = { { 0, -2 }, { -1, -1 }, { 0, -1 },
309                              { 1, -1 }, { -2, 0 },  { -1, 0 } };
310   EXPECT_EQ(kNumCoords, model.n);
311   for (int i = 0; i < kNumCoords; ++i) {
312     const int *coord = kCoords[i];
313     EXPECT_EQ(coord[0], model.coords[i][0]);
314     EXPECT_EQ(coord[1], model.coords[i][1]);
315   }
316   aom_noise_model_free(&model);
317 }
318 
TEST(NoiseModel,InitFailsWithTooLargeLag)319 TEST(NoiseModel, InitFailsWithTooLargeLag) {
320   aom_noise_model_t model;
321   aom_noise_model_params_t params = { AOM_NOISE_SHAPE_SQUARE, 10, 8, 0 };
322   EXPECT_FALSE(aom_noise_model_init(&model, params));
323   aom_noise_model_free(&model);
324 }
325 
TEST(NoiseModel,InitFailsWithTooSmallLag)326 TEST(NoiseModel, InitFailsWithTooSmallLag) {
327   aom_noise_model_t model;
328   aom_noise_model_params_t params = { AOM_NOISE_SHAPE_SQUARE, 0, 8, 0 };
329   EXPECT_FALSE(aom_noise_model_init(&model, params));
330   aom_noise_model_free(&model);
331 }
332 
TEST(NoiseModel,InitFailsWithInvalidShape)333 TEST(NoiseModel, InitFailsWithInvalidShape) {
334   aom_noise_model_t model;
335   aom_noise_model_params_t params = { aom_noise_shape(100), 3, 8, 0 };
336   EXPECT_FALSE(aom_noise_model_init(&model, params));
337   aom_noise_model_free(&model);
338 }
339 
TEST(NoiseModel,InitFailsWithInvalidBitdepth)340 TEST(NoiseModel, InitFailsWithInvalidBitdepth) {
341   aom_noise_model_t model;
342   aom_noise_model_params_t params = { AOM_NOISE_SHAPE_SQUARE, 2, 8, 0 };
343   for (int i = 0; i <= 32; ++i) {
344     params.bit_depth = i;
345     if (i == 8 || i == 10 || i == 12) {
346       EXPECT_TRUE(aom_noise_model_init(&model, params)) << "bit_depth: " << i;
347       aom_noise_model_free(&model);
348     } else {
349       EXPECT_FALSE(aom_noise_model_init(&model, params)) << "bit_depth: " << i;
350     }
351   }
352   params.bit_depth = INT_MAX;
353   EXPECT_FALSE(aom_noise_model_init(&model, params));
354 }
355 
356 // A container template class to hold a data type and extra arguments.
357 // All of these args are bundled into one struct so that we can use
358 // parameterized tests on combinations of supported data types
359 // (uint8_t and uint16_t) and bit depths (8, 10, 12).
360 template <typename T, int bit_depth, bool use_highbd>
361 struct BitDepthParams {
362   typedef T data_type_t;
363   static const int kBitDepth = bit_depth;
364   static const bool kUseHighBD = use_highbd;
365 };
366 
367 template <typename T>
368 class FlatBlockEstimatorTest : public ::testing::Test, public T {
369  public:
SetUp()370   void SetUp() override { random_.Reset(171); }
371   typedef std::vector<typename T::data_type_t> VecType;
372   VecType data_;
373   libaom_test::ACMRandom random_;
374 };
375 
376 TYPED_TEST_SUITE_P(FlatBlockEstimatorTest);
377 
TYPED_TEST_P(FlatBlockEstimatorTest,ExtractBlock)378 TYPED_TEST_P(FlatBlockEstimatorTest, ExtractBlock) {
379   const int kBlockSize = 16;
380   aom_flat_block_finder_t flat_block_finder;
381   ASSERT_EQ(1, aom_flat_block_finder_init(&flat_block_finder, kBlockSize,
382                                           this->kBitDepth, this->kUseHighBD));
383   const double normalization = flat_block_finder.normalization;
384 
385   // Test with an image of more than one block.
386   const int h = 2 * kBlockSize;
387   const int w = 2 * kBlockSize;
388   const int stride = 2 * kBlockSize;
389   this->data_.resize(h * stride, 128);
390 
391   // Set up the (0,0) block to be a plane and the (0,1) block to be a
392   // checkerboard
393   const int shift = this->kBitDepth - 8;
394   for (int y = 0; y < kBlockSize; ++y) {
395     for (int x = 0; x < kBlockSize; ++x) {
396       this->data_[y * stride + x] = (-y + x + 128) << shift;
397       this->data_[y * stride + x + kBlockSize] =
398           ((x % 2 + y % 2) % 2 ? 128 - 20 : 128 + 20) << shift;
399     }
400   }
401   std::vector<double> block(kBlockSize * kBlockSize, 1);
402   std::vector<double> plane(kBlockSize * kBlockSize, 1);
403 
404   // The block data should be a constant (zero) and the rest of the plane
405   // trend is covered in the plane data.
406   aom_flat_block_finder_extract_block(&flat_block_finder,
407                                       (uint8_t *)&this->data_[0], w, h, stride,
408                                       0, 0, &plane[0], &block[0]);
409   for (int y = 0; y < kBlockSize; ++y) {
410     for (int x = 0; x < kBlockSize; ++x) {
411       EXPECT_NEAR(0, block[y * kBlockSize + x], 1e-5);
412       EXPECT_NEAR((double)(this->data_[y * stride + x]) / normalization,
413                   plane[y * kBlockSize + x], 1e-5);
414     }
415   }
416 
417   // The plane trend is a constant, and the block is a zero mean checkerboard.
418   aom_flat_block_finder_extract_block(&flat_block_finder,
419                                       (uint8_t *)&this->data_[0], w, h, stride,
420                                       kBlockSize, 0, &plane[0], &block[0]);
421   const int mid = 128 << shift;
422   for (int y = 0; y < kBlockSize; ++y) {
423     for (int x = 0; x < kBlockSize; ++x) {
424       EXPECT_NEAR(((double)this->data_[y * stride + x + kBlockSize] - mid) /
425                       normalization,
426                   block[y * kBlockSize + x], 1e-5);
427       EXPECT_NEAR(mid / normalization, plane[y * kBlockSize + x], 1e-5);
428     }
429   }
430   aom_flat_block_finder_free(&flat_block_finder);
431 }
432 
TYPED_TEST_P(FlatBlockEstimatorTest,FindFlatBlocks)433 TYPED_TEST_P(FlatBlockEstimatorTest, FindFlatBlocks) {
434   const int kBlockSize = 32;
435   aom_flat_block_finder_t flat_block_finder;
436   ASSERT_EQ(1, aom_flat_block_finder_init(&flat_block_finder, kBlockSize,
437                                           this->kBitDepth, this->kUseHighBD));
438 
439   const int num_blocks_w = 8;
440   const int h = kBlockSize;
441   const int w = kBlockSize * num_blocks_w;
442   const int stride = w;
443   this->data_.resize(h * stride, 128);
444   std::vector<uint8_t> flat_blocks(num_blocks_w, 0);
445 
446   const int shift = this->kBitDepth - 8;
447   for (int y = 0; y < kBlockSize; ++y) {
448     for (int x = 0; x < kBlockSize; ++x) {
449       // Block 0 (not flat): constant doesn't have enough variance to qualify
450       this->data_[y * stride + x + 0 * kBlockSize] = 128 << shift;
451 
452       // Block 1 (not flat): too high of variance is hard to validate as flat
453       this->data_[y * stride + x + 1 * kBlockSize] =
454           ((uint8_t)(128 + randn(&this->random_, 5))) << shift;
455 
456       // Block 2 (flat): slight checkerboard added to constant
457       const int check = (x % 2 + y % 2) % 2 ? -2 : 2;
458       this->data_[y * stride + x + 2 * kBlockSize] = (128 + check) << shift;
459 
460       // Block 3 (flat): planar block with checkerboard pattern is also flat
461       this->data_[y * stride + x + 3 * kBlockSize] =
462           (y * 2 - x / 2 + 128 + check) << shift;
463 
464       // Block 4 (flat): gaussian random with standard deviation 1.
465       this->data_[y * stride + x + 4 * kBlockSize] =
466           ((uint8_t)(randn(&this->random_, 1) + x + 128.0)) << shift;
467 
468       // Block 5 (flat): gaussian random with standard deviation 2.
469       this->data_[y * stride + x + 5 * kBlockSize] =
470           ((uint8_t)(randn(&this->random_, 2) + y + 128.0)) << shift;
471 
472       // Block 6 (not flat): too high of directional gradient.
473       const int strong_edge = x > kBlockSize / 2 ? 64 : 0;
474       this->data_[y * stride + x + 6 * kBlockSize] =
475           ((uint8_t)(randn(&this->random_, 1) + strong_edge + 128.0)) << shift;
476 
477       // Block 7 (not flat): too high gradient.
478       const int big_check = ((x >> 2) % 2 + (y >> 2) % 2) % 2 ? -16 : 16;
479       this->data_[y * stride + x + 7 * kBlockSize] =
480           ((uint8_t)(randn(&this->random_, 1) + big_check + 128.0)) << shift;
481     }
482   }
483 
484   EXPECT_EQ(4, aom_flat_block_finder_run(&flat_block_finder,
485                                          (uint8_t *)&this->data_[0], w, h,
486                                          stride, &flat_blocks[0]));
487 
488   // First two blocks are not flat
489   EXPECT_EQ(0, flat_blocks[0]);
490   EXPECT_EQ(0, flat_blocks[1]);
491 
492   // Next 4 blocks are flat.
493   EXPECT_EQ(255, flat_blocks[2]);
494   EXPECT_EQ(255, flat_blocks[3]);
495   EXPECT_EQ(255, flat_blocks[4]);
496   EXPECT_EQ(255, flat_blocks[5]);
497 
498   // Last 2 are not flat by threshold
499   EXPECT_EQ(0, flat_blocks[6]);
500   EXPECT_EQ(0, flat_blocks[7]);
501 
502   // Add the noise from non-flat block 1 to every block.
503   for (int y = 0; y < kBlockSize; ++y) {
504     for (int x = 0; x < kBlockSize * num_blocks_w; ++x) {
505       this->data_[y * stride + x] +=
506           (this->data_[y * stride + x % kBlockSize + kBlockSize] -
507            (128 << shift));
508     }
509   }
510   // Now the scored selection will pick the one that is most likely flat (block
511   // 0)
512   EXPECT_EQ(1, aom_flat_block_finder_run(&flat_block_finder,
513                                          (uint8_t *)&this->data_[0], w, h,
514                                          stride, &flat_blocks[0]));
515   EXPECT_EQ(1, flat_blocks[0]);
516   EXPECT_EQ(0, flat_blocks[1]);
517   EXPECT_EQ(0, flat_blocks[2]);
518   EXPECT_EQ(0, flat_blocks[3]);
519   EXPECT_EQ(0, flat_blocks[4]);
520   EXPECT_EQ(0, flat_blocks[5]);
521   EXPECT_EQ(0, flat_blocks[6]);
522   EXPECT_EQ(0, flat_blocks[7]);
523 
524   aom_flat_block_finder_free(&flat_block_finder);
525 }
526 
527 REGISTER_TYPED_TEST_SUITE_P(FlatBlockEstimatorTest, ExtractBlock,
528                             FindFlatBlocks);
529 
530 typedef ::testing::Types<BitDepthParams<uint8_t, 8, false>,   // lowbd
531                          BitDepthParams<uint16_t, 8, true>,   // lowbd in 16-bit
532                          BitDepthParams<uint16_t, 10, true>,  // highbd data
533                          BitDepthParams<uint16_t, 12, true> >
534     AllBitDepthParams;
535 // Note the empty final argument can be removed if C++20 is made the minimum
536 // requirement.
537 INSTANTIATE_TYPED_TEST_SUITE_P(FlatBlockInstatiation, FlatBlockEstimatorTest,
538                                AllBitDepthParams, );
539 
540 template <typename T>
541 class NoiseModelUpdateTest : public ::testing::Test, public T {
542  public:
543   static const int kWidth = 128;
544   static const int kHeight = 128;
545   static const int kBlockSize = 16;
546   static const int kNumBlocksX = kWidth / kBlockSize;
547   static const int kNumBlocksY = kHeight / kBlockSize;
548 
SetUp()549   void SetUp() override {
550     const aom_noise_model_params_t params = { AOM_NOISE_SHAPE_SQUARE, 3,
551                                               T::kBitDepth, T::kUseHighBD };
552     ASSERT_TRUE(aom_noise_model_init(&model_, params));
553 
554     random_.Reset(100171);
555 
556     data_.resize(kWidth * kHeight * 3);
557     denoised_.resize(kWidth * kHeight * 3);
558     noise_.resize(kWidth * kHeight * 3);
559     renoise_.resize(kWidth * kHeight);
560     flat_blocks_.resize(kNumBlocksX * kNumBlocksY);
561 
562     for (int c = 0, offset = 0; c < 3; ++c, offset += kWidth * kHeight) {
563       data_ptr_[c] = &data_[offset];
564       noise_ptr_[c] = &noise_[offset];
565       denoised_ptr_[c] = &denoised_[offset];
566       strides_[c] = kWidth;
567 
568       data_ptr_raw_[c] = (uint8_t *)&data_[offset];
569       denoised_ptr_raw_[c] = (uint8_t *)&denoised_[offset];
570     }
571     chroma_sub_[0] = 0;
572     chroma_sub_[1] = 0;
573   }
574 
NoiseModelUpdate(int block_size=kBlockSize)575   int NoiseModelUpdate(int block_size = kBlockSize) {
576     return aom_noise_model_update(&model_, data_ptr_raw_, denoised_ptr_raw_,
577                                   kWidth, kHeight, strides_, chroma_sub_,
578                                   &flat_blocks_[0], block_size);
579   }
580 
TearDown()581   void TearDown() override { aom_noise_model_free(&model_); }
582 
583  protected:
584   aom_noise_model_t model_;
585   std::vector<typename T::data_type_t> data_;
586   std::vector<typename T::data_type_t> denoised_;
587 
588   std::vector<double> noise_;
589   std::vector<double> renoise_;
590   std::vector<uint8_t> flat_blocks_;
591 
592   typename T::data_type_t *data_ptr_[3];
593   typename T::data_type_t *denoised_ptr_[3];
594 
595   double *noise_ptr_[3];
596   int strides_[3];
597   int chroma_sub_[2];
598   libaom_test::ACMRandom random_;
599 
600  private:
601   uint8_t *data_ptr_raw_[3];
602   uint8_t *denoised_ptr_raw_[3];
603 };
604 
605 TYPED_TEST_SUITE_P(NoiseModelUpdateTest);
606 
TYPED_TEST_P(NoiseModelUpdateTest,UpdateFailsNoFlatBlocks)607 TYPED_TEST_P(NoiseModelUpdateTest, UpdateFailsNoFlatBlocks) {
608   EXPECT_EQ(AOM_NOISE_STATUS_INSUFFICIENT_FLAT_BLOCKS,
609             this->NoiseModelUpdate());
610 }
611 
TYPED_TEST_P(NoiseModelUpdateTest,UpdateSuccessForZeroNoiseAllFlat)612 TYPED_TEST_P(NoiseModelUpdateTest, UpdateSuccessForZeroNoiseAllFlat) {
613   this->flat_blocks_.assign(this->flat_blocks_.size(), 1);
614   this->denoised_.assign(this->denoised_.size(), 128);
615   this->data_.assign(this->denoised_.size(), 128);
616   EXPECT_EQ(AOM_NOISE_STATUS_INTERNAL_ERROR, this->NoiseModelUpdate());
617 }
618 
TYPED_TEST_P(NoiseModelUpdateTest,UpdateFailsBlockSizeTooSmall)619 TYPED_TEST_P(NoiseModelUpdateTest, UpdateFailsBlockSizeTooSmall) {
620   this->flat_blocks_.assign(this->flat_blocks_.size(), 1);
621   this->denoised_.assign(this->denoised_.size(), 128);
622   this->data_.assign(this->denoised_.size(), 128);
623   EXPECT_EQ(AOM_NOISE_STATUS_INVALID_ARGUMENT,
624             this->NoiseModelUpdate(6 /* block_size=6 is too small*/));
625 }
626 
TYPED_TEST_P(NoiseModelUpdateTest,UpdateSuccessForWhiteRandomNoise)627 TYPED_TEST_P(NoiseModelUpdateTest, UpdateSuccessForWhiteRandomNoise) {
628   aom_noise_model_t &model = this->model_;
629   const int width = this->kWidth;
630   const int height = this->kHeight;
631 
632   const int shift = this->kBitDepth - 8;
633   for (int y = 0; y < height; ++y) {
634     for (int x = 0; x < width; ++x) {
635       this->data_ptr_[0][y * width + x] = int(64 + y + randn(&this->random_, 1))
636                                           << shift;
637       this->denoised_ptr_[0][y * width + x] = (64 + y) << shift;
638       // Make the chroma planes completely correlated with the Y plane
639       for (int c = 1; c < 3; ++c) {
640         this->data_ptr_[c][y * width + x] = this->data_ptr_[0][y * width + x];
641         this->denoised_ptr_[c][y * width + x] =
642             this->denoised_ptr_[0][y * width + x];
643       }
644     }
645   }
646   this->flat_blocks_.assign(this->flat_blocks_.size(), 1);
647   EXPECT_EQ(AOM_NOISE_STATUS_OK, this->NoiseModelUpdate());
648 
649   const double kCoeffEps = 0.075;
650   const int n = model.n;
651   for (int c = 0; c < 3; ++c) {
652     for (int i = 0; i < n; ++i) {
653       EXPECT_NEAR(0, model.latest_state[c].eqns.x[i], kCoeffEps);
654       EXPECT_NEAR(0, model.combined_state[c].eqns.x[i], kCoeffEps);
655     }
656     // The second and third channels are highly correlated with the first.
657     if (c > 0) {
658       ASSERT_EQ(n + 1, model.latest_state[c].eqns.n);
659       ASSERT_EQ(n + 1, model.combined_state[c].eqns.n);
660 
661       EXPECT_NEAR(1, model.latest_state[c].eqns.x[n], kCoeffEps);
662       EXPECT_NEAR(1, model.combined_state[c].eqns.x[n], kCoeffEps);
663     }
664   }
665 
666   // The fitted noise strength should be close to the standard deviation
667   // for all intensity bins.
668   const double kStdEps = 0.1;
669   const double normalize = 1 << shift;
670 
671   for (int i = 0; i < model.latest_state[0].strength_solver.eqns.n; ++i) {
672     EXPECT_NEAR(1.0,
673                 model.latest_state[0].strength_solver.eqns.x[i] / normalize,
674                 kStdEps);
675     EXPECT_NEAR(1.0,
676                 model.combined_state[0].strength_solver.eqns.x[i] / normalize,
677                 kStdEps);
678   }
679 
680   aom_noise_strength_lut_t lut;
681   aom_noise_strength_solver_fit_piecewise(
682       &model.latest_state[0].strength_solver, -1, &lut);
683   ASSERT_EQ(2, lut.num_points);
684   EXPECT_NEAR(0.0, lut.points[0][0], 1e-5);
685   EXPECT_NEAR(1.0, lut.points[0][1] / normalize, kStdEps);
686   EXPECT_NEAR((1 << this->kBitDepth) - 1, lut.points[1][0], 1e-5);
687   EXPECT_NEAR(1.0, lut.points[1][1] / normalize, kStdEps);
688   aom_noise_strength_lut_free(&lut);
689 }
690 
TYPED_TEST_P(NoiseModelUpdateTest,UpdateSuccessForScaledWhiteNoise)691 TYPED_TEST_P(NoiseModelUpdateTest, UpdateSuccessForScaledWhiteNoise) {
692   aom_noise_model_t &model = this->model_;
693   const int width = this->kWidth;
694   const int height = this->kHeight;
695 
696   const double kCoeffEps = 0.055;
697   const double kLowStd = 1;
698   const double kHighStd = 4;
699   const int shift = this->kBitDepth - 8;
700   for (int y = 0; y < height; ++y) {
701     for (int x = 0; x < width; ++x) {
702       for (int c = 0; c < 3; ++c) {
703         // The image data is bimodal:
704         // Bottom half has low intensity and low noise strength
705         // Top half has high intensity and high noise strength
706         const int avg = (y < height / 2) ? 4 : 245;
707         const double std = (y < height / 2) ? kLowStd : kHighStd;
708         this->data_ptr_[c][y * width + x] =
709             ((uint8_t)std::min((int)255,
710                                (int)(2 + avg + randn(&this->random_, std))))
711             << shift;
712         this->denoised_ptr_[c][y * width + x] = (2 + avg) << shift;
713       }
714     }
715   }
716   // Label all blocks as flat for the update
717   this->flat_blocks_.assign(this->flat_blocks_.size(), 1);
718   EXPECT_EQ(AOM_NOISE_STATUS_OK, this->NoiseModelUpdate());
719 
720   const int n = model.n;
721   // The noise is uncorrelated spatially and with the y channel.
722   // All coefficients should be reasonably close to zero.
723   for (int c = 0; c < 3; ++c) {
724     for (int i = 0; i < n; ++i) {
725       EXPECT_NEAR(0, model.latest_state[c].eqns.x[i], kCoeffEps);
726       EXPECT_NEAR(0, model.combined_state[c].eqns.x[i], kCoeffEps);
727     }
728     if (c > 0) {
729       ASSERT_EQ(n + 1, model.latest_state[c].eqns.n);
730       ASSERT_EQ(n + 1, model.combined_state[c].eqns.n);
731 
732       // The correlation to the y channel should be low (near zero)
733       EXPECT_NEAR(0, model.latest_state[c].eqns.x[n], kCoeffEps);
734       EXPECT_NEAR(0, model.combined_state[c].eqns.x[n], kCoeffEps);
735     }
736   }
737 
738   // Noise strength should vary between kLowStd and kHighStd.
739   const double kStdEps = 0.15;
740   // We have to normalize fitted standard deviation based on bit depth.
741   const double normalize = (1 << shift);
742 
743   ASSERT_EQ(20, model.latest_state[0].strength_solver.eqns.n);
744   for (int i = 0; i < model.latest_state[0].strength_solver.eqns.n; ++i) {
745     const double a = i / 19.0;
746     const double expected = (kLowStd * (1.0 - a) + kHighStd * a);
747     EXPECT_NEAR(expected,
748                 model.latest_state[0].strength_solver.eqns.x[i] / normalize,
749                 kStdEps);
750     EXPECT_NEAR(expected,
751                 model.combined_state[0].strength_solver.eqns.x[i] / normalize,
752                 kStdEps);
753   }
754 
755   // If we fit a piecewise linear model, there should be two points:
756   // one near kLowStd at 0, and the other near kHighStd and 255.
757   aom_noise_strength_lut_t lut;
758   aom_noise_strength_solver_fit_piecewise(
759       &model.latest_state[0].strength_solver, 2, &lut);
760   ASSERT_EQ(2, lut.num_points);
761   EXPECT_NEAR(0, lut.points[0][0], 1e-4);
762   EXPECT_NEAR(kLowStd, lut.points[0][1] / normalize, kStdEps);
763   EXPECT_NEAR((1 << this->kBitDepth) - 1, lut.points[1][0], 1e-5);
764   EXPECT_NEAR(kHighStd, lut.points[1][1] / normalize, kStdEps);
765   aom_noise_strength_lut_free(&lut);
766 }
767 
TYPED_TEST_P(NoiseModelUpdateTest,UpdateSuccessForCorrelatedNoise)768 TYPED_TEST_P(NoiseModelUpdateTest, UpdateSuccessForCorrelatedNoise) {
769   aom_noise_model_t &model = this->model_;
770   const int width = this->kWidth;
771   const int height = this->kHeight;
772   const int kNumCoeffs = 24;
773   const double kStd = 4;
774   const double kStdEps = 0.3;
775   const double kCoeffEps = 0.065;
776   // Use different coefficients for each channel
777   const double kCoeffs[3][24] = {
778     { 0.02884, -0.03356, 0.00633,  0.01757,  0.02849,  -0.04620,
779       0.02833, -0.07178, 0.07076,  -0.11603, -0.10413, -0.16571,
780       0.05158, -0.07969, 0.02640,  -0.07191, 0.02530,  0.41968,
781       0.21450, -0.00702, -0.01401, -0.03676, -0.08713, 0.44196 },
782     { 0.00269, -0.01291, -0.01513, 0.07234,  0.03208,   0.00477,
783       0.00226, -0.00254, 0.03533,  0.12841,  -0.25970,  -0.06336,
784       0.05238, -0.00845, -0.03118, 0.09043,  -0.36558,  0.48903,
785       0.00595, -0.11938, 0.02106,  0.095956, -0.350139, 0.59305 },
786     { -0.00643, -0.01080, -0.01466, 0.06951, 0.03707,  -0.00482,
787       0.00817,  -0.00909, 0.02949,  0.12181, -0.25210, -0.07886,
788       0.06083,  -0.01210, -0.03108, 0.08944, -0.35875, 0.49150,
789       0.00415,  -0.12905, 0.02870,  0.09740, -0.34610, 0.58824 },
790   };
791 
792   ASSERT_EQ(model.n, kNumCoeffs);
793   this->chroma_sub_[0] = this->chroma_sub_[1] = 1;
794 
795   this->flat_blocks_.assign(this->flat_blocks_.size(), 1);
796 
797   // Add different noise onto each plane
798   const int shift = this->kBitDepth - 8;
799   for (int c = 0; c < 3; ++c) {
800     noise_synth(&this->random_, model.params.lag, model.n, model.coords,
801                 kCoeffs[c], this->noise_ptr_[c], width, height);
802     const int x_shift = c > 0 ? this->chroma_sub_[0] : 0;
803     const int y_shift = c > 0 ? this->chroma_sub_[1] : 0;
804     for (int y = 0; y < (height >> y_shift); ++y) {
805       for (int x = 0; x < (width >> x_shift); ++x) {
806         const uint8_t value = 64 + x / 2 + y / 4;
807         this->data_ptr_[c][y * width + x] =
808             (uint8_t(value + this->noise_ptr_[c][y * width + x] * kStd))
809             << shift;
810         this->denoised_ptr_[c][y * width + x] = value << shift;
811       }
812     }
813   }
814   EXPECT_EQ(AOM_NOISE_STATUS_OK, this->NoiseModelUpdate());
815 
816   // For the Y plane, the solved coefficients should be close to the original
817   const int n = model.n;
818   for (int c = 0; c < 3; ++c) {
819     for (int i = 0; i < n; ++i) {
820       EXPECT_NEAR(kCoeffs[c][i], model.latest_state[c].eqns.x[i], kCoeffEps);
821       EXPECT_NEAR(kCoeffs[c][i], model.combined_state[c].eqns.x[i], kCoeffEps);
822     }
823     // The chroma planes should be uncorrelated with the luma plane
824     if (c > 0) {
825       EXPECT_NEAR(0, model.latest_state[c].eqns.x[n], kCoeffEps);
826       EXPECT_NEAR(0, model.combined_state[c].eqns.x[n], kCoeffEps);
827     }
828     // Correlation between the coefficient vector and the fitted coefficients
829     // should be close to 1.
830     EXPECT_LT(0.98, aom_normalized_cross_correlation(
831                         model.latest_state[c].eqns.x, kCoeffs[c], kNumCoeffs));
832 
833     noise_synth(&this->random_, model.params.lag, model.n, model.coords,
834                 model.latest_state[c].eqns.x, &this->renoise_[0], width,
835                 height);
836 
837     EXPECT_TRUE(aom_noise_data_validate(&this->renoise_[0], width, height));
838   }
839 
840   // Check fitted noise strength
841   const double normalize = 1 << shift;
842   for (int c = 0; c < 3; ++c) {
843     for (int i = 0; i < model.latest_state[c].strength_solver.eqns.n; ++i) {
844       EXPECT_NEAR(kStd,
845                   model.latest_state[c].strength_solver.eqns.x[i] / normalize,
846                   kStdEps);
847     }
848   }
849 }
850 
TYPED_TEST_P(NoiseModelUpdateTest,NoiseStrengthChangeSignalsDifferentNoiseType)851 TYPED_TEST_P(NoiseModelUpdateTest,
852              NoiseStrengthChangeSignalsDifferentNoiseType) {
853   aom_noise_model_t &model = this->model_;
854   const int width = this->kWidth;
855   const int height = this->kHeight;
856   const int block_size = this->kBlockSize;
857   // Create a gradient image with std = 2 uncorrelated noise
858   const double kStd = 2;
859   const int shift = this->kBitDepth - 8;
860 
861   for (int i = 0; i < width * height; ++i) {
862     const uint8_t val = (i % width) < width / 2 ? 64 : 192;
863     for (int c = 0; c < 3; ++c) {
864       this->noise_ptr_[c][i] = randn(&this->random_, 1);
865       this->data_ptr_[c][i] = ((uint8_t)(this->noise_ptr_[c][i] * kStd + val))
866                               << shift;
867       this->denoised_ptr_[c][i] = val << shift;
868     }
869   }
870   this->flat_blocks_.assign(this->flat_blocks_.size(), 1);
871   EXPECT_EQ(AOM_NOISE_STATUS_OK, this->NoiseModelUpdate());
872 
873   const int kNumBlocks = width * height / block_size / block_size;
874   EXPECT_EQ(kNumBlocks, model.latest_state[0].strength_solver.num_equations);
875   EXPECT_EQ(kNumBlocks, model.latest_state[1].strength_solver.num_equations);
876   EXPECT_EQ(kNumBlocks, model.latest_state[2].strength_solver.num_equations);
877   EXPECT_EQ(kNumBlocks, model.combined_state[0].strength_solver.num_equations);
878   EXPECT_EQ(kNumBlocks, model.combined_state[1].strength_solver.num_equations);
879   EXPECT_EQ(kNumBlocks, model.combined_state[2].strength_solver.num_equations);
880 
881   // Bump up noise by an insignificant amount
882   for (int i = 0; i < width * height; ++i) {
883     const uint8_t val = (i % width) < width / 2 ? 64 : 192;
884     this->data_ptr_[0][i] =
885         ((uint8_t)(this->noise_ptr_[0][i] * (kStd + 0.085) + val)) << shift;
886   }
887   EXPECT_EQ(AOM_NOISE_STATUS_OK, this->NoiseModelUpdate());
888 
889   const double kARGainTolerance = 0.02;
890   for (int c = 0; c < 3; ++c) {
891     EXPECT_EQ(kNumBlocks, model.latest_state[c].strength_solver.num_equations);
892     EXPECT_EQ(15250, model.latest_state[c].num_observations);
893     EXPECT_NEAR(1, model.latest_state[c].ar_gain, kARGainTolerance);
894 
895     EXPECT_EQ(2 * kNumBlocks,
896               model.combined_state[c].strength_solver.num_equations);
897     EXPECT_EQ(2 * 15250, model.combined_state[c].num_observations);
898     EXPECT_NEAR(1, model.combined_state[c].ar_gain, kARGainTolerance);
899   }
900 
901   // Bump up the noise strength on half the image for one channel by a
902   // significant amount.
903   for (int i = 0; i < width * height; ++i) {
904     const uint8_t val = (i % width) < width / 2 ? 64 : 128;
905     if (i % width < width / 2) {
906       this->data_ptr_[0][i] =
907           ((uint8_t)(randn(&this->random_, kStd + 0.5) + val)) << shift;
908     }
909   }
910   EXPECT_EQ(AOM_NOISE_STATUS_DIFFERENT_NOISE_TYPE, this->NoiseModelUpdate());
911 
912   // Since we didn't update the combined state, it should still be at 2 *
913   // num_blocks
914   EXPECT_EQ(kNumBlocks, model.latest_state[0].strength_solver.num_equations);
915   EXPECT_EQ(2 * kNumBlocks,
916             model.combined_state[0].strength_solver.num_equations);
917 
918   // In normal operation, the "latest" estimate can be saved to the "combined"
919   // state for continued updates.
920   aom_noise_model_save_latest(&model);
921   for (int c = 0; c < 3; ++c) {
922     EXPECT_EQ(kNumBlocks, model.latest_state[c].strength_solver.num_equations);
923     EXPECT_EQ(15250, model.latest_state[c].num_observations);
924     EXPECT_NEAR(1, model.latest_state[c].ar_gain, kARGainTolerance);
925 
926     EXPECT_EQ(kNumBlocks,
927               model.combined_state[c].strength_solver.num_equations);
928     EXPECT_EQ(15250, model.combined_state[c].num_observations);
929     EXPECT_NEAR(1, model.combined_state[c].ar_gain, kARGainTolerance);
930   }
931 }
932 
TYPED_TEST_P(NoiseModelUpdateTest,NoiseCoeffsSignalsDifferentNoiseType)933 TYPED_TEST_P(NoiseModelUpdateTest, NoiseCoeffsSignalsDifferentNoiseType) {
934   aom_noise_model_t &model = this->model_;
935   const int width = this->kWidth;
936   const int height = this->kHeight;
937   const double kCoeffs[2][24] = {
938     { 0.02884, -0.03356, 0.00633,  0.01757,  0.02849,  -0.04620,
939       0.02833, -0.07178, 0.07076,  -0.11603, -0.10413, -0.16571,
940       0.05158, -0.07969, 0.02640,  -0.07191, 0.02530,  0.41968,
941       0.21450, -0.00702, -0.01401, -0.03676, -0.08713, 0.44196 },
942     { 0.00269, -0.01291, -0.01513, 0.07234,  0.03208,   0.00477,
943       0.00226, -0.00254, 0.03533,  0.12841,  -0.25970,  -0.06336,
944       0.05238, -0.00845, -0.03118, 0.09043,  -0.36558,  0.48903,
945       0.00595, -0.11938, 0.02106,  0.095956, -0.350139, 0.59305 }
946   };
947 
948   noise_synth(&this->random_, model.params.lag, model.n, model.coords,
949               kCoeffs[0], this->noise_ptr_[0], width, height);
950   for (int i = 0; i < width * height; ++i) {
951     this->data_ptr_[0][i] = (uint8_t)(128 + this->noise_ptr_[0][i]);
952   }
953   this->flat_blocks_.assign(this->flat_blocks_.size(), 1);
954   EXPECT_EQ(AOM_NOISE_STATUS_OK, this->NoiseModelUpdate());
955 
956   // Now try with the second set of AR coefficients
957   noise_synth(&this->random_, model.params.lag, model.n, model.coords,
958               kCoeffs[1], this->noise_ptr_[0], width, height);
959   for (int i = 0; i < width * height; ++i) {
960     this->data_ptr_[0][i] = (uint8_t)(128 + this->noise_ptr_[0][i]);
961   }
962   EXPECT_EQ(AOM_NOISE_STATUS_DIFFERENT_NOISE_TYPE, this->NoiseModelUpdate());
963 }
964 REGISTER_TYPED_TEST_SUITE_P(NoiseModelUpdateTest, UpdateFailsNoFlatBlocks,
965                             UpdateSuccessForZeroNoiseAllFlat,
966                             UpdateFailsBlockSizeTooSmall,
967                             UpdateSuccessForWhiteRandomNoise,
968                             UpdateSuccessForScaledWhiteNoise,
969                             UpdateSuccessForCorrelatedNoise,
970                             NoiseStrengthChangeSignalsDifferentNoiseType,
971                             NoiseCoeffsSignalsDifferentNoiseType);
972 
973 // Note the empty final argument can be removed if C++20 is made the minimum
974 // requirement.
975 INSTANTIATE_TYPED_TEST_SUITE_P(NoiseModelUpdateTestInstatiation,
976                                NoiseModelUpdateTest, AllBitDepthParams, );
977 
TEST(NoiseModelGetGrainParameters,TestLagSize)978 TEST(NoiseModelGetGrainParameters, TestLagSize) {
979   aom_film_grain_t film_grain;
980   for (int lag = 1; lag <= 3; ++lag) {
981     aom_noise_model_params_t params = { AOM_NOISE_SHAPE_SQUARE, lag, 8, 0 };
982     aom_noise_model_t model;
983     EXPECT_TRUE(aom_noise_model_init(&model, params));
984     EXPECT_TRUE(aom_noise_model_get_grain_parameters(&model, &film_grain));
985     EXPECT_EQ(lag, film_grain.ar_coeff_lag);
986     aom_noise_model_free(&model);
987   }
988 
989   aom_noise_model_params_t params = { AOM_NOISE_SHAPE_SQUARE, 4, 8, 0 };
990   aom_noise_model_t model;
991   EXPECT_TRUE(aom_noise_model_init(&model, params));
992   EXPECT_FALSE(aom_noise_model_get_grain_parameters(&model, &film_grain));
993   aom_noise_model_free(&model);
994 }
995 
TEST(NoiseModelGetGrainParameters,TestARCoeffShiftBounds)996 TEST(NoiseModelGetGrainParameters, TestARCoeffShiftBounds) {
997   struct TestCase {
998     double max_input_value;
999     int expected_ar_coeff_shift;
1000     int expected_value;
1001   };
1002   const int lag = 1;
1003   const int kNumTestCases = 19;
1004   const TestCase test_cases[] = {
1005     // Test cases for ar_coeff_shift = 9
1006     { 0, 9, 0 },
1007     { 0.125, 9, 64 },
1008     { -0.125, 9, -64 },
1009     { 0.2499, 9, 127 },
1010     { -0.25, 9, -128 },
1011     // Test cases for ar_coeff_shift = 8
1012     { 0.25, 8, 64 },
1013     { -0.2501, 8, -64 },
1014     { 0.499, 8, 127 },
1015     { -0.5, 8, -128 },
1016     // Test cases for ar_coeff_shift = 7
1017     { 0.5, 7, 64 },
1018     { -0.5001, 7, -64 },
1019     { 0.999, 7, 127 },
1020     { -1, 7, -128 },
1021     // Test cases for ar_coeff_shift = 6
1022     { 1.0, 6, 64 },
1023     { -1.0001, 6, -64 },
1024     { 2.0, 6, 127 },
1025     { -2.0, 6, -128 },
1026     { 4, 6, 127 },
1027     { -4, 6, -128 },
1028   };
1029   aom_noise_model_params_t params = { AOM_NOISE_SHAPE_SQUARE, lag, 8, 0 };
1030   aom_noise_model_t model;
1031   EXPECT_TRUE(aom_noise_model_init(&model, params));
1032 
1033   for (int i = 0; i < kNumTestCases; ++i) {
1034     const TestCase &test_case = test_cases[i];
1035     model.combined_state[0].eqns.x[0] = test_case.max_input_value;
1036 
1037     aom_film_grain_t film_grain;
1038     EXPECT_TRUE(aom_noise_model_get_grain_parameters(&model, &film_grain));
1039     EXPECT_EQ(1, film_grain.ar_coeff_lag);
1040     EXPECT_EQ(test_case.expected_ar_coeff_shift, film_grain.ar_coeff_shift);
1041     EXPECT_EQ(test_case.expected_value, film_grain.ar_coeffs_y[0]);
1042   }
1043   aom_noise_model_free(&model);
1044 }
1045 
TEST(NoiseModelGetGrainParameters,TestNoiseStrengthShiftBounds)1046 TEST(NoiseModelGetGrainParameters, TestNoiseStrengthShiftBounds) {
1047   struct TestCase {
1048     double max_input_value;
1049     int expected_scaling_shift;
1050     int expected_value;
1051   };
1052   const int kNumTestCases = 10;
1053   const TestCase test_cases[] = {
1054     { 0, 11, 0 },      { 1, 11, 64 },     { 2, 11, 128 }, { 3.99, 11, 255 },
1055     { 4, 10, 128 },    { 7.99, 10, 255 }, { 8, 9, 128 },  { 16, 8, 128 },
1056     { 31.99, 8, 255 }, { 64, 8, 255 },  // clipped
1057   };
1058   const int lag = 1;
1059   aom_noise_model_params_t params = { AOM_NOISE_SHAPE_SQUARE, lag, 8, 0 };
1060   aom_noise_model_t model;
1061   EXPECT_TRUE(aom_noise_model_init(&model, params));
1062 
1063   for (int i = 0; i < kNumTestCases; ++i) {
1064     const TestCase &test_case = test_cases[i];
1065     aom_equation_system_t &eqns = model.combined_state[0].strength_solver.eqns;
1066     // Set the fitted scale parameters to be a constant value.
1067     for (int j = 0; j < eqns.n; ++j) {
1068       eqns.x[j] = test_case.max_input_value;
1069     }
1070     aom_film_grain_t film_grain;
1071     EXPECT_TRUE(aom_noise_model_get_grain_parameters(&model, &film_grain));
1072     // We expect a single constant segemnt
1073     EXPECT_EQ(test_case.expected_scaling_shift, film_grain.scaling_shift);
1074     EXPECT_EQ(test_case.expected_value, film_grain.scaling_points_y[0][1]);
1075     EXPECT_EQ(test_case.expected_value, film_grain.scaling_points_y[1][1]);
1076   }
1077   aom_noise_model_free(&model);
1078 }
1079 
1080 // The AR coefficients are the same inputs used to generate "Test 2" in the test
1081 // vectors
TEST(NoiseModelGetGrainParameters,GetGrainParametersReal)1082 TEST(NoiseModelGetGrainParameters, GetGrainParametersReal) {
1083   const double kInputCoeffsY[] = { 0.0315,  0.0073,  0.0218,  0.00235, 0.00511,
1084                                    -0.0222, 0.0627,  -0.022,  0.05575, -0.1816,
1085                                    0.0107,  -0.1966, 0.00065, -0.0809, 0.04934,
1086                                    -0.1349, -0.0352, 0.41772, 0.27973, 0.04207,
1087                                    -0.0429, -0.1372, 0.06193, 0.52032 };
1088   const double kInputCoeffsCB[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  0,
1089                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5 };
1090   const double kInputCoeffsCR[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   0,
1091                                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.5 };
1092   const int kExpectedARCoeffsY[] = { 4,  1,   3,  0,   1,  -3,  8, -3,
1093                                      7,  -23, 1,  -25, 0,  -10, 6, -17,
1094                                      -5, 53,  36, 5,   -5, -18, 8, 67 };
1095   const int kExpectedARCoeffsCB[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
1096                                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 84 };
1097   const int kExpectedARCoeffsCR[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   0,
1098                                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -126 };
1099   // Scaling function is initialized analytically with a sqrt function.
1100   const int kNumScalingPointsY = 12;
1101   const int kExpectedScalingPointsY[][2] = {
1102     { 0, 0 },     { 13, 44 },   { 27, 62 },   { 40, 76 },
1103     { 54, 88 },   { 67, 98 },   { 94, 117 },  { 121, 132 },
1104     { 148, 146 }, { 174, 159 }, { 201, 171 }, { 255, 192 },
1105   };
1106 
1107   const int lag = 3;
1108   aom_noise_model_params_t params = { AOM_NOISE_SHAPE_SQUARE, lag, 8, 0 };
1109   aom_noise_model_t model;
1110   EXPECT_TRUE(aom_noise_model_init(&model, params));
1111 
1112   // Setup the AR coeffs
1113   memcpy(model.combined_state[0].eqns.x, kInputCoeffsY, sizeof(kInputCoeffsY));
1114   memcpy(model.combined_state[1].eqns.x, kInputCoeffsCB,
1115          sizeof(kInputCoeffsCB));
1116   memcpy(model.combined_state[2].eqns.x, kInputCoeffsCR,
1117          sizeof(kInputCoeffsCR));
1118   for (int i = 0; i < model.combined_state[0].strength_solver.num_bins; ++i) {
1119     const double x =
1120         ((double)i) / (model.combined_state[0].strength_solver.num_bins - 1.0);
1121     model.combined_state[0].strength_solver.eqns.x[i] = 6 * sqrt(x);
1122     model.combined_state[1].strength_solver.eqns.x[i] = 3;
1123     model.combined_state[2].strength_solver.eqns.x[i] = 2;
1124 
1125     // Inject some observations into the strength solver, as during film grain
1126     // parameter extraction an estimate of the average strength will be used to
1127     // adjust correlation.
1128     const int n = model.combined_state[0].strength_solver.num_bins;
1129     for (int j = 0; j < model.combined_state[0].strength_solver.num_bins; ++j) {
1130       model.combined_state[0].strength_solver.eqns.A[i * n + j] = 1;
1131       model.combined_state[1].strength_solver.eqns.A[i * n + j] = 1;
1132       model.combined_state[2].strength_solver.eqns.A[i * n + j] = 1;
1133     }
1134   }
1135 
1136   aom_film_grain_t film_grain;
1137   EXPECT_TRUE(aom_noise_model_get_grain_parameters(&model, &film_grain));
1138   EXPECT_EQ(lag, film_grain.ar_coeff_lag);
1139   EXPECT_EQ(3, film_grain.ar_coeff_lag);
1140   EXPECT_EQ(7, film_grain.ar_coeff_shift);
1141   EXPECT_EQ(10, film_grain.scaling_shift);
1142   EXPECT_EQ(kNumScalingPointsY, film_grain.num_y_points);
1143   EXPECT_EQ(1, film_grain.update_parameters);
1144   EXPECT_EQ(1, film_grain.apply_grain);
1145 
1146   const int kNumARCoeffs = 24;
1147   for (int i = 0; i < kNumARCoeffs; ++i) {
1148     EXPECT_EQ(kExpectedARCoeffsY[i], film_grain.ar_coeffs_y[i]);
1149   }
1150   for (int i = 0; i < kNumARCoeffs + 1; ++i) {
1151     EXPECT_EQ(kExpectedARCoeffsCB[i], film_grain.ar_coeffs_cb[i]);
1152   }
1153   for (int i = 0; i < kNumARCoeffs + 1; ++i) {
1154     EXPECT_EQ(kExpectedARCoeffsCR[i], film_grain.ar_coeffs_cr[i]);
1155   }
1156   for (int i = 0; i < kNumScalingPointsY; ++i) {
1157     EXPECT_EQ(kExpectedScalingPointsY[i][0], film_grain.scaling_points_y[i][0]);
1158     EXPECT_EQ(kExpectedScalingPointsY[i][1], film_grain.scaling_points_y[i][1]);
1159   }
1160 
1161   // CB strength should just be a piecewise segment
1162   EXPECT_EQ(2, film_grain.num_cb_points);
1163   EXPECT_EQ(0, film_grain.scaling_points_cb[0][0]);
1164   EXPECT_EQ(255, film_grain.scaling_points_cb[1][0]);
1165   EXPECT_EQ(96, film_grain.scaling_points_cb[0][1]);
1166   EXPECT_EQ(96, film_grain.scaling_points_cb[1][1]);
1167 
1168   // CR strength should just be a piecewise segment
1169   EXPECT_EQ(2, film_grain.num_cr_points);
1170   EXPECT_EQ(0, film_grain.scaling_points_cr[0][0]);
1171   EXPECT_EQ(255, film_grain.scaling_points_cr[1][0]);
1172   EXPECT_EQ(64, film_grain.scaling_points_cr[0][1]);
1173   EXPECT_EQ(64, film_grain.scaling_points_cr[1][1]);
1174 
1175   EXPECT_EQ(128, film_grain.cb_mult);
1176   EXPECT_EQ(192, film_grain.cb_luma_mult);
1177   EXPECT_EQ(256, film_grain.cb_offset);
1178   EXPECT_EQ(128, film_grain.cr_mult);
1179   EXPECT_EQ(192, film_grain.cr_luma_mult);
1180   EXPECT_EQ(256, film_grain.cr_offset);
1181   EXPECT_EQ(0, film_grain.chroma_scaling_from_luma);
1182   EXPECT_EQ(0, film_grain.grain_scale_shift);
1183 
1184   aom_noise_model_free(&model);
1185 }
1186 
1187 template <typename T>
1188 class WienerDenoiseTest : public ::testing::Test, public T {
1189  public:
SetUpTestSuite()1190   static void SetUpTestSuite() { aom_dsp_rtcd(); }
1191 
1192  protected:
SetUp()1193   void SetUp() override {
1194     static const float kNoiseLevel = 5.f;
1195     static const float kStd = 4.0;
1196     static const double kMaxValue = (1 << T::kBitDepth) - 1;
1197 
1198     chroma_sub_[0] = 1;
1199     chroma_sub_[1] = 1;
1200     stride_[0] = kWidth;
1201     stride_[1] = kWidth / 2;
1202     stride_[2] = kWidth / 2;
1203     for (int k = 0; k < 3; ++k) {
1204       data_[k].resize(kWidth * kHeight);
1205       denoised_[k].resize(kWidth * kHeight);
1206       noise_psd_[k].resize(kBlockSize * kBlockSize);
1207     }
1208 
1209     const double kCoeffsY[] = { 0.0406, -0.116, -0.078, -0.152, 0.0033, -0.093,
1210                                 0.048,  0.404,  0.2353, -0.035, -0.093, 0.441 };
1211     const int kCoords[12][2] = {
1212       { -2, -2 }, { -1, -2 }, { 0, -2 }, { 1, -2 }, { 2, -2 }, { -2, -1 },
1213       { -1, -1 }, { 0, -1 },  { 1, -1 }, { 2, -1 }, { -2, 0 }, { -1, 0 }
1214     };
1215     const int kLag = 2;
1216     const int kLength = 12;
1217     libaom_test::ACMRandom random;
1218     std::vector<double> noise(kWidth * kHeight);
1219     noise_synth(&random, kLag, kLength, kCoords, kCoeffsY, &noise[0], kWidth,
1220                 kHeight);
1221     noise_psd_[0] = get_noise_psd(&noise[0], kWidth, kHeight, kBlockSize);
1222     for (int i = 0; i < kBlockSize * kBlockSize; ++i) {
1223       noise_psd_[0][i] = (float)(noise_psd_[0][i] * kStd * kStd * kScaleNoise *
1224                                  kScaleNoise / (kMaxValue * kMaxValue));
1225     }
1226 
1227     float psd_value =
1228         aom_noise_psd_get_default_value(kBlockSizeChroma, kNoiseLevel);
1229     for (int i = 0; i < kBlockSizeChroma * kBlockSizeChroma; ++i) {
1230       noise_psd_[1][i] = psd_value;
1231       noise_psd_[2][i] = psd_value;
1232     }
1233     for (int y = 0; y < kHeight; ++y) {
1234       for (int x = 0; x < kWidth; ++x) {
1235         data_[0][y * stride_[0] + x] = (typename T::data_type_t)fclamp(
1236             (x + noise[y * stride_[0] + x] * kStd) * kScaleNoise, 0, kMaxValue);
1237       }
1238     }
1239 
1240     for (int c = 1; c < 3; ++c) {
1241       for (int y = 0; y < (kHeight >> 1); ++y) {
1242         for (int x = 0; x < (kWidth >> 1); ++x) {
1243           data_[c][y * stride_[c] + x] = (typename T::data_type_t)fclamp(
1244               (x + randn(&random, kStd)) * kScaleNoise, 0, kMaxValue);
1245         }
1246       }
1247     }
1248     for (int k = 0; k < 3; ++k) {
1249       noise_psd_ptrs_[k] = &noise_psd_[k][0];
1250     }
1251   }
1252   static const int kBlockSize = 32;
1253   static const int kBlockSizeChroma = 16;
1254   static const int kWidth = 256;
1255   static const int kHeight = 256;
1256   static const int kScaleNoise = 1 << (T::kBitDepth - 8);
1257 
1258   std::vector<typename T::data_type_t> data_[3];
1259   std::vector<typename T::data_type_t> denoised_[3];
1260   std::vector<float> noise_psd_[3];
1261   int chroma_sub_[2];
1262   float *noise_psd_ptrs_[3];
1263   int stride_[3];
1264 };
1265 
1266 TYPED_TEST_SUITE_P(WienerDenoiseTest);
1267 
TYPED_TEST_P(WienerDenoiseTest,InvalidBlockSize)1268 TYPED_TEST_P(WienerDenoiseTest, InvalidBlockSize) {
1269   const uint8_t *const data_ptrs[3] = {
1270     reinterpret_cast<uint8_t *>(&this->data_[0][0]),
1271     reinterpret_cast<uint8_t *>(&this->data_[1][0]),
1272     reinterpret_cast<uint8_t *>(&this->data_[2][0]),
1273   };
1274   uint8_t *denoised_ptrs[3] = {
1275     reinterpret_cast<uint8_t *>(&this->denoised_[0][0]),
1276     reinterpret_cast<uint8_t *>(&this->denoised_[1][0]),
1277     reinterpret_cast<uint8_t *>(&this->denoised_[2][0]),
1278   };
1279   EXPECT_EQ(0, aom_wiener_denoise_2d(data_ptrs, denoised_ptrs, this->kWidth,
1280                                      this->kHeight, this->stride_,
1281                                      this->chroma_sub_, this->noise_psd_ptrs_,
1282                                      18, this->kBitDepth, this->kUseHighBD));
1283   EXPECT_EQ(0, aom_wiener_denoise_2d(data_ptrs, denoised_ptrs, this->kWidth,
1284                                      this->kHeight, this->stride_,
1285                                      this->chroma_sub_, this->noise_psd_ptrs_,
1286                                      48, this->kBitDepth, this->kUseHighBD));
1287   EXPECT_EQ(0, aom_wiener_denoise_2d(data_ptrs, denoised_ptrs, this->kWidth,
1288                                      this->kHeight, this->stride_,
1289                                      this->chroma_sub_, this->noise_psd_ptrs_,
1290                                      64, this->kBitDepth, this->kUseHighBD));
1291 }
1292 
TYPED_TEST_P(WienerDenoiseTest,InvalidChromaSubsampling)1293 TYPED_TEST_P(WienerDenoiseTest, InvalidChromaSubsampling) {
1294   const uint8_t *const data_ptrs[3] = {
1295     reinterpret_cast<uint8_t *>(&this->data_[0][0]),
1296     reinterpret_cast<uint8_t *>(&this->data_[1][0]),
1297     reinterpret_cast<uint8_t *>(&this->data_[2][0]),
1298   };
1299   uint8_t *denoised_ptrs[3] = {
1300     reinterpret_cast<uint8_t *>(&this->denoised_[0][0]),
1301     reinterpret_cast<uint8_t *>(&this->denoised_[1][0]),
1302     reinterpret_cast<uint8_t *>(&this->denoised_[2][0]),
1303   };
1304   int chroma_sub[2] = { 1, 0 };
1305   EXPECT_EQ(0, aom_wiener_denoise_2d(data_ptrs, denoised_ptrs, this->kWidth,
1306                                      this->kHeight, this->stride_, chroma_sub,
1307                                      this->noise_psd_ptrs_, 32, this->kBitDepth,
1308                                      this->kUseHighBD));
1309 
1310   chroma_sub[0] = 0;
1311   chroma_sub[1] = 1;
1312   EXPECT_EQ(0, aom_wiener_denoise_2d(data_ptrs, denoised_ptrs, this->kWidth,
1313                                      this->kHeight, this->stride_, chroma_sub,
1314                                      this->noise_psd_ptrs_, 32, this->kBitDepth,
1315                                      this->kUseHighBD));
1316 }
1317 
TYPED_TEST_P(WienerDenoiseTest,GradientTest)1318 TYPED_TEST_P(WienerDenoiseTest, GradientTest) {
1319   const int width = this->kWidth;
1320   const int height = this->kHeight;
1321   const int block_size = this->kBlockSize;
1322   const uint8_t *const data_ptrs[3] = {
1323     reinterpret_cast<uint8_t *>(&this->data_[0][0]),
1324     reinterpret_cast<uint8_t *>(&this->data_[1][0]),
1325     reinterpret_cast<uint8_t *>(&this->data_[2][0]),
1326   };
1327   uint8_t *denoised_ptrs[3] = {
1328     reinterpret_cast<uint8_t *>(&this->denoised_[0][0]),
1329     reinterpret_cast<uint8_t *>(&this->denoised_[1][0]),
1330     reinterpret_cast<uint8_t *>(&this->denoised_[2][0]),
1331   };
1332   const int ret = aom_wiener_denoise_2d(
1333       data_ptrs, denoised_ptrs, width, height, this->stride_, this->chroma_sub_,
1334       this->noise_psd_ptrs_, block_size, this->kBitDepth, this->kUseHighBD);
1335   EXPECT_EQ(1, ret);
1336 
1337   // Check the noise on the denoised image (from the analytical gradient)
1338   // and make sure that it is less than what we added.
1339   for (int c = 0; c < 3; ++c) {
1340     std::vector<double> measured_noise(width * height);
1341 
1342     double var = 0;
1343     const int shift = (c > 0);
1344     for (int x = 0; x < (width >> shift); ++x) {
1345       for (int y = 0; y < (height >> shift); ++y) {
1346         const double diff = this->denoised_[c][y * this->stride_[c] + x] -
1347                             x * this->kScaleNoise;
1348         var += diff * diff;
1349         measured_noise[y * width + x] = diff;
1350       }
1351     }
1352     var /= (width * height);
1353     const double std = sqrt(std::max(0.0, var));
1354     EXPECT_LE(std, 1.25f * this->kScaleNoise);
1355     if (c == 0) {
1356       std::vector<float> measured_psd =
1357           get_noise_psd(&measured_noise[0], width, height, block_size);
1358       std::vector<double> measured_psd_d(block_size * block_size);
1359       std::vector<double> noise_psd_d(block_size * block_size);
1360       std::copy(measured_psd.begin(), measured_psd.end(),
1361                 measured_psd_d.begin());
1362       std::copy(this->noise_psd_[0].begin(), this->noise_psd_[0].end(),
1363                 noise_psd_d.begin());
1364       EXPECT_LT(
1365           aom_normalized_cross_correlation(&measured_psd_d[0], &noise_psd_d[0],
1366                                            (int)(noise_psd_d.size())),
1367           0.35);
1368     }
1369   }
1370 }
1371 
1372 REGISTER_TYPED_TEST_SUITE_P(WienerDenoiseTest, InvalidBlockSize,
1373                             InvalidChromaSubsampling, GradientTest);
1374 
1375 // Note the empty final argument can be removed if C++20 is made the minimum
1376 // requirement.
1377 INSTANTIATE_TYPED_TEST_SUITE_P(WienerDenoiseTestInstatiation, WienerDenoiseTest,
1378                                AllBitDepthParams, );
1379