/***********************************************************************/
/* 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"          */
/***********************************************************************/
/*
 * $Header: /home/gda/dxcvs/dx/src/exec/dxmods/_normals.c,v 1.3 1999/05/10 15:45:20 gda Exp $
 */

#include <dxconfig.h>


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

#define TRUE 			1
#define FALSE 			0

#define ONE_THIRD		0.33333333333333333333
#define ONE_FOURTH		0.25

#define ABS(x)			(((x)< 0.0)?(-(x)):(x))

Object _dxfNormalsObject(Object, char *);

Object
_dxfNormals(Object object, char *method)
{
    Object copy = NULL;
    int  n;

    if (method == NULL)
    {
	DXSetError(ERROR_BAD_PARAMETER, "#10000", "target dependency");
	goto error;
    }

    copy = DXCopy(object, COPY_STRUCTURE);
    if (! copy)
	goto error;

    if (! _dxfNormalsObject(copy, method))
	goto error;
    
    return copy;

error:
    DXDelete(copy);
    return NULL;
}

static Error DetermineFunction(Object, char *, int *, float *);

#define NO_DEPENDENCY	   1
#define DEP_ON_POSITIONS   2
#define DEP_ON_CONNECTIONS 3
#define DEP_ON_FACES       4

static Error _dxfNormalsTraverse(Object, char *);

Object
_dxfNormalsObject(Object o, char *m)
{
    int n;
    char *gm = NULL;

    /*
     * make sure method string is in global memory
     */
    n = strlen(m)+1;
    gm = (char *)DXAllocate(n);
    if (! gm)
	goto error;
    strcpy(gm, m);
    gm[n-1] = '\0';

    if (! DXCreateTaskGroup())
	goto error;
    
    if (! _dxfNormalsTraverse(o, gm))
    {
	DXAbortTaskGroup();
	goto error;
    }

    if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE)
	goto error;
    
    DXFree((Pointer)gm);

    return o;

error:
    DXFree((Pointer)gm);
    return NULL;
}

static Error  _dxfPositionNormals(Pointer);
static Error  _dxfConnectionNormals(Pointer);
static Error  _dxfFaceNormals(Pointer);
static Error  _dxfAreaNormals(Pointer);

static int CompareDependency(Object, char *);

       Array _dxf_FLE_Normals(int *, int, int *, int, int *, int, float *, int nd);
    
typedef struct
{
    Object obj;
    float  rad;
} AreaTask;

static CompositeField GridNormalsCompositeField(CompositeField, int);
static Field GridNormalsField(Field, int);


static Error
_dxfNormalsTraverse(Object object, char *method)
{
    Class    class;
    int      i;
    Object   child;
    int      dep;
    float    rad;
    Object   comp;

    class = DXGetObjectClass(object);
    if (class == CLASS_GROUP)
	class = DXGetGroupClass((Group)object);

    switch(class)
    {
	case CLASS_FIELD:
	    if (DXEmptyField((Field)object))
		return OK;
	    
	    if (NULL != (comp = DXGetComponentValue((Field)object, "connections")))
	    {
		char *str;
		Object eType = DXGetAttribute(comp, "element type");

		if (! eType || DXGetObjectClass(eType) != CLASS_STRING)
		{
		    DXSetError(ERROR_INVALID_DATA,
			"invalid or missing element type attribute");
		    return ERROR;
		}

		str = DXGetString((String)eType);

		if (strcmp(str, "quads") && strcmp(str, "triangles"))
		    return OK;
	    }
	    else if (! DXGetComponentValue((Field)object, "faces"))
		return OK;

	    i = CompareDependency(object, method);
	    if (i == -1)
		goto error;
	    else if (i == 1)
		return OK;

	    if (! DetermineFunction(object, method, &dep, &rad))
		goto error;
	    
	    if (dep == NO_DEPENDENCY)
	    {
		DXSetError(ERROR_INVALID_DATA, "#10450", method);
		goto error;
	    }

	    if (! GridNormalsField((Field)object, dep))
	    {
		if (dep == DEP_ON_POSITIONS && rad <= 0.0)
		{
		    return DXAddTask(_dxfPositionNormals,
			(Pointer)object, 0, 1.0);
		}
		else if (dep == DEP_ON_POSITIONS && rad > 0.0)
		{
		    AreaTask task;

		    task.obj = object;
		    task.rad = rad;

		    return DXAddTask(_dxfAreaNormals,
				    (Pointer)&task, sizeof(task), 1.0);
		}
		else if (dep == DEP_ON_CONNECTIONS)
		{
		    return DXAddTask(_dxfConnectionNormals,(Pointer)object,0,1.0);
		}
		else 
		{
		    return DXAddTask(_dxfFaceNormals,(Pointer)object,0,1.0);
		}
	    }
	    
	    break;
	
	case CLASS_COMPOSITEFIELD:
	    
	    i = CompareDependency(object, method);
	    if (i == -1)
		goto error;
	    else if (i == 1)
		return OK;

	    if (! DetermineFunction(object, method, &dep, &rad))
		goto error;
	    
	    if (GridNormalsCompositeField((CompositeField)object, dep))
		return OK;
	    
	    if (dep == NO_DEPENDENCY)
	    {
		DXSetError(ERROR_INVALID_DATA, "#10450", "method");
		goto error;
	    }

	    /*
	     * positions and area require  special stuff (eg. Grow) at
	     * the CompositeField level
	     */

	    if (dep == DEP_ON_POSITIONS && rad <= 0.0)
	    {
		return DXAddTask(_dxfPositionNormals,(Pointer)object,0,1.0);
	    }
	    else if (dep == DEP_ON_POSITIONS && rad > 0.0)
	    {
		AreaTask task;

		task.obj = object;
		task.rad = rad;

		return DXAddTask(_dxfAreaNormals,
				(Pointer)&task, sizeof(task), 1.0);
	    }

	    /*
	     * else connections or faces, and since we don't do anything
	     * special at the composite field level, so fall through to
	     * the generic group case.
	     */
	
	case CLASS_MULTIGRID:
	case CLASS_SERIES:
	case CLASS_GROUP:

	    i = 0;
	    for (;;)
	    {
	    	child = DXGetEnumeratedMember((Group)object, i++, NULL);
		if (! child)
		    break;

		if (! _dxfNormalsTraverse(child, method))
		    goto error;
	    }
	    
	    break;

	case CLASS_XFORM:

	    if (! DXGetXformInfo((Xform)object, &child, 0))
		goto error;
	    
	    if (! _dxfNormalsObject(child, method))
		goto error;
	    
	    break;

	case CLASS_CLIPPED:

	    if (! DXGetClippedInfo((Clipped)object, &child, 0))
		goto error;
	    
	    if (! _dxfNormalsObject(child, method))
		goto error;
	    
	    break;
	
	default:
	    break;
    }

    return OK;

error:
    return ERROR;
}

