/******************************************************************************
* SBsp-Int.c - Bspline surface interpolation/approximation schemes.	      *
*******************************************************************************
* Written by Gershon Elber, Feb. 94.					      *
******************************************************************************/

#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include "cagd_loc.h"
#include "extra_fn.h"

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a set of points, PtList, computes a Bspline surface of order UOrder  M
* by VOrder that interpolates or least square approximates the given set of  M
* points.	                                                             M
*   PtList is a NULL terminated array of linked lists of CagdPtStruct        M
* structs.                                                 	             M
*   All linked lists in PtList must have the same length.		     M
*   U direction of surface is associated with array, V with the linked       M
* lists.		                                   	             M
*   The size of the control mesh of the resulting Bspline surface defaults   M
* to the number of points in PtList (if SrfUSize = SrfVSize = 0).	     M
*   However, either numbers can smaller to yield a least square              M
* approximation of the gievn data set.                     	             M
*   The created curve can be parametrized as specified by ParamType.         M
*                                                                            *
* PARAMETERS:                                                                M
*   PtList:      A NULL terminating array of linked list of points.          M
*   UOrder:      Of the to be created surface.                               M
*   VOrder:      Of the to be created surface.                               M
*   SrfUSize:    U size of the to be created surface. Must be at least as    M
*                large as the array PtList.			             M
*   SrfVSize:    V size of the to be created surface. Must be at least as    M
*                large as the length of each list in PtList.                 M
*   ParamType:    Type of parametrization.                                   M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   Constructed interpolating/approximating surface.      M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfInterpPts, interpolation, least square approximation               M
*****************************************************************************/
CagdSrfStruct *BspSrfInterpPts(CagdPtStruct **PtList,
			       int UOrder,
			       int VOrder,
			       int SrfUSize,
			       int SrfVSize,
			       CagdParametrizationType ParamType)
{
    int i, NumLinkedLists, NumPtsInList;
    CagdPtStruct *Pt, **PtPtr;
    CagdRType *UKV, *VKV, *PtUKnots, *PtVKnots, *AveKV, *RU, *RV, *R2;
    CagdSrfStruct *Srf;
    CagdCtlPtStruct
	*CtlPt = NULL,
	*CtlPtList = NULL;

    for (NumLinkedLists = 0, PtPtr = PtList;
	 *PtPtr != NULL;
	 NumLinkedLists++, PtPtr++);
    for (NumPtsInList = 0, Pt = *PtList;
	 Pt != NULL;
	 NumPtsInList++, Pt = Pt -> Pnext);

    if (SrfUSize == 0)
	SrfUSize = NumLinkedLists;
    if (SrfVSize == 0)
	SrfVSize = NumPtsInList;
    if (NumLinkedLists < 3 ||
	NumLinkedLists < UOrder ||
	NumLinkedLists < SrfUSize ||
	SrfUSize < UOrder ||
	NumPtsInList < 3 ||
	NumPtsInList < VOrder ||
	NumPtsInList < SrfVSize ||
	SrfVSize < VOrder)
	return NULL;

    RU = PtUKnots =
	(CagdRType *) IritMalloc(sizeof(CagdRType) * NumLinkedLists);
    RV = PtVKnots =
	(CagdRType *) IritMalloc(sizeof(CagdRType) * NumPtsInList);

    /* Compute parameter values at which surface will interpolate PtList. */
    switch (ParamType) {
	case CAGD_CHORD_LEN_PARAM:
        case CAGD_CENTRIPETAL_PARAM:
	    /* No real support for chord length although we might be able to */
	    /* do something useful better that uniform parametrization.      */
	case CAGD_UNIFORM_PARAM:
	default:
	    for (i = 0; i < NumLinkedLists; i++)
		*RU++ = ((RealType) i) / (NumLinkedLists - 1);
	    for (i = 0; i < NumPtsInList; i++)
		*RV++ = ((RealType) i) / (NumPtsInList - 1);
	    break;
    }

    /* Construct the two knot vectors of the Bspline surface. */
    UKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * (SrfUSize + UOrder));
    VKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * (SrfVSize + VOrder));

    AveKV = BspKnotAverage(PtUKnots, NumLinkedLists,
			   NumLinkedLists + UOrder - SrfUSize - 1);
    BspKnotAffineTrans(AveKV, SrfUSize - UOrder + 2, PtUKnots[0] - AveKV[0],
		       (PtUKnots[NumLinkedLists - 1] - PtUKnots[0]) /
			   (AveKV[SrfUSize - UOrder + 1] - AveKV[0]));
    for (i = 0, RU = UKV, R2 = AveKV; i < UOrder; i++)
	*RU++ = *R2;
    for (i = 0; i < SrfUSize - UOrder; i++)
	*RU++ = *++R2;
    for (i = 0, R2++; i < UOrder; i++)
	*RU++ = *R2;
    IritFree((VoidPtr) AveKV);

    AveKV = BspKnotAverage(PtVKnots, NumPtsInList,
			   NumPtsInList + VOrder - SrfVSize - 1);
    BspKnotAffineTrans(AveKV, SrfVSize - VOrder + 2, PtVKnots[0] - AveKV[0],
		       (PtVKnots[NumPtsInList - 1] - PtVKnots[0]) /
			   (AveKV[SrfVSize - VOrder + 1] - AveKV[0]));
    for (i = 0, RV = VKV, R2 = AveKV; i < VOrder; i++)
	*RV++ = *R2;
    for (i = 0; i < SrfVSize - VOrder; i++)
	*RV++ = *++R2;
    for (i = 0, R2++; i < VOrder; i++)
	*RV++ = *R2;
    IritFree((VoidPtr) AveKV);

    /* Convert to control points in a linear list */
    for (PtPtr = PtList; *PtPtr != NULL; PtPtr++) {
	for (Pt = *PtPtr, i = 0; Pt != NULL; Pt = Pt -> Pnext, i++) {
	    int j;

	    if (CtlPtList == NULL)
		CtlPtList = CtlPt = CagdCtlPtNew(CAGD_PT_E3_TYPE);
	    else {
		CtlPt ->Pnext = CagdCtlPtNew(CAGD_PT_E3_TYPE);
		CtlPt = CtlPt -> Pnext;
	    }
	    for (j = 0; j < 3; j++)
		CtlPt -> Coords[j + 1] = Pt -> Pt[j];
	}
	if (i != NumPtsInList) {
	    CagdCtlPtFreeList(CtlPtList);

	    IritFree((VoidPtr *) PtUKnots);
	    IritFree((VoidPtr *) PtVKnots);
	    IritFree((VoidPtr *) UKV);
	    IritFree((VoidPtr *) VKV);
    
	    CAGD_FATAL_ERROR(CAGD_ERR_PT_OR_LEN_MISMATCH);
	    return NULL;
	}
    }

    Srf = BspSrfInterpolate(CtlPtList, NumLinkedLists, NumPtsInList,
			    PtUKnots, PtVKnots, UKV, VKV,
			    SrfUSize, SrfVSize, UOrder, VOrder);
    CagdCtlPtFreeList(CtlPtList);

    IritFree((VoidPtr *) PtUKnots);
    IritFree((VoidPtr *) PtVKnots);
    IritFree((VoidPtr *) UKV);
    IritFree((VoidPtr *) VKV);

    return Srf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a set of points on a rectangular grid, PtList, parameter values the  M
