/****************************************************************************** * SBsp-Int.c - Bspline surface interpolation/approximation schemes. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Feb. 94. * ******************************************************************************/ #include #include #include #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; }