static Error _DetermineFunction(Object, char *, int *);

static Error
DetermineFunction(Object o, char *method, int *dep, float *radius)
{
    int i;

    if (method == NULL || !strcmp(method, "manhattan"))
    {
	*radius = -1.0;
	*dep = DEP_ON_POSITIONS;
	return OK;
    }

    for (i=0; method[i]; i++)
	if ((method[i] < '0' || method[i] > '9') && method[i] != '.')
	    break;
    
    if (method[i])
    {
	*radius = 0;
	return _DetermineFunction(o, method, dep);
    }
    else
    {
	sscanf(method, "%f", radius);
	*dep = DEP_ON_POSITIONS;
	return OK;
    }
}

static Error
_DetermineFunction(Object o, char *name, int *dep)
{
    Class class = DXGetObjectClass(o);
    if (class == CLASS_GROUP)
	class = DXGetGroupClass((Group)o);

    switch(class)
    {
	case CLASS_FIELD:
	{
	    Object attr; char *str; Array a;

	    if (DXEmptyField((Field)o))
	    { 
		*dep = NO_DEPENDENCY;
		return OK;
	    }

	    /*
	     * make sure the component is there
	     */
	    a = (Array)DXGetComponentValue((Field)o, name);
	    if (a == NULL)
	    {
		DXSetError(ERROR_MISSING_DATA, "#10240", name);
		*dep = NO_DEPENDENCY;
		return ERROR;
	    }

	    if (! strcmp(name, "connections"))
	    {
		*dep = DEP_ON_CONNECTIONS;
		return OK;
	    }

	    if (! strcmp(name, "faces"))
	    {
		*dep = DEP_ON_FACES;
		return OK;
	    }

	    if (! DXGetStringAttribute((Object)a, "dep", &str))
	    {
		DXSetError(ERROR_INVALID_DATA, "#10255", name, "dep");
		*dep = NO_DEPENDENCY;
		return ERROR;
	    }

	    if (! strcmp(str, "positions"))
	    {
		*dep = DEP_ON_POSITIONS;
		return OK;
	    }
	    else if (! strcmp(str, "connections"))
	    {
		*dep = DEP_ON_CONNECTIONS;
		return OK;
	    }
	    else if (! strcmp(str, "faces"))
	    {
		*dep = DEP_ON_CONNECTIONS;
		return OK;
	    }
	    else
	    {
		DXSetError(ERROR_INVALID_DATA, "#11250", name);
		*dep = NO_DEPENDENCY;
		return ERROR;
	    }
	}

	case CLASS_COMPOSITEFIELD:
	{
	    Object c;
	    int i;

	    i = 0;
	    for (*dep = NO_DEPENDENCY; *dep == NO_DEPENDENCY; )
	    {
		c = DXGetEnumeratedMember((Group)o, i++, NULL);
		if (c == NULL)
		    return OK;
		
		if (! _DetermineFunction(c, name, dep))
		    return ERROR;
	    }

	    return OK;
	}

	default:
	{
	    DXSetError(ERROR_INVALID_DATA, "#11381");
	    return ERROR;
	}
    }
}

static Error _dxfConnectionNormals_Field(Field);

static Error
connectionsNormalsTask(Pointer ptr)
{
    Field field = (Field)ptr;

    return _dxfConnectionNormals_Field(field);
}

static Error
_dxfConnectionNormals(Pointer ptr)
{
     return _dxfConnectionNormals_Field((Field)ptr);
}


static Array  _dxf_CDep_triangles(Field);
static Array  _dxf_CDep_quads(Field);

#define NORMAL(nv, nd, norm)						    \
{									    \
    float length;							    \
    int i, j;								    \
									    \
    norm->x = 0.0;							    \
    norm->y = 0.0;							    \
    norm->z = 0.0;							    \
									    \
    if (nd == 2)							    \
	for (i = 0; i < nv; i++)					    \
	{								    \
	    int j = (i == nv-1) ? 0 : i + 1;				    \
	    norm->z += (ptrs[i][0] - ptrs[j][0]) *			    \
				    (ptrs[i][1] + ptrs[j][1]);		    \
	}								    \
    else								    \
	for (i = 0; i < nv; i++)					    \
	{								    \
	    int j = (i == nv-1) ? 0 : i + 1;				    \
	    norm->x += (ptrs[i][1] - ptrs[j][1]) *			    \
				    (ptrs[i][2] + ptrs[j][2]);		    \
	    norm->y += (ptrs[i][2] - ptrs[j][2]) *			    \
				    (ptrs[i][0] + ptrs[j][0]);		    \
	    norm->z += (ptrs[i][0] - ptrs[j][0]) *			    \
				    (ptrs[i][1] + ptrs[j][1]);		    \
	}								    \
									    \
    length = vector_length(norm);					    \
    if (length == 0.0)							    \
    {									    \
	norm->x = 0.0;							    \
	norm->y = 0.0;							    \
	norm->z = 1.0;							    \
	goto degenerate;						    \
    }									    \
    else								    \
    {									    \
	float d = 1.0 / length;						    \
	vector_scale(norm, d, norm);					    \
    }									    \
}

static Error
_dxfConnectionNormals_Field(Field field)
{
    Array		array;
    int			rank, shape[10];
    Type		type;
    Category		cat;
    int			n;
    char		*primitive;
    Object 		attr;

    /*
     * Check first whether the input field is empty.
     */

    if (DXEmptyField(field))
	return OK;

    array = (Array)DXGetComponentValue(field, "connections");
    if (! array)
	return OK;
    
    DXGetArrayInfo(array, &n, NULL, NULL, NULL, NULL);
    
    DXDeleteComponent(field, "normals");
    
    array = (Array)DXGetComponentValue(field, "positions");
    DXGetArrayInfo(array, NULL, &type, &cat, &rank, shape);

    if (rank != 1 || shape[0] > 3 || shape[0] < 2)
    {
	DXSetError(ERROR_INVALID_DATA, "positions must be 2 or 3-D");
	return ERROR;
    }
    else
    {
	attr = DXGetComponentAttribute(field, "connections", "element type");
	if (!attr)
	{
	    DXSetError(ERROR_INVALID_DATA, "#10255", 
				"connections", "element type");
	    return ERROR;
	}

	if (DXGetObjectClass(attr) != CLASS_STRING)
	{
	    DXSetError(ERROR_INVALID_DATA, "#10200", "element type attribute");
	    return ERROR;
	}

	/*
	 * Check the topological primitive.
	 */

	primitive = DXGetString((String)attr);

	if (! strcmp(primitive, "triangles"))
	    array = _dxf_CDep_triangles(field);
	else if (! strcmp(primitive, "quads"))
	    array = _dxf_CDep_quads(field);
	else
	    return OK;
    }

    if (! array)
	return ERROR;

    /*
     * Set the components of the output fields.
     */
    if (! DXSetComponentValue(field, "normals",(Object) array))
	return ERROR;

    if (! DXSetComponentAttribute(field, "normals", "dep",
				(Object) DXNewString("connections")))
        return ERROR;

    if (! DXEndField(field))
	return ERROR;

    return OK;
}

