/******************************************************************************
* Cagd_gen.c - General routines used by all modules of CAGD_lib.	      *
*******************************************************************************
* Written by Gershon Elber, Mar. 90.					      *
******************************************************************************/

#ifdef __MSDOS__
#include <stdlib.h>
#include <alloc.h>
#endif /* __MSDOS__ */

#include <string.h>
#include "cagd_loc.h"

#ifdef __TURBOC__			     /* Only for TC++ 1.0 and above. */
#define __DEBUG_MALLOC__
#endif /* __TURBOC__ */

#ifdef __DEBUG_FREE_BOGUS__	      /* Detected freeing of bogus pointers. */
#define FREE_TABLE_SIZE 10000
static VoidPtr AllocPtrTable[FREE_TABLE_SIZE];
static int AllocNumPtrs = 0;
#endif /* __DEBUG_FREE_BOGUS__ */

/* Control of the linear order surface convertion to polygon: */
CagdLin2PolyType _CagdLin2Poly = CAGD_ONE_POLY_PER_COLIN;

static void CagdTransform(CagdRType **Points, int Len, int MaxCoord,
	     CagdBType IsNotRational, CagdRType Translate[3], CagdRType Scale);
static void CagdMatTransform(CagdRType **Points, int Len, int MaxCoord,
				  CagdBType IsNotRational, CagdMatStruct *Mat);

/******************************************************************************
* Allocate and reset all slots of a curve structure:			      *
******************************************************************************/
CagdCrvStruct *CagdCrvNew(CagdGeomType GType, CagdPointType PType, int Length)
{
    int i, MaxAxis = CAGD_NUM_OF_PT_COORD(PType);
    CagdCrvStruct *NewCrv = (CagdCrvStruct *) CagdMalloc(sizeof(CagdCrvStruct));

    NewCrv -> GType = GType;
    NewCrv -> PType = PType;
    NewCrv -> Length = Length;
    NewCrv -> Order = 0;
    NewCrv -> KnotVector = NULL;
    NewCrv -> Pnext = NULL;
    NewCrv -> Points[0] = NULL;			    /* The rational element. */

    for (i = !CAGD_IS_RATIONAL_PT(PType); i <= MaxAxis; i++)
	NewCrv -> Points[i] = (CagdRType *) CagdMalloc(sizeof(CagdRType) *
								 Length);

    for (i = MaxAxis + 1; i <= CAGD_MAX_PT_COORD; i++)
	NewCrv -> Points[i] = NULL;

    return NewCrv;
}

/******************************************************************************
* Allocate and reset all slots of a surface structure:			      *
******************************************************************************/
CagdSrfStruct *CagdSrfNew(CagdGeomType GType, CagdPointType PType, int ULength,
								   int VLength)
{
    int i, MaxAxis = CAGD_NUM_OF_PT_COORD(PType);
    CagdSrfStruct *NewSrf = (CagdSrfStruct *) CagdMalloc(sizeof(CagdSrfStruct));

    NewSrf -> GType = GType;
    NewSrf -> PType = PType;
    NewSrf -> ULength = ULength;
    NewSrf -> VLength = VLength;
    NewSrf -> UOrder = 0;
    NewSrf -> VOrder = 0;
    NewSrf -> UKnotVector = NULL;
    NewSrf -> VKnotVector = NULL;
    NewSrf -> Pnext = NULL;
    NewSrf -> Points[0] = NULL;			    /* The rational element. */

    for (i = !CAGD_IS_RATIONAL_PT(PType); i <= MaxAxis; i++)
	NewSrf -> Points[i] = (CagdRType *) CagdMalloc(sizeof(CagdRType) *
							ULength * VLength);

    for (i = MaxAxis + 1; i <= CAGD_MAX_PT_COORD; i++)
	NewSrf -> Points[i] = NULL;

    return NewSrf;
}

/******************************************************************************
* Allocate and reset all slots of a polygon structure:			      *
******************************************************************************/
CagdPolygonStruct *CagdPolygonNew(void)
{
    CagdPolygonStruct *NewPoly = (CagdPolygonStruct *)
    					CagdMalloc(sizeof(CagdPolygonStruct));

    NewPoly -> Pnext = NULL;

    return NewPoly;
}


/******************************************************************************
* Allocate and reset all slots of a polyline structure:			      *
******************************************************************************/
CagdPolylineStruct *CagdPolylineNew(int Length)
{
    CagdPolylineStruct *NewPoly = (CagdPolylineStruct *)
					CagdMalloc(sizeof(CagdPolylineStruct));

    NewPoly -> Pnext = NULL;
    NewPoly -> Polyline = (CagdPtStruct *) CagdMalloc(sizeof(CagdPtStruct) *
								       Length);
    NewPoly -> Length = Length;

    return NewPoly;
}

