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

typedef struct _gridlevel	*GridLevel;
typedef struct _item    	*Item;


Grid
_dxfMakeSearchGrid(float *max, float *min, int nPrims, int nDim)
{
    int 	i, j, k;
    float	span[MAXDIM], len, vol;
    int         counts[MAXDIM];
    Grid	grid;

    grid = (Grid)DXAllocate(sizeof(struct grid));
    if (! grid)
	goto error;

    /*
     * A few error conditions
     */
    if (nDim <= 0 || nDim > MAXDIM || nPrims == 0)
	goto error;

    grid->nDim = nDim;

    /*
     * Determine the "volume" of the bounding box to get a 
     * initial grid approximately the size of a primitive
     * From the average primitive "volume" guess an edge length
     * Note: add a bit to span so that we include the upper boundary.
     */
    vol = 1.0;
    for (i = 0; i < nDim; i++)
    {
	span[i] = 1.0001 * (max[i] - min[i]);
	vol *= span[i];
    }

    vol /= nPrims;
    len = (float)pow((double)vol, (double)(1.0/nDim));

    /*
     * Initial counts appropriate for grid[0].  Power of 2
     */
    for (i = 0; i < nDim; i++)
    {
	j = (span[i] / len) + 1;
	for (k = 1; k < j; k <<= 1);
	grid->counts[i] = counts[i] = k;
    }

    /*
     * _dxfInitialize grid levels.  As the grid level rises, the coarseness of
     * the grid increases: if the i-th level contains n bins along an axis,
     * the (i+1)-th will contain (n/2)+1.  Except the topmost, which contains
     * only one.
     */
    grid->itemArray   = NULL;
    for (i = 0; i < MAXGRID; i++)
    {
	grid->gridlevels[i].bucketArray = NULL;
	grid->gridlevels[i].current     = -1;
    }

    for (j = 0; j < nDim; j++)
    {
	grid->gridlevels[0].counts[j] = counts[j];
	grid->gridlevels[MAXGRID-1].counts[j] = 1;
    }

    for (i = 1; i < (MAXGRID - 1); i++)
	for (j = 0; j < nDim; j++)
	    grid->gridlevels[i].counts[j] = 
			(grid->gridlevels[i-1].counts[j] >> 1) + 1;


    /*
     * offset and scale for input points -> grid
     *		index = (input - min) * scale 
     * Include a fudge factor ensuring boundary conditions
     */
    for (i = 0; i < nDim; i++)
    {
	grid->min[i] = min[i];
	grid->scale[i] = 1.0 / span[i];
    }

    /*
     * Create allocation structure for bucket elements
     */
    grid->itemArray = DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 2);
    if (! grid->itemArray)
	goto error;

    DXReference((Object)grid->itemArray);

    if (! DXAddArrayData(grid->itemArray, 0, nPrims, NULL))
	goto error;
	
    grid->items = (struct _item *)DXGetArrayData(grid->itemArray);

    grid->next = 0;

    return grid;

error:
    if (grid)
    {
	for (i = 0; i < MAXGRID; i++)
	    DXDelete((Object)grid->gridlevels[i].bucketArray);

	DXDelete((Object)grid->itemArray);
    }

    return NULL;
}

void
_dxfFreeSearchGrid(Grid grid)
{
    GridLevel gl;
    int  i;

    if (grid)
    {
	for (i = 0, gl = grid->gridlevels; i < MAXGRID; i++, gl++)
	{
	    if (gl->bucketArray)
	    {
		DXDelete((Object)gl->bucketArray);
		gl->bucketArray = NULL;
	    }
	}

	if (grid->itemArray)
	{
	    DXDelete((Object)grid->itemArray);
	    grid->itemArray  = NULL;
	}

	DXFree((Pointer)grid);
    }
}

Grid
_dxfCopySearchGrid(Grid src)
{
    int  i;
    Grid dst;

    dst = (Grid)DXAllocate(sizeof(struct grid));
    if (! dst)
	return NULL;

    memcpy((char *)dst, (char *)src, sizeof(struct grid));

    for (i = 0; i < MAXGRID; i++)
	if (dst->gridlevels[i].bucketArray)
	    DXReference((Object)dst->gridlevels[i].bucketArray);

    DXReference((Object)dst->itemArray);

    return dst;
}