static Array
_dxf_CDep_triangles(Field field)
{
    Array		a_connections, a_positions;
    Array		a_normals = NULL;
    int			n_connections;
    Vector		*normals;
    int			i, j, k;
    int			degenerate_count;
    ArrayHandle		cHandle = NULL, pHandle = NULL;
    int			nd;

    /*
     * Do quick type checks to ensure that the components are correct.
     */
    a_connections = (Array)DXGetComponentValue(field, "connections");
    if (!a_connections)
        goto error;

    a_positions = (Array)DXGetComponentValue(field, "positions");
    if (!a_positions)
    {
	DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
	goto error;
    }
    
    if (! DXGetArrayInfo(a_positions, &n_connections, NULL, NULL, NULL, &nd))
	goto error;

    if (! DXGetArrayInfo(a_connections, &n_connections, NULL, NULL, NULL, NULL))
	goto error;

    a_normals = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
    if (!a_normals)
	goto error;

    /*
     * We now assume that there is at least one connection.
     */

    if (! DXTypeCheck(a_connections, TYPE_INT, CATEGORY_REAL, 1, 3))
	goto error;

    /*
     * Setup the final array information.
     */

    if (! DXAddArrayData(a_normals, 0, n_connections,(Pointer) NULL))
        goto error;

    /*
     * Setup and initialize the components.
     */

    if (NULL ==(normals =(Vector *)DXGetArrayData(a_normals)))
	goto error;

    if (NULL ==(cHandle = DXCreateArrayHandle(a_connections)))
	goto error;
    
    if (NULL ==(pHandle = DXCreateArrayHandle(a_positions)))
	goto error;


    /*
     * For each connection...
     */
    degenerate_count = 0;
    for(i = 0; i < n_connections; i++, normals++)
    {
	int     *element, ebuf[3];
	float	*ptrs[3], buf[9];

	element =(int *)DXCalculateArrayEntry(cHandle, i,(Pointer)ebuf);

	for(j = 0; j < 3; j++)						   
	    ptrs[j] = (float *)DXGetArrayEntry(pHandle, element[j], buf+j*nd); 

	NORMAL(3, nd, normals);

	continue;

degenerate:
	degenerate_count ++;
    }

    /*
     * Set the components of the output fields.
     */
    if (! DXSetComponentValue(field, "normals",(Object) a_normals))
	goto error;

    if (! DXSetComponentAttribute(field, "normals", "dep",
				(Object) DXNewString("connections")))
        goto error;

    if (! DXEndField(field))
	goto error;

    if (degenerate_count > 0)
	DXMessage("Number of degenerate triangles: %d", degenerate_count);

    DXFreeArrayHandle(cHandle);
    DXFreeArrayHandle(pHandle);

    return a_normals;

error:
    if (cHandle)
	DXFreeArrayHandle(cHandle);

    if (pHandle)
	DXFreeArrayHandle(pHandle);

    if (a_normals)
	DXDelete((Object)a_normals);

    return NULL;
	
}

static Array
_dxf_CDep_quads(Field field)
{
    Array		a_connections, a_positions;
    Array		a_normals = NULL;
    int			n_connections;
    Vector		*normals;
    int			i, j, nd;
    int			degenerate_count;
    ArrayHandle		cHandle = NULL, pHandle = NULL;

    /*
     * Do quick type checks to ensure that the components are correct.
     */
    a_connections = (Array)DXGetComponentValue(field, "connections");
    if (!a_connections)
        goto error;

    a_positions = (Array)DXGetComponentValue(field, "positions");
    if (!a_positions)
    {
	DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
	goto error;
    }
    
    if (! DXGetArrayInfo(a_positions, NULL, NULL, NULL, NULL, &nd))
	goto error;

    if (! DXGetArrayInfo(a_connections, &n_connections, NULL, NULL, NULL, NULL))
	goto error;

    a_normals = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
    if (!a_normals)
	goto error;

    /*
     * We now assume that there is at least one connection.
     */

    if (! DXTypeCheck(a_connections, TYPE_INT, CATEGORY_REAL, 1, 4))
	goto error;

    /*
     * Setup the final array information.
     */

    if (! DXAddArrayData(a_normals, 0, n_connections,(Pointer) NULL))
        goto error;

    /*
     * Setup and initialize the components.
     */

    if (NULL ==(normals =(Vector *)DXGetArrayData(a_normals)))
	goto error;

    if (NULL ==(cHandle = DXCreateArrayHandle(a_connections)))
	goto error;
    
    if (NULL ==(pHandle = DXCreateArrayHandle(a_positions)))
	goto error;

    /*
     * For each connection...
     */
    degenerate_count = 0;
    for(i = 0; i < n_connections; i++, normals++)
    {
	int *element, ebuf[4];
	float *ptrs[4], buf[12];

	element =(int *)DXCalculateArrayEntry(cHandle, i,(Pointer)ebuf);

	ptrs[0] = (float *)DXGetArrayEntry(pHandle, element[0], buf+0*nd); 
	ptrs[1] = (float *)DXGetArrayEntry(pHandle, element[2], buf+1*nd); 
	ptrs[2] = (float *)DXGetArrayEntry(pHandle, element[3], buf+2*nd); 
	ptrs[3] = (float *)DXGetArrayEntry(pHandle, element[1], buf+3*nd); 

	NORMAL(4, nd, normals);

	continue;

degenerate:
	degenerate_count ++;
    }

    if (degenerate_count > 0)
	DXMessage("Number of degenerate triangles: %d", degenerate_count);

    DXFreeArrayHandle(cHandle);
    DXFreeArrayHandle(pHandle);

    return a_normals;

error:
    if (cHandle)
	DXFreeArrayHandle(cHandle);

    if (pHandle)
	DXFreeArrayHandle(pHandle);

    if (a_normals)
	DXDelete((Object)a_normals);

    return NULL;
	
}

static Error _dxfPositionNormals_Field(Field);