/******************************************************************************
* Allocate and copy all slots of a curve structure:			      *
******************************************************************************/
CagdCrvStruct *CagdCrvCopy(CagdCrvStruct *Crv)
{
    int i, MaxAxis = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdCrvStruct *NewCrv = (CagdCrvStruct *) CagdMalloc(sizeof(CagdCrvStruct));

    NewCrv -> GType = Crv -> GType;
    NewCrv -> PType = Crv -> PType;
    NewCrv -> Length = Crv -> Length;
    NewCrv -> Order = Crv -> Order;
    if (Crv -> KnotVector != NULL)
	NewCrv -> KnotVector = BspKnotCopy(Crv -> KnotVector,
					   Crv -> Length + Crv -> Order);
    else
	NewCrv -> KnotVector = NULL;
    NewCrv -> Pnext = NULL;
    NewCrv -> Points[0] = NULL;			    /* The rational element. */

    for (i = !CAGD_IS_RATIONAL_PT(Crv -> PType); i <= MaxAxis; i++) {
	NewCrv -> Points[i] = (CagdRType *) CagdMalloc(sizeof(CagdRType) *
							       Crv -> Length);
	GEN_COPY(NewCrv -> Points[i], Crv -> Points[i],
		 sizeof(CagdRType) * Crv -> Length);
    }

    for (i = MaxAxis + 1; i <= CAGD_MAX_PT_COORD; i++)
	NewCrv -> Points[i] = NULL;

    return NewCrv;
}

/******************************************************************************
* Allocate and copy all slots of a surface structure:			      *
******************************************************************************/
CagdSrfStruct *CagdSrfCopy(CagdSrfStruct *Srf)
{
    int i, MaxAxis = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdSrfStruct *NewSrf = (CagdSrfStruct *) CagdMalloc(sizeof(CagdSrfStruct));

    NewSrf -> GType = Srf -> GType;
    NewSrf -> PType = Srf -> PType;
    NewSrf -> ULength = Srf -> ULength;
    NewSrf -> VLength = Srf -> VLength;
    NewSrf -> UOrder = Srf -> UOrder;
    NewSrf -> VOrder = Srf -> VOrder;
    if (Srf -> UKnotVector != NULL)
	NewSrf -> UKnotVector = BspKnotCopy(Srf -> UKnotVector,
					    Srf -> ULength + Srf -> UOrder);
    else
	NewSrf -> UKnotVector = NULL;
    if (Srf -> VKnotVector != NULL)
	NewSrf -> VKnotVector = BspKnotCopy(Srf -> VKnotVector,
					    Srf -> VLength + Srf -> VOrder);
    else
	NewSrf -> VKnotVector = NULL;
    NewSrf -> Pnext = NULL;
    NewSrf -> Points[0] = NULL;			    /* The rational element. */

    for (i = !CAGD_IS_RATIONAL_PT(Srf -> PType); i <= MaxAxis; i++) {
	NewSrf -> Points[i] = (CagdRType *) CagdMalloc(sizeof(CagdRType) *
					 Srf -> ULength * Srf -> VLength);
	GEN_COPY(NewSrf -> Points[i], Srf -> Points[i],
		 sizeof(CagdRType) * Srf -> ULength * Srf -> VLength);
    }

    for (i = MaxAxis + 1; i <= CAGD_MAX_PT_COORD; i++)
	NewSrf -> Points[i] = NULL;

    return NewSrf;
}

/******************************************************************************
* Allocate and copy all slots of a polygon structure:			      *
******************************************************************************/
CagdPolygonStruct *CagdPolygonCopy(CagdPolygonStruct *Poly)
{
    CagdPolygonStruct *NewPoly = (CagdPolygonStruct *)
    					CagdMalloc(sizeof(CagdPolygonStruct));

    GEN_COPY(NewPoly, Poly, sizeof(CagdPolygonStruct));
    NewPoly -> Pnext = NULL;

    return NewPoly;
}


