/******************************************************************************
* CagdCMrg.c - Surface/Surface merging routine.				      *
*******************************************************************************
* Written by Gershon Elber, Sep 92.					      *
******************************************************************************/

#include "cagd_loc.h"

static void InterpolateLinearSeg(CagdRType *V1, CagdRType *V2, int Len,
								int Step);

/******************************************************************************
* Merge two surfaces in the provided direction Dir. If SameEdge, it is        *
* assumed last edge of Srf1 is identical to first edge of Srf2 and one row    *
* is dropped from new mesh. Otherwise a ruled surface is fit between the two  *
* edges.								      *
******************************************************************************/
CagdSrfStruct *CagdMergeSrfSrf(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2,
			       CagdSrfDirType Dir, CagdBType SameEdge)
{
    CagdBType IsNotRational;
    int i, j, UOrder, VOrder, ULen1, VLen1, ULen2, VLen2, MaxCoord, Length;
    CagdRType **Points1, **Points2, **Points;
    CagdPointType SrfPType;
    CagdSrfStruct *Srf;

    /* Make surfaces compatible: */
    Srf1 = CagdSrfCopy(Srf1);
    Srf2 = CagdSrfCopy(Srf2);

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	case CAGD_CONST_V_DIR:
	    if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, TRUE, TRUE,
					Dir == CAGD_CONST_V_DIR,
					Dir == CAGD_CONST_U_DIR))
		FATAL_ERROR(CAGD_ERR_SRF_FAIL_CMPT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    /* Verify surface geometric types. */
    switch (Srf1 -> GType) {
	case CAGD_SBEZIER_TYPE:
	    Srf = CnvrtBezier2BsplineSrf(Srf1);
	    CagdSrfFree(Srf1);
	    Srf1 = Srf;
	    Srf = CnvrtBezier2BsplineSrf(Srf2);
	    CagdSrfFree(Srf2);
	    Srf2 = Srf;
	    break;
	case CAGD_SBSPLINE_TYPE:
	    break;
	case CAGD_SPOWER_TYPE:
	    FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_SRF);
	    break;
    }

    UOrder = Srf1 -> UOrder;
    VOrder = Srf1 -> VOrder;
    ULen1 = Srf1 -> ULength;
    VLen1 = Srf1 -> VLength;
    ULen2 = Srf2 -> ULength;
    VLen2 = Srf2 -> VLength;
    Points1 = Srf1 -> Points;
    Points2 = Srf2 -> Points;
    IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf1);
    MaxCoord = CAGD_NUM_OF_PT_COORD(Srf1 -> PType),
    SrfPType = CAGD_MAKE_PT_TYPE(!IsNotRational, MaxCoord);

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    Length = SameEdge ? ULen1 + ULen2 - 1 : ULen1 + ULen2 + UOrder - 2;
	    Srf = BspSrfNew(Length, VLen1, UOrder, VOrder, SrfPType);

	    /* Update knot vectors. We assume open end condition here... */
	    CAGD_GEN_COPY(Srf -> UKnotVector, Srf1 -> UKnotVector,
			  (ULen1 + UOrder - 1) * sizeof(CagdRType));
	    CAGD_GEN_COPY(Srf -> VKnotVector, Srf1 -> VKnotVector,
			  (VLen1 + VOrder) * sizeof(CagdRType));

	    if (SameEdge) {
		/* Copy kv of second surface immediately after. */
		CAGD_GEN_COPY(&Srf -> UKnotVector[ULen1 + UOrder - 1],
			      &Srf2 -> UKnotVector[UOrder],
			      ULen2 * sizeof(CagdRType));
		BspKnotAffineTrans(&Srf -> UKnotVector[ULen1 + UOrder - 1],
				   ULen2,
				   Srf -> UKnotVector[ULen1 + UOrder - 2] -
				   Srf2 -> UKnotVector[0],
				   1.0);
	    }
	    else {
		/* Copy kv of second surface order after. */
		CAGD_GEN_COPY(&Srf -> UKnotVector[ULen1 + UOrder - 1],
			      &Srf2 -> UKnotVector[1],
			      (ULen2 + UOrder - 1) * sizeof(CagdRType));
		BspKnotAffineTrans(&Srf -> UKnotVector[ULen1 + UOrder - 1],
				   ULen2 + UOrder - 1,
				   Srf -> UKnotVector[ULen1 + UOrder - 2] -
				   Srf -> UKnotVector[ULen1 + UOrder - 1] + 1.0,
				   1.0);
	    }

	    Points = Srf -> Points;

	    for (i = 0; i < VLen1; i++) {
		for (j = IsNotRational; j <= MaxCoord; j++) {
		    CAGD_GEN_COPY(&Points[j][CAGD_MESH_UV(Srf, 0, i)],
				  &Points1[j][CAGD_MESH_UV(Srf1, 0, i)],
				  ULen1 * sizeof(CagdRType));
		    if (SameEdge) {
			/* Copy row of second surface immediately after. */
			CAGD_GEN_COPY(
			    &Points[j][CAGD_MESH_UV(Srf, ULen1 - 1, i)],
			    &Points2[j][CAGD_MESH_UV(Srf2, 0, i)],
			    ULen2 * sizeof(CagdRType));
		    }
		    else {
			/* Copy row of 2nd srf order after and lin. interp. */
			CAGD_GEN_COPY(
			    &Points[j][CAGD_MESH_UV(Srf, ULen1 + UOrder - 2, i)],
			    &Points2[j][CAGD_MESH_UV(Srf2, 0, i)],
			    ULen2 * sizeof(CagdRType));
			InterpolateLinearSeg(
			    &Points[j][CAGD_MESH_UV(Srf, ULen1 - 1, i)],
			    &Points[j][CAGD_MESH_UV(Srf, ULen1 + UOrder - 2, i)],
			    UOrder, 1);
		    }
		}
	    }
	    break;
	case CAGD_CONST_V_DIR:
	    Length = SameEdge ? VLen1 + VLen2 - 1 : VLen1 + VLen2 + VOrder - 2;
	    Srf = BspSrfNew(ULen1, Length, UOrder, VOrder, SrfPType);

	    /* Update knot vectors. We assume open end condition here... */
	    CAGD_GEN_COPY(Srf -> UKnotVector, Srf1 -> UKnotVector,
			  (ULen1 + UOrder) * sizeof(CagdRType));
	    CAGD_GEN_COPY(Srf -> VKnotVector, Srf1 -> VKnotVector,
			  (VLen1 + VOrder - 1) * sizeof(CagdRType));

	    if (SameEdge) {
		/* Copy kv of second surface immediately after. */
		CAGD_GEN_COPY(&Srf -> VKnotVector[VLen1 + VOrder - 1],
			      &Srf2 -> VKnotVector[VOrder],
			      VLen2 * sizeof(CagdRType));
		BspKnotAffineTrans(&Srf -> VKnotVector[VLen1 + VOrder - 1],
				   VLen2,
				   Srf -> VKnotVector[VLen1 + VOrder - 2] -
				   Srf2 -> VKnotVector[0],
				   1.0);
	    }
	    else {
		/* Copy kv of second surface order after. */
		CAGD_GEN_COPY(&Srf -> VKnotVector[VLen1 + VOrder - 1],
			      &Srf2 -> VKnotVector[1],
			      (VLen2 + VOrder - 1) * sizeof(CagdRType));
		BspKnotAffineTrans(&Srf -> VKnotVector[VLen1 + VOrder - 1],
				   VLen2 + VOrder - 1,
				   Srf -> VKnotVector[VLen1 + VOrder - 2] -
				   Srf -> VKnotVector[VLen1 + VOrder - 1] + 1.0,
				   1.0);
	    }

	    Points = Srf -> Points;

	    for (i = 0; i < ULen1; i++) {
		for (j = IsNotRational; j <= MaxCoord; j++) {
		    CAGD_GEN_COPY_STEP(&Points[j][CAGD_MESH_UV(Srf, i, 0)],
				       &Points1[j][CAGD_MESH_UV(Srf1, i, 0)],
				       VLen1, ULen1, ULen1, CagdRType);
		    if (SameEdge) {
			/* Copy col of second surface immediately after. */
			CAGD_GEN_COPY_STEP(
			    &Points[j][CAGD_MESH_UV(Srf, i, VLen1 - 1)],
			    &Points2[j][CAGD_MESH_UV(Srf2, i, 0)],
			    VLen2, ULen1, ULen1, CagdRType);
		    }
		    else {
			/* Copy col of 2nd srf order after and lin. interp. */
			CAGD_GEN_COPY_STEP(
			    &Points[j][CAGD_MESH_UV(Srf, i, VLen1 + VOrder - 2)],
			    &Points2[j][CAGD_MESH_UV(Srf2, i, 0)],
			    VLen2, ULen1, ULen1, CagdRType);
			InterpolateLinearSeg(
			    &Points[j][CAGD_MESH_UV(Srf, i, VLen1 - 1)],
			    &Points[j][CAGD_MESH_UV(Srf, i, VLen1 + VOrder - 2)],
			    VOrder, ULen1);
		    }
		}
	    }
	    break;
	default:
	    Srf = NULL;
	    FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    CagdSrfFree(Srf1);
    CagdSrfFree(Srf2);

    return Srf;
}

