Terrain Height map resampling with BiCubic filtering

See the previous post for the TerrainStruct.h header file.  Somewhat redundant post but setup to compile on my linux compiler with the given gnu version that I am evidently using which for some odd reason is fussy around passing by values arrays (or something funny like this).

//*****************************BiCubicResample.cpp****************************************
#define __BicubicResample_CPP
#include <vector>
/*
double cubicInterpolate (std::vector<double> * p, double x) {
 return (*p)[1] + 0.5 * x*((*p)[2] - (*p)[0] + x*(2.0*(*p)[0] - 5.0*(*p)[1] + 4.0*(*p)[2] - (*p)[3] + x*(3.0*((*p)[1] - (*p)[2]) + (*p)[3] - (*p)[0])));
}
*/
///*
double cubicInterpolate (std::vector<double> * p, double x) {
 return (*p)[1] + .5f*x*(-1.0f*(*p)[0]+(*p)[2]) + x*x*((*p)[0] - 5.0f/2.0f * (*p)[1] +2.0f*(*p)[2] -.5f*(*p)[3]) + x*x*x*(-.5f*(*p)[0] + 1.5f*(*p)[1] - 1.5f*(*p)[2] +.5f*(*p)[3]);
}
//*/

//trying a forward differencing as opposed to central difference method. Same deal.
//
/*
double cubicInterpolate (std::vector<double> * p, double x) {
 return (*p)[1] + x*(-1.0f*(*p)[0]+(*p)[1]) + x*x*(2.0f*(*p)[0] - 5.0f * (*p)[1] +4.0f*(*p)[2] -1.0f*(*p)[3]) + x*x*x*(-1.0f*(*p)[0] + 3.0f*(*p)[1] - 3.0f*(*p)[2] +1.0f*(*p)[3]);
}
*/
double bicubicInterpolate (std::vector<std::vector<double> > * p, double x, double y) {
 std::vector<double>  arr(4,0);
 for (int i = 0; i < 4; i ++){
    std::vector<double> v(4,0);
    v[0] = (*p)[i][0]; v[1] = (*p)[i][1]; v[2] = (*p)[i][2]; v[3] = (*p)[i][3];
    /*
    arr[0] = cubicInterpolate(&((*p)[0]), y);
    arr[1] = cubicInterpolate(&((*p)[1]), y);
    arr[2] = cubicInterpolate(&((*p)[2]), y);
    arr[3] = cubicInterpolate(&((*p)[3]), y);
    */
    arr[i] = cubicInterpolate(&v,y);
 }
 return cubicInterpolate(&arr, x);
}
#endif 

//**************************BuildBiCubicResample.cpp*****************************************
#ifndef __BuildBicubicResample_CPP
#define __BuildBicubicResample_CPP
#include "TerrainStruct.h"
#include "BicubicResample.cpp"
#include <vector>

terr::CPointsMap * BuildBicubicResampl( double size, double RSize, terr::CPointsMap * heightmap){
 terr::CPointsMap * Rheightmap = new terr::CPointsMap();
  
 for (int i = 0; i < RSize; i++){
  for(int j = 0; j < RSize; j++){
   
   double x = ((double) i) * (size-2.0f)/(RSize-1.0f);
   double y = ((double) j) * (size-2.0f)/(RSize-1.0f); 
   
   int p1x = (int)x; int p1y = (int)y;
   double locx = x -p1x; double locy = y-p1y; 
   int p2x = (int)x+1; int p2y = (int)y+1;
   int p0x = (x <= 0 ? p1x : p1x-1); int p0y = (y <= 0 ? p1y : p1y-1);
   int p3x = (x >= (size-2.0f) ? p2x : p2x+1); int p3y = (y >= (size-2.0f) ? p2y : p2y+1); 
   terr::Coordpair * p00 = new terr::Coordpair(p0x,p0y);
   terr::Coordpair * p01 = new terr::Coordpair(p0x,p1y);
   terr::Coordpair * p02 = new terr::Coordpair(p0x,p2y);
   terr::Coordpair * p03 = new terr::Coordpair(p0x,p3y);
   terr::Coordpair * p10 = new terr::Coordpair(p1x,p0y);
   terr::Coordpair * p11 = new terr::Coordpair(p1x,p1y);
   terr::Coordpair * p12 = new terr::Coordpair(p1x,p2y);
   terr::Coordpair * p13 = new terr::Coordpair(p1x,p3y);
   terr::Coordpair * p20 = new terr::Coordpair(p2x,p0y);
   terr::Coordpair * p21 = new terr::Coordpair(p2x,p1y);
   terr::Coordpair * p22 = new terr::Coordpair(p2x,p2y);
   terr::Coordpair * p23 = new terr::Coordpair(p2x,p3y);
   terr::Coordpair * p30 = new terr::Coordpair(p3x,p0y);
   terr::Coordpair * p31 = new terr::Coordpair(p3x,p1y);
   terr::Coordpair * p32 = new terr::Coordpair(p3x,p2y);
   terr::Coordpair * p33 = new terr::Coordpair(p3x,p3y);
   
   std::vector<double> ay(4,0);
   std::vector<std::vector<double> >  a(4,ay);
   a[0][0] = (*heightmap)[(*p00)]; a[0][1] = (*heightmap)[(*p01)]; 
   a[0][2] = (*heightmap)[(*p02)]; a[0][3] = (*heightmap)[(*p03)]; 
   a[1][0] = (*heightmap)[(*p10)]; a[1][1] = (*heightmap)[(*p11)]; 
   a[1][2] = (*heightmap)[(*p12)]; a[1][3] = (*heightmap)[(*p13)];
   a[2][0] = (*heightmap)[(*p20)]; a[2][1] = (*heightmap)[(*p21)]; 
   a[2][2] = (*heightmap)[(*p22)]; a[2][3] = (*heightmap)[(*p23)];
   a[3][0] = (*heightmap)[(*p30)]; a[3][1] = (*heightmap)[(*p31)]; 
   a[3][2] = (*heightmap)[(*p32)]; a[3][3] = (*heightmap)[(*p33)];   
   
   double height = bicubicInterpolate (&a, locx, locy); //slower matrix compute method
   terr::Coordpair * coord = new terr::Coordpair(i,j);
   (*Rheightmap)[(*coord)] = height;
  }
 }
 return Rheightmap;
}
#endif