diff --git a/color.h b/color.h
index 7b0ac50..ffac6c1 100644
--- a/color.h
+++ b/color.h
@@ -2,6 +2,8 @@
 
 #include <cstdint>
 
+constexpr uint32_t kNumColors = (1 << 16);
+
 struct Color {
   // 32-bit for compiler convenience, but values are 16-bit
   uint32_t r;
diff --git a/image.h b/image.h
index 571b81b..9059079 100644
--- a/image.h
+++ b/image.h
@@ -76,7 +76,7 @@ std::string Image<X, Y>::ToPng() {
   assert(info_ptr);
 
   png_set_write_fn(png_ptr, &ret, &WriteCallback, nullptr);
-  png_set_IHDR(png_ptr, info_ptr, X / 2, Y / 2,
+  png_set_IHDR(png_ptr, info_ptr, X, Y,
     16, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
     PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
 
diff --git a/lut.h b/lut.h
new file mode 100644
index 0000000..63adf1d
--- /dev/null
+++ b/lut.h
@@ -0,0 +1,100 @@
+#pragma once
+
+#include "color.h"
+
+struct Coord3d {
+  uint32_t r;
+  uint32_t g;
+  uint32_t b;
+};
+
+template <uint32_t R, uint32_t G, uint32_t B>
+class Lut3d : public std::array<std::array<std::array<Color, B>, G>, R> {
+ public:
+  Color MapColor(const Color& 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<uint32_t, uint32_t> FindChannelRoot(uint32_t value, uint32_t points);
+
+  constexpr static uint32_t BlockSize(uint32_t points);
+};
+
+template <uint32_t R, uint32_t G, uint32_t B>
+Color Lut3d<R, G, B>::MapColor(const Color& 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 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 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 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 inter0 = InterpolateColor(inter00, inter10, rem.g, BlockSize(G));
+  auto inter1 = InterpolateColor(inter01, inter11, rem.g, BlockSize(G));
+
+  return InterpolateColor(inter0, inter1, rem.b, BlockSize(B));
+}
+
+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),
+  };
+}
+
+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) {
+  return i0 + ((mul * 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.r, G);
+  auto root_b = FindChannelRoot(in.r, B);
+  return {
+    {root_r.first, root_g.first, root_b.first},
+    {root_r.second, root_g.second, root_b.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) {
+  // 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.
+  return std::make_pair(std::min(points - 2, value / BlockSize(points)), value % BlockSize(points));;
+}
+
+template <uint32_t R, uint32_t G, uint32_t B>
+constexpr uint32_t BlockSize(uint32_t points) {
+  return kNumColors / (points - 1);
+}
diff --git a/piphoto.cc b/piphoto.cc
index 35895fd..ae3fbab 100644
--- a/piphoto.cc
+++ b/piphoto.cc
@@ -5,6 +5,7 @@
 #include "colors.h"
 #include "coord.h"
 #include "image.h"
+#include "lut.h"
 #include "piraw.h"
 #include "util.h"
 
diff --git a/piraw.h b/piraw.h
index b0105dc..937897a 100644
--- a/piraw.h
+++ b/piraw.h
@@ -17,8 +17,8 @@ class PiRaw {
   PiRaw(const PiRaw&) = delete;
   PiRaw(PiRaw&&) = delete;
 
-  static std::unique_ptr<Image<X, Y>> FromJpeg(const std::string_view& jpeg);
-  static std::unique_ptr<Image<X, Y>> FromRaw(const std::string_view& raw);
+  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);
 
  private:
   static constexpr uint32_t kJpegHeaderBytes = 32768;
@@ -42,21 +42,21 @@ class PiRaw {
 typedef PiRaw<3280, 2464, 10, 16, 2> PiRaw2;
 
 template <uint32_t X, uint32_t Y, uint32_t D, uint32_t A, uint32_t P>
-typename std::unique_ptr<Image<X, Y>> PiRaw<X, Y, D, A, P>::FromJpeg(const std::string_view& jpeg) {
+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, Y>> PiRaw<X, Y, D, A, P>::FromRaw(const std::string_view& raw) {
+typename std::unique_ptr<Image<X / 2, Y / 2>> PiRaw<X, Y, D, A, P>::FromRaw(const std::string_view& raw) {
   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, Y>>();
+  auto image = std::make_unique<Image<X / 2, Y / 2>>();
 
   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) {