static Error
positionsNormalsTask(Pointer ptr)
{
    Field field = (Field)ptr;

    return _dxfPositionNormals_Field(field);
}

static Error
_dxfPositionNormals(Pointer ptr)
{
    Object object = (Object)ptr;
    Object child;
    int i;
    Class class = DXGetObjectClass(object);

    if (class == CLASS_GROUP)
    {
	if (! DXGrow(object, 1, NULL, NULL))
	    return ERROR;

	if (! DXCreateTaskGroup())
	    return ERROR;

	for (i = 0; ;i++)
	{
	    child = DXGetEnumeratedMember((Group)object, i, NULL);
	    if (! child)
		break;
	    
	    if (! DXAddTask(positionsNormalsTask, (Pointer)child, 0, 1.0))
	    {
		DXAbortTaskGroup();
		return ERROR;
	    }
	}

	if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE)
	    return ERROR;

	if (! DXShrink(object))
	    return ERROR;
	
	return OK;
    }
    else 
	 return _dxfPositionNormals_Field((Field)object);
}

static Array  _dxf_PDep_triangles(Field);
static Array  _dxf_PDep_quads(Field);

static Error
_dxfPositionNormals_Field(Field field)
{
    Array		array;
    int			rank, shape[10];
    Type		type;
    Category		cat;
    int			n;
    char		*primitive;
    Object 		attr;

    /*
     * Check first whether the input field is empty.
     */

    if (DXEmptyField(field))
	return OK;

    array = (Array)DXGetComponentValue(field, "connections");
    if (! array)
	return OK;
    
    DXGetArrayInfo(array, &n, NULL, NULL, NULL, NULL);
    if (n <= 0)
	return OK;
    
    DXDeleteComponent(field, "normals");
    
    array = (Array)DXGetComponentValue(field, "positions");
    if (! array)
    {
	DXSetError(ERROR_INVALID_DATA, "#10200", "positions component");
	return ERROR;
    }

    DXGetArrayInfo(array, &n, &type, &cat, &rank, shape);

    if (rank != 1 || shape[0] < 2 || shape[0] > 3)
    {
	DXSetError(ERROR_INVALID_DATA, "positions must be 2 or 3-D");
	return ERROR;
    }
    else
    {
	attr = DXGetComponentAttribute(field, "connections", "element type");
	if (!attr)
	{
	    DXSetError(ERROR_INVALID_DATA, "#10255",
					"connections", "element type");
	    return ERROR;
	}

	if (DXGetObjectClass(attr) != CLASS_STRING)
	{
	    DXSetError(ERROR_INVALID_DATA, "#10200", "element type attribute");
	    return ERROR;
	}

	/*
	 * Check the topological primitive.
	 */

	primitive = DXGetString((String)attr);

	if (! strcmp(primitive, "triangles"))
	    array =_dxf_PDep_triangles(field);
	else if (! strcmp(primitive, "quads"))
	    array =_dxf_PDep_quads(field);
	else
	    return OK;
    }

    if (! array)
	return ERROR;

    /*
     * Set the components of the output fields.
     */

    if (DXGetComponentValue(field, "original normals"))
	DXDeleteComponent(field, "original normals");

    if (! (DXSetComponentValue(field, "normals",(Object)array) &&
	   DXChangedComponentValues(field, "normals")	       &&
	   DXEndField(field)))
	return ERROR;

    return OK;
}

static Array
_dxf_PDep_triangles(Field field)
{
    Array		nArray;
    Array		a_connections, a_positions;
    int			n_connections, n_positions;
    float		*pts;
    int		        *tri;
    float		centroid[3];
    float		length;
    float		dist;
    float		invdist;
    Vector		*normals, *n;
    int			i, nd;
    int			degenerate_count;
    ArrayHandle		pHandle = NULL;
    ArrayHandle		cHandle = NULL;

    degenerate_count = 0;

    /* 
     * Initialize outputs.
     */

    nArray = NULL;
    if (NULL == (nArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3)))
	goto error;

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

    DXGetArrayInfo(a_connections, &n_connections, NULL, NULL, NULL, NULL);
    DXGetArrayInfo(a_positions, &n_positions, NULL, NULL, NULL, &nd);

    if (! DXTypeCheck(a_connections, TYPE_INT, CATEGORY_REAL, 1, 3))
	goto error;

    /*
     * Setup the final array information.
     */
    if (! DXAddArrayData(nArray, 0, n_positions, (Pointer)NULL))
        goto error;

    /*
     * Setup and initialize the components.
     */

    normals = (Vector *)DXGetArrayData(nArray);

    for(i = 0, n = normals; i < n_positions; i++, n++)
	n->x = n->y = n->z = 0.0;

    /*
     * For each connection...
     */

    pHandle = DXCreateArrayHandle(a_positions);
    cHandle = DXCreateArrayHandle(a_connections);
    if (! pHandle || ! cHandle)
	goto error;

    for(i = 0; i < n_connections; i++)
    {
	Vector enorm;
	int    j, *element, ebuf[3];
	float  *ptrs[3], buf[9];

	element =(int *)DXCalculateArrayEntry(cHandle, i,(Pointer)ebuf);

	for(j = 0; j < 3; j++)						   
	    ptrs[j] = (float *)DXGetArrayEntry(pHandle, element[j], buf+j*nd); 

	NORMAL(3, nd, (&enorm));

        /*
	 * Determine the coordinates of the centroid of the triangle.
	 */

	if (nd == 3)
	{
	    centroid[0] = (ptrs[0][0]+ptrs[1][0]+ptrs[2][0])*ONE_THIRD;
	    centroid[1] = (ptrs[0][1]+ptrs[1][1]+ptrs[2][1])*ONE_THIRD;
	    centroid[2] = (ptrs[0][2]+ptrs[1][2]+ptrs[2][2])*ONE_THIRD;

	    for (j = 0; j < 3; j++)
	    {
		float x, y, z;
		x = ptrs[j][0] - centroid[0];
		y = ptrs[j][1] - centroid[1];
		z = ptrs[j][2] - centroid[2];
		dist = ABS(x)+ ABS(y)+ ABS(z);
		if (dist > 0.0)
		{
		    invdist = 1.0 / dist;
		    n = normals + element[j];
		    n->x += enorm.x * invdist;
		    n->y += enorm.y * invdist;
		    n->z += enorm.z * invdist;
		}
	    }
	}
	else
	{
	    centroid[0] = (ptrs[0][0]+ptrs[1][0]+ptrs[2][0])*ONE_THIRD;
	    centroid[1] = (ptrs[0][1]+ptrs[1][1]+ptrs[2][1])*ONE_THIRD;

	    for (j = 0; j < 3; j++)
	    {
		float x, y, z;
		x = ptrs[j][0] - centroid[0];
		y = ptrs[j][1] - centroid[1];
		dist = ABS(x)+ ABS(y);
		if (dist > 0.0)
		{
		    invdist = 1.0 / dist;
		    n = normals + element[j];
		    n->z += enorm.z * invdist;
		}
	    }
	}

	continue;

degenerate:
	degenerate_count ++;						    
    }

    DXFreeArrayHandle(pHandle);
    DXFreeArrayHandle(cHandle);

    /*
     * Now normalize all the normals
     */

    for(i = 0; i < n_positions; i++)
    {
	if (vector_zero((Vector *)normals))
	{
	    normals->x = (float)0.0;
	    normals->y = (float)0.0;
	    normals->z = (float)1.0;
	}

	if (! vector_normalize((Vector *)normals, (Vector *)normals))
	    goto error;
	
	normals ++;
    }

    if (degenerate_count > 0)
	DXMessage("Number of degenerate triangles: %d", degenerate_count);

    return nArray;

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