Error
_dxfAddItemToSearchGrid(Grid grid, float **vertices, int nVerts, int prim)
{
    int		 itemNum, bucketOffset;
    int		 level;
    GridLevel 	 gridlevel;
    Item	 itemPtr;
    register int i, j, nDim;
    int		 g;
    int		 same;
    int		 *stride;
    float 	 fMin[MAXDIM], fMax[MAXDIM];
    int		 iMin[MAXDIM], iMax[MAXDIM];

    nDim    = grid->nDim;

    /*
     * Get minmax of element in grid space
     */
    for (i = 0; i < nDim; i++)
    {
	fMax[i] = -DXD_MAX_FLOAT;
	fMin[i] =  DXD_MAX_FLOAT;
    }

    for (i = 0; i < nVerts; i++)
    {
	float *min   = grid->min;
	float *scale = grid->scale;

	for (j = 0; j < nDim; j++)
	{
	    float t = (vertices[i][j] - *min++) * *scale++;

	    if (t < fMin[j])
		fMin[j] = t;

	    if (t > fMax[j])
		fMax[j] = t;
	}
    }

    /*
     * For each grid level in turn, map the minmax into the integer index
     * space of the grid level.  Exit when the integer minmax falls into one
     * bin.  Note that the topmost level contains a single bin.
     */
    gridlevel = grid->gridlevels;
    for (level = 0; level < (MAXGRID - 1); level ++, gridlevel ++)
    {
	for (i = 0; i < nDim; i++)
	{
	    iMin[i] = (int)(fMin[i] * (gridlevel->counts[i] - 1));
	    if (iMin[i] < 0)
		iMin[i] = 0;

	    iMax[i] = (int)(fMax[i] * (gridlevel->counts[i] - 1));
	    if (iMax[i] >= gridlevel->counts[i])
		iMax[i] = gridlevel->counts[i] - 1;

	    if (iMin[i] != iMax[i])
		break;
	}

	if (i == nDim)
	    break;
    }

    /*
     * If we didn't succeed, stick it into the top bin
     */
    if (level == MAXGRID-1)
	for (i = 0; i < nDim; i++)
	    iMin[i] = 0;

    /*
     * If necessary, create the desired grid level.
     */
    if (gridlevel->bucketArray == NULL)
    {
	int nBins, *gc, *st, lc;

	/* 
	 * grid counts along the axes, and how many bins in grid?
	 * the following code is a bit obscure in order to avoid a
	 * optimizer bug on the 860
	 */
	gc = gridlevel->counts;
	st = gridlevel->stride;

	/*
	 * Initial setup: stride for first dimension is 1, number of 
	 * bins along that axis is the grid number reduced by g.
	 */
	lc = *gc++;
	nBins = lc;
	*st++ = 1;

	/*
	 * st and gc point to the second slots in the 
	 * counts and strides arrays respectively.  Note the 
	 * loop executes (nDim-1) times for the remaining axes.
	 */
	for (i = 1; i < nDim; i++)
	{
	    /*
	     * stride for the next axis is the length of the last
	     * axis times the stride for the last axis.
	     */
	    *st = *(st-1) * lc;

	    /*
	     * the number of bins in the current axis is...
	     */
	    lc = *gc;
	    
	    /*
	     * And multiply them in to the total bin count
	     */
	    nBins *= lc;

	    /*
	     * now bump the pointers
	     */
	    st ++;
	    gc ++;
	}

	/*
	 * DXAllocate grid's bins
	 */
	gridlevel->bucketArray = DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 1);
	if (! gridlevel->bucketArray)
	    goto error;
    
	DXReference((Object)gridlevel->bucketArray);
    
	if (! DXAddArrayData(gridlevel->bucketArray, 0, nBins, NULL))
	    goto error;

	gridlevel->buckets = (int *)DXGetArrayData(gridlevel->bucketArray);

	memset(gridlevel->buckets, -1, nBins * sizeof(int));

	gridlevel->count = 0;
    }

    /*
     * Now stick the item into that bucket in that grid.
     */
    stride    = gridlevel->stride;

    bucketOffset = iMin[0];
    for (i = 1; i < nDim; i++)
       bucketOffset += iMin[i] * stride[i];

    /*
     * Get a new item
     */
    itemNum = grid->next++;
    itemPtr = grid->items + itemNum;

    /*
     * Set up pointer to item
     */
    itemPtr->primNum = prim;

    /*
     * DXAdd to linked list
     */
    itemPtr->next = gridlevel->buckets[bucketOffset];
    gridlevel->buckets[bucketOffset] = itemNum;

    gridlevel->count ++;
		
    return 1;

