/***********************************************************************/
/* Open Visualization Data Explorer                                    */
/* (C) Copyright IBM Corp. 1989,1999                                   */
/* ALL RIGHTS RESERVED                                                 */
/* This code licensed under the                                        */
/*    "IBM PUBLIC LICENSE - Open Visualization Data Explorer"          */
/***********************************************************************/

#include <dxconfig.h>



#include <stdio.h>
#include <math.h>
#include <string.h>
#include <dx/dx.h>
#include "quadsRR2DClass.h"

static void  _dxfCleanup(QuadsRR2DInterpolator);
static Error _dxfInitializeTask(Pointer);
static Error _dxfInitialize(QuadsRR2DInterpolator qi);

#define DIAGNOSTIC(str) \
	DXMessage("quads interpolator: %s\n", (str));

int
_dxfRecognizeQuadsRR2D(Field field)
{
    Array 	array;
    Array 	topology;
    int	  	nDim;
    int	  	status;
    Type	type;
    Category	cat;
    int		rank, shape[32];

    CHECK(field, CLASS_FIELD);

    ELT_TYPECHECK(field, "quads");

    status = OK;

    array = (Array)DXGetComponentValue(field, "positions");
    if (!array)
	return 0;
    
    if (!DXQueryGridPositions(array, &nDim, NULL, NULL, NULL))
	return 0;

    if (nDim != 2)
	return 0;

    if (! DXGetArrayInfo(array, NULL, &type, &cat, &rank, shape))
	return 0;

    if (cat != CATEGORY_REAL)
	return 0;

    if (rank != 1 || shape[0] != 2)
	return 0;

    array = (Array)DXGetComponentValue(field, "connections");
    if (!array)
	return 0;
    
    if (!DXQueryGridConnections(array, &nDim, NULL))
	return 0;
    
    if (nDim != 2)
	return 0;

    return status;
}

QuadsRR2DInterpolator 
_dxfNewQuadsRR2DInterpolator(Field field, 
		enum interp_init initType, double fuzz, Matrix *m)
{
    return (QuadsRR2DInterpolator)_dxf_NewQuadsRR2DInterpolator(field, 
			initType, fuzz, m, &_dxdquadsrr2dinterpolator_class);
}


QuadsRR2DInterpolator 
_dxf_NewQuadsRR2DInterpolator(Field field, 
		enum interp_init initType, float fuzz, Matrix *m,
		struct quadsrr2dinterpolator_class *class)
{
    QuadsRR2DInterpolator qi;

    qi = (QuadsRR2DInterpolator)_dxf_NewFieldInterpolator(field, fuzz, m,
			(struct fieldinterpolator_class *)class);
    if (! qi)
	return NULL;

    qi->pointsArray = NULL;
    qi->dataArray   = NULL;

    /*
     * The grid should be regular
     */
    {
	float localOrigin[2], displacement[2];
	Array connections, positions;
	float deltas[4], d;

	positions = (Array)DXGetComponentValue(field, "positions");
	connections = (Array)DXGetComponentValue(field, "connections");

	if (!positions || !connections)
	{
	    DXDelete((Object)(qi));
	    return NULL;
	}

	DXQueryGridPositions(positions, NULL, NULL, localOrigin, deltas);

	d = 1.0 / (deltas[0]*deltas[3] - deltas[1]*deltas[2]);

	((float *)((Interpolator)qi)->matrix.A)[0] =  deltas[3] * d;
	((float *)((Interpolator)qi)->matrix.A)[1] = -deltas[1] * d;
	((float *)((Interpolator)qi)->matrix.A)[2] = -deltas[2] * d;
	((float *)((Interpolator)qi)->matrix.A)[3] =  deltas[0] * d;

	DXGetMeshOffsets((MeshArray)connections, qi->meshOffsets);

	((Interpolator)qi)->matrix.b[0] = localOrigin[0]
			- (qi->meshOffsets[0]*deltas[0] + qi->meshOffsets[1]*deltas[1]);
	((Interpolator)qi)->matrix.b[1] = localOrigin[1]
			- (qi->meshOffsets[0]*deltas[2] + qi->meshOffsets[1]*deltas[3]);
    }

    if (initType == INTERP_INIT_PARALLEL)
    {
	if (! DXAddTask(_dxfInitializeTask, (Pointer)&qi, sizeof(qi), 1.0))
	{
	    DXDelete((Object)qi);
	    return NULL;
	}
    }
    else if (initType == INTERP_INIT_IMMEDIATE)
    {
	if (! _dxfInitialize(qi))
	{
	    DXDelete((Object)qi);
	    return NULL;
	}
    }

    return qi;
}

static Error
_dxfInitializeTask(Pointer p)
{
    return _dxfInitialize(*(QuadsRR2DInterpolator *)p);
}


