/******************************************************************************
* TriInterp.c - Interpolation of trivariate data using trivariates.           *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, Oct. 94.					      *
******************************************************************************/

#include "geomat3d.h"
#include "triv_loc.h"

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Interpolates control points of given trivariate, preserving the order    M
* and continuity of the original trivariate.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   TV:   Trivariate to interpolate its control mesh.                        M
*                                                                            *
* RETURN VALUE:                                                              M
*   TrivTVStruct *: The interpolating trivariate.                            M
*                                                                            M
* KEYWORDS:                                                                  M
*   TrivInterpTrivar, interpolation                                          M
*****************************************************************************/
TrivTVStruct *TrivInterpTrivar(TrivTVStruct *TV)
{
    CagdPointType
	PType = TV -> PType;
    CagdBType
	IsRational = CAGD_IS_RATIONAL_PT(PType);
    int i, j, k, l, m,
	MaxCoord = CAGD_NUM_OF_PT_COORD(PType),
	ULength = TV -> ULength,
	VLength = TV -> VLength,
	WLength = TV -> WLength,
	UOrder = TV -> UOrder,
	VOrder = TV -> VOrder,
	WOrder = TV -> WOrder,
	UVLength = ULength * VLength;
    TrivTVStruct
	*InterpTV = TrivTVCopy(TV);
    CagdRType
        *UKnotVector = TRIV_IS_BSPLINE_TV(TV) ? TV -> UKnotVector
				: BspKnotUniformOpen(ULength, ULength, NULL),
        *VKnotVector = TRIV_IS_BSPLINE_TV(TV) ? TV -> VKnotVector
				: BspKnotUniformOpen(VLength, VLength, NULL),
        *WKnotVector = TRIV_IS_BSPLINE_TV(TV) ? TV -> WKnotVector
				: BspKnotUniformOpen(WLength, WLength, NULL),
	*UNodes = BspKnotNodes(UKnotVector, ULength + UOrder, UOrder),
	*VNodes = BspKnotNodes(VKnotVector, VLength + VOrder, VOrder),
	*WNodes = BspKnotNodes(WKnotVector, WLength + WOrder, WOrder),
	**InterpPoints = InterpTV -> Points,
        **OrigPoints = TV -> Points;
    CagdSrfStruct
	**Srfs = (CagdSrfStruct **) IritMalloc(sizeof(CagdSrfStruct *) *
								     WLength);

    /* Interpolate the WLength surfaces: */
    for (k = 0; k < WLength; k++) {
	int Index = k * UVLength;
	CagdCtlPtStruct *Pt,
	     *PtList = NULL;

	for (l = 0; l < UVLength; l++, Index++) {
	    Pt = CagdCtlPtNew(TV -> PType);
	    for (m = !IsRational; m <= MaxCoord; m++)
		Pt -> Coords[m] = OrigPoints[m][Index];

	    LIST_PUSH(Pt, PtList);
	}
	PtList = CagdListReverse(PtList);

	Srfs[k] = BspSrfInterpolate(PtList, ULength, VLength, UNodes, VNodes,
				   UKnotVector, VKnotVector, ULength, VLength,
				   UOrder, VOrder);

	CagdCtlPtFreeList(PtList);
    }

    /* Interpolate the third, W, dimension. */
    for (i = 0; i < ULength; i++) {
	for (j = 0; j< VLength; j++) {
	    int Index = TRIV_MESH_UVW(TV, i, j, 0);
	    CagdRType **CrvPoints;
	    CagdCrvStruct *Crv;
	    CagdCtlPtStruct *Pt,
	        *PtList = NULL;

	    for (k = 0; k < WLength; k++) {
		CagdRType
		    **SrfPoints = Srfs[k] -> Points;

		Pt = CagdCtlPtNew(TV -> PType);
		for (m = !IsRational; m <= MaxCoord; m++)
		  Pt -> Coords[m] = SrfPoints[m][Index];

		LIST_PUSH(Pt, PtList);
	    }
	    PtList = CagdListReverse(PtList);

	    Crv = BspCrvInterpolate(PtList, WLength, WNodes, WKnotVector,
				    WLength, WOrder, FALSE);
	    CrvPoints = Crv -> Points;

	    for (k = 0; k < WLength; k++) {
		for (m = !IsRational; m <= MaxCoord; m++)
		   InterpPoints[m][Index + k * UVLength] = CrvPoints[m][k];
	    }

	    CagdCrvFree(Crv);
	}
    }

    IritFree(UNodes);
    IritFree(VNodes);
    IritFree(WNodes);
    if (!TRIV_IS_BSPLINE_TV(TV)) {
	IritFree(UKnotVector);
	IritFree(VKnotVector);
	IritFree(WKnotVector);
    }
    for (i = 0; i < WLength; i++)
	CagdSrfFree(Srfs[i]);
    IritFree(Srfs);

    return InterpTV;    
}