#include <float.h>
#include <tcl.h>
#include <Unirect.h>
#include "RpField1D.h"
#include "RpFieldRect3D.h"

extern int GetFloatFromObj(Tcl_Interp *interp, Tcl_Obj *objPtr, 
	float *valuePtr);
extern int GetAxisFromObj(Tcl_Interp *interp, Tcl_Obj *objPtr, int *indexPtr);

static INLINE char *
skipspaces(char *string) 
{
    while (isspace(*string)) {
	string++;
    }
    return string;
}

static INLINE char *
getline(char **stringPtr, char *endPtr) 
{
    char *line, *p;

    line = skipspaces(*stringPtr);
    for (p = line; p < endPtr; p++) {
	if (*p == '\n') {
	    *p++ = '\0';
	    *stringPtr = p;
	    return line;
	}
    }
    *stringPtr = p;
    return line;
}

int 
Rappture::Unirect3d::LoadData(Tcl_Interp *interp, int objc, 
			      Tcl_Obj *const *objv)
{
    int num[3], nValues;
    float min[3], max[3];
    float *values;
    const char *units[4], *order;


    if ((objc & 0x01) == 0) {
        Tcl_AppendResult(interp, Tcl_GetString(objv[0]), ": ",
                "wrong number of arguments: should be key-value pairs",
                (char *)NULL);
        return TCL_ERROR;
    }

    /* Default order is  z, y, x. */
    int axis1, axis2, axis3;
    axis1 = 2; 			/* Z-axis */
    axis2 = 1; 			/* Y-axis */
    axis3 = 0;			/* X-axis */

    values = NULL;
    num[0] = num[1] = num[2] = nValues = 0;
    min[0] = min[1] = min[2] = max[0] = max[1] = max[2] = 0.0f;
    order = units[0] = units[1] = units[2] = units[3] = NULL;

    int i;
    for (i = 1; i < objc; i += 2) {
        const char *string;
        char c;

        string = Tcl_GetString(objv[i]);
        c = string[0];
        if ((c == 'x') && (strcmp(string, "xmin") == 0)) {
            if (GetFloatFromObj(interp, objv[i+1], min+2) != TCL_OK) {
                return TCL_ERROR;
            }
        } else if ((c == 'x') && (strcmp(string, "xmax") == 0)) {
            if (GetFloatFromObj(interp, objv[i+1], max+2) != TCL_OK) {
                return TCL_ERROR;
            }
        } else if ((c == 'x') && (strcmp(string, "xnum") == 0)) {
            if (Tcl_GetIntFromObj(interp, objv[i+1], num+2) != TCL_OK) {
                return TCL_ERROR;
            }
            if (num[2] <= 0) {
                Tcl_AppendResult(interp, "bad xnum value: must be > 0",
                     (char *)NULL);
                return TCL_ERROR;
            }
        } else if ((c == 'x') && (strcmp(string, "xunits") == 0)) {
            units[2] = Tcl_GetString(objv[i+1]);
        } else if ((c == 'y') && (strcmp(string, "ymin") == 0)) {
            if (GetFloatFromObj(interp, objv[i+1], min+1) != TCL_OK) {
                return TCL_ERROR;
            }
        } else if ((c == 'y') && (strcmp(string, "ymax") == 0)) {
            if (GetFloatFromObj(interp, objv[i+1], max+1) != TCL_OK) {
                return TCL_ERROR;
            }
        } else if ((c == 'y') && (strcmp(string, "ynum") == 0)) {
            if (Tcl_GetIntFromObj(interp, objv[i+1], num+1) != TCL_OK) {
                return TCL_ERROR;
            }
            if (num[1] <= 0) {
                Tcl_AppendResult(interp, "bad ynum value: must be > 0",
                                 (char *)NULL);
                return TCL_ERROR;
            }
        } else if ((c == 'y') && (strcmp(string, "yunits") == 0)) {
            units[1] = Tcl_GetString(objv[i+1]);
        } else if ((c == 'z') && (strcmp(string, "zmin") == 0)) {
            if (GetFloatFromObj(interp, objv[i+1], min) != TCL_OK) {
                return TCL_ERROR;
            }
        } else if ((c == 'z') && (strcmp(string, "zmax") == 0)) {
            if (GetFloatFromObj(interp, objv[i+1], max) != TCL_OK) {
                return TCL_ERROR;
            }
        } else if ((c == 'z') && (strcmp(string, "znum") == 0)) {
            if (Tcl_GetIntFromObj(interp, objv[i+1], num) != TCL_OK) {
                return TCL_ERROR;
            }
            if (num[0] <= 0) {
                Tcl_AppendResult(interp, "bad znum value: must be > 0",
                                 (char *)NULL);
                return TCL_ERROR;
            }
        } else if ((c == 'z') && (strcmp(string, "zunits") == 0)) {
            units[0] = Tcl_GetString(objv[i+1]);
        } else if ((c == 'v') && (strcmp(string, "values") == 0)) {
            Tcl_Obj **vobj;

            if (Tcl_ListObjGetElements(interp, objv[i+1], &nValues, &vobj)
		!= TCL_OK) {
                return TCL_ERROR;
            }
            values = new float[nValues];
            int j;
            for (j = 0; j < nValues; j++) {
                if (GetFloatFromObj(interp, vobj[j], values + j)!=TCL_OK) {
                    return TCL_ERROR;
		}
	    }
        } else if ((c == 'u') && (strcmp(string, "units") == 0)) {
            _vUnits = strdup(Tcl_GetString(objv[i+1]));
        } else if ((c == 'o') && (strcmp(string, "order") == 0)) {
	    Tcl_Obj **axes;
	    int n;

	    if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &axes) 
		!= TCL_OK) {
                return TCL_ERROR;
            }
	    if (n != 3) {
                return TCL_ERROR;
	    }
	    if ((GetAxisFromObj(interp, axes[0], &axis1) != TCL_OK) ||
		(GetAxisFromObj(interp, axes[1], &axis2) != TCL_OK) ||
		(GetAxisFromObj(interp, axes[2], &axis3) != TCL_OK)) {
		return TCL_ERROR;
	    }
        } else {
            Tcl_AppendResult(interp, "unknown key \"", string,
                (char *)NULL);
            return TCL_ERROR;
        }
    }
    if (values == NULL) {
        Tcl_AppendResult(interp, "missing \"values\" key", (char *)NULL);
        return TCL_ERROR;
    }
    if ((size_t)nValues != (num[0] * num[1] * num[2] * _nComponents)) {
        Tcl_AppendResult(interp, 
		"wrong number of values: must be xnum*ynum*znum*extents", 
			 (char *)NULL);
        return TCL_ERROR;
    }
    
    if ((axis1 != 2) || (axis2 != 1) || (axis3 != 0)) {
	// Reorder the data into x, y, z where x varies fastest and so on.
	int z;
	float *data, *dp;

	dp = data = new float[nValues];
	for (z = 0; z < num[0]; z++) {
	    int y;

	    for (y = 0; y < num[1]; y++) {
		int x;

		for (x = 0; x < num[2]; x++) {
		    int i;
		    
		    /* Compute the index from the data's described ordering. */
		    i = ((z*num[axis2]*num[axis3]) + (y*num[axis3]) + x) * 3;
		    for(size_t v = 0; v < _nComponents; v++) {
			dp[v] = values[i+v];
		    }
		    dp += _nComponents;
		}
	    }
	}
	delete [] values;
	values = data;
    }

    _values = values;
    _nValues = nValues;
    if (units[3] != NULL) {
	_vUnits = strdup(units[3]);
    }
    _xMin = min[axis3];
    _xMax = max[axis3];
    _xNum = num[axis3];
    if (units[axis3] != NULL) {
	_xUnits = strdup(units[axis3]);
    }
    _yMin = min[axis2];
    _yMax = max[axis2];
    _yNum = num[axis2];
    if (units[axis2] != NULL) {
	_yUnits = strdup(units[axis2]);
    }
    _zMin = min[axis1];
    _zMax = max[axis1];
    _zNum = num[axis1];
    if (units[axis1] != NULL) {
	_zUnits = strdup(units[axis1]);
    }
    _initialized = true;
    return TCL_OK;
}