static Error
_dxfInitialize(QuadsRR2DInterpolator qi)
{
    Field		field;
    int			nDim;
    int			i;
    float		deltas[4], d;
    Type		dataType;
    Category		dataCategory;
    float		edge0, edge1;
    Point		diagonal;

    field = (Field)((Interpolator)qi)->dataObject;

    qi->fieldInterpolator.initialized = 1;

    /*
     * De-reference data
     */
    qi->dataArray   = (Array)DXGetComponentValue(field, "data");
    if (!qi->dataArray) 
    {
	DXSetError(ERROR_MISSING_DATA, "#10240", "data");
	return ERROR;
    }
    DXReference((Object)qi->dataArray);

    DXGetArrayInfo(qi->dataArray, NULL, &((Interpolator)qi)->type,
	&((Interpolator)qi)->category, &((Interpolator)qi)->rank, 
	((Interpolator)qi)->shape);

    /*
     * Get the grid. 
     */
    qi->pointsArray = (Array)DXGetComponentValue(field, "positions");
    if (!qi->pointsArray)
    {
	DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
	return ERROR;
    }
    DXReference((Object)qi->pointsArray);

    /*
     * get info about data values
     */
    DXGetArrayInfo(qi->dataArray, NULL, &dataType, &dataCategory, NULL, NULL);

    /*
     * Don't worry about maintaining shape of input; just determine how 
     * many values to interpolate.
     */
    qi->nElements = DXGetItemSize(qi->dataArray) / 
				DXTypeSize(((Interpolator)qi)->type);

    /*
     * The grid should be regular
     */
    DXQueryGridPositions(qi->pointsArray, &nDim, qi->counts, NULL, deltas);

    /*
     * Figure out how big to increment pointers along each dimension
     * Use either number of positions or number of primitive elements along
     * each axis depending on data dependency.
     */
    if (qi->fieldInterpolator.data_dependency == DATA_POSITIONS_DEPENDENT)
    {
	qi->size[1] = 1;
	for (i = 0; i >= 0; i--)
	    qi->size[i] = qi->counts[i+1] * qi->size[i+1];
    }
    else
    {
	qi->size[1] = 1;
	for (i = 0; i >= 0; i--)
	    qi->size[i] = (qi->counts[i+1] - 1) * qi->size[i+1];
    }

    qi->dHandle = DXCreateArrayHandle(qi->dataArray);
    if (! qi->dHandle)
	return ERROR;

    return OK;
}

Error
_dxfQuadsRR2DInterpolator_Delete(QuadsRR2DInterpolator qi)
{
    _dxfCleanup(qi);
    _dxfFieldInterpolator_Delete((FieldInterpolator) qi);
}

