home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_300 / 300_01 / contour.c < prev    next >
C/C++ Source or Header  |  1989-12-30  |  24KB  |  457 lines
  1. /*   HEADER:   CUG300;
  2.       TITLE:   Contour ploting demonstration;
  3.        DATE:   4/20/1989;
  4. DESCRIPTION:   "Plot contours using mat.lib functions";
  5.     VERSION:   2.04;
  6.    FILENAME:   CONTOUR.C;
  7.    SEE-ALSO:   MAT_V2D.H;
  8. */
  9.  
  10. /*===================================================================*/
  11. /* ROUTINES TO GENERATE CADD CONTOUR LINES FROM INPUT GRID VALUES    */
  12. /*===================================================================*/
  13. /*      0,1       1,1                                                */
  14. /*   m4.-----------.m3            ^ 0,jub                            */
  15. /*     | .       . |              |                                  */
  16. /*     |   .   .   |              |                                  */
  17. /*     |     .m0   |              |                                  */
  18. /*     |    .  .   |              |                                  */
  19. /*     |  .      . |              |                                  */
  20. /*   m1|___________|m2           0,0 --------------> iub,0           */
  21. /*     0,0        1,0                                                */
  22. /*                                                                   */
  23. /*   grid square layout            coordinate system layout          */
  24. /*   -------------------           ------------------------          */
  25. /*===================================================================*/
  26.  
  27. #include <stdio.h>
  28. #include "mat_v2d.h"
  29. #include <math.h>
  30. #include <stdlib.h>
  31. #include <conio.h>
  32.  
  33.  
  34. #define  FALSE       0
  35. #define  TRUE        1
  36. #define  NO_ROWS     24
  37. #define  NO_COLS     72
  38.  
  39. #define  xchk(x)     if(x<0)puts("x<0");if(x>NO_COLS)puts("x>NO_COLS")
  40. #define  ychk(y)     if(y<0)puts("y<0");if(y>NO_COLS)puts("y>NO_COLS")
  41.  
  42. void  box          (float x1,float y1,float x2,float y2 );
  43. void  dda          (float x1,float y1,float x2,float y2,char symbol);
  44. void  display      (void);
  45. void  erase        (void);
  46. void  GridSelect   (int iub, int jub, int nc, struct fmat *x,
  47.                     struct fmat *y,struct fmat *d, struct fmat *z);
  48. void ContourSelect (int i,int j, int nc,struct fmat *x,
  49.                           struct fmat *y,struct fmat *d,struct fmat *z,
  50.                           int iub, int jub);
  51. void drawit        (float x1, float x2, float y1, float y2,
  52.                     struct fmat *x, struct fmat *y,int iub,int jub );
  53. void set_graphics  (void);
  54. void file_output   (void);
  55.  
  56.  
  57. char  frame[NO_COLS][NO_ROWS];
  58.  
  59.  
  60. /*============================================================================*/
  61.  
  62. main()
  63.  
  64. /*============================================================================*/
  65.  
  66. {
  67. int i,i1,j,j1,k,iub,jub,nc,prmerr,num_rows,num_cols,token_size;
  68. struct   fmat *x,*y,*d,*z;
  69. struct   tmat *input_data;
  70. char     filename[20];
  71.  
  72. /*       Input arrays and indexes                                             */
  73.  
  74.     printf("Enter contour data filename => ");
  75.    scanf("%s",filename);
  76.    mtcnt(filename,&num_rows,&num_cols,&token_size);    /* Get input file size */
  77.    tdim(input_data,num_rows,num_cols,token_size);      /* Dimension input     */
  78.    mtget(filename,input_data);                         /* array               */
  79.  
  80.     iub = no_toks(input_data,0);                        /* Get indices from    */
  81.     jub = no_toks(input_data,1);                        /* the number of tokens*/
  82.     nc  = no_toks(input_data,2);                        /* stored in the       */
  83.                                                        /* appropriate input   */
  84.    fdim(d,iub,jub+1);                                  /* data line header    */
  85.    fdim(x,iub,DCLVCT);
  86.    fdim(y,jub,DCLVCT);                                 /* Dimension needed    */
  87.    fdim(z,nc,DCLVCT);                                  /* arrays              */
  88.  
  89.     for(i=0; i<iub; i++) {                              /* Store x grid loca-  */
  90.       fck(x,i,VCT);                                    /* tion.               */
  91.         f(x,i,VCT) = tf(input_data,0,i,iub);
  92.    }
  93.     for(j=0; j<jub; j++) {                              /* Store y grid loca-  */
  94.       fck(y,j,VCT);                                    /* tions               */
  95.         f(y,j,VCT) = tf(input_data,1,j,jub);
  96.    }
  97.    for(i=0; i<nc; i++) {                               /* Store contour       */
  98.       fck(z,i,VCT);                                    /* levels to be        */
  99.       f(z,i,VCT) = tf(input_data,2,i,nc);              /* plotted             */
  100.    }
  101.    for (j=3,k=jub-1; j<num_rows; j++,k--)  {
  102.       for (i=0; i<iub; i++) {                          /* Store f(x,y) for    */
  103.             fck(d,i,k);                                   /* each grid line      */
  104.             f(d,i,k) = tf(input_data,j,i,iub);            /* intersection        */
  105.       }
  106.    }
  107.    iub--;jub--;
  108.  
  109. /*       Check input parameters for validity                                  */
  110.  
  111.    prmerr = FALSE;                                     /* Test for no grid    */
  112.     if (iub<=0 || jub<=0)                               /* intersections.      */
  113.       prmerr = TRUE;                                   /*                     */
  114.    if (nc <=0)                                         /* No contours         */
  115.       prmerr = TRUE;                                   /* requested           */
  116.    for (k = 1; k <= (nc-1); k++)                       /*                     */
  117.    {                                                   /*                     */
  118.       if (f(z,k,VCT) <= f(z,k-1,VCT) )                 /* Duplicate contours  */
  119.          prmerr = TRUE;                                /*                     */
  120.    }                                                   /*                     */
  121.    if (prmerr)                                         /*                     */
  122.    {                                                   /*                     */
  123.       printf ("Error in input parameters");            /* Message if any      */
  124.       return;                                          /* errors found        */
  125.    }                                                   /*                     */
  126.    set_graphics();                                     /* Initialize graphic  */
  127.                                                        /* frame buffer        */
  128.    GridSelect (iub, jub, nc, x, y, d, z );             /* Plot contours       */
  129.                                                        /*                     */
  130.    display();                                          /* Display all         */
  131.                                                        /* contours plotted    */
  132.    file_output();                                      /*                     */
  133.                                                        /* Place contour in    */
  134.                                                        /* an ASCII file       */
  135. }                                                      /*                     */
  136.                                                                               
  137.  
  138. /*============================================================================*/
  139.  
  140. void GridSelect  (   int iub, int jub, int nc, struct fmat *x,
  141.                      struct fmat *y,struct fmat *d, struct fmat *z )
  142.  
  143. /*============================================================================*/
  144. {
  145. int i, j;
  146.  
  147.  
  148. /*    Scan the input array top down & left to right                           */
  149.  
  150.  
  151.  
  152.    for  (j = (jub-1); j >= 0; j--)
  153.    {
  154.       for (i = 0; i <= (iub-1); i++)
  155.       {
  156.         ContourSelect (i,j,nc,x,y,d,z,iub,jub);          /* Select contour  to  */
  157.       }                                                /* plot                */
  158.    }
  159. }
  160.  
  161. /*============================================================================*/
  162.  
  163. void ContourSelect ( int i,int j, int nc,struct fmat *x,
  164.                      struct fmat *y, struct fmat *d,
  165.                      struct fmat *z, int iub, int jub )
  166.  
  167. /*============================================================================*/
  168.  
  169. {
  170.  
  171. /*   Declare & initialize variables                                           */
  172.  
  173. float     h [5],    xh [5],     yh[5];
  174.  
  175. int      ish[5];
  176. int    caseval,k,i1,j1,m1, m2, m3;
  177. int            m;
  178. static int  im[4] =   { 0, 1, 1, 0 };
  179. static int  jm[4] =   { 0, 0, 1, 1 };
  180. static int  castab [3][3][3] =
  181.                          { { { 0, 0, 8 },
  182.                              { 0, 2, 5 },              /*   (k,i,j)       */
  183.                              { 7, 6, 9 }, },           /*   -- i ---->    */
  184.                                                        /*  | \            */  
  185.                                  { { 0, 0, 4 },        /*      k          */  
  186.                                    { 1, 0, 1 },        /*  j              */  
  187.                                    { 4, 0, 0 }, },     /*  |       \      */  
  188.                                                        /*  v         @    */  
  189.                                         { { 9, 6, 7 },                        
  190.                                           { 5, 2, 0 },                        
  191.                                           { 8, 0, 0 }  }  };                  
  192.                                                                               
  193. float dmin, dmax, x1, y1, x2, y2;                                             
  194.                                                                               
  195. /*   Find the lowest vertex                                                   */
  196.                                                                               
  197.    if ( f(d,i,j) < f(d,i,j+1) )
  198.       dmin = f(d,i,j);                                                        
  199.    else dmin = f(d,i,j+1);
  200.    if ( f(d,i+1,j)  < dmin )
  201.       dmin = f(d,i+1,j);                                                      
  202.    if ( f(d,i+1,j+1) < dmin )
  203.       dmin = f(d,i+1,j+1);                                                    
  204.                                                                               
  205. /*   Find highest vertex                                                      */
  206.                                                                               
  207.    if ( f(d,i,j) > f(d,i,j+1) )
  208.       dmax = f(d,i,j);
  209.       else dmax = f(d,i,j+1);                                                 
  210.    if ( f(d,i+1,j) > dmax)
  211.       dmax= f(d,i+1,j);
  212.    if (f(d,i+1,j+1) > dmax)
  213.        dmax = f(d,i+1,j+1);                                                   
  214.    if ( dmax < f(z,0,VCT) || dmin > f(z,nc-1,VCT) )
  215.        return;  /* No contours to be plotted in grid */
  216.  
  217. /*   Draw each contour in grid square                                         */
  218.  
  219.    for (k = 0; k <= (nc-1); k++)
  220.    {
  221.       if (!( f(z,k,VCT) < dmin || f(z,k,VCT) > dmax ))
  222.       {   /* Contour in grid */
  223.  
  224. /*   Determine line segment case                                              */
  225.  
  226.          for (m = 4; m >= 0; m--)
  227.          {
  228.             if ( m > 0 )
  229.             {
  230.                h[m] = f(d,(i+im[m-1]),(j+jm[m-1])) - f(z,k,VCT);
  231.                xh[m]= f(x,i+im[m-1],VCT);
  232.                yh[m]= f(y,j+jm[m-1],VCT);
  233.             }
  234.             if ( m==0 )
  235.             {
  236.                h[0] = (h[1]+h[2]+h[3]+h[4])/4;
  237.                xh[0] = ( f(x,i,VCT) + f(x,i+1,VCT) )/2;
  238.                yh[0] = ( f(y,j,VCT) + f(y,j+1,VCT) )/2;
  239.             }
  240.  
  241.             if      ( h[m] > 0 )  ish[m]= 2;
  242.             else if ( h[m] < 0 )  ish[m]= 0;
  243.             else                  ish[m]= 1;
  244.  
  245.          } /* end m loop */
  246.  
  247.          for (m = 1; m <= 4; m++)
  248.          {
  249.             m1 =m;
  250.             m2 = 0;
  251.             m3 = m+1;
  252.             if (m3 == 5) m3 = 1;
  253.             caseval = castab [ish[m1]] [ish[m2]] [ish[m3]];
  254.                                                                               
  255. /*    Scan each triangle in grid square                                       */
  256.                                                                               
  257.             switch (caseval)
  258.             {                                                                 
  259.                case 0 : break;      /* No intersection */
  260.                                                                               
  261.                case 1 : x1 =xh[m1]; /* Line between vertices m1 and m2        */
  262.                         y1=yh[m1];                                            
  263.                         x2=xh[m2];                                            
  264.                         y2=yh[m2];                                            
  265.                         break;                                                
  266.                case 2 : x1 =xh[m2]; /* Line between vertices m2 and m3        */
  267.                         y1=yh[m2];                                            
  268.                         x2=xh[m3];                                            
  269.                         y2=yh[m3];                                            
  270.                         break;                                                
  271.                case 3 : x1 =xh[m3];/* Line between vertices m3 and m1         */
  272.                         y1=yh[m3];                                            
  273.                         x2=xh[m1];                                            
  274.                         y2=yh[m1];                                              
  275.                         break;                                                
  276.                case 4 : x1=xh[m1]; /* Line between vertices m1 and side m2-m3 */
  277.                         y1=yh[m1];                                            
  278.                         x2=(h[m3]*xh[m2]-h[m2]*xh[m3])/(h[m3]-h[m2]);         
  279.                         y2=(h[m3]*yh[m2]-h[m2]*yh[m3])/(h[m3]-h[m2]);         
  280.                         break;                                                
  281.                case 5 : x1=xh[m2]; /* Line between vertices m2 and side m3-m1 */
  282.                         y1=yh[m2];                                            
  283.                         x2=(h[m1]*xh[m3]-h[m3]*xh[m1])/(h[m1]-h[m3]);         
  284.                         y2=(h[m1]*yh[m3]-h[m3]*yh[m1])/(h[m1]-h[m3]);         
  285.                         break;                                                
  286.                case 6 : x1=xh[m3]; /* Line between vertices m3 and side m1-m2 */
  287.                         y1=yh[m3];                                            
  288.                         x2=(h[m2]*xh[m1]-h[m1]*xh[m2])/(h[m2]-h[m1]);         
  289.                         y2=(h[m2]*yh[m1]-h[m1]*yh[m2])/(h[m2]-h[m1]);         
  290.                         break;                                                
  291.                                    /*  Line between sides m1-m2 & m2-m3       */
  292.                                                                               
  293.                case 7 : x1 =(h[m2]*xh[m1]-h[m1]*xh[m2])/(h[m2]-h[m1]);        
  294.                         y1 =(h[m2]*yh[m1]-h[m1]*yh[m2])/(h[m2]-h[m1]);        
  295.                         x2 =(h[m3]*xh[m2]-h[m2]*xh[m3])/(h[m3]-h[m2]);        
  296.                         y2 =(h[m3]*yh[m2]-h[m2]*yh[m3])/(h[m3]-h[m2]);        
  297.                         break;                                                
  298.                                    /*  Line between sides m2-m3 & m3-m1       */
  299.                                                                               
  300.                case 8 : x1 =(h[m3]*xh[m2]-h[m2]*xh[m3])/(h[m3]-h[m2]);        
  301.                         y1 =(h[m3]*yh[m2]-h[m2]*yh[m3])/(h[m3]-h[m2]);        
  302.                         x2 =(h[m1]*xh[m3]-h[m3]*xh[m1])/(h[m1]-h[m3]);        
  303.                         y2 =(h[m1]*yh[m3]-h[m3]*yh[m1])/(h[m1]-h[m3]);        
  304.                         break;                                                
  305.                                    /*  Line between sides m3-m1 & m1-m2       */
  306.                                                                               
  307.                case 9 : x1 =(h[m1]*xh[m3]-h[m3]*xh[m1])/(h[m1]-h[m3]);        
  308.                         y1 =(h[m1]*yh[m3]-h[m3]*yh[m1])/(h[m1]-h[m3]);        
  309.                         x2 =(h[m2]*xh[m1]-h[m1]*xh[m2])/(h[m2]-h[m1]);        
  310.                         y2 =(h[m2]*yh[m1]-h[m1]*yh[m2])/(h[m2]-h[m1]);        
  311.                         break;                                               
  312.  
  313.              default : printf("The defined case for vertex %d, contour %d\n",m,k);
  314.                     printf("grid %d,%d does not exist\n", k,i,j);
  315.             }                                                                 
  316.             if (caseval!=0)
  317.             {                                                                 
  318.                drawit(x1,y1,x2,y2,x,y,iub,jub);                               
  319.             }                                                                 
  320.          } /* Loop m */                                                       
  321.       } /* if */                                                              
  322.    } /* Loop k */                                                             
  323. }
  324.                                                                               
  325. /*============================================================================*/
  326.                                                                               
  327. void drawit (float x1, float y1, float x2, float y2,                          
  328.              struct fmat *x, struct fmat *y,int iub, int jub )
  329.                                                                               
  330. /*============================================================================*/
  331. {                                                                             
  332.    x1 = (NO_COLS-3)*x1/(f(x,iub,VCT))+1;               /* Draw the selected   */
  333.    y1 = (NO_ROWS-3)*y1/(f(y,jub,VCT))+1;               /* line segment in     */
  334.    x2 = (NO_COLS-3)*x2/(f(x,iub,VCT))+1;               /* the frame buffer    */
  335.    y2 = (NO_ROWS-3)*y2/(f(y,jub,VCT))+1;                                      
  336.                                                                               
  337.    dda(x1,y1,x2,y2,'+');                                                      
  338. }                                                                             
  339.                                                                               
  340. /*============================================================================*/
  341.                                                                               
  342. void set_graphics(void)
  343.                                                                               
  344. /*============================================================================*/
  345. {                                                                             
  346. erase();                                               /* Clear the fuffer    */
  347. box(0,0,NO_COLS-1,NO_ROWS-1);                          /* Draw box in frame   */
  348. }                                                      /* buffer              */
  349.                                                                               
  350. /*===========================================================================*/
  351.                                                                               
  352. void  box(float x1,float y1,float x2,float y2 )                               
  353.                                                                               
  354. /*===========================================================================*/
  355.                                                                               
  356. {                                                                             
  357.    xchk(x1);xchk(x2);ychk(y1);ychk(y2);                /* Draw box with ASCII */
  358.    frame[x1][y1]       ='*';                           /* characters in       */
  359.    dda(x1,y1+1, x1,y2-1,'|');                          /* the graphic frame   */
  360.    frame[x1][y2]       ='*';                                                  
  361.    dda(x1+1,y2, x2-1,y2,'-');                                                 
  362.    frame[x2][y2]       ='*';                                                  
  363.    dda(x2,y2-1, x2,y1+1,'|');                                                 
  364.    frame[x2][y1]       ='*';                                                  
  365.    dda(x2-1,y1, x1+1,y1,'-');                                                 
  366. }                                                                             
  367.                                                                               
  368. /*===========================================================================*/
  369.                                                                               
  370. void  dda (float x1,float y1,float x2,float y2,char symbol)                   
  371.                                                                               
  372. /*===========================================================================*/
  373. {                                                                             
  374. float    x,y,dx_step,dy_step;                           /* Draw a line in    */
  375. int steps,i,dx,dy;                                 /* the bufffer       */
  376.                                                                               
  377.    dx=x2-x1;                              /* Calc deltas */                   
  378.    dy=y2-y1;                                                                  
  379.    steps=max(abs(dx),abs(dy))+1;          /* Calc no steps */
  380.    dx_step=(float)dx/steps;               /* Step size */
  381.    dy_step=(float)dy/steps;
  382.    x=x1+.5;                               /* Offsets for rounding */
  383.    y=y1+.5;
  384.    for (i=0; i<steps; i++) {              /* Set pixels */
  385.       frame[floor(x)][floor(y)]=symbol;
  386.       x+=dx_step;
  387.       y+=dy_step;
  388.    }
  389.    frame[floor(x)][floor(y)]=symbol;
  390.                                           /* Last pixel */
  391. }
  392.  
  393. /*===========================================================================*/
  394.  
  395. void  display (void)
  396.  
  397. /*===========================================================================*/
  398. {
  399. int        x,y;                                        /* Display the frame  */
  400.                                                        /* buffer on the std  */
  401.    clrscr();                                           /* output device      */
  402.    for (y=NO_ROWS-1; y>=0; y--) {
  403.       for (x=0; x<NO_COLS; x++) {
  404.          putch(frame[x][y]);
  405.       }
  406.       putch('\n');
  407.       putch('\r');
  408.    }
  409. }
  410.                                                                               
  411. /*===========================================================================*/
  412.                                                                               
  413. void  erase (void)
  414.  
  415. /*===========================================================================*/
  416. {
  417. int   x,y;
  418.                                                        /* Clear the frame    */
  419.    for(y=NO_ROWS; y>0; y--) {                          /* buffer             */
  420.       for (x=1; x<=NO_COLS; x++) {
  421.          frame[x][y]='\x20';
  422.       }
  423.    }
  424. }
  425.  
  426. /*===========================================================================*/
  427.  
  428. void  file_output (void)
  429.  
  430. /*===========================================================================*/
  431. {
  432. char ch,filename[24];
  433. int   x,y;
  434. FILE  *FN;
  435.  
  436.  
  437.    printf("ASCII file copy desired ? (Y/N) => ");ch=getche();
  438.  
  439.    if (ch=='Y'||ch=='y') {                             /* Output the frame   */
  440.                                                        /* buffer to a text-  */
  441.       clrscr();                                        /* file               */
  442.       printf("Enter the output filename desired:");
  443.       scanf("%s",filename);
  444.  
  445.       FN=fopen(filename,"w");
  446.  
  447.       for (y=NO_ROWS-1; y>=0; y--) {
  448.          for (x=0; x<NO_COLS; x++) {
  449.             putc(frame[x][y],FN);
  450.          }
  451.       putc('\n',FN);
  452.       }
  453.       fclose(FN);
  454.    }
  455. }
  456.  
  457.