diff --git a/Makefile b/Makefile
index cecfca7..88111f5 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
 all: piphoto
 
-objects = piphoto.o color.o coord.o util.o
+objects = piphoto.o util.o
 
 piphoto: $(objects) Makefile
 	clang-3.9 -O3 -g -Weverything -Werror --std=c++1z --stdlib=libc++ -o piphoto $(objects) -lc++ -lunwind -lpng
diff --git a/color.cc b/color.cc
deleted file mode 100644
index 42f5dcb..0000000
--- a/color.cc
+++ /dev/null
@@ -1,13 +0,0 @@
-#include "color.h"
-
-#include <iomanip>
-
-std::ostream& operator<<(std::ostream& os, const Color& color) {
-  return os
-    << std::hex << std::setfill('0')
-    << "rgb("
-    << "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 b19a4f5..4bd057c 100644
--- a/color.h
+++ b/color.h
@@ -2,23 +2,48 @@
 
 #include <array>
 #include <cstdint>
+#include <iomanip>
 #include <iostream>
 
 #include "intmath.h"
 
-constexpr uint32_t kNumColors = UINT16_MAX;
+constexpr uint32_t kMaxColor = 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;
+template <uint32_t C>
+struct Color : public std::array<uint32_t, C> {
+  constexpr uint32_t Difference(const Color<C>& other) const;
+  constexpr Color<C> Interpolate(const Color<C>& other, uint32_t mul, uint32_t div) const;
 };
 
-constexpr uint32_t Color::Difference(const Color& other) const {
-  return (
-    AbsDiff(this->at(0), other.at(0)) +
-    AbsDiff(this->at(1), other.at(1)) +
-    AbsDiff(this->at(2), other.at(2))
-  );
+struct RgbColor : public Color<3> {};
+
+template <uint32_t C>
+constexpr uint32_t Color<C>::Difference(const Color<C>& other) const {
+  uint32_t diff = 0;
+  for (uint32_t c = 0; c < C; ++c) {
+    diff += AbsDiff(this->at(c), other.at(c));
+  }
+  return diff;
 }
 
-std::ostream& operator<<(std::ostream& os, const Color& color);
+template <uint32_t C>
+constexpr Color<C> Color<C>::Interpolate(const Color<C>& other, uint32_t mul, uint32_t div) const {
+  Color<C> ret;
+  for (uint32_t c = 0; c < C; ++c) {
+    ret.at(c) = ::Interpolate(this->at(c), other.at(c), mul, div);
+  }
+  return ret;
+}
+
+template <uint32_t C>
+std::ostream& operator<<(std::ostream& os, const Color<C>& color) {
+  os << std::hex << std::setfill('0') << "Color(";
+  for (uint32_t c = 0; c < C; ++c) {
+    os << "0x" << std::setw(4) << color.at(0);
+    if (c < C - 1) {
+      os << ", ";
+    }
+  }
+  return os << ")" << std::dec;
+}
diff --git a/colorchecker.h b/colorchecker.h
index aea22bc..7d4a2a8 100644
--- a/colorchecker.h
+++ b/colorchecker.h
@@ -14,36 +14,38 @@
 // 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}}},
+constexpr std::array<RgbColor, 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}}}},
 }};
 
