/***********************************************************************/
/* 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 <string.h>
#include "objectClass.h"
#include "internals.h"

/*
 * DXBoundingBox returns a bounding box of an object.  The box is not
 * necessarily aligned to world coordinates if a transform has been applied.
 *
 * The top level bounding box routine just calls _dxf_BoundingBox with
 * a null transform; the transform argument is for propagating the
 * transforms downward.
 */


#define VALID_POSITIONS_ONLY 1
#define ALL_POSITIONS        0

Object
DXBoundingBox(Object o, Point *box)
{
    if (DXGetError()!=ERROR_NONE)
	DXErrorReturn(ERROR_INTERNAL,
		    "DXBoundingBox called with error code set");
    return _dxfBoundingBox(o, box, NULL, ALL_POSITIONS);
}

Object
DXValidPositionsBoundingBox(Object o, Point *box)
{
    if (DXGetError()!=ERROR_NONE)
	DXErrorReturn(ERROR_INTERNAL,
		    "DXBoundingBox called with error code set");
    return _dxfBoundingBox(o, box, NULL, VALID_POSITIONS_ONLY);
}

/*
 * compute a bounding box from a set of points
 */

static
expand(float *out, int dim, int nboxpts, float *min, float *max)
{
    int i, j, k, m, n;
    n = dim * nboxpts;
    for (i=0, m=n/2; i<dim; i++, m/=2) {
	    for (j=0; j<n; j+=2*m) {
	    for (k=0; k<m; k+=dim) {
		out[i+j+k] = min[i];
		out[m+i+j+k] = max[i];
	    }
	}
    }
}


static Error
bbox(float *out, float *in, int dim, int nboxpts, int npoints,
	InvalidComponentHandle ich)
{
    float min[100], max[100];
    int i, j, k, m, n;

    /* callers should have checked this */
    ASSERT(npoints!=0);
    ASSERT(dim<100);

    /* initialize */
    for (i=0; i<dim; i++) {
        min[i] = DXD_MAX_FLOAT;
        max[i] = -DXD_MAX_FLOAT;
    }

    /* compute min/max */
    if (ich)
    {
	n = dim * npoints;
	for (i=k=0; i<n; k++, i+=dim) {
	    if (! DXIsElementInvalid(ich, k))
	    {
		for (j=0; j<dim; j++) {
		    if (in[i+j] < min[j])
			min[j] = in[i+j];
		    if (in[i+j] > max[j])
			max[j] = in[i+j];
		}
	    }
	}
    }
    else
    {
	n = dim * npoints;
	for (i=0; i<n; i+=dim) {
	    for (j=0; j<dim; j++) {
		if (in[i+j] < min[j])
		    min[j] = in[i+j];
		if (in[i+j] > max[j])
		    max[j] = in[i+j];
	    }
	}
    }

    /* expand min/max out to nboxpts corners */
    expand(out, dim, nboxpts, min, max);
}


/*
 * compute a bounding box from a group
 */

Object
_dxfGroup_BoundingBox(Group g, Point *box, Matrix *m, int validity)
{
    Object o;
    int i, j, k, n=0, align=0;
    Point min, max, p;
    Matrix mat, inv;

    /* compute the collective min/max of components */
    for (i=0; o=DXGetEnumeratedMember(g, i, NULL); i++) {

	/* get the component bounding box */
	if (!_dxfBoundingBox(o, box, m, validity)) {
	    if (DXGetError()!=ERROR_NONE)
		return NULL;		/* a real error occured */
	    continue;			/* or, the field had no bounding box */
	}

	/*
	 * Align group box to first box.  We do this by considering
	 * box[0] to be the origin of a coordinate system and box[1]-box[0],
	 * box[2]-box[0] and box[4]-box[0] to be the axes.  We transform
	 * all the bounding boxes into this coordinate system and form
	 * a bounding box there, and then transform that bounding box back
	 * to the original coordinate system when we're done.
	 *
	 * Modification: only do this if all dimensions are non-zero.
	 */
	if (box) {
	    if (n==0) {
		Vector v1, v2, v4;
		float l1, l2, l4;
		v1 = DXSub(box[1], box[0]), l1 = fabs(DXLength(v1));
		v2 = DXSub(box[2], box[0]), l2 = fabs(DXLength(v2));
		v4 = DXSub(box[4], box[0]), l4 = fabs(DXLength(v4));
		/* XXX - use an epsilon here */
		if (l1!=0.0 && l2!=0.0 && l4!=0.0) {
		    *(Point*) mat.A[0] = v1;
		    *(Point*) mat.A[1] = v2;
		    *(Point*) mat.A[2] = v4;
		    *(Point*) mat.b    = box[0];
		    inv = DXInvert(mat);
		    align = 1;
		}
		min.x = min.y = min.z = DXD_MAX_FLOAT;
		max.x = max.y = max.z = -DXD_MAX_FLOAT;
	    }
	    for (j=0; j<8; j++) {
		p = align? DXApply(box[j],inv) : box[j];
		if (p.x<min.x) min.x = p.x;  if (p.x>max.x) max.x = p.x;
		if (p.y<min.y) min.y = p.y;  if (p.y>max.y) max.y = p.y;
		if (p.z<min.z) min.z = p.z;  if (p.z>max.z) max.z = p.z;
	    }
	    n++;
	}
    }

    /*
     * DXTransform the bounding box back to the original
     * coordinate system (see above).
     */
    if (box) {
	if (n==0)
	    return NULL;
	for (i=0; i<8; i++) {
	    p = DXPt(i&1?max.x:min.x, i&2?max.y:min.y, i&4?max.z:min.z);
	    box[i] = align? DXApply(p, mat) : p;
	}
    }

    return (Object) g;
}