error:
    return 0;
}

void
_dxfInitGetSearchGrid(Grid grid, float *pt) 
{ 
    int i;
    float f;

    for (i = 0; i < grid->nDim; i++)
    {
	f = (pt[i] - grid->min[i]) * grid->scale[i];
	if (f < 0)
	    f = 0.0;
	else if (f >= 1.0) 
	    f = 1.0;

	grid->gridPoint[i] = f;
    }

    grid->currentGrid = -1;
}

int
_dxfGetNextSearchGrid(Grid grid)
{
    Item	item;
    int		g, i, offset;
    GridLevel	gridlevel;
    int		index[MAXDIM];

    /*
     * Get pointer to current grid.  Note that the current grid
     * number is initialized to -1 for the search.
     */
    g = grid->currentGrid;
    gridlevel = grid->gridlevels + g;

    /*
     * If this is the initial call, or if we have exhausted the elements at
     * this grid level, move up.
     */
    if (g == -1 || gridlevel->current == -1)
    {
	do
	{
	    g++;
	    gridlevel ++;

	    /*
	     * If there ain't no more buckets, quit
	     */
	    if (g == MAXGRID)
		return -1;
	    
	    /*
	     * If this grid level is for real...
	     */
	    if (gridlevel->bucketArray)
	    {
		/*
		 * shift the indices to correspond to the larger 
		 * gridlevel
		 */
		for (i = 0; i < grid->nDim; i++)
		    index[i] = grid->gridPoint[i] * (gridlevel->counts[i] - 1);

		/*
		 * Get the bucket offset
		 */
		offset = index[0];
		for (i = 1; i < grid->nDim; i++)
		    offset += index[i] * gridlevel->stride[i];
	
		gridlevel->current = gridlevel->buckets[offset];
	    }
	}
	while(gridlevel->current == -1);

	grid->currentGrid = g;
    }

    if (gridlevel->current == -1)
    {
	return -1;
    }
    else
    {
	item = grid->items + gridlevel->current;
	gridlevel->current = item->next;
	return item->primNum;
    }

}

#if 0
#define SAMPLE_SIZE	32

Error
SizeEstimate(Field field, float *size)
{
    Array p, e;
    float *points;
    int *elems, *ePtr;
    int nElems, nDim, nVerts, nSamples;
    int i, j, k, random;
    float l, s;

    p = (Array)DXGetComponentValue(field, "positions");
    if (! p)
	return ERROR;

    e = (Array)DXGetComponentValue(field, "connections");
    if (! e)
	return ERROR;

    points = (float *)DXGetArrayData(p);
    DXGetArrayInfo(p, NULL, NULL, NULL, NULL, &nDim);

    elems = (int *)DXGetArrayData(e);
    DXGetArrayInfo(e, &nElems, NULL, NULL, NULL, &nVerts);

    if (nElems > SAMPLE_SIZE)
    {
	nSamples = SAMPLE_SIZE;
	random = 1;
    }
    else
    {
	nSamples = nElems;
	random = 0;
    }
    
    s = 999999.0;
    for (i = 0; i < nSamples; i++)
    {
 	float *p0, *p1, m;

	if (random)
	    k = rand() % nElems;
	else
	    k = i;
    
	ePtr = elems + k*nVerts;

	p0 = points + nDim*ePtr[0];
	p1 = points + nDim*ePtr[nVerts-1];

	m = 0.0;
	for (j = 0; j < nDim; j++)
	{
	    l = *p0++ - *p1++;
	    if (l < 0.0)
		m += -l;
	    else
		m += l;
	}

	if (m < s)
	    s = m;
    }

    *size = s;

    return OK;
}
#endif