/******************************************************************************
* Merge a list of surfaces, in order.					      *
******************************************************************************/
CagdSrfStruct *CagdMergeSrfList(CagdSrfStruct *SrfList, CagdSrfDirType Dir,
							CagdBType SameEdge)
{
    if (SrfList != NULL && SrfList -> Pnext != NULL) {
	CagdSrfStruct
	    *MergedSrf = CagdSrfCopy(SrfList);

	for (SrfList = SrfList -> Pnext;
	     SrfList != NULL;
	     SrfList = SrfList -> Pnext) {
	    CagdSrfStruct
		*TmpSrf = CagdMergeSrfSrf(MergedSrf, SrfList, Dir, SameEdge);

	    CagdSrfFree(MergedSrf);
	    MergedSrf = TmpSrf;
	}
	return MergedSrf;
    }
    else
	return SrfList ? CagdSrfCopy(SrfList) : NULL;
}

/******************************************************************************
* Linearly interpolate between V1 and V2 values Len times (Len includes V1    *
* and V2) and step the array using Step.				      *
******************************************************************************/
static void InterpolateLinearSeg(CagdRType *V1, CagdRType *V2, int Len,
								int Step)
{
    int i;
    CagdRType
	*V = V1 + Step;

    if (Len-- <= 2)
	return;				     /* No middle points to interp. */

    for (i = 1; i < Len; i++) {
	*V = (i * (*V2) + (Len - i) * (*V1)) / Len;

	V += Step;
    }
}