int 
Rappture::Unirect2d::LoadData(Tcl_Interp *interp, int objc, 
			      Tcl_Obj *const *objv)
{

    if ((objc & 0x01) == 0) {
        Tcl_AppendResult(interp, Tcl_GetString(objv[0]), ": ",
                "wrong number of arguments: should be key-value pairs",
                (char *)NULL);
        return TCL_ERROR;
    }

    int axis[2];
    axis[0] = 1; 			/* X-axis */
    axis[1] = 0;			/* Y-axis */

    _xNum = _yNum = _nValues = 0;
    _xMin = _yMin = _xMax = _yMax = 0.0f;
    if (_xUnits != NULL) {
	free(_xUnits);
    }
    if (_yUnits != NULL) {
	free(_yUnits);
    }
    if (_vUnits != NULL) {
	free(_vUnits);
    }
    _xUnits = _yUnits = _vUnits = NULL;
    if (_values != NULL) {
	delete [] _values;
    }
    _values = NULL;

    int i;
    for (i = 1; i < objc; i += 2) {
        const char *string;
        char c;

        string = Tcl_GetString(objv[i]);
        c = string[0];
        if ((c == 'x') && (strcmp(string, "xmin") == 0)) {
            if (GetFloatFromObj(interp, objv[i+1], &_xMin) != TCL_OK) {
                return TCL_ERROR;
            }
        } else if ((c == 'x') && (strcmp(string, "xmax") == 0)) {
            if (GetFloatFromObj(interp, objv[i+1], &_xMax) != TCL_OK) {
                return TCL_ERROR;
            }
        } else if ((c == 'x') && (strcmp(string, "xnum") == 0)) {
	    int n;
            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
                return TCL_ERROR;
            }
            if (n <= 0) {
                Tcl_AppendResult(interp, "bad xnum value: must be > 0",
                     (char *)NULL);
                return TCL_ERROR;
            }
	    _xNum = n;
        } else if ((c == 'x') && (strcmp(string, "xunits") == 0)) {
            _xUnits = strdup(Tcl_GetString(objv[i+1]));
        } else if ((c == 'y') && (strcmp(string, "ymin") == 0)) {
            if (GetFloatFromObj(interp, objv[i+1], &_yMin) != TCL_OK) {
                return TCL_ERROR;
            }
        } else if ((c == 'y') && (strcmp(string, "ymax") == 0)) {
            if (GetFloatFromObj(interp, objv[i+1], &_yMax) != TCL_OK) {
                return TCL_ERROR;
            }
        } else if ((c == 'y') && (strcmp(string, "ynum") == 0)) {
	    int n;
            if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
                return TCL_ERROR;
            }
            if (n <= 0) {
                Tcl_AppendResult(interp, "bad ynum value: must be > 0",
                                 (char *)NULL);
                return TCL_ERROR;
            }
	    _yNum = n;
        } else if ((c == 'y') && (strcmp(string, "yunits") == 0)) {
            _yUnits = strdup(Tcl_GetString(objv[i+1]));
        } else if ((c == 'v') && (strcmp(string, "values") == 0)) {
            Tcl_Obj **vobj;
	    int n;

            if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &vobj) != TCL_OK){
                return TCL_ERROR;
            }
	    if (n <= 0) {
                Tcl_AppendResult(interp, "empty values list : must be > 0",
                                 (char *)NULL);
                return TCL_ERROR;
	    }
	    _nValues = n;
            _values = new float[_nValues];
            size_t j;
            for (j = 0; j < _nValues; j++) {
                if (GetFloatFromObj(interp, vobj[j], _values + j)!=TCL_OK) {
                    return TCL_ERROR;
		}
	    }
        } else if ((c == 'u') && (strcmp(string, "units") == 0)) {
            _vUnits = strdup(Tcl_GetString(objv[i+1]));
        } else if ((c == 'e') && (strcmp(string, "extents") == 0)) {
	    int n;

	    if (Tcl_GetIntFromObj(interp, objv[i+1], &n) != TCL_OK) {
                return TCL_ERROR;
            }
	    if (n <= 0) {
                Tcl_AppendResult(interp, "bad extents value: must be > 0",
                                 (char *)NULL);
                return TCL_ERROR;
            }
	    _nComponents = n;
        } else if ((c == 'a') && (strcmp(string, "axisorder") == 0)) {
	    Tcl_Obj **order;
	    int n;

	    if (Tcl_ListObjGetElements(interp, objv[i+1], &n, &order) 
		!= TCL_OK) {
                return TCL_ERROR;
            }
	    if (n != 2) {
		Tcl_AppendResult(interp, 
			"wrong # of axes defined for ordering the data",
			(char *)NULL);
                return TCL_ERROR;
	    }
	    if ((GetAxisFromObj(interp, order[0], axis) != TCL_OK) ||
		(GetAxisFromObj(interp, order[1], axis+1) != TCL_OK)) {
		return TCL_ERROR;
	    }
        } else {
            Tcl_AppendResult(interp, "unknown key \"", string,
                (char *)NULL);
            return TCL_ERROR;
        }
    }
    if (_values == NULL) {
        Tcl_AppendResult(interp, "missing \"values\" key", (char *)NULL);
        return TCL_ERROR;
    }
    if (_nValues != (_xNum * _yNum * _nComponents)) {
        Tcl_AppendResult(interp, 
		"wrong number of values: must be xnum*ynum*extents", 
			 (char *)NULL);
        return TCL_ERROR;
    }
    
