/**************************************************************** * This program seeks the equation of a conical passing through * * 5 given points. * * ------------------------------------------------------------- * * SAMPLE RUN: * * * * SEEK A CONICAL PASSING THROUGH 5 POINTS * * * * X1 Y1 = 0 2 * * X2 Y2 = 1 5 * * X3 Y3 = 5 77 * * X4 Y4 = -4 50 * * X5 Y5 = 2 14 * * * * * * The coefficients of the conical ax2+bxy+cy2+dx+ey+f = 0 are: * * * * a = 1.500000 * * b = 0.000000 * * c = -0.000000 * * d = -0.000000 * * e = -0.500000 * * f = 1.000000 * * * * C++ version by J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***************************************************************** Explanations: ------------ We seek the cartesian equation ax^2+bxy+cy^2+dx+ey+f=0 of the conical passing through 5 points: (x1,y1),(x2,y2)... to (x5,y5). The unknowns are a,b,c,d,e,f. Since the equation can be multiplied by any real number, the actual number of independant unknowns is 5. That is the reason why only five different points are necassary to determine the conical. The linear system to solve is a system of five homogeneous equations with 6 unknowns. The resolution of the system will be made by successive steps: Let us suppose that f<>0 and let us fix f=1, then the other parameters a,b,c,d,e are solutions of the linear system: x1^2 a + x1y1 b + y1^2 c + x1 d + y1 e = -1 x2^2 a + x2y2 b + y2^2 c + x2 d + y2 e = -1 ------------------------------------------- x5^2 a + x5y5 b + y5^2 c + x5 d + y5 e = -1 This linear system of size 5 is solved by a classical Gauss-Jordan method found in procedure SolveSystem (see also program MATRICES/sysmat.cpp). if the system has a solution (det <> 0) the problem is finished, else the hypothesis f<>0 was wrong and we then try e<>0 (that is e=1). The linear system to solve is then: x1^2 a + x1y1 b + y1^2 c + x1 d + f = -y1 x2^2 a + x2y2 b + y2^2 c + x2 d + f = -y2 ------------------------------------------- x5^2 a + x5y5 b + y5^2 c + x5 d + f = -y5 if the system has a solution (det <> 0) the problem is finished, else the hypothesis e<>0 was again wrong and we successively explore d<>0, c<>0... to a<>0. So in the most unfavourable case, that is if only a degenerated conical answers the problem, we can have up to six linear systems of size 5 to solve. Usually (f<>0 or e<>0), we only solve one or two linear systems. ----------------------------------------------------------------------------------*/ #include #include #define EPSMACH 2e-16 #define FALSE 0 #define TRUE 1 typedef double MAT5[6][6]; //index 0 not used here typedef double VEC5[6][2]; // " " double x1,yy1,x2,y2,x3,y3,x4,y4,x5,y5; //y1 already defined by Visual C++ double a,b,c,d,e,f; double det; /****************************************** * SOLVING A LINEAR MATRIX SYSTEM AX = B * * with Gauss-Jordan method using full * * pivoting at each step. During the pro- * * cess, original AA and BB matrices are * * destroyed to spare storage location. * * --------------------------------------- * * INPUTS: AA MATRIX N*N * * BB MATRIX N*M * * --------------------------------------- * * OUTPUS: AA INVERSE OF AA N*N * * DET DETERMINANT OF AA * * BB SOLUTION MATRIX N*M * * --------------------------------------- * * NOTA - If M=0 inversion of AA matrix * * only. * * * * C++ version by J-P Moreau * * (here limited to N=5 and M=1). * ******************************************/ void SolveSystem(int N, int M,MAT5 AA,VEC5 BB, double *DET) { double PC[6],PL[6],CS[6]; double PV,PAV,temp,TT; int I,IK,J,JK,K; //LABELS: e50, fin, fin1, fin2 // Initializations *DET = 1.0; for (I=1; I<=N; I++) { PC[I] = 0.0; PL[I] = 0.0; CS[I] = 0.0; } // Main loop for (K=1; K<=N; K++) { // Searching greatest pivot PV=AA[K][K]; IK=K; JK=K; PAV=fabs(PV); for (I=K; I<=N; I++) for (J=K; J<=N; J++) { temp = fabs(AA[I][J]); if (temp > PAV) { PV=AA[I][J]; PAV=fabs(PV); IK=I; JK=J; } } // Search terminated, the pivot is in location I=IK, J=JK. // Memorizing pivot location: PC[K]=JK; PL[K]=IK; // DETERMINANT DET is actualised // If DET=0, ERROR MESSAGE and STOP // Machine dependant EPSMACH equals here 2E-16 if (IK!=K) *DET=-*DET; if (JK!=K) *DET=-*DET; *DET=*DET*PV; temp= fabs(*DET); if (temp < EPSMACH) goto e50; // POSITIONNING PIVOT IN K,K: if (IK!=K) for (I=1; I<=N; I++) { // EXCHANGE LINES IK and K of matrix AA: TT=AA[IK][I]; AA[IK][I]=AA[K][I]; AA[K][I]=TT; } if (M!=0) for (I=1; I<=M; I++) { TT=BB[IK][I]; BB[IK][I]=BB[K][I]; BB[K][I]=TT; } // OT is at correct line if (JK!=K) for (I=1; I<=N; I++) { // EXCHANGE COLUMNS JK and K of matrix AA: TT=AA[I][JK]; AA[I][JK]=AA[I][K]; AA[I][K]=TT; } // The PIVOT is at correct column. // and is located in K,K. // Column K of matrix AA is stored in CS vector // then column K is set to zero. for (I=1; I<=N; I++) { CS[I]=AA[I][K]; AA[I][K]= 0.0; } CS[K]= 0.0; AA[K][K]= 1.0; // Line K of matrix AA is modified: temp = fabs(PV); if (temp < EPSMACH) goto e50; for (I=1; I<=N; I++) AA[K][I]=AA[K][I]/PV; if (M!=0) for (I=1; I<=M; I++) BB[K][I]=BB[K][I]/PV; // Other lines of matrix AA are modified: for (J=1; J<=N; J++) { if (J==K) goto fin; for (I=1; I<=N; I++) // Line J of matrix AA is modified: AA[J][I]=AA[J][I]-CS[J]*AA[K][I]; if (M!=0) for (I=1; I<=M; I++) BB[J][I]=BB[J][I]-CS[J]*BB[K][I]; fin: ;} // Line K is ready. } // of K loop // MATRIX AA INVERSION IS DONE - REARRANGEMENT OF MATRIX AA // EXCHANGE LINES for (I=N; I>0; I--) { IK=int(PC[I]); if (IK==I) goto fin1; // EXCHANGE LINES I and PC(I) of matrix AA: for (J=1; J<=N; J++) { TT=AA[I][J]; AA[I][J]=AA[IK][J]; AA[IK][J]=TT; } if (M!=0) for (J=1; J<=M; J++) { TT=BB[I][J]; BB[I][J]=BB[IK][J]; BB[IK][J]=TT; } // NO MORE EXCHANGE is NECESSARY // Go to next line. fin1: ;} // of loop i=N downto // EXCHANGE COLUMNS for (J=N; J>0; J--) { JK=int(PL[J]); if (JK==J) goto fin2; // EXCHANGE COLUMNS J ET PL(J) of matrix AA: for (I=1; I<=N; I++) { TT=AA[I][J]; AA[I][J]=AA[I][JK]; AA[I][JK]=TT; } // NO MORE EXCHANGE is NECESSARY // Go to next column. fin2: ;} // REARRANGEMENT TERMINATED. e50:;} // of SOLVESYSTEM() //function used by SeekConical() void Solve(int i,MAT5 tm,VEC5 tv, double *a, double *b, double *c, double *d, double *e, double *f) { SolveSystem(5,1,tm,tv,&det); if (det!=0) { switch(i) { case 1: // case f<>0 *a=tv[1][1]; *b=tv[2][1]; *c=tv[3][1]; *d=tv[4][1]; *e=tv[5][1]; *f=1; break; case 2: // case e<>0 *a=tv[1][1]; *b=tv[2][1]; *c=tv[3][1]; *d=tv[4][1]; *e=1; *f=tv[5][1]; break; case 3: // case d<>0 *a=tv[1][1]; *b=tv[2][1]; *c=tv[3][1]; *d=1; *e=tv[4][1]; *f=tv[5][1]; break; case 4: // case c<>0 *a=tv[1][1]; *b=tv[2][1]; *c=1; *d=tv[3][1]; *e=tv[4][1]; *f=tv[5][1]; break; case 5: // case b<>0 *a=tv[1][1]; *b=1; *c=tv[2][1]; *d=tv[3][1]; *e=tv[4][1]; *f=tv[5][1]; break; case 6: // case a<>0 *a=1; *b=tv[1][1]; *c=tv[2][1]; *d=tv[3][1]; *e=tv[4][1]; *f=tv[5][1]; } } } bool SeekConical(double x1,double yy1,double x2,double y2,double x3,double y3,double x4, double y4,double x5,double y5,double *a,double *b,double *c,double *d, double *e,double *f) { //Label: e100 MAT5 tm; VEC5 tv; double x[6],y[6]; int i,j; x[1]=x1; y[1]=yy1; x[2]=x2; y[2]=y2; x[3]=x3; y[3]=y3; x[4]=x4; y[4]=y4; x[5]=x5; y[5]=y5; // verify that the 5 given points are different for (i=2; i<=5; i++) for (j=1; j0 for (i=1; i<=5; i++) { tm[i][1]=x[i]*x[i]; tm[i][2]=x[i]*y[i]; tm[i][3]=y[i]*y[i]; tm[i][4]=x[i]; tm[i][5]=y[i]; tv[i][1]=-1; } Solve(1,tm,tv,a,b,c,d,e,f); if (det!=0) goto e100; // try e<>0 for (i=1; i<=5; i++) { tm[i][1]=x[i]*x[i]; tm[i][2]=x[i]*y[i]; tm[i][3]=y[i]*y[i]; tm[i][4]=x[i]; tm[i][5]=1; tv[i][1]=-y[i]; } Solve(2,tm,tv,a,b,c,d,e,f); if (det!=0) goto e100; // try d<>0 for (i=1; i<=5; i++) { tm[i][1]=x[i]*x[i]; tm[i][2]=x[i]*y[i]; tm[i][3]=y[i]*y[i]; tm[i][4]=y[i]; tm[i][5]=1; tv[i][1]=-x[i]; } Solve(3,tm,tv,a,b,c,d,e,f); if (det!=0) goto e100; // try c<>0 for (i=1; i<=5; i++) { tm[i][1]=x[i]*x[i]; tm[i][2]=x[i]*y[i]; tm[i][3]=x[i]; tm[i][4]=y[i]; tm[i][5]=1; tv[i][1]=-y[i]*y[i]; } Solve(4,tm,tv,a,b,c,d,e,f); if (det!=0) goto e100; // try b<>0 for (i=1; i<=5; i++) { tm[i][1]=x[i]*x[i]; tm[i][2]=y[i]*y[i]; tm[i][3]=x[i]; tm[i][4]=y[i]; tm[i][5]=1; tv[i][1]=-x[i]*y[i]; } Solve(5,tm,tv,a,b,c,d,e,f); if (det!=0) goto e100; // last try a<>0 for (i=1; i<=5; i++) { tm[i][1]=x[i]*y[i]; tm[i][2]=y[i]*y[i]; tm[i][3]=x[i]; tm[i][4]=y[i]; tm[i][5]=1; tv[i][1]=-x[i]*x[i]; } Solve(6,tm,tv,a,b,c,d,e,f); if (det!=0) goto e100; else return 0; e100: return TRUE; } //det<>0 means success void main() { printf("\n SEEK A CONICAL PASSING THROUGH 5 POINTS\n\n"); printf(" X1 Y1 = "); scanf("%lf %lf", &x1,&yy1); printf(" X2 Y2 = "); scanf("%lf %lf", &x2,&y2); printf(" X3 Y3 = "); scanf("%lf %lf", &x3,&y3); printf(" X4 Y4 = "); scanf("%lf %lf", &x4,&y4); printf(" X5 Y5 = "); scanf("%lf %lf", &x5,&y5); printf("\n\n\n"); if (SeekConical(x1,yy1,x2,y2,x3,y3,x4,y4,x5,y5,&a,&b,&c,&d,&e,&f)) { printf(" The coefficients of the conical ax2+bxy+cy2+dx+ey+f = 0 are:\n\n"); printf(" a=%10.6f\n", a); printf(" b=%10.6f\n", b); printf(" c=%10.6f\n", c); printf(" d=%10.6f\n", d); printf(" e=%10.6f\n", e); printf(" f=%10.6f\n", f); } else printf(" No conical found.\n"); printf("\n\n"); } // end of file conical1.cpp