Need performance tips on this 4D LUT

[onqel]

Hi guys
This is my first post here so bare with me ;)

I've made a rather simple solution for a lookup-table class for up to 4 input variables of float/double types. There's only linear interpolation at the moment. I'm looking for tips to improve the performance on this one as it is a bit crude at the moment.. I'm thinking about learning SSE2 intrinsics in the future, is this a place where it could come to benefit so I might be able to cubic/spline interpolation instead? I'm not a professional programmer, just doing this as a hobby :)

Thanks for any feedback in advance!

/*
    X = input 1 (number of results / table size)
    Y = input 2 (number of 1D LUT's)
    Z = input 3 (number of 2D LUT's)
    W = input 4 (number of 3D LUT's)
*/

template <int X, typename T=double>
class LookupTable
{
public:
    LookupTable()
    {
        for (int i=0; i<X; i++)
        {
            val[i] = 2.0;
        }
        imax = X-1;
    }
    
    void setup(T xMax, T xMin)
    {
        mX = (T) X / (xMax - xMin);
        limH = xMax;
        limL = xMin;
        
        cLimH = xMax;
        cLimL = xMin;
    }
    
    T *getTable() { return val; }
    
    inline void setX(uint &x, T &value)
    {
        val[x] = value;
    }
    
    inline const T get(T &x_val)
    {
// Keep value within limits by clipping
        if (x_val >= (limH)) x_val = limH;
        else if (x_val <= limL) x_val = limL;
        
        const T x0 = (limH - x_val) * mX;
        const uint row = (uint) x0;
        const T x = x0 - (T) row;
        const T y0 = val[row] * (1.0-x);
        const T y1 = val[row+1] * x;
        return y0+y1;
    }
    
    
private:
    uint imax;
    T mX, limH, limL, cLimH, cLimL;
    T val[X];
};


template <int X=0, int Y=0, typename T=double>
class LookupTable2D
{
public:
    LookupTable2D()
    {

    }
    
    void setup(T xMax, T xMin, T yMax, T yMin)
    {
        mY = (T) Y / (yMax - yMin);
        limH = yMax;
        limL = yMin;
        
        for (int i=0; i<Y; i++)
        {
            lut[i].setup(xMax, xMin);
        }
    }
    
    inline void setXY(uint &x, uint &y, T &value)
    {
        lut[y].setX(x, value);
    }
    
    inline const T get(T &x_val, T &y_val)
    {
        if (y_val >= limH) y_val = limH;
        else if (y_val <= limL) y_val = limL;
        
        const T x0 = (limH - y_val) * mY;
        const uint row = (uint) x0;
        const T x = x0 - (double) row;
        const T y0 = lut[row].get(x_val) * (1.0-x);
        const T y1 = lut[row+1].get(x_val) * x;
        return y0+y1;
    }
    
private:
    T mY, limH, limL;
    
    LookupTable<X, T> lut[Y];
};


template <int X, int Y=0, int Z=0, typename T=double>
class LookupTable3D
{
public:
    
    LookupTable3D()
    {

    }
    
    void setup(T xMax, T xMin, T yMax, T yMin, T zMax, T zMin)
    {
        mZ = (T) Z / (zMax - zMin);
        limH = zMax;
        limL = zMin;
        
        for (int i=0; i<Z; i++)
        {
            lut2D[i].setup(xMax, xMin, yMax, yMin);
        }
    }
    
    
    inline void setXYZ(uint &x, uint &y, uint &z, T &value)
    {
        lut2D[z].setXY(x, y, value);
    }
    
    inline const T get(T &x_val, T &y_val, T &z_val)
    {
        if (z_val >= limH) z_val = limH;
        else if (z_val <= limL) z_val = limL;
        
        const T x0 = (limH - z_val) * mZ;
        const uint row = x0;
        const T x = x0 - (T) row;
        const T y0 = lut2D[row].get(x_val, y_val) * (1.0-x);
        const T y1 = lut2D[row+1].get(x_val, y_val) * x;
        return y0+y1;
    }
    
    
private:
    T mZ;
    T limH, limL;
    
    LookupTable2D<X, Y, T> lut2D[Z];
};


template <int X, int Y=0, int Z=0, int W=0, typename T=double>
class LookupTable4D
{
public:
    LookupTable4D()
    {
        
    }
    
    void setup(T xMax, T xMin, T yMax, T yMin, T zMax, T zMin, T wMax, T wMin)
    {
        mW = (T) W / (wMax - wMin);
        limH = wMax;
        limL = wMin;
        
        for (int i=0; i<W; i++)
        {
            lut3D[i].setup(xMax, xMin, yMax, yMin, zMax, zMin);
        }
    }
    
    inline void setXYZW(uint &x, uint &y, uint &z, uint &w, T value)
    {
        lut3D[w].setXYZ(x, y, z, value);
    }
    
    inline const T get(T &x_val, T &y_val, T &z_val, T &w_val)
    {
        if (w_val >= limH) w_val = limH;
        else if (w_val <= limL) w_val = limL;
        
        const T x0 = (limH - w_val) * mW;
        const uint row = x0;
        const T x = x0 - (T) row;
        const T y0 = lut3D[row].get(x_val, y_val, z_val) * (1.0-x);
        const T y1 = lut3D[row+1].get(x_val, y_val, z_val) * x;
        return y0+y1;
    }
    
    
private:
    T mW;
    T limH, limL;
    
    LookupTable3D<X, Y, Z, T> lut3D[W];
};



using namespace std;

int main(int argc, const char * argv[])
{
// Total size = X*Y*Z*W
    const uint X = 60;
    const uint Y = 20;
    const uint Z = 20;
    const uint W = 20;
    
    LookupTable4D<X, Y, Z, W, float> lut4D;

    double xmin = -10.0;
    double xmax = 10.25;
    double x_step = (xmax - xmin) / (double) X;
    
    double ymin = -10.0;
    double ymax = 10.14;
    double y_step = (ymax - ymin) / (double) Y;
    
    double zmin = -10.0;
    double zmax = 10.25;
    double z_step = (zmax - zmin) / (double) Z;
    
    double wmin = -10.0;
    double wmax = 10.0;
    double w_step = (wmax - wmin) / (double) W;
    
    double x = xmax;
    double y = ymax;
    double z = zmax;
    double w = wmax;
    
    lut4D.setup(xmax, xmin, ymax, ymin, zmax, zmin, wmax, wmin);
    
    uint i, j, k, l;
   
// Fill LUT with test values
    for (i=0; i<W; i++)
    {
        z = zmax;
        for (j=0; j<Z; j++)
        {
            y = ymax;

            for (k=0; k<Y; k++)
            {
                x = xmax;
                for (l=0; l<X; l++)
                {
                    float val = x/2.0+y*3.0+z-1.0+w*w;
                    lut4D.setXYZW(l, k, j, i, val);
                    // X decrement
                    x -= x_step;
                }
                // Y decrement
                y -= y_step;
            }
            // Z decrement
            z -= z_step;
        }
        // W decrement
        w -= w_step;
    }
    
    x = 3.3;
    y = 5.5;
    z = 1.562;
    w = 3.45;
    
    cout.precision(12);
    cout << x/2.0+y*3.0+z-1.0+w*w << endl;
    
    float a = x;
    float b = y;
    float c = z;
    float d = w;
    
    float tmp = lut4D.get(a, b, c, d);
    cout << tmp << endl;

    return 0;
}

Sorry if the code is too messy !