/***********************************************************************/
/* 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 <math.h>
#include "cameraClass.h"
#include "internals.h"
#include "render.h"

static void doboxpts(Point *pts, int ortho,
		     float minx, float miny, float minz,
		     float max, float maxy, float maxz,
		     int *behind);
Array
_dxf_TransformArray(Matrix *, Camera, Array, Array *, Field, int *, float);


Array
_dxf_TransformBox(Matrix *, Camera, Array, int *, float);

/*
 * DXTransform a regular array.  This consists of applying the full
 * transformation to the origin and all but the translation part to the
 * deltas: if the transform is xA + b$ and each point is of the form o
 * + id$, then each transformed point is of the form (oA+b)+(idA).
 * This is only applicable if there is no perspective transformation
 * involved.
 */

static Array
TransformRegularArray(Matrix m, Camera c, Array points, Array *box)
{
    int n;
    Point origin, delta;
    Array xpoints;

    /* is regular array transformation applicable? */
    if ((c && !c->ortho) || DXGetArrayClass(points)!=CLASS_REGULARARRAY)
	return NULL;

    /* better to transform existing box in regular case */
    if (box)
	*box = NULL;

    /* type check */
    if (!DXTypeCheck(points, TYPE_FLOAT, CATEGORY_REAL, 1, 3) &&
        !DXTypeCheck(points, TYPE_FLOAT, CATEGORY_REAL, 1, 2) &&
        !DXTypeCheck(points, TYPE_FLOAT, CATEGORY_REAL, 1, 1))
	DXErrorReturn(ERROR_BAD_TYPE,
	    "positions data must be 1, 2, or 3-D float vectors");

    /* set z to 0 in case this is 2-d or 1-d */
    origin = delta = DXPt(0,0,0);

    /* get regular array info */
    if (!DXGetRegularArrayInfo((RegularArray)points, &n,
			     (Pointer)&origin, (Pointer)&delta))
	return NULL;

    /* apply full transformation to the origin */
    origin = DXApply(origin, m);

    /* apply transformation less translation part to the delta */
    m.b[0] = m.b[1] = m.b[2] = 0;
    delta = DXApply(delta, m);

    /* construct the regular array result */
    xpoints = (Array) DXNewRegularArray(TYPE_FLOAT, 3, n,
				      (Pointer)&origin, (Pointer)&delta);

    return xpoints;
}


/*
 * DXTransform a product geometry array.  This consists of applying the
 * full transformation to the first term and all but the translation part
 * to the remaining terms: if the transformation is xA+b$ and each point
 * is of the form t_0+t_1+t_2+\cdots$, where each t_i$ is a point from
 * a different term, then the transformed points are of the form
 * (t_0A+b)+(t_1A)+(t_2A)+\cdots$.  This is only applicable if there is
 * no perspective transformation involved.
 *
 * Note the hack involving FUDGE: for regular partitioned grids that
 * exactly coincide with the pixel grid, this method of computing the
 * transform can introduce roundoff errors causing origin+n*deltax
 * of one partition to overlap the origin of the next partition, causing
 * pixels to appear to be in more than one partition, resulting in
 * volume rendering artifacts.  We "solve" this by moving the grid a
 * tiny bit, thus avoiding exact pixel alignment of the grid.  Note that
 * all this does is make the artifact much less likely; it doesn't
 * eliminate this possibility.
 */

#define FUDGE (1.0/93.0)

static Array
TransformProductArray(Matrix m, Camera c, Array points, Array *box,
		      Field f, int *behind)
{
    int i, n;
    Array xpoints, terms[100];

    /* is product array transformation applicable? */
    if ((c && !c->ortho) || DXGetArrayClass(points)!=CLASS_PRODUCTARRAY)
	return NULL;

    /* type check */
    if (!DXTypeCheck(points, TYPE_FLOAT, CATEGORY_REAL, 1, 3) &&
        !DXTypeCheck(points, TYPE_FLOAT, CATEGORY_REAL, 1, 2) &&
        !DXTypeCheck(points, TYPE_FLOAT, CATEGORY_REAL, 1, 1))
	DXErrorReturn(ERROR_BAD_TYPE,
	    "positions data must be 1, 2, or 3-D float vectors");
    /* get terms */
    if (!DXGetProductArrayInfo((ProductArray)points, &n, terms))
	return NULL;
    ASSERT(n<100);

    /* apply whole transformation to first term */
    if (c) {				/* fudge only final camera transform */
	m.b[0] += FUDGE;		/* fudge the x - see comment above */
	m.b[1] += FUDGE;		/* fudge the y - see comment above */
    }
    terms[0] = _dxf_TransformArray(&m, NULL, terms[0], NULL, f, NULL, 0.0);
    if (!terms[0])
	return NULL;

    /* transform the box with the fudged matrix */
    if (box) {
	Array b = (Array) DXGetComponentValue(f, BOX);
	*box = _dxf_TransformBox(&m, NULL, b, behind, 0.0);
	if (!*box)
	    return NULL;
    }

    /* now zero the translation and apply rotation part to rest of terms */
    m.b[0] = m.b[1] = m.b[2] = 0;
    for (i=1; i<n; i++) {
	terms[i] = _dxf_TransformArray(&m, NULL, terms[i], NULL, f, NULL, 0.0);
	if (!terms[i])
	    return NULL;
    }

    /* construct product array result */
    return (Array) DXNewProductArrayV(n, terms);
}


/*
 * Transforms the input array of points or normals by transformation
 * {\tt t} followed by the camera transformation specified by {\tt c}.
 * The transformation is a point or normal transformation according to the
 * setting of the {\tt normals} flag.
 * XXX - vectorize this
 * XXX - handle different data types
*/

