diff --git a/Makefile b/Makefile
index 546a24c..cecfca7 100644
--- a/Makefile
+++ b/Makefile
@@ -3,10 +3,10 @@ all: piphoto
 objects = piphoto.o color.o coord.o util.o
 
 piphoto: $(objects) Makefile
-	clang-3.9 -O3 -g -Weverything -Werror -Wno-c++98-compat -Wno-c++98-c++11-compat-pedantic --std=c++1z --stdlib=libc++ -o piphoto $(objects) -lc++ -lunwind -lpng
+	clang-3.9 -O3 -g -Weverything -Werror --std=c++1z --stdlib=libc++ -o piphoto $(objects) -lc++ -lunwind -lpng
 
 %.o: %.cc *.h Makefile
-	clang-3.9 -O3 -g -Weverything -Werror -Wno-c++98-compat -Wno-c++98-c++11-compat-pedantic --std=c++1z --stdlib=libc++ -c -o $@ $<
+	clang-3.9 -O3 -g -Weverything -Werror -Wno-padded -Wno-c++98-compat -Wno-c++98-c++11-compat-pedantic --std=c++1z --stdlib=libc++ -c -o $@ $<
 
 run: piphoto
 	./piphoto
diff --git a/color.cc b/color.cc
index e47de3b..42f5dcb 100644
--- a/color.cc
+++ b/color.cc
@@ -6,8 +6,8 @@ std::ostream& operator<<(std::ostream& os, const Color& color) {
   return os
     << std::hex << std::setfill('0')
     << "rgb("
-    << "0x" << std::setw(4) << color.r << ", "
-    << "0x" << std::setw(4) << color.g << ", "
-    << "0x" << std::setw(4) << color.b
+    << "0x" << std::setw(4) << color.at(0) << ", "
+    << "0x" << std::setw(4) << color.at(1) << ", "
+    << "0x" << std::setw(4) << color.at(2)
     << ")" << std::dec;
 }
diff --git a/color.h b/color.h
index e085a67..b19a4f5 100644
--- a/color.h
+++ b/color.h
@@ -1,24 +1,23 @@
 #pragma once
 
+#include <array>
 #include <cstdint>
 #include <iostream>
 
-constexpr uint32_t kNumColors = (1 << 16);
+#include "intmath.h"
 