int
_dxfQuadsRR2DInterpolator_PrimitiveInterpolate(QuadsRR2DInterpolator qi,
			int *n, float **points, Pointer *values, int fuzzFlag)
{
    float	    x, y;
    int 	    ix, iy;
    float	    dx, dy;
    float 	    xMax, xMin;
    float 	    yMax, yMin;
    int		    ixMax;
    int		    iyMax;
    int		    sz0, sz1;
    float	    *org;
    float	    *iD;
    float	    value00, value01, value10, value11;
    float    	    weight00, weight01, weight10, weight11;
    float	    dxdy;
    int		    baseIndex, i, knt;
    float	    *p;
    Pointer	    v;
    float	    *A, *b;
    float	    fuzz;
    int		    dep;
    int		    itemSize;
    int		    typeSize;
    int 	    *mo;
    float	    ptx, pty;
    InvalidComponentHandle icH = ((FieldInterpolator)qi)->invCon;
    char	    *dbuf = NULL;
    Matrix	    *xform;

    if (! qi->fieldInterpolator.initialized)
    {
	if (! _dxfInitialize(qi))
	{
	    _dxfCleanup(qi);
	    return 0;
	}

	qi->fieldInterpolator.initialized = 1;
    }

    /*
     * De-reference indexing info from interpolator
     */
    A  = (float *)((Interpolator)qi)->matrix.A;
    b  = (float *)((Interpolator)qi)->matrix.b;
    mo = qi->meshOffsets;

    sz0 = qi->size[0];
    sz1 = qi->size[1];

    dep = qi->fieldInterpolator.data_dependency;
    typeSize = DXTypeSize(((Interpolator)qi)->type);
    itemSize = qi->nElements * typeSize;

    fuzz = qi->fuzz;

    ixMax = qi->counts[0] - 1;
    iyMax = qi->counts[1] - 1;

    /*
     * De-reference bounding box
     */
    xMax = ((Interpolator)qi)->max[0] + fuzz;
    xMin = ((Interpolator)qi)->min[0] - fuzz;
    yMax = ((Interpolator)qi)->max[1] + fuzz;
    yMin = ((Interpolator)qi)->min[1] - fuzz;

    if (((FieldInterpolator)qi)->xflag)
	xform = &(((FieldInterpolator)qi)->xform);
    else
	xform = NULL;


    p = *points;
    v = *values;

    dbuf = (char *)DXAllocate(4*itemSize);
    if (! dbuf)
	return ERROR;

    while (*n > 0)
    {
	float xpt[2];
	float *pPtr;

	if (xform)
	{
	    xpt[0] = p[0]*xform->A[0][0] +
		     p[1]*xform->A[1][0] +
		          xform->b[0];

	    xpt[1] = p[0]*xform->A[0][1] +
		     p[1]*xform->A[1][1] +
		          xform->b[1];

	    pPtr = xpt;
	}
	else
	    pPtr = p;


	ptx = pPtr[0] - b[0];
	pty = pPtr[1] - b[1];

	x = ptx*A[0] + pty*A[2] - mo[0];
	y = ptx*A[1] + pty*A[3] - mo[1];

	if (fuzzFlag == FUZZ_ON)
	{
	    if ((x < fuzz        || 
		 x > ixMax+fuzz) ||
		(y < -fuzz       || 
		 y > iyMax+fuzz)) break;
	}
	else
	{
	    if ((x < 0      || 
		 x > ixMax) ||
		(y < 0      || 
		 y > iyMax)) break;
	}

	/*
	 * Note: assuming we passed the above test, we really will
	 * be interpolating.
	 */

#define INTERPOLATE(type, round)				\
{								\
    type *v00, *v10, *v01, *v11, *r;				\
								\
    baseIndex = ix*sz0 + iy*sz1;				\
								\
    v00 = (type *)DXGetArrayEntry(qi->dHandle, 			\
			baseIndex, (Pointer)dbuf);		\
								\
    v10 = (type *)DXGetArrayEntry(qi->dHandle, 			\
		baseIndex+sz0, (Pointer)(dbuf+itemSize));	\
								\
    v01 = (type *)DXGetArrayEntry(qi->dHandle, 			\
		baseIndex+sz1, (Pointer)(dbuf+2*itemSize));	\
								\
    v11 = (type *)DXGetArrayEntry(qi->dHandle, 			\
		baseIndex+sz1+sz0, (Pointer)(dbuf+3*itemSize));	\
								\
    r = (type *)v;						\
								\
    for (i = 0; i < qi->nElements; i++)				\
	*r++ = weight00*(*v00++) + weight01*(*v01++) +		\
	       weight10*(*v10++) + weight11*(*v11++) +		\
	       round;						\
    								\
    v = (Pointer)r;						\
}

	if (((FieldInterpolator)qi)->cstData)
	{
	    memcpy(v, ((FieldInterpolator)qi)->cstData, itemSize);
	    v = (Pointer)(((char *)v) + itemSize);
	}
	else if (dep == DATA_POSITIONS_DEPENDENT)
	{
            Type dataType;

	    ix = x;
	    dx = x - ix;
	    if (fuzzFlag != FUZZ_ON && (ix>ixMax || (ix==ixMax && dx>0.0) || ix<0))
		goto not_found;

	    if (ix >= ixMax)
	    {
		ix = ixMax - 1;
		dx = 1.0;
	    }
	    else if (ix < 0 || dx < 0.0)
	    {
		ix = 0;
		dx = 0.0;
	    }

	    iy = y;
	    dy = y - iy;
	    if (fuzzFlag != FUZZ_ON && (iy>iyMax || (iy==iyMax && dy>0.0) || iy<0))
		goto not_found;

	    if (iy >= iyMax)
	    {
		iy = iyMax - 1;
		dy = 1.0;
	    }
	    else if (iy < 0 || dy < 0.0)
	    {
		iy = 0;
		dy = 0.0;
	    }

	    if (icH)
	    {
		baseIndex = ix*(qi->counts[1] - 1) + iy;
		if (DXIsElementInvalid(icH, baseIndex))
		    goto not_found;
	    }

	    dxdy = dx * dy;

	    weight11 = dxdy;
	    weight10 = dx - dxdy;
	    weight01 = dy - dxdy;
	    weight00 = 1 - dx - dy + dxdy;

            if ((dataType = ((Interpolator)qi)->type) == TYPE_FLOAT)
            {
                INTERPOLATE(float, 0.0);
            }
            else if (dataType == TYPE_DOUBLE)
            {
                INTERPOLATE(double, 0.0);
            }
            else if (dataType == TYPE_INT)
            {
                INTERPOLATE(int, 0.5);
            }
            else if (dataType == TYPE_SHORT)
            {
                INTERPOLATE(short, 0.5);
            }
            else if (dataType == TYPE_USHORT)
            {
                INTERPOLATE(ushort, 0.5);
            }
            else if (dataType == TYPE_UINT)
            {
                INTERPOLATE(uint, 0.5);
            }
            else if (dataType == TYPE_BYTE)
            {
                INTERPOLATE(byte, 0.5);
            }
            else if (dataType == TYPE_UBYTE)
            {
                INTERPOLATE(ubyte, 0.5);
            }
            else
            {
                INTERPOLATE(unsigned char, 0.5);
            }
	}
	else
	{
	    ix = x;
	    if (fuzzFlag != FUZZ_ON &&
		(ix>ixMax || (ix==ixMax && dx>0.0) || ix<0))
		goto not_found;

	    if (ix >= ixMax)
		ix = ixMax - 1;
	    else if (ix < 0)
		ix = 0;

	    iy = y;
	    if (fuzzFlag != FUZZ_ON &&
		(iy>iyMax || (iy==iyMax && dy>0.0) || iy<0))
		goto not_found;

	    if (iy >= iyMax)
		iy = iyMax - 1;
	    else if (iy < 0)
		iy = 0;

	    if (icH)
	    {
		baseIndex = ix*(qi->counts[1] - 1) + iy;
		if (DXIsElementInvalid(icH, baseIndex))
		    goto not_found;
	    }

	    baseIndex = ix*sz0 + iy*sz1;

	    memcpy(v, DXGetArrayEntry(qi->dHandle,
		baseIndex, (Pointer)dbuf), itemSize);

	    v = (Pointer)(((char *)v) + itemSize);
	}

	/*
	 * Only use fuzz on first sample
	 */
	fuzzFlag = FUZZ_OFF;
	fuzz = 0.0;

	/*
	 * Only use fuzz on first sample
	 */
	p += 2;
	(*n)--;
    }

not_found:

    *points = (float *)p;
    *values = (Pointer)v;

    if (dbuf)
	DXFree((Pointer)dbuf);

    return OK;
}

