home *** CD-ROM | disk | FTP | other *** search
/ vis-ftp.cs.umass.edu / vis-ftp.cs.umass.edu.tar / vis-ftp.cs.umass.edu / pub / Software / ASCENDER / poly2curves.c < prev    next >
C/C++ Source or Header  |  1996-04-24  |  16KB  |  713 lines
  1. #include "../polygons.h"
  2.  
  3.  
  4. #define EXPECTED_CONVERGENCE    250
  5. #define EPSILON         0.0001
  6.     
  7. int singular(float **M, int n);
  8.  
  9. #define TINY 1.0e-8
  10.  
  11.  
  12.  
  13. /******************************************************************************
  14.  
  15.             Christopher Jaynes
  16.               Feb 7, 1994
  17.  
  18.  
  19.     These routines take the final polygon set from the Feature
  20.     Relation Graph and draw the polygons selected as the maximal
  21.     independent subset of cycles and display them as an overlay on
  22.     the image.
  23.  
  24.     The polygons
  25.  
  26. ***************************************************************************/
  27. void constraints_3d(vertexList *vert);
  28.  
  29.  
  30. void poly2curves(c_handle pane,  MGvertexList *MetaGraph, c_handle im)
  31. {
  32.  
  33.     MGedgeList *MG;
  34.     vertexList  *vert;
  35.     vertexList  *ptr;
  36.     int  index=0;
  37.     int i=0;
  38.  
  39.     c_handle world;
  40.     c_handle RCDEcurve;
  41.     c_handle polygons_fs;
  42.  
  43.  
  44.  
  45.     MG = MetaGraph->eList;
  46.  
  47.     if (MG == NULL)
  48.         fprintf(stderr,"ERRORRRRRRR!\n");
  49.  
  50.     /** Make the 2d Feature Set to contain the rooptop polygons **/
  51.     world = get_2d_world(im);
  52.     polygons_fs = get_fs_named("UMass 2d Buildings",world);
  53.  
  54.  
  55.     while (MG != NULL) {
  56.         vert = MG->vertex->vList;
  57.         if (vert == NULL)
  58.             fprintf(stderr,"Polygon with no Vertecies !\n");
  59.         /** Impose constraints on this polygon **/
  60.         if (CONSTRAIN) 
  61.             constraints_3d(vert); 
  62.           i=0;
  63.           RCDEcurve = make_2d_building_curve(TRUE,world);
  64.           while (vert != NULL)  {
  65.             if (CONSTRAIN) 
  66.                add_vertex_internal(RCDEcurve,i,
  67.                     vert->vertex->xpos,
  68.                     vert->vertex->ypos,0.0);
  69.             else
  70.                add_vertex_internal(RCDEcurve,i,
  71.                 (double)vert->vertex->x,
  72.                 (double)vert->vertex->y,0.0);
  73.             vert = vert->next;
  74.             i++;
  75.           }
  76.           while (number_of_vertices(RCDEcurve) > i) 
  77.             delete_vertex_internal(RCDEcurve,i+1);
  78.  
  79.           add_object(RCDEcurve,polygons_fs);
  80.           index++;
  81.           MG = MG->next;
  82.     } 
  83.     expose_feature_set(polygons_fs,top_view(1,in_pane));
  84.  
  85. }
  86.  
  87.  
  88.  
  89. /* 
  90.  * Constraints:
  91.  *
  92.  * Use world knowledge about the shape of the polygons to choose a
  93.  * least squares solution for the position of the vertecies according
  94.  * to orthogonality, parallelism.
  95.  *
  96.  */
  97.  
  98. #define MAXVERT        20
  99.  
  100. float V[2*MAXVERT];           /* Vector of variables   (x,y) */
  101. float M[2*MAXVERT];           /* Vector of measurments (u,v) */
  102. float C[MAXVERT];             /* Constraint equations (evaluated at M) */
  103. float L[MAXVERT];             /* Lambda                        */
  104. float J[MAXVERT][2*MAXVERT];  /* Jacobian of C, with respect to V      */
  105. float R[MAXVERT];             /* J x M */
  106. float Theta[MAXVERT];         /* Angle Constraints imposed on the polygon */
  107. float Z[MAXVERT];
  108.  
  109. float **A;
  110. float *B;
  111.  
  112. /** Function Prototypes **/
  113. float nextX(int i, int n);
  114. float prevX(int i, int n);
  115. float nextY(int i, int n);
  116. float prevY(int i, int n);
  117. float norm(float x1,float y1,float x2,float y2);
  118. float Mx(int i);
  119. float My(int i);
  120. float error(int n);
  121. float normN(int i, int n);
  122. float normP(int i, int n);
  123.  
  124. void computeR(int n);
  125. void buildA(int n);
  126. void buildB(int n);
  127. void loadC(int n);
  128. void loadJ(int n);
  129. void initA(int n);
  130. void init(int n);
  131.  
  132.  
  133.  
  134. float **As, **vs;
  135. float *ws;
  136.  
  137. void constraints_3d(vertexList *vert)
  138. /* 
  139.  * Backproject the polygon edges into a horizonal plane, 
  140.  * impose the constraints on the final 3D shape of the polygon
  141.  * using a least-squares method.
  142.  * 
  143.  * Re-project the solution to capture the final polygon shape
  144.  * in the image coordinates.
  145.  *
  146.  */
  147. {
  148.  
  149.     int N;    /* Number of datapoints */
  150.      int *Index;
  151.         float d;
  152.         float e;
  153.     mat_4x1 point, result;
  154.         int i,j,k;
  155.     int count=0;
  156.         
  157.  
  158.  
  159.      N = loadM_in3D(vert);
  160.  
  161.     j=0;
  162.  
  163.         Index = (int *)vector(1,3*N);
  164.         A = matrix(1,3*N,1,3*N);
  165.         B = vector(1,3*N);
  166.  
  167.         e = 1.0;
  168.  
  169.         while (e > EPSILON) {
  170.                 init(N);
  171.                 loadC(N);
  172.  
  173.                 loadJ(N);
  174.  
  175.                 computeR(N);
  176.  
  177.                 initA(N);
  178.                 buildA(N);
  179.  
  180.         if (singular(A,3*N)) {
  181.             while (vert != NULL) {
  182.                 vert->vertex->xpos = (double)vert->vertex->x;
  183.                 vert->vertex->ypos = (double)vert->vertex->y;
  184.                 vert = vert->next;
  185.             }
  186.             return;
  187.         }
  188.                 ludcmp(A,(3*N),Index,&d);
  189.                 for (j=1; j <= N; j++)
  190.                         d *= A[j][j];
  191.                 if (fabs(d) <= EPSILON){/* If so, return leave polygon as is. */
  192.              while (vert != NULL) {
  193.                                 vert->vertex->xpos = (double)vert->vertex->x;
  194.                                 vert->vertex->ypos = (double)vert->vertex->y;
  195.                                 vert = vert->next;
  196.                         }
  197.                         return;
  198.         }
  199.  
  200.                 buildB(N);
  201.  
  202.                 lubksb(A,(3*N),Index,B);
  203.                 
  204.                 e = error(N);
  205.                 loadM_in3D(NULL,N);
  206.  
  207.         if (count > EXPECTED_CONVERGENCE)
  208.             break;
  209.         count++;
  210.  
  211.          
  212.         }
  213.  
  214.     /** NOTE: The solution has been computed in a 3D, horizontal plane.
  215.         reproject the points into the image coordinate system. **/
  216.     j = 1;
  217.     i=0;
  218.     while (vert != NULL) {
  219.         point[0] = B[j++];
  220.         point[1] = B[j++];
  221.         point[2] = Z[i++];
  222.         point[3] = 1.0;
  223.         project_3d_to_2d(IMAGE_NAME,point,result);
  224.         vert->vertex->xpos = result[0];
  225.         vert->vertex->ypos = result[1];
  226. /*
  227.         vert->vertex->x = (int) floor(result[0]);
  228.         vert->vertex->y = (int) floor(result[1]);
  229. */
  230.         vert = vert->next;
  231.     }
  232.  
  233.     free_matrix(A,1,3*N,1,3*N);
  234.     free_vector(B,1,3*N);
  235.     free_matrix(As,1,3*N,1,3*N);
  236.     free_matrix(vs,1,3*N,1,3*N);
  237.     free_vector(ws,1,3*N);
  238. }
  239.  
  240.  
  241.  
  242. int loadM_in3D(vertexList *v, int n)
  243. {
  244.  
  245.     int i = 0;
  246.     int j;
  247.     float temp;
  248.     mat_4x1 point, result;
  249.  
  250.     /* When v is NULL, we should copy the current solution
  251.        from B, into M */
  252.  
  253.     i=0;
  254.     if (v == NULL) { 
  255.        for (i=0; i < 2*n; i++) 
  256.         M[i] = B[i+1];
  257.     } else {
  258.        i=0; j = 0;    
  259.        while (v != NULL) {
  260.         point[0] = (double) v->vertex->x;
  261.         point[1] = (double) v->vertex->y;
  262.         point[2] = 0.0;
  263.         point[3] = 1.0;
  264.         project_2d_to_KP(IMAGE_NAME,point,result);
  265.         M[j++] = (float) result[0];
  266.         M[j++] = (float) result[1];
  267.         Z[i] = result[2];
  268.         Theta[i++] = M_PI_2;
  269.         v = v->next;
  270.        }
  271.     }
  272.     return(i);
  273. }
  274.  
  275.  
  276. float prevX(int i, int n)
  277. {
  278.         if (i == 0)
  279.                 return(M[(2*n)-2]);
  280.         else
  281.                 return(M[(2*i)-2]);
  282. }
  283.  
  284. float nextX(int i, int n)
  285. {
  286.         if (i == (n-1))
  287.                 return(M[0]);
  288.         else
  289.                 return(M[(2*i)+2]);
  290. }
  291.  
  292. float prevY(int i, int n)
  293. {
  294.         if (i == 0)
  295.                 return(M[(2*n)-1]);
  296.         else
  297.                 return(M[(2*i)-1]);
  298. }
  299.  
  300. float nextY(int i, int n)
  301. {
  302.         if (i == (n-1))
  303.                 return(M[1]);
  304.         else {
  305.                 return(M[(2*i)+3]);
  306.         }
  307. }
  308.  
  309.  
  310. void loadC(int n)
  311. /**  
  312.         Compute 'C' the Constraint equations at 'M'
  313.  
  314.         Intersecting lines on the polygon are constrained to
  315.         meet at the correct angle( right now 90 degrees !).  That is:
  316.                 L1 -- L2 meet at angle Theta[1].
  317. **/
  318. {
  319.  
  320.         int i, j;
  321.  
  322.         for (i=0; i < n; i++) {
  323.  
  324.                 C[i] = ((nextX(i,n)-Mx(i))*(Mx(i)-prevX(i,n))) +
  325.                         ((nextY(i,n)-My(i))*(My(i)-prevY(i,n))) -
  326.                         (normN(i,n)*normP(i,n)*cos(Theta[i]));
  327.         }
  328.  
  329. }
  330.  
  331. void loadJ(int n)
  332. /**
  333.         J is is Jacobian of C with respect to the matrix V (at M)
  334. **/
  335. {        
  336.  
  337.         int i, j;
  338.  
  339.         J[0][0] =  -2.0*Mx(0) + prevX(0,n) + nextX(0,n)
  340.                 -0.5*((normP(0,n)*Theta[0]*(-2.0*nextX(0,n) + 2.0*Mx(0))) /
  341.                      (normN(0,n)))
  342.                 -0.5*((normN(0,n)*Theta[0]*(-2.0*prevX(0,n) + 2.0*Mx(0))) /
  343.                      (normP(0,n)));
  344.  
  345.         J[0][1] =  -2.0*My(0) + prevY(0,n) + nextY(0,n)
  346.                 -0.5*((normP(0,n)*Theta[0]* (-2.0*nextY(0,n) + 2.0*My(0))) /
  347.                      (normN(0,n)))
  348.                 -0.5*((normP(0,n)*Theta[0]* (-2.0*prevY(0,n) + 2.0*My(0))) /
  349.                      (normN(0,n)));
  350.  
  351.         J[0][2] = -prevX(0,n) + Mx(0) - 0.5*((normP(0,n)*Theta[0]*
  352.                         (2.0*nextX(0,n) - 2.0*Mx(0)) ) / normN(0,n));
  353.  
  354.  
  355.         J[0][3] = -prevY(0,n) + My(0) - 0.5*((normP(0,n) * Theta[0]*
  356.                         (2.0*nextY(0,n) - 2.0*My(0)) ) / normN(0,n));
  357.  
  358.         J[0][(2*n)-2] = -nextX(0,n) + Mx(0) -
  359.                  0.5*((normN(0,n)*Theta[0]* (2.0*prevX(0,n) - 2.0*Mx(0))) /
  360.                         normP(0,n));
  361.         J[0][(2*n)-1] = -nextY(0,n) + My(0) -
  362.                  0.5*((normN(0,n)*Theta[0]* (2.0*prevY(0,n) - 2.0*My(0))) /
  363.                         normP(0,n));
  364.     for (i=1,j=0; i < n-1; i++,j= j + 2) {
  365.                 J[i][j]  = -nextX(i,n) + Mx(i) -
  366.                  0.5*((normN(i,n)*Theta[i]* (2.0*prevX(i,n) - 2.0*Mx(i))) /
  367.                         normP(i,n));
  368.                 J[i][j+1] = -nextY(i,n) + My(i) -
  369.                  0.5*((normN(i,n)*Theta[i]* (2.0*prevY(i,n) - 2.0*My(i))) /
  370.                         normP(i,n));
  371.  
  372.                 J[i][j+2] =  -2.0*Mx(i) + prevX(i,n) + nextX(i,n)
  373.                 -0.5*((normP(i,n)*Theta[i]* (-2.0*nextX(i,n) + 2.0*Mx(i))) /
  374.                      (normN(i,n)))
  375.                 -0.5*((normN(i,n)*Theta[i]* (-2.0*prevX(i,n) + 2.0*Mx(i))) /
  376.                      (normP(i,n)));
  377.                 J[i][j+3] =  -2.0*My(i) + prevY(i,n) + nextY(i,n)
  378.                 -0.5*((normP(i,n)*Theta[i]* (-2.0*nextY(i,n) + 2.0*My(i))) /
  379.                      (normN(i,n)))
  380.                 -0.5*((normN(i,n)*Theta[i]* (-2.0*prevY(i,n) + 2.0*My(i))) /
  381.                      (normP(i,n)));
  382.  
  383.                 J[i][j+4] = -prevX(i,n) + Mx(i) -
  384.                          0.5*((normP(i,n)* Theta[i]*(2.0*nextX(i,n) - 2.0*Mx(i))) /
  385.                                 normN(i,n));
  386.  
  387.                 J[i][j+5] = -prevY(i,n) + My(i) -
  388.                          0.5*((normP(i,n)* Theta[i]*(2.0*nextY(i,n) - 2.0*My(i))) /
  389.                                 normN(i,n));
  390.         }
  391.       J[i][0] = -prevX(i,n) + Mx(i) - 0.5*((normP(i,n)*Theta[i]*
  392.                         (2.0*nextX(i,n) - 2.0*Mx(i)) ) /
  393.                                 normN(i,n));
  394.         J[i][1] = -prevY(i,n) + My(i) - 0.5*((normP(i,n)*Theta[i]*
  395.                         (2.0*nextY(i,n) - 2.0*My(i)) ) /
  396.                                 normN(i,n));
  397.  
  398.         J[i][j] = -nextX(i,n) + Mx(i) - 0.5*((normN(i,n) *Theta[i]*
  399.                          (2.0*prevX(i,n) - 2.0*Mx(i))) /
  400.                         normP(i,n));
  401.  
  402.         J[i][j+1] = -nextY(i,n) + My(i) - 0.5*((normN(i,n) *Theta[i]*
  403.                          (2.0*prevY(i,n) - 2.0*My(i))) /
  404.                         normP(i,n));
  405.  
  406.         J[i][j+2] =  -2.0*Mx(i) + prevX(i,n) + nextX(i,n)
  407.                 -0.5*((normP(i,n)*Theta[i]* (-2.0*nextX(i,n) + 2.0*Mx(i))) /
  408.                      (normN(i,n)))
  409.                 -0.5*((normN(i,n)*Theta[i]* (-2.0*prevX(i,n) + 2.0*Mx(i))) /
  410.                      (normP(i,n)));
  411.  
  412.  
  413.         J[i][j+3] =  -2.0*My(i) + prevY(i,n) + nextY(i,n)
  414.                 -0.5*((normP(i,n)*Theta[i]* (-2.0*nextY(i,n) + 2.0*My(i))) /
  415.                      (normN(i,n)))
  416.                 -0.5*((normN(i,n)*Theta[i]* (-2.0*prevY(i,n) + 2.0*My(i))) /
  417.                      (normP(i,n)));
  418.  
  419. }
  420.  
  421. void buildB(int n)
  422. /* -2m
  423.                 Cm -Jm *M
  424.  
  425.         Here, Jm*M is stored in the matrix R,  result of the
  426.         Jacobian(nx2n) * M (2nx1).
  427.  
  428.         The matrix B, is used to solve the final linearized system
  429.         for V.  The solution is then used as another inital estimate
  430.         for the linearization of C.
  431. **/
  432.  
  433. {
  434.         int i,j;
  435.  
  436.         for (j=0,i=1; i <= 2*n; j++,i++)
  437.                 B[i] = 2.0*M[j];
  438.         for(j=0,i=2*n+1; i <= 3*n; i++,j++) {
  439.                 B[i] =  R[j] - C[j];
  440.         }
  441. }
  442.  
  443.  
  444.  
  445. void buildA(int n)
  446. {
  447.  
  448.         int i, j, k;
  449.         int l;
  450.  
  451.         /** Upper left is 2*I **/
  452.         for (i=1; i <= (2*n); i++)
  453.                 A[i][i] = 2.0;
  454.  
  455.         /** Lower left is the jacobian J **/
  456.         for (k=0,i=(2*n)+1; i <= 3*n; k++,i++)
  457.            for(l=0,j=1; j <= (2*n); l++,j++)
  458.                 A[i][j] = J[k][l];
  459.  
  460.         /** Lower right is zero **/
  461.         for (i=(2*n)+1; i <= 3*n; i++)
  462.            for (j=(2*n)+1; j <= 3*n; j++)
  463.                 A[i][j] = 0.0;
  464.  
  465.         /**  Upper right is transpose of J **/
  466.         for (k=0,i=(2*n)+1; i <= 3*n; k++,i++)
  467.            for (l=0,j=1; j <= 2*n; l++,j++)
  468.                 A[j][i] = J[k][l];
  469. }
  470.  
  471. void
  472. computeR(int n)
  473. {
  474.  
  475. /* R is the Matrix - Vector Product of the Jacobian J and the parameter
  476.    vector V */
  477.      
  478.  
  479.         int i;
  480.         int j;
  481.  
  482.         for (i=0; i < n; i++)
  483.                 for (j=0; j < (2*n); j++)
  484.                         R[i] += J[i][j] * M[j];
  485. }
  486.  
  487. void initA(int n)
  488. {
  489.  
  490.         int i, j;
  491.  
  492.         for (i=1; i <=(3*n); i++)
  493.            for (j=1; j <=(3*n); j++)
  494.                 A[i][j] = 0.0;
  495. }
  496.  
  497. void init(int n)
  498. {
  499.          
  500.         int i,j;
  501.  
  502.         for (i=0; i < n; i++)
  503.            for (j = 0; j < 2*n; j++)
  504.                 J[i][j] = 0.0;
  505.  
  506.         for (i=0; i < n; i ++)
  507.                 R[i] = 0.0;
  508.  
  509. }
  510.  
  511.  
  512. float Mx(int i)
  513. {
  514.         return(M[2*i]);
  515. }
  516.  
  517. float My(int i)
  518. {
  519.         return(M[(2*i)+1]);
  520. }
  521.  
  522.  
  523. float norm(float x1,float y1,float x2,float y2)
  524. {
  525.  
  526.         return( (float)sqrt( ((x1-x2)*(x1-x2)) + ((y1-y2)*(y1-y2))));
  527. }
  528.  
  529. float error(int n)
  530. {
  531.  
  532.         int     i;
  533.         float e;
  534.         
  535.     e = 0.0;
  536.  
  537.         for (i = 0; i < (2*n); i++) {
  538.                 e += (M[i]-B[i+1])*(M[i] - B[i+1]);
  539.         }
  540.  
  541.         return(sqrt(e));
  542. }
  543.  
  544.  
  545. float normP(int i, int n)
  546. {
  547.  
  548.         return( norm(prevX(i,n),prevY(i,n),Mx(i),My(i)) );
  549. }
  550. float normN(int i, int n)
  551. {
  552.  
  553.         return( norm(nextX(i,n),nextY(i,n),Mx(i),My(i)) );
  554. }
  555.  
  556.  
  557.  
  558. /*************************************
  559.     Get 3d stuff -----
  560.  
  561.  
  562. ****/
  563.  
  564. void poly2_3dcurves(c_handle pane, MGvertexList *MetaGraph)
  565. {
  566.  
  567.  
  568.     c_handle world3d;
  569.     c_handle *features;
  570.     c_handle RCDEcurve;
  571.     c_handle im;
  572.     c_handle view;
  573.     c_handle obj_matrix;
  574.     c_handle obj_trans;
  575.     c_handle extrudedCurve;
  576.     mat_4x1 p, result;
  577.  
  578.     MGedgeList *MG;
  579.         vertexList  *vert;
  580.         vertexList  *tempPtr;
  581.         int  index=0;
  582.         int i=0;
  583.     int closed=1;
  584.     c_handle poly_fs;
  585.     double avgHeight=0.0;
  586.     int count;
  587.  
  588.  
  589.         world3d = get_3d_world_named(IMAGE_NAME_3D);
  590.     view = top_view(1,in_pane);
  591.     poly_fs = get_fs_named("UMass 3d Buildings",world3d);
  592.  
  593.  
  594.     MG = MetaGraph->eList;
  595.  
  596.     if (MG == NULL)
  597.         fprintf(stderr,"ERRORRRRRRR!\n");
  598.  
  599.     while (MG != NULL) {
  600.         vert = MG->vertex->vList;
  601.         if (vert == NULL)
  602.             fprintf(stderr,"Polygon with no Vertecies !\n");
  603.         /** Impose constraints on this polygon  **/
  604.         if (CONSTRAIN)
  605.             constraints_3d(vert); 
  606.          
  607.  
  608.  
  609.  
  610.             tempPtr = vert;
  611.         avgHeight = 0.0;
  612.         count = 0;
  613.           while (tempPtr != NULL) {
  614.             avgHeight += vert->vertex->height;
  615.             count++;
  616.             tempPtr = tempPtr->next;
  617.               } 
  618.           avgHeight = avgHeight / (double)count;
  619.  
  620.  
  621.           i=0;
  622.           RCDEcurve = make_3d_closed_curve(6,&world3d);
  623.           while (vert != NULL)  {
  624.             p[0] = (double)vert->vertex->x;
  625.             p[1] = (double)vert->vertex->y;
  626.             p[2] = avgHeight;
  627.             p[3] = 1.0;
  628.             
  629.             project_2d_to_KP(IMAGE_NAME,p,result);
  630.             add_vertex_internal(RCDEcurve,i,
  631.                     result[0],result[1],
  632.                     avgHeight);
  633.             vert = vert->next;
  634.             i++;
  635.           }
  636.           while (number_of_vertices(RCDEcurve) > i) 
  637.             delete_vertex_internal(RCDEcurve,i+1);
  638.  
  639.         select_fs(world3d,poly_fs);
  640.         extrudedCurve = extrude_roof_curve(RCDEcurve);
  641.         index++;
  642.         MG = MG->next;
  643.      } 
  644.     expose_feature_set(poly_fs,top_view(1,in_pane));
  645.  
  646. }
  647.  
  648.  
  649. c_handle get_fs_named(char *Fname, c_handle wld)
  650. {
  651.  
  652.     c_handle *features;
  653.     int    done=FALSE;
  654.  
  655.     
  656.     features = (c_handle *)get_feature_sets(wld,0);
  657.  
  658.     while (!done && *features != -1) {
  659.         if (strcmp((char *)name(*features),Fname) == 0)
  660.             return(*features);
  661.         else 
  662.             features++;
  663.     }
  664. }
  665.  
  666.  
  667.  
  668. int singular(float **M, int n)
  669. /*
  670.  * Compute a determinant of matrix M, size n.
  671.  *
  672.  */
  673. {
  674.  
  675.    int i,j;
  676.    float big,temp;
  677.  
  678.    As = matrix(1,n,1,n);
  679.    vs = matrix(1,n,1,n);
  680.    ws = vector(1,n);
  681.  
  682.    for (i=1; i <=n; i++) 
  683.      for (j=1; j <=n; j++)
  684.         As[i][j] = M[i][j];
  685.  
  686.    svdcmp(As,n,n,ws,vs);
  687.  
  688.    for (i=1; i <n; i++) {
  689.     if (ws[i] < TINY) {
  690.         return(TRUE);
  691.     }
  692.  
  693.    }
  694.  
  695.  
  696.    return(FALSE);
  697.  
  698. /**
  699.    
  700.    for (i=1; i < n; i++) {
  701.     big = 0.0;
  702.     for (j=1; j<=n; j++)
  703.         if ((temp=fabs(M[i][j])) > big) big=temp;
  704.     if (big < TINY) {
  705.         return(TRUE);
  706.     }
  707.    }
  708.  
  709.    return(FALSE);
  710. ***/
  711. }
  712.  
  713.