#define XFORM {							\
    ox = mp->A[0][0]*ix+mp->A[1][0]*iy+mp->A[2][0]*iz+mp->b[0];	\
    oy = mp->A[0][1]*ix+mp->A[1][1]*iy+mp->A[2][1]*iz+mp->b[1];	\
    oz = mp->A[0][2]*ix+mp->A[1][2]*iy+mp->A[2][2]*iz+mp->b[2];	\
    out[i].x = ox, out[i].y = oy, out[i].z = oz;		\
}

#define BBOX {							\
    if (ox > maxx) maxx = ox; if (ox < minx) minx = ox;		\
    if (oy > maxy) maxy = oy; if (oy < miny) miny = oy;		\
}

#define PBOX {							\
    if (ox > maxx) maxx = ox; if (ox < minx) minx = ox;		\
    if (oy > maxy) maxy = oy; if (oy < miny) miny = oy;		\
    if (oz > maxz) maxz = oz; if (oz < minz) minz = oz;		\
}

static Array
TransformIrregular(Matrix *mp, Camera c, Array points, Point *box, int *behind)
{
    int i, n, rank, shape[100], dim;
    Array xpoints;
    Point *in, *out;
    float ix, iy, iz, ox, oy, oz, minx, miny, minz, maxx, maxy, maxz;

    /* get input data to local memory */
    in = (Point *) DXGetArrayData(points);
    if (!DXGetArrayInfo(points, &n, NULL, NULL, &rank, shape))
	DXErrorReturn(ERROR_INTERNAL, "DXGetArrayInfo failed!");
    if (rank!=1) {
	DXSetError(ERROR_BAD_PARAMETER,
	    "positions data must be 1, 2, or 3-D float vectors (not scalar)");
	return NULL;
    }
    dim = shape[0];

    /* create output array, allocate storage for it */
    xpoints = DXAllocateArray(DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3), n);
    out = (Point *) DXGetArrayData(DXAddArrayData(xpoints, 0, n, NULL));
    if (!out)
	return NULL;

    if (box) {
	minx = miny = minz = DXD_MAX_FLOAT;
	maxx = maxy = maxz = -DXD_MAX_FLOAT;
    }

    /* transform points */
    if (dim==3) {
	if (box) {
	    if (!c || c->ortho) {
		for (i=0; i<n; i++) {
		    ix=in[i].x, iy=in[i].y, iz=in[i].z;
		    XFORM;
		    BBOX;
		}
	    } else {
		for (i=0; i<n; i++) {
		    ix=in[i].x, iy=in[i].y, iz=in[i].z;
		    XFORM;
		    PBOX;
		}
	    }
	} else {
	    for (i=0; i<n; i++) {
		ix=in[i].x, iy=in[i].y, iz=in[i].z;
		XFORM;
	    }
	}
    } else if (dim==2) {
	struct in2 {float x, y;} *in2 = (struct in2 *)in;
	for (i=0; i<n; i++) {
	    ix=in2[i].x, iy=in2[i].y, iz=0;
	    XFORM;
	    if (box) {
		PBOX;
	    }
	}
    } else if (dim==1) {
	struct in1 {float x;} *in1 = (struct in1 *)in;
	for (i=0; i<n; i++) {
	    ix=in1[i].x, iy=0, iz=0;
	    XFORM;
	    if (box) {
		PBOX;
	    }
	}
    } else {
	DXSetError(ERROR_BAD_PARAMETER,
		 "array is %d-d, expecting 1-d, 2-d, or 3-d", dim);
	return NULL;
    }

    /* expand to eight points */
    if (box)
	doboxpts(box, !c || c->ortho,
		 minx, miny, minz, maxx, maxy, maxz, behind);

    return xpoints;
}

Array
_dxf_TransformNormals(Matrix *mp, Array points)
{
    int i, n, dim, bad = 0;
    Array xpoints;
    Point *in, *out;
    Matrix m;
    Point cstIn, cstOut;

    /* normals unchanged? */
    if (!mp)
	return points;

    /* type check */
    if (!DXTypeCheck(points, TYPE_FLOAT, CATEGORY_REAL, 1, 3))
	DXErrorReturn(ERROR_BAD_PARAMETER, "normals not 3-d");

    /* matrix is different for normals */
    m = DXAdjointTranspose(*mp);
    
    if (DXQueryConstantArray(points, &n, &cstIn))
    {
	float x, y, z, d;
	x=m.A[0][0]*cstIn.x+m.A[1][0]*cstIn.y+m.A[2][0]*cstIn.z+m.b[0];
	y=m.A[0][1]*cstIn.x+m.A[1][1]*cstIn.y+m.A[2][1]*cstIn.z+m.b[1];
	z=m.A[0][2]*cstIn.x+m.A[1][2]*cstIn.y+m.A[2][2]*cstIn.z+m.b[2];
	d = sqrt(x*x+y*y+z*z);
	if (d) d = (float)1.0 / d;
	else bad = n;
	cstOut.x = x*d;
	cstOut.y = y*d;
	cstOut.z = z*d;

	xpoints = (Array)DXNewConstantArray(n, (Pointer)&cstOut,
				TYPE_FLOAT, CATEGORY_REAL, 1, 3);
	if (! xpoints)
	    return NULL;
    }
    else 
    {
	/* get input data */
	in = (Point *) DXGetArrayData(points);
	if (!DXGetArrayInfo(points, &n, NULL, NULL, NULL, &dim))
	    DXErrorReturn(ERROR_INTERNAL, "DXGetArrayInfo failed!");

	/* create output array, allocate storage for it */
	xpoints = DXAllocateArray(DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3), n);
	out = (Point *) DXGetArrayData(DXAddArrayData(xpoints, 0, n, NULL));
	if (!out)
	    return NULL;


	/* transform points */
	for (i=0; i<n; i++) {
	    float x, y, z, d;
	    x=m.A[0][0]*in[i].x+m.A[1][0]*in[i].y+m.A[2][0]*in[i].z+m.b[0];
	    y=m.A[0][1]*in[i].x+m.A[1][1]*in[i].y+m.A[2][1]*in[i].z+m.b[1];
	    z=m.A[0][2]*in[i].x+m.A[1][2]*in[i].y+m.A[2][2]*in[i].z+m.b[2];
	    d = sqrt(x*x+y*y+z*z);
	    if (d) d = (float)1.0 / d;
	    else bad++;
	    out[i].x = x*d;
	    out[i].y = y*d;
	    out[i].z = z*d;
	}
    }

    if (bad)
	DXWarning("%d zero-length %s detected",
		bad, bad>1? "normals were" : "normal was");

    return xpoints;
}