#ifndef notdef
    if ((axis[0] != 1) || (axis[1] != 0)) {
	fprintf(stderr, "reordering data\n");
	// Reorder the data into x, y where x varies fastest and so on.
	size_t y;
	float *data, *dp;

	dp = data = new float[_nValues];
	for (y = 0; y < _yNum; y++) {
	    size_t x;

	    for (x = 0; x < _xNum; x++) {
		size_t i, v;
		    
		/* Compute the index from the data's described ordering. */
		i = (y + (_yNum * x)) * _nComponents;
		for(v = 0; v < _nComponents; v++) {
		    dp[v] = _values[i+v];
		}
		dp += _nComponents;
	    }
	}
	delete [] _values;
	_values = data;
    }
#endif
    _initialized = true;
    return TCL_OK;
}


bool
Rappture::Unirect3d::ImportDx(Rappture::Outcome &result, int nComponents,
	size_t length, char *string) 
{
    size_t nx, ny, nz, npts;
    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
    char *p, *endPtr;

    dx = dy = dz = 0.0;         // Suppress compiler warning.
    x0 = y0 = z0 = 0.0;		// May not have an origin line.
    for (p = string, endPtr = p + length; p < endPtr; /*empty*/) {
	char *line;

	line = getline(&p, endPtr);
        if (line == endPtr) {
	    break;
	}
        if ((line[0] == '#') || (line == '\0')) {
	    continue;		// Skip blank or comment lines.
	}
	if (sscanf(line, "object %*d class gridpositions counts %d %d %d", 
		   &nx, &ny, &nz) == 3) {
	    if ((nx < 0) || (ny < 0) || (nz < 0)) {
		result.addError("invalid grid size: x=%d, y=%d, z=%d", 
			nx, ny, nz);
		return false;
	    }
	} else if (sscanf(line, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
	    // found origin
	} else if (sscanf(line, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
	    // found one of the delta lines
	    if (ddx != 0.0) {
		dx = ddx;
	    } else if (ddy != 0.0) {
		dy = ddy;
	    } else if (ddz != 0.0) {
		dz = ddz;
	    }
	} else if (sscanf(line, "object %*d class array type %*s shape 3"
		" rank 1 items %d data follows", &npts) == 1) {
	    printf("#points=%d\n", npts);
	    if (npts != nx*ny*nz) {
		result.addError("inconsistent data: expected %d points"
				" but found %d points", nx*ny*nz, npts);
		return false;
	    }
	    break;
	} else if (sscanf(line, "object %*d class array type %*s rank 0"
		" times %d data follows", &npts) == 1) {
	    if (npts != nx*ny*nz) {
		result.addError("inconsistent data: expected %d points"
				" but found %d points", nx*ny*nz, npts);
		return false;
	    }
	    break;
	}
    }
    if (npts != nx*ny*nz) {
	result.addError("inconsistent data: expected %d points"
			" but found %d points", nx*ny*nz, npts);
	return false;
    }

    _initialized = false;
    _xValueMin = _yValueMin = _zValueMin = FLT_MAX;
    _xValueMax = _yValueMax = _zValueMax = -FLT_MAX;
    _xMin = x0, _yMin = y0, _zMin = z0;
    _xNum = nx, _yNum = ny, _zNum = nz;
    _xMax = _xMin + dx * _xNum;
    _yMax = _yMin + dy * _yNum;
    _zMax = _zMin + dz * _zNum;
    _nComponents = nComponents;

    if (_values != NULL) {
	delete [] _values;
    }
    _values = new float[npts * _nComponents];
    _nValues = 0;
    for (size_t ix = 0; ix < _xNum; ix++) {
	for (size_t iy = 0; iy < _yNum; iy++) {
	    for (size_t iz = 0; iz < _zNum; iz++) {
		char *line;
		if ((p == endPtr) || (_nValues > (size_t)npts)) {
		    break;
		}
		line = getline(&p, endPtr);
		if ((line[0] == '#') || (line[0] == '\0')) {
		    continue;	// Skip blank or comment lines.
		}
		double vx, vy, vz;
		if (sscanf(line, "%lg %lg %lg", &vx, &vy, &vz) == 3) {
		    int nindex = (iz*nx*ny + iy*nx + ix) * 3;
		    if (vx < _xValueMin) {
			_xValueMin = vx;
		    } else if (vx > _xValueMax) {
			_xValueMax = vx;
		    }
		    if (vy < _yValueMin) {
			_yValueMin = vy;
		    } else if (vy > _yValueMax) {
			_yValueMax = vy;
		    }
		    if (vz < _zValueMin) {
			_zValueMin = vz;
		    } else if (vz > _zValueMax) {
			_zValueMax = vz;
		    }
		    _values[nindex] = vx;
		    _values[nindex+1] = vy;
		    _values[nindex+2] = vz;
		    _nValues++;
		}
	    }
	}
    }
    // make sure that we read all of the expected points
    if (_nValues != npts) {
	result.addError("inconsistent data: expected %d points"
			" but found %d points", npts, _nValues);
	delete []  _values;
	_values = NULL;
	return false;
    }
    _nValues *= _nComponents;
    _initialized = true;
    return true;
}


bool
Rappture::Unirect3d::Resample(Rappture::Outcome &result, int nSamples)
{
    Rappture::Mesh1D xgrid(_xMin, _xMax, _xNum);
    Rappture::Mesh1D ygrid(_yMin, _yMax, _yNum);
    Rappture::Mesh1D zgrid(_zMin, _zMax, _zNum);
    Rappture::FieldRect3D xfield(xgrid, ygrid, zgrid);
    Rappture::FieldRect3D yfield(xgrid, ygrid, zgrid);
    Rappture::FieldRect3D zfield(xgrid, ygrid, zgrid);

    size_t i, j;
    for (i = 0, j = 0; i < _nValues; i += _nComponents, j++) {
	double vx, vy, vz;

	vx = _values[i];
	vy = _values[i+1];
	vz = _values[i+2];
	
	xfield.define(j, vx);
	yfield.define(j, vy);
	zfield.define(j, vz);
    }
    // Figure out a good mesh spacing
    double dx, dy, dz;
    dx = xfield.rangeMax(Rappture::xaxis) - xfield.rangeMin(Rappture::xaxis);
    dy = xfield.rangeMax(Rappture::yaxis) - xfield.rangeMin(Rappture::yaxis);
    dz = xfield.rangeMax(Rappture::zaxis) - xfield.rangeMin(Rappture::zaxis);

    double dmin;
    dmin = pow((dx*dy*dz)/(nSamples*nSamples*nSamples), 0.333);
    
    printf("dx:%lf dy:%lf dz:%lf dmin:%lf\n", dx, dy, dz, dmin);

    /* Recompute new number of points for each axis. */
    _xNum = (size_t)ceil(dx/dmin);
    _yNum = (size_t)ceil(dy/dmin);
    _zNum = (size_t)ceil(dz/dmin);
    
#ifndef NV40
    // must be an even power of 2 for older cards
    _xNum = (int)pow(2.0, ceil(log10((double)_xNum)/log10(2.0)));
    _yNum = (int)pow(2.0, ceil(log10((double)_yNum)/log10(2.0)));
    _zNum = (int)pow(2.0, ceil(log10((double)_zNum)/log10(2.0)));
#endif

    size_t n = _nComponents * _xNum * _yNum * _zNum;
    float *data = new float[n];
    memset(data, 0, sizeof(float) * n);
    
    // Generate the uniformly sampled rectangle that we need for a volume
    float *destPtr = data;
    for (size_t i = 0; i < _zNum; i++) {
	double z;

	z = _zMin + (i * dmin);
	for (size_t j = 0; j < _yNum; j++) {
	    double y;
		
	    y = _yMin + (j * dmin);
	    for (size_t k = 0; k < _xNum; k++) {
		double x;

		x = _xMin + (k * dmin);
		destPtr[0] = xfield.value(x, y, z);
		destPtr[1] = yfield.value(x, y, z);
		destPtr[2] = zfield.value(x, y, z);
	    }
	}
    }
    delete [] _values;
    _values = data;
    _nValues = _xNum * _yNum * _zNum * _nComponents;
    return true;
}


void
Rappture::Unirect3d::GetVectorRange(void)
{
    assert(_nComponents == 3);
    _magMax = -DBL_MAX, _magMin = DBL_MAX;
    size_t i;
    for (i = 0; i < _nValues; i += _nComponents) {
	double vx, vy, vz, vm;

	vx = _values[i];
	vy = _values[i+1];
	vz = _values[i+2];
		    
	vm = sqrt(vx*vx + vy*vy + vz*vz);
	if (vm > _magMax) {
	    _magMax = vm;
	}
	if (vm < _magMin) {
	    _magMin = vm;
	}
    }
}