static Array
_dxf_PDep_quads(Field field)
{
    Array		nArray;
    Array		a_connections, a_positions;
    int			n_connections, n_positions;
    float		*pts;
    int		        *quad;
    float		centroid[3];
    float		length;
    float		dist;
    float		invdist;
    Vector		*normals, *n;
    int			i, nd;
    int			degenerate_count;
    ArrayHandle		pHandle = NULL;
    ArrayHandle		cHandle = NULL;

    degenerate_count = 0;

    /* 
     * Initialize outputs.
     */

    nArray = NULL;

    if (NULL == (nArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3)))
	goto error;

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

    DXGetArrayInfo(a_connections, &n_connections, NULL, NULL, NULL, NULL);
    DXGetArrayInfo(a_positions, &n_positions, NULL, NULL, NULL, &nd);

    if (! DXTypeCheck(a_connections, TYPE_INT, CATEGORY_REAL, 1, 4))
	goto error;

    /*
     * Setup the final array information.
     */
    if (! DXAddArrayData(nArray, 0, n_positions, (Pointer)NULL))
        goto error;

    /*
     * Setup and initialize the components.
     */

    normals = (Vector *)DXGetArrayData(nArray);

    for(i = 0, n = normals; i < n_positions; i++, n++)
	n->x = n->y = n->z = 0.0;

    pHandle = DXCreateArrayHandle(a_positions);
    cHandle = DXCreateArrayHandle(a_connections);
    if (! pHandle || ! cHandle)
	goto error;

    for(i = 0; i < n_connections; i++)
    {
	Vector enorm;
	int    j, *element, ebuf[4];
	float  *ptrs[4], buf[12];

	element = (int *)DXCalculateArrayEntry(cHandle, i,(Pointer)ebuf);

	ptrs[0] = (float *)DXGetArrayEntry(pHandle, element[0], buf+0*nd); 
	ptrs[1] = (float *)DXGetArrayEntry(pHandle, element[2], buf+1*nd); 
	ptrs[2] = (float *)DXGetArrayEntry(pHandle, element[3], buf+2*nd); 
	ptrs[3] = (float *)DXGetArrayEntry(pHandle, element[1], buf+3*nd); 

	NORMAL(4, nd, (&enorm));

	/*
	 * Determine the addresses of the normals at the 
	 * three points.
	 */
	if (nd == 3)
	{
	    centroid[0] = (ptrs[0][0]+ptrs[1][0]+ptrs[2][0]+ptrs[3][0])*ONE_FOURTH;
	    centroid[1] = (ptrs[0][1]+ptrs[1][1]+ptrs[2][1]+ptrs[3][1])*ONE_FOURTH;
	    centroid[2] = (ptrs[0][2]+ptrs[1][2]+ptrs[2][2]+ptrs[3][2])*ONE_FOURTH;

	    for (j = 0; j < 4; j++)
	    {
		float x, y, z;
		x = ptrs[j][0] - centroid[0];
		y = ptrs[j][1] - centroid[1];
		z = ptrs[j][2] - centroid[2];
		dist = ABS(x)+ ABS(y)+ ABS(z);
		if (dist > 0.0)
		{
		    invdist = 1.0 / dist;
		    n = normals + element[j];
		    n->x += enorm.x * invdist;
		    n->y += enorm.y * invdist;
		    n->z += enorm.z * invdist;
		}
	    }
	}
	else
	{
	    centroid[0] = (ptrs[0][0]+ptrs[1][0]+ptrs[2][0]+ptrs[3][0])*ONE_FOURTH;
	    centroid[1] = (ptrs[0][1]+ptrs[1][1]+ptrs[2][1]+ptrs[3][1])*ONE_FOURTH;

	    n->x = 0.0;
	    n->y = 0.0;

	    for (j = 0; j < 4; j++)
	    {
		float x, y, z;
		x = ptrs[j][0] - centroid[0];
		y = ptrs[j][1] - centroid[1];
		dist = ABS(x)+ ABS(y);
		if (dist > 0.0)
		{
		    invdist = 1.0 / dist;
		    n = normals + element[j];
		    n->z += enorm.z * invdist;
		}
	    }
	}

	continue;

degenerate:
	degenerate_count ++;						   
    }

    DXFreeArrayHandle(pHandle);
    DXFreeArrayHandle(cHandle);

    /*
     * Now normalize all the normals
     */

    for(i = 0; i < n_positions; i++)
    {
	if (vector_zero((Vector *)normals))
	{
	    normals->x = (float)0.0;
	    normals->y = (float)0.0;
	    normals->z = (float)1.0;
	}

	if (! vector_normalize((Vector *)normals, (Vector *)normals))
	    goto error;
	
	normals ++;
    }

    if (degenerate_count > 0)
	DXMessage("Number of degenerate triangles: %d", degenerate_count);

    return nArray;

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

typedef struct
{
    Field field;
    float radius;
} AreaArgs;

static Error _dxfAreaNormals_Field(Field, float);

static Error
areaNormalsTask(Pointer ptr)
{
    AreaArgs *args = (AreaArgs *)ptr;

    return _dxfAreaNormals_Field(args->field, args->radius);
}