/*
 * Transforms the input array of points or normals by transformation
 * {\tt t} followed by the camera transformation specified by {\tt c}.
 * The transformation is a point or normal transformation according to the
 * setting of the {\tt normals} flag.
 * Returns a pointer to a new array holding the transformed points or null to
 * indicate an error.
 */

Array
_dxf_TransformArray(Matrix *mp, Camera c, Array points,
		Array *box, Field f, int *behind, float fuzz)
{
    Array xpoints, a;
    Matrix m;
    Point *boxpts;

    /*
     * if both mp and c are null, the transform requested
     * is identity, so return without doing anything.
     * Because of this, you can call DXTransform(o, NULL, NULL)
     * to make sure everything is in world coordinates as often as
     * you want without penalty.
     */
    if (!mp && !c && DXTypeCheck(points,TYPE_FLOAT,CATEGORY_REAL,1,3))
	return points;

    /* get the transform */
    if (mp)
	m = DXConcatenate(*mp, DXGetCameraMatrixWithFuzz(c, fuzz));
    else
	m = DXGetCameraMatrixWithFuzz(c, fuzz);

    /* regular array? */
    if (xpoints = TransformRegularArray(m, c, points, box))
	return xpoints;

    /* product array? */
    if (xpoints = TransformProductArray(m, c, points, box, f, behind))
	return xpoints;

    /* irregular array */
    else {
	if (box) {
	    *box = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
	    if (!*box) return NULL;
	    boxpts = (Point *) DXGetArrayData(DXAddArrayData(*box, 0, 8, NULL));
	    if (!boxpts) return NULL;
	} else
	    boxpts = NULL;
	return TransformIrregular(&m, c, points, boxpts, behind);
    }
}

/*
 * DXTransform a box, rectifying to the new coordinate system,
 * and expanding to 8 points if necessary
 */


Array
_dxf_TransformBox(Matrix *mp, Camera c, Array points, int *behind, float fuzz)
{
    float minx, miny, minz, maxx, maxy, maxz;
    Point *pts;
    Array box;
    int i, n;

    /* transform */
    box = _dxf_TransformArray(mp, c, points, NULL, NULL, NULL, fuzz);
    if (!box) return NULL;
    DXGetArrayInfo(box, &n, NULL, NULL, NULL, NULL);
    pts = (Point *) DXGetArrayData(box);
    if (!pts) return NULL;

    /* compute min/max */
    minx = miny = minz = DXD_MAX_FLOAT;
    maxx = maxy = maxz = -DXD_MAX_FLOAT;
    for (i=0; i<n; i++) {
	if (pts[i].x < minx) minx = pts[i].x;
	if (pts[i].y < miny) miny = pts[i].y;
	if (pts[i].z < minz) minz = pts[i].z;
	if (pts[i].x > maxx) maxx = pts[i].x;
	if (pts[i].y > maxy) maxy = pts[i].y;
	if (pts[i].z > maxz) maxz = pts[i].z;
    }

    /* expand to 8 points if necessary */
    if (n!=8) {
	Array a = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
	if (!DXAddArrayData(a, 0, 8, NULL)) return NULL;
	DXDelete((Object)box);
	box = a;
    }	

    /* rectify */
    pts = (Point *) DXGetArrayData(box);
    if (!pts) return NULL;
    doboxpts(pts, !c || c->ortho, minx, miny, minz, maxx, maxy, maxz, behind);

    return box;
}


/*
 * While the perspective divide for positions is deferred until
 * the object is tiled to allow clipping, the perspective
 * divide for the bounding box is done during the transformation
 * stage.  This routine computes the screen-aligned bounding box,
 * taking into account perspective.
 *
 * The result that we leave in the "bounding box" is rather strange,
 * in that 1) the min/max in x have had the perspective divide done
 * while the min/max in z have not, and 2) we clip the bounding box
 * at the nearPlane plane to be only the portion BEHIND the viewer; this
 * means box[0].z is near, and box[7].z can be used to determine whether
 * a clipping pass is needed.
 *
 * MAX_XY is the maximum x,y that we will allow (by choice of near)
 * after the perspective divide.  This must be large enough to
 * allow nearby triangles that are partly in front of and partly
 * behind the camera to remain on screen as long as possible, but
 * not so large as to allow overflow (in float->int conversion) or
 * insufficient precision in the triangle routines.
 *
 * We also limit near by SMALL, because even if x,y are uniformly
 * 0, we can't allow 1/near to overflow.
 */

#define MAX_XY 100000.0
#define SMALL 1e-35