/******************************************************************************
* Allocate and copy all slots of a polyline structure:			      *
******************************************************************************/
CagdPolylineStruct *CagdPolylineCopy(CagdPolylineStruct *Poly)
{
    CagdPolylineStruct *NewPoly = (CagdPolylineStruct *)
    					CagdMalloc(sizeof(CagdPolylineStruct));

    NewPoly -> Polyline = (CagdPtStruct *) CagdMalloc(sizeof(CagdPtStruct) *
							Poly -> Length);
    GEN_COPY(NewPoly -> Polyline, Poly -> Polyline,
					sizeof(CagdPtStruct) * Poly -> Length);
    NewPoly -> Pnext = NULL;
    NewPoly -> Length = Poly -> Length;

    return NewPoly;
}

/******************************************************************************
* Deallocate and free all slots of a curve structure:			      *
******************************************************************************/
void CagdCrvFree(CagdCrvStruct *Crv)
{
    int i, MaxAxis = CAGD_NUM_OF_PT_COORD(Crv -> PType);

    for (i = !CAGD_IS_RATIONAL_SRF(Crv); i <= MaxAxis; i++)
	CagdFree((VoidPtr) Crv -> Points[i]);

    if (Crv -> KnotVector != NULL)
	CagdFree((VoidPtr) Crv -> KnotVector);

    CagdFree((VoidPtr) Crv);
}

/******************************************************************************
* Deallocate and free a curve structure list:				      *
******************************************************************************/
void CagdCrvFreeList(CagdCrvStruct *CrvList)
{
    CagdCrvStruct *CrvTemp;

    while (CrvList) {
	CrvTemp = CrvList -> Pnext;
	CagdCrvFree(CrvList);
	CrvList = CrvTemp;
    }
}

/******************************************************************************
* Deallocate and free all slots of a surface structure:			      *
******************************************************************************/
void CagdSrfFree(CagdSrfStruct *Srf)
{
    int i, MaxAxis = CAGD_NUM_OF_PT_COORD(Srf -> PType);

    for (i = !CAGD_IS_RATIONAL_SRF(Srf); i <= MaxAxis; i++)
	CagdFree((VoidPtr) Srf -> Points[i]);

    if (Srf -> UKnotVector != NULL)
	CagdFree((VoidPtr) Srf -> UKnotVector);
    if (Srf -> VKnotVector != NULL)
	CagdFree((VoidPtr) Srf -> VKnotVector);

    CagdFree((VoidPtr) Srf);
}

/******************************************************************************
* Deallocate and free a curve structure list:				      *
******************************************************************************/
void CagdSrfFreeList(CagdSrfStruct *SrfList)
{
    CagdSrfStruct *SrfTemp;

    while (SrfList) {
	SrfTemp = SrfList -> Pnext;
	CagdSrfFree(SrfList);
	SrfList = SrfTemp;
    }
}

/******************************************************************************
* Deallocate and free a curve structure list:				      *
******************************************************************************/
void CagdPolylineFree(CagdPolylineStruct *Poly)
{
    CagdFree((VoidPtr) Poly -> Polyline);
    CagdFree((VoidPtr) Poly);
}

/******************************************************************************
* Deallocate and free a curve structure list:				      *
******************************************************************************/
void CagdPolylineFreeList(CagdPolylineStruct *PolyList)
{
    CagdPolylineStruct *PolyTemp;

    while (PolyList) {
	PolyTemp = PolyList -> Pnext;
	CagdPolylineFree(PolyList);
	PolyList = PolyTemp;
    }
}

/******************************************************************************
* Deallocate and free a curve structure list:				      *
******************************************************************************/
void CagdPolygonFree(CagdPolygonStruct *Poly)
{
    CagdFree((VoidPtr) Poly);
}

/******************************************************************************
* Deallocate and free a curve structure list:				      *
******************************************************************************/
void CagdPolygonFreeList(CagdPolygonStruct *PolyList)
{
    CagdPolygonStruct *PolyTemp;

    while (PolyList) {
	PolyTemp = PolyList -> Pnext;
	CagdPolygonFree(PolyList);
	PolyList = PolyTemp;
    }
}

#if defined(__DEBUG_MALLOC__) || defined(__DEBUG_FREE_BOGUS__)
/*****************************************************************************
* My Routine to	allocate dynamic memory. All program requests must call this *
* routine (no direct call to malloc). Dies if no memory.		     *
*****************************************************************************/
static void AllocError(const char *Msg, VoidPtr *p)
{
    char s[128];

    /* Note s is practically ignored, but you may want to print it. */
#ifdef __MSDOS__
    sprintf(s, "%s, Ptr = %p\n", Msg, p);
#else
    sprintf(s, "%s, Ptr = %lx\n", Msg, p);
#endif /* __MSDOS__ */

    FATAL_ERROR(CAGD_ERR_ALLOC_ERR);
}
#endif /* __DEBUG_MALLOC__ || __DEBUG_FREE_BOGUS__ */