static void
_dxfCleanup(QuadsRR2DInterpolator qi)
{								
    if (qi->dHandle)
    {
	DXFreeArrayHandle(qi->dHandle);
	qi->dHandle = NULL;
    }

    if (qi->pointsArray)
    {
	DXDelete((Object)qi->pointsArray);
	qi->pointsArray = NULL;
    }

    if (qi->dataArray)
    {
	DXDelete((Object)qi->dataArray);
	qi->dataArray = NULL;
    }
}

Object
_dxfQuadsRR2DInterpolator_Copy(QuadsRR2DInterpolator old, enum copy copy)
{
    QuadsRR2DInterpolator new;

    new = (QuadsRR2DInterpolator)
	    _dxf_NewObject((struct object_class *)&_dxdquadsrr2dinterpolator_class);

    if (!(_dxf_CopyQuadsRR2DInterpolator(new, old, copy)))
    {
	DXDelete((Object)new);
	return NULL;
    }
    else
	return (Object)new;
}

QuadsRR2DInterpolator
_dxf_CopyQuadsRR2DInterpolator(QuadsRR2DInterpolator new, 
				QuadsRR2DInterpolator old, enum copy copy)
{

    if (! _dxf_CopyFieldInterpolator((FieldInterpolator)new,
					(FieldInterpolator)old, copy))
	return NULL;
    
    memcpy((char *)new->size,   (char *)old->size,   2*sizeof(int));
    memcpy((char *)new->counts, (char *)old->counts, 2*sizeof(float));
    memcpy((char *)new->meshOffsets, (char *)old->meshOffsets, sizeof(old->meshOffsets));

    new->fuzz       = old->fuzz;
    new->nElements  = old->nElements;

    if (new->fieldInterpolator.initialized)
    {
	new->pointsArray    = (Array)DXReference((Object)old->pointsArray);
	new->dataArray      = (Array)DXReference((Object)old->dataArray);

	new->dHandle = DXCreateArrayHandle(new->dataArray);
	if (! new->dHandle)
	    return ERROR;
    }
    else
    {
	new->pointsArray    = NULL;
	new->dataArray      = NULL;
	new->dHandle	    = NULL;
    }

    if (DXGetError())
	return NULL;

    return new;
}

Interpolator
_dxfQuadsRR2DInterpolator_LocalizeInterpolator(QuadsRR2DInterpolator qi)
{
    if (qi->fieldInterpolator.localized)
	return (Interpolator)qi;

    qi->fieldInterpolator.localized = 1;

    if (qi->fieldInterpolator.initialized)
	qi->dHandle = DXCreateArrayHandle(qi->dataArray);

    if (DXGetError())
        return NULL;
    else
        return (Interpolator)qi;
}