static void
doboxpts(Point *pts, int ortho,
	 float minx, float miny, float minz,
	 float maxx, float maxy, float maxz,
	 int *behind
) {
    float nearPlane;

    if (!ortho) {

	float xy, max = 0;
	float mxz, mnz;

	/* extreme x,y */
	xy = minx;  if (xy<0) xy = -xy;  if (xy>max) max = xy;
	xy = maxx;  if (xy<0) xy = -xy;  if (xy>max) max = xy;
	xy = miny;  if (xy<0) xy = -xy;  if (xy>max) max = xy;
	xy = maxy;  if (xy<0) xy = -xy;  if (xy>max) max = xy;

	/* compute nearPlane to avoid overflow and loss of precision */
	nearPlane = max / MAX_XY;
	if (nearPlane<SMALL) nearPlane = SMALL;
	nearPlane = -nearPlane;

	/* x,y bounds after perspective divide */
	mxz = maxz>nearPlane? nearPlane : maxz;
	mnz = minz>nearPlane? nearPlane : minz;
	minx /= (minx<0? -mxz : -mnz);
	maxx /= (maxx>0? -mxz : -mnz);
	miny /= (miny<0? -mxz : -mnz);
	maxy /= (maxy>0? -mxz : -mnz);

	/* is object behind camera? */
	if (behind && minz>nearPlane)
	    *behind = 1;

    } else
	nearPlane = 0;

    /* generate the bounding box */
    if (maxz<nearPlane) maxz = nearPlane;
    pts[0].x = minx,  pts[0].y = miny,  pts[0].z = nearPlane;
    pts[1].x = minx,  pts[1].y = miny,  pts[1].z = maxz;
    pts[2].x = minx,  pts[2].y = maxy,  pts[2].z = nearPlane;
    pts[3].x = minx,  pts[3].y = maxy,  pts[3].z = maxz;
    pts[4].x = maxx,  pts[4].y = miny,  pts[4].z = nearPlane;
    pts[5].x = maxx,  pts[5].y = miny,  pts[5].z = maxz;
    pts[6].x = maxx,  pts[6].y = maxy,  pts[6].z = nearPlane;
    pts[7].x = maxx,  pts[7].y = maxy,  pts[7].z = maxz;

}

static Array ApplyTransform_GenericArray(Array, Matrix *, int);
static Array ApplyTransform_RegularArray(Array, Matrix *, int);
static Array ApplyTransform_ProductArray(Array, Matrix *, int);
static Array ApplyTransform_ConstantArray(Array, Matrix *, int);

static Field ApplyTransform_Field(Field, Matrix *);
static Array ApplyTransform_Array(Array, Matrix *, int);

static Object AT_traverse(Object, Matrix *);
static Error  AT_field_task(Pointer);

static Array SetDimensionality(Array, int);

Object
DXApplyTransform(Object object, Matrix *stack)
{
    if (! DXCreateTaskGroup())
	goto error;
    
    object = AT_traverse(object, stack);
    if (! object)
    {
	DXAbortTaskGroup();
	goto error;
    }

    if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE)
	goto error;
    
    return object;

error:
    return NULL;
}

static Object
AT_traverse(Object object, Matrix *stack)
{
    switch (DXGetObjectClass(object))
    {
	case CLASS_ARRAY:
	{
	    Object attr = DXGetAttribute(object, "geometric");
	    if (attr)
	    {
		if (! strcmp(DXGetString((String)attr), "vector"))
		{
		    Matrix ax;
		    ax = DXAdjointTranspose(*stack);
		    return (Object)ApplyTransform_Array((Array)object, &ax, 0);
		}
		else
		    return (Object)ApplyTransform_Array((Array)object,stack,0);
	    }
	    else
		return object;
	}

	case CLASS_FIELD:
	    return (Object)ApplyTransform_Field((Field)object, stack);
	
	case CLASS_XFORM:
	{
	    Matrix newmatrix, product;
	    Object child;

	    if (! DXGetXformInfo((Xform)object, &child, &newmatrix))
		return NULL;
	    
	    if (stack)
	    {
		product = DXConcatenate(newmatrix, *stack);
		return AT_traverse(child, &product);
	    }
	    else
		return AT_traverse(child, &newmatrix);
	}

	case CLASS_CLIPPED:
	{
	    Object clipper, clipped;

	    if (! DXGetClippedInfo((Clipped)object, &clipped, &clipper))
		return NULL;
	    
	    if (clipped && NULL == (clipped = AT_traverse(clipped, stack)))
		return NULL;

	    if (clipper && NULL == (clipper = AT_traverse(clipper, stack)))
	    {
		DXDelete(clipped);
		return NULL;
	    }

	    return (Object)DXNewClipped(clipped, clipper);
	}

	case CLASS_GROUP:
	{
	    int i;
	    Object child;
	    Group  result;

	    result = (Group)DXCopy(object, COPY_HEADER);

	    i = 0;
	    while (NULL !=
		(child = DXGetEnumeratedMember((Group)object, i, NULL)))
	    {
		child = AT_traverse(child, stack);
		if (! child)
		{
		    DXDelete((Object)result);
		    return NULL;
		}
		if (! DXSetEnumeratedMember(result, i++, child))
		{
		    DXDelete((Object)result);
		    return NULL;
		}
	    }

	    return (Object)result;
	}

	default:
	    return DXCopy(object, COPY_STRUCTURE);
    }
}

#define NOT_GEOMETRIC		0
#define GEOMETRIC_VECTORS	1
#define GEOMETRIC_POSITIONS	2

typedef struct {
    Field field;
    Matrix matrix;
} AT_field_args; 