-struct Color {
-  // 32-bit for compiler convenience, but values are 16-bit
-  uint32_t r;
-  uint32_t g;
-  uint32_t b;
+constexpr uint32_t kNumColors = UINT16_MAX;
 
+// 32-bit for compiler convenience, but values are 16-bit
+struct Color : public std::array<uint32_t, 3> {
   constexpr uint32_t Difference(const Color& other) const;
 };
 
 constexpr uint32_t Color::Difference(const Color& other) const {
   return (
-    ((r > other.r) ? (r - other.r) : (other.r - r)) +
-    ((g > other.g) ? (g - other.g) : (other.g - g)) +
-    ((b > other.b) ? (b - other.b) : (other.b - b))
+    AbsDiff(this->at(0), other.at(0)) +
+    AbsDiff(this->at(1), other.at(1)) +
+    AbsDiff(this->at(2), other.at(2))
   );
 }
 
diff --git a/colorchecker.h b/colorchecker.h
index 133154b..aea22bc 100644
--- a/colorchecker.h
+++ b/colorchecker.h
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <array>
+#include <cstdint>
 #include <numeric>
 
 #include "color.h"
@@ -8,35 +9,36 @@
 #include "coord.h"
 #include "image.h"
 #include "lut.h"
+#include "minimum.h"
 
 // Maximum LUT size that has each point adjacent to at least one ColorChecker color.
 typedef Lut3d<4, 3, 3> ColorCheckerLut3d;
 
 constexpr std::array<Color, 24> kColorCheckerSrgb = {{
-  {0x7300, 0x5200, 0x4400},
-  {0xc200, 0x9600, 0x8200},
-  {0x6200, 0x7a00, 0x9d00},
-  {0x5700, 0x6c00, 0x4300},
-  {0x8500, 0x8000, 0xb100},
-  {0x6700, 0xbd00, 0xaa00},
-  {0xd600, 0x7e00, 0x2c00},
-  {0x5000, 0x5b00, 0xa600},
-  {0xc100, 0x5a00, 0x6300},
-  {0x5e00, 0x3c00, 0x6c00},
-  {0x9d00, 0xbc00, 0x4000},
-  {0xe000, 0xa300, 0x2e00},
-  {0x3800, 0x3d00, 0x9600},
-  {0x4600, 0x9400, 0x4900},
-  {0xaf00, 0x3600, 0x3c00},
-  {0xe700, 0xc700, 0x1f00},
-  {0xbb00, 0x5600, 0x9500},
-  {0x0800, 0x8500, 0xa100},
-  {0xf300, 0xf300, 0xf200},
-  {0xc800, 0xc800, 0xc800},
-  {0xa000, 0xa000, 0xa000},
-  {0x7a00, 0x7a00, 0x7900},
-  {0x5500, 0x5500, 0x5500},
-  {0x3400, 0x3400, 0x3400},
+  {{{0x7300, 0x5200, 0x4400}}},
+  {{{0xc200, 0x9600, 0x8200}}},
+  {{{0x6200, 0x7a00, 0x9d00}}},
+  {{{0x5700, 0x6c00, 0x4300}}},
+  {{{0x8500, 0x8000, 0xb100}}},
+  {{{0x6700, 0xbd00, 0xaa00}}},
+  {{{0xd600, 0x7e00, 0x2c00}}},
+  {{{0x5000, 0x5b00, 0xa600}}},
+  {{{0xc100, 0x5a00, 0x6300}}},
+  {{{0x5e00, 0x3c00, 0x6c00}}},
+  {{{0x9d00, 0xbc00, 0x4000}}},
+  {{{0xe000, 0xa300, 0x2e00}}},
+  {{{0x3800, 0x3d00, 0x9600}}},
+  {{{0x4600, 0x9400, 0x4900}}},
+  {{{0xaf00, 0x3600, 0x3c00}}},
+  {{{0xe700, 0xc700, 0x1f00}}},
+  {{{0xbb00, 0x5600, 0x9500}}},
+  {{{0x0800, 0x8500, 0xa100}}},
+  {{{0xf300, 0xf300, 0xf200}}},
+  {{{0xc800, 0xc800, 0xc800}}},
+  {{{0xa000, 0xa000, 0xa000}}},
+  {{{0x7a00, 0x7a00, 0x7900}}},
+  {{{0x5500, 0x5500, 0x5500}}},
+  {{{0x3400, 0x3400, 0x3400}}},
 }};
 
 template <uint32_t X, uint32_t Y>
@@ -103,3 +105,37 @@ std::unique_ptr<Image<X, Y>> HighlightClosest(const Image<X, Y>& image) {
   }
   return out;
 }
+
+template <uint32_t P, uint32_t R, uint32_t G, uint32_t B, uint32_t X, uint32_t Y>
+uint32_t OptimizeLut(const Image<X, Y>& image, Lut3d<R, G, B>* lut) {
+  uint32_t diff = 0;
+
+  for (uint32_t r = 0; r < R; ++r) {
+    auto& rect = lut->at(r);
+
+    for (uint32_t g = 0; g < G; ++g) {
+      auto& row = rect.at(g);
+
+      for (uint32_t b = 0; b < B; ++b) {
+        auto& color = row.at(b);
+
+        std::cout << "(" << r << ", " << g << ", " << b << ")" << std::endl;
+
+        for (uint32_t c = 0; c < color.size(); ++c) {
+          auto min = FindPossibleMinimum<uint32_t, uint32_t, 4>(
+            0, UINT16_MAX,
+            [&image, &lut, r, g, b, c](uint32_t val) {
+              auto test_lut = *lut;
+              test_lut[r][g][b].at(c) = val;
+              return ScoreImage(*test_lut.MapImage(image));
+            });
+          std::cout << "\tC" << c << ": " << color.at(c) << " -> " << min << std::endl;
+          diff += AbsDiff(color.at(c), min);
+          color.at(c) = min;
+        }
+      }
+    }
+  }
+
+  return diff;
+}
diff --git a/colors.h b/colors.h
index 0ac1671..85db66d 100644
--- a/colors.h
+++ b/colors.h
@@ -2,5 +2,5 @@
 
 #include "color.h"
 