-template <uint32_t X, uint32_t Y>
-std::array<Coord, kColorCheckerSrgb.size()> FindClosest(const Image<X, Y>& image) {
-  std::array<Coord, kColorCheckerSrgb.size()> closest;
+template <uint32_t X, uint32_t Y, uint32_t C>
+std::array<Coord<2>, kColorCheckerSrgb.size()> FindClosest(const Image<X, Y, C>& image) {
+  static_assert(C == 3);
+
+  std::array<Coord<2>, kColorCheckerSrgb.size()> closest;
   std::array<uint32_t, kColorCheckerSrgb.size()> diff;
   diff.fill(UINT32_MAX);
 
@@ -57,7 +59,7 @@ std::array<Coord, kColorCheckerSrgb.size()> FindClosest(const Image<X, Y>& image
         auto pixel_diff = pixel.Difference(kColorCheckerSrgb.at(cc));
         if (pixel_diff < diff.at(cc)) {
           diff.at(cc) = pixel_diff;
-          closest.at(cc) = {x, y};
+          closest.at(cc) = {{{x, y}}};
         }
       }
     }
@@ -66,8 +68,10 @@ std::array<Coord, kColorCheckerSrgb.size()> FindClosest(const Image<X, Y>& image
   return closest;
 }
 
-template <uint32_t X, uint32_t Y>
-uint32_t ScoreImage(const Image<X, Y>& image) {
+template <uint32_t X, uint32_t Y, uint32_t C>
+uint32_t ScoreImage(const Image<X, Y, C>& image) {
+  static_assert(C == 3);
+
   std::array<uint32_t, kColorCheckerSrgb.size()> diff;
   diff.fill(UINT32_MAX);
 
@@ -89,44 +93,48 @@ uint32_t ScoreImage(const Image<X, Y>& image) {
   return std::accumulate(diff.begin(), diff.end(), UINT32_C(0));
 }
 
-template <uint32_t X, uint32_t Y>
-std::unique_ptr<Image<X, Y>> HighlightClosest(const Image<X, Y>& image) {
-  auto out = std::make_unique<Image<X, Y>>(image);
+template <uint32_t X, uint32_t Y, uint32_t C>
+std::unique_ptr<Image<X, Y, C>> HighlightClosest(const Image<X, Y, C>& image) {
+  static_assert(C == 3);
+
+  auto out = std::make_unique<Image<X, Y, C>>(image);
 
   auto closest = FindClosest(*out);
   for (uint32_t cc = 0; cc < kColorCheckerSrgb.size(); ++cc) {
     const auto& coord = closest.at(cc);
     const auto& color = kColorCheckerSrgb.at(cc);
-    out->DrawSquare({std::max(5U, coord.x) - 5, std::max(5U, coord.y) - 5}, kBlack, 10);
-    out->DrawSquare({std::max(6U, coord.x) - 6, std::max(6U, coord.y) - 6}, color, 12);
-    out->DrawSquare({std::max(7U, coord.x) - 7, std::max(7U, coord.y) - 7}, color, 14);
-    out->DrawSquare({std::max(8U, coord.x) - 8, std::max(8U, coord.y) - 8}, color, 16);
-    out->DrawSquare({std::max(9U, coord.x) - 9, std::max(9U, coord.y) - 9}, kWhite, 18);
+    out->DrawSquare({{{std::max(5U, coord.at(0)) - 5, std::max(5U, coord.at(1)) - 5}}}, kBlack, 10);
+    out->DrawSquare({{{std::max(6U, coord.at(0)) - 6, std::max(6U, coord.at(1)) - 6}}}, color, 12);
+    out->DrawSquare({{{std::max(7U, coord.at(0)) - 7, std::max(7U, coord.at(1)) - 7}}}, color, 14);
+    out->DrawSquare({{{std::max(8U, coord.at(0)) - 8, std::max(8U, coord.at(1)) - 8}}}, color, 16);
+    out->DrawSquare({{{std::max(9U, coord.at(0)) - 9, std::max(9U, coord.at(1)) - 9}}}, kWhite, 18);
   }
   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) {
+template <uint32_t P, uint32_t LUT_X, uint32_t LUT_Y, uint32_t LUT_Z, uint32_t IMG_X, uint32_t IMG_Y, uint32_t C>
+uint32_t OptimizeLut(const Image<IMG_X, IMG_Y, C>& image, Lut3d<LUT_X, LUT_Y, LUT_Z>* lut) {
+  static_assert(C == 3);
+
   uint32_t diff = 0;
 
-  for (uint32_t r = 0; r < R; ++r) {
-    auto& rect = lut->at(r);
+  for (uint32_t x = 0; x < LUT_X; ++x) {
+    auto& rect = lut->at(x);
 
-    for (uint32_t g = 0; g < G; ++g) {
-      auto& row = rect.at(g);
+    for (uint32_t y = 0; y < LUT_Y; ++y) {
+      auto& row = rect.at(y);
 
-      for (uint32_t b = 0; b < B; ++b) {
-        auto& color = row.at(b);
+      for (uint32_t z = 0; z < LUT_Z; ++z) {
+        auto& color = row.at(z);
 
-        std::cout << "(" << r << ", " << g << ", " << b << ")" << std::endl;
+        std::cout << Coord<3>{{{x, y, z}}} << std::endl;
 
-        for (uint32_t c = 0; c < color.size(); ++c) {
+        for (uint32_t c = 0; c < C; ++c) {
           auto min = FindPossibleMinimum<uint32_t, uint32_t, 4>(
             0, UINT16_MAX,
-            [&image, &lut, r, g, b, c](uint32_t val) {
+            [&image, &lut, x, y, z, c](uint32_t val) {
               auto test_lut = *lut;
-              test_lut[r][g][b].at(c) = val;
+              test_lut.at(x).at(y).at(z).at(c) = val;
               return ScoreImage(*test_lut.MapImage(image));
             });
           std::cout << "\tC" << c << ": " << color.at(c) << " -> " << min << std::endl;
diff --git a/colors.h b/colors.h
index 85db66d..5d22f86 100644
--- a/colors.h
+++ b/colors.h
@@ -2,5 +2,5 @@
 
 #include "color.h"
 
-constexpr Color kBlack = {{{0x0000, 0x0000, 0x0000}}};
-constexpr Color kWhite = {{{0xffff, 0xffff, 0xffff}}};
+constexpr RgbColor kBlack = {{{{0x0000, 0x0000, 0x0000}}}};
+constexpr RgbColor kWhite = {{{{0xffff, 0xffff, 0xffff}}}};
diff --git a/coord.cc b/coord.cc
deleted file mode 100644
index 94fa9d6..0000000
--- a/coord.cc
+++ /dev/null
@@ -1,6 +0,0 @@
-#include "coord.h"
-
-std::ostream& operator<<(std::ostream& os, const Coord& coord) {
-  return os << "(" << coord.x << ", " << coord.y << ")";
-}
-
diff --git a/coord.h b/coord.h
index 53b8125..8196821 100644
--- a/coord.h
+++ b/coord.h
@@ -3,15 +3,17 @@
 #include <cstdint>
 #include <ostream>
 
-struct Coord {
-  uint32_t x;
-  uint32_t y;
-};
+template <uint32_t D>
+struct Coord : public std::array<uint32_t, D> {};
 
-struct Coord3d {
-  uint32_t r;
-  uint32_t g;
-  uint32_t b;
-};
-
-std::ostream& operator<<(std::ostream& os, const Coord& coord);
+template <uint32_t D>
+std::ostream& operator<<(std::ostream& os, const Coord<D>& coord) {
+  os << "(";
+  for (uint32_t d = 0; d < D; ++d) {
+    os << coord.at(d);
+    if (d < D - 1) {
+      os << ", ";
+    }
+  }
+  return os << ")";
+}
diff --git a/image.h b/image.h
index 8e8de25..888db41 100644
--- a/image.h
+++ b/image.h
@@ -8,56 +8,56 @@
 #include "color.h"
 #include "coord.h"
 
-template <uint32_t X, uint32_t Y>
-class Image : public std::array<std::array<Color, X>, Y> {
+template <uint32_t X, uint32_t Y, uint32_t C>
+class Image : public std::array<std::array<Color<C>, X>, Y> {
  public:
-  constexpr const Color& GetPixel(const Coord& coord) const;
+  constexpr const Color<C>& GetPixel(const Coord<2>& coord) const;
 
-  void SetPixel(const Coord& coord, const Color& color);
-  void DrawXLine(const Coord& start, const Color& color, uint32_t length);
-  void DrawYLine(const Coord& start, const Color& color, uint32_t length);
-  void DrawRectangle(const Coord& start, const Color& color, uint32_t x_length, uint32_t y_length);
-  void DrawSquare(const Coord& start, const Color& color, uint32_t length);
+  void SetPixel(const Coord<2>& coord, const Color<C>& color);
+  void DrawXLine(const Coord<2>& start, const Color<C>& color, uint32_t length);
+  void DrawYLine(const Coord<2>& start, const Color<C>& color, uint32_t length);
+  void DrawRectangle(const Coord<2>& start, const Color<C>& color, uint32_t x_length, uint32_t y_length);
+  void DrawSquare(const Coord<2>& start, const Color<C>& color, uint32_t length);
 
   std::string ToPng();
 };
 
-template <uint32_t X, uint32_t Y>
-constexpr const Color& Image<X, Y>::GetPixel(const Coord& coord) const {
-  return this->at(coord.y).at(coord.x);
+template <uint32_t X, uint32_t Y, uint32_t C>
+constexpr const Color<C>& Image<X, Y, C>::GetPixel(const Coord<2>& coord) const {
+  return this->at(coord.at(1)).at(coord.at(0));
 }
 
-template <uint32_t X, uint32_t Y>
-void Image<X, Y>::SetPixel(const Coord& coord, const Color& color) {
-  this->at(coord.y).at(coord.x) = color;
+template <uint32_t X, uint32_t Y, uint32_t C>
+void Image<X, Y, C>::SetPixel(const Coord<2>& coord, const Color<C>& color) {
+  this->at(coord.at(1)).at(coord.at(0)) = color;
 }
 
-template <uint32_t X, uint32_t Y>
-void Image<X, Y>::DrawXLine(const Coord& coord, const Color& color, uint32_t length) {
-  auto& row = this->at(coord.y);
+template <uint32_t X, uint32_t Y, uint32_t C>
+void Image<X, Y, C>::DrawXLine(const Coord<2>& coord, const Color<C>& color, uint32_t length) {
+  auto& row = this->at(coord.at(1));
 
-  for (uint32_t x = coord.x; x < std::min(X, coord.x + length); ++x) {
+  for (uint32_t x = coord.at(0); x < std::min(X, coord.at(0) + length); ++x) {
     row.at(x) = color;
   }
 }
 
-template <uint32_t X, uint32_t Y>
-void Image<X, Y>::DrawYLine(const Coord& coord, const Color& color, uint32_t length) {
-  for (uint32_t y = coord.y; y <= std::min(Y, coord.y + length); ++y) {
-    SetPixel({coord.x, y}, color);
+template <uint32_t X, uint32_t Y, uint32_t C>
+void Image<X, Y, C>::DrawYLine(const Coord<2>& coord, const Color<C>& color, uint32_t length) {
+  for (uint32_t y = coord.at(1); y <= std::min(Y, coord.at(1) + length); ++y) {
+    SetPixel({{{coord.at(0), y}}}, color);
   }
 }
 
-template <uint32_t X, uint32_t Y>
-void Image<X, Y>::DrawRectangle(const Coord& start, const Color& color, uint32_t x_length, uint32_t y_length) {
+template <uint32_t X, uint32_t Y, uint32_t C>
+void Image<X, Y, C>::DrawRectangle(const Coord<2>& start, const Color<C>& color, uint32_t x_length, uint32_t y_length) {
   DrawXLine(start, color, x_length);
-  DrawXLine({start.x, start.y + y_length}, color, x_length);
+  DrawXLine({{{start.at(0), start.at(1) + y_length}}}, color, x_length);
   DrawYLine(start, color, y_length);
-  DrawYLine({start.x + x_length, start.y}, color, y_length);
+  DrawYLine({{{start.at(0) + x_length, start.at(1)}}}, color, y_length);
 }
 
-template <uint32_t X, uint32_t Y>
-void Image<X, Y>::DrawSquare(const Coord& start, const Color& color, uint32_t length) {
+template <uint32_t X, uint32_t Y, uint32_t C>
+void Image<X, Y, C>::DrawSquare(const Coord<2>& start, const Color<C>& color, uint32_t length) {
   DrawRectangle(start, color, length, length);
 }
 
@@ -66,8 +66,10 @@ static inline void WriteCallback(png_structp png_ptr, png_bytep data, png_size_t
   dest->append(reinterpret_cast<char*>(data), length);
 }
 
-template <uint32_t X, uint32_t Y>
-std::string Image<X, Y>::ToPng() {
+template <uint32_t X, uint32_t Y, uint32_t C>
+std::string Image<X, Y, C>::ToPng() {
+  static_assert(C == 3);  // PNG only supports RGB
+
   std::string ret;
 
   auto png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
diff --git a/intmath.h b/intmath.h
index 79cd229..59eb3ae 100644
--- a/intmath.h
+++ b/intmath.h
@@ -4,3 +4,12 @@ template <typename T>
 constexpr T AbsDiff(T a, T b) {
   return (a > b) ? (a - b) : (b - a);
 }
+
+template <typename T>
+constexpr T Interpolate(T val0, T val1, T mul, T div) {
+  if (val1 > val0) {
+    return val0 + ((mul * (val1 - val0)) / div);
+  } else {
+    return val1 - (((div - mul) * (val0 - val1)) / div);
+  }
+}
diff --git a/lut.h b/lut.h
index dd2f1f7..ee6c48c 100644
--- a/lut.h
+++ b/lut.h
@@ -3,22 +3,20 @@
 #include "color.h"
 #include "coord.h"
 
-template <uint32_t R, uint32_t G, uint32_t B>
-class Lut3d : public std::array<std::array<std::array<Color, B>, G>, R> {
+// Hardcoded to Color<3>, so color dimensions == LUT dimensions
+template <uint32_t X, uint32_t Y, uint32_t Z>
+class Lut3d : public std::array<std::array<std::array<Color<3>, X>, Y>, Z> {
  public:
-  static Lut3d<R, G, B> Identity();
+  static Lut3d<X, Y, Z> Identity();
 
-  Color MapColor(const Color& in) const;
+  Color<3> MapColor(const Color<3>& in) const;
 
-  template <uint32_t X, uint32_t Y>
-  std::unique_ptr<Image<X, Y>> MapImage(const Image<X, Y>& in) const;
+  template <uint32_t IMG_X, uint32_t IMG_Y, uint32_t C>
+  std::unique_ptr<Image<IMG_X, IMG_Y, C>> MapImage(const Image<IMG_X, IMG_Y, C>& in) const;
 
  private:
-  constexpr static Color InterpolateColor(const Color& i0, const Color& i1, uint32_t mul, uint32_t div);
-  constexpr static uint32_t Interpolate(uint32_t i0, uint32_t i1, uint32_t mul, uint32_t div);
-
   // Return value is (root_indices, remainders)
-  constexpr static std::pair<Coord3d, Coord3d> FindRoot(const Color& in);
+  constexpr static std::pair<Coord<3>, Coord<3>> FindRoot(const Color<3>& in);
   constexpr static std::pair<uint32_t, uint32_t> FindChannelRoot(uint32_t value, uint32_t points);
 
   constexpr static uint32_t BlockSize(uint32_t points);
@@ -27,22 +25,22 @@ class Lut3d : public std::array<std::array<std::array<Color, B>, G>, R> {
 // Minimum size LUT
 typedef Lut3d<2, 2, 2> MinimalLut3d;
 
-template <uint32_t R, uint32_t G, uint32_t B>
-Lut3d<R, G, B> Lut3d<R, G, B>::Identity() {
-  Lut3d<R, G, B> ret;
+template <uint32_t X, uint32_t Y, uint32_t Z>
+Lut3d<X, Y, Z> Lut3d<X, Y, Z>::Identity() {
+  Lut3d<X, Y, Z> ret;
 
-  Color color;
-  for (uint32_t r = 0; r < R; ++r) {
-    auto& rect = ret.at(r);
-    color.at(0) = std::min(kNumColors - 1, BlockSize(R) * r);
+  Color<3> color;
+  for (uint32_t x = 0; x < X; ++x) {
+    auto& rect = ret.at(x);
+    color.at(0) = std::min(kMaxColor, BlockSize(X) * x);
 
-    for (uint32_t g = 0; g < G; ++g) {
-      auto& row = rect.at(g);
-      color.at(1) = std::min(kNumColors - 1, BlockSize(G) * g);
+    for (uint32_t y = 0; y < Y; ++y) {
+      auto& row = rect.at(y);
+      color.at(1) = std::min(kMaxColor, BlockSize(Y) * y);
 
-      for (uint32_t b = 0; b < B; ++b) {
-        color.at(2) = std::min(kNumColors - 1, BlockSize(B) * b);
-        row.at(b) = color;
+      for (uint32_t z = 0; z < Z; ++z) {
+        color.at(2) = std::min(kMaxColor, BlockSize(Z) * z);
+        row.at(z) = color;
       }
     }
   }
@@ -50,51 +48,51 @@ Lut3d<R, G, B> Lut3d<R, G, B>::Identity() {
   return ret;
 }
 
-template <uint32_t R, uint32_t G, uint32_t B>
-Color Lut3d<R, G, B>::MapColor(const Color& in) const {
+template <uint32_t X, uint32_t Y, uint32_t Z>
+Color<3> Lut3d<X, Y, Z>::MapColor(const Color<3>& in) const {
   const auto root_rem = FindRoot(in);
   const auto& root = root_rem.first;
   const auto& rem = root_rem.second;
 
   // https://en.wikipedia.org/wiki/Trilinear_interpolation
-  auto inter00 = InterpolateColor(
-    this->at(root.r + 0).at(root.g + 0).at(root.b + 0),
-    this->at(root.r + 1).at(root.g + 0).at(root.b + 0),
-    rem.r,
-    BlockSize(R));
+  auto inter00 =
+    this->at(root.at(0) + 0).at(root.at(1) + 0).at(root.at(2) + 0).Interpolate(
+    this->at(root.at(0) + 1).at(root.at(1) + 0).at(root.at(2) + 0),
+    rem.at(0),
+    BlockSize(X));
 
-  auto inter01 = InterpolateColor(
-    this->at(root.r + 0).at(root.g + 0).at(root.b + 1),
-    this->at(root.r + 1).at(root.g + 0).at(root.b + 1),
-    rem.r,
-    BlockSize(R));
+  auto inter01 =
+    this->at(root.at(0) + 0).at(root.at(1) + 0).at(root.at(2) + 1).Interpolate(
+    this->at(root.at(0) + 1).at(root.at(1) + 0).at(root.at(2) + 1),
+    rem.at(0),
+    BlockSize(X));
 
-  auto inter10 = InterpolateColor(
-    this->at(root.r + 0).at(root.g + 1).at(root.b + 0),
-    this->at(root.r + 1).at(root.g + 1).at(root.b + 0),
-    rem.r,
-    BlockSize(R));
+  auto inter10 =
+    this->at(root.at(0) + 0).at(root.at(1) + 1).at(root.at(2) + 0).Interpolate(
+    this->at(root.at(0) + 1).at(root.at(1) + 1).at(root.at(2) + 0),
+    rem.at(0),
+    BlockSize(X));
 
-  auto inter11 = InterpolateColor(
-    this->at(root.r + 0).at(root.g + 1).at(root.b + 1),
-    this->at(root.r + 1).at(root.g + 1).at(root.b + 1),
-    rem.r,
-    BlockSize(R));
+  auto inter11 =
+    this->at(root.at(0) + 0).at(root.at(1) + 1).at(root.at(2) + 1).Interpolate(
+    this->at(root.at(0) + 1).at(root.at(1) + 1).at(root.at(2) + 1),
+    rem.at(0),
+    BlockSize(X));
 
-  auto inter0 = InterpolateColor(inter00, inter10, rem.g, BlockSize(G));
-  auto inter1 = InterpolateColor(inter01, inter11, rem.g, BlockSize(G));
+  auto inter0 = inter00.Interpolate(inter10, rem.at(1), BlockSize(Y));
+  auto inter1 = inter01.Interpolate(inter11, rem.at(1), BlockSize(Y));
 
-  return InterpolateColor(inter0, inter1, rem.b, BlockSize(B));
+  return inter0.Interpolate(inter1, rem.at(2), BlockSize(Z));
 }
 
-template <uint32_t R, uint32_t G, uint32_t B>
-template <uint32_t X, uint32_t Y>
-std::unique_ptr<Image<X, Y>> Lut3d<R, G, B>::MapImage(const Image<X, Y>& in) const {
-  auto out = std::make_unique<Image<X, Y>>();
+template <uint32_t X, uint32_t Y, uint32_t Z>
+template <uint32_t IMG_X, uint32_t IMG_Y, uint32_t C>
+std::unique_ptr<Image<IMG_X, IMG_Y, C>> Lut3d<X, Y, Z>::MapImage(const Image<IMG_X, IMG_Y, C>& in) const {
+  auto out = std::make_unique<Image<IMG_X, IMG_Y, C>>();
 
-  for (uint32_t y = 0; y < Y; ++y) {
-    for (uint32_t x = 0; x < X; ++x) {
-      Coord coord = {x, y};
+  for (uint32_t y = 0; y < IMG_Y; ++y) {
+    for (uint32_t x = 0; x < IMG_X; ++x) {
+      Coord<2> coord = {{{x, y}}};
       out->SetPixel(coord, MapColor(in.GetPixel(coord)));
     }
   }
@@ -102,37 +100,19 @@ std::unique_ptr<Image<X, Y>> Lut3d<R, G, B>::MapImage(const Image<X, Y>& in) con
   return out;
 }
 
-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.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) {
-  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.at(0), R);
-  auto root_g = FindChannelRoot(in.at(1), G);
-  auto root_b = FindChannelRoot(in.at(2), B);
+template <uint32_t X, uint32_t Y, uint32_t Z>
+constexpr std::pair<Coord<3>, Coord<3>> Lut3d<X, Y, Z>::FindRoot(const Color<3>& in) {
+  auto root_x = FindChannelRoot(in.at(0), X);
+  auto root_y = FindChannelRoot(in.at(1), Y);
+  auto root_z = FindChannelRoot(in.at(2), Z);
   return {
-    {root_r.first, root_g.first, root_b.first},
-    {root_r.second, root_g.second, root_b.second},
+    {{{root_x.first, root_y.first, root_z.first}}},
+    {{{root_x.second, root_y.second, root_z.second}}},
   };
 }
 
-template <uint32_t R, uint32_t G, uint32_t B>
-constexpr std::pair<uint32_t, uint32_t> Lut3d<R, G, B>::FindChannelRoot(const uint32_t value, const uint32_t points) {
+template <uint32_t X, uint32_t Y, uint32_t Z>
+constexpr std::pair<uint32_t, uint32_t> Lut3d<X, Y, Z>::FindChannelRoot(const uint32_t value, const uint32_t points) {
   // points - 1 is the last point index. Since we're going to fidn the cube
   // around this point by adding to the root, we need to be at least 1 less
   // than that.
@@ -140,7 +120,7 @@ constexpr std::pair<uint32_t, uint32_t> Lut3d<R, G, B>::FindChannelRoot(const ui
   return std::make_pair(index, value - (index * BlockSize(points)));
 }
 
-template <uint32_t R, uint32_t G, uint32_t B>
-constexpr uint32_t Lut3d<R, G, B>::BlockSize(uint32_t points) {
-  return kNumColors / (points - 1);
+template <uint32_t X, uint32_t Y, uint32_t Z>
+constexpr uint32_t Lut3d<X, Y, Z>::BlockSize(uint32_t points) {
+  return (kMaxColor + 1) / (points - 1);
 }
diff --git a/piraw.h b/piraw.h
index 268177c..706306c 100644
--- a/piraw.h
+++ b/piraw.h
@@ -10,15 +10,15 @@ namespace std {
 using string_view = experimental::string_view;
 }
 
-template <uint32_t X, uint32_t Y, uint32_t D, uint32_t A, uint32_t P>
+template <uint32_t X, uint32_t Y, uint32_t C, uint32_t D, uint32_t A, uint32_t P>
 class PiRaw {
  public:
   PiRaw() = delete;
   PiRaw(const PiRaw&) = delete;
   PiRaw(PiRaw&&) = delete;
 
-  static std::unique_ptr<Image<X / 2, Y / 2>> FromJpeg(const std::string_view& jpeg);
-  static std::unique_ptr<Image<X / 2, Y / 2>> FromRaw(const std::string_view& raw);
+  static std::unique_ptr<Image<X / 2, Y / 2, C>> FromJpeg(const std::string_view& jpeg);
+  static std::unique_ptr<Image<X / 2, Y / 2, C>> FromRaw(const std::string_view& raw);
 
  private:
   static constexpr uint32_t kJpegHeaderBytes = 32768;
@@ -36,27 +36,30 @@ class PiRaw {
   typedef std::array<uint32_t, kPixelsPerChunk> Chunk;
 
   static constexpr Chunk GetChunk(const std::string_view& raw, const uint32_t x_chunk, const uint32_t y);
-  static constexpr Color CombineRaw(uint32_t y0x0, uint32_t y0x1, uint32_t y1x0, uint32_t y1x1);
+  static constexpr Color<C> CombineRaw(uint32_t y0x0, uint32_t y0x1, uint32_t y1x0, uint32_t y1x1);
 };
 
-typedef PiRaw<3280, 2464, 10, 16, 2> PiRaw2;
+typedef PiRaw<3280, 2464, 3, 10, 16, 2> PiRaw2;
+
+template <uint32_t X, uint32_t Y, uint32_t C, uint32_t D, uint32_t A, uint32_t P>
+typename std::unique_ptr<Image<X / 2, Y / 2, C>> PiRaw<X, Y, C, D, A, P>::FromJpeg(const std::string_view& jpeg) {
+  static_assert(C == 3);
 
-template <uint32_t X, uint32_t Y, uint32_t D, uint32_t A, uint32_t P>
-typename std::unique_ptr<Image<X / 2, Y / 2>> PiRaw<X, Y, D, A, P>::FromJpeg(const std::string_view& jpeg) {
   auto container_len = GetRawBytes() + kJpegHeaderBytes;
   assert(jpeg.substr(jpeg.size() - container_len, 4) == kJpegHeaderMagic);
   return FromRaw(jpeg.substr(jpeg.size() - GetRawBytes(), GetRawBytes()));
 }
 
-template <uint32_t X, uint32_t Y, uint32_t D, uint32_t A, uint32_t P>
-typename std::unique_ptr<Image<X / 2, Y / 2>> PiRaw<X, Y, D, A, P>::FromRaw(const std::string_view& raw) {
+template <uint32_t X, uint32_t Y, uint32_t C, uint32_t D, uint32_t A, uint32_t P>
+typename std::unique_ptr<Image<X / 2, Y / 2, C>> PiRaw<X, Y, C, D, A, P>::FromRaw(const std::string_view& raw) {
+  static_assert(C == 3);
   static_assert(X % 2 == 0);
   static_assert(Y % 2 == 0);
   static_assert(kPixelsPerChunk == 4);
 
   assert(raw.size() == GetRawBytes());
 
-  auto image = std::make_unique<Image<X / 2, Y / 2>>();
+  auto image = std::make_unique<Image<X / 2, Y / 2, C>>();
 
   for (uint32_t y = 0, out_y = 0; y < Y; y += 2, ++out_y) {
     for (uint32_t x_chunk = 0, out_x = 0; x_chunk < X / kPixelsPerChunk; ++x_chunk, out_x += kPixelsPerChunk / 2) {
@@ -69,34 +72,35 @@ typename std::unique_ptr<Image<X / 2, Y / 2>> PiRaw<X, Y, D, A, P>::FromRaw(cons
   return image;
 }
 
-template <uint32_t X, uint32_t Y, uint32_t D, uint32_t A, uint32_t P>
-constexpr uint32_t PiRaw<X, Y, D, A, P>::GetRawBytes() {
+template <uint32_t X, uint32_t Y, uint32_t C, uint32_t D, uint32_t A, uint32_t P>
+constexpr uint32_t PiRaw<X, Y, C, D, A, P>::GetRawBytes() {
   return GetRowBytes() * GetNumRows();
 }
 
-template <uint32_t X, uint32_t Y, uint32_t D, uint32_t A, uint32_t P>
-constexpr uint32_t PiRaw<X, Y, D, A, P>::GetRowBytes() {
+template <uint32_t X, uint32_t Y, uint32_t C, uint32_t D, uint32_t A, uint32_t P>
+constexpr uint32_t PiRaw<X, Y, C, D, A, P>::GetRowBytes() {
   return Align(Align(X + P) * D / kBitsPerByte);
 }
 
-template <uint32_t X, uint32_t Y, uint32_t D, uint32_t A, uint32_t P>
-constexpr uint32_t PiRaw<X, Y, D, A, P>::GetNumRows() {
+template <uint32_t X, uint32_t Y, uint32_t C, uint32_t D, uint32_t A, uint32_t P>
+constexpr uint32_t PiRaw<X, Y, C, D, A, P>::GetNumRows() {
   return Align(Y + P);
 }
 
-template <uint32_t X, uint32_t Y, uint32_t D, uint32_t A, uint32_t P>
-constexpr uint32_t PiRaw<X, Y, D, A, P>::GetChunkBytes() {
+template <uint32_t X, uint32_t Y, uint32_t C, uint32_t D, uint32_t A, uint32_t P>
+constexpr uint32_t PiRaw<X, Y, C, D, A, P>::GetChunkBytes() {
   return D * kPixelsPerChunk / kBitsPerByte;
 }
 
-template <uint32_t X, uint32_t Y, uint32_t D, uint32_t A, uint32_t P>
-constexpr uint32_t PiRaw<X, Y, D, A, P>::Align(uint32_t val) {
+template <uint32_t X, uint32_t Y, uint32_t C, uint32_t D, uint32_t A, uint32_t P>
+constexpr uint32_t PiRaw<X, Y, C, D, A, P>::Align(uint32_t val) {
   return (~(A - 1)) & ((val) + (A - 1));
 }
 
-template <uint32_t X, uint32_t Y, uint32_t D, uint32_t A, uint32_t P>
-constexpr typename PiRaw<X, Y, D, A, P>::Chunk PiRaw<X, Y, D, A, P>::GetChunk(const std::string_view& raw, const uint32_t x_chunk, const uint32_t y) {
+template <uint32_t X, uint32_t Y, uint32_t C, uint32_t D, uint32_t A, uint32_t P>
+constexpr typename PiRaw<X, Y, C, D, A, P>::Chunk PiRaw<X, Y, C, D, A, P>::GetChunk(const std::string_view& raw, const uint32_t x_chunk, const uint32_t y) {
   // Function is bit depth & layout specific
+  static_assert(C == 3);
   static_assert(D == 10);
 
   uint32_t start = y * GetRowBytes() + x_chunk * GetChunkBytes();
@@ -114,10 +118,12 @@ constexpr typename PiRaw<X, Y, D, A, P>::Chunk PiRaw<X, Y, D, A, P>::GetChunk(co
   return ret;
 }
 
-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) {
+template <uint32_t X, uint32_t Y, uint32_t C, uint32_t D, uint32_t A, uint32_t P>
+constexpr Color<C> PiRaw<X, Y, C, D, A, P>::CombineRaw(uint32_t y0x0, uint32_t y0x1, uint32_t y1x0, uint32_t y1x1) {
   // Function is bit layout specific
-  Color ret;
+  static_assert(C == 3);
+
+  Color<C> ret;
   ret.at(0) = y1x1;
   ret.at(1) = (y0x1 + y1x0) / 2;
   ret.at(2) = y0x0;