static Field
ApplyTransform_Field(Field field, Matrix *matrix)
{
    AT_field_args args;

    field = (Field)DXCopy((Object)field, COPY_STRUCTURE);
    if (! field) 
	return NULL;

    if (! matrix)
	return field;

    args.field = field;
    args.matrix = *matrix;

    if (! DXAddTask(AT_field_task, (Pointer)&args, sizeof(args), 1.0))
	goto error;
    
    return field;

error:
    DXDelete((Object)field);
    return NULL;
}

static Error
AT_field_task(Pointer ptr)
{
    Array a;
    int   i, doit;
    char  *name;
    Matrix adjxpose;
    int    adjxpose_done = 0;
    Field  field;
    Matrix *matrix;

    field  = ((AT_field_args *)ptr)->field;
    matrix = &((AT_field_args *)ptr)->matrix;

    i = 0;
    while (NULL != (a = (Array)DXGetEnumeratedComponentValue(field, i++, &name)))
    {
	doit = NOT_GEOMETRIC;

	if (! strcmp(name, "positions"))
	{
	    a = (Array)ApplyTransform_Array(a, matrix, 0);
	}
	else if (! strcmp(name, "normals")   ||
	         ! strcmp(name, "gradient")  ||
	         ! strcmp(name, "binormals") ||
	         ! strcmp(name, "tangents"))
	{
	    if (! adjxpose_done)
	    {
		adjxpose = DXAdjointTranspose(*matrix);
		adjxpose_done = 1;
	    }

	    a = (Array)ApplyTransform_Array(a, &adjxpose, 1);
	}
	else
	{
	    Object attr = DXGetComponentAttribute(field, name, "geometric");
	    if (attr)
	    {
		if (! strcmp(DXGetString((String)attr), "positional"))
		    a = ApplyTransform_Array(a, matrix, 0);
		else if (! strcmp(DXGetString((String)attr), "vector"))
		{
		    if (! adjxpose_done)
		    {
			adjxpose = DXAdjointTranspose(*matrix);
			adjxpose_done = 1;
		    }

		    a = ApplyTransform_Array(a, &adjxpose, 0);
		}
	    }
	}

	if (! a)
	    goto error;

	if (! DXSetComponentValue(field, name, (Object)a))
	    goto error;
    }

    DXDeleteComponent(field, "box");
    if (! DXEndField(field))
	goto error;

    return OK;

error:
    return ERROR;
}

static Array
ApplyTransform_Array(Array array, Matrix *m, int normalize)
{
    Array a;
    switch (DXGetArrayClass(array))
    {
	case CLASS_REGULARARRAY:
	    if (DXQueryConstantArray(array, NULL, NULL))
		a = ApplyTransform_ConstantArray(array, m, normalize);
	    else
		a = ApplyTransform_RegularArray(array, m, normalize);
	    break;
	case CLASS_PRODUCTARRAY:
	    a = ApplyTransform_ProductArray(array, m, normalize);
	    break;
	case CLASS_CONSTANTARRAY:
	    a = ApplyTransform_ConstantArray(array, m, normalize);
	    break;
	case CLASS_ARRAY:
	    a = ApplyTransform_GenericArray(array, m, normalize);
	    break;
	default:
	    DXSetError(ERROR_INVALID_DATA, 
			"Inappropriate array class for transformation");
	    return NULL;
    }

    return a;
}
	    
static Array
ApplyTransform_ProductArray(Array array, Matrix *m, int norm)
{
    ProductArray result = NULL;
    int          i, n;
    Array        terms[100];
    Matrix	 m0;
    int		 dim;

    if (! DXGetProductArrayInfo((ProductArray)array, &n, terms))
	goto error;
    
    if (NULL == 
	(terms[0] = (Array)ApplyTransform_Array((Array)terms[0], m, norm)))
	goto error;

    if (n > 1)
    {
	memcpy((char *)&m0, (char *)m, sizeof(Matrix));
	m0.b[0] = m0.b[1] = m0.b[2] = 0.0;

	for (i = 1; i < n; i++)
	{
	    terms[i] = (Array)ApplyTransform_Array((Array)terms[i], &m0, norm);
	    if (! terms[i])
	    {
		for (i--; i >= 0; i--)
		    DXDelete((Object)terms[i]);
		goto error;
	    }
	}
    }

    /*
     * Find max. dimensionality term and pad out remainder to match
     */
    dim = 0;
    for (i = 0; i < n; i++)
    {   
	int d;

	DXGetArrayInfo(terms[i], NULL, NULL, NULL, NULL, &d);

	if (d > dim)
	    dim = d;
    }

    for (i = 0; i < n; i++)
    {  
	Array new = SetDimensionality(terms[i], dim);

	if (! new)
	    goto error;

	if (new != terms[i])
	{
	    DXDelete((Object)terms[i]);
	    terms[i] = new;
	}
    }

    result = DXNewProductArrayV(n, terms);
    if (! result)
    {
	for (i--; i >= 0; i--)
	    DXDelete((Object)terms[i]);
	goto error;
    }

    return (Array)result;

error:
    DXDelete((Object)result);
    return NULL;
}

#define TRANSFORM_REGULARARRAY(type, round) 		\
{							\
    type *io = (type *)inOrigin;			\
    type *id = (type *)inDeltas;			\
    type *oo = (type *)outOrigin;			\
    type *od = (type *)outDeltas;			\
    int  i, j;						\
							\
    for (i = 0; i < nDout; i++)				\
    {							\
	double fo, fd;					\
							\
	fd = 0.0;					\
	fo = matrix->b[i];				\
							\
	for (j = 0; j < nDin; j++)			\
	{						\
	    fo += io[j]*matrix->A[j][i];		\
	    fd += id[j]*matrix->A[j][i];		\
	}						\
							\
	oo[i] = (type)(fo + round);			\
	od[i] = (type)(fd + round);			\
    }							\
}

