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 <algorithm>
 #include <iostream>
+#include <string>
 #include <stdio.h>
 #include <sys/stat.h>
 #include <vector>
 
 namespace {
 
-inline uint8_t BlendOnBlack(const uint8_t val, const uint8_t alpha) {
-  return (static_cast<int>(val) * static_cast<int>(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<typename T>
 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<const char*>(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<uint8_t> in_data_temp;
-  in_data_temp.reserve(rgb->size);
-
-  auto* rgb_data = static_cast<uint8_t*>(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<const char*>(in_data->data), in_data->size);
-  std::vector<uint8_t> rgb;
+sapi::StatusOr<std::string> 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<char[]> 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<GuetzliInitData> 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<float>(
+    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<int>(val) * static_cast<int>(alpha) + 128) / 255;
+}
+
+bool ReadPNG(const std::string& data, int* xsize, int* ysize,
+             std::vector<uint8_t>* 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<const char*>(in_data->data), in_data->size);
-  guetzli::JPEGData jpg;
-
-  auto result = guetzli::ReadJpeg(data, 
-    static_cast<guetzli::JpegReadMode>(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<double>(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<char[]> buf(new char[fsize]);
-  status = read(fd, buf.get(), fsize);
+  int xsize, ysize;
+  std::vector<uint8_t> 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<std::pair<std::string, int>> 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::string, int>{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 <iostream>
+#include <memory>
 
 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<ImageType> 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>* 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<double>(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>* 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<double>(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> 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<ProcessingParams> 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, &params, &input, &output);
-
-    if (!process_status.ok()) {
-      return process_status;
-    }
-  } else {
-    auto process_status = ProcessJpeg(&api, &params, &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<ImageType> 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>* params, 
-                          sapi::v::LenVal* input, 
-                          sapi::v::LenVal* output) const;
-
-  absl::Status ProcessJpeg(GuetzliApi* api, 
-                          sapi::v::Struct<Params>* params, 
-                          sapi::v::LenVal* input, 
-                          sapi::v::LenVal* output) const;
+  sapi::StatusOr<ImageType> 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>* 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<double>(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>* 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<double>(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> 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, &params, &input, &output);
-
-    if (!process_status.ok()) {
-      return process_status;
-    }
-  } else {
-    auto process_status = ProcessJpeg(&api, &params, &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 <assert.h>
+#include <math.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <string.h>
+
+#include <algorithm>
+#include <array>
+
+
+// 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<char *>(malloc(bytes + kCacheLineSize));
+  if (allocated == nullptr) {
+    return nullptr;
+  }
+  const uintptr_t misalignment =
+      reinterpret_cast<uintptr_t>(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<char *>(aligned_pointer);
+  assert(reinterpret_cast<uintptr_t>(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<float> 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<int>(1, m * fabs(sigma));
+  std::vector<float> 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<float>& 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<int>(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<float>& 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<float> 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<float> 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<int>(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<ImageF> OpsinDynamicsImage(const std::vector<ImageF>& rgb) {
+  PROFILER_FUNC;
+  std::vector<ImageF> xyb(3);
+  std::vector<ImageF> 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 <class V>
+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<ImageF>& 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<ImageF>* BUTTERAUGLI_RESTRICT mask,
+                     std::vector<ImageF>* BUTTERAUGLI_RESTRICT mask_dc) {
+  std::vector<ImageF> mask_xyb0 = CreatePlanes<float>(xsize, ysize, 3);
+  std::vector<ImageF> mask_xyb1 = CreatePlanes<float>(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<ImageF>& rgb0)
+    : xsize_(rgb0[0].xsize()),
+      ysize_(rgb0[0].ysize()),
+      num_pixels_(xsize_ * ysize_) {
+  if (xsize_ < 8 || ysize_ < 8) return;
+  std::vector<ImageF> xyb0 = OpsinDynamicsImage(rgb0);
+  SeparateFrequencies(xsize_, ysize_, xyb0, pi0_);
+}
+
+void ButteraugliComparator::Mask(
+    std::vector<ImageF>* BUTTERAUGLI_RESTRICT mask,
+    std::vector<ImageF>* BUTTERAUGLI_RESTRICT mask_dc) const {
+  MaskPsychoImage(pi0_, pi0_, xsize_, ysize_, mask, mask_dc);
+}
+
+void ButteraugliComparator::Diffmap(const std::vector<ImageF>& rgb1,
+                                    ImageF &result) const {
+  PROFILER_FUNC;
+  if (xsize_ < 8 || ysize_ < 8) return;
+  DiffmapOpsinDynamicsImage(OpsinDynamicsImage(rgb1), result);
+}
+
+void ButteraugliComparator::DiffmapOpsinDynamicsImage(
+    const std::vector<ImageF>& 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<ImageF> block_diff_dc(3);
+  std::vector<ImageF> 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<ImageF> mask_xyb;
+  std::vector<ImageF> 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 <bool fastMode, class Tag>
+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 <class Tag>
+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<float> 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<float>(norm1) + absval);
+
+      // Primary symmetric quadratic objective.
+      diffs[ix] = scaler * diff;
+
+      const float scaler2 = norm2_0lt1 / (static_cast<float>(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<false, Tag>(&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<false, Tag>(&diffs[0], x0, y0, xsize_, ysize_);
+    }
+    for (; x0 < xsize_ - 4; ++x0) {
+      row_diff[x0] +=
+          PaddedMaltaUnit<true, Tag>(&diffs[0], x0, y0, xsize_, ysize_);
+    }
+
+    for (; x0 < xsize_; ++x0) {
+      row_diff[x0] +=
+          PaddedMaltaUnit<false, Tag>(&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<false, Tag>(&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<MaltaTag>(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<MaltaTagLF>(lum0, lum1, xsize_, ysize_,
+                               w_0gt1, w_0lt1,
+                               norm1, len,
+                               mulli, block_diff_ac);
+}
+
+ImageF ButteraugliComparator::CombineChannels(
+    const std::vector<ImageF>& mask_xyb,
+    const std::vector<ImageF>& mask_xyb_dc,
+    const std::vector<ImageF>& block_diff_dc,
+    const std::vector<ImageF>& 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 <stdio.h>
+
+// ===== Functions used by Mask only =====
+static std::array<double, 512> MakeMask(
+    double extmul, double extoff,
+    double mul, double offset,
+    double scaler) {
+  std::array<double, 512> 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<double, 512> 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<double, 512> 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<double, 512> 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<double, 512> 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<ImageF>& xyb0,
+          const std::vector<ImageF>& xyb1,
+          std::vector<ImageF>* BUTTERAUGLI_RESTRICT mask,
+          std::vector<ImageF>* 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<float>(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<double>(
+      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<ImageF> &rgb0_image,
+                        const std::vector<ImageF> &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<size_t>(kMax, xsize);
+    size_t yscaled = std::max<size_t>(kMax, ysize);
+    std::vector<ImageF> scaled0 = CreatePlanes<float>(xscaled, yscaled, 3);
+    std::vector<ImageF> scaled1 = CreatePlanes<float>(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<size_t>(xsize - 1, std::max(0, x - xborder));
+          size_t y2 = std::min<size_t>(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<ImageF> &rgb0,
+                          const std::vector<ImageF> &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<std::vector<float> > &rgb, std::vector<float> &quant) {
+  if (xsize < 16 || ysize < 16) {
+    return false;  // Butteraugli is undefined for small images.
+  }
+  size_t size = xsize * ysize;
+
+  std::vector<ImageF> rgb_planes = PlanesFromPacked(xsize, ysize, rgb);
+  std::vector<ImageF> scale_xyb;
+  std::vector<ImageF> 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<double>(std::max<double>(score * (kTableSize - 1), 0.0),
+                           kTableSize - 2);
+  int ix = static_cast<int>(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<uint8_t>(255 * pow(v, 0.5) + 0.5);
+  }
+}
+
+}  // namespace
+
+void CreateHeatMapImage(const std::vector<float>& distmap,
+                        double good_threshold, double bad_threshold,
+                        size_t xsize, size_t ysize,
+                        std::vector<uint8_t>* 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 <cassert>
+#include <cmath>
+#include <cstdint>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <memory>
+#include <vector>
+
+#define BUTTERAUGLI_ENABLE_CHECKS 0
+
+// This is the main interface to butteraugli image similarity
+// analysis function.
+
+namespace butteraugli {
+
+template<typename T>
+class Image;
+
+using Image8 = Image<uint8_t>;
+using ImageF = Image<float>;
+
+// 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<ImageF> &rgb0,
+                          const std::vector<ImageF> &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<std::vector<float> > &rgb, std::vector<float> &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 <typename T>
+using CacheAlignedUniquePtrT = std::unique_ptr<T[], void (*)(void *)>;
+
+using CacheAlignedUniquePtr = CacheAlignedUniquePtrT<uint8_t>;
+
+template <typename T = uint8_t>
+static inline CacheAlignedUniquePtrT<T> Allocate(const size_t entries) {
+  return CacheAlignedUniquePtrT<T>(
+      static_cast<T * const BUTTERAUGLI_RESTRICT>(
+          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 <size_t multiple>
+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 <typename ComponentType>
+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<uint8_t*>(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<T *>(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<const T *>(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<intptr_t>(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 <typename T>
+static inline std::vector<Image<T>> CreatePlanes(const size_t xsize,
+                                                 const size_t ysize,
+                                                 const size_t num_planes) {
+  std::vector<Image<T>> 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 <typename T>
+static inline Image<T> CopyPixels(const Image<T> &other) {
+  Image<T> 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 <typename T>
+static inline std::vector<Image<T>> CopyPlanes(
+    const std::vector<Image<T>> &planes) {
+  std::vector<Image<T>> copy;
+  copy.reserve(planes.size());
+  for (const Image<T> &plane : planes) {
+    copy.push_back(CopyPixels(plane));
+  }
+  return copy;
+}
+
+// Compacts a padded image into a preallocated packed vector.
+template <typename T>
+static inline void CopyToPacked(const Image<T> &from, std::vector<T> *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 <typename T>
+static inline void CopyFromPacked(const std::vector<T> &from, Image<T> *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 <typename T>
+static inline std::vector<Image<T>> PlanesFromPacked(
+    const size_t xsize, const size_t ysize,
+    const std::vector<std::vector<T>> &packed) {
+  std::vector<Image<T>> planes;
+  planes.reserve(packed.size());
+  for (const std::vector<T> &p : packed) {
+    planes.push_back(Image<T>(xsize, ysize));
+    CopyFromPacked(p, &planes.back());
+  }
+  return planes;
+}
+
+template <typename T>
+static inline std::vector<std::vector<T>> PackedFromPlanes(
+    const std::vector<Image<T>> &planes) {
+  assert(!planes.empty());
+  const size_t num_pixels = planes[0].xsize() * planes[0].ysize();
+  std::vector<std::vector<T>> packed;
+  packed.reserve(planes.size());
+  for (const Image<T> &image : planes) {
+    packed.push_back(std::vector<T>(num_pixels));
+    CopyToPacked(image, &packed.back());
+  }
+  return packed;
+}
+
+struct PsychoImage {
+  std::vector<ImageF> uhf;
+  std::vector<ImageF> hf;
+  std::vector<ImageF> mf;
+  std::vector<ImageF> lf;
+};
+
+class ButteraugliComparator {
+ public:
+  ButteraugliComparator(const std::vector<ImageF>& rgb0);
+
+  // Computes the butteraugli map between the original image given in the
+  // constructor and the distorted image give here.
+  void Diffmap(const std::vector<ImageF>& rgb1, ImageF& result) const;
+
+  // Same as above, but OpsinDynamicsImage() was already applied.
+  void DiffmapOpsinDynamicsImage(const std::vector<ImageF>& 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<ImageF>* BUTTERAUGLI_RESTRICT mask,
+            std::vector<ImageF>* 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<ImageF>& scale_xyb,
+                         const std::vector<ImageF>& scale_xyb_dc,
+                         const std::vector<ImageF>& block_diff_dc,
+                         const std::vector<ImageF>& block_diff_ac) const;
+
+  const size_t xsize_;
+  const size_t ysize_;
+  const size_t num_pixels_;
+  PsychoImage pi0_;
+};
+
+void ButteraugliDiffmap(const std::vector<ImageF> &rgb0,
+                        const std::vector<ImageF> &rgb1,
+                        ImageF &diffmap);
+
+double ButteraugliScoreFromDiffmap(const ImageF& distmap);
+
+// Generate rgb-representation of the distance between two images.
+void CreateHeatMapImage(const std::vector<float> &distmap,
+                        double good_threshold, double bad_threshold,
+                        size_t xsize, size_t ysize,
+                        std::vector<uint8_t> *heatmap);
+
+// Compute values of local frequency and dc masking based on the activity
+// in the two images.
+void Mask(const std::vector<ImageF>& xyb0,
+          const std::vector<ImageF>& xyb1,
+          std::vector<ImageF>* BUTTERAUGLI_RESTRICT mask,
+          std::vector<ImageF>* BUTTERAUGLI_RESTRICT mask_dc);
+
+template <class V>
+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 <class V>
+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<ImageF> OpsinDynamicsImage(const std::vector<ImageF>& 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 <int INDEX>
+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<INDEX - 1>(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 <int N>
+  static double EvaluatePolynomial(const double x,
+                                   const double (&coefficients)[N]) {
+    double b1 = 0.0;
+    double b2 = 0.0;
+    ClenshawRecursion<N - 1>(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<float>(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 <cmath>
+#include <cstdint>
+#include <cstdio>
+#include <vector>
+#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<Image8>* 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<uint8_t>(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<Image8>* 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<uint8_t>(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<Image8>& rgb,
+                      std::vector<ImageF>& 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<Image8> ReadImageOrDie(const char* filename) {
+  std::vector<Image8> 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<double>(std::max<double>(
+      score * (kTableSize - 1), 0.0), kTableSize - 2);
+  int ix = static_cast<int>(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<uint8_t>(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<uint8_t>* 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<Image8> rgb1 = ReadImageOrDie(argv[1]);
+  std::vector<Image8> 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<ImageF> 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<uint8_t> 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 <<EOF >$@\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 = ["."],
+)