Object
_dxfCompositeField_BoundingBox(Group g, Point *box, Matrix *m, int validity)
{
    Object o;
    int i, j, k, n=0;
    Point min, max, *p;

    min.x = min.y = min.z = DXD_MAX_FLOAT;
    max.x = max.y = max.z = -DXD_MAX_FLOAT;

    /* compute the collective min/max of components */
    for (i=0; o=DXGetEnumeratedMember(g, i, NULL); i++) {

	/* get the component bounding box */
	if (!_dxfBoundingBox(o, box, NULL, validity)) {
	    if (DXGetError()!=ERROR_NONE)
		return NULL;		/* a real error occured */
	    continue;			/* or, the field had no bounding box */
	}

	/*
	 * Align group box to first box.  We do this by considering
	 * box[0] to be the origin of a coordinate system and box[1]-box[0],
	 * box[2]-box[0] and box[4]-box[0] to be the axes.  We transform
	 * all the bounding boxes into this coordinate system and form
	 * a bounding box there, and then transform that bounding box back
	 * to the original coordinate system when we're done.
	 *
	 * Modification: only do this if all dimensions are non-zero.
	 */
	if (box)
	{
	    for (j=0, p=box; j<8; j++, p++)
	    {
		if (p->x<min.x) min.x = p->x;  if (p->x>max.x) max.x = p->x;
		if (p->y<min.y) min.y = p->y;  if (p->y>max.y) max.y = p->y;
		if (p->z<min.z) min.z = p->z;  if (p->z>max.z) max.z = p->z;
	    }
	    n++;
	}
    }

    /*
     * DXTransform the bounding box back to the original
     * coordinate system (see above).
     */
    if (box) {
	if (n==0)
	    return NULL;
	for (i=0; i<8; i++)
	{
	    if (m)
	    {
		box[i] = DXApply(DXPt(i&1?max.x:min.x,
				      i&2?max.y:min.y,
				      i&4?max.z:min.z), *m);
	    }
	    else
	    {
		box[i] = DXPt(i&1?max.x:min.x,
			      i&2?max.y:min.y,
			      i&4?max.z:min.z);
	    }
	}
    }

    return (Object) g;
}


/*
 * compute a bounding box from a field
 */

static
Error
_copyout(Point *out, Matrix *m, Array a)
{
    int i, rank, dim;
    float *in;

    if (!DXGetArrayInfo(a, NULL, NULL, NULL, &rank, &dim))
	return NULL;
    if (rank==0)
	dim = 1;
    in = (float *) DXGetArrayData(a);
    if (!in)
	return NULL;

    if (dim==3) {
	if (m)
	    for (i=0; i<8; i++)
		out[i] = DXApply(((Point*)in)[i], *m);
	else
	    for (i=0; i<8; i++)
		out[i] = ((Point*)in)[i];
    } else if (dim==2) {
	if (m)
	    for (i=0; i<4; i++)
		out[2*i] = out[2*i+1] = DXApply(DXPt(in[2*i], in[2*i+1], 0), *m);
	else
	    for (i=0; i<4; i++)
		out[2*i] = out[2*i+1] = DXPt(in[2*i], in[2*i+1], 0);
    } else if (dim==1) {
	if (m)
	    for (i=0; i<2; i++)
		out[4*i] = out[4*i+1] = out[4*i+2] = out[4*i+3]
		    = DXApply(DXPt(in[i],0,0), *m);
	else
	    for (i=0; i<2; i++)
		out[4*i] = out[4*i+1] = out[4*i+2] = out[4*i+3]
		    = DXPt(in[i],0,0);
    } else {
	DXSetError(ERROR_BAD_PARAMETER, "positions must be 3-dimensional");
	return NULL;
    }

    return OK;
}


/*
 * Check or figure out dimensionality
 */