static Array
ApplyTransform_RegularArray(Array array, Matrix *matrix, int norm)
{
    Type     type;
    Category cat;
    Pointer  inOrigin = NULL, outOrigin = NULL;
    Pointer  inDeltas = NULL, outDeltas = NULL;
    int      rank, shape[32], size, n, count, nDin, nDout, i, j;
    Array    result = NULL;

    if (! DXGetArrayInfo(array, &count, &type, &cat, &rank, shape))
	return NULL;
    
    if (rank != 1 || shape[0] > 3)
    {
	DXSetError(ERROR_INVALID_DATA, 
	    "data to be transformed must be vectors of dimensionality <= 3");
	return NULL;
    }

    if (norm)
    {
	DXSetError(ERROR_INVALID_DATA, 
	    "cannot normalize transformed non-constant regular vectors");
	return NULL;
    }

    nDin = shape[0];

    for (nDout = 3; nDout >= 1; nDout --)
    {
	int zero = (matrix->b[nDout-1] == 0.0);
	
	for (i = 0; i < nDin && zero; i++)
	    zero = (matrix->A[i][nDout-1] == 0.0);

	if (! zero)
	    break;
    }

    inOrigin  = DXAllocate(nDin*DXTypeSize(type));
    inDeltas  = DXAllocate(nDin*DXTypeSize(type));
    outOrigin = DXAllocate(nDout*DXTypeSize(type));
    outDeltas = DXAllocate(nDout*DXTypeSize(type));

    if (! inOrigin || ! inDeltas || ! outOrigin || ! outDeltas)
    {
	DXFree(inOrigin);
	DXFree(inDeltas);
	DXFree(outOrigin);
	DXFree(outDeltas);
	return NULL;
    }

    DXGetRegularArrayInfo((RegularArray)array, &n, inOrigin, inDeltas);

    switch(type)
    {
	case TYPE_DOUBLE:
		TRANSFORM_REGULARARRAY(double, 0.0);
		break;

	case TYPE_FLOAT:
		TRANSFORM_REGULARARRAY(float, 0.0);
		break;

	case TYPE_INT:
		TRANSFORM_REGULARARRAY(int, 0.5);
		break;

	case TYPE_UINT:
		TRANSFORM_REGULARARRAY(uint, 0.5);
		break;

	case TYPE_SHORT:
		TRANSFORM_REGULARARRAY(short, 0.5);
		break;

	case TYPE_USHORT:
		TRANSFORM_REGULARARRAY(ushort, 0.5);
		break;

	case TYPE_BYTE:
		TRANSFORM_REGULARARRAY(byte, 0.5);
		break;

	case TYPE_UBYTE:
		TRANSFORM_REGULARARRAY(ubyte, 0.5);
		break;
    }

    result=(Array)DXNewRegularArray(type, nDout, n, outOrigin, outDeltas);

    DXFree(outOrigin);
    DXFree(outDeltas);
    DXFree(inOrigin);
    DXFree(inDeltas);

    return result;
}

#define TRANSFORM_CONSTANTARRAY(type, round) 		\
{							\
    type *io = (type *)inVector;			\
    type *oo = (type *)outVector;			\
    int  i, j;						\
							\
    for (i = 0; i < nDout; i++)				\
    {							\
	double fo;					\
							\
	fo = matrix->b[i];				\
							\
	for (j = 0; j < nDin; j++)			\
	    fo += io[j]*matrix->A[j][i];		\
							\
	oo[i] = (type)(fo + round);			\
    }							\
}

#define NORMALIZE_CONSTANTARRAY(type)			\
{							\
    type *oo = (type *)outVector;			\
    double foo = 0.0;					\
    for (i = 0; i < nDout; i++) foo += oo[i]*oo[i];     \
    if (foo != 0.0)					\
	for (i = 0; i < nDout; i++) oo[i] /= foo;	\
}

static Array
ApplyTransform_ConstantArray(Array array, Matrix *matrix, int norm)
{
    Type     type;
    Category cat;
    int      rank, shape[32], size, n, count, nDin, nDout, i, j;
    Array    result = NULL;
    Pointer  inVector;
    Pointer  outVector;

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

    if (norm && type != TYPE_FLOAT && type != TYPE_DOUBLE)
    {
	DXSetError(ERROR_INVALID_DATA, 
		"cannot normalize non-floating point vectors");
	return NULL;
    }
    
    if (rank != 1 || shape[0] > 3)
    {
	DXSetError(ERROR_INVALID_DATA, 
	    "data to be transformed must be vectors of dimensionality <= 3");
	return NULL;
    }

    nDin = shape[0];

    for (nDout = 3; nDout >= 1; nDout --)
    {
	int zero = (matrix->b[nDout-1] == 0.0);
	
	for (i = 0; i < nDin && zero; i++)
	    zero = (matrix->A[i][nDout-1] == 0.0);

	if (! zero)
	    break;
    }

    if (norm && type != TYPE_FLOAT && type != TYPE_DOUBLE)
    {
	DXSetError(ERROR_INVALID_DATA, 
		"cannot normalize non-floating point vectors");
	return NULL;
    }

    inVector  = DXGetConstantArrayData(array);
    outVector = DXAllocate(nDin*DXTypeSize(type));
    if (! outVector)
	return NULL;
    
    switch(type)
    {
	case TYPE_DOUBLE:
		TRANSFORM_CONSTANTARRAY(double, 0.0);
		if (norm)
		    NORMALIZE_CONSTANTARRAY(double);
		break;

	case TYPE_FLOAT:
		TRANSFORM_CONSTANTARRAY(float, 0.0);
		if (norm)
		    NORMALIZE_CONSTANTARRAY(float);
		break;

	case TYPE_INT:
		TRANSFORM_CONSTANTARRAY(int, 0.5);
		if (norm)
		    NORMALIZE_CONSTANTARRAY(int);
		break;

	case TYPE_UINT:
		TRANSFORM_CONSTANTARRAY(uint, 0.5);
		if (norm)
		    NORMALIZE_CONSTANTARRAY(uint);
		break;

	case TYPE_SHORT:
		TRANSFORM_CONSTANTARRAY(short, 0.5);
		if (norm)
		    NORMALIZE_CONSTANTARRAY(short);
		break;

	case TYPE_USHORT:
		TRANSFORM_CONSTANTARRAY(ushort, 0.5);
		if (norm)
		    NORMALIZE_CONSTANTARRAY(ushort);
		break;

	case TYPE_BYTE:
		TRANSFORM_CONSTANTARRAY(byte, 0.5);
		if (norm)
		    NORMALIZE_CONSTANTARRAY(byte);
		break;

	case TYPE_UBYTE:
		TRANSFORM_CONSTANTARRAY(ubyte, 0.5);
		if (norm)
		    NORMALIZE_CONSTANTARRAY(ubyte);
		break;
    }

    result=(Array)DXNewConstantArray(count, outVector, type, cat, rank, nDout);
    DXFree(outVector);

    return result;
}

