/***********************************************************************/
/* 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 "fle2DClass.h"

#define FALSE 0
#define TRUE 1

static Error _dxfInitializeTask(Pointer);
static Error _dxfInitialize(Fle2DInterpolator);
static int   _dxfCleanup(Fle2DInterpolator);

static int CountIntersections(Fle2DInterpolator, float *, int);
static int InsideFace(Fle2DInterpolator, float *, int);
static int FindFace(Fle2DInterpolator, float *);

int
_dxfRecognizeFLE2D(Field field)
{
    Array    array;
    int      nDim;
    Type     t;
    Category c;
    int      r, s[32];
    Object   depAttr;

    CHECK(field, CLASS_FIELD);

    array = (Array)DXGetComponentValue(field, "faces");
    if (! array)
	return 0;

    DXGetArrayInfo(array, NULL, &t, &c, &r, s);

    if (t != TYPE_INT || c != CATEGORY_REAL ||
			!(r == 0 || (r == 1 && s[0] == 1)))
    {
	DXSetError(ERROR_INVALID_DATA, "invalid faces component");
	return 0;
    }

    array = (Array)DXGetComponentValue(field, "loops");
    if (! array)
	return 0;

    DXGetArrayInfo(array, NULL, &t, &c, &r, s);

    if (t != TYPE_INT || c != CATEGORY_REAL ||
			!(r == 0 || (r == 1 && s[0] == 1)))
    {
	DXSetError(ERROR_INVALID_DATA, "invalid loops component");
	return 0;
    }

    array = (Array)DXGetComponentValue(field, "edges");
    if (! array)
	return 0;

    DXGetArrayInfo(array, NULL, &t, &c, &r, s);

    if (t != TYPE_INT || c != CATEGORY_REAL ||
			!(r == 0 || (r == 1 && s[0] == 1)))
    {
	DXSetError(ERROR_INVALID_DATA, "invalid edges component");
	return 0;
    }

    array = (Array)DXGetComponentValue(field, "positions");
    if (!array)
    {
	DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
	return 0;
    }

    DXGetArrayInfo(array, NULL, &t, &c, &r, s);
	
    if (t != TYPE_FLOAT || c != CATEGORY_REAL || r != 1 || s[0] != 2)
    {
	DXSetError(ERROR_INVALID_DATA, 
		"fle positions must be float and 2-D");
	return 0;
    }

    depAttr = DXGetComponentAttribute(field, "data", "dep");
    if (! depAttr)
    {
	DXSetError(ERROR_INVALID_DATA, "missing data dependency");
	return 0;
    }

    if (DXGetObjectClass(depAttr) != CLASS_STRING)
    {
	DXSetError(ERROR_INVALID_DATA, "invalid data dependency attribute");
	return 0;
    }

    if (strcmp(DXGetString((String)depAttr), "faces"))
    {
	DXSetError(ERROR_INVALID_DATA, 
		"invalid data dependence: fle data must be dep faces");
	return 0;
    }
    
    return 1;
}

Fle2DInterpolator
_dxfNewFle2DInterpolator(Field field,
		enum interp_init initType, double fuzz, Matrix *m) 
{
    return (Fle2DInterpolator)_dxf_NewFle2DInterpolator(field, 
			    initType, fuzz, m, &_dxdfle2dinterpolator_class);
}

Fle2DInterpolator
_dxf_NewFle2DInterpolator(Field field,
		enum interp_init initType, float fuzz, Matrix *m,
		struct fle2dinterpolator_class *class)
{
    Fle2DInterpolator 	fle;
    float	        *mm, *MM;

    fle = (Fle2DInterpolator)_dxf_NewFieldInterpolator(field, fuzz, m,
	    (struct fieldinterpolator_class *)class);

    if (! fle)
	return NULL;

    fle->fArray = NULL;
    fle->lArray = NULL;
    fle->eArray = NULL;
    fle->pArray = NULL;
    fle->dArray = NULL;
    fle->mmfArray = NULL;
    fle->mmlArray = NULL;

    mm = ((Interpolator)fle)->min;
    MM = ((Interpolator)fle)->max;

    if (((MM[0] - mm[0]) * (MM[1] - mm[1])) == 0.0)
    {
	DXDelete((Object)fle);
	return NULL;
    }

    fle->hint = -1;

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

    }
    else if (initType == INTERP_INIT_IMMEDIATE)
    {
	if (! _dxfInitialize(fle))
	{
	    DXDelete((Object)fle);
	    return NULL;
	}
    }

    return fle;
}

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

static Error
_dxfInitialize(Fle2DInterpolator fle)
{
    Field	field;
    int		nDim;
    Type	dataType;
    Category	dataCategory;
    int		i, j, k;
    int 	nPoints;
    float	lx, ly, lX, lY;
    float	fx, fy, fX, fY;
    float	fuzz;
    float	len, area;
    int		*faces, *loops, *edges, *faceEdges = NULL;
    float	*positions;
    MinMax	*lmm, *fmm;

    fle->fieldInterpolator.initialized = 1;

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

    fle->fArray = (Array)DXGetComponentValue(field, "faces");
    fle->lArray = (Array)DXGetComponentValue(field, "loops");
    fle->eArray = (Array)DXGetComponentValue(field, "edges");
    fle->pArray = (Array)DXGetComponentValue(field, "positions");
    fle->dArray = (Array)DXGetComponentValue(field, "data");

    DXReference((Object)fle->fArray);
    DXReference((Object)fle->lArray);
    DXReference((Object)fle->eArray);
    DXReference((Object)fle->pArray);
    DXReference((Object)fle->dArray);

    fle->faces     = (int   *)DXGetArrayData(fle->fArray);
    fle->loops     = (int   *)DXGetArrayData(fle->lArray);
    fle->edges     = (int   *)DXGetArrayData(fle->eArray);
    fle->positions = (float *)DXGetArrayData(fle->pArray);

    DXGetArrayInfo(fle->fArray, &fle->nFaces, NULL, NULL, NULL, NULL);
    DXGetArrayInfo(fle->lArray, &fle->nLoops, NULL, NULL, NULL, NULL);
    DXGetArrayInfo(fle->eArray, &fle->nEdges, NULL, NULL, NULL, NULL);

    fle->mmfArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 4);
    fle->mmlArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 4);
    if (! fle->mmfArray || ! fle->mmlArray)
	goto error;

    DXReference((Object)fle->mmfArray);
    DXReference((Object)fle->mmlArray);

    if (! DXAddArrayData(fle->mmfArray, 0, fle->nFaces, NULL))
	goto error;

    if (! DXAddArrayData(fle->mmlArray, 0, fle->nLoops, NULL))
	goto error;

    fmm = fle->minmax_faces = (MinMax *)DXGetArrayData(fle->mmfArray);
    if (! fmm)
	goto error;

    lmm = fle->minmax_loops = (MinMax *)DXGetArrayData(fle->mmlArray);
    if (! lmm)
	goto error;

    for (i = 0; i < fle->nFaces; i++, fmm++)
    {
	int sl = fle->faces[i];
	int el = (i == fle->nFaces-1) ? fle->nLoops : fle->faces[i+1];

	fx = fy =  DXD_MAX_FLOAT;
	fX = fY = -DXD_MAX_FLOAT;

	for (j = sl; j < el; j++, lmm++)
	{
	    int se = fle->loops[j];
	    int ee = (j == fle->nLoops-1) ? fle->nEdges : fle->loops[j+1];

	    lx = ly =  DXD_MAX_FLOAT;
	    lX = lY = -DXD_MAX_FLOAT;
	
	    for (k = se; k < ee; k++)
	    {
		float *xy = fle->positions+(2*fle->edges[k]);

		if (*xy < lx) lx = *xy;
		if (*xy > lX) lX = *xy;
		xy++;
		if (*xy < ly) ly = *xy;
		if (*xy > lY) lY = *xy;
	    }

            if (lx < fx) fx = lx;
	    if (lX > fX) fX = lX;
	    if (ly < fy) fy = ly;
	    if (lY > fY) fY = lY;

	    lmm->x = lx;
	    lmm->y = ly;
	    lmm->X = lX;
	    lmm->Y = lY;
	}

	fmm->x = fx;
	fmm->y = fy;
	fmm->X = fX;
	fmm->Y = fY;
    }

    DXGetArrayInfo(fle->dArray, NULL, &((Interpolator)fle)->type,
		&((Interpolator)fle)->category,
		&((Interpolator)fle)->rank, ((Interpolator)fle)->shape);

    fle->dHandle = DXCreateArrayHandle(fle->dArray);
    if (! fle->dHandle)
	return ERROR;


    return OK;

error:

    DXDelete((Pointer)fle->mmfArray);
    fle->mmfArray = NULL;

    DXDelete((Pointer)fle->mmlArray);
    fle->mmlArray = NULL;

    return ERROR;
}

Error
_dxfFle2DInterpolator_Delete(Fle2DInterpolator fle)
{
    _dxfCleanup(fle);
    _dxfFieldInterpolator_Delete((FieldInterpolator) fle);
}

static int
CountIntersections(Fle2DInterpolator fle, float *xy, int loop)
{
    int   se = fle->loops[loop];
    int   le = (loop == fle->nLoops-1) ? fle->nEdges : fle->loops[loop+1];
    int   e, p, q, knt;
    float t, x, *ppos, *qpos, *top, *bottom;

    p    = fle->edges[se];
    ppos = fle->positions + (p<<1);
    for (e = se, knt = 0; e < le; e++, p = q, ppos = qpos)
    {
	if (e == le-1)
	    q = fle->edges[se];
	else
	    q = fle->edges[e+1];
    
	qpos = fle->positions + (q<<1);

	if (qpos[1] > ppos[1])
	{
	    top = qpos; bottom = ppos;
	}
	else if (qpos[1] < ppos[1])
	{
	    top = ppos; bottom = qpos;
	}
	else
	    continue;
	
	if (top[1] < xy[1])
	    continue;
	
	if (bottom[1] >= xy[1])
	    continue;

	if (top[0] < xy[0] && bottom[0] < xy[0])
	    continue;
	
	if (top[0] > xy[0] && bottom[0] > xy[0])
	{
	    knt ++;
	    continue;
	}
	
	t = (xy[1] - top[1]) / (bottom[1] - top[1]);
	x = top[0] + t*(bottom[0] - top[0]);

	if (x >= xy[0])
	{
	    knt ++;
	    continue;
	}
    }

    return knt;
}

static int 
InsideFace(Fle2DInterpolator fle, float *xy, int face)
{
    int    sl = fle->faces[face];
    int    el = (face == fle->nFaces-1) ? fle->nLoops : fle->faces[face+1];
    int    knt;
    int    l;
    MinMax *mm = fle->minmax_loops + sl;

    for (l = sl, knt = 0; l < el; l++, mm++)
    {
	if (mm->x < xy[0] && mm->X > xy[0] &&
	    mm->y < xy[1] && mm->Y > xy[1])
	    knt += CountIntersections(fle, xy, l);
    }
	
    if (knt & 0x1)
	return face;

    return -1;
}


static int 
FindFace(Fle2DInterpolator fle, float *xy)
{
    int i;

    if (fle->hint >= 0)
    {
	MinMax *mm = fle->minmax_faces + fle->hint;

	if (mm->x < xy[0] && mm->X > xy[0] &&
			    mm->y < xy[1] && mm->Y > xy[1])
	{
	    if (-1 != InsideFace(fle, xy, fle->hint))
		return fle->hint;
	}
    }

    for (i = 0; i < fle->nFaces; i++)
    {
	MinMax *mm = fle->minmax_faces + i;

	if (mm->x < xy[0] && mm->X > xy[0] &&
			    mm->y < xy[1] && mm->Y > xy[1])
	{
	    if (-1 != InsideFace(fle, xy, i))
	    {
		fle->hint = i;
		return i;
	    }
	}
    }

    fle->hint = -1;
    return -1;
}

int
_dxfFle2DInterpolator_PrimitiveInterpolate(Fle2DInterpolator fle,
		    int *n, float **points, Pointer *values, int fuzzFlag)
{
    int found;
    Pointer v;
    float *p;
    int itemSize;
    Matrix *xform;
    char *dbuf = NULL;

    if (! fle->fieldInterpolator.initialized)
    {
	if (! _dxfInitialize(fle))
	{
	    _dxfCleanup(fle);
	    return ERROR;
	}

	fle->fieldInterpolator.initialized = 1;
    }

    itemSize = DXGetItemSize(fle->dArray);

    dbuf = (char *)DXAllocate(itemSize);
    if (! dbuf)
	return ERROR;
    
    if (((FieldInterpolator)fle)->xflag)
	xform = &(((FieldInterpolator)fle)->xform);
    else
	xform = NULL;

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

    /*
     * For each point in the input, attempt to interpolate the point.
     * When a point cannot be interpolated, quit.
     */
    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;


	found = FindFace(fle, pPtr);

	if (found == -1)
	    break;
	
	memcpy(v, DXGetArrayEntry(fle->dHandle, found, dbuf), itemSize);

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

	/*
	 * Only use fuzz on first point
	 */
	fuzzFlag = FUZZ_OFF;

	p += 2;
	*n -= 1;
    }

    DXFree((Pointer)dbuf);

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

    return OK;

