/****************************************************************************** * CagdOslo.c - an implementation of the oslo algorithm (1). * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Aug. 92. * ******************************************************************************/ #include "cagd_loc.h" #define OSLO_IRIT_EPS 1e-20 #define OSLO_APX_EQ(x, y) (ABS((x) - (y)) < OSLO_IRIT_EPS) #ifdef DEBUG_OSLO static void PrintAlphaMat(BspKnotAlphaCoeffType *A); #endif /* DEBUG_OSLO */ /***************************************************************************** * DESCRIPTION: M * Computes the values of the alpha coefficients, Ai,k(j) of order k: M * M * n n m m n V * _ _ _ _ _ V * C(t) = \ Pi Bi,k(t) = \ Pi \ Ai,k(j) Nj,k(t) = \ ( \ Pi Ai,k(j) ) Nj,k(t)V * / / / / / V * - - - - - V * i=0 i=0 j=0 j=0 i=0 V * M * Let T be the original knot vector and let t be the refined one, i.e. T is M * a subset of t. M * The Ai,k(j) are computed from the following recursive definition: M * M * 1, T(i) <= t(i) < T(i+1) V * / V * Ai,1(j) = V * \ V * 0, otherwise. V * V * V * T(j+k-1) - T(i) T(i+k) - T(j+k-1) V * Ai,k(j) = --------------- Ai,k-1(j) + --------------- Ai+1,k-1(j) V * T(i+k-1) - T(i) T(i+k) - T(i+1) V * M * LengthKVT + k is the length of KVT and similarly LengthKVt + k is the M * length of KVt. In other words, LengthKVT and LengthKVt are the control M * points len... M * M * The output matrix has LengthKVT rows and LengthKVt columns (#cols > #rows) M * ColIndex/Length hold LengthKVt pairs of first non zero scalar and length ofM * non zero values in that column, so not all LengthKVT scalars are blended. M * * * PARAMETERS: M * k: Order of geometry. M * KVT: Original knot vector. M * LengthKVT: Length of original control polygon with KVT knot vector. M * KVt: Refined knot vector. Must contain all knots of KVT. M * LengthKVt: Length of refined control polygon with KVt knot vector. M * * * RETURN VALUE: M * BspKnotAlphaCoeffType *: A matrix to multiply the coefficients of the M * geometry using KVT, in order to get the M * coefficients under the space defined using M * KVt that represent the same geometry. M * * * SEE ALSO: M * BspKnotFreeAlphaCoef, BspKnotEvalAlphaCoefMerge, BspCrvKnotInsert, M * BspSrfKnotInsert M * * * KEYWORDS: M * BspKnotEvalAlphaCoef, alpha matrix, refinement M *****************************************************************************/ BspKnotAlphaCoeffType *BspKnotEvalAlphaCoef(int k, CagdRType *KVT, int LengthKVT, CagdRType *KVt, int LengthKVt) { int i, j, o; CagdRType *m, **Rows; BspKnotAlphaCoeffType *A = (BspKnotAlphaCoeffType *) IritMalloc(sizeof(BspKnotAlphaCoeffType)); A -> Order = k; A -> Length = LengthKVT; A -> RefLength = LengthKVt; A -> Matrix = (CagdRType *) IritMalloc(sizeof(CagdRType) * (LengthKVT + 1) * LengthKVt); Rows = A -> Rows = (CagdRType **) IritMalloc(sizeof(CagdRType *) * (LengthKVT + 1)); A -> ColIndex = (int *) IritMalloc(sizeof(int) * LengthKVt); A -> ColLength = (int *) IritMalloc(sizeof(int) * LengthKVt); /* Update the row pointers to point onto the matrix rows. */ for (i = 0, j = 0; i <= LengthKVT; i++, j += LengthKVt) Rows[i] = &A -> Matrix[j]; /* Initialize the matrix with according to order 1: */ for (m = A -> Matrix, i = (LengthKVT + 1) * LengthKVt; i > 0; i--) *m++ = 0.0; for (i = 0; i < LengthKVT; i++) { CagdRType *RowI = Rows[i], KVTI = KVT[i], KVTI1 = KVT[i + 1]; for (j = 0; j < LengthKVt; j++, RowI++) { if (KVTI <= KVt[j] && KVt[j] < KVTI1) *RowI = 1.0; } } /* Iterate upto order k: */ for (o = 2; o <= k; o++) { for (i = 0; i < LengthKVT; i++) { CagdRType *RowI = Rows[i], *RowI1 = Rows[i + 1], KVTI = KVT[i], KVTIO = KVT[i + o], t1 = KVT[i + o - 1] - KVTI, t2 = KVTIO - KVT[i + 1]; /* If ti == 0, the whole term should be ignored. */ t1 = t1 < OSLO_IRIT_EPS ? 0.0 : 1.0 / t1; t2 = t2 < OSLO_IRIT_EPS ? 0.0 : 1.0 / t2; for (j = 0; j < LengthKVt; j++, RowI++, RowI1++) { int jo1 = j + o - 1; *RowI = *RowI * (KVt[jo1] - KVTI) * t1 + *RowI1 * (KVTIO - KVt[jo1]) * t2; } } } /* Update the row non zero indices. */ for (i = LengthKVt - 1; i >= 0; i--) { for (j = 0; j < LengthKVT && OSLO_APX_EQ(Rows[j][i], 0.0); j++); A -> ColIndex[i] = j; for (j = LengthKVT - 1; j >= 0 && OSLO_APX_EQ(Rows[j][i], 0.0); j--); if ((A -> ColLength[i] = j + 1 - A -> ColIndex[i]) <= 0) CAGD_FATAL_ERROR(CAGD_ERR_DEGEN_ALPHA); } return A; } /***************************************************************************** * DESCRIPTION: M * Frees the BspKnotAlphaCoeffType data structrure. M * * * PARAMETERS: M * A: Alpha matrix to free. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * BspKnotEvalAlphaCoef, BspKnotEvalAlphaCoefMerge, BspCrvKnotInsert, M * BspSrfKnotInsert M * * * KEYWORDS: M * BspKnotFreeAlphaCoef, alpha matrix, refinement M *****************************************************************************/ void BspKnotFreeAlphaCoef(BspKnotAlphaCoeffType *A) { IritFree((VoidPtr) A -> Matrix); IritFree((VoidPtr) A -> Rows); IritFree((VoidPtr) A -> ColIndex); IritFree((VoidPtr) A -> ColLength); IritFree((VoidPtr) A); } /***************************************************************************** * DESCRIPTION: M * Same as EvalAlphaCoef but the new knot set NewKV is merged with KVT to M * form the new knot vector KVt. M * * * PARAMETERS: M * k: Order of geometry. M * KVT: Original knot vector. M * LengthKVT: Length of original knot vector. M * NewKV: A sequence of new knots to introduce into KVT. M * LengthNewKV: Length of new knot sequence. M * * * RETURN VALUE: M * BspKnotAlphaCoeffType *: A matrix to multiply the coefficients of the M * geometry using KVT, in order to get the M * coefficients under the space defined using M * KVt that represent the same geometry. M * * * SEE ALSO: M * BspKnotFreeAlphaCoef, BspKnotEvalAlphaCoef, BspCrvKnotInsert, M * BspSrfKnotInsert M * * * KEYWORDS: M * BspKnotEvalAlphaCoefMerge, alpha matrix, refinement M *****************************************************************************/ BspKnotAlphaCoeffType *BspKnotEvalAlphaCoefMerge(int k, CagdRType *KVT, int LengthKVT, CagdRType *NewKV, int LengthNewKV) { BspKnotAlphaCoeffType *A; if (NewKV == NULL || LengthNewKV == 0) { A = BspKnotEvalAlphaCoef(k, KVT, LengthKVT, KVT, LengthKVT); } else { int LengthKVt; CagdRType *KVt = BspKnotMergeTwo(KVT, LengthKVT + k, NewKV, LengthNewKV, 0, &LengthKVt); A = BspKnotEvalAlphaCoef(k, KVT, LengthKVT, KVt, LengthKVt - k); IritFree(KVt); } return A; } /***************************************************************************** * DESCRIPTION: M * Prepares a refinement vector for the given knot vector domain with n M * inserted knots equally spaced. M * * * PARAMETERS: M * n: Number of knots to introduce. M * Tmin: Minimum domain to introduce knots. M * Tmax: Maximum domain to introduce knots. M * * * RETURN VALUE: M * CagdRType *: A vector of n knots uniformly spaced between TMin and TMax. M * * * KEYWORDS: M * BspKnotPrepEquallySpaced, refinement M *****************************************************************************/ CagdRType *BspKnotPrepEquallySpaced(int n, CagdRType Tmin, CagdRType Tmax) { int i; CagdRType dt, t, *RefKV; if (n <= 0) { CAGD_FATAL_ERROR(CAGD_ERR_WRONG_INDEX); return NULL; } dt = (Tmax - Tmin) / (n + 1), t = Tmin + dt, RefKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * n); for (i = 0; i < n; i++, t += dt) { RefKV[i] = t; } return RefKV; } #ifdef DEBUG_OSLO /***************************************************************************** * DESCRIPTION: * * Prints the content of the alpha matrix. * * * * PARAMETERS: * * A: Alpha matrix to print. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void PrintAlphaMat(BspKnotAlphaCoeffType *A) { int i, j; fprintf(stderr, "Order = %d, Length = %d\n", A -> Order, A -> Length); for (i = 0; i < A -> Length; i++) { for (j = 0; j < A -> RefLength; j++) fprintf(stderr, " %9.5lg", A -> Rows[i][j]); fprintf(stderr, "\n"); } fprintf(stderr, " "); for (j = 0; j < A -> RefLength; j++) fprintf(stderr, " %3d %3d |", A -> ColIndex[j], A -> ColLength[j]); fprintf(stderr, "\n"); } #endif /* DEBUG_OSLO */