/******************** GRAD ******************************
  grad.c
  Daniel Geist - Mallinckrodt Institue of Radiology 1987.
   This program reads a set of ct scans for one
   patient. It outputs two sets of files containing predefined 
   views of bone surfaces. One set contains distance values and
   the other gradient values.

   Danny Geist & Mike Vannier

   Contact for further information:
	Michael W. Vannier
	Mallinckrodt Institute of Radiology
	Washington University School of Medicine
	510 South Kingshighway Blvd.
	St. Louis, Mo. 63110 USA

   Ake Wallin has made the following changes - 1988

   -Command line input
   -Threshold of gradient for differencing
   -Views to choose: BO (bottom), TO (top), RL (right lateral),
    LL (left lateral), RE (rear), FR (front), NO (none, just leave distance)
   -Clipping planes: RL (right lateral), LL (lateral left), RE (rear),
    FR (front)
   -Scan of data in "correct" order, assuming:
    object:
    low z - feet, high z - head
    low x - right, high x - left
    low y - back, high y - front
    image:
    low x - left, high x - right
    low y - down, high y - up
		    front
    view BO : right      left
		    back
		    front
    view TO : left       right
		    back
		    head
    view RL : back      front
		    feet
		    head
    view LL : front     back
		    feet
		    head
    view RE : left      right
		    feet
		    head
    view FR : right     left
		    feet
*********************************************************/
#include <stdio.h>
#include <math.h>
#include <ctype.h>

/* some global variables*/
int FIRSTSLICE,LASTSLICE,THRESHOLD,NLINES;
float ZOOM;
char DR;
int  buffer[2][6][256];          /* input buffer */
float fxbuf[5][256],fybuf[5][256]; /* X,Y Views floating buffers */

float huge fzbuf[256][256];        /* Z view floating buffer */
char fnamein[50]="ctbild.000";        /* input file name  */
float GRAD_THRESHOLD = 10.0;
unsigned int   views = 0xFFFF;
int   clipx[2], clipy[2];      /* clipping planes */
int   dispmode;                /* distance and/or gradient */
int   image_or[2];            /* image oriention in x and y */
int   object_or[3];           /* image oriention in x, y and z */
int   header_blocks;          /* number of header blocks in CT's */
char *usestr1 = "Usage: grad [file] [-z] [-f] [-l] [-t] [-d] [-n(g|d)] [-g]";
char *usestr2 = "            [-h] [-v(bo|to|rl|ll|re|fr|no)] [-c(rl|ll|re|fr)]";
char *usestr3 = "            [-i(x(r|l)|y(u|d))] [-o(x(r|l)|y(f|b)|z(h|f))]";

succ(i)   /* (i+1) modolus 3 */
int i;
{
    return(i==2?0:i+1);
}

prev(i)   /* (i-1) modolus 3 */
int i;
{
    return(i==0?2:i-1);
}

/* Set input file - add slice number to extension */
setfilename(filenum)
int filenum;
{
    int i;

    for (i = strlen(fnamein)-1; i > 0; i--)
      if (fnamein[i] == '.') break;
    if (i == 0) {
      printf("CT file names inconsistend!\n");
      exit(1);
    }
    fnamein[++i]=filenum/100+'0';
    fnamein[++i]=(filenum%100)/10+'0';
    fnamein[++i]='0'+filenum%10;
}

/* interpolate from bottom line to top line n-1 lines  */
interpolate(line,bot,top,n)
int line,bot,top,n;
{
    int next,i,j,x,inc;
    inc=top>bot?1:(-1);    /* interpolate backward or forwards ? */
    next=bot+inc;
    for(i=1,j=n-1;i<n;i++,j--){   /* do for each next line of interpolation */
        for(x=0;x<256;x++) buffer[line][next][x]=
            (buffer[line][bot][x]*j+buffer[line][top][x]*i)/n;
        next+=inc;
    }
}

/* midpoint - fraction part of threshold transition distance */ 
float midpoint(b,a)
float b,a;
{
    return( (THRESHOLD-a) / (b-a) );
}

/* get floating point distance values  */
getdistances(xstart,xstop,xdir,
             ystart,ystop,ydir,
             start_slice,end_slice,zdir,
             xbufxstart, xbufxdir, xbufystart, xbufydir,
             ybufxstart, ybufxdir, ybufystart, ybufydir,
             zbufxstart, zbufxdir, zbufystart, zbufydir,
             pass)