/*****************************************************************************
* My Routine to	allocate dynamic memory. All program requests must call this *
* routine (no direct call to malloc). Dies if no memory.		     *
*****************************************************************************/
VoidPtr CagdMalloc(unsigned size)
{
    VoidPtr p;

    if ((p = malloc(size)) != NULL) {
#	ifdef __DEBUG_FREE_BOGUS__
	    int i;

	    /* Save the allocated pointer so we can search for it when freed. */
	    for (i = 0; i < AllocNumPtrs; i++)
	        if (AllocPtrTable[i] == NULL) {
		    AllocPtrTable[i] = p;
		    break;
		}

	    if (i >= AllocNumPtrs) AllocPtrTable[AllocNumPtrs++] = p;
#	endif /* __DEBUG_FREE_BOGUS__ */

	return p;
    }

    CagdFatalError(CAGD_ERR_NOT_ENOUGH_MEM);

    return NULL;				    /* Make warnings silent. */
}

/*****************************************************************************
* My Routine to	free dynamic memory. All program requests must call this     *
* routine (no direct call to free).					     *
*****************************************************************************/
void CagdFree(VoidPtr p)
{
#ifdef __DEBUG_FREE_BOGUS__
    int i;

    /* Compare the freed pointer with the list of allocated ones. */
    for (i = 0; i < AllocNumPtrs; i++) if (AllocPtrTable[i] == p) break;
    if (i >= AllocNumPtrs) AllocError("Free unallocated pointer", p);
#endif /* __DEBUG_FREE_BOGUS__ */

#ifdef __DEBUG_MALLOC__
    switch (heapcheck()) {
	case _HEAPCORRUPT:
	    AllocError("Heap is corrupted", p);
	    break;
	case _BADNODE:
	    AllocError("Attempt to free a bogus pointer", p);
	    break;
	case _FREEENTRY:
	    AllocError("Attempt to free an already freed pointer", p);
	    break;
	case _HEAPOK:
	case _HEAPEMPTY:
	case _USEDENTRY:
	    break;
	default:
	    AllocError("Allocation error", p);
	    break;

    }
#endif /* __DEBUG_MALLOC__ */

    free(p);
}

/******************************************************************************
* Linear transform in place given Crv as specified by Translate and Scale.    *
******************************************************************************/
void CagdCrvTransform(CagdCrvStruct *Crv, CagdRType Translate[3],
							     CagdRType Scale)
{
    switch (Crv -> GType) {
	case CAGD_CBEZIER_TYPE:
	case CAGD_CBSPLINE_TYPE:
    	    CagdTransform(Crv -> Points,
			  Crv -> Length,
    			  CAGD_NUM_OF_PT_COORD(Crv -> PType),
    			  !CAGD_IS_RATIONAL_CRV(Crv),
			  Translate,
			  Scale);

	    break;
	case CAGD_CPOWER_TYPE:
	    FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
	    break;
    }
}

/******************************************************************************
* Linear transform in place given Srf as specified by Translate and Scale.    *
******************************************************************************/
void CagdSrfTransform(CagdSrfStruct *Srf, CagdRType Translate[3],
							     CagdRType Scale)
{
    switch (Srf -> GType) {
	case CAGD_SBEZIER_TYPE:
	case CAGD_SBSPLINE_TYPE:
	    CagdTransform(Srf -> Points,
	    		  Srf -> ULength * Srf -> VLength,
	                  CAGD_NUM_OF_PT_COORD(Srf -> PType),
			  !CAGD_IS_RATIONAL_SRF(Srf),
		          Translate,
        	          Scale);
	    break;
	case CAGD_SPOWER_TYPE:
	    FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_SRF);
	    break;
    }
}

/******************************************************************************
* Linear transform in place given Data as specified by Translate and Scale.   *
******************************************************************************/
static void CagdTransform(CagdRType **Points, int Len, int MaxCoord,
	     CagdBType IsNotRational, CagdRType Translate[3], CagdRType Scale)
{
    int i, j;

    if (IsNotRational)
	for (i = 1; i <= MaxCoord; i++)
	    for (j = 0; j < Len; j++)
		Points[i][j] = (Points[i][j] + Translate[i - 1]) * Scale;
    else
	for (i = 1; i <= MaxCoord; i++)
	    for (j = 0; j < Len; j++)
		Points[i][j] = (Points[i][j] +
				Translate[i - 1] * Points[W][j]) * Scale;
}