static
Error
check(Array a, int *dim, int *nboxpts)
{
    Type type;
    Category category;
    int rank, i;

    /* type check */
    if (!DXGetArrayInfo(a, NULL, &type, &category, &rank, NULL))
	return ERROR;
    if (type!=TYPE_FLOAT || category!=CATEGORY_REAL || (rank!=0 && rank!=1))
	DXErrorReturn(ERROR_BAD_TYPE, "bad positions/box type in bounding box");

    /* get info */
    if (rank==0)
	*dim = 1;
    else if (!DXGetArrayInfo(a, NULL, NULL, NULL, NULL, dim))
	return NULL;
    *nboxpts = 1;
    for (i=0; i<*dim; i++)
	*nboxpts *= 2;

    return OK;
}

/* XXX */
#define MAX_DIMS 25

Object
_dxfField_BoundingBox(Field f, Point *box, Matrix *m, int validity)
{
    Array points_array, box_array;
    float *in, *out;
    int i, j, k, dim, nboxpts, npoints, n, counts[MAX_DIMS];
    float origins[MAX_DIMS], deltas[MAX_DIMS*MAX_DIMS];
    float min[MAX_DIMS], max[MAX_DIMS];
    InvalidComponentHandle ich = NULL;
    int do_invalids = 0;

    /*
     * If its an empty field, then return NULL but no error.
     */
    if (DXEmptyField(f))
	goto null_return;

    /* points_array XXX return NULL if no points, but not an error */
    points_array = (Array) DXGetComponentValue(f, POSITIONS);
    if (!points_array)
	goto null_return;

    if (validity == VALID_POSITIONS_ONLY)
    {
	int invalids_present = (int)DXGetComponentValue(f, "invalid positions");

	/*
	 * if theres a valid_box component, or there are no invalids and
	 * there is a box component, then just return it, applying
	 * transform if any
	 */
	box_array = (Array) DXGetComponentValue(f, VBOX);
	if (! box_array && ! invalids_present)
	    box_array = (Array)DXGetComponentValue(f, BOX);

	if (box_array && box)
	    goto done;
    
	if (invalids_present)
	{
	    ich = DXCreateInvalidComponentHandle((Object)f, NULL, "positions");
	    if (! ich)
		goto null_return;

	    do_invalids = 1;
	}
	else
	{
	    validity = ALL_POSITIONS;
	    do_invalids = 0;
	}
    }
    else if (NULL != (box_array = (Array)DXGetComponentValue(f, BOX)))
	goto done;

    /* get dim, nboxpts */
    if (!check(points_array, &dim, &nboxpts))
	goto null_return;

    if (dim>MAX_DIMS)
    {
	DXSetError(ERROR_BAD_PARAMETER, "dimensionality of object is too high");
	goto null_return;
    }

    DXGetArrayInfo(points_array, &npoints, NULL, NULL, NULL, NULL);
    if (npoints==0 || (do_invalids && npoints == DXGetInvalidCount(ich)))
	goto null_return;

    /*
     * allocate array
     */
    box_array = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, dim);
    if (!DXAddArrayData(box_array, 0, nboxpts, NULL))
	goto null_return;

    if (do_invalids)
    {
	DXSetComponentValue(f, VBOX, (Object)box_array);
	DXSetAttribute((Object)box_array, DER, 
		  (Object)DXMakeStringList(2, POSITIONS, INVALID_POSITIONS));
    }
    else
    {
	DXSetComponentValue(f, BOX, (Object)box_array);
	DXSetAttribute((Object)box_array, DER, O_POSITIONS);
    }

    out = (float *) DXGetArrayData(box_array);
    if (!out)
	return ERROR;

    /* is it regular or irregular? */
    if (DXQueryGridPositions(points_array, NULL, counts, origins, deltas))
    {
	/*
	 * Check to see if the corners of the box are
	 * valid.  If so, we can just use the fast mechanism
	 */
	if (do_invalids)
	{
	    /* 
	     * quick out if the corners are valid
	     */
	    int index, i, j, strides[MAX_DIMS];

	    for (i = 0; i < dim; i++)
	    {
		strides[i] = counts[i] - 1;

		for (j = i+1; j < dim; j++)
		   strides[i] *= counts[j];
	    }
	    
	    do_invalids = 0;
	    for (i = 0; !do_invalids && i < 1<<dim; i++)
	    {
		for (j = index = 0; j < dim; j++)
		    if (i & (1<<j)) index += strides[j];
		
		do_invalids = DXIsElementInvalid(ich, index);
	    }
	}

	/*
	 * If all the corners were valid, just look at them
	 */
	if (! do_invalids)
	{
	    /*
	     * compute min/max
	     */
	    for (i=0; i<dim; i++) {
		float mn=0, mx=0, o;
		for (j=0; j<dim; j++) {
		    float m = deltas[j*dim+i] * (counts[j]-1);
		    if (m>0) mx += m;
		    else     mn += m;
		}
		o = origins[i];
		min[i] = o + mn;
		max[i] = o + mx;
	    }

	    /*
	     * expand min/max out to nboxpts corners
	     */
	    expand(out, dim, nboxpts, min, max);

	    goto done;
	}

    }
    
    if (DXGetArrayClass(points_array) == CLASS_ARRAY)
    {
	in = (float *) DXGetArrayData(points_array);
	if (!in)
	    return NULL;

	DXGetArrayInfo(points_array, &npoints, NULL, NULL, NULL, NULL);
	if (npoints==0)
	    return NULL;

	/* compute bounding box from points */
	bbox(out, in, dim, nboxpts, npoints, ich);
    }
    else
    {
	int	    i, j;
	float       min[MAX_DIMS], max[MAX_DIMS], *pt;
	ArrayHandle handle = DXCreateArrayHandle(points_array);
	Pointer     scratch = DXAllocate(DXGetItemSize(points_array));
	if (! handle || ! scratch)
	{
	    DXFreeArrayHandle(handle);
	    DXFree(scratch);
	    return ERROR;
	}

	DXGetArrayInfo(points_array, &npoints, NULL, NULL, NULL, NULL);
	if (npoints==0)
	{
	    DXFreeArrayHandle(handle);
	    DXFree(scratch);
	    return NULL;
	}

	for (i = 0; i < dim; i++)
	{
	    min[i] =  DXD_MAX_FLOAT;
	    max[i] = -DXD_MAX_FLOAT;
	}

	if (ich)
	    for (i = 0; i < npoints; i++)
	    {
		if (! DXIsElementInvalid(ich, i))
		{
		    pt = (float *)DXCalculateArrayEntry(handle, i, scratch);
		    for (j = 0; j < dim; j++)
		    {
			if (pt[j] < min[j]) min[j] = pt[j];
			if (pt[j] > max[j]) max[j] = pt[j];
		    }
		}
	    }
	else
	    for (i = 0; i < npoints; i++)
	    {
		pt = (float *)DXIterateArray(handle, i, (Pointer)pt, scratch);
		for (j = 0; j < dim; j++)
		{
		    if (pt[j] < min[j]) min[j] = pt[j];
		    if (pt[j] > max[j]) max[j] = pt[j];
		}
	    }
	

	DXFreeArrayHandle(handle);
	DXFree(scratch);

	/* expand min/max out to nboxpts corners */
	expand(out, dim, nboxpts, min, max);
    }