static Error
_dxfAreaNormals(Pointer ptr)
{
    Object object = ((AreaTask *)ptr)->obj;
    float  radius = ((AreaTask *)ptr)->rad;
    int i;
    Object child;
    
    if (DXGetObjectClass(object) == CLASS_COMPOSITEFIELD)
    {
	if (! DXGrow(object, 1, NULL, NULL))
	    return ERROR;

	if (! DXCreateTaskGroup())
	    return ERROR;

	for (i = 0; ;i++)
	{
	    AreaArgs a;

	    child = DXGetEnumeratedMember((Group)object, i, NULL);
	    if (! child)
		break;

	    a.radius = radius;
	    a.field  = (Field)child;
	    
	    if (! DXAddTask(areaNormalsTask, (Pointer)&a, sizeof(a), 1.0))
	    {
		DXAbortTaskGroup();
		return ERROR;
	    }
	}

	if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE)
	    return ERROR;

	if (! DXShrink(object))
	    return ERROR;
	
	return OK;
    }
    else 
	 return _dxfAreaNormals_Field((Field)object, radius);
}

static float  _dxfTriangleArea(float *, float *, float *);
static int    _dxfInside(float *, float, float *);
static float  _dxfDistance(float *, float *);
static Error  _dxfCheckTriangle(int, float, int, int *,
		    int *, int *, float *, float *, float *);
static float  _dxfIntersect(float *, float, float *, float *);
static float  _dxfTriangleArea(float *, float *, float *);

static float scale = 1.0;

static Error
_dxfAreaNormals_Field(Field field, float radius)
{
    Array	triArray, nbrArray, ptArray, nrmArray;
    int		nTriangles;
    int		nPoints;
    float	*pts, *normals, *facets;
    int		*tris, *nbrs;
    int		*pointTris;
    int 	*pointFlags;
    float	*f;
    int		*t;
    int		p, q, r, i, j;
    float	v0[3], v1[3];
    float	totEdgeLen;

    facets = NULL;
    nrmArray = NULL;

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

    DXMarkTimeLocal("nbrs start");

    nbrArray = DXNeighbors(field);
    if (! nbrArray)
	goto error;
    
    nbrs = (int *)DXGetArrayData(nbrArray);
    pts = (float *)DXGetArrayData(ptArray);
    
    if (! DXQueryOriginalSizes(field, &nPoints, NULL))
	DXGetArrayInfo(ptArray, &nPoints, NULL, NULL, NULL, NULL);

    DXGetArrayInfo(triArray, &nTriangles, NULL, NULL, NULL, NULL);

    pointTris = (int *)DXAllocate(nPoints * sizeof(int));
    if (! pointTris)
	goto error;
    
    memset((int *)pointTris, -1, nPoints * sizeof(int));

    tris = (int *)DXGetArrayData(triArray);

    facets = (float *)DXAllocate(nTriangles*3*sizeof(float));
    if (! facets)
	goto error;

    memset((int *)facets, 0, nTriangles*3*sizeof(float));

    f = facets;
    t = tris;
    totEdgeLen = 0.0;
    for(i = 0; i < nTriangles; i++)
    {
	p = t[0];
	q = t[1];
	r = t[2];

	totEdgeLen += _dxfDistance(pts+(3*p), pts+(3*q));
	totEdgeLen += _dxfDistance(pts+(3*q), pts+(3*r));
	totEdgeLen += _dxfDistance(pts+(3*r), pts+(3*p));

	v0[0] = pts[3*p + 0] - pts[3*q + 0];
	v0[1] = pts[3*p + 1] - pts[3*q + 1];
	v0[2] = pts[3*p + 2] - pts[3*q + 2];

	v1[0] = pts[3*p + 0] - pts[3*r + 0];
	v1[1] = pts[3*p + 1] - pts[3*r + 1];
	v1[2] = pts[3*p + 2] - pts[3*r + 2];

	f[0] = v0[1] * v1[2] - v1[1] * v0[2];
	f[1] = v0[2] * v1[0] - v1[2] * v0[0];
	f[2] = v0[0] * v1[1] - v1[0] * v0[1];

	vector_normalize((Vector *)f, (Vector *)f);

	f += 3;

	for (j = 0; j < 3; j++)
	{
	    if (*t < nPoints)
		pointTris[*t] = i;
	    t++;
	}
    }

    radius *= (totEdgeLen /(3*nTriangles));

    nrmArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
    if (! nrmArray)
	goto error;

    if (! DXAddArrayData(nrmArray, 0, nPoints, NULL))
	goto error;

    normals = (float *)DXGetArrayData(nrmArray);

    memset((int *)normals, 0, nPoints*3*sizeof(float));

    pointFlags = (int *)DXAllocate(nTriangles * sizeof(int));
    if (! pointFlags)
	goto error;
    
    memset((int *)pointFlags, -1, nTriangles * sizeof(int));

    for(i = 0; i < nPoints; i++)
    {
	if (pointTris[i] >= 0)
	{
	    _dxfCheckTriangle(i, radius, pointTris[i], pointFlags,
		    tris, nbrs, pts, facets, normals);
	    if (! vector_normalize((Vector *)(normals+3*i), 
					(Vector *)(normals+3*i)))
	    {
		normals[3*i+0] = 0.0;
		normals[3*i+1] = 1.0;
		normals[3*i+2] = 0.0;
	    }
	}
	else
	{
	    normals[3*i+0] = 0.0;
	    normals[3*i+1] = 1.0;
	    normals[3*i+2] = 0.0;
	}
    }

    DXFree((Pointer)facets);
    DXFree((Pointer)pointTris);
    DXFree((Pointer)pointFlags);

    /*
     * Set the components of the output fields.
     */

    if (DXGetComponentValue(field, "original normals"))
	DXDeleteComponent(field, "original normals");

    if (! (DXSetComponentValue(field, "normals",(Object)nrmArray) &&
	   DXSetComponentAttribute(field, "normals", "dep", 	
					(Object)DXNewString("positions")) &&
	   DXChangedComponentValues(field, "normals")		        &&
	   DXEndField(field)))
    {
	goto error;
    }

    return OK;

error:
    DXFree((Pointer)facets);
    DXFree((Pointer)pointTris);
    DXFree((Pointer)pointFlags);
    DXDelete((Object)nrmArray);

    return ERROR;
}

