diff --git a/internal/ceres/cubic_interpolation_test.cc b/internal/ceres/cubic_interpolation_test.cc
index 0b4bb12..04a7b33 100644
--- a/internal/ceres/cubic_interpolation_test.cc
+++ b/internal/ceres/cubic_interpolation_test.cc
@@ -31,31 +31,147 @@
 #include "ceres/cubic_interpolation.h"
 
 #include "ceres/jet.h"
+#include "ceres/internal/scoped_ptr.h"
 #include "glog/logging.h"
 #include "gtest/gtest.h"
 
 namespace ceres {
 namespace internal {
 
-TEST(CubicInterpolator, NeedsAtleastTwoValues) {
-  double x[] = {1};
-  EXPECT_DEATH_IF_SUPPORTED(CubicInterpolator c(x, 0), "num_values > 1");
-  EXPECT_DEATH_IF_SUPPORTED(CubicInterpolator c(x, 1), "num_values > 1");
+static const double kTolerance = 1e-12;
+
+TEST(Array1D, OneDataDimension) {
+  int x[] = {1, 2, 3};
+  Array1D<int, 1> array(x, 3);
+  for (int i = 0; i < 3; ++i) {
+    double value;
+    array.GetValue(i, &value);
+    EXPECT_EQ(value, static_cast<double>(i + 1));
+  }
 }
 
-static const double kTolerance = 1e-12;
+TEST(Array1D, TwoDataDimensionIntegerDataInterleaved) {
+  int x[] = {1, 5,
+             2, 6,
+             3, 7};
+
+  Array1D<int, 2, true> array(x, 3);
+  for (int i = 0; i < 3; ++i) {
+    double value[2];
+    array.GetValue(i, value);
+    EXPECT_EQ(value[0], static_cast<double>(i + 1));
+    EXPECT_EQ(value[1], static_cast<double>(i + 5));
+  }
+}
+
+TEST(Array1D, TwoDataDimensionIntegerDataStacked) {
+  int x[] = {1, 2, 3,
+             5, 6, 7};
+
+  Array1D<int, 2, false> array(x, 3);
+  for (int i = 0; i < 3; ++i) {
+    double value[2];
+    array.GetValue(i, value);
+    EXPECT_EQ(value[0], static_cast<double>(i + 1));
+    EXPECT_EQ(value[1], static_cast<double>(i + 5));
+  }
+}
+
+TEST(Array2D, OneDataDimensionRowMajor) {
+  int x[] = {1, 2, 3,
+             2, 3, 4};
+  Array2D<int, 1, true, true> array(x, 2, 3);
+  for (int r = 0; r < 2; ++r) {
+    for (int c = 0; c < 3; ++c) {
+      double value;
+      array.GetValue(r, c, &value);
+      EXPECT_EQ(value, static_cast<double>(r + c + 1));
+    }
+  }
+}
+
+TEST(Array2D, TwoDataDimensionRowMajorInterleaved) {
+  int x[] = {1, 4, 2, 8, 3, 12,
+             2, 8, 3, 12, 4, 16};
+  Array2D<int, 2, true, true> array(x, 2, 3);
+  for (int r = 0; r < 2; ++r) {
+    for (int c = 0; c < 3; ++c) {
+      double value[2];
+      array.GetValue(r, c, value);
+      EXPECT_EQ(value[0], static_cast<double>(r + c + 1));
+      EXPECT_EQ(value[1], static_cast<double>(4 *(r + c + 1)));
+    }
+  }
+}
+
+TEST(Array2D, TwoDataDimensionRowMajorStacked) {
+  int x[] = {1,  2,  3,
+             2,  3,  4,
+             4,  8, 12,
+             8, 12, 16};
+  Array2D<int, 2, true, false> array(x, 2, 3);
+  for (int r = 0; r < 2; ++r) {
+    for (int c = 0; c < 3; ++c) {
+      double value[2];
+      array.GetValue(r, c, value);
+      EXPECT_EQ(value[0], static_cast<double>(r + c + 1));
+      EXPECT_EQ(value[1], static_cast<double>(4 *(r + c + 1)));
+    }
+  }
+}
+
+TEST(Array2D, TwoDataDimensionColMajorInterleaved) {
+  int x[] = { 1,  4, 2,  8,
+              2,  8, 3, 12,
+              3, 12, 4, 16};
+  Array2D<int, 2, false, true> array(x, 2, 3);
+  for (int r = 0; r < 2; ++r) {
+    for (int c = 0; c < 3; ++c) {
+      double value[2];
+      array.GetValue(r, c, value);
+      EXPECT_EQ(value[0], static_cast<double>(r + c + 1));
+      EXPECT_EQ(value[1], static_cast<double>(4 *(r + c + 1)));
+    }
+  }
+}
+
+TEST(Array2D, TwoDataDimensionColMajorStacked) {
+  int x[] = {1,   2,
+             2,   3,
+             3,   4,
+             4,   8,
+             8,  12,
+             12, 16};
+  Array2D<int, 2, false, false> array(x, 2, 3);
+  for (int r = 0; r < 2; ++r) {
+    for (int c = 0; c < 3; ++c) {
+      double value[2];
+      array.GetValue(r, c, value);
+      EXPECT_EQ(value[0], static_cast<double>(r + c + 1));
+      EXPECT_EQ(value[1], static_cast<double>(4 *(r + c + 1)));
+    }
+  }
+}
+
 
 class CubicInterpolatorTest : public ::testing::Test {
  public:
+  template <int kDataDimension>
   void RunPolynomialInterpolationTest(const double a,
                                       const double b,
                                       const double c,
                                       const double d) {
+    values_.reset(new double[kDataDimension * kNumSamples]);
+
     for (int x = 0; x < kNumSamples; ++x) {
-      values_[x] = a * x * x * x + b * x * x + c * x + d;
+      for (int dim = 0; dim < kDataDimension; ++dim) {
+      values_[x * kDataDimension + dim] =
+          (dim * dim  + 1) * (a  * x * x * x + b * x * x + c * x + d);
+      }
     }
 
-    CubicInterpolator interpolator(values_, kNumSamples);
+    Array1D<double, kDataDimension> array(values_.get(), kNumSamples);
+    CubicInterpolator<Array1D<double, kDataDimension> > interpolator(array);
 
     // Check values in the all the cells but the first and the last
     // ones. In these cells, the interpolated function values should
@@ -66,46 +182,63 @@
     // function values and its derivatives not to match.
     for (int j = 0; j < kNumTestSamples; ++j) {
       const double x = 1.0 + 7.0 / (kNumTestSamples - 1) * j;
-      const double expected_f = a * x * x * x + b * x * x + c * x + d;
-      const double expected_dfdx = 3.0 * a * x * x + 2.0 * b * x + c;
-      double f, dfdx;
+      double expected_f[kDataDimension], expected_dfdx[kDataDimension];
+      double f[kDataDimension], dfdx[kDataDimension];
 
-      EXPECT_TRUE(interpolator.Evaluate(x, &f, &dfdx));
-      EXPECT_NEAR(f, expected_f, kTolerance)
-          << "x: " << x
-          << " actual f(x): " << expected_f
-          << " estimated f(x): " << f;
-      EXPECT_NEAR(dfdx, expected_dfdx, kTolerance)
-          << "x: " << x
-          << " actual df(x)/dx: " << expected_dfdx
-          << " estimated df(x)/dx: " << dfdx;
+      for (int dim = 0; dim < kDataDimension; ++dim) {
+        expected_f[dim] =
+            (dim * dim  + 1) * (a  * x * x * x + b * x * x + c * x + d);
+        expected_dfdx[dim] = (dim * dim + 1) * (3.0 * a * x * x + 2.0 * b * x + c);
+      }
+
+      EXPECT_TRUE(interpolator.Evaluate(x, f, dfdx));
+      for (int dim = 0; dim < kDataDimension; ++dim) {
+        EXPECT_NEAR(f[dim], expected_f[dim], kTolerance)
+            << "x: " << x << " dim: " << dim
+            << " actual f(x): " << expected_f[dim]
+            << " estimated f(x): " << f[dim];
+        EXPECT_NEAR(dfdx[dim], expected_dfdx[dim], kTolerance)
+            << "x: " << x << " dim: " << dim
+            << " actual df(x)/dx: " << expected_dfdx[dim]
+            << " estimated df(x)/dx: " << dfdx[dim];
+      }
     }
   }
 
  private:
   static const int kNumSamples = 10;
   static const int kNumTestSamples = 100;
-  double values_[kNumSamples];
+  scoped_array<double> values_;
 };
 
 TEST_F(CubicInterpolatorTest, ConstantFunction) {
-  RunPolynomialInterpolationTest(0.0, 0.0, 0.0, 0.5);
+  RunPolynomialInterpolationTest<1>(0.0, 0.0, 0.0, 0.5);
+  RunPolynomialInterpolationTest<2>(0.0, 0.0, 0.0, 0.5);
+  RunPolynomialInterpolationTest<3>(0.0, 0.0, 0.0, 0.5);
 }
 
 TEST_F(CubicInterpolatorTest, LinearFunction) {
-  RunPolynomialInterpolationTest(0.0, 0.0, 1.0, 0.5);
+  RunPolynomialInterpolationTest<1>(0.0, 0.0, 1.0, 0.5);
+  RunPolynomialInterpolationTest<2>(0.0, 0.0, 1.0, 0.5);
+  RunPolynomialInterpolationTest<3>(0.0, 0.0, 1.0, 0.5);
 }
 
 TEST_F(CubicInterpolatorTest, QuadraticFunction) {
-  RunPolynomialInterpolationTest(0.0, 0.4, 1.0, 0.5);
+  RunPolynomialInterpolationTest<1>(0.0, 0.4, 1.0, 0.5);
+  RunPolynomialInterpolationTest<2>(0.0, 0.4, 1.0, 0.5);
+  RunPolynomialInterpolationTest<3>(0.0, 0.4, 1.0, 0.5);
 }
 
+
 TEST(CubicInterpolator, JetEvaluation) {
-  const double values[] = {1.0, 2.0, 2.0, 3.0};
-  CubicInterpolator interpolator(values, 4);
-  double f, dfdx;
+  const double values[] = {1.0, 2.0, 2.0, 5.0, 3.0, 9.0, 2.0, 7.0};
+
+  Array1D<double, 2, true> array(values, 4);
+  CubicInterpolator<Array1D<double, 2, true> > interpolator(array);
+
+  double f[2], dfdx[2];
   const double x = 2.5;
-  EXPECT_TRUE(interpolator.Evaluate(x, &f, &dfdx));
+  EXPECT_TRUE(interpolator.Evaluate(x, f, dfdx));
 
   // Create a Jet with the same scalar part as x, so that the output
   // Jet will be evaluated at x.
@@ -116,42 +249,48 @@
   x_jet.v(2) = 1.2;
   x_jet.v(3) = 1.3;
 
-  Jet<double, 4> f_jet;
-  EXPECT_TRUE(interpolator.Evaluate(x_jet, &f_jet));
+  Jet<double, 4> f_jets[2];
+  EXPECT_TRUE(interpolator.Evaluate(x_jet, f_jets));
 
   // Check that the scalar part of the Jet is f(x).
-  EXPECT_EQ(f_jet.a, f);
+  EXPECT_EQ(f_jets[0].a, f[0]);
+  EXPECT_EQ(f_jets[1].a, f[1]);
 
   // Check that the derivative part of the Jet is dfdx * x_jet.v
   // by the chain rule.
-  EXPECT_EQ((f_jet.v - dfdx * x_jet.v).norm(), 0.0);
+  EXPECT_NEAR((f_jets[0].v - dfdx[0] * x_jet.v).norm(), 0.0, kTolerance);
+  EXPECT_NEAR((f_jets[1].v - dfdx[1] * x_jet.v).norm(), 0.0, kTolerance);
 }
 
 class BiCubicInterpolatorTest : public ::testing::Test {
  public:
+  template <int kDataDimension>
   void RunPolynomialInterpolationTest(const Eigen::Matrix3d& coeff) {
+    values_.reset(new double[kNumRows * kNumCols * kDataDimension]);
     coeff_ = coeff;
-    double* v = values_;
+    double* v = values_.get();
     for (int r = 0; r < kNumRows; ++r) {
       for (int c = 0; c < kNumCols; ++c) {
-        *v++ = EvaluateF(r, c);
+        for (int dim = 0; dim < kDataDimension; ++dim) {
+          *v++ = (dim * dim + 1) * EvaluateF(r, c);
+        }
       }
     }
-    BiCubicInterpolator interpolator(values_, kNumRows, kNumCols);
+
+    Array2D<double, kDataDimension> array(values_.get(), kNumRows, kNumCols);
+    BiCubicInterpolator<Array2D<double, kDataDimension> > interpolator(array);
 
     for (int j = 0; j < kNumRowSamples; ++j) {
       const double r = 1.0 + 7.0 / (kNumRowSamples - 1) * j;
       for (int k = 0; k < kNumColSamples; ++k) {
         const double c = 1.0 + 7.0 / (kNumColSamples - 1) * k;
-        const double expected_f = EvaluateF(r, c);
-        const double expected_dfdr = EvaluatedFdr(r, c);
-        const double expected_dfdc = EvaluatedFdc(r, c);
-        double f, dfdr, dfdc;
-
-        EXPECT_TRUE(interpolator.Evaluate(r, c, &f, &dfdr, &dfdc));
-        EXPECT_NEAR(f, expected_f, kTolerance);
-        EXPECT_NEAR(dfdr, expected_dfdr, kTolerance);
-        EXPECT_NEAR(dfdc, expected_dfdc, kTolerance);
+        double f[kDataDimension], dfdr[kDataDimension], dfdc[kDataDimension];
+        EXPECT_TRUE(interpolator.Evaluate(r, c, f, dfdr, dfdc));
+        for (int dim = 0; dim < kDataDimension; ++dim) {
+          EXPECT_NEAR(f[dim], (dim * dim + 1) * EvaluateF(r, c), kTolerance);
+          EXPECT_NEAR(dfdr[dim], (dim * dim + 1) * EvaluatedFdr(r, c), kTolerance);
+          EXPECT_NEAR(dfdc[dim], (dim * dim + 1) * EvaluatedFdc(r, c), kTolerance);
+        }
       }
     }
   }
@@ -187,18 +326,22 @@
   static const int kNumCols = 10;
   static const int kNumRowSamples = 100;
   static const int kNumColSamples = 100;
-  double values_[kNumRows * kNumCols];
+  scoped_array<double> values_;
 };
 
 TEST_F(BiCubicInterpolatorTest, ZeroFunction) {
   Eigen::Matrix3d coeff = Eigen::Matrix3d::Zero();
-  RunPolynomialInterpolationTest(coeff);
+  RunPolynomialInterpolationTest<1>(coeff);
+  RunPolynomialInterpolationTest<2>(coeff);
+  RunPolynomialInterpolationTest<3>(coeff);
 }
 
 TEST_F(BiCubicInterpolatorTest, Degree00Function) {
   Eigen::Matrix3d coeff = Eigen::Matrix3d::Zero();
   coeff(2, 2) = 1.0;
-  RunPolynomialInterpolationTest(coeff);
+  RunPolynomialInterpolationTest<1>(coeff);
+  RunPolynomialInterpolationTest<2>(coeff);
+  RunPolynomialInterpolationTest<3>(coeff);
 }
 
 TEST_F(BiCubicInterpolatorTest, Degree01Function) {
@@ -206,7 +349,9 @@
   coeff(2, 2) = 1.0;
   coeff(0, 2) = 0.1;
   coeff(2, 0) = 0.1;
-  RunPolynomialInterpolationTest(coeff);
+  RunPolynomialInterpolationTest<1>(coeff);
+  RunPolynomialInterpolationTest<2>(coeff);
+  RunPolynomialInterpolationTest<3>(coeff);
 }
 
 TEST_F(BiCubicInterpolatorTest, Degree10Function) {
@@ -214,7 +359,9 @@
   coeff(2, 2) = 1.0;
   coeff(0, 1) = 0.1;
   coeff(1, 0) = 0.1;
-  RunPolynomialInterpolationTest(coeff);
+  RunPolynomialInterpolationTest<1>(coeff);
+  RunPolynomialInterpolationTest<2>(coeff);
+  RunPolynomialInterpolationTest<3>(coeff);
 }
 
 TEST_F(BiCubicInterpolatorTest, Degree11Function) {
@@ -224,7 +371,9 @@
   coeff(1, 0) = 0.1;
   coeff(0, 2) = 0.2;
   coeff(2, 0) = 0.2;
-  RunPolynomialInterpolationTest(coeff);
+  RunPolynomialInterpolationTest<1>(coeff);
+  RunPolynomialInterpolationTest<2>(coeff);
+  RunPolynomialInterpolationTest<3>(coeff);
 }
 
 TEST_F(BiCubicInterpolatorTest, Degree12Function) {
@@ -235,7 +384,9 @@
   coeff(0, 2) = 0.2;
   coeff(2, 0) = 0.2;
   coeff(1, 1) = 0.3;
-  RunPolynomialInterpolationTest(coeff);
+  RunPolynomialInterpolationTest<1>(coeff);
+  RunPolynomialInterpolationTest<2>(coeff);
+  RunPolynomialInterpolationTest<3>(coeff);
 }
 
 TEST_F(BiCubicInterpolatorTest, Degree21Function) {
@@ -246,7 +397,9 @@
   coeff(0, 2) = 0.2;
   coeff(2, 0) = 0.2;
   coeff(0, 0) = 0.3;
-  RunPolynomialInterpolationTest(coeff);
+  RunPolynomialInterpolationTest<1>(coeff);
+  RunPolynomialInterpolationTest<2>(coeff);
+  RunPolynomialInterpolationTest<3>(coeff);
 }
 
 TEST_F(BiCubicInterpolatorTest, Degree22Function) {
@@ -259,17 +412,22 @@
   coeff(0, 0) = 0.3;
   coeff(0, 1) = -0.4;
   coeff(1, 0) = -0.4;
-  RunPolynomialInterpolationTest(coeff);
+  RunPolynomialInterpolationTest<1>(coeff);
+  RunPolynomialInterpolationTest<2>(coeff);
+  RunPolynomialInterpolationTest<3>(coeff);
 }
 
 TEST(BiCubicInterpolator, JetEvaluation) {
-  const double values[] = {1.0, 2.0, 2.0, 3.0,
-                           1.0, 2.0, 2.0, 3.0};
-  BiCubicInterpolator interpolator(values, 2, 4);
-  double f, dfdr, dfdc;
+  const double values[] = {1.0, 5.0, 2.0, 10.0, 2.0, 6.0, 3.0, 5.0,
+                           1.0, 2.0, 2.0,  2.0, 2.0, 2.0, 3.0, 1.0};
+
+  Array2D<double, 2> array(values, 2, 4);
+  BiCubicInterpolator<Array2D<double, 2> > interpolator(array);
+
+  double f[2], dfdr[2], dfdc[2];
   const double r = 0.5;
   const double c = 2.5;
-  EXPECT_TRUE(interpolator.Evaluate(r, c, &f, &dfdr, &dfdc));
+  EXPECT_TRUE(interpolator.Evaluate(r, c, f, dfdr, dfdc));
 
   // Create a Jet with the same scalar part as x, so that the output
   // Jet will be evaluated at x.
@@ -287,10 +445,16 @@
   c_jet.v(2) = 4.2;
   c_jet.v(3) = 5.3;
 
-  Jet<double, 4> f_jet;
-  EXPECT_TRUE(interpolator.Evaluate(r_jet, c_jet, &f_jet));
-  EXPECT_EQ(f_jet.a, f);
-  EXPECT_EQ((f_jet.v - dfdr * r_jet.v - dfdc * c_jet.v).norm(), 0.0);
+  Jet<double, 4> f_jets[2];
+  EXPECT_TRUE(interpolator.Evaluate(r_jet, c_jet, f_jets));
+  EXPECT_EQ(f_jets[0].a, f[0]);
+  EXPECT_EQ(f_jets[1].a, f[1]);
+  EXPECT_NEAR((f_jets[0].v - dfdr[0] * r_jet.v - dfdc[0] * c_jet.v).norm(),
+              0.0,
+              kTolerance);
+  EXPECT_NEAR((f_jets[1].v - dfdr[1] * r_jet.v - dfdc[1] * c_jet.v).norm(),
+              0.0,
+              kTolerance);
 }
 
 }  // namespace internal