diff --git a/docs/source/bibliography.rst b/docs/source/bibliography.rst
index bd99758..bab6e61 100644
--- a/docs/source/bibliography.rst
+++ b/docs/source/bibliography.rst
@@ -44,6 +44,10 @@
    with indeterminacy**, *IEEE Transactions on Information Theory*
    47(5):2017-2028, 2001.
 
+.. [Keys] R. G. Keys, **Cubic convolution interpolation for digital
+   image processing**, *IEEE Trans. on Acoustics, Speech, and Signal
+   Processing*, 29(6), 1981.
+
 .. [KushalAgarwal] A. Kushal and S. Agarwal, **Visibility based
    preconditioning for bundle adjustment**, *In Proceedings of the
    IEEE Conference on Computer Vision and Pattern Recognition*, 2012.
diff --git a/docs/source/nnls_modeling.rst b/docs/source/nnls_modeling.rst
index 37f7c59..f9e9de6 100644
--- a/docs/source/nnls_modeling.rst
+++ b/docs/source/nnls_modeling.rst
@@ -67,8 +67,8 @@
 :math:`f_i\left(x_{i_1},...,x_{i_k}\right)` and the matrices
 
  .. math:: J_{ij} = \frac{\partial}{\partial
-	   x_{i_j}}f_i\left(x_{i_1},...,x_{i_k}\right),\quad \forall j
-	   \in \{1, \ldots, k\}
+           x_{i_j}}f_i\left(x_{i_1},...,x_{i_k}\right),\quad \forall j
+           \in \{1, \ldots, k\}
 
 .. class:: CostFunction
 
@@ -495,6 +495,39 @@
    argument. Please be careful when setting the size parameters.
 
 
+Numeric Differentiation & LocalParameterization
+-----------------------------------------------
+
+   If your cost function depends on a parameter block that must lie on
+   a manifold and the functor cannot be evaluated for values of that
+   parameter block not on the manifold then you may have problems
+   numerically differentiating such functors.
+
+   This is because numeric differentiation in Ceres is performed by
+   perturbing the individual coordinates of the parameter blocks that
+   a cost functor depends on. In doing so, we assume that the
+   parameter blocks live in an Euclidean space and ignore the
+   structure of manifold that they live As a result some of the
+   perturbations may not lie on the manifold corresponding to the
+   parameter block.
+
+   For example consider a four dimensional parameter block that is
+   interpreted as a unit Quaternion. Perturbing the coordinates of
+   this parameter block will violate the unit norm property of the
+   parameter block.
+
+   Fixing this problem requires that :class:`NumericDiffCostFunction`
+   be aware of the :class:`LocalParameterization` associated with each
+   parameter block and only generate perturbations in the local
+   tangent space of each parameter block.
+
+   For now this is not considered to be a serious enough problem to
+   warrant changing the :class:`NumericDiffCostFunction` API. Further,
+   in most cases it is relatively straightforward to project a point
+   off the manifold back onto the manifold before using it in the
+   functor. For example in case of the Quaternion, normalizing the
+   4-vector before using it does the trick.
+
    **Alternate Interface**
 
    For a variety of reason, including compatibility with legacy code,
@@ -571,6 +604,10 @@
    As a rule of thumb, try using :class:`NumericDiffCostFunction` before
    you use :class:`DynamicNumericDiffCostFunction`.
 