* surface should interpolate or approximate these grid points, U/VParams,    M
* the expected two knot vectors of the surface, U/VKV, the expected lengths  M
* U/VLength and orders U/VOrder of the Bspline surface, computes the Bspline M
* surface's coefficients.						     M
*   All points in PtList are assumed of the same type.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   PtList:     A long linked list (NumUPts * NumVPts) of points to          M
*               interpolated or least square approximate.                    M
*   NumUPts:    Number of points in PtList in the U direction.               M
*   NumVPts:    Number of points in PtList in the V direction.               M
*   UParams:    Parameter at which surface should interpolate or             M
*               approxximate PtList in the U direction.                      M
*   VParams:    Parameter at which surface should interpolate or             M
*               approxximate PtList in the V direction.                      M
*   UKV:        Requested knot vector form the surface in the U direction.   M
*   VKV:        Requested knot vector form the surface in the V direction.   M
*   ULength:    Requested length of control mesh of surface in U direction.  M
*   VLength:    Requested length of control mesh of surface in V direction.  M
*   UOrder:     Requested order of surface in U direction.                   M
*   VOrder:     Requested order of surface in U direction.                   M
*									     *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   Constructed interpolating/approximating surface.      M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfInterpolate, interpolation, least square approximation             M
*****************************************************************************/
CagdSrfStruct *BspSrfInterpolate(CagdCtlPtStruct *PtList,
				 int NumUPts,
				 int NumVPts,
				 CagdRType *UParams,
				 CagdRType *VParams,
				 CagdRType *UKV,
				 CagdRType *VKV,
				 int ULength,
				 int VLength,
				 int UOrder,
				 int VOrder)
{
    CagdPointType
	PtType = PtList -> PtType;
    CagdBType
	IsRational = CAGD_IS_RATIONAL_PT(PtType);
    int i, j,
	NumCoords = CAGD_NUM_OF_PT_COORD(PtType);
    CagdCtlPtStruct *Pt;
    CagdRType **SPoints;
    CagdSrfStruct *Srf;
    CagdCrvStruct **Crvs;

    /* Construct the Bspline surface and its two knot vectors. */
    Srf = BspSrfNew(ULength, VLength, UOrder, VOrder, PtType);
    SPoints = Srf -> Points;
    CAGD_GEN_COPY(Srf -> UKnotVector, UKV,
		  (ULength + UOrder) * sizeof(CagdRType));
    CAGD_GEN_COPY(Srf -> VKnotVector, VKV,
		  (VLength + VOrder) * sizeof(CagdRType));

    /* Interpolate the rows as set of curves. */
    Crvs = (CagdCrvStruct **) IritMalloc(sizeof(CagdCrvStruct *) * NumVPts);
    for (i = 0, Pt = PtList; i < NumVPts; i++) {
	Crvs[i] = BspCrvInterpolate(Pt, NumUPts, UParams, UKV, ULength,
				    UOrder, FALSE);

	for (j = 0; j < NumUPts; j++)		     /* Advance to next row. */
	    Pt = Pt -> Pnext;
    }

    /* Interpolate the columns from the curves of the rows. */
    for (i = 0; i < ULength; i++) {
	int k;
	CagdCrvStruct *Crv;
	CagdCtlPtStruct
	    *CtlPtList = NULL,
	    *CtlPt = NULL;
	CagdRType **CPoints;

	for (j = 0; j < NumVPts; j++) {
	    CagdRType
		**Points = Crvs[j] -> Points;

	    if (CtlPtList == NULL)
		CtlPtList = CtlPt = CagdCtlPtNew(CAGD_PT_E3_TYPE);
	    else {
		CtlPt ->Pnext = CagdCtlPtNew(CAGD_PT_E3_TYPE);
		CtlPt = CtlPt -> Pnext;
	    }

	    for (k = !IsRational; k <= NumCoords; k++)
		CtlPt -> Coords[k] = Points[k][i];
	}

	/* Interpolate the column, copy to mesh, and free the curve. */
	Crv = BspCrvInterpolate(CtlPtList, NumVPts, VParams, VKV, VLength,
				VOrder, FALSE);
	CPoints = Crv -> Points;
	CagdCtlPtFreeList(CtlPtList);

	for (j = 0; j < VLength; j++)
	    for (k = !IsRational; k <= NumCoords; k++)
		SPoints[k][i + j * ULength] = CPoints[k][j];

	CagdCrvFree(Crv);
    }

    for (i = 0; i < NumVPts; i++)
	CagdCrvFree(Crvs[i]);
    IritFree((VoidPtr) Crvs);

    return Srf;
}
