Bicubic Interpolation
Here pretty good explanation.

I'll start by consider the 1D case of cubic interpolation (because its easier to explain) and then go onto the 2D case.

The basic idea of cubic interpolation is to estimate the values between any two points by as a cubic function, $f(x)$, with first derivative, $f'(x)$

$$f(x) = ax^3 + bx^2 + cx + d$$
$$f'(x) = 3ax^2 + 2bx + c$$

so that, when spliced to the adjacent cubic splines on the left and on the right, the interpolated function (the "zeroth derivative") and its first derivative are continuous.  That means, at the left boundary of some given spline, the zeroth and first derivatives are equal to the those of the right boundary of the adjacent spline immediately to the left.

Without loss of generality, we can assume that we are interpolating between $x=0$ and $x=1$.  On the left and right boundaries, we see that

$$f(0) = d$$
$$f(1) = a+b+c+d$$
$$f'(0) = c$$
$$f'(1) = 3a+2b+c$$

Solving for the polynomial coefficients

$$a = 2f(0) - 2f(1) + f'(0) + f'(1)$$
$$b = -3f(0) + 3f(1) - 2f'(0) - f'(1)$$
$$c = f'(0)$$
$$d = f(0)$$

Now consider the four points around the region being estimated, which are $(-1,y_{-1})$, $(0,y_0)$, $(1,y_1)$, $(2,y_2)$, and impose boundary conditions that 

$$f(0) = y_0$$
$$f(1) = y_1$$

and define the two boundary derivatives as

$$f'(0) = \frac{y_1-y_{-1}}{2}$$
$$f'(1) = \frac{y_2-y_0}{2}$$

which is consistent to how they will be defined for the adjacent splines.

Then 

$$a = -\frac{y_{-1}}{2} + \frac{3y_0}{2} - \frac{3y_1}{2} + \frac{y_2}{2}$$
$$b = y_{-1} - \frac{5y_0}{2} + 2y_1 - \frac{y_2}{2}$$
$$c = -\frac{y_{-1}}{2} + \frac{y_1}{2}$$
$$d = y_0$$

once the coefficients are known you can interpolate $f(x)$ for any point between the central 2 points, $(0,y_0)$ and $(1,y_1)$.

For bicubic interpolation the principle is pretty much the same but you estimate a surface using 16 points (4x4 grid) rather than just a curve. One simple way to do this is to first interpolate the columns and then interpolate the resulting rows.