From 1927dfafabc06a7eced67d1d3fd60f2c0c88c152 Mon Sep 17 00:00:00 2001 From: Magnus Jedvert Date: Tue, 11 Sep 2018 12:56:06 +0200 Subject: [PATCH] Add tool for aligning color space of video files MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This class adds logic for aligning color space of a test video compared to a reference video. If there is a color space mismatch, it typically does not have much impact on human perception, but it has a big impact on PSNR and SSIM calculations. For example, aligning a test run with VP8 improves PSNR and SSIM from: Average PSNR: 29.142818, average SSIM: 0.946026 to: Average PSNR: 38.146229, average SSIM: 0.965388. The optiomal color transformation between the two videos were: 0.86 0.01 0.00 14.37 0.00 0.88 0.00 15.32 0.00 0.00 0.88 15.74 which is converting YUV full range to YUV limited range. There is already a CL out for fixing this discrepancy here: https://webrtc-review.googlesource.com/c/src/+/94543 After that, hopefully there is no color space mismatch when saving the raw YUV values. It's good that the video quality tool is color space agnostic anyway, and can compensate for differences when the test video is obtained by e.g. filming a physical device screen. Also, the linear least square logic will be used for compensating geometric distorisions in a follow-up CL. Bug: webrtc:9642 Change-Id: I499713960a0544d8e45c5d09886e68ec829b28a7 Reviewed-on: https://webrtc-review.googlesource.com/c/95950 Reviewed-by: Sami Kalliomäki Reviewed-by: Patrik Höglund Commit-Queue: Magnus Jedvert Cr-Commit-Position: refs/heads/master@{#25193} --- rtc_tools/BUILD.gn | 8 + .../frame_analyzer/linear_least_squares.cc | 200 +++++++++++++++ .../frame_analyzer/linear_least_squares.h | 53 ++++ .../linear_least_squares_unittest.cc | 91 +++++++ .../frame_analyzer/video_color_aligner.cc | 239 ++++++++++++++++++ .../frame_analyzer/video_color_aligner.h | 49 ++++ .../video_color_aligner_unittest.cc | 174 +++++++++++++ 7 files changed, 814 insertions(+) create mode 100644 rtc_tools/frame_analyzer/linear_least_squares.cc create mode 100644 rtc_tools/frame_analyzer/linear_least_squares.h create mode 100644 rtc_tools/frame_analyzer/linear_least_squares_unittest.cc create mode 100644 rtc_tools/frame_analyzer/video_color_aligner.cc create mode 100644 rtc_tools/frame_analyzer/video_color_aligner.h create mode 100644 rtc_tools/frame_analyzer/video_color_aligner_unittest.cc 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