/******************************************************************************
* Linear transform in place the given Crv as specified by Mat.		      *
******************************************************************************/
void CagdCrvMatTransform(CagdCrvStruct *Crv, CagdMatStruct *Mat)
{
    switch (Crv -> GType) {
	case CAGD_CBEZIER_TYPE:
	case CAGD_CBSPLINE_TYPE:
	    CagdMatTransform(Crv -> Points,
    			     Crv -> Length,
                	     CAGD_NUM_OF_PT_COORD(Crv -> PType),
			     !CAGD_IS_RATIONAL_CRV(Crv),
		             Mat);
	    break;
	case CAGD_CPOWER_TYPE:
	    FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_CRV);
	    break;
    }
}

/******************************************************************************
* Linear transform in place the given Srf as specified by Mat.		      *
******************************************************************************/
void CagdSrfMatTransform(CagdSrfStruct *Srf, CagdMatStruct *Mat)
{
    switch (Srf -> GType) {
	case CAGD_SBEZIER_TYPE:
	case CAGD_SBSPLINE_TYPE:
	    CagdMatTransform(Srf -> Points,
    			     Srf -> ULength * Srf -> VLength,
        	             CAGD_NUM_OF_PT_COORD(Srf -> PType),
			     !CAGD_IS_RATIONAL_SRF(Srf),
		             Mat);
	    break;
	case CAGD_SPOWER_TYPE:
	    FATAL_ERROR(CAGD_ERR_POWER_NO_SUPPORT);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_UNDEF_SRF);
	    break;
    }
}

/******************************************************************************
* Linear transform in place the given Data as specified by Mat.		      *
******************************************************************************/
static void CagdMatTransform(CagdRType **Points, int Len, int MaxCoord,
				  CagdBType IsNotRational, CagdMatStruct *Mat)
{
    int i, j;
    CagdRType P[4], Q[4];

    if (IsNotRational)
	for (i = 0; i < Len; i++) {
	    for (j = 1; j <= MaxCoord; j++) P[j - 1] = Points[j][i];
	    P[3] = 1.0;

	    /* Zero unused coords. */
	    for (j = MaxCoord + 1; j < 3; i++) P[j - 1] = 0.0;

            CagdMultMatVec(Mat, P, Q);

	    for (j = 1; j <= MaxCoord; j++) Points[j][i] = Q[j - 1];
        }
    else
	for (i = 0; i < Len; i++) {
	    for (j = 1; j <= MaxCoord; j++) P[j - 1] = Points[j][i];
            P[3] = Points[W][i];

	    /* Zero unused coords. */
	    for (j = MaxCoord + 1; j < 3; i++) P[j - 1] = 0.0;

            CagdMultMatVec(Mat, P, Q);

	    for (j = 1; j <= MaxCoord; j++) Points[j][i] = Q[j - 1];
	    Points[W][i] = Q[3];
        }
}

/*****************************************************************************
* Routine to create one polygon given its vertices.			     *
*****************************************************************************/
CagdPolygonStruct *_CagdMakePolygon(CagdBType ComputeNormals,
				    CagdPtStruct *Pt1,
				    CagdPtStruct *Pt2,
				    CagdPtStruct *Pt3,
				    CagdVecStruct *Nl1,
				    CagdVecStruct *Nl2,
				    CagdVecStruct *Nl3)
{
    CagdPolygonStruct *Poly = CagdPolygonNew();

    CAGD_COPY_POINT(Poly -> Polygon[0], *Pt1);
    CAGD_COPY_POINT(Poly -> Polygon[1], *Pt2);
    CAGD_COPY_POINT(Poly -> Polygon[2], *Pt3);

    if (ComputeNormals) {
	CAGD_COPY_VECTOR(Poly -> Normal[0], *Nl1);
	CAGD_COPY_VECTOR(Poly -> Normal[1], *Nl2);
	CAGD_COPY_VECTOR(Poly -> Normal[2], *Nl3);
    }
    else {
	CAGD_RESET_VECTOR(Poly -> Normal[0]);
	CAGD_RESET_VECTOR(Poly -> Normal[1]);
	CAGD_RESET_VECTOR(Poly -> Normal[2]);
    }

    return Poly;
}

/*****************************************************************************
* Set the way (co)linear surfaces are converted into polygons.		     *
*****************************************************************************/
void CagdSetLinear2Poly(CagdLin2PolyType Lin2Poly)
{
    _CagdLin2Poly = Lin2Poly;
}