home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 293_02 / grad1.c < prev    next >
C/C++ Source or Header  |  1989-08-25  |  33KB  |  1,012 lines
  1. /******************** GRAD ******************************
  2.   grad.c
  3.   Daniel Geist - Mallinckrodt Institue of Radiology 1987.
  4.    This program reads a set of ct scans for one
  5.    patient. It outputs two sets of files containing predefined 
  6.    views of bone surfaces. One set contains distance values and
  7.    the other gradient values.
  8.  
  9.    Danny Geist & Mike Vannier
  10.  
  11.    Contact for further information:
  12.     Michael W. Vannier
  13.     Mallinckrodt Institute of Radiology
  14.     Washington University School of Medicine
  15.     510 South Kingshighway Blvd.
  16.     St. Louis, Mo. 63110 USA
  17.  
  18.    Ake Wallin has made the following changes - 1988
  19.  
  20.    -Command line input
  21.    -Threshold of gradient for differencing
  22.    -Views to choose: BO (bottom), TO (top), RL (right lateral),
  23.     LL (left lateral), RE (rear), FR (front), NO (none, just leave distance)
  24.    -Clipping planes: RL (right lateral), LL (lateral left), RE (rear),
  25.     FR (front)
  26.    -Scan of data in "correct" order, assuming:
  27.     object:
  28.     low z - feet, high z - head
  29.     low x - right, high x - left
  30.     low y - back, high y - front
  31.     image:
  32.     low x - left, high x - right
  33.     low y - down, high y - up
  34.             front
  35.     view BO : right      left
  36.             back
  37.             front
  38.     view TO : left       right
  39.             back
  40.             head
  41.     view RL : back      front
  42.             feet
  43.             head
  44.     view LL : front     back
  45.             feet
  46.             head
  47.     view RE : left      right
  48.             feet
  49.             head
  50.     view FR : right     left
  51.             feet
  52. *********************************************************/
  53. #include <stdio.h>
  54. #include <math.h>
  55. #include <ctype.h>
  56.  
  57. /* some global variables*/
  58. int FIRSTSLICE,LASTSLICE,THRESHOLD,NLINES;
  59. float ZOOM;
  60. char DR;
  61. int  buffer[2][6][256];          /* input buffer */
  62. float fxbuf[5][256],fybuf[5][256]; /* X,Y Views floating buffers */
  63.  
  64. float huge fzbuf[256][256];        /* Z view floating buffer */
  65. char fnamein[50]="ctbild.000";        /* input file name  */
  66. float GRAD_THRESHOLD = 10.0;
  67. unsigned int   views = 0xFFFF;
  68. int   clipx[2], clipy[2];      /* clipping planes */
  69. int   dispmode;                /* distance and/or gradient */
  70. int   image_or[2];            /* image oriention in x and y */
  71. int   object_or[3];           /* image oriention in x, y and z */
  72. int   header_blocks;          /* number of header blocks in CT's */
  73. char *usestr1 = "Usage: grad [file] [-z] [-f] [-l] [-t] [-d] [-n(g|d)] [-g]";
  74. char *usestr2 = "            [-h] [-v(bo|to|rl|ll|re|fr|no)] [-c(rl|ll|re|fr)]";
  75. char *usestr3 = "            [-i(x(r|l)|y(u|d))] [-o(x(r|l)|y(f|b)|z(h|f))]";
  76.  
  77. succ(i)   /* (i+1) modolus 3 */
  78. int i;
  79. {
  80.     return(i==2?0:i+1);
  81. }
  82.  
  83. prev(i)   /* (i-1) modolus 3 */
  84. int i;
  85. {
  86.     return(i==0?2:i-1);
  87. }
  88.  
  89. /* Set input file - add slice number to extension */
  90. setfilename(filenum)
  91. int filenum;
  92. {
  93.     int i;
  94.  
  95.     for (i = strlen(fnamein)-1; i > 0; i--)
  96.       if (fnamein[i] == '.') break;
  97.     if (i == 0) {
  98.       printf("CT file names inconsistend!\n");
  99.       exit(1);
  100.     }
  101.     fnamein[++i]=filenum/100+'0';
  102.     fnamein[++i]=(filenum%100)/10+'0';
  103.     fnamein[++i]='0'+filenum%10;
  104. }
  105.  
  106. /* interpolate from bottom line to top line n-1 lines  */
  107. interpolate(line,bot,top,n)
  108. int line,bot,top,n;
  109. {
  110.     int next,i,j,x,inc;
  111.     inc=top>bot?1:(-1);    /* interpolate backward or forwards ? */
  112.     next=bot+inc;
  113.     for(i=1,j=n-1;i<n;i++,j--){   /* do for each next line of interpolation */
  114.         for(x=0;x<256;x++) buffer[line][next][x]=
  115.             (buffer[line][bot][x]*j+buffer[line][top][x]*i)/n;
  116.         next+=inc;
  117.     }
  118. }
  119.  
  120. /* midpoint - fraction part of threshold transition distance */ 
  121. float midpoint(b,a)
  122. float b,a;
  123. {
  124.     return( (THRESHOLD-a) / (b-a) );
  125. }
  126.  
  127. /* get floating point distance values  */
  128. getdistances(xstart,xstop,xdir,
  129.              ystart,ystop,ydir,
  130.              start_slice,end_slice,zdir,
  131.              xbufxstart, xbufxdir, xbufystart, xbufydir,
  132.              ybufxstart, ybufxdir, ybufystart, ybufydir,
  133.              zbufxstart, zbufxdir, zbufystart, zbufydir,
  134.              pass)
  135. int xstart,xstop,xdir,ystart,ystop,ydir,start_slice,end_slice,zdir,pass;
  136. int  xbufxstart, xbufxdir, xbufystart, xbufydir;
  137. int  ybufxstart, ybufxdir, ybufystart, ybufydir;
  138. int  zbufxstart, zbufxdir, zbufystart, zbufydir;
  139. /* also used are global variables: fxbuf, fybuf, fzbuf, buffer,
  140.    ZOOM, DR, image_or, object_or. */
  141. {
  142.     int z,x,y,i,j,start,stop,inc,line,rzoom,inter;
  143.     float remain;
  144.     FILE *fxfloat,*fyfloat,*fzfloat,*fn[2];
  145.     char filename[13];
  146.     int  filegap;
  147.     int  xbfx, ybfx, zbfx, zbfy;
  148.  
  149.     NLINES=0;         /* number of output lines in X,Y view directions */
  150.     remain=0;         /* remainder of interpolation after roundoff */
  151.     rzoom=ZOOM+0.5;   /* rounded zoom factor */
  152.     /* X view floating output file */
  153.     sprintf(filename,"%c:xdis%d.dat",DR,pass);
  154.     if ((fxfloat=fopen(filename,"wb")) == NULL) {
  155.       printf("Error creating %s!\n", filename);
  156.       exit(1);
  157.     }
  158.     /* Y view floating output file */
  159.     sprintf(filename,"%c:ydis%d.dat",DR,pass);
  160.     if ((fyfloat=fopen(filename,"wb")) == NULL) {
  161.       printf("Error creating %s!\n", filename);
  162.       exit(1);
  163.     }
  164.     for(i=0;i<256;i++)
  165.         for(j=0;j<256;j++)fzbuf[i][j]=256; /* clear Z view buffer */
  166.     
  167.     for(z=start_slice; z!=end_slice; z += zdir){  /* For each Slice */
  168.       /* open next two slice files */
  169.       setfilename(z);
  170.       if ((fn[0]=fopen(fnamein,"rb")) == NULL) continue;
  171.     for (filegap = zdir; filegap + z != end_slice+zdir; filegap+=zdir) {
  172.           setfilename(z+filegap);
  173.           if ((fn[1]=fopen(fnamein,"rb")) != NULL) break;
  174.       }
  175.       if (fn[1] == NULL) continue;
  176.  
  177.       inter=rzoom; /* interpolation factor assumed rounded zoom */
  178.       /* correct interpolation factor according to floating remainder */
  179.       remain+=rzoom-ZOOM;
  180.       if(remain>=1){
  181.         inter-=1;
  182.         remain-=1;
  183.       }
  184.       else if(remain<=(-1)){
  185.         inter+=1;
  186.         remain+=1;
  187.       }
  188.  
  189.       line=0;               /* current input buffer line */
  190.       for(j=0;j<inter;j++)  /* clear X,Y floating buffers */
  191.         for(i=0;i<256;i++) fxbuf[j][i]=fybuf[j][i]=256;
  192.       /* set index for zbuf */
  193.       zbfy = zbufystart;
  194.       /* set index for xbuf */
  195.       xbfx = xbufxstart;
  196.       for(y=ystart; y!=ystop; y+=ydir){              /* For each line */
  197.     /* skip to line*/
  198.     for(i=0;i<2;i++)fseek(fn[i],(long)512*(y+header_blocks),SEEK_SET);
  199.         fread(buffer[line][0],1,512,fn[0]);   
  200.         fread(buffer[line][inter],1,512,fn[1]);
  201.         interpolate(line,0,inter,inter); /* interpolate in_between */
  202.     for(i=0;i<inter;i++){   /* For each interpolation line */
  203.           /* set index for zbuf */
  204.       zbfx = zbufxstart;
  205.           /* set index for ybuf */
  206.           ybfx = ybufxstart;
  207.       for(x=xstart; x!=xstop; x+=xdir) { /* For each Voxel value */
  208.             /* find threshold transition*/
  209.         if (buffer[line][i+1][x] >= THRESHOLD) {
  210.               /* if first transition in X direction get floating
  211.                  distance */
  212.               if(fxbuf[i][xbfx]==256.0) fxbuf[i][xbfx]=(x==xstart)?0:
  213.                                     (x-xstart)*xdir-1+
  214.                      midpoint((float)buffer[line][i+1][x],
  215.                             (float)buffer[line][i+1][x-xdir]);
  216.  
  217.               /* if first transition in Y direction get floating
  218.                    distance */
  219.               if(fybuf[i][ybfx]==256.0) fybuf[i][ybfx]=(y==ystart)?0:
  220.                     (y-ystart)*ydir-1+
  221.                      midpoint((float)buffer[line][i+1][x],
  222.                               (float)buffer[1-line][i+1][x]);
  223.  
  224.               /* if first transition in Z direction get floating
  225.                    distance */
  226.               if(fzbuf[zbfy][zbfx]==256.0) fzbuf[zbfy][zbfx]=
  227.                    (i==0) && (buffer[line][i][x]>=THRESHOLD) ?
  228.                     NLINES : NLINES+i+
  229.                     midpoint((float)buffer[line][i+1][x],
  230.                              (float)buffer[line][i][x]);
  231.             }
  232.             /* change index of ybuf and zbuf */
  233.             ybfx += ybufxdir;
  234.         zbfx += zbufxdir;
  235.         if ((ybfx < 0 || ybfx > 255) && x+xdir != xstop) {
  236.               printf("Error in index YBFX!\n");
  237.               exit(1);
  238.             }
  239.         if ((zbfx < 0 || zbfx > 255) && x+xdir != xstop) {
  240.               printf("Error in index ZBFX!\n");
  241.               exit(1);
  242.             }
  243.           }
  244.         }
  245.         line=1-line;  /* change current input buffer line */
  246.         /* update index of xbuf and zbuf */
  247.         xbfx += xbufxdir;
  248.         zbfy += zbufydir;
  249.     if ((zbfy < 0 || zbfy > 255) && y+ydir != ystop) {
  250.           printf("Error in index ZBFY!\n");
  251.           exit(1);
  252.         }
  253.     if ((xbfx < 0 || xbfx > 255) && y+ydir != ystop) {
  254.           printf("Error in index XBFX!\n");
  255.           exit(1);
  256.         }
  257.       }
  258.       NLINES+=inter;           /* increment number of output lines */
  259.  
  260.       fclose(fn[0]);           /* close slice files */
  261.       fclose(fn[1]);
  262.       /* write output lines to files */
  263.       if (fwrite(fxbuf,1,sizeof(float)*inter*256,fxfloat) <
  264.         inter*256*sizeof(float))
  265.     printf("Error writing to FXBUF!\n");
  266.       if (fwrite(fybuf,1,sizeof(float)*inter*256,fyfloat) <
  267.         inter*256*sizeof(float))
  268.     printf("Error writing to FYBUF!\n");
  269.       printf(" line %d\r", z);
  270.     }
  271.     printf("\n");
  272.     fclose(fxfloat); /* close X,Y floating files */
  273.     fclose(fyfloat);
  274.     /* write Z floating file */
  275.     sprintf(filename,"%c:zdis%d.dat",DR,pass);
  276.     fzfloat=fopen(filename,"wb");
  277.     for(i=0;i<256;i++)fwrite(fzbuf[i],sizeof(float),256,fzfloat);
  278.     fclose(fzfloat);
  279.  
  280.     /* should we reverse the order of the lines? */
  281.     if (xbufydir < 0) {
  282.       sprintf(filename,"%c:xdis%d.dat",DR,pass);
  283.       fxfloat=fopen(filename,"rb");
  284.       for (line = 0; line < NLINES; line++)
  285.         fread(fzbuf[line], 1, sizeof(float)*256, fxfloat);
  286.       fclose (fxfloat);
  287.       fxfloat=fopen(filename,"wb");
  288.       for (line = 1; line <= NLINES; line++)
  289.         fwrite(fzbuf[NLINES-line], sizeof(float), 256, fxfloat);
  290.       fclose (fxfloat);
  291.     }
  292.     if (ybufydir < 0) {
  293.       sprintf(filename,"%c:ydis%d.dat",DR,pass);
  294.       fyfloat=fopen(filename,"rb");
  295.       for (line = 0; line < NLINES; line++)
  296.         fread(fzbuf[line], 1, sizeof(float)*256, fyfloat);
  297.       fclose (fyfloat);
  298.       fyfloat=fopen(filename,"wb");
  299.       for (line = 1; line <= NLINES; line++)
  300.         fwrite(fzbuf[NLINES-line], sizeof(float), 256, fyfloat);
  301.       fclose (fyfloat);
  302.     }
  303.  
  304. }
  305.  
  306. unsigned char grad(y1,y2,z1,z2,y_factor,z_factor)
  307. float y1,y2,z1,z2;
  308. int y_factor,z_factor;
  309. {
  310.     float gx,gy,gz;
  311.     unsigned char gxint;
  312.     /* get z and y components of gradient */
  313.     gz=(z1-z2)/(float) z_factor;
  314.     gy=(y1-y2)/(float) y_factor;
  315.     /*compute gx - normalized x component of gradient */
  316.     gx=1/sqrt(1+gz*gz+gy*gy);
  317.     gxint=255*gx+0.5;      /*scale gx by 256 and round */
  318.     return(gxint);
  319. }
  320.  
  321. /* get gradient and depth shades for one image line */
  322. doline(lineg,lined,line,prevline,succline,z_factor,fg,fd)
  323. unsigned char lineg[],lined[];  /*output buffers */
  324. int line,prevline,succline;  /* input buffer configuration */
  325. int z_factor; /* for choosing forward, backward or central differences */
  326. FILE *fg,*fd;  /* output files */
  327. {
  328.     int i;
  329.     float y1, y2, z1, z2;
  330.     int factor;
  331.     /* do first pixel in line */
  332.     i = 0;
  333.     if(fxbuf[line][0]==256)lineg[0]=lined[0]=0;
  334.     else{
  335.         if (dispmode & 1)
  336.           lined[0]=255-fxbuf[line][0]; /*distance shade */
  337.         if (dispmode & 2)
  338.           lineg[0]=grad(fxbuf[line][0],fxbuf[line][1],fxbuf[prevline][0],
  339.                         fxbuf[succline][0],1,z_factor);
  340.     }
  341.     /* do rest of pixels inside line */
  342.     for(i=1;i<255;i++)
  343.       if(fxbuf[line][i]==256)lineg[i]=lined[i]=0;
  344.     else{
  345.         if (dispmode & 1)
  346.           lined[i]=255-fxbuf[line][i]; /*distance shade */
  347.         if (dispmode & 2) {
  348.           y1 = fxbuf[line][i-1];
  349.           y2 = fxbuf[line][i+1];
  350.           factor = 2;
  351.           if (fabs(y1 - y2) > GRAD_THRESHOLD) {
  352.             if (fabs(y1 - fxbuf[line][i]) < GRAD_THRESHOLD/2.) {
  353.               y2 = fxbuf[line][i];
  354.               factor = 1;
  355.             }
  356.             else if (fabs(fxbuf[line][i] - y2) < GRAD_THRESHOLD/2.) {
  357.               y1 = fxbuf[line][i];
  358.               factor = 1;
  359.             }
  360.           }
  361.           z1 = fxbuf[prevline][i];
  362.           z2 = fxbuf[succline][i];
  363.           if (z_factor == 2 && fabs(z1 - z2) > GRAD_THRESHOLD) {
  364.             if (fabs(z1 - fxbuf[line][i]) < GRAD_THRESHOLD/2.) {
  365.               z2 = fxbuf[line][i];
  366.               z_factor = 1;
  367.             }
  368.             else if (fabs(fxbuf[line][i] - z2) < GRAD_THRESHOLD/2.) {
  369.               z1 = fxbuf[line][i];
  370.               z_factor = 1;
  371.             }
  372.           }
  373.  
  374. /*        lineg[i]=grad(fxbuf[line][i-1],fxbuf[line][i+1],fxbuf[prevline][i],
  375.         fxbuf[succline][i],2,z_factor); */
  376.           lineg[i]=grad(y1, y2, z1, z2, factor, z_factor);
  377.         }
  378.     }
  379.     /* do last pixel in line */
  380.     if(fxbuf[line][255]==256) lineg[255]=lined[255]=0;
  381.     else{
  382.         if (dispmode & 1)
  383.           lined[255]=255-fxbuf[line][255]; /*distance shade */
  384.         if (dispmode & 2)
  385.           lineg[255]=grad(fxbuf[line][255],fxbuf[line][254],
  386.                           fxbuf[prevline][255],
  387.                           fxbuf[succline][255],1,z_factor);
  388.     }
  389.     if (dispmode & 2)
  390.       fwrite(lineg,1,256,fg); /* write to output files */
  391.     if (dispmode & 1)
  392.       fwrite(lined,1,256,fd);
  393. }
  394.  
  395. /* create an gradient and distance shaded view */
  396. doviews(namedis,nameg,named,nlines)
  397. char *namedis,*nameg,*named; /* floating file , gradiet and distance files*/
  398. int nlines;                  /* number of lines in image */
  399. {
  400.     FILE *fg,*fd,*ffloat;
  401.     int z,i,j,k,midline;
  402.     char lined[256],lineg[256]; /* gradient and distance value buffers */
  403.     midline=1;                  /* middle line in input buffer */
  404.     /* open output and input files */
  405.     fd = fg = NULL;
  406.     if (dispmode & 1)
  407.       fd=fopen(named,"wb");
  408.     if (dispmode & 2)
  409.       fg=fopen(nameg,"wb");
  410.     ffloat=fopen(namedis,"rb");
  411.     fread(fxbuf,1,3*256*sizeof(float),ffloat); /* read first three floating lines */
  412.     /* do first line */
  413.     doline(lineg,lined,0,0,1,1,fg,fd);
  414.     /* do rest of lines */
  415.     for(z=0;z<(nlines-2);z++){      /*for each inside line */
  416.         doline(lineg,lined,midline,prev(midline),succ(midline),2,fg,fd);
  417.         fread(fxbuf[prev(midline)],1,256*sizeof(float),ffloat); /*read next floating line */
  418.         midline=succ(midline);
  419.     }
  420.     /* do last line */
  421.     doline(lineg,lined,midline,prev(midline),midline,1,fg,fd);
  422.     /* close all files */
  423.     if (fg != NULL) fclose(fg);
  424.     if (fd != NULL) fclose(fd);
  425.     fclose(ffloat);
  426. }
  427.  
  428. /**********************************************************/
  429. /**** MAIN ***** MAIN ***** MAIN ***** MAIN ***** MAIN ****/
  430. /**********************************************************/
  431. /* Usage: grad [filename] [-f] [-l] [-z] [-t] [-d] [-n(d|g)] [-g]
  432.               [-h] [-v(bo|to|rl|ll|re|fr)] [-c(rl|ll|re|fr)]
  433.               [-i(x(r|l)|y(u|d))] [-o(x(r|l)|y(f|b)|z(h|f))]
  434.  */
  435. main(argc, argv)
  436. int argc;
  437. char *argv[];
  438. {   char  filename[13];
  439.     FILE *par;
  440.     FILE *fp;
  441.     int   i, n;
  442.     char  cmd[80];
  443.     int   scanxdir[3], scanydir[3], scanzdir[3];  /* scanning directions */
  444.     int   imagex[4], imagey[4], imagez[4];        /* image/scan direction */
  445.     int   temp_nlines;                            /* temporary line numbers */
  446.  
  447.     /* set clipping planes */
  448.     clipx[0] = clipy[0] = 0;
  449.     clipx[1] = clipy[1] = 256;
  450.     /* set image orientation to left->right, down->up */
  451.     image_or[0] = image_or[1] = 1;
  452.     /* set object orientation to right->left, back->front, feet->head */
  453.     object_or[0] = object_or[1] = object_or[2] = 1;
  454.     /* get parameters from the command line */
  455.     /* set other defaults */
  456.     ZOOM = 1.0;
  457.     FIRSTSLICE = 1;
  458.     LASTSLICE = 256;
  459.     THRESHOLD = 175;
  460.     DR = 'D';
  461.     dispmode = 3;
  462.     header_blocks = 1;
  463.     if (argc > 1) {
  464.       /* get parameters */
  465.       for (argc--, argv++; argc > 0; argc--, argv++) {
  466.         if (**argv == '-' || **argv == '/') {
  467.           switch (tolower(*(*argv+1))) {
  468.             case 'z':
  469.               ZOOM = atof(*argv+2);
  470.               break;
  471.             case 'f':
  472.               FIRSTSLICE = atoi(*argv+2);
  473.               break;
  474.             case 'l':
  475.               LASTSLICE = atoi(*argv+2);
  476.               break;
  477.             case 't':
  478.               THRESHOLD = atoi(*argv+2);
  479.               break;
  480.             case 'g':
  481.               GRAD_THRESHOLD = atof(*argv+2);
  482.               break;
  483.         case 'd':
  484.           DR = *(*argv+2);
  485.           break;
  486.         case 'h':
  487.           header_blocks = atoi(*argv+2);
  488.           break;
  489.             case 'n':
  490.               switch (tolower(*(*argv+2))) {
  491.                 case 'd':
  492.                   dispmode &= 0xfe;
  493.                   break;
  494.                 case 'g':
  495.                   dispmode &= 0xfd;
  496.                   break;
  497.         default:
  498.           usage();
  499.               }
  500.               break;
  501.         case 'v':
  502.           /* viewing directions */
  503.           if (views == 0xFFFF) views = 0;
  504.           if (tolower(*(*argv+2)) == 'b' &&
  505.           tolower(*(*argv+3)) == 'o')
  506.         views |= 1;
  507.           else if (tolower(*(*argv+2)) == 't' &&
  508.                tolower(*(*argv+3)) == 'o')
  509.         views |= 1 << 1;
  510.           else if (tolower(*(*argv+2)) == 'r' &&
  511.                tolower(*(*argv+3)) == 'l')
  512.         views |= 1 << 2;
  513.           else if (tolower(*(*argv+2)) == 'l' &&
  514.                tolower(*(*argv+3)) == 'l')
  515.         views |= 1 << 3;
  516.           else if (tolower(*(*argv+2)) == 'r' &&
  517.                tolower(*(*argv+3)) == 'e')
  518.         views |= 1 << 4;
  519.           else if (tolower(*(*argv+2)) == 'f' &&
  520.                tolower(*(*argv+3)) == 'r')
  521.         views |= 1 << 5;
  522.           else if (tolower(*(*argv+2)) == 'n' &&
  523.                tolower(*(*argv+3)) == 'o')
  524.             views = 0;
  525.           else
  526.         usage();
  527.           break;
  528.         case 'c':
  529.           /* clipping planes */
  530.           if (tolower(*(*argv+2)) == 'r' &&
  531.           tolower(*(*argv+3)) == 'l')
  532.         clipx[0] = atoi(*argv + 4);
  533.           else if (tolower(*(*argv+2)) == 'l' &&
  534.                tolower(*(*argv+3)) == 'l')
  535.         clipx[1] = atoi(*argv + 4);
  536.           else if (tolower(*(*argv+2)) == 'r' &&
  537.                tolower(*(*argv+3)) == 'e')
  538.         clipy[0] = atoi(*argv + 4);
  539.           else if (tolower(*(*argv+2)) == 'f' &&
  540.                tolower(*(*argv+3)) == 'r')
  541.         clipy[1] = atoi(*argv + 4);
  542.           else
  543.         usage();
  544.           break;
  545.         case 'i':
  546.           /* image orientation */
  547.           if (tolower(*(*argv+2)) == 'x' &&
  548.           tolower(*(*argv+3)) == 'r')
  549.         image_or[0] = 1;
  550.           else if (tolower(*(*argv+2)) == 'x' &&
  551.           tolower(*(*argv+3)) == 'l')
  552.         image_or[0] = -1;
  553.           else if (tolower(*(*argv+2)) == 'y' &&
  554.                tolower(*(*argv+3)) == 'u')
  555.         image_or[1] = 1;
  556.           else if (tolower(*(*argv+2)) == 'y' &&
  557.                tolower(*(*argv+3)) == 'd')
  558.         image_or[1] = -1;
  559.           else
  560.         usage();
  561.           break;
  562.         case 'o':
  563.           /* object orientation */
  564.           if (tolower(*(*argv+2)) == 'x' &&
  565.           tolower(*(*argv+3)) == 'r')
  566.         object_or[0] = -1;
  567.           else if (tolower(*(*argv+2)) == 'x' &&
  568.           tolower(*(*argv+3)) == 'l')
  569.         object_or[0] = 1;
  570.           else if (tolower(*(*argv+2)) == 'y' &&
  571.                tolower(*(*argv+3)) == 'f')
  572.         object_or[1] = 1;
  573.           else if (tolower(*(*argv+2)) == 'y' &&
  574.                tolower(*(*argv+3)) == 'b')
  575.         object_or[1] = -1;
  576.           else if (tolower(*(*argv+2)) == 'z' &&
  577.                tolower(*(*argv+3)) == 'h')
  578.         object_or[2] = 1;
  579.           else if (tolower(*(*argv+2)) == 'z' &&
  580.                tolower(*(*argv+3)) == 'f')
  581.         object_or[2] = -1;
  582.           else
  583.         usage();
  584.           break;
  585.         default:
  586.           usage();
  587.           }
  588.         }
  589.         else {
  590.           strncpy(fnamein, *argv, 45);
  591.       strcat(fnamein, ".000");
  592.         }
  593.       }
  594.     }
  595.     else {
  596.       cmd[0] = 1;
  597.       while(cmd[0]) {
  598.     /* display defaults */
  599.     printf("1 : filename          - %s\n", fnamein);
  600.     printf("2 : first slice       - %d, last slice - %d\n", FIRSTSLICE,
  601.                                 LASTSLICE);
  602.     printf("3 : zoom              - %f\n", ZOOM);
  603.     printf("4 : dens. threshold   - %d\n", THRESHOLD);
  604.     printf("5 : temp. drive       - %c\n", toupper(DR));
  605.     printf("6 : display           - (%s) (%s)\n",
  606.            (dispmode & 1) ? "distance" : "",
  607.            (dispmode & 2) ? "gradient" : "");
  608.     printf("7 : grad. threshold   - %f\n", GRAD_THRESHOLD);
  609.     printf("8 : header blocks     - %d\n", header_blocks);
  610.     printf("9 : views             - %s %s %s %s %s %s %s\n",
  611.            (views & 1) ? "BO" : "",
  612.            (views & 1<<1) ? "TO" : "",
  613.            (views & 1<<2) ? "RL" : "",
  614.            (views & 1<<3) ? "LL" : "",
  615.            (views & 1<<4) ? "RE" : "",
  616.            (views & 1<<5) ? "FR" : "",
  617.            (views == 0) ? "NO" : "");
  618.     printf("10: clipping planes   - RL %d LL %d RE %d FR %d\n",
  619.            clipx[0], clipx[1], clipy[0], clipy[1]);
  620.     printf("11: image orientation - X %d Y %d\n",
  621.            image_or[0], image_or[1]);
  622.     printf("12: objct orientation - X %d Y %d Z %d\n",
  623.            object_or[0], object_or[1], object_or[2]);
  624.     printf("\n Enter number to change any item, or ENTER to start\n");
  625.     if (gets(cmd) == NULL) exit(0);
  626.     if (cmd[0]) {
  627.       for (i = 0; !isdigit(cmd[i]) && cmd[i] != '\0'; i++) ;
  628.       if (isdigit(cmd[i])) {
  629.         sscanf(&cmd[i], "%d", &i);
  630.         switch (i) {
  631.           case 1:
  632.         printf("Enter file name : ");
  633.         if (gets(cmd) == NULL) exit(0);
  634.         sscanf(cmd, "%s", fnamein);
  635.         strcat(fnamein, ".000");
  636.         break;
  637.           case 2:
  638.         printf("Enter first scan : ");
  639.         if (gets(cmd) == NULL) exit(0);
  640.         sscanf(cmd, "%d", &FIRSTSLICE);
  641.         printf("Enter last scan : ");
  642.         if (gets(cmd) == NULL) exit(0);
  643.         sscanf(cmd, "%d", &LASTSLICE);
  644.         break;
  645.           case 3:
  646.         printf("Enter zoom factor : ");
  647.         if (gets(cmd) == NULL) exit(0);
  648.         sscanf(cmd, "%f", &ZOOM);
  649.         break;
  650.           case 4:
  651.         printf("Enter threshold : ");
  652.         if (gets(cmd) == NULL) exit(0);
  653.         sscanf(cmd, "%d", &THRESHOLD);
  654.         break;
  655.           case 5:
  656.         printf("Enter temporary drive : ");
  657.         if (gets(cmd) == NULL) exit(0);
  658.         for (i = 0; !isalpha(DR = cmd[i]) && DR != '\0'; i++) ;
  659.         if (!isalpha(DR)) DR = 'd';
  660.         break;
  661.           case 6:
  662.         printf("1 to turn distance %s, and 2 gradient %s\n",
  663.                (dispmode & 1) ? "off" : "on",
  664.                (dispmode & 2) ? "off" : "on");
  665.         if (gets(cmd) == NULL) exit(0);
  666.         sscanf(cmd, "%d", &i);
  667.         dispmode ^= 1 << --i;
  668.         break;
  669.           case 7:
  670.         printf("Enter gradient threshold : ");
  671.         if (gets(cmd) == NULL) exit(0);
  672.         sscanf(cmd, "%f", &GRAD_THRESHOLD);
  673.         break;
  674.           case 8:
  675.         printf("Enter number of header blocks : ");
  676.         if (gets(cmd) == NULL) exit(0);
  677.         sscanf(cmd, "%d", &header_blocks);
  678.         break;
  679.           case 9:
  680.         printf("1 BO %s, 2 TO %s, 3 RL %s, 4 LL %s, 5 RE %s, 6 FR %s, 7 NO %s ",
  681.                (views & 1) ? "off" : "on",
  682.                (views & 1<<1) ? "off" : "on",
  683.                (views & 1<<2) ? "off" : "on",
  684.                (views & 1<<3) ? "off" : "on",
  685.                (views & 1<<4) ? "off" : "on",
  686.                (views & 1<<5) ? "off" : "on",
  687.                (views == 0) ? "off" : "on");
  688.         if (gets(cmd) == NULL) exit(0);
  689.         sscanf(cmd, "%d", &i);
  690.         if (i == 7)
  691.           views == 0;
  692.         else
  693.           views ^= 1 << --i;
  694.         break;
  695.           case 10:
  696.         printf("Enter for clipping: 1 RL, 2 LL, 3 RE, 4 FR : ");
  697.         if (gets(cmd) == NULL) exit(0);
  698.         sscanf(cmd, "%d", &i);
  699.         if (i > 0 && i < 5) {
  700.           printf("Enter clipping : ");
  701.           if (gets(cmd) == NULL) exit(0);
  702.           sscanf(cmd, "%d", &n);
  703.           switch (i) {
  704.             case 1:
  705.               clipx[0] = n;
  706.               break;
  707.             case 2:
  708.               clipx[1] = n;
  709.               break;
  710.             case 3:
  711.               clipy[0] = n;
  712.               break;
  713.             case 4:
  714.               clipy[1] = n;
  715.               break;
  716.           }
  717.         }
  718.         break;
  719.           case 11:
  720.         printf("Enter for image orienation: 1 X, 2 Y : ");
  721.         if (gets(cmd) == NULL) exit(0);
  722.         sscanf(cmd, "%d", &i);
  723.         if (i > 0 && i < 3)
  724.           image_or[i-1] *= -1;
  725.         break;
  726.           case 12:
  727.         printf("Enter for object orienation: 1 X, 2 Y, 3 Z : ");
  728.         if (gets(cmd) == NULL) exit(0);
  729.         sscanf(cmd, "%d", &i);
  730.         if (i > 0 && i < 4)
  731.           object_or[i-1] *= -1;
  732.         break;
  733.           default:
  734.         printf("Unvalid command!\n");
  735.         break;
  736.         }
  737.       }
  738.  
  739.     }
  740.       }
  741.     }
  742.     THRESHOLD+=1024;
  743.  
  744.   /* test parameters */
  745.   if (ZOOM > 5.0 || ZOOM < 1.0) {
  746.     printf("ZOOM must not be smaller than 1 or larger than 5!\n");
  747.     exit(1);
  748.   }
  749.   if (FIRSTSLICE < 1 || LASTSLICE <= FIRSTSLICE) {
  750.     printf("FIRSTSLICE and LASTSLICE are not within valid range!\n");
  751.     exit(1);
  752.   }
  753.   for (i = FIRSTSLICE; i <= LASTSLICE; i++) {
  754.     setfilename(i);
  755.     if ((fp = fopen(fnamein, "rb")) == NULL) continue;
  756.     fclose(fp);
  757.     break;
  758.   }
  759.   if (i > LASTSLICE) {
  760.     setfilename(FIRSTSLICE);
  761.     printf("No files with a name between %s and", fnamein);
  762.     setfilename(LASTSLICE);
  763.     printf(" %s!\n", fnamein);
  764.     exit(1);
  765.   }
  766.   else
  767.     FIRSTSLICE = i;
  768.   for (i = LASTSLICE; i >= FIRSTSLICE; i--) {
  769.     setfilename(i);
  770.     if ((fp = fopen(fnamein, "rb")) == NULL) continue;
  771.     fclose(fp);
  772.     break;
  773.   }
  774.   LASTSLICE = i;
  775.   if (LASTSLICE <= FIRSTSLICE) {
  776.     printf("FIRSTSLICE and LASTSLICE are not within valid range!\n");
  777.     exit(1);
  778.   }
  779.   if (GRAD_THRESHOLD < 1.0) {
  780.     printf("GRAD_THRESHOLD is too small!\n");
  781.     exit(1);
  782.   }
  783.   if (!isalpha(DR)) {
  784.     printf("Temporary drive has to be a letter!\n");
  785.     exit(1);
  786.   }
  787.   if (header_blocks < 0) {
  788.     printf("Number of header blocks has to be positive!\n");
  789.     exit(1);
  790.   }
  791.   if (clipx[0] < 0 || clipx[1] > 256 || clipx[1] <= clipx[0]) {
  792.     printf("Lateral clipping planes are not within valid range!\n");
  793.   }
  794.   if (clipy[0] < 0 || clipy[1] > 256 || clipy[1] <= clipy[0]) {
  795.     printf("Frontal clipping planes are not within valid range!\n");
  796.   }
  797.  
  798.   /* temporary number of lines in output is used in finding
  799.      traversal of slices */
  800.   temp_nlines = ((float)(LASTSLICE-FIRSTSLICE)*ZOOM+.5);
  801.  
  802.   /* get floating distance values (first pass over data)*/
  803.     if (!views || (views & 1) || (views & 1<<2) || (views & 1<<4)) {
  804.       /* set scan start/end/directions depending on object/image
  805.          orientations */
  806.       scanxdir[0] = clipx[0];
  807.       scanxdir[1] = clipx[1];
  808.       scanxdir[2] = 1;
  809.       scanydir[0] = clipy[0];
  810.       scanydir[1] = clipy[1];
  811.       scanydir[2] = 1;
  812.       scanzdir[0] = FIRSTSLICE;
  813.       scanzdir[1] = LASTSLICE;
  814.       scanzdir[2] = 1;
  815.       imagex[0] = clipy[0];
  816.       imagex[1] = 1;
  817.       imagex[2] = 0;
  818.       imagex[3] = 1;
  819.       imagey[0] = 255 - clipx[0];
  820.       imagey[1] = -1;
  821.       imagey[2] = 0;
  822.       imagey[3] = 1;
  823.       imagez[0] = clipx[0];
  824.       imagez[1] = 1;
  825.       imagez[2] = clipy[0];
  826.       imagez[3] = 1;
  827.       if (object_or[0] == -1) {
  828.         scanxdir[0] = clipx[1]-1;
  829.         scanxdir[1] = clipx[0]-1;
  830.         scanxdir[2] = -1;
  831.         imagey[0] = clipx[1]-1;
  832.         imagez[0] = 256 - clipx[1];
  833.       }
  834.       if (object_or[1] == -1) {
  835.         scanydir[0] = clipy[1]-1;
  836.         scanydir[1] = clipy[0]-1;
  837.         scanydir[2] = -1;
  838.         imagex[0] = 256 - clipy[1];
  839.         imagez[2] = 256 - clipy[1];
  840.       }
  841.       if (object_or[2] == -1) {
  842.         scanzdir[0] = LASTSLICE;
  843.         scanzdir[1] = FIRSTSLICE;
  844.         scanzdir[2] = -1;
  845.       }
  846.       if (image_or[0] == -1) {
  847.         imagex[0] = 255-imagex[0];
  848.         imagex[1] *= -1;
  849.         imagey[0] = 255-imagey[0];
  850.         imagey[1] *= -1;
  851.         imagez[0] = 255-imagez[0];
  852.         imagez[1] *= -1;
  853.       }
  854.       if (image_or[1] == -1) {
  855.         imagex[2] = temp_nlines-imagex[2];
  856.         imagex[3] *= -1;
  857.         imagey[2] = temp_nlines-imagey[2];
  858.         imagey[3] *= -1;
  859.         imagez[2] = 255-imagez[2];
  860.         imagez[3] *= -1;
  861.       }
  862.  
  863.       printf("starting forward pass on data\n");
  864.       getdistances(scanxdir[0], scanxdir[1], scanxdir[2],
  865.                    scanydir[0], scanydir[1], scanydir[2],
  866.                    scanzdir[0], scanzdir[1], scanzdir[2],
  867.                    imagex[0], imagex[1], imagex[2], imagex[3],
  868.                    imagey[0], imagey[1], imagey[2], imagey[3],
  869.                    imagez[0], imagez[1], imagez[2], imagez[3],
  870.                    1);
  871.       /* create images */
  872.       if (views & 1) {
  873.     printf("doing bottom (Z) view(s)\n");
  874.         sprintf(filename,"%c:zdis%d.dat",DR,1);
  875.         doviews(filename,"gbo.out","dbo.out",256);
  876.       }
  877.       if (views & 1<<2) {
  878.     printf("doing right lateral (X) view(s)\n");
  879.         sprintf(filename,"%c:xdis%d.dat",DR,1);
  880.         doviews(filename,"grl.out","drl.out",NLINES);
  881.       }
  882.       if (views & 1<<4) {
  883.     printf("doing rear (Y) view(s)\n");
  884.         sprintf(filename,"%c:ydis%d.dat",DR,1);
  885.         doviews(filename,"gre.out","dre.out",NLINES);
  886.       }
  887.     }
  888.  /* get floating distance values (second pass over data)*/
  889.     if (!views || (views & 1<<1) || (views & 1<<3) || (views & 1<<5)) {
  890.       /* set scan start/end/directions depending on object/image
  891.          orientations */
  892.       scanxdir[0] = clipx[1]-1;
  893.       scanxdir[1] = clipx[0]-1;
  894.       scanxdir[2] = -1;
  895.       scanydir[0] = clipy[1]-1;
  896.       scanydir[1] = clipy[0]-1;
  897.       scanydir[2] = -1;
  898.       scanzdir[0] = LASTSLICE;
  899.       scanzdir[1] = FIRSTSLICE;
  900.       scanzdir[2] = -1;
  901.       imagex[0] = 256-clipy[1];
  902.       imagex[1] = 1;
  903.       imagex[2] = temp_nlines;
  904.       imagex[3] = -1;
  905.       imagey[0] = clipx[1]-1;
  906.       imagey[1] = -1;
  907.       imagey[2] = temp_nlines;
  908.       imagey[3] = -1;
  909.       imagez[0] = 256-clipx[1];
  910.       imagez[1] = 1;
  911.       imagez[2] = clipy[1]-1;
  912.       imagez[3] = -1;
  913.       if (object_or[0] == -1) {
  914.         scanxdir[0] = clipx[0];
  915.         scanxdir[1] = clipx[1];
  916.         scanxdir[2] = 1;
  917.         imagey[0] = 255 - clipx[0];
  918.         imagez[0] = clipx[0];
  919.       }
  920.       if (object_or[1] == -1) {
  921.         scanydir[0] = clipy[0];
  922.         scanydir[1] = clipy[1];
  923.         scanydir[2] = 1;
  924.         imagex[0] = clipy[0];
  925.         imagez[2] = 255 - clipy[0];
  926.       }
  927.       if (object_or[2] == -1) {
  928.         scanzdir[0] = FIRSTSLICE;
  929.         scanzdir[1] = LASTSLICE;
  930.         scanzdir[2] = 1;
  931.       }
  932.       if (image_or[0] == -1) {
  933.         imagex[0] = 255-imagex[0];
  934.         imagex[1] *= -1;
  935.         imagey[0] = 255-imagey[0];
  936.         imagey[1] *= -1;
  937.         imagez[0] = 255-imagez[0];
  938.         imagez[1] *= -1;
  939.       }
  940.       if (image_or[1] == -1) {
  941.         imagex[2] = temp_nlines-imagex[2];
  942.         imagex[3] *= -1;
  943.         imagey[2] = temp_nlines-imagey[2];
  944.         imagey[3] *= -1;
  945.         imagez[2] = 255-imagez[2];
  946.         imagez[3] *= -1;
  947.       }
  948.       printf("starting backward pass on data\n");
  949.       getdistances(scanxdir[0], scanxdir[1], scanxdir[2],
  950.                    scanydir[0], scanydir[1], scanydir[2],
  951.                    scanzdir[0], scanzdir[1], scanzdir[2],
  952.                    imagex[0], imagex[1], imagex[2], imagex[3],
  953.                    imagey[0], imagey[1], imagey[2], imagey[3],
  954.                    imagez[0], imagez[1], imagez[2], imagez[3],
  955.                    2);
  956.       /* create images */
  957.       if (views & 1<<1) {
  958.     printf("doing top (Z) view(s)\n");
  959.         sprintf(filename,"%c:zdis%d.dat",DR,2);
  960.         doviews(filename,"gto.out","dto.out",256);
  961.       }
  962.       if (views & 1<<3) {
  963.     printf("doing left lateral (X) view(s)\n");
  964.         sprintf(filename,"%c:xdis%d.dat",DR,2);
  965.         doviews(filename,"gll.out","dll.out",NLINES);
  966.       }
  967.       if (views & 1<<5) {
  968.     printf("doing front (Y) view(s)\n");
  969.         sprintf(filename,"%c:ydis%d.dat",DR,2);
  970.         doviews(filename,"gfr.out","dfr.out",NLINES);
  971.       }
  972.     }
  973.     printf("number of lines = %d\n",NLINES);
  974.  
  975.     /* write temporary drive and number of lines */ 
  976.     par=fopen("param.dat","w");
  977.     fprintf(par,"%c %d\n",DR,NLINES);
  978.     /* write views */
  979.     if (views & 1)
  980.       fprintf(par, "BO ");
  981.     if (views & (1 << 1))
  982.       fprintf(par, "TO ");
  983.     if (views & (1 << 2))
  984.       fprintf(par, "RL ");
  985.     if (views & (1 << 3))
  986.       fprintf(par, "LL ");
  987.     if (views & (1 << 4))
  988.       fprintf(par, "RE ");
  989.     if (views & (1 << 5))
  990.       fprintf(par, "FR ");
  991.     if (!views)
  992.       fprintf(par, "NO ");
  993.     fprintf(par, "\n");
  994.     /* write clipping planes */
  995.     fprintf(par, "RL %d LL %d RE %d FR %d\n", clipx[0], clipx[1],
  996.                         clipy[0], clipy[1]);
  997.     /* write image and object orientation */
  998.     fprintf(par, "image: x %2d  y %2d\n", image_or[0], image_or[1]);
  999.     fprintf(par, "object: x %2d  y %2d  z %2d\n",
  1000.         object_or[0], object_or[1], object_or[2]);
  1001.     fclose(par);
  1002. }
  1003.  
  1004.  
  1005. usage()
  1006. {
  1007.   printf("%s\n", usestr1);
  1008.   printf("%s\n", usestr2);
  1009.   printf("%s\n", usestr3);
  1010.   exit(1);
  1011. }
  1012.