Add interpolation routines

@cdeil @keflavich - here are some initial reprojection routines based on my PR to astropy. This now includes a cube reprojection routine if the third axis is independent.

@keflavich - I was going to implement your bad pixel preservation logic from FITS_tools, but map_coordinates seems to do the right thing by default (preserve nans and infs). Do you think this is because scipy implemented it in the mean time? Otherwise, any feature from FITS_tools that I'm missing here?

This needs:
* [ ] Tests
* [ ] Docs
* [x] A full n-d reprojection function
Non-integers in the `.shape` tuple should be impossible - the `int(round(` operation is therefore unnecessary.  `mul` should also respect the types.  It might be worth instead `assert`ing that all the inputs are ints.