int xstart,xstop,xdir,ystart,ystop,ydir,start_slice,end_slice,zdir,pass;
int  xbufxstart, xbufxdir, xbufystart, xbufydir;
int  ybufxstart, ybufxdir, ybufystart, ybufydir;
int  zbufxstart, zbufxdir, zbufystart, zbufydir;
/* also used are global variables: fxbuf, fybuf, fzbuf, buffer,
   ZOOM, DR, image_or, object_or. */
{
    int z,x,y,i,j,start,stop,inc,line,rzoom,inter;
    float remain;
    FILE *fxfloat,*fyfloat,*fzfloat,*fn[2];
    char filename[13];
    int  filegap;
    int  xbfx, ybfx, zbfx, zbfy;

    NLINES=0;         /* number of output lines in X,Y view directions */
    remain=0;         /* remainder of interpolation after roundoff */
    rzoom=ZOOM+0.5;   /* rounded zoom factor */
    /* X view floating output file */
    sprintf(filename,"%c:xdis%d.dat",DR,pass);
    if ((fxfloat=fopen(filename,"wb")) == NULL) {
      printf("Error creating %s!\n", filename);
      exit(1);
    }
    /* Y view floating output file */
    sprintf(filename,"%c:ydis%d.dat",DR,pass);
    if ((fyfloat=fopen(filename,"wb")) == NULL) {
      printf("Error creating %s!\n", filename);
      exit(1);
    }
    for(i=0;i<256;i++)
        for(j=0;j<256;j++)fzbuf[i][j]=256; /* clear Z view buffer */
    
    for(z=start_slice; z!=end_slice; z += zdir){  /* For each Slice */
      /* open next two slice files */
	  setfilename(z);
      if ((fn[0]=fopen(fnamein,"rb")) == NULL) continue;
	for (filegap = zdir; filegap + z != end_slice+zdir; filegap+=zdir) {
          setfilename(z+filegap);
          if ((fn[1]=fopen(fnamein,"rb")) != NULL) break;
	  }
      if (fn[1] == NULL) continue;

      inter=rzoom; /* interpolation factor assumed rounded zoom */
      /* correct interpolation factor according to floating remainder */
      remain+=rzoom-ZOOM;
      if(remain>=1){
        inter-=1;
        remain-=1;
      }
      else if(remain<=(-1)){
        inter+=1;
        remain+=1;
      }

      line=0;               /* current input buffer line */
      for(j=0;j<inter;j++)  /* clear X,Y floating buffers */
	    for(i=0;i<256;i++) fxbuf[j][i]=fybuf[j][i]=256;
      /* set index for zbuf */
      zbfy = zbufystart;
      /* set index for xbuf */
      xbfx = xbufxstart;
      for(y=ystart; y!=ystop; y+=ydir){              /* For each line */
	/* skip to line*/
	for(i=0;i<2;i++)fseek(fn[i],(long)512*(y+header_blocks),SEEK_SET);
        fread(buffer[line][0],1,512,fn[0]);   
        fread(buffer[line][inter],1,512,fn[1]);
        interpolate(line,0,inter,inter); /* interpolate in_between */
	for(i=0;i<inter;i++){   /* For each interpolation line */
          /* set index for zbuf */
	  zbfx = zbufxstart;
          /* set index for ybuf */
          ybfx = ybufxstart;
	  for(x=xstart; x!=xstop; x+=xdir) { /* For each Voxel value */
            /* find threshold transition*/
	    if (buffer[line][i+1][x] >= THRESHOLD) {
              /* if first transition in X direction get floating
	             distance */
              if(fxbuf[i][xbfx]==256.0) fxbuf[i][xbfx]=(x==xstart)?0:
                                    (x-xstart)*xdir-1+
                     midpoint((float)buffer[line][i+1][x],
                            (float)buffer[line][i+1][x-xdir]);

		      /* if first transition in Y direction get floating
                   distance */
              if(fybuf[i][ybfx]==256.0) fybuf[i][ybfx]=(y==ystart)?0:
                    (y-ystart)*ydir-1+
                     midpoint((float)buffer[line][i+1][x],
                              (float)buffer[1-line][i+1][x]);
 
	          /* if first transition in Z direction get floating
                   distance */
              if(fzbuf[zbfy][zbfx]==256.0) fzbuf[zbfy][zbfx]=
                   (i==0) && (buffer[line][i][x]>=THRESHOLD) ?
                    NLINES : NLINES+i+
                    midpoint((float)buffer[line][i+1][x],
                             (float)buffer[line][i][x]);
            }
            /* change index of ybuf and zbuf */
            ybfx += ybufxdir;
	    zbfx += zbufxdir;
	    if ((ybfx < 0 || ybfx > 255) && x+xdir != xstop) {
              printf("Error in index YBFX!\n");
              exit(1);
            }
	    if ((zbfx < 0 || zbfx > 255) && x+xdir != xstop) {
              printf("Error in index ZBFX!\n");
              exit(1);
            }
          }
        }
        line=1-line;  /* change current input buffer line */
        /* update index of xbuf and zbuf */
        xbfx += xbufxdir;
        zbfy += zbufydir;
	if ((zbfy < 0 || zbfy > 255) && y+ydir != ystop) {
          printf("Error in index ZBFY!\n");
          exit(1);
        }
	if ((xbfx < 0 || xbfx > 255) && y+ydir != ystop) {
          printf("Error in index XBFX!\n");
          exit(1);
        }
      }
      NLINES+=inter;           /* increment number of output lines */

      fclose(fn[0]);           /* close slice files */
      fclose(fn[1]);
      /* write output lines to files */
      if (fwrite(fxbuf,1,sizeof(float)*inter*256,fxfloat) <
		inter*256*sizeof(float))
	printf("Error writing to FXBUF!\n");
      if (fwrite(fybuf,1,sizeof(float)*inter*256,fyfloat) <
		inter*256*sizeof(float))
	printf("Error writing to FYBUF!\n");
      printf(" line %d\r", z);
    }
    printf("\n");
    fclose(fxfloat); /* close X,Y floating files */
    fclose(fyfloat);
    /* write Z floating file */
    sprintf(filename,"%c:zdis%d.dat",DR,pass);
    fzfloat=fopen(filename,"wb");
    for(i=0;i<256;i++)fwrite(fzbuf[i],sizeof(float),256,fzfloat);
    fclose(fzfloat);

    /* should we reverse the order of the lines? */
    if (xbufydir < 0) {
      sprintf(filename,"%c:xdis%d.dat",DR,pass);
      fxfloat=fopen(filename,"rb");
      for (line = 0; line < NLINES; line++)
        fread(fzbuf[line], 1, sizeof(float)*256, fxfloat);
      fclose (fxfloat);
      fxfloat=fopen(filename,"wb");
      for (line = 1; line <= NLINES; line++)
        fwrite(fzbuf[NLINES-line], sizeof(float), 256, fxfloat);
      fclose (fxfloat);
    }
    if (ybufydir < 0) {
      sprintf(filename,"%c:ydis%d.dat",DR,pass);
      fyfloat=fopen(filename,"rb");
      for (line = 0; line < NLINES; line++)
        fread(fzbuf[line], 1, sizeof(float)*256, fyfloat);
      fclose (fyfloat);
      fyfloat=fopen(filename,"wb");
      for (line = 1; line <= NLINES; line++)
        fwrite(fzbuf[NLINES-line], sizeof(float), 256, fyfloat);
      fclose (fyfloat);
    }

}

