From d41a3d6d16bccb057aa78ffc68a9a65219c9ce7e Mon Sep 17 00:00:00 2001 From: Bohdan Tyshchenko Date: Tue, 11 Aug 2020 15:48:48 -0700 Subject: [PATCH] First version of guetzli sandbox --- oss-internship-2020/guetzli/.bazelrc | 2 + oss-internship-2020/guetzli/BUILD.bazel | 31 +- oss-internship-2020/guetzli/WORKSPACE | 31 +- .../guetzli/external/guetzli.BUILD | 18 + .../guetzli/guetzli_entry_points.cc | 237 +- .../guetzli/guetzli_entry_points.h | 31 +- .../guetzli/guetzli_sandboxed.cc | 44 +- .../guetzli/guetzli_transaction.cc | 180 +- .../guetzli/guetzli_transaction.h | 190 +- .../guetzli/third_party/butteraugli/BUILD | 24 + .../guetzli/third_party/butteraugli/LICENSE | 201 ++ .../guetzli/third_party/butteraugli/README.md | 68 + .../guetzli/third_party/butteraugli/WORKSPACE | 25 + .../butteraugli/butteraugli/Makefile | 7 + .../butteraugli/butteraugli/butteraugli.cc | 1994 +++++++++++++++++ .../butteraugli/butteraugli/butteraugli.h | 619 +++++ .../butteraugli/butteraugli_main.cc | 457 ++++ .../third_party/butteraugli/jpeg.BUILD | 89 + .../guetzli/third_party/butteraugli/png.BUILD | 33 + .../third_party/butteraugli/zlib.BUILD | 36 + 20 files changed, 3835 insertions(+), 482 deletions(-) create mode 100644 oss-internship-2020/guetzli/.bazelrc create mode 100644 oss-internship-2020/guetzli/external/guetzli.BUILD create mode 100755 oss-internship-2020/guetzli/third_party/butteraugli/BUILD create mode 100755 oss-internship-2020/guetzli/third_party/butteraugli/LICENSE create mode 100755 oss-internship-2020/guetzli/third_party/butteraugli/README.md create mode 100755 oss-internship-2020/guetzli/third_party/butteraugli/WORKSPACE create mode 100755 oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/Makefile create mode 100755 oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli.cc create mode 100755 oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli.h create mode 100755 oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli_main.cc create mode 100755 oss-internship-2020/guetzli/third_party/butteraugli/jpeg.BUILD create mode 100755 oss-internship-2020/guetzli/third_party/butteraugli/png.BUILD create mode 100755 oss-internship-2020/guetzli/third_party/butteraugli/zlib.BUILD diff --git a/oss-internship-2020/guetzli/.bazelrc b/oss-internship-2020/guetzli/.bazelrc new file mode 100644 index 0000000..a68e070 --- /dev/null +++ b/oss-internship-2020/guetzli/.bazelrc @@ -0,0 +1,2 @@ +# Build in C++17 mode without a custom CROSSTOOL +build --cxxopt=-std=c++17 diff --git a/oss-internship-2020/guetzli/BUILD.bazel b/oss-internship-2020/guetzli/BUILD.bazel index 345f879..6e3332a 100644 --- a/oss-internship-2020/guetzli/BUILD.bazel +++ b/oss-internship-2020/guetzli/BUILD.bazel @@ -1,9 +1,3 @@ -load("@rules_cc//cc:defs.bzl", "cc_proto_library") -load("@rules_proto//proto:defs.bzl", "proto_library") -load( - "@com_google_sandboxed_api//sandboxed_api/bazel:proto.bzl", - "sapi_proto_library", -) load( "@com_google_sandboxed_api//sandboxed_api/bazel:sapi.bzl", "sapi_library", @@ -17,23 +11,17 @@ cc_library( "@guetzli//:guetzli_lib", "@com_google_sandboxed_api//sandboxed_api:lenval_core", "@com_google_sandboxed_api//sandboxed_api:vars", - #"@com_google_sandboxed_api//sandboxed_api/sandbox2/util:temp_file", visibility error "@png_archive//:png" ], - visibility = ["//visibility:public"] ) sapi_library( name = "guetzli_sapi", - #srcs = ["guetzli_transaction.cc"], // Error when try to place definitions insde .cc file + srcs = ["guetzli_transaction.cc"], hdrs = ["guetzli_sandbox.h", "guetzli_transaction.h"], functions = [ - "ProcessJPEGString", - "ProcessRGBData", - "ButteraugliScoreQuality", - "ReadPng", - "ReadJpegData", - "ReadDataFromFd", + "ProcessJpeg", + "ProcessRgb", "WriteDataToFd" ], input_files = ["guetzli_entry_points.h"], @@ -43,23 +31,10 @@ sapi_library( namespace = "guetzli::sandbox" ) -# cc_library( -# name = "guetzli_sapi_transaction", -# #srcs = ["guetzli_transaction.cc"], -# hdrs = ["guetzli_transaction.h"], -# deps = [ -# ":guetzli_sapi" -# ], -# visibility = ["//visibility:public"] -# ) - cc_binary( name="guetzli_sandboxed", srcs=["guetzli_sandboxed.cc"], - includes = ["."], - visibility= [ "//visibility:public" ], deps = [ - #":guetzli_sapi_transaction" ":guetzli_sapi" ] ) \ No newline at end of file diff --git a/oss-internship-2020/guetzli/WORKSPACE b/oss-internship-2020/guetzli/WORKSPACE index 17887de..f7ec64c 100644 --- a/oss-internship-2020/guetzli/WORKSPACE +++ b/oss-internship-2020/guetzli/WORKSPACE @@ -15,6 +15,7 @@ workspace(name = "guetzli_sandboxed") load("@bazel_tools//tools/build_defs/repo:git.bzl", "git_repository") +load("@bazel_tools//tools/build_defs/repo:http.bzl", "http_archive") load("@bazel_tools//tools/build_defs/repo:utils.bzl", "maybe") # Include the Sandboxed API dependency if it does not already exist in this @@ -46,22 +47,19 @@ load("@com_google_protobuf//:protobuf_deps.bzl", "protobuf_deps") protobuf_deps() -maybe( - git_repository, +local_repository( + name = "butteraugli", + path = "third_party/butteraugli/" +) + +http_archive( name = "guetzli", - remote = "https://github.com/google/guetzli.git", - branch = "master" + build_file = "guetzli.BUILD", + sha256 = "39632357e49db83d9560bf0de560ad833352f36d23b109b0e995b01a37bddb57", + strip_prefix = "guetzli-master", + url = "https://github.com/google/guetzli/archive/master.zip" ) -maybe( - git_repository, - name = "googletest", - remote = "https://github.com/google/googletest", - tag = "release-1.10.0" -) - -load("@bazel_tools//tools/build_defs/repo:http.bzl", "http_archive") - http_archive( name = "png_archive", build_file = "png.BUILD", @@ -76,4 +74,11 @@ http_archive( sha256 = "8d7e9f698ce48787b6e1c67e6bff79e487303e66077e25cb9784ac8835978017", strip_prefix = "zlib-1.2.10", url = "http://zlib.net/fossils/zlib-1.2.10.tar.gz", +) + +maybe( + git_repository, + name = "googletest", + remote = "https://github.com/google/googletest", + tag = "release-1.10.0" ) \ No newline at end of file diff --git a/oss-internship-2020/guetzli/external/guetzli.BUILD b/oss-internship-2020/guetzli/external/guetzli.BUILD new file mode 100644 index 0000000..dec29a1 --- /dev/null +++ b/oss-internship-2020/guetzli/external/guetzli.BUILD @@ -0,0 +1,18 @@ +package(default_visibility = ["//visibility:public"]) + +cc_library( + name = "guetzli_lib", + srcs = glob( + [ + "guetzli/*.h", + "guetzli/*.cc", + "guetzli/*.inc", + ], + exclude = ["guetzli/guetzli.cc"], + ), + copts = [ "-Wno-sign-compare" ], + visibility= [ "//visibility:public" ], + deps = [ + "@butteraugli//:butteraugli_lib", + ], +) \ No newline at end of file diff --git a/oss-internship-2020/guetzli/guetzli_entry_points.cc b/oss-internship-2020/guetzli/guetzli_entry_points.cc index 446150c..a815673 100644 --- a/oss-internship-2020/guetzli/guetzli_entry_points.cc +++ b/oss-internship-2020/guetzli/guetzli_entry_points.cc @@ -3,18 +3,25 @@ #include "guetzli_entry_points.h" #include "png.h" #include "sandboxed_api/sandbox2/util/fileops.h" +#include "sandboxed_api/util/statusor.h" #include #include +#include #include #include #include namespace { -inline uint8_t BlendOnBlack(const uint8_t val, const uint8_t alpha) { - return (static_cast(val) * static_cast(alpha) + 128) / 255; -} +constexpr int kBytesPerPixel = 350; +constexpr int kLowestMemusageMB = 100; // in MB + +struct GuetzliInitData { + std::string in_data; + guetzli::Params params; + guetzli::ProcessStats stats; +}; template void CopyMemoryToLenVal(const T* data, size_t size, @@ -26,68 +33,63 @@ void CopyMemoryToLenVal(const T* data, size_t size, out_data->data = new_out; } -} // namespace - -extern "C" bool ProcessJPEGString(const guetzli::Params* params, - int verbose, - sapi::LenValStruct* in_data, - sapi::LenValStruct* out_data) -{ - std::string in_data_temp(static_cast(in_data->data), - in_data->size); - - guetzli::ProcessStats stats; - if (verbose > 0) { - stats.debug_output_file = stderr; - } - - std::string temp_out = ""; - auto result = guetzli::Process(*params, &stats, in_data_temp, &temp_out); - - if (result) { - CopyMemoryToLenVal(temp_out.data(), temp_out.size(), out_data); - } - - return result; -} - -extern "C" bool ProcessRGBData(const guetzli::Params* params, - int verbose, - sapi::LenValStruct* rgb, - int w, int h, - sapi::LenValStruct* out_data) -{ - std::vector in_data_temp; - in_data_temp.reserve(rgb->size); - - auto* rgb_data = static_cast(rgb->data); - std::copy(rgb_data, rgb_data + rgb->size, std::back_inserter(in_data_temp)); - - guetzli::ProcessStats stats; - if (verbose > 0) { - stats.debug_output_file = stderr; - } - - std::string temp_out = ""; - auto result = - guetzli::Process(*params, &stats, in_data_temp, w, h, &temp_out); - - //TODO: Move shared part of the code to another function - if (result) { - CopyMemoryToLenVal(temp_out.data(), temp_out.size(), out_data); - } - - return result; -} - -extern "C" bool ReadPng(sapi::LenValStruct* in_data, - int* xsize, int* ysize, - sapi::LenValStruct* rgb_out) -{ - std::string data(static_cast(in_data->data), in_data->size); - std::vector rgb; +sapi::StatusOr ReadFromFd(int fd) { + struct stat file_data; + auto status = fstat(fd, &file_data); - png_structp png_ptr = + if (status < 0) { + return absl::FailedPreconditionError( + "Error reading input from fd" + ); + } + + auto fsize = file_data.st_size; + + std::unique_ptr buf(new char[fsize]); + status = read(fd, buf.get(), fsize); + + if (status < 0) { + lseek(fd, 0, SEEK_SET); + return absl::FailedPreconditionError( + "Error reading input from fd" + ); + } + + return std::string(buf.get(), fsize); +} + +sapi::StatusOr PrepareDataForProcessing( + const ProcessingParams* processing_params) { + auto input_status = ReadFromFd(processing_params->remote_fd); + + if (!input_status.ok()) { + return input_status.status(); + } + + guetzli::Params guetzli_params; + guetzli_params.butteraugli_target = static_cast( + guetzli::ButteraugliScoreForQuality(processing_params->quality)); + + guetzli::ProcessStats stats; + + if (processing_params->verbose) { + stats.debug_output_file = stderr; + } + + return GuetzliInitData{ + std::move(input_status.value()), + guetzli_params, + stats + }; +} + +inline uint8_t BlendOnBlack(const uint8_t val, const uint8_t alpha) { + return (static_cast(val) * static_cast(alpha) + 128) / 255; +} + +bool ReadPNG(const std::string& data, int* xsize, int* ysize, + std::vector* rgb) { + png_structp png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING, nullptr, nullptr, nullptr); if (!png_ptr) { return false; @@ -129,7 +131,7 @@ extern "C" bool ReadPng(sapi::LenValStruct* in_data, *xsize = png_get_image_width(png_ptr, info_ptr); *ysize = png_get_image_height(png_ptr, info_ptr); - rgb.resize(3 * (*xsize) * (*ysize)); + rgb->resize(3 * (*xsize) * (*ysize)); const int components = png_get_channels(png_ptr, info_ptr); switch (components) { @@ -137,7 +139,7 @@ extern "C" bool ReadPng(sapi::LenValStruct* in_data, // GRAYSCALE for (int y = 0; y < *ysize; ++y) { const uint8_t* row_in = row_pointers[y]; - uint8_t* row_out = &(rgb)[3 * y * (*xsize)]; + uint8_t* row_out = &(*rgb)[3 * y * (*xsize)]; for (int x = 0; x < *xsize; ++x) { const uint8_t gray = row_in[x]; row_out[3 * x + 0] = gray; @@ -151,7 +153,7 @@ extern "C" bool ReadPng(sapi::LenValStruct* in_data, // GRAYSCALE + ALPHA for (int y = 0; y < *ysize; ++y) { const uint8_t* row_in = row_pointers[y]; - uint8_t* row_out = &(rgb)[3 * y * (*xsize)]; + uint8_t* row_out = &(*rgb)[3 * y * (*xsize)]; for (int x = 0; x < *xsize; ++x) { const uint8_t gray = BlendOnBlack(row_in[2 * x], row_in[2 * x + 1]); row_out[3 * x + 0] = gray; @@ -165,7 +167,7 @@ extern "C" bool ReadPng(sapi::LenValStruct* in_data, // RGB for (int y = 0; y < *ysize; ++y) { const uint8_t* row_in = row_pointers[y]; - uint8_t* row_out = &(rgb)[3 * y * (*xsize)]; + uint8_t* row_out = &(*rgb)[3 * y * (*xsize)]; memcpy(row_out, row_in, 3 * (*xsize)); } break; @@ -174,7 +176,7 @@ extern "C" bool ReadPng(sapi::LenValStruct* in_data, // RGBA for (int y = 0; y < *ysize; ++y) { const uint8_t* row_in = row_pointers[y]; - uint8_t* row_out = &(rgb)[3 * y * (*xsize)]; + uint8_t* row_out = &(*rgb)[3 * y * (*xsize)]; for (int x = 0; x < *xsize; ++x) { const uint8_t alpha = row_in[4 * x + 3]; row_out[3 * x + 0] = BlendOnBlack(row_in[4 * x + 0], alpha); @@ -189,53 +191,84 @@ extern "C" bool ReadPng(sapi::LenValStruct* in_data, return false; } png_destroy_read_struct(&png_ptr, &info_ptr, nullptr); - - CopyMemoryToLenVal(rgb.data(), rgb.size(), rgb_out); - return true; } -extern "C" bool ReadJpegData(sapi::LenValStruct* in_data, - int mode, - int* xsize, int* ysize) -{ - std::string data(static_cast(in_data->data), in_data->size); - guetzli::JPEGData jpg; - - auto result = guetzli::ReadJpeg(data, - static_cast(mode), &jpg); - - if (result) { - *xsize = jpg.width; - *ysize = jpg.height; - } - - return result; +bool CheckMemoryLimitExceeded(int memlimit_mb, int xsize, int ysize) { + double pixels = static_cast(xsize) * ysize; + return memlimit_mb != -1 + && (pixels * kBytesPerPixel / (1 << 20) > memlimit_mb + || memlimit_mb < kLowestMemusageMB); } -extern "C" double ButteraugliScoreQuality(double quality) { - return guetzli::ButteraugliScoreForQuality(quality); +} // namespace + +extern "C" bool ProcessJpeg(const ProcessingParams* processing_params, + sapi::LenValStruct* output) { + auto processing_data_status = PrepareDataForProcessing(processing_params); + + if (!processing_data_status.status().ok()) { + fprintf(stderr, "%s\n", processing_data_status.status().ToString().c_str()); + return false; + } + + guetzli::JPEGData jpg_header; + if (!guetzli::ReadJpeg(processing_data_status.value().in_data, + guetzli::JPEG_READ_HEADER, &jpg_header)) { + fprintf(stderr, "Error reading JPG data from input file\n"); + return false; + } + + if (CheckMemoryLimitExceeded(processing_params->memlimit_mb, + jpg_header.width, jpg_header.height)) { + fprintf(stderr, "Memory limit would be exceeded.\n"); + return false; + } + + std::string out_data; + if (!guetzli::Process(processing_data_status.value().params, + &processing_data_status.value().stats, + processing_data_status.value().in_data, + &out_data)) { + fprintf(stderr, "Guezli processing failed\n"); + return false; + } + + CopyMemoryToLenVal(out_data.data(), out_data.size(), output); + return true; } -extern "C" bool ReadDataFromFd(int fd, sapi::LenValStruct* out_data) { - struct stat file_data; - auto status = fstat(fd, &file_data); - - if (status < 0) { +extern "C" bool ProcessRgb(const ProcessingParams* processing_params, + sapi::LenValStruct* output) { + auto processing_data_status = PrepareDataForProcessing(processing_params); + + if (!processing_data_status.status().ok()) { + fprintf(stderr, "%s\n", processing_data_status.status().ToString().c_str()); return false; } - - auto fsize = file_data.st_size; - std::unique_ptr buf(new char[fsize]); - status = read(fd, buf.get(), fsize); + int xsize, ysize; + std::vector rgb; - if (status < 0) { + if (!ReadPNG(processing_data_status.value().in_data, &xsize, &ysize, &rgb)) { + fprintf(stderr, "Error reading PNG data from input file\n"); return false; } - - CopyMemoryToLenVal(buf.get(), fsize, out_data); + if (CheckMemoryLimitExceeded(processing_params->memlimit_mb, xsize, ysize)) { + fprintf(stderr, "Memory limit would be exceeded.\n"); + return false; + } + + std::string out_data; + if (!guetzli::Process(processing_data_status.value().params, + &processing_data_status.value().stats, + rgb, xsize, ysize, &out_data)) { + fprintf(stderr, "Guetzli processing failed\n"); + return false; + } + + CopyMemoryToLenVal(out_data.data(), out_data.size(), output); return true; } diff --git a/oss-internship-2020/guetzli/guetzli_entry_points.h b/oss-internship-2020/guetzli/guetzli_entry_points.h index 6fd0f11..90b0d97 100644 --- a/oss-internship-2020/guetzli/guetzli_entry_points.h +++ b/oss-internship-2020/guetzli/guetzli_entry_points.h @@ -4,26 +4,15 @@ #include "sandboxed_api/lenval_core.h" #include "sandboxed_api/vars.h" -extern "C" bool ProcessJPEGString(const guetzli::Params* params, - int verbose, - sapi::LenValStruct* in_data, - sapi::LenValStruct* out_data); - -extern "C" bool ProcessRGBData(const guetzli::Params* params, - int verbose, - sapi::LenValStruct* rgb, - int w, int h, - sapi::LenValStruct* out_data); - -extern "C" bool ReadPng(sapi::LenValStruct* in_data, - int* xsize, int* ysize, - sapi::LenValStruct* rgb_out); - -extern "C" bool ReadJpegData(sapi::LenValStruct* in_data, - int mode, int* xsize, int* ysize); - -extern "C" double ButteraugliScoreQuality(double quality); - -extern "C" bool ReadDataFromFd(int fd, sapi::LenValStruct* out_data); +struct ProcessingParams { + int remote_fd = -1; + int verbose = 0; + int quality = 0; + int memlimit_mb = 0; +}; +extern "C" bool ProcessJpeg(const ProcessingParams* processing_params, + sapi::LenValStruct* output); +extern "C" bool ProcessRgb(const ProcessingParams* processing_params, + sapi::LenValStruct* output); extern "C" bool WriteDataToFd(int fd, sapi::LenValStruct* data); \ No newline at end of file diff --git a/oss-internship-2020/guetzli/guetzli_sandboxed.cc b/oss-internship-2020/guetzli/guetzli_sandboxed.cc index a62f249..3393b8f 100644 --- a/oss-internship-2020/guetzli/guetzli_sandboxed.cc +++ b/oss-internship-2020/guetzli/guetzli_sandboxed.cc @@ -14,17 +14,6 @@ namespace { constexpr int kDefaultJPEGQuality = 95; constexpr int kDefaultMemlimitMB = 6000; // in MB -//constexpr absl::string_view kMktempSuffix = "XXXXXX"; - -// sapi::StatusOr> CreateNamedTempFile( -// absl::string_view prefix) { -// std::string name_template = absl::StrCat(prefix, kMktempSuffix); -// int fd = mkstemp(&name_template[0]); -// if (fd < 0) { -// return absl::UnknownError("Error creating temp file"); -// } -// return std::pair{std::move(name_template), fd}; -// } void TerminateHandler() { fprintf(stderr, "Unhandled exception. Most likely insufficient memory available.\n" @@ -96,21 +85,6 @@ int main(int argc, const char** argv) { return 1; } - // auto out_temp_file = CreateNamedTempFile("/tmp/" + std::string(argv[opt_idx + 1])); - // if (!out_temp_file.ok()) { - // fprintf(stderr, "Can't create temporary output file: %s\n", - // argv[opt_idx + 1]); - // return 1; - // } - // sandbox2::file_util::fileops::FDCloser out_fd_closer( - // out_temp_file.value().second); - - // if (unlink(out_temp_file.value().first.c_str()) < 0) { - // fprintf(stderr, "Error unlinking temp out file: %s\n", - // out_temp_file.value().first.c_str()); - // return 1; - // } - sandbox2::file_util::fileops::FDCloser out_fd_closer( open(".", O_TMPFILE | O_RDWR, S_IRUSR | S_IWUSR)); @@ -119,14 +93,6 @@ int main(int argc, const char** argv) { return 1; } - // sandbox2::file_util::fileops::FDCloser out_fd_closer(open(argv[opt_idx + 1], - // O_RDWR | O_CREAT, S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP)); - - // if (out_fd_closer.get() < 0) { - // fprintf(stderr, "Can't open output file: %s\n", argv[opt_idx + 1]); - // return 1; - // } - guetzli::sandbox::TransactionParams params = { in_fd_closer.get(), out_fd_closer.get(), @@ -143,13 +109,19 @@ int main(int argc, const char** argv) { if (remove(argv[opt_idx + 1]) < 0) { fprintf(stderr, "Error deleting existing output file: %s\n", argv[opt_idx + 1]); + return 1; } } std::stringstream path; path << "/proc/self/fd/" << out_fd_closer.get(); - linkat(AT_FDCWD, path.str().c_str(), AT_FDCWD, argv[opt_idx + 1], - AT_SYMLINK_FOLLOW); + + if (linkat(AT_FDCWD, path.str().c_str(), AT_FDCWD, argv[opt_idx + 1], + AT_SYMLINK_FOLLOW) < 0) { + fprintf(stderr, "Error linking %s\n", + argv[opt_idx + 1]); + return 1; + } } else { fprintf(stderr, "%s\n", result.ToString().c_str()); // Use cerr instead ? diff --git a/oss-internship-2020/guetzli/guetzli_transaction.cc b/oss-internship-2020/guetzli/guetzli_transaction.cc index 9403783..0f7e1bd 100644 --- a/oss-internship-2020/guetzli/guetzli_transaction.cc +++ b/oss-internship-2020/guetzli/guetzli_transaction.cc @@ -1,11 +1,38 @@ #include "guetzli_transaction.h" #include +#include namespace guetzli { namespace sandbox { absl::Status GuetzliTransaction::Init() { + // Close remote fd if transaction is repeated + if (in_fd_.GetRemoteFd() != -1) { + SAPI_RETURN_IF_ERROR(in_fd_.CloseRemoteFd(sandbox()->GetRpcChannel())); + } + if (out_fd_.GetRemoteFd() != -1) { + SAPI_RETURN_IF_ERROR(out_fd_.CloseRemoteFd(sandbox()->GetRpcChannel())); + } + + // Reposition back to the beginning of file + if (lseek(in_fd_.GetValue(), 0, SEEK_CUR) != 0) { + if (lseek(in_fd_.GetValue(), 0, SEEK_SET) != 0) { + return absl::FailedPreconditionError( + "Error returnig cursor to the beginning" + ); + } + } + + // Choosing between jpg and png modes + sapi::StatusOr image_type = GetImageTypeFromFd(in_fd_.GetValue()); + + if (!image_type.ok()) { + return image_type.status(); + } + + image_type_ = image_type.value(); + SAPI_RETURN_IF_ERROR(sandbox()->TransferToSandboxee(&in_fd_)); SAPI_RETURN_IF_ERROR(sandbox()->TransferToSandboxee(&out_fd_)); @@ -18,125 +45,36 @@ absl::Status GuetzliTransaction::Init() { "Error receiving remote FD: remote output fd is set to -1"); } + in_fd_.OwnLocalFd(false); // FDCloser will close local fd + out_fd_.OwnLocalFd(false); // FDCloser will close local fd + return absl::OkStatus(); } - absl::Status GuetzliTransaction::ProcessPng(GuetzliAPi* api, - sapi::v::Struct* params, - sapi::v::LenVal* input, - sapi::v::LenVal* output) const { - sapi::v::Int xsize; - sapi::v::Int ysize; - sapi::v::LenVal rgb_in(0); - - auto read_result = api->ReadPng(input->PtrBefore(), xsize.PtrBoth(), - ysize.PtrBoth(), rgb_in.PtrBoth()); - - if (!read_result.value_or(false)) { - return absl::FailedPreconditionError( - "Error reading PNG data from input file" - ); - } - - double pixels = static_cast(xsize.GetValue()) * ysize.GetValue(); - if (params_.memlimit_mb != -1 - && (pixels * kBytesPerPixel / (1 << 20) > params_.memlimit_mb - || params_.memlimit_mb < kLowestMemusageMB)) { - return absl::FailedPreconditionError( - "Memory limit would be exceeded" - ); - } - - auto result = api->ProcessRGBData(params->PtrBefore(), params_.verbose, - rgb_in.PtrBefore(), xsize.GetValue(), - ysize.GetValue(), output->PtrBoth()); - if (!result.value_or(false)) { - return absl::FailedPreconditionError( - "Guetzli processing failed" - ); - } - - return absl::OkStatus(); - } - - absl::Status GuetzliTransaction::ProcessJpeg(GuetzliApi* api, - sapi::v::Struct* params, - sapi::v::LenVal* input, - sapi::v::LenVal* output) const { - ::sapi::v::Int xsize; - ::sapi::v::Int ysize; - auto read_result = api->ReadJpegData(input->PtrBefore(), 0, xsize.PtrBoth(), - ysize.PtrBoth()); - - if (!read_result.value_or(false)) { - return absl::FailedPreconditionError( - "Error reading JPG data from input file" - ); - } - - double pixels = static_cast(xsize.GetValue()) * ysize.GetValue(); - if (params_.memlimit_mb != -1 - && (pixels * kBytesPerPixel / (1 << 20) > params_.memlimit_mb - || params_.memlimit_mb < kLowestMemusageMB)) { - return absl::FailedPreconditionError( - "Memory limit would be exceeded" - ); - } - - auto result = api->ProcessJPEGString(params->PtrBefore(), params_.verbose, - input->PtrBefore(), output->PtrBoth()); - - if (!result.value_or(false)) { - return absl::FailedPreconditionError( - "Guetzli processing failed" - ); - } - - return absl::OkStatus(); - } - absl::Status GuetzliTransaction::Main() { GuetzliApi api(sandbox()); - - sapi::v::LenVal input(0); sapi::v::LenVal output(0); - sapi::v::Struct params; - - auto read_result = api.ReadDataFromFd(in_fd_.GetRemoteFd(), input.PtrBoth()); - if (!read_result.value_or(false)) { - return absl::FailedPreconditionError( - "Error reading data inside sandbox" - ); - } - - auto score_quality_result = api.ButteraugliScoreQuality(params_.quality); - - if (!score_quality_result.ok()) { - return absl::FailedPreconditionError( - "Error calculating butteraugli score" - ); - } - - params.mutable_data()->butteraugli_target = score_quality_result.value(); - - static const unsigned char kPNGMagicBytes[] = { - 0x89, 'P', 'N', 'G', '\r', '\n', 0x1a, '\n', + sapi::v::Struct processing_params; + *processing_params.mutable_data() = {in_fd_.GetRemoteFd(), + params_.verbose, + params_.quality, + params_.memlimit_mb }; + + auto result_status = image_type_ == ImageType::JPEG ? + api.ProcessJpeg(processing_params.PtrBefore(), output.PtrBoth()) : + api.ProcessRgb(processing_params.PtrBefore(), output.PtrBoth()); - if (input.GetDataSize() >= 8 && - memcmp(input.GetData(), kPNGMagicBytes, sizeof(kPNGMagicBytes)) == 0) { - auto process_status = ProcessPng(&api, ¶ms, &input, &output); - - if (!process_status.ok()) { - return process_status; - } - } else { - auto process_status = ProcessJpeg(&api, ¶ms, &input, &output); - - if (!process_status.ok()) { - return process_status; - } + if (!result_status.value_or(false)) { + std::stringstream error_stream; + error_stream << "Error processing " + << (image_type_ == ImageType::JPEG ? "jpeg" : "rgb") << " data" + << std::endl; + + return absl::FailedPreconditionError( + error_stream.str() + ); } auto write_result = api.WriteDataToFd(out_fd_.GetRemoteFd(), @@ -156,5 +94,27 @@ time_t GuetzliTransaction::CalculateTimeLimitFromImageSize( return (pixels / kMpixPixels + 5) * 60; } +sapi::StatusOr GuetzliTransaction::GetImageTypeFromFd(int fd) const { + static const unsigned char kPNGMagicBytes[] = { + 0x89, 'P', 'N', 'G', '\r', '\n', 0x1a, '\n', + }; + char read_buf[8]; + + if (read(fd, read_buf, 8) != 8) { + return absl::FailedPreconditionError( + "Error determining type of the input file" + ); + } + + if (lseek(fd, 0, SEEK_SET) != 0) { + return absl::FailedPreconditionError( + "Error returnig cursor to the beginning" + ); + } + + return memcmp(read_buf, kPNGMagicBytes, sizeof(kPNGMagicBytes)) == 0 ? + ImageType::PNG : ImageType::JPEG; +} + } // namespace sandbox } // namespace guetzli \ No newline at end of file diff --git a/oss-internship-2020/guetzli/guetzli_transaction.h b/oss-internship-2020/guetzli/guetzli_transaction.h index 5d6e83a..8dbd7ec 100644 --- a/oss-internship-2020/guetzli/guetzli_transaction.h +++ b/oss-internship-2020/guetzli/guetzli_transaction.h @@ -10,22 +10,24 @@ namespace guetzli { namespace sandbox { -constexpr int kDefaultTransactionRetryCount = 1; +constexpr int kDefaultTransactionRetryCount = 0; constexpr uint64_t kMpixPixels = 1'000'000; -constexpr int kBytesPerPixel = 350; -constexpr int kLowestMemusageMB = 100; // in MB - -struct TransactionParams { - int in_fd; - int out_fd; - int verbose; - int quality; - int memlimit_mb; +enum class ImageType { + JPEG, + PNG }; -//Add optional time limit/retry count as a constructors arguments -//Use differenet status errors +struct TransactionParams { + int in_fd = -1; + int out_fd = -1; + int verbose = 0; + int quality = 0; + int memlimit_mb = 0; +}; + +// Instance of this transaction shouldn't be reused +// Create a new one for each processing operation class GuetzliTransaction : public sapi::Transaction { public: GuetzliTransaction(TransactionParams&& params) @@ -34,7 +36,9 @@ class GuetzliTransaction : public sapi::Transaction { , in_fd_(params_.in_fd) , out_fd_(params_.out_fd) { + //TODO: Add retry count as a parameter sapi::Transaction::set_retry_count(kDefaultTransactionRetryCount); + //TODO: Try to use sandbox().set_wall_limit instead of infinite time limit sapi::Transaction::SetTimeLimit(0); } @@ -42,15 +46,7 @@ class GuetzliTransaction : public sapi::Transaction { absl::Status Init() override; absl::Status Main() final; - absl::Status ProcessPng(GuetzliApi* api, - sapi::v::Struct* params, - sapi::v::LenVal* input, - sapi::v::LenVal* output) const; - - absl::Status ProcessJpeg(GuetzliApi* api, - sapi::v::Struct* params, - sapi::v::LenVal* input, - sapi::v::LenVal* output) const; + sapi::StatusOr GetImageTypeFromFd(int fd) const; // As guetzli takes roughly 1 minute of CPU per 1 MPix we need to calculate // approximate time for transaction to complete @@ -59,158 +55,8 @@ class GuetzliTransaction : public sapi::Transaction { const TransactionParams params_; sapi::v::Fd in_fd_; sapi::v::Fd out_fd_; + ImageType image_type_ = ImageType::JPEG; }; -absl::Status GuetzliTransaction::Init() { - SAPI_RETURN_IF_ERROR(sandbox()->TransferToSandboxee(&in_fd_)); - SAPI_RETURN_IF_ERROR(sandbox()->TransferToSandboxee(&out_fd_)); - - if (in_fd_.GetRemoteFd() < 0) { - return absl::FailedPreconditionError( - "Error receiving remote FD: remote input fd is set to -1"); - } - if (out_fd_.GetRemoteFd() < 0) { - return absl::FailedPreconditionError( - "Error receiving remote FD: remote output fd is set to -1"); - } - - return absl::OkStatus(); -} - - absl::Status GuetzliTransaction::ProcessPng(GuetzliApi* api, - sapi::v::Struct* params, - sapi::v::LenVal* input, - sapi::v::LenVal* output) const { - sapi::v::Int xsize; - sapi::v::Int ysize; - sapi::v::LenVal rgb_in(0); - - auto read_result = api->ReadPng(input->PtrBefore(), xsize.PtrBoth(), - ysize.PtrBoth(), rgb_in.PtrBoth()); - - if (!read_result.value_or(false)) { - return absl::FailedPreconditionError( - "Error reading PNG data from input file" - ); - } - - double pixels = static_cast(xsize.GetValue()) * ysize.GetValue(); - if (params_.memlimit_mb != -1 - && (pixels * kBytesPerPixel / (1 << 20) > params_.memlimit_mb - || params_.memlimit_mb < kLowestMemusageMB)) { - return absl::FailedPreconditionError( - "Memory limit would be exceeded" - ); - } - - auto result = api->ProcessRGBData(params->PtrBefore(), params_.verbose, - rgb_in.PtrBefore(), xsize.GetValue(), - ysize.GetValue(), output->PtrBoth()); - if (!result.value_or(false)) { - return absl::FailedPreconditionError( - "Guetzli processing failed" - ); - } - - return absl::OkStatus(); - } - - absl::Status GuetzliTransaction::ProcessJpeg(GuetzliApi* api, - sapi::v::Struct* params, - sapi::v::LenVal* input, - sapi::v::LenVal* output) const { - sapi::v::Int xsize; - sapi::v::Int ysize; - auto read_result = api->ReadJpegData(input->PtrBefore(), 0, xsize.PtrBoth(), - ysize.PtrBoth()); - - if (!read_result.value_or(false)) { - return absl::FailedPreconditionError( - "Error reading JPG data from input file" - ); - } - - double pixels = static_cast(xsize.GetValue()) * ysize.GetValue(); - if (params_.memlimit_mb != -1 - && (pixels * kBytesPerPixel / (1 << 20) > params_.memlimit_mb - || params_.memlimit_mb < kLowestMemusageMB)) { - return absl::FailedPreconditionError( - "Memory limit would be exceeded" - ); - } - - auto result = api->ProcessJPEGString(params->PtrBefore(), params_.verbose, - input->PtrBefore(), output->PtrBoth()); - - if (!result.value_or(false)) { - return absl::FailedPreconditionError( - "Guetzli processing failed" - ); - } - - return absl::OkStatus(); - } - -absl::Status GuetzliTransaction::Main() { - GuetzliApi api(sandbox()); - - sapi::v::LenVal input(0); - sapi::v::LenVal output(0); - sapi::v::Struct params; - - auto read_result = api.ReadDataFromFd(in_fd_.GetRemoteFd(), input.PtrBoth()); - - if (!read_result.value_or(false)) { - return absl::FailedPreconditionError( - "Error reading data inside sandbox" - ); - } - - auto score_quality_result = api.ButteraugliScoreQuality(params_.quality); - - if (!score_quality_result.ok()) { - return absl::FailedPreconditionError( - "Error calculating butteraugli score" - ); - } - - params.mutable_data()->butteraugli_target = score_quality_result.value(); - - static const unsigned char kPNGMagicBytes[] = { - 0x89, 'P', 'N', 'G', '\r', '\n', 0x1a, '\n', - }; - - if (input.GetDataSize() >= 8 && - memcmp(input.GetData(), kPNGMagicBytes, sizeof(kPNGMagicBytes)) == 0) { - auto process_status = ProcessPng(&api, ¶ms, &input, &output); - - if (!process_status.ok()) { - return process_status; - } - } else { - auto process_status = ProcessJpeg(&api, ¶ms, &input, &output); - - if (!process_status.ok()) { - return process_status; - } - } - - auto write_result = api.WriteDataToFd(out_fd_.GetRemoteFd(), - output.PtrBefore()); - - if (!write_result.value_or(false)) { - return absl::FailedPreconditionError( - "Error writing file inside sandbox" - ); - } - - return absl::OkStatus(); -} - -time_t GuetzliTransaction::CalculateTimeLimitFromImageSize( - uint64_t pixels) const { - return (pixels / kMpixPixels + 5) * 60; -} - } // namespace sandbox } // namespace guetzli diff --git a/oss-internship-2020/guetzli/third_party/butteraugli/BUILD b/oss-internship-2020/guetzli/third_party/butteraugli/BUILD new file mode 100755 index 0000000..f553c0b --- /dev/null +++ b/oss-internship-2020/guetzli/third_party/butteraugli/BUILD @@ -0,0 +1,24 @@ +cc_library( + name = "butteraugli_lib", + srcs = [ + "butteraugli/butteraugli.cc", + "butteraugli/butteraugli.h", + ], + hdrs = [ + "butteraugli/butteraugli.h", + ], + copts = ["-Wno-sign-compare"], + visibility = ["//visibility:public"], +) + +cc_binary( + name = "butteraugli", + srcs = ["butteraugli/butteraugli_main.cc"], + copts = ["-Wno-sign-compare"], + visibility = ["//visibility:public"], + deps = [ + ":butteraugli_lib", + "@jpeg_archive//:jpeg", + "@png_archive//:png", + ], +) \ No newline at end of file diff --git a/oss-internship-2020/guetzli/third_party/butteraugli/LICENSE b/oss-internship-2020/guetzli/third_party/butteraugli/LICENSE new file mode 100755 index 0000000..261eeb9 --- /dev/null +++ b/oss-internship-2020/guetzli/third_party/butteraugli/LICENSE @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/oss-internship-2020/guetzli/third_party/butteraugli/README.md b/oss-internship-2020/guetzli/third_party/butteraugli/README.md new file mode 100755 index 0000000..4623442 --- /dev/null +++ b/oss-internship-2020/guetzli/third_party/butteraugli/README.md @@ -0,0 +1,68 @@ +# butteraugli + +> A tool for measuring perceived differences between images + +## Introduction + +Butteraugli is a project that estimates the psychovisual similarity of two +images. It gives a score for the images that is reliable in the domain of barely +noticeable differences. Butteraugli not only gives a scalar score, but also +computes a spatial map of the level of differences. + +One of the main motivations for this project is the statistical differences in +location and density of different color receptors, particularly the low density +of blue cones in the fovea. Another motivation comes from more accurate modeling +of ganglion cells, particularly the frequency space inhibition. + +## Use + +Butteraugli can work as a quality metric for lossy image and video compression. +On our small test corpus butteraugli performs better than our implementations of +the reference methods, psnrhsv-m, ssim, and our yuv-color-space variant of ssim. +One possible use is to define the quality level setting used in a jpeg +compressor, or to compare two or more compression methods at the same level of +psychovisual differences. + +Butteraugli is intended to be a research tool more than a practical tool for +choosing compression formats. We don't know how well butteraugli performs with +major deformations -- we have mostly tuned it within a small range of quality, +roughly corresponding to jpeg qualities 90 to 95. + +## Interface + +Only a C++ interface is provided. The interface takes two images and outputs a +map together with a scalar value defining the difference. The scalar value can +be compared to two reference values that divide the value space into three +experience classes: 'great', 'acceptable' and 'not acceptable'. + +## Build instructions + +Install [Bazel](http://bazel.build) by following the +[instructions](https://www.bazel.build/docs/install.html). Run `bazel build -c opt +//:butteraugli` in the directory that contains this README file to build the +[command-line utility](#cmdline-tool). If you want to use Butteraugli as a +library, depend on the `//:butteraugli_lib` target. + +Alternatively, you can use the Makefile provided in the `butteraugli` directory, +after ensuring that [libpng](http://www.libpng.org/) and +[libjpeg](http://ijg.org/) are installed. On some systems you might need to also +install corresponding `-dev` packages. + +The code is portable and also compiles on Windows after defining +`_CRT_SECURE_NO_WARNINGS` in the project settings. + +## Command-line utility {#cmdline-tool} + +Butteraugli, apart from the library, comes bundled with a comparison tool. The +comparison tool supports PNG and JPG images as inputs. To compare images, run: + +``` +butteraugli image1.{png|jpg} image2.{png|jpg} +``` + +The tool can also produce a heatmap of differences between images. The heatmap +will be output as a PNM image. To produce one, run: + +``` +butteraugli image1.{png|jpg} image2.{png|jpg} heatmap.pnm +``` diff --git a/oss-internship-2020/guetzli/third_party/butteraugli/WORKSPACE b/oss-internship-2020/guetzli/third_party/butteraugli/WORKSPACE new file mode 100755 index 0000000..4d6ed65 --- /dev/null +++ b/oss-internship-2020/guetzli/third_party/butteraugli/WORKSPACE @@ -0,0 +1,25 @@ +workspace(name = "butteraugli") + +new_http_archive( + name = "png_archive", + url = "http://github.com/glennrp/libpng/archive/v1.2.57.zip", + sha256 = "a941dc09ca00148fe7aaf4ecdd6a67579c293678ed1e1cf633b5ffc02f4f8cf7", + strip_prefix = "libpng-1.2.57", + build_file = "png.BUILD", +) + +new_http_archive( + name = "zlib_archive", + url = "http://zlib.net/fossils/zlib-1.2.10.tar.gz", + sha256 = "8d7e9f698ce48787b6e1c67e6bff79e487303e66077e25cb9784ac8835978017", + strip_prefix = "zlib-1.2.10", + build_file = "zlib.BUILD", +) + +new_http_archive( + name = "jpeg_archive", + url = "http://www.ijg.org/files/jpegsrc.v9b.tar.gz", + sha256 = "240fd398da741669bf3c90366f58452ea59041cacc741a489b99f2f6a0bad052", + strip_prefix = "jpeg-9b", + build_file = "jpeg.BUILD", +) diff --git a/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/Makefile b/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/Makefile new file mode 100755 index 0000000..76b3a9b --- /dev/null +++ b/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/Makefile @@ -0,0 +1,7 @@ +LDLIBS += -lpng -ljpeg +CXXFLAGS += -std=c++11 -I.. +LINK.o = $(LINK.cc) + +all: butteraugli.o butteraugli_main.o butteraugli + +butteraugli: butteraugli.o butteraugli_main.o diff --git a/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli.cc b/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli.cc new file mode 100755 index 0000000..77c91cc --- /dev/null +++ b/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli.cc @@ -0,0 +1,1994 @@ +// Copyright 2016 Google Inc. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// +// Author: Jyrki Alakuijala (jyrki.alakuijala@gmail.com) +// +// The physical architecture of butteraugli is based on the following naming +// convention: +// * Opsin - dynamics of the photosensitive chemicals in the retina +// with their immediate electrical processing +// * Xyb - hybrid opponent/trichromatic color space +// x is roughly red-subtract-green. +// y is yellow. +// b is blue. +// Xyb values are computed from Opsin mixing, not directly from rgb. +// * Mask - for visual masking +// * Hf - color modeling for spatially high-frequency features +// * Lf - color modeling for spatially low-frequency features +// * Diffmap - to cluster and build an image of error between the images +// * Blur - to hold the smoothing code + +#include "butteraugli/butteraugli.h" + +#include +#include +#include +#include +#include + +#include +#include + + +// Restricted pointers speed up Convolution(); MSVC uses a different keyword. +#ifdef _MSC_VER +#define __restrict__ __restrict +#endif + +#ifndef PROFILER_ENABLED +#define PROFILER_ENABLED 0 +#endif +#if PROFILER_ENABLED +#else +#define PROFILER_FUNC +#define PROFILER_ZONE(name) +#endif + +namespace butteraugli { + +void *CacheAligned::Allocate(const size_t bytes) { + char *const allocated = static_cast(malloc(bytes + kCacheLineSize)); + if (allocated == nullptr) { + return nullptr; + } + const uintptr_t misalignment = + reinterpret_cast(allocated) & (kCacheLineSize - 1); + // malloc is at least kPointerSize aligned, so we can store the "allocated" + // pointer immediately before the aligned memory. + assert(misalignment % kPointerSize == 0); + char *const aligned = allocated + kCacheLineSize - misalignment; + memcpy(aligned - kPointerSize, &allocated, kPointerSize); + return BUTTERAUGLI_ASSUME_ALIGNED(aligned, 64); +} + +void CacheAligned::Free(void *aligned_pointer) { + if (aligned_pointer == nullptr) { + return; + } + char *const aligned = static_cast(aligned_pointer); + assert(reinterpret_cast(aligned) % kCacheLineSize == 0); + char *allocated; + memcpy(&allocated, aligned - kPointerSize, kPointerSize); + assert(allocated <= aligned - kPointerSize); + assert(allocated >= aligned - kCacheLineSize); + free(allocated); +} + +static inline bool IsNan(const float x) { + uint32_t bits; + memcpy(&bits, &x, sizeof(bits)); + const uint32_t bitmask_exp = 0x7F800000; + return (bits & bitmask_exp) == bitmask_exp && (bits & 0x7FFFFF); +} + +static inline bool IsNan(const double x) { + uint64_t bits; + memcpy(&bits, &x, sizeof(bits)); + return (0x7ff0000000000001ULL <= bits && bits <= 0x7fffffffffffffffULL) || + (0xfff0000000000001ULL <= bits && bits <= 0xffffffffffffffffULL); +} + +static inline void CheckImage(const ImageF &image, const char *name) { + for (size_t y = 0; y < image.ysize(); ++y) { + const float * const BUTTERAUGLI_RESTRICT row = image.Row(y); + for (size_t x = 0; x < image.xsize(); ++x) { + if (IsNan(row[x])) { + printf("Image %s @ %lu,%lu (of %lu,%lu)\n", name, x, y, image.xsize(), + image.ysize()); + exit(1); + } + } + } +} + +#if BUTTERAUGLI_ENABLE_CHECKS + +#define CHECK_NAN(x, str) \ + do { \ + if (IsNan(x)) { \ + printf("%d: %s\n", __LINE__, str); \ + abort(); \ + } \ + } while (0) + +#define CHECK_IMAGE(image, name) CheckImage(image, name) + +#else + +#define CHECK_NAN(x, str) +#define CHECK_IMAGE(image, name) + +#endif + + +// Purpose of kInternalGoodQualityThreshold: +// Normalize 'ok' image degradation to 1.0 across different versions of +// butteraugli. +static const double kInternalGoodQualityThreshold = 20.35; +static const double kGlobalScale = 1.0 / kInternalGoodQualityThreshold; + +inline float DotProduct(const float u[3], const float v[3]) { + return u[0] * v[0] + u[1] * v[1] + u[2] * v[2]; +} + +std::vector ComputeKernel(float sigma) { + const float m = 2.25; // Accuracy increases when m is increased. + const float scaler = -1.0 / (2 * sigma * sigma); + const int diff = std::max(1, m * fabs(sigma)); + std::vector kernel(2 * diff + 1); + for (int i = -diff; i <= diff; ++i) { + kernel[i + diff] = exp(scaler * i * i); + } + return kernel; +} + +void ConvolveBorderColumn( + const ImageF& in, + const std::vector& kernel, + const float weight_no_border, + const float border_ratio, + const size_t x, + float* const BUTTERAUGLI_RESTRICT row_out) { + const int offset = kernel.size() / 2; + int minx = x < offset ? 0 : x - offset; + int maxx = std::min(in.xsize() - 1, x + offset); + float weight = 0.0f; + for (int j = minx; j <= maxx; ++j) { + weight += kernel[j - x + offset]; + } + // Interpolate linearly between the no-border scaling and border scaling. + weight = (1.0f - border_ratio) * weight + border_ratio * weight_no_border; + float scale = 1.0f / weight; + for (size_t y = 0; y < in.ysize(); ++y) { + const float* const BUTTERAUGLI_RESTRICT row_in = in.Row(y); + float sum = 0.0f; + for (int j = minx; j <= maxx; ++j) { + sum += row_in[j] * kernel[j - x + offset]; + } + row_out[y] = sum * scale; + } +} + +// Computes a horizontal convolution and transposes the result. +ImageF Convolution(const ImageF& in, + const std::vector& kernel, + const float border_ratio) { + ImageF out(in.ysize(), in.xsize()); + const int len = kernel.size(); + const int offset = kernel.size() / 2; + float weight_no_border = 0.0f; + for (int j = 0; j < len; ++j) { + weight_no_border += kernel[j]; + } + float scale_no_border = 1.0f / weight_no_border; + const int border1 = in.xsize() <= offset ? in.xsize() : offset; + const int border2 = in.xsize() - offset; + std::vector scaled_kernel = kernel; + for (int i = 0; i < scaled_kernel.size(); ++i) { + scaled_kernel[i] *= scale_no_border; + } + // left border + for (int x = 0; x < border1; ++x) { + ConvolveBorderColumn(in, kernel, weight_no_border, border_ratio, x, + out.Row(x)); + } + // middle + for (size_t y = 0; y < in.ysize(); ++y) { + const float* const BUTTERAUGLI_RESTRICT row_in = in.Row(y); + for (int x = border1; x < border2; ++x) { + const int d = x - offset; + float* const BUTTERAUGLI_RESTRICT row_out = out.Row(x); + float sum = 0.0f; + for (int j = 0; j < len; ++j) { + sum += row_in[d + j] * scaled_kernel[j]; + } + row_out[y] = sum; + } + } + // right border + for (int x = border2; x < in.xsize(); ++x) { + ConvolveBorderColumn(in, kernel, weight_no_border, border_ratio, x, + out.Row(x)); + } + return out; +} + +// A blur somewhat similar to a 2D Gaussian blur. +// See: https://en.wikipedia.org/wiki/Gaussian_blur +ImageF Blur(const ImageF& in, float sigma, float border_ratio) { + std::vector kernel = ComputeKernel(sigma); + return Convolution(Convolution(in, kernel, border_ratio), + kernel, border_ratio); +} + +// Clamping linear interpolator. +inline double InterpolateClampNegative(const double *array, + int size, double ix) { + if (ix < 0) { + ix = 0; + } + int baseix = static_cast(ix); + double res; + if (baseix >= size - 1) { + res = array[size - 1]; + } else { + double mix = ix - baseix; + int nextix = baseix + 1; + res = array[baseix] + mix * (array[nextix] - array[baseix]); + } + return res; +} + +double GammaMinArg() { + double out0, out1, out2; + OpsinAbsorbance(0.0, 0.0, 0.0, &out0, &out1, &out2); + return std::min(out0, std::min(out1, out2)); +} + +double GammaMaxArg() { + double out0, out1, out2; + OpsinAbsorbance(255.0, 255.0, 255.0, &out0, &out1, &out2); + return std::max(out0, std::max(out1, out2)); +} + +double SimpleGamma(double v) { + static const double kGamma = 0.372322653176; + static const double limit = 37.8000499603; + double bright = v - limit; + if (bright >= 0) { + static const double mul = 0.0950819040934; + v -= bright * mul; + } + { + static const double limit2 = 74.6154406429; + double bright2 = v - limit2; + if (bright2 >= 0) { + static const double mul = 0.01; + v -= bright2 * mul; + } + } + { + static const double limit2 = 82.8505938033; + double bright2 = v - limit2; + if (bright2 >= 0) { + static const double mul = 0.0316722592629; + v -= bright2 * mul; + } + } + { + static const double limit2 = 92.8505938033; + double bright2 = v - limit2; + if (bright2 >= 0) { + static const double mul = 0.221249885752; + v -= bright2 * mul; + } + } + { + static const double limit2 = 102.8505938033; + double bright2 = v - limit2; + if (bright2 >= 0) { + static const double mul = 0.0402547853939; + v -= bright2 * mul; + } + } + { + static const double limit2 = 112.8505938033; + double bright2 = v - limit2; + if (bright2 >= 0) { + static const double mul = 0.021471798711500003; + v -= bright2 * mul; + } + } + static const double offset = 0.106544447664; + static const double scale = 10.7950943969; + double retval = scale * (offset + pow(v, kGamma)); + return retval; +} + +static inline double Gamma(double v) { + //return SimpleGamma(v); + return GammaPolynomial(v); +} + +std::vector OpsinDynamicsImage(const std::vector& rgb) { + PROFILER_FUNC; + std::vector xyb(3); + std::vector blurred(3); + const double kSigma = 1.2; + for (int i = 0; i < 3; ++i) { + xyb[i] = ImageF(rgb[i].xsize(), rgb[i].ysize()); + blurred[i] = Blur(rgb[i], kSigma, 0.0f); + } + for (size_t y = 0; y < rgb[0].ysize(); ++y) { + const float* const BUTTERAUGLI_RESTRICT row_r = rgb[0].Row(y); + const float* const BUTTERAUGLI_RESTRICT row_g = rgb[1].Row(y); + const float* const BUTTERAUGLI_RESTRICT row_b = rgb[2].Row(y); + const float* const BUTTERAUGLI_RESTRICT row_blurred_r = blurred[0].Row(y); + const float* const BUTTERAUGLI_RESTRICT row_blurred_g = blurred[1].Row(y); + const float* const BUTTERAUGLI_RESTRICT row_blurred_b = blurred[2].Row(y); + float* const BUTTERAUGLI_RESTRICT row_out_x = xyb[0].Row(y); + float* const BUTTERAUGLI_RESTRICT row_out_y = xyb[1].Row(y); + float* const BUTTERAUGLI_RESTRICT row_out_b = xyb[2].Row(y); + for (size_t x = 0; x < rgb[0].xsize(); ++x) { + float sensitivity[3]; + { + // Calculate sensitivity based on the smoothed image gamma derivative. + float pre_mixed0, pre_mixed1, pre_mixed2; + OpsinAbsorbance(row_blurred_r[x], row_blurred_g[x], row_blurred_b[x], + &pre_mixed0, &pre_mixed1, &pre_mixed2); + // TODO: use new polynomial to compute Gamma(x)/x derivative. + sensitivity[0] = Gamma(pre_mixed0) / pre_mixed0; + sensitivity[1] = Gamma(pre_mixed1) / pre_mixed1; + sensitivity[2] = Gamma(pre_mixed2) / pre_mixed2; + } + float cur_mixed0, cur_mixed1, cur_mixed2; + OpsinAbsorbance(row_r[x], row_g[x], row_b[x], + &cur_mixed0, &cur_mixed1, &cur_mixed2); + cur_mixed0 *= sensitivity[0]; + cur_mixed1 *= sensitivity[1]; + cur_mixed2 *= sensitivity[2]; + RgbToXyb(cur_mixed0, cur_mixed1, cur_mixed2, + &row_out_x[x], &row_out_y[x], &row_out_b[x]); + } + } + return xyb; +} + +// Make area around zero less important (remove it). +static BUTTERAUGLI_INLINE float RemoveRangeAroundZero(float w, float x) { + return x > w ? x - w : x < -w ? x + w : 0.0f; +} + +// Make area around zero more important (2x it until the limit). +static BUTTERAUGLI_INLINE float AmplifyRangeAroundZero(float w, float x) { + return x > w ? x + w : x < -w ? x - w : 2.0f * x; +} + +// XybLowFreqToVals converts from low-frequency XYB space to the 'vals' space. +// Vals space can be converted to L2-norm space (Euclidean and normalized) +// through visual masking. +template +BUTTERAUGLI_INLINE void XybLowFreqToVals(const V &x, const V &y, const V &b_arg, + V *BUTTERAUGLI_RESTRICT valx, + V *BUTTERAUGLI_RESTRICT valy, + V *BUTTERAUGLI_RESTRICT valb) { + static const double xmuli = 5.57547552483; + static const double ymuli = 1.20828034498; + static const double bmuli = 6.08319517575; + static const double y_to_b_muli = -0.628811683685; + + const V xmul(xmuli); + const V ymul(ymuli); + const V bmul(bmuli); + const V y_to_b_mul(y_to_b_muli); + const V b = b_arg + y_to_b_mul * y; + *valb = b * bmul; + *valx = x * xmul; + *valy = y * ymul; +} + +static ImageF SuppressInBrightAreas(size_t xsize, size_t ysize, + double mul, double mul2, double reg, + const ImageF& hf, + const ImageF& brightness) { + ImageF inew(xsize, ysize); + for (size_t y = 0; y < ysize; ++y) { + const float* const rowhf = hf.Row(y); + const float* const rowbr = brightness.Row(y); + float* const rownew = inew.Row(y); + for (size_t x = 0; x < xsize; ++x) { + float v = rowhf[x]; + float scaler = mul * reg / (reg + rowbr[x]); + rownew[x] = scaler * v; + } + } + return inew; +} + + +static float SuppressHfInBrightAreas(float hf, float brightness, + float mul, float reg) { + float scaler = mul * reg / (reg + brightness); + return scaler * hf; +} + +static float SuppressUhfInBrightAreas(float hf, float brightness, + float mul, float reg) { + float scaler = mul * reg / (reg + brightness); + return scaler * hf; +} + +static float MaximumClamp(float v, float maxval) { + static const double kMul = 0.688059627878; + if (v >= maxval) { + v -= maxval; + v *= kMul; + v += maxval; + } else if (v < -maxval) { + v += maxval; + v *= kMul; + v -= maxval; + } + return v; +} + +static ImageF MaximumClamping(size_t xsize, size_t ysize, const ImageF& ix, + double yw) { + static const double kMul = 0.688059627878; + ImageF inew(xsize, ysize); + for (size_t y = 0; y < ysize; ++y) { + const float* const rowx = ix.Row(y); + float* const rownew = inew.Row(y); + for (size_t x = 0; x < xsize; ++x) { + double v = rowx[x]; + if (v >= yw) { + v -= yw; + v *= kMul; + v += yw; + } else if (v < -yw) { + v += yw; + v *= kMul; + v -= yw; + } + rownew[x] = v; + } + } + return inew; +} + +static ImageF SuppressXByY(size_t xsize, size_t ysize, + const ImageF& ix, const ImageF& iy, + const double yw) { + static const double s = 0.745954517135; + ImageF inew(xsize, ysize); + for (size_t y = 0; y < ysize; ++y) { + const float* const rowx = ix.Row(y); + const float* const rowy = iy.Row(y); + float* const rownew = inew.Row(y); + for (size_t x = 0; x < xsize; ++x) { + const double xval = rowx[x]; + const double yval = rowy[x]; + const double scaler = s + (yw * (1.0 - s)) / (yw + yval * yval); + rownew[x] = scaler * xval; + } + } + return inew; +} + +static void SeparateFrequencies( + size_t xsize, size_t ysize, + const std::vector& xyb, + PsychoImage &ps) { + PROFILER_FUNC; + ps.lf.resize(3); // XYB + ps.mf.resize(3); // XYB + ps.hf.resize(2); // XY + ps.uhf.resize(2); // XY + // Extract lf ... + static const double kSigmaLf = 7.46953768697; + static const double kSigmaHf = 3.734768843485; + static const double kSigmaUhf = 1.8673844217425; + // At borders we move some more of the energy to the high frequency + // parts, because there can be unfortunate continuations in tiling + // background color etc. So we want to represent the borders with + // some more accuracy. + static double border_lf = -0.00457628248637; + static double border_mf = -0.271277366628; + static double border_hf = 0.147068973249; + for (int i = 0; i < 3; ++i) { + ps.lf[i] = Blur(xyb[i], kSigmaLf, border_lf); + // ... and keep everything else in mf. + ps.mf[i] = ImageF(xsize, ysize); + for (size_t y = 0; y < ysize; ++y) { + for (size_t x = 0; x < xsize; ++x) { + ps.mf[i].Row(y)[x] = xyb[i].Row(y)[x] - ps.lf[i].Row(y)[x]; + } + } + if (i == 2) { + ps.mf[i] = Blur(ps.mf[i], kSigmaHf, border_mf); + break; + } + // Divide mf into mf and hf. + ps.hf[i] = ImageF(xsize, ysize); + for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_mf = ps.mf[i].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[i].Row(y); + for (size_t x = 0; x < xsize; ++x) { + row_hf[x] = row_mf[x]; + } + } + ps.mf[i] = Blur(ps.mf[i], kSigmaHf, border_mf); + static const double w0 = 0.120079806822; + static const double w1 = 0.03430529365; + if (i == 0) { + for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_mf = ps.mf[0].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[0].Row(y); + for (size_t x = 0; x < xsize; ++x) { + row_hf[x] -= row_mf[x]; + row_mf[x] = RemoveRangeAroundZero(w0, row_mf[x]); + } + } + } else { + for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_mf = ps.mf[1].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[1].Row(y); + for (size_t x = 0; x < xsize; ++x) { + row_hf[x] -= row_mf[x]; + row_mf[x] = AmplifyRangeAroundZero(w1, row_mf[x]); + } + } + } + } + // Suppress red-green by intensity change in the high freq channels. + static const double suppress = 2.96534974403; + ps.hf[0] = SuppressXByY(xsize, ysize, ps.hf[0], ps.hf[1], suppress); + + for (int i = 0; i < 2; ++i) { + // Divide hf into hf and uhf. + ps.uhf[i] = ImageF(xsize, ysize); + for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_uhf = ps.uhf[i].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[i].Row(y); + for (size_t x = 0; x < xsize; ++x) { + row_uhf[x] = row_hf[x]; + } + } + ps.hf[i] = Blur(ps.hf[i], kSigmaUhf, border_hf); + static const double kRemoveHfRange = 0.0287615200377; + static const double kMaxclampHf = 78.8223237675; + static const double kMaxclampUhf = 5.8907152736; + static const float kMulSuppressHf = 1.10684769012; + static const float kMulRegHf = 0.478741530298; + static const float kRegHf = 2000 * kMulRegHf; + static const float kMulSuppressUhf = 1.76905001176; + static const float kMulRegUhf = 0.310148420674; + static const float kRegUhf = 2000 * kMulRegUhf; + + if (i == 0) { + for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_uhf = ps.uhf[0].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[0].Row(y); + for (size_t x = 0; x < xsize; ++x) { + row_uhf[x] -= row_hf[x]; + row_hf[x] = RemoveRangeAroundZero(kRemoveHfRange, row_hf[x]); + } + } + } else { + for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_uhf = ps.uhf[1].Row(y); + float* BUTTERAUGLI_RESTRICT const row_hf = ps.hf[1].Row(y); + float* BUTTERAUGLI_RESTRICT const row_lf = ps.lf[1].Row(y); + for (size_t x = 0; x < xsize; ++x) { + row_uhf[x] -= row_hf[x]; + row_hf[x] = MaximumClamp(row_hf[x], kMaxclampHf); + row_uhf[x] = MaximumClamp(row_uhf[x], kMaxclampUhf); + row_uhf[x] = SuppressUhfInBrightAreas(row_uhf[x], row_lf[x], + kMulSuppressUhf, kRegUhf); + row_hf[x] = SuppressHfInBrightAreas(row_hf[x], row_lf[x], + kMulSuppressHf, kRegHf); + + } + } + } + } + // Modify range around zero code only concerns the high frequency + // planes and only the X and Y channels. + // Convert low freq xyb to vals space so that we can do a simple squared sum + // diff on the low frequencies later. + for (size_t y = 0; y < ysize; ++y) { + float* BUTTERAUGLI_RESTRICT const row_x = ps.lf[0].Row(y); + float* BUTTERAUGLI_RESTRICT const row_y = ps.lf[1].Row(y); + float* BUTTERAUGLI_RESTRICT const row_b = ps.lf[2].Row(y); + for (size_t x = 0; x < xsize; ++x) { + float valx, valy, valb; + XybLowFreqToVals(row_x[x], row_y[x], row_b[x], &valx, &valy, &valb); + row_x[x] = valx; + row_y[x] = valy; + row_b[x] = valb; + } + } +} + +static void SameNoiseLevels(const ImageF& i0, const ImageF& i1, + const double kSigma, + const double w, + const double maxclamp, + ImageF* BUTTERAUGLI_RESTRICT diffmap) { + ImageF blurred(i0.xsize(), i0.ysize()); + for (size_t y = 0; y < i0.ysize(); ++y) { + const float* BUTTERAUGLI_RESTRICT const row0 = i0.Row(y); + const float* BUTTERAUGLI_RESTRICT const row1 = i1.Row(y); + float* BUTTERAUGLI_RESTRICT const to = blurred.Row(y); + for (size_t x = 0; x < i0.xsize(); ++x) { + double v0 = fabs(row0[x]); + double v1 = fabs(row1[x]); + if (v0 > maxclamp) v0 = maxclamp; + if (v1 > maxclamp) v1 = maxclamp; + to[x] = v0 - v1; + } + + } + blurred = Blur(blurred, kSigma, 0.0); + for (size_t y = 0; y < i0.ysize(); ++y) { + const float* BUTTERAUGLI_RESTRICT const row = blurred.Row(y); + float* BUTTERAUGLI_RESTRICT const row_diff = diffmap->Row(y); + for (size_t x = 0; x < i0.xsize(); ++x) { + double diff = row[x]; + row_diff[x] += w * diff * diff; + } + } +} + +static void L2Diff(const ImageF& i0, const ImageF& i1, const double w, + ImageF* BUTTERAUGLI_RESTRICT diffmap) { + if (w == 0) { + return; + } + for (size_t y = 0; y < i0.ysize(); ++y) { + const float* BUTTERAUGLI_RESTRICT const row0 = i0.Row(y); + const float* BUTTERAUGLI_RESTRICT const row1 = i1.Row(y); + float* BUTTERAUGLI_RESTRICT const row_diff = diffmap->Row(y); + for (size_t x = 0; x < i0.xsize(); ++x) { + double diff = row0[x] - row1[x]; + row_diff[x] += w * diff * diff; + } + } +} + +// i0 is the original image. +// i1 is the deformed copy. +static void L2DiffAsymmetric(const ImageF& i0, const ImageF& i1, + double w_0gt1, + double w_0lt1, + ImageF* BUTTERAUGLI_RESTRICT diffmap) { + if (w_0gt1 == 0 && w_0lt1 == 0) { + return; + } + w_0gt1 *= 0.8; + w_0lt1 *= 0.8; + for (size_t y = 0; y < i0.ysize(); ++y) { + const float* BUTTERAUGLI_RESTRICT const row0 = i0.Row(y); + const float* BUTTERAUGLI_RESTRICT const row1 = i1.Row(y); + float* BUTTERAUGLI_RESTRICT const row_diff = diffmap->Row(y); + for (size_t x = 0; x < i0.xsize(); ++x) { + // Primary symmetric quadratic objective. + double diff = row0[x] - row1[x]; + row_diff[x] += w_0gt1 * diff * diff; + + // Secondary half-open quadratic objectives. + const double fabs0 = fabs(row0[x]); + const double too_small = 0.4 * fabs0; + const double too_big = 1.0 * fabs0; + + if (row0[x] < 0) { + if (row1[x] > -too_small) { + double v = row1[x] + too_small; + row_diff[x] += w_0lt1 * v * v; + } else if (row1[x] < -too_big) { + double v = -row1[x] - too_big; + row_diff[x] += w_0lt1 * v * v; + } + } else { + if (row1[x] < too_small) { + double v = too_small - row1[x]; + row_diff[x] += w_0lt1 * v * v; + } else if (row1[x] > too_big) { + double v = row1[x] - too_big; + row_diff[x] += w_0lt1 * v * v; + } + } + } + } +} + +// Making a cluster of local errors to be more impactful than +// just a single error. +ImageF CalculateDiffmap(const ImageF& diffmap_in) { + PROFILER_FUNC; + // Take square root. + ImageF diffmap(diffmap_in.xsize(), diffmap_in.ysize()); + static const float kInitialSlope = 100.0f; + for (size_t y = 0; y < diffmap.ysize(); ++y) { + const float* const BUTTERAUGLI_RESTRICT row_in = diffmap_in.Row(y); + float* const BUTTERAUGLI_RESTRICT row_out = diffmap.Row(y); + for (size_t x = 0; x < diffmap.xsize(); ++x) { + const float orig_val = row_in[x]; + // TODO(b/29974893): Until that is fixed do not call sqrt on very small + // numbers. + row_out[x] = (orig_val < (1.0f / (kInitialSlope * kInitialSlope)) + ? kInitialSlope * orig_val + : std::sqrt(orig_val)); + } + } + { + static const double kSigma = 1.72547472444; + static const double mul1 = 0.458794906198; + static const float scale = 1.0f / (1.0f + mul1); + static const double border_ratio = 1.0; // 2.01209066992; + ImageF blurred = Blur(diffmap, kSigma, border_ratio); + for (int y = 0; y < diffmap.ysize(); ++y) { + const float* const BUTTERAUGLI_RESTRICT row_blurred = blurred.Row(y); + float* const BUTTERAUGLI_RESTRICT row = diffmap.Row(y); + for (int x = 0; x < diffmap.xsize(); ++x) { + row[x] += mul1 * row_blurred[x]; + row[x] *= scale; + } + } + } + return diffmap; +} + +void MaskPsychoImage(const PsychoImage& pi0, const PsychoImage& pi1, + const size_t xsize, const size_t ysize, + std::vector* BUTTERAUGLI_RESTRICT mask, + std::vector* BUTTERAUGLI_RESTRICT mask_dc) { + std::vector mask_xyb0 = CreatePlanes(xsize, ysize, 3); + std::vector mask_xyb1 = CreatePlanes(xsize, ysize, 3); + static const double muls[4] = { + 0, + 1.64178305129, + 0.831081703362, + 3.23680933546, + }; + for (int i = 0; i < 2; ++i) { + double a = muls[2 * i]; + double b = muls[2 * i + 1]; + for (size_t y = 0; y < ysize; ++y) { + const float* const BUTTERAUGLI_RESTRICT row_hf0 = pi0.hf[i].Row(y); + const float* const BUTTERAUGLI_RESTRICT row_hf1 = pi1.hf[i].Row(y); + const float* const BUTTERAUGLI_RESTRICT row_uhf0 = pi0.uhf[i].Row(y); + const float* const BUTTERAUGLI_RESTRICT row_uhf1 = pi1.uhf[i].Row(y); + float* const BUTTERAUGLI_RESTRICT row0 = mask_xyb0[i].Row(y); + float* const BUTTERAUGLI_RESTRICT row1 = mask_xyb1[i].Row(y); + for (size_t x = 0; x < xsize; ++x) { + row0[x] = a * row_uhf0[x] + b * row_hf0[x]; + row1[x] = a * row_uhf1[x] + b * row_hf1[x]; + } + } + } + Mask(mask_xyb0, mask_xyb1, mask, mask_dc); +} + +ButteraugliComparator::ButteraugliComparator(const std::vector& rgb0) + : xsize_(rgb0[0].xsize()), + ysize_(rgb0[0].ysize()), + num_pixels_(xsize_ * ysize_) { + if (xsize_ < 8 || ysize_ < 8) return; + std::vector xyb0 = OpsinDynamicsImage(rgb0); + SeparateFrequencies(xsize_, ysize_, xyb0, pi0_); +} + +void ButteraugliComparator::Mask( + std::vector* BUTTERAUGLI_RESTRICT mask, + std::vector* BUTTERAUGLI_RESTRICT mask_dc) const { + MaskPsychoImage(pi0_, pi0_, xsize_, ysize_, mask, mask_dc); +} + +void ButteraugliComparator::Diffmap(const std::vector& rgb1, + ImageF &result) const { + PROFILER_FUNC; + if (xsize_ < 8 || ysize_ < 8) return; + DiffmapOpsinDynamicsImage(OpsinDynamicsImage(rgb1), result); +} + +void ButteraugliComparator::DiffmapOpsinDynamicsImage( + const std::vector& xyb1, + ImageF &result) const { + PROFILER_FUNC; + if (xsize_ < 8 || ysize_ < 8) return; + PsychoImage pi1; + SeparateFrequencies(xsize_, ysize_, xyb1, pi1); + result = ImageF(xsize_, ysize_); + DiffmapPsychoImage(pi1, result); +} + +void ButteraugliComparator::DiffmapPsychoImage(const PsychoImage& pi1, + ImageF& result) const { + PROFILER_FUNC; + const float hf_asymmetry_ = 0.8f; + if (xsize_ < 8 || ysize_ < 8) { + return; + } + std::vector block_diff_dc(3); + std::vector block_diff_ac(3); + for (int c = 0; c < 3; ++c) { + block_diff_dc[c] = ImageF(xsize_, ysize_, 0.0); + block_diff_ac[c] = ImageF(xsize_, ysize_, 0.0); + } + + static const double wUhfMalta = 5.1409625726; + static const double norm1Uhf = 58.5001247061; + MaltaDiffMap(pi0_.uhf[1], pi1.uhf[1], + wUhfMalta * hf_asymmetry_, + wUhfMalta / hf_asymmetry_, + norm1Uhf, + &block_diff_ac[1]); + + static const double wUhfMaltaX = 4.91743441556; + static const double norm1UhfX = 687196.39002; + MaltaDiffMap(pi0_.uhf[0], pi1.uhf[0], + wUhfMaltaX * hf_asymmetry_, + wUhfMaltaX / hf_asymmetry_, + norm1UhfX, + &block_diff_ac[0]); + + static const double wHfMalta = 153.671655716; + static const double norm1Hf = 83150785.9592; + MaltaDiffMapLF(pi0_.hf[1], pi1.hf[1], + wHfMalta * sqrt(hf_asymmetry_), + wHfMalta / sqrt(hf_asymmetry_), + norm1Hf, + &block_diff_ac[1]); + + static const double wHfMaltaX = 668.358918152; + static const double norm1HfX = 0.882954368025; + MaltaDiffMapLF(pi0_.hf[0], pi1.hf[0], + wHfMaltaX * sqrt(hf_asymmetry_), + wHfMaltaX / sqrt(hf_asymmetry_), + norm1HfX, + &block_diff_ac[0]); + + static const double wMfMalta = 6841.81248144; + static const double norm1Mf = 0.0135134962487; + MaltaDiffMapLF(pi0_.mf[1], pi1.mf[1], wMfMalta, wMfMalta, norm1Mf, + &block_diff_ac[1]); + + static const double wMfMaltaX = 813.901703816; + static const double norm1MfX = 16792.9322251; + MaltaDiffMapLF(pi0_.mf[0], pi1.mf[0], wMfMaltaX, wMfMaltaX, norm1MfX, + &block_diff_ac[0]); + + static const double wmul[9] = { + 0, + 32.4449876135, + 0, + 0, + 0, + 0, + 1.01370836411, + 0, + 1.74566011615, + }; + + static const double maxclamp = 85.7047444518; + static const double kSigmaHfX = 10.6666499623; + static const double w = 884.809801415; + SameNoiseLevels(pi0_.hf[1], pi1.hf[1], kSigmaHfX, w, maxclamp, + &block_diff_ac[1]); + + for (int c = 0; c < 3; ++c) { + if (c < 2) { + L2DiffAsymmetric(pi0_.hf[c], pi1.hf[c], + wmul[c] * hf_asymmetry_, + wmul[c] / hf_asymmetry_, + &block_diff_ac[c]); + } + L2Diff(pi0_.mf[c], pi1.mf[c], wmul[3 + c], &block_diff_ac[c]); + L2Diff(pi0_.lf[c], pi1.lf[c], wmul[6 + c], &block_diff_dc[c]); + } + + std::vector mask_xyb; + std::vector mask_xyb_dc; + MaskPsychoImage(pi0_, pi1, xsize_, ysize_, &mask_xyb, &mask_xyb_dc); + + result = CalculateDiffmap( + CombineChannels(mask_xyb, mask_xyb_dc, block_diff_dc, block_diff_ac)); +} + +// Allows PaddedMaltaUnit to call either function via overloading. +struct MaltaTagLF {}; +struct MaltaTag {}; + +static float MaltaUnit(MaltaTagLF, const float* BUTTERAUGLI_RESTRICT d, + const int xs) { + const int xs3 = 3 * xs; + float retval = 0; + { + // x grows, y constant + float sum = + d[-4] + + d[-2] + + d[0] + + d[2] + + d[4]; + retval += sum * sum; + } + { + // y grows, x constant + float sum = + d[-xs3 - xs] + + d[-xs - xs] + + d[0] + + d[xs + xs] + + d[xs3 + xs]; + retval += sum * sum; + } + { + // both grow + float sum = + d[-xs3 - 3] + + d[-xs - xs - 2] + + d[0] + + d[xs + xs + 2] + + d[xs3 + 3]; + retval += sum * sum; + } + { + // y grows, x shrinks + float sum = + d[-xs3 + 3] + + d[-xs - xs + 2] + + d[0] + + d[xs + xs - 2] + + d[xs3 - 3]; + retval += sum * sum; + } + { + // y grows -4 to 4, x shrinks 1 -> -1 + float sum = + d[-xs3 - xs + 1] + + d[-xs - xs + 1] + + d[0] + + d[xs + xs - 1] + + d[xs3 + xs - 1]; + retval += sum * sum; + } + { + // y grows -4 to 4, x grows -1 -> 1 + float sum = + d[-xs3 - xs - 1] + + d[-xs - xs - 1] + + d[0] + + d[xs + xs + 1] + + d[xs3 + xs + 1]; + retval += sum * sum; + } + { + // x grows -4 to 4, y grows -1 to 1 + float sum = + d[-4 - xs] + + d[-2 - xs] + + d[0] + + d[2 + xs] + + d[4 + xs]; + retval += sum * sum; + } + { + // x grows -4 to 4, y shrinks 1 to -1 + float sum = + d[-4 + xs] + + d[-2 + xs] + + d[0] + + d[2 - xs] + + d[4 - xs]; + retval += sum * sum; + } + { + /* 0_________ + 1__*______ + 2___*_____ + 3_________ + 4____0____ + 5_________ + 6_____*___ + 7______*__ + 8_________ */ + float sum = + d[-xs3 - 2] + + d[-xs - xs - 1] + + d[0] + + d[xs + xs + 1] + + d[xs3 + 2]; + retval += sum * sum; + } + { + /* 0_________ + 1______*__ + 2_____*___ + 3_________ + 4____0____ + 5_________ + 6___*_____ + 7__*______ + 8_________ */ + float sum = + d[-xs3 + 2] + + d[-xs - xs + 1] + + d[0] + + d[xs + xs - 1] + + d[xs3 - 2]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2_*_______ + 3__*______ + 4____0____ + 5______*__ + 6_______*_ + 7_________ + 8_________ */ + float sum = + d[-xs - xs - 3] + + d[-xs - 2] + + d[0] + + d[xs + 2] + + d[xs + xs + 3]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2_______*_ + 3______*__ + 4____0____ + 5__*______ + 6_*_______ + 7_________ + 8_________ */ + float sum = + d[-xs - xs + 3] + + d[-xs + 2] + + d[0] + + d[xs - 2] + + d[xs + xs - 3]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2________* + 3______*__ + 4____0____ + 5__*______ + 6*________ + 7_________ + 8_________ */ + + float sum = + d[xs + xs - 4] + + d[xs - 2] + + d[0] + + d[-xs + 2] + + d[-xs - xs + 4]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2*________ + 3__*______ + 4____0____ + 5______*__ + 6________* + 7_________ + 8_________ */ + float sum = + d[-xs - xs - 4] + + d[-xs - 2] + + d[0] + + d[xs + 2] + + d[xs + xs + 4]; + retval += sum * sum; + } + { + /* 0__*______ + 1_________ + 2___*_____ + 3_________ + 4____0____ + 5_________ + 6_____*___ + 7_________ + 8______*__ */ + float sum = + d[-xs3 - xs - 2] + + d[-xs - xs - 1] + + d[0] + + d[xs + xs + 1] + + d[xs3 + xs + 2]; + retval += sum * sum; + } + { + /* 0______*__ + 1_________ + 2_____*___ + 3_________ + 4____0____ + 5_________ + 6___*_____ + 7_________ + 8__*______ */ + float sum = + d[-xs3 - xs + 2] + + d[-xs - xs + 1] + + d[0] + + d[xs + xs - 1] + + d[xs3 + xs - 2]; + retval += sum * sum; + } + return retval; +} + +static float MaltaUnit(MaltaTag, const float* BUTTERAUGLI_RESTRICT d, + const int xs) { + const int xs3 = 3 * xs; + float retval = 0; + { + // x grows, y constant + float sum = + d[-4] + + d[-3] + + d[-2] + + d[-1] + + d[0] + + d[1] + + d[2] + + d[3] + + d[4]; + retval += sum * sum; + } + { + // y grows, x constant + float sum = + d[-xs3 - xs] + + d[-xs3] + + d[-xs - xs] + + d[-xs] + + d[0] + + d[xs] + + d[xs + xs] + + d[xs3] + + d[xs3 + xs]; + retval += sum * sum; + } + { + // both grow + float sum = + d[-xs3 - 3] + + d[-xs - xs - 2] + + d[-xs - 1] + + d[0] + + d[xs + 1] + + d[xs + xs + 2] + + d[xs3 + 3]; + retval += sum * sum; + } + { + // y grows, x shrinks + float sum = + d[-xs3 + 3] + + d[-xs - xs + 2] + + d[-xs + 1] + + d[0] + + d[xs - 1] + + d[xs + xs - 2] + + d[xs3 - 3]; + retval += sum * sum; + } + { + // y grows -4 to 4, x shrinks 1 -> -1 + float sum = + d[-xs3 - xs + 1] + + d[-xs3 + 1] + + d[-xs - xs + 1] + + d[-xs] + + d[0] + + d[xs] + + d[xs + xs - 1] + + d[xs3 - 1] + + d[xs3 + xs - 1]; + retval += sum * sum; + } + { + // y grows -4 to 4, x grows -1 -> 1 + float sum = + d[-xs3 - xs - 1] + + d[-xs3 - 1] + + d[-xs - xs - 1] + + d[-xs] + + d[0] + + d[xs] + + d[xs + xs + 1] + + d[xs3 + 1] + + d[xs3 + xs + 1]; + retval += sum * sum; + } + { + // x grows -4 to 4, y grows -1 to 1 + float sum = + d[-4 - xs] + + d[-3 - xs] + + d[-2 - xs] + + d[-1] + + d[0] + + d[1] + + d[2 + xs] + + d[3 + xs] + + d[4 + xs]; + retval += sum * sum; + } + { + // x grows -4 to 4, y shrinks 1 to -1 + float sum = + d[-4 + xs] + + d[-3 + xs] + + d[-2 + xs] + + d[-1] + + d[0] + + d[1] + + d[2 - xs] + + d[3 - xs] + + d[4 - xs]; + retval += sum * sum; + } + { + /* 0_________ + 1__*______ + 2___*_____ + 3___*_____ + 4____0____ + 5_____*___ + 6_____*___ + 7______*__ + 8_________ */ + float sum = + d[-xs3 - 2] + + d[-xs - xs - 1] + + d[-xs - 1] + + d[0] + + d[xs + 1] + + d[xs + xs + 1] + + d[xs3 + 2]; + retval += sum * sum; + } + { + /* 0_________ + 1______*__ + 2_____*___ + 3_____*___ + 4____0____ + 5___*_____ + 6___*_____ + 7__*______ + 8_________ */ + float sum = + d[-xs3 + 2] + + d[-xs - xs + 1] + + d[-xs + 1] + + d[0] + + d[xs - 1] + + d[xs + xs - 1] + + d[xs3 - 2]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2_*_______ + 3__**_____ + 4____0____ + 5_____**__ + 6_______*_ + 7_________ + 8_________ */ + float sum = + d[-xs - xs - 3] + + d[-xs - 2] + + d[-xs - 1] + + d[0] + + d[xs + 1] + + d[xs + 2] + + d[xs + xs + 3]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2_______*_ + 3_____**__ + 4____0____ + 5__**_____ + 6_*_______ + 7_________ + 8_________ */ + float sum = + d[-xs - xs + 3] + + d[-xs + 2] + + d[-xs + 1] + + d[0] + + d[xs - 1] + + d[xs - 2] + + d[xs + xs - 3]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2_________ + 3______**_ + 4____0*___ + 5__**_____ + 6**_______ + 7_________ + 8_________ */ + + float sum = + d[xs + xs - 4] + + d[xs + xs - 3] + + d[xs - 2] + + d[xs - 1] + + d[0] + + d[1] + + d[-xs + 2] + + d[-xs + 3]; + retval += sum * sum; + } + { + /* 0_________ + 1_________ + 2**_______ + 3__**_____ + 4____0*___ + 5______**_ + 6_________ + 7_________ + 8_________ */ + float sum = + d[-xs - xs - 4] + + d[-xs - xs - 3] + + d[-xs - 2] + + d[-xs - 1] + + d[0] + + d[1] + + d[xs + 2] + + d[xs + 3]; + retval += sum * sum; + } + { + /* 0__*______ + 1__*______ + 2___*_____ + 3___*_____ + 4____0____ + 5____*____ + 6_____*___ + 7_____*___ + 8_________ */ + float sum = + d[-xs3 - xs - 2] + + d[-xs3 - 2] + + d[-xs - xs - 1] + + d[-xs - 1] + + d[0] + + d[xs] + + d[xs + xs + 1] + + d[xs3 + 1]; + retval += sum * sum; + } + { + /* 0______*__ + 1______*__ + 2_____*___ + 3_____*___ + 4____0____ + 5____*____ + 6___*_____ + 7___*_____ + 8_________ */ + float sum = + d[-xs3 - xs + 2] + + d[-xs3 + 2] + + d[-xs - xs + 1] + + d[-xs + 1] + + d[0] + + d[xs] + + d[xs + xs - 1] + + d[xs3 - 1]; + retval += sum * sum; + } + return retval; +} + +// Returns MaltaUnit. "fastMode" avoids bounds-checks when x0 and y0 are known +// to be far enough from the image borders. +template +static BUTTERAUGLI_INLINE float PaddedMaltaUnit( + float* const BUTTERAUGLI_RESTRICT diffs, const size_t x0, const size_t y0, + const size_t xsize_, const size_t ysize_) { + int ix0 = y0 * xsize_ + x0; + const float* BUTTERAUGLI_RESTRICT d = &diffs[ix0]; + if (fastMode || + (x0 >= 4 && y0 >= 4 && x0 < (xsize_ - 4) && y0 < (ysize_ - 4))) { + return MaltaUnit(Tag(), d, xsize_); + } + + float borderimage[9 * 9]; + for (int dy = 0; dy < 9; ++dy) { + int y = y0 + dy - 4; + if (y < 0 || y >= ysize_) { + for (int dx = 0; dx < 9; ++dx) { + borderimage[dy * 9 + dx] = 0.0f; + } + } else { + for (int dx = 0; dx < 9; ++dx) { + int x = x0 + dx - 4; + if (x < 0 || x >= xsize_) { + borderimage[dy * 9 + dx] = 0.0f; + } else { + borderimage[dy * 9 + dx] = diffs[y * xsize_ + x]; + } + } + } + } + return MaltaUnit(Tag(), &borderimage[4 * 9 + 4], 9); +} + +template +static void MaltaDiffMapImpl(const ImageF& lum0, const ImageF& lum1, + const size_t xsize_, const size_t ysize_, + const double w_0gt1, + const double w_0lt1, + double norm1, + const double len, const double mulli, + ImageF* block_diff_ac) { + const float kWeight0 = 0.5; + const float kWeight1 = 0.33; + + const double w_pre0gt1 = mulli * sqrt(kWeight0 * w_0gt1) / (len * 2 + 1); + const double w_pre0lt1 = mulli * sqrt(kWeight1 * w_0lt1) / (len * 2 + 1); + const float norm2_0gt1 = w_pre0gt1 * norm1; + const float norm2_0lt1 = w_pre0lt1 * norm1; + + std::vector diffs(ysize_ * xsize_); + for (size_t y = 0, ix = 0; y < ysize_; ++y) { + const float* BUTTERAUGLI_RESTRICT const row0 = lum0.Row(y); + const float* BUTTERAUGLI_RESTRICT const row1 = lum1.Row(y); + for (size_t x = 0; x < xsize_; ++x, ++ix) { + const float absval = 0.5 * std::abs(row0[x]) + 0.5 * std::abs(row1[x]); + const float diff = row0[x] - row1[x]; + const float scaler = norm2_0gt1 / (static_cast(norm1) + absval); + + // Primary symmetric quadratic objective. + diffs[ix] = scaler * diff; + + const float scaler2 = norm2_0lt1 / (static_cast(norm1) + absval); + const double fabs0 = fabs(row0[x]); + + // Secondary half-open quadratic objectives. + const double too_small = 0.55 * fabs0; + const double too_big = 1.05 * fabs0; + + if (row0[x] < 0) { + if (row1[x] > -too_small) { + double impact = scaler2 * (row1[x] + too_small); + if (diff < 0) { + diffs[ix] -= impact; + } else { + diffs[ix] += impact; + } + } else if (row1[x] < -too_big) { + double impact = scaler2 * (-row1[x] - too_big); + if (diff < 0) { + diffs[ix] -= impact; + } else { + diffs[ix] += impact; + } + } + } else { + if (row1[x] < too_small) { + double impact = scaler2 * (too_small - row1[x]); + if (diff < 0) { + diffs[ix] -= impact; + } else { + diffs[ix] += impact; + } + } else if (row1[x] > too_big) { + double impact = scaler2 * (row1[x] - too_big); + if (diff < 0) { + diffs[ix] -= impact; + } else { + diffs[ix] += impact; + } + } + } + } + } + + size_t y0 = 0; + // Top + for (; y0 < 4; ++y0) { + float* const BUTTERAUGLI_RESTRICT row_diff = block_diff_ac->Row(y0); + for (size_t x0 = 0; x0 < xsize_; ++x0) { + row_diff[x0] += + PaddedMaltaUnit(&diffs[0], x0, y0, xsize_, ysize_); + } + } + + // Middle + for (; y0 < ysize_ - 4; ++y0) { + float* const BUTTERAUGLI_RESTRICT row_diff = block_diff_ac->Row(y0); + size_t x0 = 0; + for (; x0 < 4; ++x0) { + row_diff[x0] += + PaddedMaltaUnit(&diffs[0], x0, y0, xsize_, ysize_); + } + for (; x0 < xsize_ - 4; ++x0) { + row_diff[x0] += + PaddedMaltaUnit(&diffs[0], x0, y0, xsize_, ysize_); + } + + for (; x0 < xsize_; ++x0) { + row_diff[x0] += + PaddedMaltaUnit(&diffs[0], x0, y0, xsize_, ysize_); + } + } + + // Bottom + for (; y0 < ysize_; ++y0) { + float* const BUTTERAUGLI_RESTRICT row_diff = block_diff_ac->Row(y0); + for (size_t x0 = 0; x0 < xsize_; ++x0) { + row_diff[x0] += + PaddedMaltaUnit(&diffs[0], x0, y0, xsize_, ysize_); + } + } +} + +void ButteraugliComparator::MaltaDiffMap( + const ImageF& lum0, const ImageF& lum1, + const double w_0gt1, + const double w_0lt1, + const double norm1, ImageF* BUTTERAUGLI_RESTRICT block_diff_ac) const { + PROFILER_FUNC; + const double len = 3.75; + static const double mulli = 0.354191303559; + MaltaDiffMapImpl(lum0, lum1, xsize_, ysize_, w_0gt1, w_0lt1, + norm1, len, + mulli, block_diff_ac); +} + +void ButteraugliComparator::MaltaDiffMapLF( + const ImageF& lum0, const ImageF& lum1, + const double w_0gt1, + const double w_0lt1, + const double norm1, ImageF* BUTTERAUGLI_RESTRICT block_diff_ac) const { + PROFILER_FUNC; + const double len = 3.75; + static const double mulli = 0.405371989604; + MaltaDiffMapImpl(lum0, lum1, xsize_, ysize_, + w_0gt1, w_0lt1, + norm1, len, + mulli, block_diff_ac); +} + +ImageF ButteraugliComparator::CombineChannels( + const std::vector& mask_xyb, + const std::vector& mask_xyb_dc, + const std::vector& block_diff_dc, + const std::vector& block_diff_ac) const { + PROFILER_FUNC; + ImageF result(xsize_, ysize_); + for (size_t y = 0; y < ysize_; ++y) { + float* const BUTTERAUGLI_RESTRICT row_out = result.Row(y); + for (size_t x = 0; x < xsize_; ++x) { + float mask[3]; + float dc_mask[3]; + float diff_dc[3]; + float diff_ac[3]; + for (int i = 0; i < 3; ++i) { + mask[i] = mask_xyb[i].Row(y)[x]; + dc_mask[i] = mask_xyb_dc[i].Row(y)[x]; + diff_dc[i] = block_diff_dc[i].Row(y)[x]; + diff_ac[i] = block_diff_ac[i].Row(y)[x]; + } + row_out[x] = (DotProduct(diff_dc, dc_mask) + DotProduct(diff_ac, mask)); + } + } + return result; +} + +double ButteraugliScoreFromDiffmap(const ImageF& diffmap) { + PROFILER_FUNC; + float retval = 0.0f; + for (size_t y = 0; y < diffmap.ysize(); ++y) { + const float * const BUTTERAUGLI_RESTRICT row = diffmap.Row(y); + for (size_t x = 0; x < diffmap.xsize(); ++x) { + retval = std::max(retval, row[x]); + } + } + return retval; +} + +#include + +// ===== Functions used by Mask only ===== +static std::array MakeMask( + double extmul, double extoff, + double mul, double offset, + double scaler) { + std::array lut; + for (int i = 0; i < lut.size(); ++i) { + const double c = mul / ((0.01 * scaler * i) + offset); + lut[i] = kGlobalScale * (1.0 + extmul * (c + extoff)); + if (lut[i] < 1e-5) { + lut[i] = 1e-5; + } + assert(lut[i] >= 0.0); + lut[i] *= lut[i]; + } + return lut; +} + +double MaskX(double delta) { + static const double extmul = 2.59885507073; + static const double extoff = 3.08805636789; + static const double offset = 0.315424196682; + static const double scaler = 16.2770141832; + static const double mul = 5.62939030582; + static const std::array lut = + MakeMask(extmul, extoff, mul, offset, scaler); + return InterpolateClampNegative(lut.data(), lut.size(), delta); +} + +double MaskY(double delta) { + static const double extmul = 0.9613705131; + static const double extoff = -0.581933100068; + static const double offset = 1.00846207765; + static const double scaler = 2.2342321176; + static const double mul = 6.64307621174; + static const std::array lut = + MakeMask(extmul, extoff, mul, offset, scaler); + return InterpolateClampNegative(lut.data(), lut.size(), delta); +} + +double MaskDcX(double delta) { + static const double extmul = 10.0470705878; + static const double extoff = 3.18472654033; + static const double offset = 0.0551512255218; + static const double scaler = 70.0; + static const double mul = 0.373092999662; + static const std::array lut = + MakeMask(extmul, extoff, mul, offset, scaler); + return InterpolateClampNegative(lut.data(), lut.size(), delta); +} + +double MaskDcY(double delta) { + static const double extmul = 0.0115640939227; + static const double extoff = 45.9483175519; + static const double offset = 0.0142290066313; + static const double scaler = 5.0; + static const double mul = 2.52611324247; + static const std::array lut = + MakeMask(extmul, extoff, mul, offset, scaler); + return InterpolateClampNegative(lut.data(), lut.size(), delta); +} + +ImageF DiffPrecompute(const ImageF& xyb0, const ImageF& xyb1) { + PROFILER_FUNC; + const size_t xsize = xyb0.xsize(); + const size_t ysize = xyb0.ysize(); + ImageF result(xsize, ysize); + size_t x2, y2; + for (size_t y = 0; y < ysize; ++y) { + if (y + 1 < ysize) { + y2 = y + 1; + } else if (y > 0) { + y2 = y - 1; + } else { + y2 = y; + } + const float* const BUTTERAUGLI_RESTRICT row0_in = xyb0.Row(y); + const float* const BUTTERAUGLI_RESTRICT row1_in = xyb1.Row(y); + const float* const BUTTERAUGLI_RESTRICT row0_in2 = xyb0.Row(y2); + const float* const BUTTERAUGLI_RESTRICT row1_in2 = xyb1.Row(y2); + float* const BUTTERAUGLI_RESTRICT row_out = result.Row(y); + for (size_t x = 0; x < xsize; ++x) { + if (x + 1 < xsize) { + x2 = x + 1; + } else if (x > 0) { + x2 = x - 1; + } else { + x2 = x; + } + double sup0 = (fabs(row0_in[x] - row0_in[x2]) + + fabs(row0_in[x] - row0_in2[x])); + double sup1 = (fabs(row1_in[x] - row1_in[x2]) + + fabs(row1_in[x] - row1_in2[x])); + static const double mul0 = 0.918416534734; + row_out[x] = mul0 * std::min(sup0, sup1); + static const double cutoff = 55.0184555849; + if (row_out[x] >= cutoff) { + row_out[x] = cutoff; + } + } + } + return result; +} + +void Mask(const std::vector& xyb0, + const std::vector& xyb1, + std::vector* BUTTERAUGLI_RESTRICT mask, + std::vector* BUTTERAUGLI_RESTRICT mask_dc) { + PROFILER_FUNC; + const size_t xsize = xyb0[0].xsize(); + const size_t ysize = xyb0[0].ysize(); + mask->resize(3); + *mask_dc = CreatePlanes(xsize, ysize, 3); + double muls[2] = { + 0.207017089891, + 0.267138152891, + }; + double normalizer = { + 1.0 / (muls[0] + muls[1]), + }; + static const double r0 = 2.3770330432; + static const double r1 = 9.04353323561; + static const double r2 = 9.24456601467; + static const double border_ratio = -0.0724948220913; + + { + // X component + ImageF diff = DiffPrecompute(xyb0[0], xyb1[0]); + ImageF blurred = Blur(diff, r2, border_ratio); + (*mask)[0] = ImageF(xsize, ysize); + for (size_t y = 0; y < ysize; ++y) { + for (size_t x = 0; x < xsize; ++x) { + (*mask)[0].Row(y)[x] = blurred.Row(y)[x]; + } + } + } + { + // Y component + (*mask)[1] = ImageF(xsize, ysize); + ImageF diff = DiffPrecompute(xyb0[1], xyb1[1]); + ImageF blurred1 = Blur(diff, r0, border_ratio); + ImageF blurred2 = Blur(diff, r1, border_ratio); + for (size_t y = 0; y < ysize; ++y) { + for (size_t x = 0; x < xsize; ++x) { + const double val = normalizer * ( + muls[0] * blurred1.Row(y)[x] + + muls[1] * blurred2.Row(y)[x]); + (*mask)[1].Row(y)[x] = val; + } + } + } + // B component + (*mask)[2] = ImageF(xsize, ysize); + static const double mul[2] = { + 16.6963293877, + 2.1364621982, + }; + static const double w00 = 36.4671237619; + static const double w11 = 2.1887170895; + static const double w_ytob_hf = std::max( + 0.086624184478, + 0.0); + static const double w_ytob_lf = 21.6804277046; + static const double p1_to_p0 = 0.0513061271723; + + for (size_t y = 0; y < ysize; ++y) { + for (size_t x = 0; x < xsize; ++x) { + const double s0 = (*mask)[0].Row(y)[x]; + const double s1 = (*mask)[1].Row(y)[x]; + const double p1 = mul[1] * w11 * s1; + const double p0 = mul[0] * w00 * s0 + p1_to_p0 * p1; + + (*mask)[0].Row(y)[x] = MaskX(p0); + (*mask)[1].Row(y)[x] = MaskY(p1); + (*mask)[2].Row(y)[x] = w_ytob_hf * MaskY(p1); + (*mask_dc)[0].Row(y)[x] = MaskDcX(p0); + (*mask_dc)[1].Row(y)[x] = MaskDcY(p1); + (*mask_dc)[2].Row(y)[x] = w_ytob_lf * MaskDcY(p1); + } + } +} + +void ButteraugliDiffmap(const std::vector &rgb0_image, + const std::vector &rgb1_image, + ImageF &result_image) { + const size_t xsize = rgb0_image[0].xsize(); + const size_t ysize = rgb0_image[0].ysize(); + static const int kMax = 8; + if (xsize < kMax || ysize < kMax) { + // Butteraugli values for small (where xsize or ysize is smaller + // than 8 pixels) images are non-sensical, but most likely it is + // less disruptive to try to compute something than just give up. + // Temporarily extend the borders of the image to fit 8 x 8 size. + int xborder = xsize < kMax ? (kMax - xsize) / 2 : 0; + int yborder = ysize < kMax ? (kMax - ysize) / 2 : 0; + size_t xscaled = std::max(kMax, xsize); + size_t yscaled = std::max(kMax, ysize); + std::vector scaled0 = CreatePlanes(xscaled, yscaled, 3); + std::vector scaled1 = CreatePlanes(xscaled, yscaled, 3); + for (int i = 0; i < 3; ++i) { + for (int y = 0; y < yscaled; ++y) { + for (int x = 0; x < xscaled; ++x) { + size_t x2 = std::min(xsize - 1, std::max(0, x - xborder)); + size_t y2 = std::min(ysize - 1, std::max(0, y - yborder)); + scaled0[i].Row(y)[x] = rgb0_image[i].Row(y2)[x2]; + scaled1[i].Row(y)[x] = rgb1_image[i].Row(y2)[x2]; + } + } + } + ImageF diffmap_scaled; + ButteraugliDiffmap(scaled0, scaled1, diffmap_scaled); + result_image = ImageF(xsize, ysize); + for (int y = 0; y < ysize; ++y) { + for (int x = 0; x < xsize; ++x) { + result_image.Row(y)[x] = diffmap_scaled.Row(y + yborder)[x + xborder]; + } + } + return; + } + ButteraugliComparator butteraugli(rgb0_image); + butteraugli.Diffmap(rgb1_image, result_image); +} + +bool ButteraugliInterface(const std::vector &rgb0, + const std::vector &rgb1, + ImageF &diffmap, + double &diffvalue) { + const size_t xsize = rgb0[0].xsize(); + const size_t ysize = rgb0[0].ysize(); + if (xsize < 1 || ysize < 1) { + return false; // No image. + } + for (int i = 1; i < 3; i++) { + if (rgb0[i].xsize() != xsize || rgb0[i].ysize() != ysize || + rgb1[i].xsize() != xsize || rgb1[i].ysize() != ysize) { + return false; // Image planes must have same dimensions. + } + } + ButteraugliDiffmap(rgb0, rgb1, diffmap); + diffvalue = ButteraugliScoreFromDiffmap(diffmap); + return true; +} + +bool ButteraugliAdaptiveQuantization(size_t xsize, size_t ysize, + const std::vector > &rgb, std::vector &quant) { + if (xsize < 16 || ysize < 16) { + return false; // Butteraugli is undefined for small images. + } + size_t size = xsize * ysize; + + std::vector rgb_planes = PlanesFromPacked(xsize, ysize, rgb); + std::vector scale_xyb; + std::vector scale_xyb_dc; + Mask(rgb_planes, rgb_planes, &scale_xyb, &scale_xyb_dc); + quant.reserve(size); + + // Mask gives us values in 3 color channels, but for now we take only + // the intensity channel. + for (size_t y = 0; y < ysize; ++y) { + for (size_t x = 0; x < xsize; ++x) { + quant.push_back(scale_xyb[1].Row(y)[x]); + } + } + return true; +} + +double ButteraugliFuzzyClass(double score) { + static const double fuzzy_width_up = 6.07887388532; + static const double fuzzy_width_down = 5.50793514384; + static const double m0 = 2.0; + static const double scaler = 0.840253347958; + double val; + if (score < 1.0) { + // val in [scaler .. 2.0] + val = m0 / (1.0 + exp((score - 1.0) * fuzzy_width_down)); + val -= 1.0; // from [1 .. 2] to [0 .. 1] + val *= 2.0 - scaler; // from [0 .. 1] to [0 .. 2.0 - scaler] + val += scaler; // from [0 .. 2.0 - scaler] to [scaler .. 2.0] + } else { + // val in [0 .. scaler] + val = m0 / (1.0 + exp((score - 1.0) * fuzzy_width_up)); + val *= scaler; + } + return val; +} + +double ButteraugliFuzzyInverse(double seek) { + double pos = 0; + for (double range = 1.0; range >= 1e-10; range *= 0.5) { + double cur = ButteraugliFuzzyClass(pos); + if (cur < seek) { + pos -= range; + } else { + pos += range; + } + } + return pos; +} + +namespace { + +void ScoreToRgb(double score, double good_threshold, double bad_threshold, + uint8_t rgb[3]) { + double heatmap[12][3] = { + {0, 0, 0}, + {0, 0, 1}, + {0, 1, 1}, + {0, 1, 0}, // Good level + {1, 1, 0}, + {1, 0, 0}, // Bad level + {1, 0, 1}, + {0.5, 0.5, 1.0}, + {1.0, 0.5, 0.5}, // Pastel colors for the very bad quality range. + {1.0, 1.0, 0.5}, + { + 1, 1, 1, + }, + { + 1, 1, 1, + }, + }; + if (score < good_threshold) { + score = (score / good_threshold) * 0.3; + } else if (score < bad_threshold) { + score = 0.3 + + (score - good_threshold) / (bad_threshold - good_threshold) * 0.15; + } else { + score = 0.45 + (score - bad_threshold) / (bad_threshold * 12) * 0.5; + } + static const int kTableSize = sizeof(heatmap) / sizeof(heatmap[0]); + score = std::min(std::max(score * (kTableSize - 1), 0.0), + kTableSize - 2); + int ix = static_cast(score); + double mix = score - ix; + for (int i = 0; i < 3; ++i) { + double v = mix * heatmap[ix + 1][i] + (1 - mix) * heatmap[ix][i]; + rgb[i] = static_cast(255 * pow(v, 0.5) + 0.5); + } +} + +} // namespace + +void CreateHeatMapImage(const std::vector& distmap, + double good_threshold, double bad_threshold, + size_t xsize, size_t ysize, + std::vector* heatmap) { + heatmap->resize(3 * xsize * ysize); + for (size_t y = 0; y < ysize; ++y) { + for (size_t x = 0; x < xsize; ++x) { + int px = xsize * y + x; + double d = distmap[px]; + uint8_t* rgb = &(*heatmap)[3 * px]; + ScoreToRgb(d, good_threshold, bad_threshold, rgb); + } + } +} + +} // namespace butteraugli diff --git a/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli.h b/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli.h new file mode 100755 index 0000000..2f5d938 --- /dev/null +++ b/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli.h @@ -0,0 +1,619 @@ +// Copyright 2016 Google Inc. All Rights Reserved. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// +// Disclaimer: This is not an official Google product. +// +// Author: Jyrki Alakuijala (jyrki.alakuijala@gmail.com) + +#ifndef BUTTERAUGLI_BUTTERAUGLI_H_ +#define BUTTERAUGLI_BUTTERAUGLI_H_ + +#include +#include +#include +#include +#include +#include +#include +#include + +#define BUTTERAUGLI_ENABLE_CHECKS 0 + +// This is the main interface to butteraugli image similarity +// analysis function. + +namespace butteraugli { + +template +class Image; + +using Image8 = Image; +using ImageF = Image; + +// ButteraugliInterface defines the public interface for butteraugli. +// +// It calculates the difference between rgb0 and rgb1. +// +// rgb0 and rgb1 contain the images. rgb0[c][px] and rgb1[c][px] contains +// the red image for c == 0, green for c == 1, blue for c == 2. Location index +// px is calculated as y * xsize + x. +// +// Value of pixels of images rgb0 and rgb1 need to be represented as raw +// intensity. Most image formats store gamma corrected intensity in pixel +// values. This gamma correction has to be removed, by applying the following +// function: +// butteraugli_val = 255.0 * pow(png_val / 255.0, gamma); +// A typical value of gamma is 2.2. It is usually stored in the image header. +// Take care not to confuse that value with its inverse. The gamma value should +// be always greater than one. +// Butteraugli does not work as intended if the caller does not perform +// gamma correction. +// +// diffmap will contain an image of the size xsize * ysize, containing +// localized differences for values px (indexed with the px the same as rgb0 +// and rgb1). diffvalue will give a global score of similarity. +// +// A diffvalue smaller than kButteraugliGood indicates that images can be +// observed as the same image. +// diffvalue larger than kButteraugliBad indicates that a difference between +// the images can be observed. +// A diffvalue between kButteraugliGood and kButteraugliBad indicates that +// a subtle difference can be observed between the images. +// +// Returns true on success. + +bool ButteraugliInterface(const std::vector &rgb0, + const std::vector &rgb1, + ImageF &diffmap, + double &diffvalue); + +const double kButteraugliQuantLow = 0.26; +const double kButteraugliQuantHigh = 1.454; + +// Converts the butteraugli score into fuzzy class values that are continuous +// at the class boundary. The class boundary location is based on human +// raters, but the slope is arbitrary. Particularly, it does not reflect +// the expectation value of probabilities of the human raters. It is just +// expected that a smoother class boundary will allow for higher-level +// optimization algorithms to work faster. +// +// Returns 2.0 for a perfect match, and 1.0 for 'ok', 0.0 for bad. Because the +// scoring is fuzzy, a butteraugli score of 0.96 would return a class of +// around 1.9. +double ButteraugliFuzzyClass(double score); + +// Input values should be in range 0 (bad) to 2 (good). Use +// kButteraugliNormalization as normalization. +double ButteraugliFuzzyInverse(double seek); + +// Returns a map which can be used for adaptive quantization. Values can +// typically range from kButteraugliQuantLow to kButteraugliQuantHigh. Low +// values require coarse quantization (e.g. near random noise), high values +// require fine quantization (e.g. in smooth bright areas). +bool ButteraugliAdaptiveQuantization(size_t xsize, size_t ysize, + const std::vector > &rgb, std::vector &quant); + +// Implementation details, don't use anything below or your code will +// break in the future. + +#ifdef _MSC_VER +#define BUTTERAUGLI_RESTRICT __restrict +#else +#define BUTTERAUGLI_RESTRICT __restrict__ +#endif + +#ifdef _MSC_VER +#define BUTTERAUGLI_INLINE __forceinline +#else +#define BUTTERAUGLI_INLINE inline +#endif + +#ifdef __clang__ +// Early versions of Clang did not support __builtin_assume_aligned. +#define BUTTERAUGLI_HAS_ASSUME_ALIGNED __has_builtin(__builtin_assume_aligned) +#elif defined(__GNUC__) +#define BUTTERAUGLI_HAS_ASSUME_ALIGNED 1 +#else +#define BUTTERAUGLI_HAS_ASSUME_ALIGNED 0 +#endif + +// Returns a void* pointer which the compiler then assumes is N-byte aligned. +// Example: float* PIK_RESTRICT aligned = (float*)PIK_ASSUME_ALIGNED(in, 32); +// +// The assignment semantics are required by GCC/Clang. ICC provides an in-place +// __assume_aligned, whereas MSVC's __assume appears unsuitable. +#if BUTTERAUGLI_HAS_ASSUME_ALIGNED +#define BUTTERAUGLI_ASSUME_ALIGNED(ptr, align) __builtin_assume_aligned((ptr), (align)) +#else +#define BUTTERAUGLI_ASSUME_ALIGNED(ptr, align) (ptr) +#endif // BUTTERAUGLI_HAS_ASSUME_ALIGNED + +// Functions that depend on the cache line size. +class CacheAligned { + public: + static constexpr size_t kPointerSize = sizeof(void *); + static constexpr size_t kCacheLineSize = 64; + + // The aligned-return annotation is only allowed on function declarations. + static void *Allocate(const size_t bytes); + static void Free(void *aligned_pointer); +}; + +template +using CacheAlignedUniquePtrT = std::unique_ptr; + +using CacheAlignedUniquePtr = CacheAlignedUniquePtrT; + +template +static inline CacheAlignedUniquePtrT Allocate(const size_t entries) { + return CacheAlignedUniquePtrT( + static_cast( + CacheAligned::Allocate(entries * sizeof(T))), + CacheAligned::Free); +} + +// Returns the smallest integer not less than "amount" that is divisible by +// "multiple", which must be a power of two. +template +static inline size_t Align(const size_t amount) { + static_assert(multiple != 0 && ((multiple & (multiple - 1)) == 0), + "Align<> argument must be a power of two"); + return (amount + multiple - 1) & ~(multiple - 1); +} + +// Single channel, contiguous (cache-aligned) rows separated by padding. +// T must be POD. +// +// Rationale: vectorization benefits from aligned operands - unaligned loads and +// especially stores are expensive when the address crosses cache line +// boundaries. Introducing padding after each row ensures the start of a row is +// aligned, and that row loops can process entire vectors (writes to the padding +// are allowed and ignored). +// +// We prefer a planar representation, where channels are stored as separate +// 2D arrays, because that simplifies vectorization (repeating the same +// operation on multiple adjacent components) without the complexity of a +// hybrid layout (8 R, 8 G, 8 B, ...). In particular, clients can easily iterate +// over all components in a row and Image requires no knowledge of the pixel +// format beyond the component type "T". The downside is that we duplicate the +// xsize/ysize members for each channel. +// +// This image layout could also be achieved with a vector and a row accessor +// function, but a class wrapper with support for "deleter" allows wrapping +// existing memory allocated by clients without copying the pixels. It also +// provides convenient accessors for xsize/ysize, which shortens function +// argument lists. Supports move-construction so it can be stored in containers. +template +class Image { + // Returns cache-aligned row stride, being careful to avoid 2K aliasing. + static size_t BytesPerRow(const size_t xsize) { + // Allow reading one extra AVX-2 vector on the right margin. + const size_t row_size = xsize * sizeof(T) + 32; + const size_t align = CacheAligned::kCacheLineSize; + size_t bytes_per_row = (row_size + align - 1) & ~(align - 1); + // During the lengthy window before writes are committed to memory, CPUs + // guard against read after write hazards by checking the address, but + // only the lower 11 bits. We avoid a false dependency between writes to + // consecutive rows by ensuring their sizes are not multiples of 2 KiB. + if (bytes_per_row % 2048 == 0) { + bytes_per_row += align; + } + return bytes_per_row; + } + + public: + using T = ComponentType; + + Image() : xsize_(0), ysize_(0), bytes_per_row_(0), + bytes_(static_cast(nullptr), Ignore) {} + + Image(const size_t xsize, const size_t ysize) + : xsize_(xsize), + ysize_(ysize), + bytes_per_row_(BytesPerRow(xsize)), + bytes_(Allocate(bytes_per_row_ * ysize)) {} + + Image(const size_t xsize, const size_t ysize, T val) + : xsize_(xsize), + ysize_(ysize), + bytes_per_row_(BytesPerRow(xsize)), + bytes_(Allocate(bytes_per_row_ * ysize)) { + for (size_t y = 0; y < ysize_; ++y) { + T* const BUTTERAUGLI_RESTRICT row = Row(y); + for (int x = 0; x < xsize_; ++x) { + row[x] = val; + } + } + } + + Image(const size_t xsize, const size_t ysize, + uint8_t * const BUTTERAUGLI_RESTRICT bytes, + const size_t bytes_per_row) + : xsize_(xsize), + ysize_(ysize), + bytes_per_row_(bytes_per_row), + bytes_(bytes, Ignore) {} + + // Move constructor (required for returning Image from function) + Image(Image &&other) + : xsize_(other.xsize_), + ysize_(other.ysize_), + bytes_per_row_(other.bytes_per_row_), + bytes_(std::move(other.bytes_)) {} + + // Move assignment (required for std::vector) + Image &operator=(Image &&other) { + xsize_ = other.xsize_; + ysize_ = other.ysize_; + bytes_per_row_ = other.bytes_per_row_; + bytes_ = std::move(other.bytes_); + return *this; + } + + void Swap(Image &other) { + std::swap(xsize_, other.xsize_); + std::swap(ysize_, other.ysize_); + std::swap(bytes_per_row_, other.bytes_per_row_); + std::swap(bytes_, other.bytes_); + } + + // How many pixels. + size_t xsize() const { return xsize_; } + size_t ysize() const { return ysize_; } + + T *const BUTTERAUGLI_RESTRICT Row(const size_t y) { +#ifdef BUTTERAUGLI_ENABLE_CHECKS + if (y >= ysize_) { + printf("Row %zu out of bounds (ysize=%zu)\n", y, ysize_); + abort(); + } +#endif + void *row = bytes_.get() + y * bytes_per_row_; + return reinterpret_cast(BUTTERAUGLI_ASSUME_ALIGNED(row, 64)); + } + + const T *const BUTTERAUGLI_RESTRICT Row(const size_t y) const { +#ifdef BUTTERAUGLI_ENABLE_CHECKS + if (y >= ysize_) { + printf("Const row %zu out of bounds (ysize=%zu)\n", y, ysize_); + abort(); + } +#endif + void *row = bytes_.get() + y * bytes_per_row_; + return reinterpret_cast(BUTTERAUGLI_ASSUME_ALIGNED(row, 64)); + } + + // Raw access to byte contents, for interfacing with other libraries. + // Unsigned char instead of char to avoid surprises (sign extension). + uint8_t * const BUTTERAUGLI_RESTRICT bytes() { return bytes_.get(); } + const uint8_t * const BUTTERAUGLI_RESTRICT bytes() const { + return bytes_.get(); + } + size_t bytes_per_row() const { return bytes_per_row_; } + + // Returns number of pixels (some of which are padding) per row. Useful for + // computing other rows via pointer arithmetic. + intptr_t PixelsPerRow() const { + static_assert(CacheAligned::kCacheLineSize % sizeof(T) == 0, + "Padding must be divisible by the pixel size."); + return static_cast(bytes_per_row_ / sizeof(T)); + } + + private: + // Deleter used when bytes are not owned. + static void Ignore(void *ptr) {} + + // (Members are non-const to enable assignment during move-assignment.) + size_t xsize_; // original intended pixels, not including any padding. + size_t ysize_; + size_t bytes_per_row_; // [bytes] including padding. + CacheAlignedUniquePtr bytes_; +}; + +// Returns newly allocated planes of the given dimensions. +template +static inline std::vector> CreatePlanes(const size_t xsize, + const size_t ysize, + const size_t num_planes) { + std::vector> planes; + planes.reserve(num_planes); + for (size_t i = 0; i < num_planes; ++i) { + planes.emplace_back(xsize, ysize); + } + return planes; +} + +// Returns a new image with the same dimensions and pixel values. +template +static inline Image CopyPixels(const Image &other) { + Image copy(other.xsize(), other.ysize()); + const void *BUTTERAUGLI_RESTRICT from = other.bytes(); + void *BUTTERAUGLI_RESTRICT to = copy.bytes(); + memcpy(to, from, other.ysize() * other.bytes_per_row()); + return copy; +} + +// Returns new planes with the same dimensions and pixel values. +template +static inline std::vector> CopyPlanes( + const std::vector> &planes) { + std::vector> copy; + copy.reserve(planes.size()); + for (const Image &plane : planes) { + copy.push_back(CopyPixels(plane)); + } + return copy; +} + +// Compacts a padded image into a preallocated packed vector. +template +static inline void CopyToPacked(const Image &from, std::vector *to) { + const size_t xsize = from.xsize(); + const size_t ysize = from.ysize(); +#if BUTTERAUGLI_ENABLE_CHECKS + if (to->size() < xsize * ysize) { + printf("%zu x %zu exceeds %zu capacity\n", xsize, ysize, to->size()); + abort(); + } +#endif + for (size_t y = 0; y < ysize; ++y) { + const float* const BUTTERAUGLI_RESTRICT row_from = from.Row(y); + float* const BUTTERAUGLI_RESTRICT row_to = to->data() + y * xsize; + memcpy(row_to, row_from, xsize * sizeof(T)); + } +} + +// Expands a packed vector into a preallocated padded image. +template +static inline void CopyFromPacked(const std::vector &from, Image *to) { + const size_t xsize = to->xsize(); + const size_t ysize = to->ysize(); + assert(from.size() == xsize * ysize); + for (size_t y = 0; y < ysize; ++y) { + const float* const BUTTERAUGLI_RESTRICT row_from = + from.data() + y * xsize; + float* const BUTTERAUGLI_RESTRICT row_to = to->Row(y); + memcpy(row_to, row_from, xsize * sizeof(T)); + } +} + +template +static inline std::vector> PlanesFromPacked( + const size_t xsize, const size_t ysize, + const std::vector> &packed) { + std::vector> planes; + planes.reserve(packed.size()); + for (const std::vector &p : packed) { + planes.push_back(Image(xsize, ysize)); + CopyFromPacked(p, &planes.back()); + } + return planes; +} + +template +static inline std::vector> PackedFromPlanes( + const std::vector> &planes) { + assert(!planes.empty()); + const size_t num_pixels = planes[0].xsize() * planes[0].ysize(); + std::vector> packed; + packed.reserve(planes.size()); + for (const Image &image : planes) { + packed.push_back(std::vector(num_pixels)); + CopyToPacked(image, &packed.back()); + } + return packed; +} + +struct PsychoImage { + std::vector uhf; + std::vector hf; + std::vector mf; + std::vector lf; +}; + +class ButteraugliComparator { + public: + ButteraugliComparator(const std::vector& rgb0); + + // Computes the butteraugli map between the original image given in the + // constructor and the distorted image give here. + void Diffmap(const std::vector& rgb1, ImageF& result) const; + + // Same as above, but OpsinDynamicsImage() was already applied. + void DiffmapOpsinDynamicsImage(const std::vector& xyb1, + ImageF& result) const; + + // Same as above, but the frequency decomposition was already applied. + void DiffmapPsychoImage(const PsychoImage& ps1, ImageF &result) const; + + void Mask(std::vector* BUTTERAUGLI_RESTRICT mask, + std::vector* BUTTERAUGLI_RESTRICT mask_dc) const; + + private: + void MaltaDiffMapLF(const ImageF& y0, + const ImageF& y1, + double w_0gt1, + double w_0lt1, + double normalization, + ImageF* BUTTERAUGLI_RESTRICT block_diff_ac) const; + + void MaltaDiffMap(const ImageF& y0, + const ImageF& y1, + double w_0gt1, + double w_0lt1, + double normalization, + ImageF* BUTTERAUGLI_RESTRICT block_diff_ac) const; + + ImageF CombineChannels(const std::vector& scale_xyb, + const std::vector& scale_xyb_dc, + const std::vector& block_diff_dc, + const std::vector& block_diff_ac) const; + + const size_t xsize_; + const size_t ysize_; + const size_t num_pixels_; + PsychoImage pi0_; +}; + +void ButteraugliDiffmap(const std::vector &rgb0, + const std::vector &rgb1, + ImageF &diffmap); + +double ButteraugliScoreFromDiffmap(const ImageF& distmap); + +// Generate rgb-representation of the distance between two images. +void CreateHeatMapImage(const std::vector &distmap, + double good_threshold, double bad_threshold, + size_t xsize, size_t ysize, + std::vector *heatmap); + +// Compute values of local frequency and dc masking based on the activity +// in the two images. +void Mask(const std::vector& xyb0, + const std::vector& xyb1, + std::vector* BUTTERAUGLI_RESTRICT mask, + std::vector* BUTTERAUGLI_RESTRICT mask_dc); + +template +BUTTERAUGLI_INLINE void RgbToXyb(const V &r, const V &g, const V &b, + V *BUTTERAUGLI_RESTRICT valx, + V *BUTTERAUGLI_RESTRICT valy, + V *BUTTERAUGLI_RESTRICT valb) { + *valx = r - g; + *valy = r + g; + *valb = b; +} + +template +BUTTERAUGLI_INLINE void OpsinAbsorbance(const V &in0, const V &in1, + const V &in2, + V *BUTTERAUGLI_RESTRICT out0, + V *BUTTERAUGLI_RESTRICT out1, + V *BUTTERAUGLI_RESTRICT out2) { + // https://en.wikipedia.org/wiki/Photopsin absorbance modeling. + static const double mixi0 = 0.254462330846; + static const double mixi1 = 0.488238255095; + static const double mixi2 = 0.0635278003854; + static const double mixi3 = 1.01681026909; + static const double mixi4 = 0.195214015766; + static const double mixi5 = 0.568019861857; + static const double mixi6 = 0.0860755536007; + static const double mixi7 = 1.1510118369; + static const double mixi8 = 0.07374607900105684; + static const double mixi9 = 0.06142425304154509; + static const double mixi10 = 0.24416850520714256; + static const double mixi11 = 1.20481945273; + + const V mix0(mixi0); + const V mix1(mixi1); + const V mix2(mixi2); + const V mix3(mixi3); + const V mix4(mixi4); + const V mix5(mixi5); + const V mix6(mixi6); + const V mix7(mixi7); + const V mix8(mixi8); + const V mix9(mixi9); + const V mix10(mixi10); + const V mix11(mixi11); + + *out0 = mix0 * in0 + mix1 * in1 + mix2 * in2 + mix3; + *out1 = mix4 * in0 + mix5 * in1 + mix6 * in2 + mix7; + *out2 = mix8 * in0 + mix9 * in1 + mix10 * in2 + mix11; +} + +std::vector OpsinDynamicsImage(const std::vector& rgb); + +ImageF Blur(const ImageF& in, float sigma, float border_ratio); + +double SimpleGamma(double v); + +double GammaMinArg(); +double GammaMaxArg(); + +// Polynomial evaluation via Clenshaw's scheme (similar to Horner's). +// Template enables compile-time unrolling of the recursion, but must reside +// outside of a class due to the specialization. +template +static inline void ClenshawRecursion(const double x, const double *coefficients, + double *b1, double *b2) { + const double x_b1 = x * (*b1); + const double t = (x_b1 + x_b1) - (*b2) + coefficients[INDEX]; + *b2 = *b1; + *b1 = t; + + ClenshawRecursion(x, coefficients, b1, b2); +} + +// Base case +template <> +inline void ClenshawRecursion<0>(const double x, const double *coefficients, + double *b1, double *b2) { + const double x_b1 = x * (*b1); + // The final iteration differs - no 2 * x_b1 here. + *b1 = x_b1 - (*b2) + coefficients[0]; +} + +// Rational polynomial := dividing two polynomial evaluations. These are easier +// to find than minimax polynomials. +struct RationalPolynomial { + template + static double EvaluatePolynomial(const double x, + const double (&coefficients)[N]) { + double b1 = 0.0; + double b2 = 0.0; + ClenshawRecursion(x, coefficients, &b1, &b2); + return b1; + } + + // Evaluates the polynomial at x (in [min_value, max_value]). + inline double operator()(const double x) const { + // First normalize to [0, 1]. + const double x01 = (x - min_value) / (max_value - min_value); + // And then to [-1, 1] domain of Chebyshev polynomials. + const double xc = 2.0 * x01 - 1.0; + + const double yp = EvaluatePolynomial(xc, p); + const double yq = EvaluatePolynomial(xc, q); + if (yq == 0.0) return 0.0; + return static_cast(yp / yq); + } + + // Domain of the polynomials; they are undefined elsewhere. + double min_value; + double max_value; + + // Coefficients of T_n (Chebyshev polynomials of the first kind). + // Degree 5/5 is a compromise between accuracy (0.1%) and numerical stability. + double p[5 + 1]; + double q[5 + 1]; +}; + +static inline double GammaPolynomial(double value) { + static const RationalPolynomial r = { + 0.971783, 590.188894, + { + 98.7821300963361, 164.273222212631, 92.948112871376, + 33.8165311212688, 6.91626704983562, 0.556380877028234 + }, + { + 1, 1.64339473427892, 0.89392405219969, 0.298947051776379, + 0.0507146002577288, 0.00226495093949756 + }}; + return r(value); +} + +} // namespace butteraugli + +#endif // BUTTERAUGLI_BUTTERAUGLI_H_ diff --git a/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli_main.cc b/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli_main.cc new file mode 100755 index 0000000..f38af1d --- /dev/null +++ b/oss-internship-2020/guetzli/third_party/butteraugli/butteraugli/butteraugli_main.cc @@ -0,0 +1,457 @@ +#include +#include +#include +#include +#include "butteraugli/butteraugli.h" + +extern "C" { +#include "png.h" +#include "jpeglib.h" +} + +namespace butteraugli { +namespace { + +// "rgb": cleared and filled with same-sized image planes (one per channel); +// either RGB, or RGBA if the PNG contains an alpha channel. +bool ReadPNG(FILE* f, std::vector* rgb) { + png_structp png_ptr = + png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); + if (!png_ptr) { + return false; + } + + png_infop info_ptr = png_create_info_struct(png_ptr); + if (!info_ptr) { + png_destroy_read_struct(&png_ptr, NULL, NULL); + return false; + } + + if (setjmp(png_jmpbuf(png_ptr)) != 0) { + // Ok we are here because of the setjmp. + png_destroy_read_struct(&png_ptr, &info_ptr, NULL); + return false; + } + + rewind(f); + png_init_io(png_ptr, f); + + // The png_transforms flags are as follows: + // packing == convert 1,2,4 bit images, + // strip == 16 -> 8 bits / channel, + // shift == use sBIT dynamics, and + // expand == palettes -> rgb, grayscale -> 8 bit images, tRNS -> alpha. + const unsigned int png_transforms = + PNG_TRANSFORM_PACKING | PNG_TRANSFORM_EXPAND | PNG_TRANSFORM_STRIP_16; + + png_read_png(png_ptr, info_ptr, png_transforms, NULL); + + png_bytep* row_pointers = png_get_rows(png_ptr, info_ptr); + + const int xsize = png_get_image_width(png_ptr, info_ptr); + const int ysize = png_get_image_height(png_ptr, info_ptr); + const int components = png_get_channels(png_ptr, info_ptr); + + *rgb = CreatePlanes(xsize, ysize, 3); + + switch (components) { + case 1: { + // GRAYSCALE + for (int y = 0; y < ysize; ++y) { + const uint8_t* const BUTTERAUGLI_RESTRICT row = row_pointers[y]; + uint8_t* const BUTTERAUGLI_RESTRICT row0 = (*rgb)[0].Row(y); + uint8_t* const BUTTERAUGLI_RESTRICT row1 = (*rgb)[1].Row(y); + uint8_t* const BUTTERAUGLI_RESTRICT row2 = (*rgb)[2].Row(y); + + for (int x = 0; x < xsize; ++x) { + const uint8_t gray = row[x]; + row0[x] = row1[x] = row2[x] = gray; + } + } + break; + } + case 2: { + // GRAYSCALE_ALPHA + rgb->push_back(Image8(xsize, ysize)); + for (int y = 0; y < ysize; ++y) { + const uint8_t* const BUTTERAUGLI_RESTRICT row = row_pointers[y]; + uint8_t* const BUTTERAUGLI_RESTRICT row0 = (*rgb)[0].Row(y); + uint8_t* const BUTTERAUGLI_RESTRICT row1 = (*rgb)[1].Row(y); + uint8_t* const BUTTERAUGLI_RESTRICT row2 = (*rgb)[2].Row(y); + uint8_t* const BUTTERAUGLI_RESTRICT row3 = (*rgb)[3].Row(y); + + for (int x = 0; x < xsize; ++x) { + const uint8_t gray = row[2 * x + 0]; + const uint8_t alpha = row[2 * x + 1]; + row0[x] = gray; + row1[x] = gray; + row2[x] = gray; + row3[x] = alpha; + } + } + break; + } + case 3: { + // RGB + for (int y = 0; y < ysize; ++y) { + const uint8_t* const BUTTERAUGLI_RESTRICT row = row_pointers[y]; + uint8_t* const BUTTERAUGLI_RESTRICT row0 = (*rgb)[0].Row(y); + uint8_t* const BUTTERAUGLI_RESTRICT row1 = (*rgb)[1].Row(y); + uint8_t* const BUTTERAUGLI_RESTRICT row2 = (*rgb)[2].Row(y); + + for (int x = 0; x < xsize; ++x) { + row0[x] = row[3 * x + 0]; + row1[x] = row[3 * x + 1]; + row2[x] = row[3 * x + 2]; + } + } + break; + } + case 4: { + // RGBA + rgb->push_back(Image8(xsize, ysize)); + for (int y = 0; y < ysize; ++y) { + const uint8_t* const BUTTERAUGLI_RESTRICT row = row_pointers[y]; + uint8_t* const BUTTERAUGLI_RESTRICT row0 = (*rgb)[0].Row(y); + uint8_t* const BUTTERAUGLI_RESTRICT row1 = (*rgb)[1].Row(y); + uint8_t* const BUTTERAUGLI_RESTRICT row2 = (*rgb)[2].Row(y); + uint8_t* const BUTTERAUGLI_RESTRICT row3 = (*rgb)[3].Row(y); + + for (int x = 0; x < xsize; ++x) { + row0[x] = row[4 * x + 0]; + row1[x] = row[4 * x + 1]; + row2[x] = row[4 * x + 2]; + row3[x] = row[4 * x + 3]; + } + } + break; + } + default: + png_destroy_read_struct(&png_ptr, &info_ptr, NULL); + return false; + } + png_destroy_read_struct(&png_ptr, &info_ptr, NULL); + return true; +} + +const double* NewSrgbToLinearTable() { + double* table = new double[256]; + for (int i = 0; i < 256; ++i) { + const double srgb = i / 255.0; + table[i] = + 255.0 * (srgb <= 0.04045 ? srgb / 12.92 + : std::pow((srgb + 0.055) / 1.055, 2.4)); + } + return table; +} + +void jpeg_catch_error(j_common_ptr cinfo) { + (*cinfo->err->output_message) (cinfo); + jmp_buf* jpeg_jmpbuf = (jmp_buf*) cinfo->client_data; + jpeg_destroy(cinfo); + longjmp(*jpeg_jmpbuf, 1); +} + +// "rgb": cleared and filled with same-sized image planes (one per channel); +// either RGB, or RGBA if the PNG contains an alpha channel. +bool ReadJPEG(FILE* f, std::vector* rgb) { + rewind(f); + + struct jpeg_decompress_struct cinfo; + struct jpeg_error_mgr jerr; + cinfo.err = jpeg_std_error(&jerr); + jmp_buf jpeg_jmpbuf; + cinfo.client_data = &jpeg_jmpbuf; + jerr.error_exit = jpeg_catch_error; + if (setjmp(jpeg_jmpbuf)) { + return false; + } + + jpeg_create_decompress(&cinfo); + + jpeg_stdio_src(&cinfo, f); + jpeg_read_header(&cinfo, TRUE); + jpeg_start_decompress(&cinfo); + + int row_stride = cinfo.output_width * cinfo.output_components; + JSAMPARRAY buffer = (*cinfo.mem->alloc_sarray) + ((j_common_ptr) &cinfo, JPOOL_IMAGE, row_stride, 1); + + const size_t xsize = cinfo.output_width; + const size_t ysize = cinfo.output_height; + + *rgb = CreatePlanes(xsize, ysize, 3); + + switch (cinfo.out_color_space) { + case JCS_GRAYSCALE: + while (cinfo.output_scanline < cinfo.output_height) { + jpeg_read_scanlines(&cinfo, buffer, 1); + + const uint8_t* const BUTTERAUGLI_RESTRICT row = buffer[0]; + uint8_t* const BUTTERAUGLI_RESTRICT row0 = + (*rgb)[0].Row(cinfo.output_scanline - 1); + uint8_t* const BUTTERAUGLI_RESTRICT row1 = + (*rgb)[1].Row(cinfo.output_scanline - 1); + uint8_t* const BUTTERAUGLI_RESTRICT row2 = + (*rgb)[2].Row(cinfo.output_scanline - 1); + + for (int x = 0; x < xsize; x++) { + const uint8_t gray = row[x]; + row0[x] = row1[x] = row2[x] = gray; + } + } + break; + + case JCS_RGB: + while (cinfo.output_scanline < cinfo.output_height) { + jpeg_read_scanlines(&cinfo, buffer, 1); + + const uint8_t* const BUTTERAUGLI_RESTRICT row = buffer[0]; + uint8_t* const BUTTERAUGLI_RESTRICT row0 = + (*rgb)[0].Row(cinfo.output_scanline - 1); + uint8_t* const BUTTERAUGLI_RESTRICT row1 = + (*rgb)[1].Row(cinfo.output_scanline - 1); + uint8_t* const BUTTERAUGLI_RESTRICT row2 = + (*rgb)[2].Row(cinfo.output_scanline - 1); + + for (int x = 0; x < xsize; x++) { + row0[x] = row[3 * x + 0]; + row1[x] = row[3 * x + 1]; + row2[x] = row[3 * x + 2]; + } + } + break; + + default: + jpeg_destroy_decompress(&cinfo); + return false; + } + + jpeg_finish_decompress(&cinfo); + jpeg_destroy_decompress(&cinfo); + return true; +} + +// Translate R, G, B channels from sRGB to linear space. If an alpha channel +// is present, overlay the image over a black or white background. Overlaying +// is done in the sRGB space; while technically incorrect, this is aligned with +// many other software (web browsers, WebP near lossless). +void FromSrgbToLinear(const std::vector& rgb, + std::vector& linear, int background) { + const size_t xsize = rgb[0].xsize(); + const size_t ysize = rgb[0].ysize(); + static const double* const kSrgbToLinearTable = NewSrgbToLinearTable(); + + if (rgb.size() == 3) { // RGB + for (int c = 0; c < 3; c++) { + linear.push_back(ImageF(xsize, ysize)); + for (int y = 0; y < ysize; ++y) { + const uint8_t* const BUTTERAUGLI_RESTRICT row_rgb = rgb[c].Row(y); + float* const BUTTERAUGLI_RESTRICT row_linear = linear[c].Row(y); + for (size_t x = 0; x < xsize; x++) { + const int value = row_rgb[x]; + row_linear[x] = kSrgbToLinearTable[value]; + } + } + } + } else { // RGBA + for (int c = 0; c < 3; c++) { + linear.push_back(ImageF(xsize, ysize)); + for (int y = 0; y < ysize; ++y) { + const uint8_t* const BUTTERAUGLI_RESTRICT row_rgb = rgb[c].Row(y); + float* const BUTTERAUGLI_RESTRICT row_linear = linear[c].Row(y); + const uint8_t* const BUTTERAUGLI_RESTRICT row_alpha = rgb[3].Row(y); + for (size_t x = 0; x < xsize; x++) { + int value; + if (row_alpha[x] == 255) { + value = row_rgb[x]; + } else if (row_alpha[x] == 0) { + value = background; + } else { + const int fg_weight = row_alpha[x]; + const int bg_weight = 255 - fg_weight; + value = + (row_rgb[x] * fg_weight + background * bg_weight + 127) / 255; + } + row_linear[x] = kSrgbToLinearTable[value]; + } + } + } + } +} + +std::vector ReadImageOrDie(const char* filename) { + std::vector rgb; + FILE* f = fopen(filename, "rb"); + if (!f) { + fprintf(stderr, "Cannot open %s\n", filename); + exit(1); + } + unsigned char magic[2]; + if (fread(magic, 1, 2, f) != 2) { + fprintf(stderr, "Cannot read from %s\n", filename); + exit(1); + } + if (magic[0] == 0xFF && magic[1] == 0xD8) { + if (!ReadJPEG(f, &rgb)) { + fprintf(stderr, "File %s is a malformed JPEG.\n", filename); + exit(1); + } + } else { + if (!ReadPNG(f, &rgb)) { + fprintf(stderr, "File %s is neither a valid JPEG nor a valid PNG.\n", + filename); + exit(1); + } + } + fclose(f); + return rgb; +} + +static void ScoreToRgb(double score, double good_threshold, + double bad_threshold, uint8_t rgb[3]) { + double heatmap[12][3] = { + { 0, 0, 0 }, + { 0, 0, 1 }, + { 0, 1, 1 }, + { 0, 1, 0 }, // Good level + { 1, 1, 0 }, + { 1, 0, 0 }, // Bad level + { 1, 0, 1 }, + { 0.5, 0.5, 1.0 }, + { 1.0, 0.5, 0.5 }, // Pastel colors for the very bad quality range. + { 1.0, 1.0, 0.5 }, + { 1, 1, 1, }, + { 1, 1, 1, }, + }; + if (score < good_threshold) { + score = (score / good_threshold) * 0.3; + } else if (score < bad_threshold) { + score = 0.3 + (score - good_threshold) / + (bad_threshold - good_threshold) * 0.15; + } else { + score = 0.45 + (score - bad_threshold) / + (bad_threshold * 12) * 0.5; + } + static const int kTableSize = sizeof(heatmap) / sizeof(heatmap[0]); + score = std::min(std::max( + score * (kTableSize - 1), 0.0), kTableSize - 2); + int ix = static_cast(score); + double mix = score - ix; + for (int i = 0; i < 3; ++i) { + double v = mix * heatmap[ix + 1][i] + (1 - mix) * heatmap[ix][i]; + rgb[i] = static_cast(255 * pow(v, 0.5) + 0.5); + } +} + +void CreateHeatMapImage(const ImageF& distmap, double good_threshold, + double bad_threshold, size_t xsize, size_t ysize, + std::vector* heatmap) { + heatmap->resize(3 * xsize * ysize); + for (size_t y = 0; y < ysize; ++y) { + for (size_t x = 0; x < xsize; ++x) { + int px = xsize * y + x; + double d = distmap.Row(y)[x]; + uint8_t* rgb = &(*heatmap)[3 * px]; + ScoreToRgb(d, good_threshold, bad_threshold, rgb); + } + } +} + +// main() function, within butteraugli namespace for convenience. +int Run(int argc, char* argv[]) { + if (argc != 3 && argc != 4) { + fprintf(stderr, + "Usage: %s {image1.(png|jpg|jpeg)} {image2.(png|jpg|jpeg)} " + "[heatmap.ppm]\n", + argv[0]); + return 1; + } + + std::vector rgb1 = ReadImageOrDie(argv[1]); + std::vector rgb2 = ReadImageOrDie(argv[2]); + + if (rgb1.size() != rgb2.size()) { + fprintf(stderr, "Different number of channels: %lu vs %lu\n", rgb1.size(), + rgb2.size()); + exit(1); + } + + for (size_t c = 0; c < rgb1.size(); ++c) { + if (rgb1[c].xsize() != rgb2[c].xsize() || + rgb1[c].ysize() != rgb2[c].ysize()) { + fprintf( + stderr, "The images are not equal in size: (%lu,%lu) vs (%lu,%lu)\n", + rgb1[c].xsize(), rgb2[c].xsize(), rgb1[c].ysize(), rgb2[c].ysize()); + return 1; + } + } + + // TODO: Figure out if it is a good idea to fetch the gamma from the image + // instead of applying sRGB conversion. + std::vector linear1, linear2; + // Overlay the image over a black background. + FromSrgbToLinear(rgb1, linear1, 0); + FromSrgbToLinear(rgb2, linear2, 0); + ImageF diff_map, diff_map_on_white; + double diff_value; + if (!butteraugli::ButteraugliInterface(linear1, linear2, diff_map, + diff_value)) { + fprintf(stderr, "Butteraugli comparison failed\n"); + return 1; + } + ImageF* diff_map_ptr = &diff_map; + if (rgb1.size() == 4 || rgb2.size() == 4) { + // If the alpha channel is present, overlay the image over a white + // background as well. + FromSrgbToLinear(rgb1, linear1, 255); + FromSrgbToLinear(rgb2, linear2, 255); + double diff_value_on_white; + if (!butteraugli::ButteraugliInterface(linear1, linear2, diff_map_on_white, + diff_value_on_white)) { + fprintf(stderr, "Butteraugli comparison failed\n"); + return 1; + } + if (diff_value_on_white > diff_value) { + diff_value = diff_value_on_white; + diff_map_ptr = &diff_map_on_white; + } + } + printf("%lf\n", diff_value); + + if (argc == 4) { + const double good_quality = ::butteraugli::ButteraugliFuzzyInverse(1.5); + const double bad_quality = ::butteraugli::ButteraugliFuzzyInverse(0.5); + std::vector rgb; + CreateHeatMapImage(*diff_map_ptr, good_quality, bad_quality, + rgb1[0].xsize(), rgb2[0].ysize(), &rgb); + FILE* const fmap = fopen(argv[3], "wb"); + if (fmap == NULL) { + fprintf(stderr, "Cannot open %s\n", argv[3]); + perror("fopen"); + return 1; + } + bool ok = true; + if (fprintf(fmap, "P6\n%lu %lu\n255\n", + rgb1[0].xsize(), rgb1[0].ysize()) < 0){ + perror("fprintf"); + ok = false; + } + if (ok && fwrite(rgb.data(), 1, rgb.size(), fmap) != rgb.size()) { + perror("fwrite"); + ok = false; + } + if (fclose(fmap) != 0) { + perror("fclose"); + ok = false; + } + if (!ok) return 1; + } + + return 0; +} + +} // namespace +} // namespace butteraugli + +int main(int argc, char** argv) { return butteraugli::Run(argc, argv); } diff --git a/oss-internship-2020/guetzli/third_party/butteraugli/jpeg.BUILD b/oss-internship-2020/guetzli/third_party/butteraugli/jpeg.BUILD new file mode 100755 index 0000000..92c9ddc --- /dev/null +++ b/oss-internship-2020/guetzli/third_party/butteraugli/jpeg.BUILD @@ -0,0 +1,89 @@ +# Description: +# The Independent JPEG Group's JPEG runtime library. + +licenses(["notice"]) # custom notice-style license, see LICENSE + +cc_library( + name = "jpeg", + srcs = [ + "cderror.h", + "cdjpeg.h", + "jaricom.c", + "jcapimin.c", + "jcapistd.c", + "jcarith.c", + "jccoefct.c", + "jccolor.c", + "jcdctmgr.c", + "jchuff.c", + "jcinit.c", + "jcmainct.c", + "jcmarker.c", + "jcmaster.c", + "jcomapi.c", + "jconfig.h", + "jcparam.c", + "jcprepct.c", + "jcsample.c", + "jctrans.c", + "jdapimin.c", + "jdapistd.c", + "jdarith.c", + "jdatadst.c", + "jdatasrc.c", + "jdcoefct.c", + "jdcolor.c", + "jdct.h", + "jddctmgr.c", + "jdhuff.c", + "jdinput.c", + "jdmainct.c", + "jdmarker.c", + "jdmaster.c", + "jdmerge.c", + "jdpostct.c", + "jdsample.c", + "jdtrans.c", + "jerror.c", + "jfdctflt.c", + "jfdctfst.c", + "jfdctint.c", + "jidctflt.c", + "jidctfst.c", + "jidctint.c", + "jinclude.h", + "jmemmgr.c", + "jmemnobs.c", + "jmemsys.h", + "jmorecfg.h", + "jquant1.c", + "jquant2.c", + "jutils.c", + "jversion.h", + "transupp.h", + ], + hdrs = [ + "jerror.h", + "jpegint.h", + "jpeglib.h", + ], + includes = ["."], + visibility = ["//visibility:public"], +) + +genrule( + name = "configure", + outs = ["jconfig.h"], + cmd = "cat <$@\n" + + "#define HAVE_PROTOTYPES 1\n" + + "#define HAVE_UNSIGNED_CHAR 1\n" + + "#define HAVE_UNSIGNED_SHORT 1\n" + + "#define HAVE_STDDEF_H 1\n" + + "#define HAVE_STDLIB_H 1\n" + + "#ifdef WIN32\n" + + "#define INLINE __inline\n" + + "#else\n" + + "#define INLINE __inline__\n" + + "#endif\n" + + "EOF\n", +) diff --git a/oss-internship-2020/guetzli/third_party/butteraugli/png.BUILD b/oss-internship-2020/guetzli/third_party/butteraugli/png.BUILD new file mode 100755 index 0000000..9ff982b --- /dev/null +++ b/oss-internship-2020/guetzli/third_party/butteraugli/png.BUILD @@ -0,0 +1,33 @@ +# Description: +# libpng is the official PNG reference library. + +licenses(["notice"]) # BSD/MIT-like license + +cc_library( + name = "png", + srcs = [ + "png.c", + "pngerror.c", + "pngget.c", + "pngmem.c", + "pngpread.c", + "pngread.c", + "pngrio.c", + "pngrtran.c", + "pngrutil.c", + "pngset.c", + "pngtrans.c", + "pngwio.c", + "pngwrite.c", + "pngwtran.c", + "pngwutil.c", + ], + hdrs = [ + "png.h", + "pngconf.h", + ], + includes = ["."], + linkopts = ["-lm"], + visibility = ["//visibility:public"], + deps = ["@zlib_archive//:zlib"], +) diff --git a/oss-internship-2020/guetzli/third_party/butteraugli/zlib.BUILD b/oss-internship-2020/guetzli/third_party/butteraugli/zlib.BUILD new file mode 100755 index 0000000..edb77fd --- /dev/null +++ b/oss-internship-2020/guetzli/third_party/butteraugli/zlib.BUILD @@ -0,0 +1,36 @@ +package(default_visibility = ["//visibility:public"]) + +licenses(["notice"]) # BSD/MIT-like license (for zlib) + +cc_library( + name = "zlib", + srcs = [ + "adler32.c", + "compress.c", + "crc32.c", + "crc32.h", + "deflate.c", + "deflate.h", + "gzclose.c", + "gzguts.h", + "gzlib.c", + "gzread.c", + "gzwrite.c", + "infback.c", + "inffast.c", + "inffast.h", + "inffixed.h", + "inflate.c", + "inflate.h", + "inftrees.c", + "inftrees.h", + "trees.c", + "trees.h", + "uncompr.c", + "zconf.h", + "zutil.c", + "zutil.h", + ], + hdrs = ["zlib.h"], + includes = ["."], +)