#define TRANSFORM_GENERICARRAY(type, round)				   \
{									   \
    type *ip = (type *)inPointer;					   \
    type *op = (type *)outPointer;					   \
									   \
    for (i = 0; i < count; i++, ip += nDin, op += nDout) 		   \
    {									   \
	for (j = 0; j < nDout; j++)					   \
	{								   \
	    double f = m->b[j];						   \
									   \
	    for (k = 0; k < nDin; k++)					   \
		f += ip[k] * m->A[k][j];				   \
	    								   \
	    op[j] = (type)(f + round);					   \
	}								   \
    }									   \
}


#define NORMALIZE_GENERICARRAY(type)					   \
{									   \
    float foo;								   \
    type *op = (type *)outPointer;		   			   \
    for (j = 0; j < count; j++, op += nDout)				   \
    {									   \
	for (i = 0, foo = 0.0; i < nDout; i++) foo += op[i]*op[i];	   \
	if (foo != 0.0)							   \
	{								   \
	    foo = 1.0 / sqrt(foo);					   \
	    for (i = 0; i < nDout; i++) op[i] *= foo;			   \
	}								   \
    }									   \
}

static Array
ApplyTransform_GenericArray(Array array, Matrix *m, int norm)
{
    Category          cat;
    int               rank, shape[32], size, n;
    Array             result = NULL;
    int               i, j, k;
    int		      count;
    int		      nDin, nDout;
    Type	      type;
    Pointer  	      inPointer, outPointer;

    if (! DXGetArrayInfo(array, &count, &type, &cat, &rank, shape))
	goto error;
    
    if (rank != 1 || shape[0] > 3)
    {
	DXSetError(ERROR_INVALID_DATA, 
	    "data to be transformed must be vectors of dimensionality <= 3");
	goto error;
    }

    nDin = shape[0];

    for (nDout = 3; nDout >= 1; nDout --)
    {
	int zero = (m->b[nDout-1] == 0.0);
	
	for (i = 0; i < nDin && zero; i++)
	    zero = (m->A[i][nDout-1] == 0.0);

	if (! zero)
	    break;
    }

    result = DXNewArray(type, cat, 1, nDout);
    if (! result)
	goto error;

    if (! DXAddArrayData(result, 0, count, NULL))
	goto error;
    
    inPointer  = DXGetArrayData(array);
    outPointer = DXGetArrayData(result);

    switch(type)
    {
	case TYPE_DOUBLE:
		TRANSFORM_GENERICARRAY(double, 0.0);
		if (norm)
		    NORMALIZE_GENERICARRAY(double);
		break;

	case TYPE_FLOAT:
		TRANSFORM_GENERICARRAY(float, 0.0);
		if (norm)
		    NORMALIZE_GENERICARRAY(float);
		break;

	case TYPE_INT:
		TRANSFORM_GENERICARRAY(int, 0.5);
		if (norm)
		    NORMALIZE_GENERICARRAY(int);
		break;

	case TYPE_UINT:
		TRANSFORM_GENERICARRAY(uint, 0.5);
		if (norm)
		    NORMALIZE_GENERICARRAY(uint);
		break;

	case TYPE_SHORT:
		TRANSFORM_GENERICARRAY(short, 0.5);
		if (norm)
		    NORMALIZE_GENERICARRAY(short);
		break;

	case TYPE_USHORT:
		TRANSFORM_GENERICARRAY(ushort, 0.5);
		if (norm)
		    NORMALIZE_GENERICARRAY(ushort);
		break;

	case TYPE_BYTE:
		TRANSFORM_GENERICARRAY(byte, 0.5);
		if (norm)
		    NORMALIZE_GENERICARRAY(byte);
		break;

	case TYPE_UBYTE:
		TRANSFORM_GENERICARRAY(ubyte, 0.5);
		if (norm)
		    NORMALIZE_GENERICARRAY(ubyte);
		break;
    }

    return result;

error:
    DXDelete((Object)result);
    return NULL;
}

	    
static Array
SetDimensionality_ProductArray(Array array, int n)
{
    ProductArray result = NULL;
    int          i;
    Array        terms[100];

    if (! DXGetProductArrayInfo((ProductArray)array, &n, terms))
	goto error;
    
    for (i = 0; i < n; i++)
    {
	terms[i] = (Array)SetDimensionality((Array)terms[i], n);
	if (! terms[i])
	{
	    for (i--; i >= 0; i--)
		DXDelete((Object)terms[i]);
	    goto error;
	}
    }

    result = DXNewProductArrayV(n, terms);
    if (! result)
    {
	for (i--; i >= 0; i--)
	    DXDelete((Object)terms[i]);
	goto error;
    }

    return (Array)result;

error:
    DXDelete((Object)result);
    return NULL;
}