+   **WARNING** The same caution about mixing local parameterizations
+   with numeric differentiation applies as is the case with
+   :class:`NumericDiffCostFunction``.
+
 :class:`CostFunctionToFunctor`
 ==============================
 
@@ -1686,3 +1723,121 @@
 .. function:: void AngleAxisRotatePoint<T>(const T angle_axis[3], const T pt[3], T result[3])
 
    .. math:: y = R(\text{angle_axis}) x
+
+
+Cubic Interpolation
+===================
+
+Optimization problems often involve functions that are given in the
+form of a table of values, for example an image. Evaluating these
+functions and their derivatives requires interpolating these
+values. Interpolating tabulated functions is a vast area of research
+and there are a lot of libraries which implement a variety of
+interpolation schemes. However, using them within the automatic
+differentiation framework in Ceres is quite painful. To this end,
+Ceres provides the ability to interpolate one dimensional and two
+dimensional tabular functions.
+
+The one dimensional interpolation is based on the Cubic Hermite
+Spline, also known as the Catmull-Rom Spline. This produces a first
+order differentiable interpolating function. The two dimensional
+interpolation scheme is a generalization of the one dimensional scheme
+where the interpolating function is assumed to be separable in the two
+dimensions,
+
+More details of the construction can be found `Linear Methods for
+Image Interpolation <http://www.ipol.im/pub/art/2011/g_lmii/>`_ by
+Pascal Getreuer.
+
+.. class:: CubicInterpolator
+
+Given as input a one dimensional array like object, which provides
+the following interface.
+
+.. code::
+
+  struct Array {
+    enum { DATA_DIMENSION = 2; };
+    void GetValue(int n, double* f) const;
+    int NumValues() const;
+  };
+
+Where, ``GetValue`` gives us the value of a function :math:`f`
+(possibly vector valued) on the integers :code:`{0, ..., NumValues() -
+1}` and the enum ``DATA_DIMENSION`` indicates the dimensionality of
+the function being interpolated. For example if you are interpolating
+a color image with three channels (Red, Green & Blue), then
+``DATA_DIMENSION = 3``.
+
+:class:`CubicInterpolator` uses Cubic Hermite splines to produce a
+smooth approximation to it that can be used to evaluate the
+:math:`f(x)` and :math:`f'(x)` at any real valued point in the
+interval :code:`[0, NumValues() - 1]`. For example, the following code
+interpolates an array of four numbers.
+
+.. code::
+
+  const double data[] = {1.0, 2.0, 5.0, 6.0};
+  Array1D<double, 1> array(x, 4);
+  CubicInterpolator interpolator(array);
+  double f, dfdx;
+  CHECK(interpolator.Evaluate(1.5, &f, &dfdx));
+
+
+In the above code we use ``Array1D`` a templated helper class that
+allows easy interfacing between ``C++`` arrays and
+:class:`CubicInterpolator`.
+
+``Array1D`` supports vector valued functions where the various
+coordinates of the function can be interleaved or stacked. It also
+allows the use of any numeric type as input, as long as it can be
+safely cast to a double.
+
+.. class:: BiCubicInterpolator
+
+Given as input a two dimensional array like object, which provides
+the following interface:
+
+.. code::
+
+  struct Array {
+    enum { DATA_DIMENSION = 1 };
+    void GetValue(int row, int col, double* f) const;
+    int NumRows() const;
+    int NumCols() const;
+  };
+
+Where, ``GetValue`` gives us the value of a function :math:`f`
+(possibly vector valued) on the integer grid :code:`{0, ...,
+NumRows() - 1} x {0, ..., NumCols() - 1}` and the enum
+``DATA_DIMENSION`` indicates the dimensionality of the function being
+interpolated. For example if you are interpolating a color image with
+three channels (Red, Green & Blue), then ``DATA_DIMENSION = 3``.
+
+:class:`BiCubicInterpolator` uses the cubic convolution interpolation
+algorithm of R. Keys [Keys]_, to produce a smooth approximation to it
+that can be used to evaluate the :math:`f(r,c)`, :math:`\frac{\partial
+f(r,c)}{\partial r}` and :math:`\frac{\partial f(r,c)}{\partial c}` at
+any real valued point in the quad :code:`[0, NumRows() - 1] x [0,
+NumCols() - 1]`.
+
+For example the following code interpolates a two dimensional array.
+
+.. code::
+
+   const double data[] = {1.0, 3.0, -1.0, 4.0,
+                          3.6, 2.1,  4.2, 2.0,
+                          2.0, 1.0,  3.1, 5.2};
+   Array2D<double, 1>  array(data, 3, 4);
+   BiCubicInterpolator interpolator(array);
+   double f, dfdr, dfdc;
+   CHECK(interpolator.Evaluate(1.2, 2.5, &f, &dfdr, &dfdc));
+
+In the above code, the templated helper class ``Array2D`` is used to
+make a ``C++`` array look like a two dimensional table to
+:class:`BiCubicInterpolator`.
+
+``Array2D`` supports row or column major layouts. It also supports
+vector valued functions where the individual coordinates of the
+function may be interleaved or stacked. It also allows the use of any
+numeric type as input, as long as it can be safely cast to double.
diff --git a/docs/source/users.rst b/docs/source/users.rst
index d04927d..20f0e18 100644
--- a/docs/source/users.rst
+++ b/docs/source/users.rst
@@ -48,29 +48,3 @@
 .. _calibrating robot-camera systems:
    http://rosindustrial.org/news/2014/9/24/industrial-calibration-library-update-and-presentation
 .. _skinned control meshes: http://research.microsoft.com/en-us/projects/handmodelingfrommonoculardepth/
