/******************************************************************************
* An interface to an implementation of the marching cube algorithm.	      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
*						Gershon Elber, Dec 1992.      *
******************************************************************************/

#include <stdio.h>
#include "irit_sm.h"
#include "cagd_lib.h"
#include "iritprsr.h"
#include "allocate.h"
#include "mrchcube.h"
#include "triv_loc.h"

typedef enum {
    INPUT_ASCII = 1,
    INPUT_INTEGER,
    INPUT_LONG,
    INPUT_BYTE,
    INPUT_IRIT_FLOAT,
    INPUT_IRIT_DOUBLE
} ScalarFormatType;

static PointType
    GlblCubeDim = { 1.0, 1.0, 1.0 };
static int
    GlblSkipInputData = 1,
    GlblDataWidth = 100,
    GlblDataDepth = 100,
    GlblDataHeight = 100;
static ScalarFormatType
    GlblInputFormat = INPUT_ASCII;
static CagdRType
    *GlblLayerOne = NULL,
    *GlblLayerTwo = NULL;

static CagdRType GetOneScalar(FILE *Fin);
static MCCubeCornerScalarStruct *GetCubeFile(FILE *f, int FirstTime);
static MCCubeCornerScalarStruct *GetCubeTV(TrivTVStruct *TV, int FirstTime);
static void EstimateGradient(MCCubeCornerScalarStruct *CCS);

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Extract a polygonal iso-surface out of volumetric data file.             M
*                                                                            *
* PARAMETERS:                                                                M
*   FileName:   Containing the volumetric data.                              M
*   DataType:   Type of scalar value in volume.  Can be one of:              M
*		   1 - Regular ascii (seperated by while spaces.             V
*		   2 - Two bytes short integer.				     V
*		   3 - Four bytes long integer.				     V
*		   4 - One byte (char) integer.				     V
*		   5 - Four bytes float.				     V
*		   6 - Eight bytes double.				     V
*   CubeDim:	Width, height, and depth of a single cube, in object space   M
*		coordinates.						     M
*   Width:      Of volumetric data set.                                      M
*   Height:     Of volumetric data set.                                      M
*   Depth:      Of volumetric data set.                                      M
*   SkipFactor: Typically 1.  For 2, only every second sample is considered  M
*		and for i, only every i'th sample is considered, in all axes.M
*   IsoVal:     At which to extract the iso-surface.                         M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *: A polygonal approximation of the iso-surface at IsoVal M
*		computed using the marching cubes algorithm.	             M
*                                                                            *
* SEE ALSO:                                                                  M
*   MCThresholdCube, MCExtractIsoSurface2                                    M
*                                                                            *
* KEYWORDS:                                                                  M
*   MCExtractIsoSurface                                                      M
*****************************************************************************/
IPObjectStruct *MCExtractIsoSurface(char *FileName,
				    int DataType,
				    PointType CubeDim,
				    int Width,
				    int Height,
				    int Depth,
				    int SkipFactor,
				    CagdRType IsoVal)
{
    FILE *f;
    MCCubeCornerScalarStruct *CCS;
    IPPolygonStruct
	*AllPolys = NULL;
    int IsCirc = IritPrsrSetPolyListCirc(FALSE);

    IritPrsrSetPolyListCirc(IsCirc); /* Restore state, now that we know it. */

    if ((f = fopen(FileName, "r")) == NULL) {
        TRIV_FATAL_ERROR(TRIV_ERR_FAIL_READ_FILE);
        return NULL;
    }

    PT_COPY(GlblCubeDim, CubeDim);
    GlblDataWidth = Width;
    GlblDataHeight = Height;
    GlblDataDepth = Depth;
    GlblSkipInputData = SkipFactor;
    GlblInputFormat = (ScalarFormatType) DataType;

    GlblLayerOne = (CagdRType *) IritMalloc(sizeof(CagdRType) * GlblDataWidth *
							       GlblDataHeight);
    GlblLayerTwo = (CagdRType *) IritMalloc(sizeof(CagdRType) * GlblDataWidth *
							       GlblDataHeight);
    GetCubeFile(f, TRUE);                          /* Initialize local info. */

    while ((CCS = GetCubeFile(f, FALSE)) != NULL) {
	MCPolygonStruct
	    *MCPolys = MCThresholdCube(CCS, IsoVal);

	while (MCPolys != NULL) {
	    int i;
	    MCPolygonStruct
		*MCPolyTmp = MCPolys -> Pnext;

	    /* Convert to regular irit polygons. */
	    for (i = 2; i < MCPolys -> NumOfVertices - 1; i++) {
	        int j;
	        IPVertexStruct
		    *V3 = IPAllocVertex(0, NULL, NULL),
		    *V2 = IPAllocVertex(0, NULL, V3),
		    *V1 = IPAllocVertex(0, NULL, V2);
		IPPolygonStruct
		    *Pl = IPAllocPolygon(0, V1, AllPolys);

		AllPolys = Pl;

		for (j = 0; j < 3; j++) {
		    V1 -> Coord[j] = MCPolys -> V[0][j];
		    V2 -> Coord[j] = MCPolys -> V[i - 1][j];
		    V3 -> Coord[j] = MCPolys -> V[i][j];

		    V1 -> Normal[j] = MCPolys -> N[0][j];
		    V2 -> Normal[j] = MCPolys -> N[i - 1][j];
		    V3 -> Normal[j] = MCPolys -> N[i][j];
		}
		IP_SET_NORMAL_VRTX(V1);
		IP_SET_NORMAL_VRTX(V2);
		IP_SET_NORMAL_VRTX(V3);

		IritPrsrUpdatePolyPlane(Pl);
		if (DOT_PROD(V1 -> Normal, Pl -> Plane) < 0) {
		    PT_SCALE(Pl -> Plane, -1.0);
		    Pl -> Plane[3] *= -1.0;
		}

		if (IsCirc)
		    V3 -> Pnext = V1;
	    }

	    IritFree((VoidPtr) MCPolys);

	    MCPolys = MCPolyTmp;
	}
    }

    IritFree((VoidPtr) GlblLayerOne);
    IritFree((VoidPtr) GlblLayerTwo);

    if (AllPolys != NULL)
        return GenPOLYObject(AllPolys);
    else
        return NULL;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Extract a polygonal iso-surface out of a trivariate function.            M
*                                                                            *
* PARAMETERS:                                                                M
*   TV:		The trivariate to compute an iso surface approximation for,  M
*   Axis:	of the trivariate to handle, 1 for X, 2 for Y, etc.	     M
*   TrivarNormals:   If TRUE normal are computed using gradient of the       M
*		trivariate, if FALSE normal are estimated using differencing.M
*   CubeDim:	Width, height, and depth of a single cube, in object space   M
*		coordinates.						     M
*   SkipFactor: Typically 1.  For 2, only every second sample is considered  M
*		and for i, only every i'th sample is considered, in all axes.M
*   IsoVal:     At which to extract the iso-surface.                         M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *: A polygonal approximation of the iso-surface at IsoVal M
*		computed using the marching cubes algorithm.	             M
*                                                                            *
* SEE ALSO:                                                                  M
*   MCThresholdCube, MCExtractIsoSurface                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   MCExtractIsoSurface2                                                     M
*****************************************************************************/
IPObjectStruct *MCExtractIsoSurface2(TrivTVStruct *TV,
				     int Axis,
				     CagdBType TrivarNormals,
				     PointType CubeDim,
				     int SkipFactor,
				     CagdRType IsoVal)
{
    MCCubeCornerScalarStruct *CCS;
    TrivTVStruct *TVGradient[3];
    IPPolygonStruct
	*AllPolys = NULL;
    int i,
	IsCirc = IritPrsrSetPolyListCirc(FALSE);
    CagdRType UMin, UMax, VMin, VMax, WMin, WMax;

    IritPrsrSetPolyListCirc(IsCirc); /* Restore state, now that we know it. */

    if (Axis < 1 || Axis > CAGD_MAX_PT_COORD || TV -> Points[Axis] == NULL) {
	TRIV_FATAL_ERROR(TRIV_ERR_INVALID_AXIS);
	return NULL;
    }

    /* Make sure the domain is consistent with the data size. */
    TV = TrivTVCopy(TV);
    TrivTVDomain(TV, &UMin, &UMax, &VMin, &VMax, &WMin, &WMax);
    if (!APX_EQ(UMin, 0.0) || !APX_EQ(UMax, TV -> ULength - 1) ||
	!APX_EQ(VMin, 0.0) || !APX_EQ(VMax, TV -> VLength - 1) ||
	!APX_EQ(WMin, 0.0) || !APX_EQ(WMax, TV -> WLength - 1)) {
	BspKnotAffineTrans2(TV -> UKnotVector, TV -> ULength + TV -> UOrder,
			    0.0, TV -> ULength - 1.0);
	BspKnotAffineTrans2(TV -> VKnotVector, TV -> VLength + TV -> VOrder,
			    0.0, TV -> VLength - 1.0);
	BspKnotAffineTrans2(TV -> WKnotVector, TV -> WLength + TV -> WOrder,
			    0.0, TV -> WLength - 1.0);
    }

    if (TrivarNormals &&
	((TVGradient[0] = TrivTVDerive(TV, TRIV_CONST_U_DIR)) == NULL ||
	 (TVGradient[1] = TrivTVDerive(TV, TRIV_CONST_V_DIR)) == NULL ||
	 (TVGradient[2] = TrivTVDerive(TV, TRIV_CONST_W_DIR)) == NULL)) {
	TrivTVFree(TV);
	return NULL;
    }

    for (i = 0; i < 3; i++) {
	if (CubeDim[i] < IRIT_UEPS)
	    CubeDim[i] = IRIT_UEPS;
    }

    PT_COPY(GlblCubeDim, CubeDim);
    GlblDataWidth = TV -> ULength;
    GlblDataHeight = TV -> VLength;
    GlblDataDepth = TV -> WLength;
    GlblSkipInputData = SkipFactor;

    GetCubeTV(TV, TRUE);                           /* Initialize local info. */

    while ((CCS = GetCubeTV(TV, FALSE)) != NULL) {
	MCPolygonStruct
	    *MCPolys = MCThresholdCube(CCS, IsoVal);

	while (MCPolys != NULL) {
	    int i;
	    MCPolygonStruct
		*MCPolyTmp = MCPolys -> Pnext;

	    /* Convert to regular irit polygons. */
	    for (i = 2; i < MCPolys -> NumOfVertices - 1; i++) {
	        int j;
	        IPVertexStruct
		    *V3 = IPAllocVertex(0, NULL, NULL),
		    *V2 = IPAllocVertex(0, NULL, V3),
		    *V1 = IPAllocVertex(0, NULL, V2);
		IPPolygonStruct
		    *Pl = IPAllocPolygon(0, V1, AllPolys);

		AllPolys = Pl;

		for (j = 0; j < 3; j++) {
		    V1 -> Coord[j] = MCPolys -> V[0][j];
		    V2 -> Coord[j] = MCPolys -> V[i - 1][j];
		    V3 -> Coord[j] = MCPolys -> V[i][j];
		}

		if (TrivarNormals) {
		    IPVertexStruct *V;

		    /* Compute the gradients for the three vertices. */
		    for (j = 0, V = V1; j < 3; j++, V = V -> Pnext) {
			int l;
			for (l = 0; l < 3; l++) {
			    CagdRType
			        *R = TrivTVEval(TVGradient[l], 
						V -> Coord[0] / CubeDim[0],
						V -> Coord[1] / CubeDim[1],
						V -> Coord[2] / CubeDim[2]);

			    V -> Normal[l] = R[1];
			}
			PT_NORMALIZE(V -> Normal);
		    }
		}
		else {
		    for (j = 0; j < 3; j++) {
			V1 -> Normal[j] = MCPolys -> N[0][j];
			V2 -> Normal[j] = MCPolys -> N[i - 1][j];
			V3 -> Normal[j] = MCPolys -> N[i][j];
		    }
		}
		IP_SET_NORMAL_VRTX(V1);
		IP_SET_NORMAL_VRTX(V2);
		IP_SET_NORMAL_VRTX(V3);

		IritPrsrUpdatePolyPlane(Pl);
		if (DOT_PROD(V1 -> Normal, Pl -> Plane) < 0) {
		    PT_SCALE(Pl -> Plane, -1.0);
		    Pl -> Plane[3] *= -1.0;
		}

		if (IsCirc)
		    V3 -> Pnext = V1;
	    }

	    IritFree((VoidPtr) MCPolys);

	    MCPolys = MCPolyTmp;
	}
    }

    if (TrivarNormals) {
	TrivTVFree(TVGradient[0]);
	TrivTVFree(TVGradient[1]);
	TrivTVFree(TVGradient[2]);
    }
    TrivTVFree(TV);

    if (AllPolys != NULL)
        return GenPOLYObject(AllPolys);
    else {
	fprintf(stderr, "Empty iso surface resulted.\n");
        return NULL;
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Get one scalar value from input stream.                                  *
*                                                                            *
* PARAMETERS:                                                                *
*   Fin:   Input file to read from.                                          *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdRType:   One read scalar value.                                       *
*****************************************************************************/
static CagdRType GetOneScalar(FILE *Fin)
{
    short i;
    long l;
    float f;
    double d;
    CagdRType r;

    switch (GlblInputFormat) {
	case INPUT_ASCII:
	    if (fscanf(Fin, "%lf", &r) != 1)
		return IRIT_INFNTY;
	    break;
	case INPUT_INTEGER:
	    i = getc(Fin);
	    r = getc(Fin) * 256.0 + i;
	    break;
	case INPUT_LONG:
	    if (fread(&l, 4, 1, Fin) != 4)
		return IRIT_INFNTY;
	    r = l;
	    break;
	case INPUT_BYTE:
	    if ((r = getc(Fin)) == EOF)
		return IRIT_INFNTY;
	    break;
	case INPUT_IRIT_FLOAT:
	    if (fread(&f, 4, 1, Fin) != 4)
		return IRIT_INFNTY;
	    r = f;
	    break;
	case INPUT_IRIT_DOUBLE:
	    if (fread(&d, 8, 1, Fin) != 8)
		return IRIT_INFNTY;
	    r = d;
	    break;
	default:
	    fprintf(stderr, "Input format requested not supported.\n");
	    return IRIT_INFNTY;
	    break;
    }

    return r;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Reads inout file and returns one cube at a time.                         *
*                                                                            *
* PARAMETERS:                                                                *
*   FirstTime:    TRUE for first time, FALSE otherwise.                      *
*                                                                            *
* RETURN VALUE:                                                              *
*   MCCubeCornerScalarStruct *:                                              *
*****************************************************************************/
static MCCubeCornerScalarStruct *GetCubeFile(FILE *f, int FirstTime)
{
    static MCCubeCornerScalarStruct CCS;
    static int
	LayerNumber = -1,
	LayerCountX = -1,
	LayerCountY = 0;
    int i, j;
    CagdRType *p;

    if (FirstTime) {
	/* Initialize the CCS constant data */
	CCS.CubeDim[0] = GlblCubeDim[0];
	CCS.CubeDim[1] = GlblCubeDim[1];
	CCS.CubeDim[2] = GlblCubeDim[2];
	LayerNumber = -GlblSkipInputData;
	LayerCountX = -1;
	LayerCountY = 0;

	for (p = GlblLayerTwo, i = 0; i < GlblDataWidth * GlblDataHeight; i++)
	    if ((*p++ = GetOneScalar(f)) == IRIT_INFNTY)
	        return NULL;

	return NULL;
    }

    if (LayerCountX == -1) {		             /* Read the next layer. */
        LayerNumber += GlblSkipInputData;
        if (LayerNumber >= GlblDataDepth - GlblSkipInputData)
	    return NULL;                                            /* Done! */

	p = GlblLayerOne;                            /* Swap the two layers. */
	GlblLayerOne = GlblLayerTwo;
	GlblLayerTwo = p;

	for (j = 0; j < GlblSkipInputData; j++) { /* And get the next layer. */
	    for (p = GlblLayerTwo, i = 0;
		 i < GlblDataWidth * GlblDataHeight;
		 i++)
	        if ((*p++ = GetOneScalar(f)) == IRIT_INFNTY)
		    return NULL;
	}

	LayerCountX = 0;
	LayerCountY = 0;
    }

    CCS.Vrtx0Lctn[0] = LayerCountX * GlblCubeDim[0] / GlblSkipInputData;
    CCS.Vrtx0Lctn[1] = LayerCountY * GlblCubeDim[1] / GlblSkipInputData;
    CCS.Vrtx0Lctn[2] = LayerNumber * GlblCubeDim[2] / GlblSkipInputData;
    CCS.Corners[0] = GlblLayerOne[LayerCountY * GlblDataWidth + LayerCountX];
    CCS.Corners[1] = GlblLayerOne[LayerCountY * GlblDataWidth +
			          LayerCountX + GlblSkipInputData];
    CCS.Corners[2] = GlblLayerOne[(LayerCountY + GlblSkipInputData) * GlblDataWidth +
			          LayerCountX + GlblSkipInputData];
    CCS.Corners[3] = GlblLayerOne[(LayerCountY + GlblSkipInputData) * GlblDataWidth +
			          LayerCountX];
    CCS.Corners[4] = GlblLayerTwo[LayerCountY * GlblDataWidth + LayerCountX];
    CCS.Corners[5] = GlblLayerTwo[LayerCountY * GlblDataWidth +
			          LayerCountX + GlblSkipInputData];
    CCS.Corners[6] = GlblLayerTwo[(LayerCountY + GlblSkipInputData) * GlblDataWidth +
			          LayerCountX + GlblSkipInputData];
    CCS.Corners[7] = GlblLayerTwo[(LayerCountY + GlblSkipInputData) * GlblDataWidth +
			          LayerCountX];

    EstimateGradient(&CCS);

    LayerCountX += GlblSkipInputData;
    if (LayerCountX >= GlblDataWidth - GlblSkipInputData) {
	LayerCountY += GlblSkipInputData;
	LayerCountX = 0;
	if (LayerCountY >= GlblDataHeight - GlblSkipInputData)
	    LayerCountX = -1;                                /* No more data */
    }

    return &CCS;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Reads inout file and returns one cube at a time.                         *
*                                                                            *
* PARAMETERS:                                                                *
*   FirstTime:    TRUE for first time, FALSE otherwise.                      *
*                                                                            *
* RETURN VALUE:                                                              *
*   MCCubeCornerScalarStruct *:                                              *
*****************************************************************************/
static MCCubeCornerScalarStruct *GetCubeTV(TrivTVStruct *TV, int FirstTime)
{
    static MCCubeCornerScalarStruct CCS;
    static int
	LayerNumber = -1,
	LayerCountX = -1,
	LayerCountY = 0;
    int j;

    if (FirstTime) {
	/* Initialize the CCS constant data */
	CCS.CubeDim[0] = GlblCubeDim[0];
	CCS.CubeDim[1] = GlblCubeDim[1];
	CCS.CubeDim[2] = GlblCubeDim[2];
	LayerNumber = -GlblSkipInputData;
	LayerCountX = -1;
	LayerCountY = 0;

	GlblLayerTwo = TV -> Points[1];

	return NULL;
    }

    if (LayerCountX == -1) {		             /* Read the next layer. */
        LayerNumber += GlblSkipInputData;
        if (LayerNumber >= GlblDataDepth - GlblSkipInputData)
	    return NULL;                                            /* Done! */

	GlblLayerOne = GlblLayerTwo;     /* Swap the two layers and advance. */
	for (j = 0; j < GlblSkipInputData; j++)   /* And get the next layer. */
	    GlblLayerTwo = GlblLayerTwo + GlblDataWidth * GlblDataHeight;

	LayerCountX = 0;
	LayerCountY = 0;
    }

    CCS.Vrtx0Lctn[0] = LayerCountX * GlblCubeDim[0] / GlblSkipInputData;
    CCS.Vrtx0Lctn[1] = LayerCountY * GlblCubeDim[1] / GlblSkipInputData;
    CCS.Vrtx0Lctn[2] = LayerNumber * GlblCubeDim[2] / GlblSkipInputData;
    CCS.Corners[0] = GlblLayerOne[LayerCountY * GlblDataWidth + LayerCountX];
    CCS.Corners[1] = GlblLayerOne[LayerCountY * GlblDataWidth +
			          LayerCountX + GlblSkipInputData];
    CCS.Corners[2] = GlblLayerOne[(LayerCountY + GlblSkipInputData) * GlblDataWidth +
			          LayerCountX + GlblSkipInputData];
    CCS.Corners[3] = GlblLayerOne[(LayerCountY + GlblSkipInputData) * GlblDataWidth +
			          LayerCountX];
    CCS.Corners[4] = GlblLayerTwo[LayerCountY * GlblDataWidth + LayerCountX];
    CCS.Corners[5] = GlblLayerTwo[LayerCountY * GlblDataWidth +
			          LayerCountX + GlblSkipInputData];
    CCS.Corners[6] = GlblLayerTwo[(LayerCountY + GlblSkipInputData) * GlblDataWidth +
			          LayerCountX + GlblSkipInputData];
    CCS.Corners[7] = GlblLayerTwo[(LayerCountY + GlblSkipInputData) * GlblDataWidth +
			          LayerCountX];

    EstimateGradient(&CCS);

    LayerCountX += GlblSkipInputData;
    if (LayerCountX >= GlblDataWidth - GlblSkipInputData) {
	LayerCountY += GlblSkipInputData;
	LayerCountX = 0;
	if (LayerCountY >= GlblDataHeight - GlblSkipInputData)
	    LayerCountX = -1;                                /* No more data */
    }

    return &CCS;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Estimates normals based on first order differencing.                     *
*                                                                            *
* PARAMETERS:                                                                *
*   CCS:       A single cube to estimate normals to its vertices.            *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void EstimateGradient(MCCubeCornerScalarStruct *CCS)
{
    CCS -> GradientX[0] =
	CCS -> GradientX[1] =
	    CCS -> Corners[1] - CCS -> Corners[0];
    CCS -> GradientX[2] =
	CCS -> GradientX[3] =
	    CCS -> Corners[2] - CCS -> Corners[3];
    CCS -> GradientX[4] =
	CCS -> GradientX[5] =
	    CCS -> Corners[5] - CCS -> Corners[4];
    CCS -> GradientX[6] =
	CCS -> GradientX[7] =
	    CCS -> Corners[6] - CCS -> Corners[7];
    
    CCS -> GradientY[0] =
	CCS -> GradientY[3] =
	    CCS -> Corners[3] - CCS -> Corners[0];
    CCS -> GradientY[1] =
	CCS -> GradientY[2] =
	    CCS -> Corners[2] - CCS -> Corners[1];
    CCS -> GradientY[4] =
	CCS -> GradientY[7] =
	    CCS -> Corners[7] - CCS -> Corners[4];
    CCS -> GradientY[5] =
	CCS -> GradientY[6] =
	    CCS -> Corners[6] - CCS -> Corners[5];
    
    CCS -> GradientZ[0] =
	CCS -> GradientZ[4] =
	    CCS -> Corners[4] - CCS -> Corners[0];
    CCS -> GradientZ[1] =
	CCS -> GradientZ[5] =
	    CCS -> Corners[5] - CCS -> Corners[1];
    CCS -> GradientZ[2] =
	CCS -> GradientZ[6] =
	    CCS -> Corners[6] - CCS -> Corners[2];
    CCS -> GradientZ[3] =
	CCS -> GradientZ[7] =
	    CCS -> Corners[7] - CCS -> Corners[3];

    CCS -> HasGradient = TRUE;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Loads a volumetric data set as a trivariate function of prescribed       M
* orders.  Uniform open end conditions are created for it.                   M
*                                                                            *
* PARAMETERS:                                                                M
*   FileName:    To load the trivariate data set from.                       M
*   DataType:    Type of scalar value in volume data file.  Can be one of:   M
*		   1 - Regular ascii (seperated by while spaces.             V
*		   2 - Two bytes short integer.				     V
*		   3 - Four bytes long integer.				     V
*		   4 - One byte (char) integer.				     V
*		   5 - Four bytes float.				     V
*		   6 - Eight bytes double.				     V
*   VolSize:     Dimensions of trivariate voluem.                            M
*   Orders:      Orders of the three axis of the volume (in U, V, and W).    M
*                                                                            *
* RETURN VALUE:                                                              M
*   TrivTVStrcut *:  Loaded trivariate, or NULL if error.                    M
*                                                                            *
* KEYWORDS:                                                                  M
*   TrivLoadVolumeIntoTV                                                     M
*****************************************************************************/
TrivTVStruct *TrivLoadVolumeIntoTV(char *FileName,
				   int DataType,
				   VectorType VolSize,
				   VectorType Orders)
{
    int i;
    FILE *f;
    TrivTVStruct
        *TV = TrivBspTVNew((int) VolSize[0],
			   (int) VolSize[1],
			   (int) VolSize[2],
			   (int) Orders[0],
			   (int) Orders[1],
			   (int) Orders[2],
			   CAGD_PT_E1_TYPE);
    CagdRType
	*R = TV -> Points[1];

    BspKnotUniformOpen(TV -> ULength, TV -> UOrder, TV -> UKnotVector);
    BspKnotUniformOpen(TV -> VLength, TV -> VOrder, TV -> VKnotVector);
    BspKnotUniformOpen(TV -> WLength, TV -> WOrder, TV -> WKnotVector);

    BspKnotAffineTrans2(TV -> UKnotVector, TV -> ULength + TV -> UOrder,
			0.0, TV -> ULength - 1.0);
    BspKnotAffineTrans2(TV -> VKnotVector, TV -> VLength + TV -> VOrder,
			0.0, TV -> VLength - 1.0);
    BspKnotAffineTrans2(TV -> WKnotVector, TV -> WLength + TV -> WOrder,
			0.0, TV -> WLength - 1.0);
    
    if ((f = fopen(FileName, "r")) == NULL) {
	TrivTVFree(TV);
	TRIV_FATAL_ERROR(TRIV_ERR_READ_FAIL);
	return NULL;
    }

    GlblInputFormat = (ScalarFormatType) DataType;

    for (i = TV -> ULength * TV -> VLength * TV -> WLength; i > 0; i--) {
	if ((*R++ = GetOneScalar(f)) == IRIT_INFNTY) {
	    TrivTVFree(TV);
	    TRIV_FATAL_ERROR(TRIV_ERR_READ_FAIL);
	    return NULL;
	}
    }

    return TV;
}