static Array
SetDimensionality_RegularArray(Array array, int n)
{
    Type     	   type;
    Category 	   cat;
    unsigned char  *origin = NULL;
    unsigned char  *delta = NULL;
    int      	   rank, shape[32],  count, nD, i, j;
    Array    	   result = NULL;

    if (! DXGetArrayInfo(array, &count, &type, &cat, &rank, shape))
	return NULL;
    
    if (rank != 1)
    { 
	DXSetError(ERROR_INTERNAL, "can only set dimensionality of vectors");
	return NULL;
    }
    
    if (shape[0] == n)
	return array;

    if (shape[0] > n)
	nD = shape[0];
    else
	nD = n;

    origin  = DXAllocate(nD*DXTypeSize(type));
    delta   = DXAllocate(nD*DXTypeSize(type));

    if (! origin || ! delta)
    {
	DXFree(origin);
	DXFree(delta);
	return NULL;
    }

    DXGetRegularArrayInfo((RegularArray)array, NULL, origin, delta);

    if (n > shape[0])
    {
	int inSize  = shape[0]*DXTypeSize(type)*DXCategorySize(cat);
	int outSize = n*DXTypeSize(type)*DXCategorySize(cat);

	memset(origin+inSize, 0, outSize-inSize);
	memset(delta+inSize, 0, outSize-inSize);
    }

    result  = (Array)DXNewRegularArray(type, n, count, origin, delta);

    DXFree(origin);
    DXFree(delta);

    return result;
}

static Array
SetDimensionality_ConstantArray(Array array, int n)
{
    int      rank, shape[32], count;
    Type     type;
    Category cat;
    Array    result = NULL;
    Pointer  buf = NULL;
    int	     size;

    if (! DXGetArrayInfo(array, &count, &type, &cat, &rank, shape))
	return NULL;
    
    if (rank != 1)
    { 
	DXSetError(ERROR_INTERNAL, "can only set dimensionality of vectors");
	return NULL;
    }

    if (shape[0] == n)
	return array;
    
    size = n*DXTypeSize(type)*DXCategorySize(cat);

    buf = DXAllocateZero(size);
    if (! buf)
	return NULL;
    
    if (shape[0] < n)
	size = shape[0]*DXTypeSize(type)*DXCategorySize(cat);
    
    memcpy(buf, DXGetConstantArrayData(array), size);

    result = (Array)DXNewConstantArray(count, buf, type, cat, 1, n);

    DXFree(buf);
    
    return result;
}

static Array
SetDimensionality_GenericArray(Array array, int n)
{
    Category          cat;
    Type	      type;
    int               rank, shape[32];
    Array             result = NULL;
    int		      count;
    unsigned char     *inPointer, *outPointer;
    int		      inSize, outSize, i;

    if (! DXGetArrayInfo(array, &count, &type, &cat, &rank, shape))
	goto error;
    
    if (rank != 1)
    { 
	DXSetError(ERROR_INTERNAL, "can only set dimensionality of vectors");
	return NULL;
    }

    if (shape[0] == n)
	return array;

    result = DXNewArray(type, cat, 1, n);
    if (! result)
	goto error;

    if (! DXAddArrayData(result, 0, count, NULL))
    {
	DXDelete((Object)result);
	goto error;
    }
    
    inPointer  = (unsigned char *)DXGetArrayData(array);
    outPointer = (unsigned char *)DXGetArrayData(result);

    inSize  = DXGetItemSize(array);
    outSize = DXGetItemSize(result);

    if (shape[0] < n)
    {
	int size1 = inSize;
	int size2 = outSize - inSize;

	for (i = 0; i < count; i++)
	{
	    memcpy(outPointer, inPointer, size1);
	    memset(outPointer+size1, 0, size2);
	    inPointer  += inSize;
	    outPointer += outSize;
	}
    }
    else
    {
	for (i = 0; i < count; i++)
	{
	    memcpy(outPointer, inPointer, inSize);
	    inPointer  += inSize;
	    outPointer += outSize;
	}
    }

    return result;

error:
    DXDelete((Object)result);
    return NULL;
}

static Array
SetDimensionality(Array array, int n)
{
    Array a;

    if (DXGetObjectClass((Object)array) != CLASS_ARRAY)
    {
	DXSetError(ERROR_INTERNAL, "SetDimensionality called on non-Array");
	return NULL;
    }

    switch (DXGetArrayClass(array))
    {
	case CLASS_REGULARARRAY:
	    if (DXQueryConstantArray(array, NULL, NULL))
		a = SetDimensionality_ConstantArray(array, n);
	    else
		a = SetDimensionality_RegularArray(array, n);
	    break;
	case CLASS_PRODUCTARRAY:
	    a = SetDimensionality_ProductArray(array, n);
	    break;
	case CLASS_CONSTANTARRAY:
	    a = SetDimensionality_ConstantArray(array, n);
	    break;
	case CLASS_ARRAY:
	    a = SetDimensionality_GenericArray(array, n);
	    break;
	default:
	    DXSetError(ERROR_INVALID_DATA, 
			"Inappropriate array class for transformation");
	    return NULL;
    }

    return a;
}