static Error
_dxfCheckTriangle(int centerIndex, float radius, int triIndex, int *flags,
		    int *triangles, int *neighbors, float *positions, 
		    float *facets, float *normals)
{
    float *p, *q, *r;
    int   *triangle;
    int	  *nbrs;
    float *normal;
    int   in[3];
    float trap[12];
    float *center;
    float *facet;
    int   i, j, next;
    int   allOut, allIn;
    float area, t;

    if (flags[triIndex] == centerIndex)
	return OK;
    
    flags[triIndex] = centerIndex;

    center   = positions +(3*centerIndex);
    normal   = normals   +(3*centerIndex);
    facet    = facets    +(3*triIndex);
    triangle = triangles +(3*triIndex);
    nbrs     = neighbors +(3*triIndex);

    allIn  = 1;
    allOut = 1;
    for(i = 0; i < 3; i++)
    {
	p = positions +(3*triangle[i]);
	in[i] = _dxfInside(center, radius, p);

	if (in[i])
	    allOut = 0;
	else
	    allIn = 0;
    }

    if (allOut)
	return OK;
    
    if (allIn)
    {
	p = positions + 3*triangle[0];
	q = positions + 3*triangle[1];
	r = positions + 3*triangle[2];

	area  = _dxfTriangleArea(p, q, r);
	facet = facets + 3*(triIndex);

	normal[0] += area*facet[0];
	normal[1] += area*facet[1];
	normal[2] += area*facet[2];
    }
    else
    {
	j = 0;
	for(i = 0; i < 3; i++)
	{
	    if (i == 2)
		next = 0;
	    else
		next = i + 1;
	    
	    p = positions +(3*triangle[i]);
	    q = positions +(3*triangle[next]);

	    if (in[i])
	    {
		trap[j++] = p[0];
		trap[j++] = p[1];
		trap[j++] = p[2];

		if (! in[next])
		{
		    t = _dxfIntersect(center, radius, p, q);
		    if (t == -1.0)
		    {
			DXMessage("Intersection error");
			t = 0.0;
		    }

		    trap[j++] = p[0] + t*(q[0]-p[0]);
		    trap[j++] = p[1] + t*(q[1]-p[1]);
		    trap[j++] = p[2] + t*(q[2]-p[2]);
		}
	    }
	    else
	    {
		if (in[next])
		{
		    t = _dxfIntersect(center, radius, p, q);
		    if (-0.001 > t || 1.0001 < t)
		    {
			DXMessage("Intersection error: %f", t);
			if (t < 0.5)t = 0.0;
			else t = 1.0;
		    }

		    trap[j++] = p[0] + t*(q[0]-p[0]);
		    trap[j++] = p[1] + t*(q[1]-p[1]);
		    trap[j++] = p[2] + t*(q[2]-p[2]);
		}
	    }
	}

	if (j > 12)
	    DXMessage("clip error");
		
	facet = facets + 3*(triIndex);

	if (j == 12)
	{
	    area  = _dxfTriangleArea(trap+0, trap+3, trap+6);
	    normal[0] += area*facet[0];
	    normal[1] += area*facet[1];
	    normal[2] += area*facet[2];

	    area  = _dxfTriangleArea(trap+0, trap+6, trap+9);
	    normal[0] += area*facet[0];
	    normal[1] += area*facet[1];
	    normal[2] += area*facet[2];
	}
	else
	{
	    area  = _dxfTriangleArea(trap+0, trap+3, trap+6);
	    normal[0] += area*facet[0];
	    normal[1] += area*facet[1];
	    normal[2] += area*facet[2];
	}
    }

    for(i = 0; i < 3; i++)
    {
	if (i == 2)
	    next = 0;
	else
	    next = i + 1;
	
	if (nbrs[i] >= 0 &&(in[i] || in[next]))
	    _dxfCheckTriangle(centerIndex, radius, nbrs[i], flags,
		    triangles, neighbors, positions, facets, normals);
    }

    return OK;
}

static int
_dxfInside(float *center, float radius, float *pt)
{
    return(_dxfDistance(center, pt)< radius);
}

static float
_dxfDistance(float *p, float *q)
{
    float a, b, c;

    a = *p++ - *q++;
    b = *p++ - *q++;
    c = *p++ - *q++;

    return sqrt(a*a + b*b + c*c);
}

static float
_dxfTriangleArea(float *p, float *q, float *r)
{
    float a, b, c, S;

    a = _dxfDistance(p, q);
    b = _dxfDistance(q, r);
    c = _dxfDistance(r, p);

    S = (a + b + c)/ 2;

    return sqrt(S*(S-a)*(S-b)*(S-c));
}
	
static float 
_dxfIntersect(float *center, float radius, float *p0, float *p1)
{
    float x10, y10, z10;
    float x0c, y0c, z0c;
    float A, B, C, D;
    float t;

    x10 = p1[0] - p0[0];
    y10 = p1[1] - p0[1];
    z10 = p1[2] - p0[2];

    x0c = p0[0] - center[0];
    y0c = p0[1] - center[1];
    z0c = p0[2] - center[2];

    A = x10*x10 + y10*y10 + z10*z10;
    B = 2 *(x10*x0c + y10*y0c + z10*z0c);
    C = x0c*x0c + y0c*y0c + z0c*z0c - radius*radius;

    D = B*B - 4*A*C;

    if (D >= 0.0)
    {
	D = sqrt(D);

	t = (-B + D)/(2 * A);
	if (t >= -0.0001 && t <= 1.0001)
	    return t;

	t = (-B - D)/(2 * A);
	if (t >= -0.0001 && t <= 1.0001)
	    return t;
    }
    
    return -1.0;
}

static int 
CompareDependency(Object o, char *target)
{
    Class class = DXGetObjectClass(o);
    if (class == CLASS_GROUP)
	class = DXGetGroupClass((Group)o);

    if (class == CLASS_FIELD)
    {
	Object attr;

	if (! DXGetComponentValue((Field)o, "normals"))
	    return 0;

	attr = DXGetComponentAttribute((Field)o, "normals", "dep");
	if (! attr)
	{
	    DXSetError(ERROR_INVALID_DATA, 
		"normals component is missing the dep attribute");
	    return -1;
	}

	if (DXGetObjectClass(attr) != CLASS_STRING)
	{
	    DXSetError(ERROR_BAD_CLASS, 
		"dep attribute must be class STRING");
	    return -1;
	}

	return !strcmp(DXGetString((String)attr), target);
    }
    else if (class == CLASS_COMPOSITEFIELD)
    {
	Object child;
	int i, j;

	i = 0;
	while(NULL != (child = DXGetEnumeratedMember((Group)o, i++, NULL)))
	{
	    j = CompareDependency(child, target);
	    if (j != 1)
		return j;
	}

	return 1;
    }
    else
    {
	DXSetError(ERROR_BAD_CLASS, "field or composite field");
	return -1;
    }
}