-
-
-Publications
-============
-
-Ceres Solver is used (and cited) in the following publications:
-
-#. **User-Specific Hand Modeling from Monocular Depth
-   Sequences**, J. Taylor, R. Stebbing, V. Ramakrishna, C. Keskin, J. Shotton, S. Izadi, A. Hertzmann,
-   and A. Fitzgibbon, CVPR 2014.
-
-#. **Global Fusion of Relative Motions for Robust, Accurate and
-   Scalable Structure from Motion**, P. Moulon, P. Monasse
-   and R. Marlet, ICCV 2013.
-
-#. **Recurrent neural networks for voice activity
-   detection**, T. Hughes and K. Mierle, ICASSP 2013.
-
-#. **Street View Motion-from-Structure-from-Motion**, B. Klingner, D. Martin
-   and J. Roseborough, ICCV 2013.
-
-#. **Adaptive Structure from Motion with a contrario model
-   estimation**, P. Moulon, P. Monasse and R. Marlet, ACCV 2012.
-
-#. **Visibility based preconditioning for bundle
-   adjustment**, A. Kushal and S. Agarwal, CVPR 2012.
diff --git a/docs/source/version_history.rst b/docs/source/version_history.rst
index f666fb5..f9e9c7b 100644
--- a/docs/source/version_history.rst
+++ b/docs/source/version_history.rst
@@ -12,15 +12,43 @@
 #. Added ``CubicInterpolator`` and ``BiCubicInterpolator`` to allow
    smooth interpolation of sampled functions and integration with
    automatic differentiation.
+#. Add method to return covariance in tangent space. (Michael Vitus &
+   Steve Hsu)
+#. Add Homogeneous vector parameterization. (Michael Vitus)
 
-Bug Fixes
----------
+Bug Fixes & Minor Changes
+-------------------------
+#. Rename macros in the public API to have a ``CERES_`` prefix.
+#. Fix ``OrderedGroup::Reverse()`` when it is empty.
+#. Update the code to point to ceres-solver.org.
+#. Update documentation to point to the GitHub issue tracker.
+#. Disable ``LAPACK`` for iOS builds. (Greg Coombe)
+#. Force use of single-thread in ``Problem::Evaluate()`` without OpenMP.
+#. Less strict check for multithreading. (Chris Sweeney)
+#. Update tolerances in small_blas_test.cc (Philipp Hubner)
+#. Documentation corrections (Steve Hsu)
+#. Fixed ``sampled_function.cc`` (Pablo Speciale)
+#. Fix example code in the documentation. (Rodney Hoskinson)
+#. Improve the error handling in Conjugate Gradients.
+#. Improve preconditioner documentation.
+#. Remove dead code from fpclassify.h.
+#. Make Android.mk threads sensitive.
+#. Changed the ``CURRENT_CONFIG_INSTALL_DIR`` to be a variable local
+   to Ceres. (Chris Sweeney)
+#. Fix typo in the comments in ``Jet.h``. (Julius Ziegler)
+#. Add the ASL at ETH Zurich to the list of users.
+#. Fixed a typo in the documentation. (Richard Stebbing)
+#. Fixed a boundary handling bug in the BiCubic interpolation
+   code. (Bernhard Zeisl)
+#. Fixed a ``MSVC`` compilation bug in the cubic interpolation code
+   (Johannes Schönberger)
 #. Add covariance related files to the Android build.
 #. Update Ubuntu 14.04 installation instructions. (Filippo Basso)
 #. Improved logging for linear solver failures.
 #. Improved crash messages in ``Problem``.
 #. Hide Homebrew related variables in CMake GUI.
-#. Add SuiteSparse link dependency for compressed_col_sparse_matrix_utils_test.
+#. Add SuiteSparse link dependency for
+   compressed_col_sparse_matrix_utils_test.
 #. Autodetect Homebrew install prefix on OSX.
 #. Lint changes from William Rucklidge and Jim Roseborough.
 #. Remove ``using namespace std:`` from ``port.h``
diff --git a/include/ceres/cubic_interpolation.h b/include/ceres/cubic_interpolation.h
index 5f3e183..9e63d01 100644
--- a/include/ceres/cubic_interpolation.h
+++ b/include/ceres/cubic_interpolation.h
@@ -122,9 +122,9 @@
 //
 // Example usage:
 //
-//  const double x[] = {1.0, 2.0, 5.0, 6.0};
-//  Array1D data(x, 4);
-//  CubicInterpolator interpolator(data);
+//  const double data[] = {1.0, 2.0, 5.0, 6.0};
+//  Array1D<double, 1> array(x, 4);
+//  CubicInterpolator<Array1D<double, 1> > interpolator(array);
 //  double f, dfdx;
 //  CHECK(interpolator.Evaluator(1.5, &f, &dfdx));
 template<typename Array>
@@ -243,6 +243,17 @@
 //
 // http://en.wikipedia.org/wiki/Cubic_Hermite_spline
 // http://en.wikipedia.org/wiki/Bicubic_interpolation
+//
+// Example usage:
+//
+// const double data[] = {1.0, 3.0, -1.0, 4.0,
+//                         3.6, 2.1,  4.2, 2.0,
+//                        2.0, 1.0,  3.1, 5.2};
+//  Array2D<double, 1>  array(data, 3, 4);
+//  BiCubicInterpolator<Array2D<double, 1> > interpolator(array);
+//  double f, dfdr, dfdc;
+//  CHECK(interpolator.Evaluate(1.2, 2.5, &f, &dfdr, &dfdc));
+
 template<typename Array>
 class CERES_EXPORT BiCubicInterpolator {
  public: