/***************************************************************************** * A matching code between two freeform curves. * ****************************************************************************** * (C) Gershon Elber, Technion, Israel Institute of Technology * ****************************************************************************** * Written by: Shmuel Cohen Ver 0.1, May 1995. * *****************************************************************************/ #include "cagd_loc.h" #include "geomat3d.h" #define SAMPLE_NUMBER 500 #define BISECTOR_MIN_DOT_PROD 0.001 typedef struct LinkStruct { int x, y; } LinkStruct; static int GlblShiftFlag = 1; static int GlblNumberOfSamples; static int GlblRotationFlag = FALSE; static int GlblAllowNegativeNorm = FALSE; static int GlblIgnoreUpperBoundLimit = FALSE; static CagdBType FindTangentsVec(CagdCrvStruct *Crv, CagdVType TangentsVector[SAMPLE_NUMBER], CagdVType CrvEval[SAMPLE_NUMBER]); static void MatchingTest(CagdVType TangentsVectorsCrv1[SAMPLE_NUMBER], CagdVType TangentsVectorsCrv2[SAMPLE_NUMBER], CagdRType TangentsDotProdRes[SAMPLE_NUMBER], int NewIndex1[2 * SAMPLE_NUMBER], int NewIndex2[2 * SAMPLE_NUMBER]); static void CopyTwoArrays(int *Dest1, int *Dest2, int *Src1, int *Src2); static void RotateVector(CagdVType SrcVec, CagdVType DestVec, CagdRType Angle); static CagdRType AngleBetweenVectors(CagdVType Vector1, CagdVType Vector2); static void FindWayInMat(int i, int j, int *Cnt, LinkStruct MatchMat[SAMPLE_NUMBER][SAMPLE_NUMBER], CagdRType MaxMat[SAMPLE_NUMBER + 1][SAMPLE_NUMBER + 1], int NewIndex1[2 * SAMPLE_NUMBER], int NewIndex2[2 * SAMPLE_NUMBER]); static CagdBType IfCrvIsClose(CagdCrvStruct *Crv); static CagdRType FindNewMatching(CagdVType TanVecCrv1[SAMPLE_NUMBER], CagdVType TanVecCrv2[SAMPLE_NUMBER], CagdVType Crv1Eval[SAMPLE_NUMBER], CagdVType Crv2Eval[SAMPLE_NUMBER], int NewIndex1[2 * SAMPLE_NUMBER], int NewIndex2[2 * SAMPLE_NUMBER], CagdMatchNormFuncType MatchNormFunc); static int FindBestMatching(CagdVType TanVecCrv1[SAMPLE_NUMBER], CagdVType TanVecCrv2[SAMPLE_NUMBER], CagdVType Crv1Eval[SAMPLE_NUMBER], CagdVType Crv2Eval[SAMPLE_NUMBER], int NewIndex1[2 * SAMPLE_NUMBER], int NewIndex2[2 * SAMPLE_NUMBER], CagdMatchNormFuncType MatchNormFunc); static CagdCrvStruct *CreateReparamCrv(int *Vect, int *Vecr, int Order, int Len, int Reduce); static void FixVector(int *OldVec, CagdRType *NewVec, int Len); static void FixCrv(CagdCrvStruct *Crv); static void PolyTransform(CagdRType **Poly, int Len, CagdRType NewBegin, CagdRType NewEnd); static void VectorTransform(CagdRType *Vec, CagdRType NewBegin, CagdRType NewEnd, int Len); /***************************************************************************** * DESCRIPTION: * * Finds the tangent vector to the curve at each parameter. * * * * PARAMETERS: * * Crv: The input curve. * * TangentsVector: An array for the tangent at each parameter. * * CrvEval: The value of the curve at each parameter. * * * * RETURN VALUE: * * CagdBType: TRUE for Sucssess. * *****************************************************************************/ static CagdBType FindTangentsVec(CagdCrvStruct *Crv, CagdVType TangentsVector[SAMPLE_NUMBER], CagdVType CrvEval[SAMPLE_NUMBER]) { int i; CagdRType StartCrvPar, EndCrvPar, Par, *Pt, DPar; CagdCrvStruct *TanCrv = CagdCrvDerive(Crv); CagdCrvDomain(Crv, &StartCrvPar, &EndCrvPar); DPar = (EndCrvPar - StartCrvPar) / (GlblNumberOfSamples - 1); for (i = 0, Par = StartCrvPar; i < GlblNumberOfSamples; i++, Par += DPar) { if (Par > EndCrvPar) /* Due to floating point roundoff errors. */ Par = EndCrvPar; Pt = CagdCrvEval(TanCrv, Par); CagdCoerceToE3(TangentsVector[i], &Pt, -1, TanCrv -> PType); PT_NORMALIZE(TangentsVector[i]); Pt = CagdCrvEval(Crv , Par); CagdCoerceToE3(CrvEval[i], &Pt, -1, TanCrv -> PType); } CagdCrvFree(TanCrv); return TRUE; } /***************************************************************************** * DESCRIPTION: * * Checks the matching with the tangents' inner product norm. * * * * PARAMETERS: * * TangentsVectorsCrv1: An array for the first curve's tangents. * * TangentsVectorsCrv2: An array for the second curve's tangents. * * TangentsDotProdRes: The inner product resulting vector. * * NewIndex1: The vector of the new matching to the first curve. * * NewIndex1: The vector of the new matching to the second curve.* * * * RETURN VALUE: * * void * *****************************************************************************/ static void MatchingTest(CagdVType TangentsVectorsCrv1[SAMPLE_NUMBER], CagdVType TangentsVectorsCrv2[SAMPLE_NUMBER], CagdRType TangentsDotProdRes[SAMPLE_NUMBER], int NewIndex1[2 * SAMPLE_NUMBER], int NewIndex2[2 * SAMPLE_NUMBER]) { int i; CagdRType Sum; Sum = 0; for (i = 0; i < 2 * GlblNumberOfSamples && NewIndex1[i] != -1; i++) { TangentsDotProdRes[i] = GMVecDotProd(TangentsVectorsCrv1[NewIndex1[i]], TangentsVectorsCrv2[NewIndex1[i]]); Sum += TangentsDotProdRes[i]; #ifdef DEBUG_MATCH_2CRV if (TangentsDotProdRes[i] <= 0) fprintf(stderr, "NO MATCHING at sample point number: %d with %d\n", NewIndex1[i], NewIndex2[i]); #endif /* DEBUG_MATCH_2CRV */ } #ifdef DEBUG_MATCH_2CRV fprintf(stderr, "\nThe sum after testing matching is: %lf\n", Sum); #endif /* DEBUG_MATCH_2CRV */ } /***************************************************************************** * DESCRIPTION: * * Copies two integer vectors. * * * * PARAMETERS: * * Dest1: The first destination vector. * * Dest2: The second destination vector. * * Src1: The first source vector. * * Src2: The second source vector. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void CopyTwoArrays(int *Dest1, int *Dest2, int *Src1, int *Src2) { GEN_COPY(Dest1, Src1, sizeof(int) * 2 * SAMPLE_NUMBER); GEN_COPY(Dest2, Src2, sizeof(int) * 2 * SAMPLE_NUMBER); } /***************************************************************************** * DESCRIPTION: * * find the angle between two vectors. * * * * PARAMETERS: * * Vector1: The first input vector. * * vector2: The second input vector. * * * * RETURN VALUE: * * CagdRType: The angle between the two inputvectors * *****************************************************************************/ static CagdRType AngleBetweenVectors(CagdVType Vector1, CagdVType Vector2) { CagdRType Cos = DOT_PROD(Vector1, Vector2); return acos(Cos); } /***************************************************************************** * DESCRIPTION: * * rotate vector in input angle. * * * * PARAMETERS: * * SrcVec: The input vector. * * DestVec: The vector after the rotation. * * Angle: the rotation angle. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void RotateVector(CagdVType SrcVec, CagdVType DestVec, CagdRType Angle) { DestVec[0] = SrcVec[0] * cos(Angle) - SrcVec[1] * sin(Angle); DestVec[1] = SrcVec[0] * sin(Angle) + SrcVec[1] * cos(Angle); DestVec[2] = SrcVec[2]; } /***************************************************************************** * DESCRIPTION: * * Finds the path in the cost matrix (MaxMat). * * * * PARAMETERS: * * i: Line. * * j: Column. * * Cnt: Place in NewIndex1 and NewIndex2. * * MatchMat: The path matrix. * * MaxMat: Contains the costs. * * NewIndex1: The vector of the new matching to the first curve. * * NewIndex2: The vector of the new matching to the second curve * * * * RETURN VALUE: * * void * *****************************************************************************/ static void FindWayInMat(int i, int j, int *Cnt, LinkStruct MatchMat[SAMPLE_NUMBER][SAMPLE_NUMBER], CagdRType MaxMat[SAMPLE_NUMBER + 1][SAMPLE_NUMBER + 1], int NewIndex1[2 * SAMPLE_NUMBER], int NewIndex2[2 * SAMPLE_NUMBER]) { int k, l; if (i != 1 || j != 1) { k = MatchMat[i - 1][j - 1].x; l = MatchMat[i - 1][j - 1].y; FindWayInMat(k, l, Cnt, MatchMat, MaxMat, NewIndex1, NewIndex2); } NewIndex1[*Cnt] = i - 1; NewIndex2[(*Cnt)++] = j - 1; NewIndex1[*Cnt] = NewIndex2[*Cnt] = -1; } /***************************************************************************** * DESCRIPTION: * * Checks if the input curve is closed. * * * * PARAMETERS: * * Crv: The input curve to check. * * * * RETURN VALUE: * * CagdBType: TRUE - closed curve, FALSE - not a closed curve. * *****************************************************************************/ static CagdBType IfCrvIsClose(CagdCrvStruct *Crv) { CagdRType *Pt, StartPara, EndPara; CagdPType Begin, End; CagdCrvDomain(Crv, &StartPara, &EndPara); Pt = CagdCrvEval(Crv, StartPara); CagdCoerceToE3(Begin, &Pt, -1, Crv -> PType); Pt = CagdCrvEval(Crv, EndPara); CagdCoerceToE3(End, &Pt, -1, Crv -> PType); return PT_PT_DIST(Begin, End) < IRIT_EPS; } /***************************************************************************** * DESCRIPTION: M * Computes the distance norm to the matching. M * * * PARAMETERS: M * T1: A pointer to tangent to the first curve at i-th point. M * T2: A pointer to tangent to the second curve at j-th point. M * P1: A pointer to value of the first curve at i-th point. M * P2: A pointer to value of the second curve at j-th point. M * * * RETURN VALUE: M * CagdRType: A numeric matching value, the smaller the better. M * * * KEYWORDS: M * CagdMatchDistNorm, correspondance, matching M *****************************************************************************/ CagdRType CagdMatchDistNorm(CagdVType T1, CagdVType T2, CagdVType P1, CagdVType P2) { CagdVType V; GlblIgnoreUpperBoundLimit = TRUE; PT_SUB(V, P1, P2); return PT_LENGTH(V); } /***************************************************************************** * DESCRIPTION: M * Computes the distance norm to the matching. M * * * PARAMETERS: M * T1: A pointer to tangent to the first curve at i-th point. M * T2: A pointer to tangent to the second curve at j-th point. M * P1: A pointer to value of the first curve at i-th point. M * P2: A pointer to value of the second curve at j-th point. M * * * RETURN VALUE: M * CagdRType: A numeric matching value, the smaller the better. M * * * KEYWORDS: M * CagdMatchDistNorm, correspondance, matching M *****************************************************************************/ CagdRType CagdMatchBisectorNorm(CagdVType T1, CagdVType T2, CagdVType P1, CagdVType P2) { CagdVType N1, N2; CagdPType ClosestPt1, ClosestPt2; CagdRType Param1, Param2; N1[0] = T1[1]; N1[1] = -T1[0]; N2[0] = -T2[1]; N2[1] = T2[0]; N2[2] = N1[2] = 0.0; /* If normals are colinear, measure the distance between the one */ /* normal ray to the other point. */ if (fabs(DOT_PROD(T1, N2)) < BISECTOR_MIN_DOT_PROD) return CGDistPointLine(P1, P2, N2) + CGDistPointLine(P2, P1, N1); else { if (CG2PointsFromLineLine(P1, N1, P2, N2, ClosestPt1, &Param1, ClosestPt2, &Param2)) { return fabs(Param1 - Param2); } else return CGDistPointLine(P1, P2, N2) + CGDistPointLine(P2, P1, N1); } } /***************************************************************************** * DESCRIPTION: M * Computes the default morphing norm to the matching. M * * * PARAMETERS: M * T1: A pointer to tangent to the first curve at i-th point. M * T2: A pointer to tangent to the second curve at j-th point. M * P1: A pointer to value of the first curve at i-th point. M * P2: A pointer to value of the second curve at j-th point. M * * * RETURN VALUE: M * CagdRType: -1 for no matching or the cost of the matching for the point M * between zero and one. M * * * KEYWORDS: M * CagdMatchMorphNorm, correspondance, matching M *****************************************************************************/ CagdRType CagdMatchMorphNorm(CagdVType T1, CagdVType T2, CagdVType P1, CagdVType P2) { CagdRType MultRes = DOT_PROD(T1, T2); GlblIgnoreUpperBoundLimit = FALSE; if (!GlblAllowNegativeNorm && MultRes < 0) return -1.0; return 1.0 - MultRes; } /***************************************************************************** * DESCRIPTION: M * Computes the default ruled norm to the matching. M * * * PARAMETERS: M * T1: A pointer to tangent to the first curve at i-th point. M * T2: A pointer to tangent to the second curve at j-th point. M * P1: A pointer to value of the first curve at i-th point. M * P2: A pointer to value of the second curve at j-th point. M * * * RETURN VALUE: M * CagdRType: -1 for no matching or the cost of the matching for the point M * between zero and one. M * * * KEYWORDS: M * CagdMatchRuledNorm, correspondance, matching M *****************************************************************************/ CagdRType CagdMatchRuledNorm(CagdVType T1, CagdVType T2, CagdVType P1, CagdVType P2) { CagdPType RuledVec, CrossVec1, CrossVec2; CagdRType MultRes; GlblIgnoreUpperBoundLimit = FALSE; MultRes = DOT_PROD(T1, T2); PT_SUB(RuledVec, P1, P2); CROSS_PROD(CrossVec1, T1, RuledVec); CROSS_PROD(CrossVec2, T2, RuledVec); if (!GlblAllowNegativeNorm && DOT_PROD(CrossVec1, CrossVec2) < 0) return -1.0; return 1.0 - MultRes; } /***************************************************************************** * DESCRIPTION: * * Finds matching cost for all samples. * * * * PARAMETERS: * * TanVecCrv1: An array for the first curve tangents. * * TanVecCrv2: An array for the second curve tangents. * * Crv1Eval: An array for the first curve values. * * Crv2Eval: An array for the second curve values. * * NewIndex1: The vector of the new matching to the first curve. * * NewIndex2: The vector of the new matching to the second curve. * * MatchNormFunc: A pointer to the matching norm function. * * * * RETURN VALUE: * * CagdRType: The matching cost for the all sample point. * *****************************************************************************/ static CagdRType FindNewMatching(CagdVType TanVecCrv1[SAMPLE_NUMBER], CagdVType TanVecCrv2[SAMPLE_NUMBER], CagdVType Crv1Eval[SAMPLE_NUMBER], CagdVType Crv2Eval[SAMPLE_NUMBER], int NewIndex1[2 * SAMPLE_NUMBER], int NewIndex2[2 * SAMPLE_NUMBER], CagdMatchNormFuncType MatchNormFunc) { /* The only reason these variables are static is to prevent stack */ /* overlow on small machines like Windows NT or OS2. */ static CagdRType DotProdMat[SAMPLE_NUMBER + 1][SAMPLE_NUMBER + 1], MaxMat[SAMPLE_NUMBER + 1][SAMPLE_NUMBER + 1]; static LinkStruct MatchMat[SAMPLE_NUMBER][SAMPLE_NUMBER]; int i, j, Cnt = 0; CagdRType Sum, Min, Res; for (i = 0; i < GlblNumberOfSamples + 1; i++) DotProdMat[i][0] = DotProdMat[0][i] = MaxMat[i][0] = MaxMat[0][i] = 2 * GlblNumberOfSamples; /* Init. the DotProd matrix. */ for (i = 0; i < GlblNumberOfSamples; i++) { CagdRType *DotProdVec = &DotProdMat[i + 1][1]; for (j = 0; j < GlblNumberOfSamples; j++) { Res = MatchNormFunc(TanVecCrv1[i], TanVecCrv2[j], Crv1Eval[i], Crv2Eval[j]); if (!GlblAllowNegativeNorm && Res == -1) *DotProdVec++ = 2 * GlblNumberOfSamples; else *DotProdVec++ = Res; } } MaxMat[0][0] = 0; for (i = 1; i <= GlblNumberOfSamples; i++) { CagdRType *DotProdVec = &DotProdMat[i][1], *MaxVec1 = &MaxMat[i - 1][1], *MaxVec = &MaxMat[i][1]; LinkStruct *MatchVec = MatchMat[i - 1]; for (j = 1; j <= GlblNumberOfSamples; j++, MatchVec++, MaxVec1++, MaxVec++) { Min = MIN(MaxVec1[-1], MIN(*MaxVec1, MaxVec[-1])); *MaxVec = *DotProdVec++ + Min; if (Min == MaxVec1[-1]) { MatchVec -> x = i - 1; MatchVec -> y = j - 1; } else if (Min == MaxVec[-1]) { MatchVec -> x = i; MatchVec -> y = j - 1; } else { MatchVec -> x = i - 1; MatchVec -> y = j; } } } Sum = MaxMat[GlblNumberOfSamples][GlblNumberOfSamples]; if (GlblIgnoreUpperBoundLimit || Sum < 2 * GlblNumberOfSamples) FindWayInMat(GlblNumberOfSamples, GlblNumberOfSamples ,&Cnt, MatchMat, MaxMat, NewIndex1, NewIndex2); return Sum; } /***************************************************************************** * DESCRIPTION: * * Finds the best cost for the matching (if exist) after checking all * * the samples. * * * * PARAMETERS: * * TanVecCrv1: An array for the first curve tangents. * * TanVecCrv2: An array for the second curve tangents. * * Crv1Eval: An array for the first curve values. * * Crv2Eval: An array for the second curve values. * * NewIndex1: The vector of the new matching to the first curve. * * NewIndex2: The vector of the new matching to the second curve. * * MatchNormFunc: A pointer to the matching norm function. * * * * RETURN VALUE: * * int: -1 for no matching between the curves. * *****************************************************************************/ static int FindBestMatching(CagdVType TanVecCrv1[SAMPLE_NUMBER], CagdVType TanVecCrv2[SAMPLE_NUMBER], CagdVType Crv1Eval[SAMPLE_NUMBER], CagdVType Crv2Eval[SAMPLE_NUMBER], int NewIndex1[2 * SAMPLE_NUMBER], int NewIndex2[2 * SAMPLE_NUMBER], CagdMatchNormFuncType MatchNormFunc) { int i, j, ShiftIndex1[2 * SAMPLE_NUMBER], ShiftIndex2[2 * SAMPLE_NUMBER], RotateIndex1[2 * SAMPLE_NUMBER], RotateIndex2[2 * SAMPLE_NUMBER], Shift1 = -1, Shift2 = -1; CagdVType ShiftTanVec[2 * SAMPLE_NUMBER], RotateTanVec[2 * SAMPLE_NUMBER]; CagdRType TmpSum, Angle, ShiftSum = 2 * GlblNumberOfSamples, RotateSum = 2 * GlblNumberOfSamples; for (i = 0; i < GlblNumberOfSamples; i++) { PT_COPY(ShiftTanVec[i], TanVecCrv2[i]); PT_COPY(ShiftTanVec[i + GlblNumberOfSamples], TanVecCrv2[i]); PT_COPY(RotateTanVec[i], TanVecCrv2[i]); PT_COPY(RotateTanVec[i + GlblNumberOfSamples], TanVecCrv2[i]); } for (i = 0; i < GlblShiftFlag; i++) { for (j = 0; j < GlblNumberOfSamples; j++) NewIndex1[j] = NewIndex2[j] = j; NewIndex1[GlblNumberOfSamples] = NewIndex2[GlblNumberOfSamples] = -1; TmpSum = FindNewMatching(TanVecCrv1, &ShiftTanVec[GlblNumberOfSamples - i], Crv1Eval, Crv2Eval, NewIndex1, NewIndex2, MatchNormFunc); if (TmpSum < ShiftSum) { ShiftSum = TmpSum; Shift1 = i; CopyTwoArrays(ShiftIndex1, ShiftIndex2, NewIndex1, NewIndex2); } /* Using rotation. */ if (GlblRotationFlag) { Angle = AngleBetweenVectors(TanVecCrv1[0], RotateTanVec[GlblNumberOfSamples - i]); for (j = 0; j < GlblNumberOfSamples; j++) RotateVector(ShiftTanVec[j], RotateTanVec[j], -Angle); for (j = 0; j < GlblNumberOfSamples; j++) NewIndex1[j] = NewIndex2[j] = j; NewIndex1[GlblNumberOfSamples] = NewIndex2[GlblNumberOfSamples] = -1; /* Check the matching after rotation. */ TmpSum = FindNewMatching(TanVecCrv1, &RotateTanVec[GlblNumberOfSamples - i], Crv1Eval, Crv2Eval, NewIndex1, NewIndex2, MatchNormFunc); if (TmpSum < RotateSum) { RotateSum = TmpSum; Shift2 = i; CopyTwoArrays(RotateIndex1, RotateIndex2, NewIndex1, NewIndex2); } } } #ifdef DEBUG_MATCH_2CRV fprintf(stderr, "\nShiftSum = :%lf at %d\n", ShiftSum, Shift1); fprintf(stderr, "\nRotateSum = :%lf at %d\n", RotateSum, Shift2); #endif /* DEBUG_MATCH_2CRV */ if (Shift1 == -1 && Shift2 == -1) return -1; if (ShiftSum <= RotateSum) CopyTwoArrays(NewIndex1, NewIndex2, ShiftIndex1, ShiftIndex2); else CopyTwoArrays(NewIndex1, NewIndex2, RotateIndex1, RotateIndex2); return ShiftSum <= RotateSum ? Shift1 : Shift2; } /***************************************************************************** * DESCRIPTION: * * Builds the reparametrization Bspline curve from Vect and Vecr, * * with order of Order and with Reduce degrees of freedom. * * * * PARAMETERS: * * Vect: The new parametrization of the first curve. * * Vecr: The new parametrization of the second curve. * * Order: Of the reparametrization curves. * * Len: The length of Vect and Vecr. * * Reduce: The degree of freedom. * * * * RETURN VALUE: * * CagdCrvStruct *: The reparametrization curve. * *****************************************************************************/ static CagdCrvStruct *CreateReparamCrv(int *Vect, int *Vecr, int Order, int Len, int Reduce) { int i, j; CagdRType Delta, r, *Kv, NewVecr[2 * SAMPLE_NUMBER], NewVect[2 * SAMPLE_NUMBER]; CagdCrvStruct *Crv; CagdCtlPtStruct *Mesh, *PtList = NULL; FixVector(Vect, NewVect, Len); FixVector(Vecr, NewVecr, Len); VectorTransform(NewVect, 0, SAMPLE_NUMBER - 1, Len); VectorTransform(NewVecr, 0, SAMPLE_NUMBER - 1, Len); for (i = Len - 1; i >= 0; i--) { Mesh = CagdCtlPtNew(CAGD_PT_E1_TYPE); Mesh -> Coords[1] = NewVecr[i]; LIST_PUSH(Mesh, PtList); } Delta = ((CagdRType) Len) / (Reduce - Order); Kv = (CagdRType *) IritMalloc((Reduce + Order) * sizeof(CagdRType)); for (i = j = 0; i < Order; i++, j++) Kv[j] = NewVect[0]; for (r = Delta / 2; j < Reduce; r += Delta) Kv[j++] = NewVect[(int) r]; for (i = 0; i < Order; i++) Kv[j++] = NewVect[Len - 1]; /* Interpolate a bspline curve from the points. */ Crv = BspCrvInterpolate(PtList, Len, NewVect, Kv, Reduce, Order, FALSE); CagdCtlPtFreeList(PtList); IritFree(Kv); /* Fix the curve to be monotone. */ FixCrv(Crv); /* Return the new parametrization curve. */ return Crv; } /***************************************************************************** * DESCRIPTION: * * Fix the input integer vector, so that NewVec is increasing monotone. * * * * PARAMETERS: * * OldVec: The input vector. * * NewVec: The output (fixed) vector * * Len: The length of the vector. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void FixVector(int *OldVec, CagdRType *NewVec, int Len) { int j, i = 1; NewVec[0] = OldVec[0]; while (i < Len) { if (OldVec[i] != OldVec[i - 1]) { NewVec[i] = OldVec[i]; i++; } else { CagdRType Delta; for (j = i; OldVec[j] == OldVec[i - 1] && j < Len; j++); if (j >= Len) Delta = 1.0 / (j - i + 1.0); else Delta = (OldVec[j] - OldVec[i]) / ((CagdRType) (j - i + 1)); j--; for ( ; i <= j; i++) NewVec[i] = NewVec[i - 1] + Delta; } } if (!APX_EQ(NewVec[Len - 1], OldVec[Len - 1])) { CagdRType R = OldVec[Len - 1] / NewVec[Len - 1]; for (i = 0; i < Len; i++) NewVec[i] *= R; } for (i = 1; i < Len; i++) if (NewVec[i - 1] > NewVec[i]) fprintf(stderr, "CrvMatch: FixVector: Resulting vector is not monotone.\n"); } /***************************************************************************** * DESCRIPTION: * * Fix the input curve to be monotone. * * * * PARAMETERS: * * Crv: The input curve to be fixed. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void FixCrv(CagdCrvStruct *Crv) { int i, j; CagdRType t; for (j = 0; j < Crv -> Length; j++) for (i = 0; i < Crv -> Length - 1; i++) if (Crv -> Points[1][i] > Crv -> Points[1][i + 1]) { t = Crv -> Points[1][i]; Crv -> Points[1][i] = Crv -> Points[1][i + 1]; Crv -> Points[1][i + 1] = t; } } /***************************************************************************** * DESCRIPTION: * * Affine transform a set of points, so its new end locations are NewBegin * * and NewEnd, respectively. * * * * PARAMETERS: * * Poly: A pointer to points to change. * * Len: The length of the input poly. * * NewBegin: The new begin. * * NewEnd: The new end. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void PolyTransform(CagdRType **Poly, int Len, CagdRType NewBegin, CagdRType NewEnd) { int i; CagdRType Begin, End; Begin = Poly[1][0]; End = Poly[1][Len - 1]; for (i = 0; i < Len; i++) Poly[1][i] = (Poly[1][i] - Begin) * (NewEnd - NewBegin) / (End - Begin) + NewBegin; } /***************************************************************************** * DESCRIPTION: * * Affine transform a set of vectors, so its new end locations are NewBegin * * and NewEnd, respectively. * * * * PARAMETERS: * * Vec: The input and the output vector. * * NewBegin: The new begin. * * NewEnd: The new end. * * Len: The length of the input vector. * * * * RETURN VALUE: * * void * *****************************************************************************/ static void VectorTransform(CagdRType *Vec, CagdRType NewBegin, CagdRType NewEnd, int Len) { int i; CagdRType Begin, End; Begin = Vec[0]; End = Vec[Len - 1]; for (i = 0; i < Len; i++) Vec[i] = (Vec[i] - Begin) * (NewEnd - NewBegin) / (End - Begin) + NewBegin; } /***************************************************************************** * DESCRIPTION: M * Gets two freeform curves, Crv1, nd Crv2, computes a new parametrization M * to Crv2 using composition between Crv2 and a computed reparametrization M * that establishes a matching correspondance between Crv1 and Crv2. M * * * PARAMETERS: M * Crv1: The first curve. M * Crv2: The second curve. M * Reduce: The degrees of freedom of the reparametrization curve. The M * larger this number is, the better the reparametrization will be M * at the cost of more computation. Must be less than SampleSet M * SampleSet: Number of samples the two curves are sampled at. The M * larger this number is, the better the reparametrization will be M * ReparamOrder: Order of reparametrization curve. M * RotateFlag: use or not use rotation in finding best matching M * AllowNegativeNorm: If TRUE, negative norms are locally allowed. M * MatchNormFunc: A pointer to the matching norm. M * * * RETURN VALUE: M * CagdCrvStruct *: The second curve, Crv2, after reparametrization that M * matches the first curve. M * * * KEYWORDS: M * CagdMatchingTwoCurves, correspondance, matching M *****************************************************************************/ CagdCrvStruct *CagdMatchingTwoCurves(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2, int Reduce, int SampleSet, int ReparamOrder, int RotateFlag, int AllowNegativeNorm, CagdMatchNormFuncType MatchNormFunc) { int i, Len, Shift, NewIndex1[2 * SAMPLE_NUMBER], NewIndex2[2 * SAMPLE_NUMBER]; CagdVType Crv1Eval[SAMPLE_NUMBER], Crv2Eval[SAMPLE_NUMBER], TanVecCrv1[SAMPLE_NUMBER], TanVecCrv2[SAMPLE_NUMBER]; CagdRType StartCrv2Par, StartCrv1Par, EndCrv2Par, EndCrv1Par, ShiftPar, TangentsDotProdRes[SAMPLE_NUMBER]; CagdCrvStruct *TCrv, *NewParaCrv; Crv2 = CagdCrvCopy(Crv2); GlblNumberOfSamples = (SampleSet > SAMPLE_NUMBER) ? SAMPLE_NUMBER : SampleSet; GlblRotationFlag = RotateFlag; GlblShiftFlag = IfCrvIsClose(Crv2) ? GlblNumberOfSamples : 1; GlblAllowNegativeNorm = AllowNegativeNorm; if (MatchNormFunc == NULL) MatchNormFunc = CagdMatchMorphNorm; /* Find the tangents at each sample point. */ if (FindTangentsVec(Crv1, TanVecCrv1, Crv1Eval) == 0 || FindTangentsVec(Crv2, TanVecCrv2, Crv2Eval) == 0) { CAGD_FATAL_ERROR(CAGD_ERR_CANNOT_COMP_VEC_FIELD); return NULL; } /* Init. the new matching vectors */ for (i = 0; i < GlblNumberOfSamples; i++) NewIndex1[i] = NewIndex2[i] = i; NewIndex1[GlblNumberOfSamples] = NewIndex2[GlblNumberOfSamples] = -1; MatchingTest(TanVecCrv1, TanVecCrv2, TangentsDotProdRes, NewIndex1, NewIndex2); if ((Shift = FindBestMatching(TanVecCrv1, TanVecCrv2, Crv1Eval, Crv2Eval, NewIndex1, NewIndex2, MatchNormFunc)) == -1) { TCrv = CagdCrvReverse(Crv2); CagdCrvFree(Crv2); Crv2 = TCrv; if (FindTangentsVec(Crv1, TanVecCrv1, Crv1Eval) == 0 || FindTangentsVec(Crv2, TanVecCrv2, Crv2Eval) == 0) { CAGD_FATAL_ERROR(CAGD_ERR_CANNOT_COMP_VEC_FIELD); return NULL; } for (i = 0; i < GlblNumberOfSamples; i++) NewIndex1[i] = NewIndex2[i] = i; NewIndex1[GlblNumberOfSamples] = NewIndex2[GlblNumberOfSamples] = -1; if ((Shift = FindBestMatching(TanVecCrv1, TanVecCrv2, Crv1Eval, Crv2Eval, NewIndex1, NewIndex2, MatchNormFunc)) == -1) { return NULL; } } if (Shift != 0) { CagdCrvStruct *H1Crv2, *H2Crv2; CagdCrvDomain(Crv2, &StartCrv2Par, &EndCrv2Par); ShiftPar = StartCrv2Par + ((CagdRType) (GlblNumberOfSamples - Shift)) * (EndCrv2Par - StartCrv2Par) / ((CagdRType) GlblNumberOfSamples); H2Crv2 = CagdCrvRegionFromCrv(Crv2, StartCrv2Par , ShiftPar); H1Crv2 = CagdCrvRegionFromCrv(Crv2, ShiftPar, EndCrv2Par); CagdCrvFree(Crv2); Crv2 = CagdMergeCrvCrv(H1Crv2, H2Crv2, TRUE); CagdCrvFree(H1Crv2); CagdCrvFree(H2Crv2); BspKnotAffineTrans2(Crv2 -> KnotVector, Crv2 -> Length + Crv2 -> Order, StartCrv2Par, EndCrv2Par); } for (Len = 0; Len < 2 * GlblNumberOfSamples && NewIndex1[Len+1] != -1; Len++); #ifdef DEBUG_MATCH_2CRV for (i = 0; i < 2 * NumberOfSample && NewIndex1[i] != -1; fprintf(stderr, "\n%d %d", NewIndex1[i], NewIndex2[i]), i++); #endif /* DEBUG_MATCH_2CRV */ NewParaCrv = CreateReparamCrv(NewIndex1, NewIndex2, ReparamOrder, Len, Reduce); /* Find the domain of the two curves. */ CagdCrvDomain(Crv1, &StartCrv1Par, &EndCrv1Par); CagdCrvDomain(Crv2, &StartCrv2Par, &EndCrv2Par); /* Scale the matching curve. */ VectorTransform(NewParaCrv -> KnotVector, StartCrv1Par, EndCrv1Par, NewParaCrv -> Length + NewParaCrv -> Order); PolyTransform(NewParaCrv -> Points, NewParaCrv -> Length, StartCrv2Par, EndCrv2Par); TCrv = SymbComposeCrvCrv(Crv2, NewParaCrv); CagdCrvFree(NewParaCrv); CagdCrvFree(Crv2); return TCrv; }