unsigned char grad(y1,y2,z1,z2,y_factor,z_factor)
float y1,y2,z1,z2;
int y_factor,z_factor;
{
    float gx,gy,gz;
    unsigned char gxint;
    /* get z and y components of gradient */
    gz=(z1-z2)/(float) z_factor;
    gy=(y1-y2)/(float) y_factor;
    /*compute gx - normalized x component of gradient */
    gx=1/sqrt(1+gz*gz+gy*gy);
    gxint=255*gx+0.5;      /*scale gx by 256 and round */
    return(gxint);
}

/* get gradient and depth shades for one image line */
doline(lineg,lined,line,prevline,succline,z_factor,fg,fd)
unsigned char lineg[],lined[];  /*output buffers */
int line,prevline,succline;  /* input buffer configuration */
int z_factor; /* for choosing forward, backward or central differences */
FILE *fg,*fd;  /* output files */
{
    int i;
    float y1, y2, z1, z2;
    int factor;
    /* do first pixel in line */
    i = 0;
    if(fxbuf[line][0]==256)lineg[0]=lined[0]=0;
    else{
        if (dispmode & 1)
          lined[0]=255-fxbuf[line][0]; /*distance shade */
        if (dispmode & 2)
          lineg[0]=grad(fxbuf[line][0],fxbuf[line][1],fxbuf[prevline][0],