done:

    if (ich)
	DXFreeInvalidComponentHandle(ich);

    if (box && !_copyout(box, m, box_array))
	goto null_return;

    return (Object) f;

null_return:

    if (ich)
	DXFreeInvalidComponentHandle(ich);
    
    return NULL;
}

/*
 * bounding box of xform node
 */

Object
_dxfXform_BoundingBox(Xform x, Point *box, Matrix *mp, int validity)
{
    Matrix m;
    Object o;

    /* get this node's matrix and object */
    if (!DXGetXformInfo(x, &o, &m))
	return NULL;

    /* concatenate this node's matrix m with incoming  matrix *mp */
    if (mp)
	m = DXConcatenate(m, *mp);

    /* compute o's bounding box */
    if (!_dxfBoundingBox(o, box, &m, validity))
	return NULL;

    return (Object) x;
}


/*
 * Bounding box of clipped node.  Note: ignore bounding box of
 * clipping object;  because it is generally based on the bounding
 * box of the transformed clipped object and then has its bbox
 * recomputed by DXEndField, it will have an unnecessarily large
 * bounding box; but this is ok, because its bounding box is of no
 * interest anyway.
 */

Object
_dxfClipped_BoundingBox(Clipped c, Point *box, Matrix *mp, int validity)
{
    Object render;

    /* use bounding box of render object, ignoring clipping object */
    if (!DXGetClippedInfo(c, &render, NULL))
	return NULL;
    if (!_dxfBoundingBox(render, box, mp, validity))
	return NULL;

    return (Object) c;
}


/*
 * Screen objects are in completely different coordinates, so
 * we don't include them in DXBoundingBox.
 */

Object
_dxfScreen_BoundingBox(Screen s, Point *box, Matrix *mp, int validity)
{
    Object o;
    int fixed, i;

    if (!DXGetScreenInfo(s, &o, &fixed, NULL))
	return NULL;
    if (fixed==SCREEN_WORLD) {
	static Point zero = {0, 0, 0}, p;
	p = mp? DXApply(zero, *mp) : zero;
	if (box)
	    for (i=0; i<8; i++)
		box[i] = p;
	return (Object) s;
    } else
	return NULL;
}