-constexpr Color kBlack = Color{0x0000, 0x0000, 0x0000};
-constexpr Color kWhite = Color{0xffff, 0xffff, 0xffff};
+constexpr Color kBlack = {{{0x0000, 0x0000, 0x0000}}};
+constexpr Color kWhite = {{{0xffff, 0xffff, 0xffff}}};
diff --git a/image.h b/image.h
index 7e958fa..8e8de25 100644
--- a/image.h
+++ b/image.h
@@ -84,9 +84,9 @@ std::string Image<X, Y>::ToPng() {
   for (auto& row : *this) {
     std::array<uint16_t, X * 3> out_row;
     for (uint32_t x = 0; x < X; ++x) {
-      out_row[x * 3 + 0] = htons(static_cast<uint16_t>(row[x].r));
-      out_row[x * 3 + 1] = htons(static_cast<uint16_t>(row[x].g));
-      out_row[x * 3 + 2] = htons(static_cast<uint16_t>(row[x].b));
+      out_row[x * 3 + 0] = htons(static_cast<uint16_t>(row[x].at(0)));
+      out_row[x * 3 + 1] = htons(static_cast<uint16_t>(row[x].at(1)));
+      out_row[x * 3 + 2] = htons(static_cast<uint16_t>(row[x].at(2)));
     }
     png_write_row(png_ptr, reinterpret_cast<unsigned char*>(out_row.data()));
   }
diff --git a/intmath.h b/intmath.h
new file mode 100644
index 0000000..79cd229
--- /dev/null
+++ b/intmath.h
@@ -0,0 +1,6 @@
+#pragma once
+
+template <typename T>
+constexpr T AbsDiff(T a, T b) {
+  return (a > b) ? (a - b) : (b - a);
+}
diff --git a/lut.h b/lut.h
index 53e96b9..dd2f1f7 100644
--- a/lut.h
+++ b/lut.h
@@ -6,7 +6,7 @@
 template <uint32_t R, uint32_t G, uint32_t B>
 class Lut3d : public std::array<std::array<std::array<Color, B>, G>, R> {
  public:
-  static std::unique_ptr<Lut3d<R, G, B>> Identity();
+  static Lut3d<R, G, B> Identity();
 
   Color MapColor(const Color& in) const;
 
@@ -28,20 +28,20 @@ class Lut3d : public std::array<std::array<std::array<Color, B>, G>, R> {
 typedef Lut3d<2, 2, 2> MinimalLut3d;
 
 template <uint32_t R, uint32_t G, uint32_t B>
-std::unique_ptr<Lut3d<R, G, B>> Lut3d<R, G, B>::Identity() {
-  auto ret = std::make_unique<Lut3d<R, G, B>>();
+Lut3d<R, G, B> Lut3d<R, G, B>::Identity() {
+  Lut3d<R, G, B> ret;
 
   Color color;
   for (uint32_t r = 0; r < R; ++r) {
-    auto& rect = ret->at(r);
-    color.r = std::min(kNumColors - 1, BlockSize(R) * r);
+    auto& rect = ret.at(r);
+    color.at(0) = std::min(kNumColors - 1, BlockSize(R) * r);
 
     for (uint32_t g = 0; g < G; ++g) {
       auto& row = rect.at(g);
-      color.g = std::min(kNumColors - 1, BlockSize(G) * g);
+      color.at(1) = std::min(kNumColors - 1, BlockSize(G) * g);
 
       for (uint32_t b = 0; b < B; ++b) {
-        color.b = std::min(kNumColors - 1, BlockSize(B) * b);
+        color.at(2) = std::min(kNumColors - 1, BlockSize(B) * b);
         row.at(b) = color;
       }
     }
@@ -104,24 +104,27 @@ std::unique_ptr<Image<X, Y>> Lut3d<R, G, B>::MapImage(const Image<X, Y>& in) con
 
 template <uint32_t R, uint32_t G, uint32_t B>
 constexpr Color Lut3d<R, G, B>::InterpolateColor(const Color& i0, const Color& i1, uint32_t mul, uint32_t div) {
-  return {
-    Interpolate(i0.r, i1.r, mul, div),
-    Interpolate(i0.g, i1.g, mul, div),
-    Interpolate(i0.b, i1.b, mul, div),
-  };
+  return {{{
+    Interpolate(i0.at(0), i1.at(0), mul, div),
+    Interpolate(i0.at(1), i1.at(1), mul, div),
+    Interpolate(i0.at(2), i1.at(2), mul, div),
+  }}};
 }
 
 template <uint32_t R, uint32_t G, uint32_t B>
 constexpr uint32_t Lut3d<R, G, B>::Interpolate(uint32_t i0, uint32_t i1, uint32_t mul, uint32_t div) {
-  assert(i1 >= i0);
-  return i0 + ((mul * (i1 - i0)) / div);
+  if (i1 > i0) {
+    return i0 + ((mul * (i1 - i0)) / div);
+  } else {
+    return i1 + ((mul * (i0 - i1)) / div);
+  }
 }
 
 template <uint32_t R, uint32_t G, uint32_t B>
 constexpr std::pair<Coord3d, Coord3d> Lut3d<R, G, B>::FindRoot(const Color& in) {
-  auto root_r = FindChannelRoot(in.r, R);
-  auto root_g = FindChannelRoot(in.g, G);
-  auto root_b = FindChannelRoot(in.b, B);
+  auto root_r = FindChannelRoot(in.at(0), R);
+  auto root_g = FindChannelRoot(in.at(1), G);
+  auto root_b = FindChannelRoot(in.at(2), B);
   return {
     {root_r.first, root_g.first, root_b.first},
     {root_r.second, root_g.second, root_b.second},
diff --git a/minimum.h b/minimum.h
new file mode 100644
index 0000000..c1fe181
--- /dev/null
+++ b/minimum.h
@@ -0,0 +1,49 @@
+#pragma once
+
+template <typename I, typename O>
+struct Range {
+  I start;
+  I end;
+  I testpoint;
+  O testpoint_value;
+};
+
+// Find the minimum value of a callback within a range, using a given
+// parallelism.
+//
+// Deterministic for a given parallelism, but not guaranteed to be correct.
+// Since it does a non-exhaustive search, can be fooled by distributions with
+// multiple peaks, especially those with the minimum in a narrow valley and
+// other wider valleys.
+template <typename I, typename O, uint32_t P>
+I FindPossibleMinimum(I min, I max, std::function<O(I)> callback) {
+  if (min == max) {
+    return min;
+  }
+
+  std::array<Range<I, O>, P> ranges;
+
+  const I step = ((max - min) / P) + 1;
+  const I offset = step / 2;
+  for (uint32_t i = 0; i < P; ++i) {
+    auto& range = ranges[i];
+    range.start = std::min(max, min + i * step);
+    range.end = std::min(max, range.start + (step - 1));
+    range.testpoint = range.start + offset;
+  }
+
+  // TODO: threads
+  for (auto& range : ranges) {
+    range.testpoint_value = callback(range.testpoint);
+  }
+
+  const auto& min_range = *std::min_element(ranges.begin(), ranges.end(), [](const Range<I, O>& a, const Range<I, O>& b) {
+    return a.testpoint_value < b.testpoint_value;
+  });
+
+  if (step == 1) {
+    return min_range.testpoint;
+  } else {
+    return FindPossibleMinimum<I, O, P>(min_range.start, min_range.end, callback);
+  }
+}
diff --git a/piphoto.cc b/piphoto.cc
index a2e2652..88bec05 100644
--- a/piphoto.cc
+++ b/piphoto.cc
@@ -7,8 +7,15 @@
 
 int main() {
   auto image = PiRaw2::FromJpeg(ReadFile("test.jpg"));
+  std::cout << "Initial score: " << ScoreImage(*image) << std::endl;
+
   auto lut = MinimalLut3d::Identity();
-  auto image2 = lut->MapImage(*image);
-  std::cout << "Score: " << ScoreImage(*image2) << std::endl;
+  uint32_t diff = 1;
+  while (diff) {
+    diff = OptimizeLut<4>(*image, &lut);
+    std::cout << "diff=" << diff << " error=" << ScoreImage(*lut.MapImage(*image)) << std::endl;
+  }
+
+  auto image2 = lut.MapImage(*image);
   WriteFile("test.png", HighlightClosest(*image2)->ToPng());
 }
diff --git a/piraw.h b/piraw.h
index 630a924..268177c 100644
--- a/piraw.h
+++ b/piraw.h
@@ -118,8 +118,8 @@ template <uint32_t X, uint32_t Y, uint32_t D, uint32_t A, uint32_t P>
 constexpr Color PiRaw<X, Y, D, A, P>::CombineRaw(uint32_t y0x0, uint32_t y0x1, uint32_t y1x0, uint32_t y1x1) {
   // Function is bit layout specific
   Color ret;
-  ret.r = y1x1;
-  ret.g = (y0x1 + y1x0) / 2;
-  ret.b = y0x0;
+  ret.at(0) = y1x1;
+  ret.at(1) = (y0x1 + y1x0) / 2;
+  ret.at(2) = y0x0;
   return ret;
 }