home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / plotting / contour / contour.lha / Contour / contour.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-08-28  |  9.5 KB  |  354 lines
  1. /*** contour.c ***/
  2.  
  3. #include <stdio.h>
  4. #include <strings.h>
  5. #include <math.h>
  6. #include "common.h"
  7. #include "contour.h"
  8. #include "plot.h"
  9. #include "Xdefs.h"
  10.  
  11. /* print out a single contour */
  12. draw_cont()
  13. {
  14.    FILE    *fopen(), *ict;
  15.    extern  double cont_level;
  16.    extern  double xmin,xmax,ymin,ymax;
  17.    extern  plotptr plot_listhead;
  18.    int     npts=0,ncurves=0;
  19.    nodeptr N;
  20.    plotptr P;
  21.  
  22.    /* print out info about the image read in */
  23.    fprintf(stdout,"   Maximum z-value    = %f\n",zmax);
  24.    fprintf(stdout,"   Minimum z-value    = %f\n",zmin);
  25.  
  26.    /* check the desired contour level */
  27.    if ((zmin > cont_level) || (cont_level > zmax)) {
  28.       fprintf(stderr,"Error : The desired contour level %.2f",cont_level);
  29.       fprintf(stderr,"is not within the surface\n");
  30.    } else {
  31.       /* now find the desired contours and print it out */
  32.       find_contour(cont_level);
  33.  
  34.       /* open up the file */
  35.       if ((ict = fopen(CTFILE,"w")) == NULL) {
  36.          fprintf(stderr,"cat: can't open %s\n",CTFILE);
  37.          exit(1);
  38.       }
  39.       /* count the number of curves */
  40.       for (P=plot_listhead; P!=NULL; P=P->next) ncurves++;
  41.  
  42.       /* output the x-y limits and the number of curves */
  43.       fprintf(ict,"%10.5f %10.5f %10.5f %10.5f\n",xmin,xmax,ymin,ymax);
  44.       fprintf(ict,"%10.5f\n",(double)ncurves);
  45.  
  46.       /* print out the data for each contour */
  47.       for (P=plot_listhead; P!=NULL; P=P->next) {
  48.          /* count the number of points first */
  49.          npts = 0;
  50.          for (N=P->nodehead; N!=NULL; N=N->next) npts++;
  51.          /* now print out the info */
  52.          fprintf(ict,"%10.5f\n",(double)npts);
  53.          for (N=P->nodehead; N!=NULL; N=N->next) 
  54.              fprintf(ict,"%10.5f %10.5f\n",N->pt.x,N->pt.y);
  55.       }
  56.  
  57.       /* finished */
  58.       fprintf(stdout,"   The x-y points of the %.2f contour ",cont_level);
  59.       fprintf(stdout,"has been printed in \"%s\"\n",CTFILE);
  60.    }
  61. }
  62.  
  63. /* draw the contours */
  64. draw_view(argc,argv)
  65. int argc;
  66. char *argv[];
  67. {
  68.    FILE   *fopen(), *ips;
  69.    char   *sprintf();
  70.    float  rdfloat();
  71.    extern double zmin, zmax;
  72.    double level,cstep;
  73.    char   c, command[MAXCHAR];
  74.    int    nsteps;
  75.    int    FOUND=FALSE;
  76.  
  77.    /* find all the contours first */
  78.    fprintf(stdout,"   Maximum z-value           = %f\n",zmax);
  79.    fprintf(stdout,"   Minimum z-value           = %f\n",zmin);
  80.  
  81.    /* ask for the contour step size */
  82.    while (!FOUND) {
  83.       fprintf(stdout,"   Enter contour step size :");
  84.       fscanf(stdin,"%lf",&cstep);
  85.       while ( (c=getc(stdin)) != '\n');
  86.       if (cstep < SMALL)
  87.          fprintf(stdout,"Error : Step size %10f is too small\n",cstep);
  88.       else {
  89.          nsteps = (zmax-zmin)/cstep;
  90.          if (nsteps < 1)
  91.             fprintf(stdout,"   Error : Too few contours %d\n",nsteps);
  92.          else if (nsteps > 50 && nsteps < 100) {
  93.             fprintf(stdout,"   Warning:%d contour steps will be made\n",nsteps);
  94.             fprintf(stdout,"   Type return (or n/q) to continue :");
  95.             if ( (c=getc(stdin)) == '\n' || c == 'y' || c == 'Y' )
  96.                FOUND = TRUE;
  97.             else if (c == 'q')
  98.                exit(1);
  99.             if (c != '\n') while ( (c=getc(stdin)) != '\n');
  100.          } else if (nsteps >= 100) 
  101.             fprintf(stdout,"   Error : Too many contours %d\n",nsteps);
  102.          else
  103.             FOUND = TRUE;
  104.       }
  105.    }
  106.    fprintf(stdout,"   Contour step size         = %f\n",cstep);
  107.    fprintf(stdout,"   Number of contours        = %d\n",nsteps);
  108.  
  109.    /* now find the contours */
  110.    level = (int)(zmin/cstep)*cstep + cstep;
  111.    while (level <= zmax+SMALL) {
  112.       find_contour(level);
  113.       level += cstep;
  114.    }
  115.  
  116.    /* now draw the plots */
  117.    if ( (strcmp("xterm",term)==0)  ||
  118.         (strcmp("xterms",term)==0) || 
  119.         (strcmp("vt200",term)==0)  || 
  120.         (strcmp("vt220",term)==0)  || 
  121.         (strcmp("vs100s",term)==0) ) {
  122.       InitWindow(argc,argv);
  123.       RunOps();
  124.    }
  125.  
  126.    /* hp2648 procedures */
  127.    if ((strcmp("hp2648",term)==0) || (strcmp("2648",term)==0)) {
  128.       axeshp();
  129.       plothp();
  130.       endgrhp();
  131.    }
  132.  
  133.    /* Postscript plot */
  134.    if (postscript == ON) {
  135.       if ((ips = fopen(PSFILE,"w")) == NULL) {
  136.          fprintf(stderr,"cat: can't open %s\n",PSFILE);
  137.          exit(1);
  138.       }
  139.       initps(ips);
  140.       axesps(ips);
  141.       plotps(ips);
  142.       endgrps(ips);
  143.    }
  144.  
  145.    /* Get postscript printout */
  146.    if (postscript == ON && printplot == ON) {
  147.       fprintf(stdout,"   Type return (or n/q) to ");
  148.       fprintf(stdout,"send the plot to the %s printer :",printer);
  149.       if ( (c=getc(stdin)) == '\n' || c == 'y' || c == 'Y' ) {
  150.          sprintf(command,"%s %s %s","lpr",printer,PSFILE);
  151.          system(command);
  152.          fprintf(stdout,"   Postscript file has been ");
  153.          fprintf(stdout,"sent to the %s printer\n",printer);
  154.       } else if (c == 'q')
  155.          exit(1);
  156.       if (c != '\n') while ( (c=getc(stdin)) != '\n');
  157.    }
  158. }
  159.  
  160. /* find all the segments that are formed when triangle intersects z=level */
  161. find_contour(level)
  162. double level;
  163. {
  164.    void    delete_segm();
  165.    plotptr insert_plot();
  166.    nodeptr insert_tailnode();
  167.    extern  triaptr tria_listhead;
  168.    extern  segmptr segm_listhead;
  169.    triaptr T;
  170.    segmptr S;
  171.    plotptr P;
  172.    data    PT1,PT2;
  173.  
  174.    for (T=tria_listhead; T!=NULL; T=T->next)
  175.       find_intsct(level,T);
  176.    
  177.    /* now sort out these segments */
  178.    /* start with the segmhead */
  179.    while ((S=segm_listhead)!=NULL) {
  180.       P = insert_plot(level);
  181.       PT1 = S->pt1;
  182.       PT2 = S->pt2;
  183.       insert_tailnode(&PT1,P);
  184.       insert_tailnode(&PT2,P);
  185.       delete_segm(S);
  186.       add_to_tailplot(&PT2,P);
  187.       add_to_headplot(&PT1,P);
  188.    }
  189. }
  190.  
  191. /* add to the tail of the plot list */
  192. add_to_tailplot(pt,P)
  193. data *pt;
  194. plotptr P;
  195. {
  196.    void    delete_segm();
  197.    nodeptr insert_tailnode();
  198.    segmptr matching_segm();
  199.    segmptr S;
  200.  
  201.    /* add the further node at every loop iteration */
  202.    while ((S=matching_segm(pt))!=NULL) {
  203.       if (longline(&(S->pt1),pt) ) {
  204.          insert_tailnode(&(S->pt1),P);
  205.          *pt = S->pt1;
  206.       } else {
  207.          insert_tailnode(&(S->pt2),P);
  208.          *pt = S->pt2;
  209.       }
  210.       delete_segm(S);
  211.    }
  212.    return;
  213. }
  214.  
  215. /* add to the head of the plot list */
  216. add_to_headplot(pt,P)
  217. data *pt;
  218. plotptr P;
  219. {
  220.    void    delete_segm();
  221.    nodeptr insert_headnode();
  222.    segmptr matching_segm();
  223.    segmptr S;
  224.  
  225.    /* add the further node at every loop iteration */
  226.    while ((S=matching_segm(pt))!=NULL) {
  227.       if (longline(&(S->pt1),pt) ) {
  228.          insert_headnode(&(S->pt1),P);
  229.          *pt = S->pt1;
  230.       } else {
  231.          insert_headnode(&(S->pt2),P);
  232.          *pt = S->pt2;
  233.       }
  234.       delete_segm(S);
  235.    }
  236.    return;
  237. }
  238.  
  239. /* Return segment containing data point that matches a given point */
  240. segmptr matching_segm(pt)
  241. data *pt;
  242. {
  243.    int     longline();
  244.    extern  segmptr segm_listhead;
  245.    segmptr S,SF;
  246.    int     FOUND = FALSE;
  247.  
  248.    /* check lengths of lines */
  249.    for (S=segm_listhead; S!=NULL && !FOUND; S=S->next) {
  250.       /* check the first point */
  251.       if (fabs(S->pt1.x - pt->x) < SMALL) 
  252.          FOUND = !(longline(&(S->pt1),pt));
  253.       if (!FOUND) {
  254.          /* now check the second point */
  255.          if (fabs(S->pt2.x - pt->x) < SMALL) 
  256.             FOUND = !(longline(&(S->pt2),pt));
  257.       }
  258.       if (FOUND) SF = S;
  259.    }
  260.    
  261.    /* return the segment */
  262.    if (!FOUND) SF = NULL;
  263.    return(SF);
  264. }
  265.  
  266. /* find the segment of intersection between a triangle and a plane */
  267. find_intsct(level,T)
  268. double level;
  269. triaptr T;
  270. {
  271.    int     longline();
  272.    int     plane_intsct_tria();
  273.    int     plane_intsct_line();
  274.    segmptr insert_segm();
  275.    data    intpts[3],pta;
  276.    int     i=0;
  277.    segmptr newptr;
  278.  
  279.    if (plane_intsct_tria(level,T)) {
  280.       if (plane_intsct_line(level,T->pt1,T->pt2,&pta)) intpts[i++] = pta;
  281.       if (plane_intsct_line(level,T->pt2,T->pt3,&pta)) intpts[i++] = pta;
  282.       if (plane_intsct_line(level,T->pt3,T->pt1,&pta)) intpts[i++] = pta;
  283.       if (i==3) {
  284.          /* only way is to have plane intersect tria at one of the nodes */
  285.          if (longline(&(intpts[0]),&(intpts[1])))
  286.             insert_segm(&(intpts[0]),&(intpts[1]));
  287.          else if (longline(&(intpts[1]),&(intpts[2])))
  288.             insert_segm(&(intpts[1]),&(intpts[2]));
  289.          else if (longline(&(intpts[2]),&(intpts[0])))
  290.             insert_segm(&(intpts[2]),&(intpts[0]));
  291.          else
  292.             fprintf(stderr,"Warning: Zero segment length\n");
  293.       } else if (i!=2)
  294.          fprintf(stderr,"Warning: Nonstandard No of intersections = %5d\n",i);
  295.       else if (longline(&(intpts[0]),&(intpts[1])))
  296.          newptr = insert_segm(&(intpts[0]),&(intpts[1]));
  297.       else
  298.          fprintf(stderr,"Warning: Zero segment length\n");
  299.    }
  300. }
  301.  
  302. /* find out if a line segment is long or short */
  303. int longline(pta,ptb)
  304. data *pta,*ptb;
  305. {
  306.    int    longln = TRUE;
  307.    double dl,dx,dy,dz;
  308.  
  309.    dx = pta->x - ptb->x;
  310.    dy = pta->y - ptb->y;
  311.    dz = pta->z - ptb->z;
  312.    dl = dx*dx + dy*dy + dz*dz;
  313.    if (dl < SMALLER) longln = FALSE;
  314.    return(longln);
  315. }
  316.  
  317. /* find out the intersection point between a line and a z-plane */
  318. int plane_intsct_line(z,pt1,pt2,pt3)
  319. double z;
  320. data pt1,pt2,*pt3;
  321. {
  322.    int    intsct = TRUE;
  323.    double t;
  324.  
  325.    if ((pt1.z > z && pt2.z > z) || 
  326.        (pt1.z < z && pt2.z < z) ||
  327.        (fabs(pt1.z-pt2.z) < SMALL) )
  328.       intsct = FALSE;
  329.    if (intsct) {
  330.       t = (z - pt1.z)/(pt2.z-pt1.z);
  331.       pt3->x = pt1.x + (pt2.x-pt1.x)*t;
  332.       pt3->y = pt1.y + (pt2.y-pt1.y)*t;
  333.       pt3->z = z;
  334.    }
  335.    return(intsct);
  336. }
  337.  
  338. /* find out if the triangle intersects the plane */
  339. int plane_intsct_tria(z,T)
  340. double z;
  341. triaptr T;
  342. {
  343.    int    intsct = TRUE;
  344.    double z1,z2,z3;
  345.  
  346.    z1 = T->pt1.z;
  347.    z2 = T->pt2.z;
  348.    z3 = T->pt3.z;
  349.    if ((z1>z && z2>z && z3>z) || (z1<z && z2<z && z3<z))
  350.       intsct = FALSE;
  351.    return(intsct);
  352. }
  353.  
  354.