error:
    DXFree((Pointer)dbuf);
    return ERROR;
}

static int
_dxfCleanup(Fle2DInterpolator fle)
{								
    if (fle->dHandle)
    {
	DXFreeArrayHandle(fle->dHandle);
	fle->dHandle = NULL;
    }

    DXDelete((Object)fle->fArray);
    DXDelete((Object)fle->lArray);
    DXDelete((Object)fle->eArray);
    DXDelete((Object)fle->pArray);
    DXDelete((Object)fle->dArray);
    DXDelete((Object)fle->mmfArray);
    DXDelete((Object)fle->mmlArray);

    fle->fArray   = NULL;
    fle->lArray   = NULL;
    fle->eArray   = NULL;
    fle->pArray   = NULL;
    fle->dArray   = NULL;
    fle->mmfArray = NULL;
    fle->mmlArray = NULL;
}

Object
_dxfFle2DInterpolator_Copy(Fle2DInterpolator old, enum copy copy)
{
    Fle2DInterpolator new;

    new = (Fle2DInterpolator)
	    _dxf_NewObject((struct object_class *)&_dxdfle2dinterpolator_class);

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

Fle2DInterpolator
_dxf_CopyFle2DInterpolator(Fle2DInterpolator new, 
				Fle2DInterpolator old, enum copy copy)
{

    if (! _dxf_CopyFieldInterpolator((FieldInterpolator)new,
					(FieldInterpolator)old, copy))
	return NULL;
    
    new->nEdges    = old->nEdges;
    new->nFaces    = old->nFaces;
    new->nLoops    = old->nLoops;
    new->hint      = old->hint;

    if (new->fieldInterpolator.initialized)
    {
	new->fArray    = (Array)DXReference((Object)old->fArray);
	new->lArray    = (Array)DXReference((Object)old->lArray);
	new->eArray    = (Array)DXReference((Object)old->eArray);
	new->pArray    = (Array)DXReference((Object)old->pArray);
	new->dArray    = (Array)DXReference((Object)old->dArray);
	new->mmfArray  = (Array)DXReference((Object)old->mmfArray);
	new->mmlArray  = (Array)DXReference((Object)old->mmlArray);

	if (new->fieldInterpolator.localized)
	{
	    new->faces        = (int    *)DXGetArrayDataLocal(new->fArray);
	    new->loops        = (int    *)DXGetArrayDataLocal(new->lArray);
	    new->edges        = (int    *)DXGetArrayDataLocal(new->eArray);
	    new->positions    = (float  *)DXGetArrayDataLocal(new->pArray);
	    new->minmax_faces = (MinMax *)DXGetArrayDataLocal(new->mmfArray);
	    new->minmax_loops = (MinMax *)DXGetArrayDataLocal(new->mmlArray);
	}
	else
	{
	    new->faces        = (int    *)DXGetArrayData(new->fArray);
	    new->loops        = (int    *)DXGetArrayData(new->lArray);
	    new->edges        = (int    *)DXGetArrayData(new->eArray);
	    new->positions    = (float  *)DXGetArrayData(new->pArray);
	    new->minmax_faces = (MinMax *)DXGetArrayData(new->mmfArray);
	    new->minmax_loops = (MinMax *)DXGetArrayData(new->mmlArray);
	}

	new->dHandle = DXCreateArrayHandle(new->dArray);
	if (! new->dHandle)
	    return NULL;
    }
    else
    {
	new->fArray    = NULL;
	new->lArray    = NULL;
	new->eArray    = NULL;
	new->pArray    = NULL;
	new->dHandle   = NULL;
    }

    if (DXGetError())
	return NULL;

    return new;
}

Interpolator
_dxfFle2DInterpolator_LocalizeInterpolator(Fle2DInterpolator fle)
{
    if (fle->fieldInterpolator.localized)
	return (Interpolator)fle;

    fle->fieldInterpolator.localized = 1;

    if (fle->fieldInterpolator.initialized)
    {
        fle->faces        = (int    *)DXGetArrayDataLocal(fle->fArray);
        fle->loops        = (int    *)DXGetArrayDataLocal(fle->lArray);
        fle->edges        = (int    *)DXGetArrayDataLocal(fle->eArray);
        fle->positions    = (float  *)DXGetArrayDataLocal(fle->pArray);
        fle->minmax_faces = (MinMax *)DXGetArrayDataLocal(fle->mmfArray);
        fle->minmax_faces = (MinMax *)DXGetArrayDataLocal(fle->mmlArray);
    }

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