diff --git a/rtc_tools/BUILD.gn b/rtc_tools/BUILD.gn index c0716b4dac..3d9804c613 100644 --- a/rtc_tools/BUILD.gn +++ b/rtc_tools/BUILD.gn @@ -86,6 +86,10 @@ rtc_static_library("video_file_writer") { rtc_static_library("video_quality_analysis") { sources = [ + "frame_analyzer/linear_least_squares.cc", + "frame_analyzer/linear_least_squares.h", + "frame_analyzer/video_color_aligner.cc", + "frame_analyzer/video_color_aligner.h", "frame_analyzer/video_quality_analysis.cc", "frame_analyzer/video_quality_analysis.h", "frame_analyzer/video_temporal_aligner.cc", @@ -93,6 +97,7 @@ rtc_static_library("video_quality_analysis") { ] deps = [ ":video_file_reader", + "../api:array_view", "../api/video:video_frame_i420", "../common_video", "../rtc_base:checks", @@ -333,7 +338,9 @@ if (rtc_include_tests) { testonly = true sources = [ + "frame_analyzer/linear_least_squares_unittest.cc", "frame_analyzer/reference_less_video_analysis_unittest.cc", + "frame_analyzer/video_color_aligner_unittest.cc", "frame_analyzer/video_quality_analysis_unittest.cc", "frame_analyzer/video_temporal_aligner_unittest.cc", "frame_editing/frame_editing_unittest.cc", @@ -363,6 +370,7 @@ if (rtc_include_tests) { "//test:test_support", "//testing/gtest", "//third_party/abseil-cpp/absl/memory", + "//third_party/libyuv", ] if (rtc_enable_protobuf) { diff --git a/rtc_tools/frame_analyzer/linear_least_squares.cc b/rtc_tools/frame_analyzer/linear_least_squares.cc new file mode 100644 index 0000000000..706d212d17 --- /dev/null +++ b/rtc_tools/frame_analyzer/linear_least_squares.cc @@ -0,0 +1,200 @@ +/* + * Copyright (c) 2018 The WebRTC project authors. All Rights Reserved. + * + * Use of this source code is governed by a BSD-style license + * that can be found in the LICENSE file in the root of the source + * tree. An additional intellectual property rights grant can be found + * in the file PATENTS. All contributing project authors may + * be found in the AUTHORS file in the root of the source tree. + */ + +#include "rtc_tools/frame_analyzer/linear_least_squares.h" + +#include +#include + +#include "rtc_base/checks.h" +#include "rtc_base/logging.h" + +namespace webrtc { +namespace test { + +template +using Matrix = std::valarray>; + +namespace { + +template +R DotProduct(const std::valarray& a, const std::valarray& b) { + RTC_CHECK_EQ(a.size(), b.size()); + return std::inner_product(std::begin(a), std::end(a), std::begin(b), R(0)); +} + +// Calculates a^T * b. +template +Matrix MatrixMultiply(const Matrix& a, const Matrix& b) { + Matrix result(std::valarray(a.size()), b.size()); + for (size_t i = 0; i < a.size(); ++i) { + for (size_t j = 0; j < b.size(); ++j) + result[j][i] = DotProduct(a[i], b[j]); + } + + return result; +} + +template +Matrix Transpose(const Matrix& matrix) { + if (matrix.size() == 0) + return Matrix(); + const size_t rows = matrix.size(); + const size_t columns = matrix[0].size(); + Matrix result(std::valarray(rows), columns); + + for (size_t i = 0; i < rows; ++i) { + for (size_t j = 0; j < columns; ++j) + result[j][i] = matrix[i][j]; + } + + return result; +} + +// Convert valarray from type T to type R. +template +std::valarray ConvertTo(const std::valarray& v) { + std::valarray result(v.size()); + for (size_t i = 0; i < v.size(); ++i) + result[i] = static_cast(v[i]); + return result; +} + +// Convert valarray Matrix from type T to type R. +template +Matrix ConvertTo(const Matrix& mat) { + Matrix result(mat.size()); + for (size_t i = 0; i < mat.size(); ++i) + result[i] = ConvertTo(mat[i]); + return result; +} + +// Convert from valarray Matrix back to the more conventional std::vector. +template +std::vector> ToVectorMatrix(const Matrix& m) { + std::vector> result; + for (const std::valarray& v : m) + result.emplace_back(std::begin(v), std::end(v)); + return result; +} + +// Create a valarray Matrix from a conventional std::vector. +template +Matrix FromVectorMatrix(const std::vector>& mat) { + Matrix result(mat.size()); + for (size_t i = 0; i < mat.size(); ++i) + result[i] = std::valarray(mat[i].data(), mat[i].size()); + return result; +} + +// Returns |matrix_to_invert|^-1 * |right_hand_matrix|. |matrix_to_invert| must +// have square size. +Matrix GaussianElimination(Matrix matrix_to_invert, + Matrix right_hand_matrix) { + // |n| is the width/height of |matrix_to_invert|. + const size_t n = matrix_to_invert.size(); + // Make sure |matrix_to_invert| has square size. + for (const std::valarray& column : matrix_to_invert) + RTC_CHECK_EQ(n, column.size()); + // Make sure |right_hand_matrix| has correct size. + for (const std::valarray& column : right_hand_matrix) + RTC_CHECK_EQ(n, column.size()); + + // Transpose the matrices before and after so that we can perform Gaussian + // elimination on the columns instead of the rows, since that is easier with + // our representation. + matrix_to_invert = Transpose(matrix_to_invert); + right_hand_matrix = Transpose(right_hand_matrix); + + // Loop over the diagonal of |matrix_to_invert| and perform column reduction. + // Column reduction is a sequence of elementary column operations that is + // performed on both |matrix_to_invert| and |right_hand_matrix| until + // |matrix_to_invert| has been transformed to the identity matrix. + for (size_t diagonal_index = 0; diagonal_index < n; ++diagonal_index) { + // Make sure the diagonal element has the highest absolute value by + // swapping columns if necessary. + for (size_t column = diagonal_index + 1; column < n; ++column) { + if (std::abs(matrix_to_invert[column][diagonal_index]) > + std::abs(matrix_to_invert[diagonal_index][diagonal_index])) { + std::swap(matrix_to_invert[column], matrix_to_invert[diagonal_index]); + std::swap(right_hand_matrix[column], right_hand_matrix[diagonal_index]); + } + } + + // Reduce the diagonal element to be 1, by dividing the column with that + // value. If the diagonal element is 0, it means the system of equations has + // many solutions, and in that case we will return an arbitrary solution. + if (matrix_to_invert[diagonal_index][diagonal_index] == 0.0) { + RTC_LOG(LS_WARNING) << "Matrix is not invertible, ignoring."; + continue; + } + const double diagonal_element = + matrix_to_invert[diagonal_index][diagonal_index]; + matrix_to_invert[diagonal_index] /= diagonal_element; + right_hand_matrix[diagonal_index] /= diagonal_element; + + // Eliminate the other entries in row |diagonal_index| by making them zero. + for (size_t column = 0; column < n; ++column) { + if (column == diagonal_index) + continue; + const double row_element = matrix_to_invert[column][diagonal_index]; + matrix_to_invert[column] -= + row_element * matrix_to_invert[diagonal_index]; + right_hand_matrix[column] -= + row_element * right_hand_matrix[diagonal_index]; + } + } + + // Transpose the result before returning it, explained in comment above. + return Transpose(right_hand_matrix); +} + +} // namespace + +IncrementalLinearLeastSquares::IncrementalLinearLeastSquares() = default; +IncrementalLinearLeastSquares::~IncrementalLinearLeastSquares() = default; + +void IncrementalLinearLeastSquares::AddObservations( + const std::vector>& x, + const std::vector>& y) { + if (x.empty() || y.empty()) + return; + // Make sure all columns are the same size. + const size_t n = x[0].size(); + for (const std::vector& column : x) + RTC_CHECK_EQ(n, column.size()); + for (const std::vector& column : y) + RTC_CHECK_EQ(n, column.size()); + + // We will multiply the uint8_t values together, so we need to expand to a + // type that can safely store those values, i.e. uint16_t. + const Matrix unpacked_x = ConvertTo(FromVectorMatrix(x)); + const Matrix unpacked_y = ConvertTo(FromVectorMatrix(y)); + + const Matrix xx = MatrixMultiply(unpacked_x, unpacked_x); + const Matrix xy = MatrixMultiply(unpacked_x, unpacked_y); + if (sum_xx && sum_xy) { + *sum_xx += xx; + *sum_xy += xy; + } else { + sum_xx = xx; + sum_xy = xy; + } +} + +std::vector> +IncrementalLinearLeastSquares::GetBestSolution() const { + RTC_CHECK(sum_xx && sum_xy) << "No observations have been added"; + return ToVectorMatrix(GaussianElimination(ConvertTo(*sum_xx), + ConvertTo(*sum_xy))); +} + +} // namespace test +} // namespace webrtc diff --git a/rtc_tools/frame_analyzer/linear_least_squares.h b/rtc_tools/frame_analyzer/linear_least_squares.h new file mode 100644 index 0000000000..1b07dc1bde --- /dev/null +++ b/rtc_tools/frame_analyzer/linear_least_squares.h @@ -0,0 +1,53 @@ +/* + * Copyright (c) 2018 The WebRTC project authors. All Rights Reserved. + * + * Use of this source code is governed by a BSD-style license + * that can be found in the LICENSE file in the root of the source + * tree. An additional intellectual property rights grant can be found + * in the file PATENTS. All contributing project authors may + * be found in the AUTHORS file in the root of the source tree. + */ + +#ifndef RTC_TOOLS_FRAME_ANALYZER_LINEAR_LEAST_SQUARES_H_ +#define RTC_TOOLS_FRAME_ANALYZER_LINEAR_LEAST_SQUARES_H_ + +#include +#include +#include "absl/types/optional.h" + +namespace webrtc { +namespace test { + +// This class is used for finding a matrix b that roughly solves the equation: +// y = x * b. This is generally impossible to do exactly, so the problem is +// rephrased as finding the matrix b that minimizes the difference: +// |y - x * b|^2. Calling multiple AddObservations() is equivalent to +// concatenating the observation vectors and calling AddObservations() once. The +// reason for doing it incrementally is that we can't store the raw YUV values +// for a whole video file in memory at once. This class has a constant memory +// footprint, regardless how may times AddObservations() is called. +class IncrementalLinearLeastSquares { + public: + IncrementalLinearLeastSquares(); + ~IncrementalLinearLeastSquares(); + + // Add a number of observations. The subvectors of x and y must have the same + // length. + void AddObservations(const std::vector>& x, + const std::vector>& y); + + // Calculate and return the best linear solution, given the observations so + // far. + std::vector> GetBestSolution() const; + + private: + // Running sum of x^T * x. + absl::optional>> sum_xx; + // Running sum of x^T * y. + absl::optional>> sum_xy; +}; + +} // namespace test +} // namespace webrtc + +#endif // RTC_TOOLS_FRAME_ANALYZER_LINEAR_LEAST_SQUARES_H_ diff --git a/rtc_tools/frame_analyzer/linear_least_squares_unittest.cc b/rtc_tools/frame_analyzer/linear_least_squares_unittest.cc new file mode 100644 index 0000000000..b074aacce8 --- /dev/null +++ b/rtc_tools/frame_analyzer/linear_least_squares_unittest.cc @@ -0,0 +1,91 @@ +/* + * Copyright (c) 2018 The WebRTC project authors. All Rights Reserved. + * + * Use of this source code is governed by a BSD-style license + * that can be found in the LICENSE file in the root of the source + * tree. An additional intellectual property rights grant can be found + * in the file PATENTS. All contributing project authors may + * be found in the AUTHORS file in the root of the source tree. + */ + +#include "rtc_tools/frame_analyzer/linear_least_squares.h" + +#include "test/gtest.h" + +namespace webrtc { +namespace test { + +TEST(LinearLeastSquares, ScalarIdentityOneObservation) { + IncrementalLinearLeastSquares lls; + lls.AddObservations({{1}}, {{1}}); + EXPECT_EQ(std::vector>({{1.0}}), lls.GetBestSolution()); +} + +TEST(LinearLeastSquares, ScalarIdentityTwoObservationsOneCall) { + IncrementalLinearLeastSquares lls; + lls.AddObservations({{1, 2}}, {{1, 2}}); + EXPECT_EQ(std::vector>({{1.0}}), lls.GetBestSolution()); +} + +TEST(LinearLeastSquares, ScalarIdentityTwoObservationsTwoCalls) { + IncrementalLinearLeastSquares lls; + lls.AddObservations({{1}}, {{1}}); + lls.AddObservations({{2}}, {{2}}); + EXPECT_EQ(std::vector>({{1.0}}), lls.GetBestSolution()); +} + +TEST(LinearLeastSquares, MatrixIdentityOneObservation) { + IncrementalLinearLeastSquares lls; + lls.AddObservations({{1, 2}, {3, 4}}, {{1, 2}, {3, 4}}); + EXPECT_EQ(std::vector>({{1.0, 0.0}, {0.0, 1.0}}), + lls.GetBestSolution()); +} + +TEST(LinearLeastSquares, MatrixManyObservations) { + IncrementalLinearLeastSquares lls; + // Test that we can find the solution of the overspecified equation system: + // [1, 2] [1, 3] = [5, 11] + // [3, 4] [2, 4] [11, 25] + // [5, 6] [17, 39] + lls.AddObservations({{1}, {2}}, {{5}, {11}}); + lls.AddObservations({{3}, {4}}, {{11}, {25}}); + lls.AddObservations({{5}, {6}}, {{17}, {39}}); + + const std::vector> result = lls.GetBestSolution(); + // We allow some numerical flexibility here. + EXPECT_DOUBLE_EQ(1.0, result[0][0]); + EXPECT_DOUBLE_EQ(2.0, result[0][1]); + EXPECT_DOUBLE_EQ(3.0, result[1][0]); + EXPECT_DOUBLE_EQ(4.0, result[1][1]); +} + +TEST(LinearLeastSquares, MatrixVectorOneObservation) { + IncrementalLinearLeastSquares lls; + // Test that we can find the solution of the overspecified equation system: + // [1, 2] [1] = [5] + // [3, 4] [2] [11] + // [5, 6] [17] + lls.AddObservations({{1, 3, 5}, {2, 4, 6}}, {{5, 11, 17}}); + + const std::vector> result = lls.GetBestSolution(); + // We allow some numerical flexibility here. + EXPECT_DOUBLE_EQ(1.0, result[0][0]); + EXPECT_DOUBLE_EQ(2.0, result[0][1]); +} + +TEST(LinearLeastSquares, LinearLeastSquaresNonPerfectSolution) { + IncrementalLinearLeastSquares lls; + // Test that we can find the non-perfect solution of the overspecified + // equation system: + // [1] [20] = [21] + // [2] [39] + // [3] [60] + // [2] [41] + // [1] [19] + lls.AddObservations({{1, 2, 3, 2, 1}}, {{21, 39, 60, 41, 19}}); + + EXPECT_DOUBLE_EQ(20.0, lls.GetBestSolution()[0][0]); +} + +} // namespace test +} // namespace webrtc diff --git a/rtc_tools/frame_analyzer/video_color_aligner.cc b/rtc_tools/frame_analyzer/video_color_aligner.cc new file mode 100644 index 0000000000..7afb1e4695 --- /dev/null +++ b/rtc_tools/frame_analyzer/video_color_aligner.cc @@ -0,0 +1,239 @@ +/* + * Copyright (c) 2018 The WebRTC project authors. All Rights Reserved. + * + * Use of this source code is governed by a BSD-style license + * that can be found in the LICENSE file in the root of the source + * tree. An additional intellectual property rights grant can be found + * in the file PATENTS. All contributing project authors may + * be found in the AUTHORS file in the root of the source tree. + */ + +#include "rtc_tools/frame_analyzer/video_color_aligner.h" + +#include +#include +#include +#include +#include + +#include "api/array_view.h" +#include "api/video/i420_buffer.h" +#include "rtc_base/checks.h" +#include "rtc_base/logging.h" +#include "rtc_base/refcountedobject.h" +#include "rtc_tools/frame_analyzer/linear_least_squares.h" +#include "rtc_tools/frame_analyzer/video_quality_analysis.h" +#include "third_party/libyuv/include/libyuv/compare.h" +#include "third_party/libyuv/include/libyuv/planar_functions.h" +#include "third_party/libyuv/include/libyuv/scale.h" + +namespace webrtc { +namespace test { + +namespace { + +// Helper function for AdjustColors(). This functions calculates a single output +// row for y with the given color coefficients. The u/v channels are assumed to +// be subsampled by a factor of 2, which is the case of I420. +void CalculateYChannel(rtc::ArrayView y_data, + rtc::ArrayView u_data, + rtc::ArrayView v_data, + const std::array& coeff, + rtc::ArrayView output) { + RTC_CHECK_EQ(y_data.size(), output.size()); + // Each u/v element represents two y elements. Make sure we have enough to + // cover the Y values. + RTC_CHECK_GE(u_data.size() * 2, y_data.size()); + RTC_CHECK_GE(v_data.size() * 2, y_data.size()); + + // Do two pixels at a time since u/v are subsampled. + for (size_t i = 0; i * 2 < y_data.size() - 1; ++i) { + const float uv_contribution = + coeff[1] * u_data[i] + coeff[2] * v_data[i] + coeff[3]; + + const float val0 = coeff[0] * y_data[i * 2 + 0] + uv_contribution; + const float val1 = coeff[0] * y_data[i * 2 + 1] + uv_contribution; + + // Clamp result to a byte. + output[i * 2 + 0] = static_cast( + std::round(std::max(0.0f, std::min(val0, 255.0f)))); + output[i * 2 + 1] = static_cast( + std::round(std::max(0.0f, std::min(val1, 255.0f)))); + } + + // Handle the last pixel for odd widths. + if (y_data.size() % 2 == 1) { + const float val = coeff[0] * y_data[y_data.size() - 1] + + coeff[1] * u_data[(y_data.size() - 1) / 2] + + coeff[2] * v_data[(y_data.size() - 1) / 2] + coeff[3]; + output[y_data.size() - 1] = + static_cast(std::round(std::max(0.0f, std::min(val, 255.0f)))); + } +} + +// Helper function for AdjustColors(). This functions calculates a single output +// row for either u or v, with the given color coefficients. Y, U, and V are +// assumed to be the same size, i.e. no subsampling. +void CalculateUVChannel(rtc::ArrayView y_data, + rtc::ArrayView u_data, + rtc::ArrayView v_data, + const std::array& coeff, + rtc::ArrayView output) { + RTC_CHECK_EQ(y_data.size(), u_data.size()); + RTC_CHECK_EQ(y_data.size(), v_data.size()); + RTC_CHECK_EQ(y_data.size(), output.size()); + + for (size_t x = 0; x < y_data.size(); ++x) { + const float val = coeff[0] * y_data[x] + coeff[1] * u_data[x] + + coeff[2] * v_data[x] + coeff[3]; + // Clamp result to a byte. + output[x] = + static_cast(std::round(std::max(0.0f, std::min(val, 255.0f)))); + } +} + +// Convert a frame to four vectors consisting of [y, u, v, 1]. +std::vector> FlattenYuvData( + const rtc::scoped_refptr& frame) { + std::vector> result( + 4, std::vector(frame->ChromaWidth() * frame->ChromaHeight())); + + // Downscale the Y plane so that all YUV planes are the same size. + libyuv::ScalePlane(frame->DataY(), frame->StrideY(), frame->width(), + frame->height(), result[0].data(), frame->ChromaWidth(), + frame->ChromaWidth(), frame->ChromaHeight(), + libyuv::kFilterBox); + + libyuv::CopyPlane(frame->DataU(), frame->StrideU(), result[1].data(), + frame->ChromaWidth(), frame->ChromaWidth(), + frame->ChromaHeight()); + + libyuv::CopyPlane(frame->DataV(), frame->StrideV(), result[2].data(), + frame->ChromaWidth(), frame->ChromaWidth(), + frame->ChromaHeight()); + + std::fill(result[3].begin(), result[3].end(), 1u); + + return result; +} + +ColorTransformationMatrix VectorToColorMatrix( + const std::vector>& v) { + ColorTransformationMatrix color_transformation; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 4; ++j) + color_transformation[i][j] = v[i][j]; + } + return color_transformation; +} + +} // namespace + +ColorTransformationMatrix CalculateColorTransformationMatrix( + const rtc::scoped_refptr& reference_frame, + const rtc::scoped_refptr& test_frame) { + IncrementalLinearLeastSquares incremental_lls; + incremental_lls.AddObservations(FlattenYuvData(test_frame), + FlattenYuvData(reference_frame)); + return VectorToColorMatrix(incremental_lls.GetBestSolution()); +} + +ColorTransformationMatrix CalculateColorTransformationMatrix( + const rtc::scoped_refptr