static Error
_dxfFaceNormals(Pointer ptr)
{
    Field field = (Field)ptr;
    Array fa, la, ea, pa, na;
    int *f, *l, *e;
    float *p;
    int nf, nl, ne, ndim, rank;

    na = (Array)DXGetComponentValue(field, "normals");
    if (na)
    {
	char *str;

	if (! DXGetStringAttribute((Object)na, "dep", &str))
	{
	    DXSetError(ERROR_INVALID_DATA, "normals dependency attribute");
	    goto error;
	}

	if (! strcmp(str, "faces"))
	    return OK;
	
	na = NULL;
    }

    fa   = (Array)DXGetComponentValue(field, "faces");
    la   = (Array)DXGetComponentValue(field, "loops");
    ea   = (Array)DXGetComponentValue(field, "edges");
    pa   = (Array)DXGetComponentValue(field, "positions");
    if (! fa || ! la || ! ea || ! pa)
    {
	DXSetError(ERROR_MISSING_DATA, "one of faces/loops/edges/positions");
	goto error;
    }

    f = (int *)DXGetArrayData(fa);
    DXGetArrayInfo(fa, &nf, NULL, NULL, NULL, NULL);

    l = (int *)DXGetArrayData(la);
    DXGetArrayInfo(la, &nl, NULL, NULL, NULL, NULL);

    e = (int *)DXGetArrayData(ea);
    DXGetArrayInfo(ea, &ne, NULL, NULL, NULL, NULL);

    p = (float *)DXGetArrayData(pa);
    DXGetArrayInfo(pa, NULL, NULL, NULL, &rank, &ndim);

    if (rank != 1 || ndim < 2 || ndim > 3)
    {
	DXSetError(ERROR_INVALID_DATA, "data must be 2- or 3-D");
	goto error;
    }
    else
    {
	na = _dxf_FLE_Normals(f, nf, l, nl, e, ne, p, ndim);
    }

    if (! na)
	goto error;

    if (! DXSetStringAttribute((Object)na, "dep", "faces"))
	goto error;

    if (! DXSetComponentValue(field, "normals", (Object)na))
	goto error;
    
    return OK;

error:
    DXDelete((Object)na);
    return ERROR;
}

/*
 * This is the interface called in triangulation code
 */
Array
_dxf_FLE_Normals(int *f, int nf, int *l, int nl, int *e, int ne, float *v, int nd)
{
    int   i, face, maxloop;
    Point *normal;
    float **ptrs = NULL;
    Array nA = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
    if (! nA)
	goto error;
    
    if (! DXAddArrayData(nA, 0, nf, NULL))
	goto error;
    
    normal = (Point *)DXGetArrayData(nA);
    
    /*
     * For each face, find the longest outer loop
     */
    maxloop = 0;
    for (face = 0; face < nf; face ++)
    {
	int loop   = f[face];
	int start  = l[loop];
	int end    = (loop == nl-1) ? ne : l[loop+1];
	int knt    = end - start;

	if (knt > maxloop)
	    maxloop = knt;
    }

    ptrs = (float **)DXAllocate(maxloop*sizeof(float *));
    if (! ptrs)
	goto error;

    /*
     * For each face...
     */
    for (face = 0; face < nf; face++, normal++)
    {
	int   loop   = f[face];
	int   start  = l[loop];
	int   end    = (loop == nl-1) ? ne : l[loop+1];
	int   knt    = end - start;
	float l;

	for (i = 0; i < knt; i++)
	    ptrs[i] = v + e[i+start]*nd;
	
	NORMAL(knt, nd, normal);

	l = DXLength(*normal);
	if (l != 0.0)
	    *normal = DXMul(*normal, l);

	continue;

degenerate:
	normal->x = 0.0;
	normal->y = 0.0;
	normal->z = 1.0;
	continue;
    }

    DXFree((Pointer)ptrs);

    return nA;

error:
    DXDelete((Object)nA);
    DXFree((Pointer)ptrs);
    return NULL;
}

static Field
GridNormalsField(Field field, int dep)
{
    Array pArray, cArray, nArray;
    int pDim, cDim, n, counts[3];
    float x, y, z, d, *vectors[2], deltas[9], norm[3];
    int i, j;
    Object attr;

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

    if (! pArray || ! cArray)
	return NULL;
    
    if (! DXQueryGridConnections(cArray, &cDim, NULL) ||
	! DXQueryGridPositions(pArray, &pDim, counts, NULL, deltas))
	return NULL;
    
    if (cDim != 2 || pDim < 2 || pDim > 3)
	return NULL;
    
    for (i = j = 0; i < pDim; i++)
	if (counts[i] > 1)
	    vectors[j++] = deltas + i*pDim;
    
    if (j > 2)
	return NULL;
    
    if (pDim == 3)
    {
	x = vectors[0][1]*vectors[1][2] - vectors[0][2]*vectors[1][1];
	y = vectors[0][2]*vectors[1][0] - vectors[0][0]*vectors[1][2];
	z = vectors[0][0]*vectors[1][1] - vectors[0][1]*vectors[1][0];
    }
    else
    {
	x = 0.0;
	y = 0.0;
	z = vectors[0][0]*vectors[1][1] - vectors[0][1]*vectors[1][0];
    }

    d = x*x + y*y + z*z;
    if (d != 0.0)
    {
	d = 1.0 / sqrt(d);

	norm[0] = x * d;
	norm[1] = y * d;
	norm[2] = z * d;
    }
    else norm[0] = norm[1] = norm[2] = 0.0;

    if (dep == DEP_ON_POSITIONS)
    {
	DXGetArrayInfo(pArray, &n, NULL, NULL, NULL, NULL);
	attr = (Object)DXNewString("positions");
    }
    else if (dep == DEP_ON_CONNECTIONS)
    {
	DXGetArrayInfo(cArray, &n, NULL, NULL, NULL, NULL);
	attr = (Object)DXNewString("connections");
    }
    else
	return NULL;

    if (! attr)
	return NULL;

    nArray = (Array)DXNewConstantArray(n, (Pointer)&norm, TYPE_FLOAT,
			CATEGORY_REAL, 1, 3);
    if (! nArray)
    {
	DXDelete((Object)attr);
	return NULL;
    }
    
    if (DXGetComponentValue(field, "normals"))
	DXDeleteComponent(field, "normals");
    
    if (! DXSetComponentValue(field, "normals", (Object)nArray))
    {
	DXDelete((Object)nArray);
	DXDelete((Object)attr);
	return NULL;
    }

    if (! DXSetComponentAttribute(field, "normals", "dep", attr))
    {
	DXDelete((Object)nArray);
	DXDelete((Object)attr);
	return NULL;
    }

    return field;
}

static CompositeField
GridNormalsCompositeField(CompositeField o, int dep)
{
    int i;
    Object c;
    
    i = 0;
    while(NULL != (c = DXGetEnumeratedMember((Group)o, i++, NULL)))
	if (! GridNormalsField((Field)c, dep))
	    return NULL;
    
    return o;
}
	 
