home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource4 / 293_01 / deb.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-08-25  |  7.8 KB  |  252 lines
  1. /************************* deb.c ************************************
  2.  
  3.       3-D Reconstruction of Medical Images
  4.  
  5.     Three Dimensional Reconstruction Of Medical
  6.     Images from Serial Slices - CT, MRI, Ultrasound
  7.  
  8.  
  9.    These programs process a set of slices images (scans) for one
  10.    patient. It outputs two sets of files containing nine predefined
  11.    views of bony surfaces. One set contains distance values and
  12.    the other gradient values.
  13.  
  14.    The distance values are used as 3-D spatial topographic surface
  15.    coordinate maps for geometrical analysis of the scanned object.
  16.  
  17.    The gradient values are used for rendering the surface maps on
  18.    CRT displays for subjective viewing where perception of small
  19.    surface details is important.
  20.  
  21.     Daniel Geist, B.S.
  22.     Michael W. Vannier, M.D.
  23.  
  24.     Mallinckrodt Institute of Radiology
  25.     Washington University School of Medicine
  26.     510 S. Kingshighway Blvd.
  27.     St. Louis, Mo. 63110
  28.  
  29.     These programs may be copied and used freely for non-commercial
  30.     purposes by developers with inclusion of this notice.
  31.  
  32.  
  33. ********************************************************************/
  34. #include <stdio.h>
  35. #include <math.h>
  36. int outfile,ZMAX,FIRSTSLICE,LASTSLICE,THRESHOLD,
  37.     RIGHTMID,LEFTMID,TOPINT,midslice,midline,NLINES;
  38. float ZOOM;
  39. int  buffer[2][6][256];
  40. float fxbuf[5][256],fybuf[5][256];
  41. float huge fzbuf[256][256];
  42. /*            standard 18 output files ( 9 views x 2) */
  43. char *fnamein="ctbild.000",*fgll="gll.out",*fdll="dll.out";
  44. succ(i)
  45. int i;
  46. {return(i==2?0:i+1);}
  47. prev(i)
  48. int i;
  49. {return(i==0?2:i-1);}
  50. setfilename(filenum)
  51. int filenum;
  52. {fnamein[7]=filenum/100+'0';
  53.  fnamein[8]=(filenum%100)/10+'0';
  54.  fnamein[9]='0'+filenum%10;
  55. }
  56. interpolate(line,bot,top,n)
  57. int line,bot,top,n;
  58. {int next,i,j,x,inc;
  59.  inc=top>bot?1:(-1);
  60.  next=bot+inc;
  61.  for(i=1,j=n-1;i<n;i++,j--){
  62.       for(x=0;x<256;x++) buffer[line][next][x]=
  63.         (buffer[line][bot][x]*j+buffer[line][top][x]*i)/n;
  64.       next+=inc;
  65.  }
  66. }
  67. float midpoint(b,a)
  68. float b,a;
  69. {return( (THRESHOLD-a) / (b-a) );}
  70. getdistances()
  71. {int z,x,y,i,j,start,stop,inc,line,rzoom,inter;
  72.  float remain;
  73.  FILE *fxfloat,*fyfloat,*fzfloat,*fn[2];
  74.   NLINES=0;
  75.   remain=0;
  76.   rzoom=ZOOM+0.5;
  77.   fxfloat=fopen("xdis.dat","wb");
  78.   fyfloat=fopen("ydis.dat","wb");
  79.   for(i=0;i<256;i++)
  80.    for(j=0;j<256;j++)fzbuf[i][j]=256;
  81.   for(z=0;z<(LASTSLICE-FIRSTSLICE);z++){
  82.       for(i=0;i<2;i++){
  83.           setfilename(LASTSLICE-(z+i));
  84.           fn[i]=fopen(fnamein,"rb");
  85.           fseek(fn[i],(long)512,SEEK_SET);
  86.       }
  87.       inter=rzoom;
  88.       remain+=rzoom-ZOOM;
  89.       if(remain>=1){
  90.           inter-=1;
  91.           remain-=1;
  92.       }
  93.       else if(remain<=(-1)){
  94.           inter+=1;
  95.           remain+=1;
  96.       }
  97.       line=0;
  98.       for(j=0;j<inter;j++)
  99.        for(i=0;i<256;i++)fxbuf[j][i]=fybuf[j][i]=256;
  100.       for(y=255;y>=0;y--){
  101.           fseek(fn[0],(long)512*(y+1),SEEK_SET);
  102.           fread(buffer[line][0],1,512,fn[0]);
  103.           fseek(fn[1],(long)512*(y+1),SEEK_SET);
  104.           fread(buffer[line][inter],1,512,fn[1]);
  105.           interpolate(line,0,inter,inter);
  106.           for(i=0;i<inter;i++){
  107.              for(x=255;x>=0;x--)if (buffer[line][i+1][x]>=THRESHOLD){
  108.                 if(fxbuf[i][y]==256.0) fxbuf[i][y]=(x==255)?0:255-(x+1)+
  109.                       midpoint((float)buffer[line][i+1][x],
  110.                               (float)buffer[line][i+1][x+1]);
  111.                 if(fzbuf[y][x]==256.0) fzbuf[y][x]=
  112.                      (z==0) && (buffer[line][0][x]>=THRESHOLD)?0:NLINES+i+
  113.                       midpoint((float)buffer[line][i+1][x],
  114.                                        (float)buffer[line][i][x]);
  115.              }
  116.           }
  117.           line=1-line;
  118.       }
  119.       NLINES+=inter;  
  120.       fclose(fn[0]);
  121.       fclose(fn[1]);
  122.       fwrite(fxbuf,1,inter*1024,fxfloat);
  123.       fwrite(fybuf,1,inter*1024,fyfloat);
  124.       printf("did %d \n",z);
  125.   }
  126.   fclose(fxfloat);
  127.   fclose(fyfloat);
  128.   fzfloat=fopen("zdis.dat","wb");
  129.   for(i=0;i<256;i++){fwrite(fzbuf[i],1,1024,fzfloat);printf("wr %d",i);}
  130.   fclose(fzfloat);
  131. }
  132. unsigned char gradx(y1,y2,z1,z2)
  133. float y1,y2,z1,z2;
  134. {float gx,gy,gz;
  135.  unsigned char gxint;
  136.      /* get z and y components of gradient */
  137.   gz=(z1-z2)/2;
  138.   gy=(y1-y2)/2;
  139.      /*compute gx - normalized x component of gradient */
  140.   gx=1/sqrt(1+gz*gz+gy*gy);
  141.   gxint=255*gx+0.5;      /*scale gx by 256 */
  142.   return(gxint);
  143. }
  144. doviews(namedis,nameg,named,nlines)
  145. char *namedis,*nameg,*named;
  146. int nlines;
  147. {FILE *fg,*fd,*ffloat;
  148.  int z,i,j,k,midline;
  149.  char lined[256],lineg[256];
  150.  midline=1;
  151.  fd=fopen(named,"wb");
  152.  fg=fopen(nameg,"wb");
  153.  ffloat=fopen(namedis,"rb");
  154.  fread(fxbuf,1,3072,ffloat);
  155.  /* do first line */
  156.  if(fxbuf[0][0]==256.0)lineg[0]=lined[0]=0;
  157.  else{
  158.      lined[0]=256-fxbuf[0][0];
  159.      lineg[0]=gradx(fxbuf[0][0]*2,fxbuf[0][1]*2,fxbuf[0][0]*2,fxbuf[1][0]*2);
  160.  }
  161.  for(i=1;i<255;i++)
  162.   if(fxbuf[0][i]==256.0) lineg[i]=lined[i]=0;
  163.   else{
  164.       lined[i]=256-fxbuf[0][i];
  165.       lineg[i]=gradx(fxbuf[0][i-1],fxbuf[0][i+1],2*fxbuf[0][i],2*fxbuf[1][i]);
  166.  }
  167.  if(fxbuf[0][255]==256.0)lineg[255]=lined[255]=0;
  168.  else{
  169.      lined[255]=256-fxbuf[0][255];
  170.      lineg[255]=
  171.        gradx(fxbuf[0][255]*2,fxbuf[0][254]*2,fxbuf[0][255]*2,fxbuf[1][255]*2);
  172.  }
  173.  fwrite(lineg,1,256,fg);
  174.  fwrite(lined,1,256,fd);
  175.  printf("Begining computation of  views\n");
  176.  for(z=0;z<(nlines-2);z++){      /*for each slice */
  177.      if(fxbuf[midline][0]==256.0)lineg[0]=lined[0]=0;
  178.      else{
  179.          lined[0]=256-fxbuf[midline][0];
  180.          lineg[0]=
  181.            gradx(fxbuf[midline][0]*2,fxbuf[midline][1]*2,
  182.                    fxbuf[prev(midline)][0],fxbuf[succ(midline)][0]);
  183.      }
  184.      for(i=1;i<255;i++)
  185.        if(fxbuf[midline][i]==256.0) lineg[i]=lined[i]=0;
  186.       else{
  187.           lined[i]=256-fxbuf[midline][i];
  188.           lineg[i]=gradx(fxbuf[midline][i-1],fxbuf[midline][i+1],
  189.            fxbuf[prev(midline)][i],fxbuf[succ(midline)][i]);
  190.      }
  191.      if(fxbuf[midline][255]==256.0)lineg[255]=lined[255]=0;
  192.      else{
  193.          lined[255]=256-fxbuf[midline][255];
  194.         lineg[255]=gradx(fxbuf[midline][255]*2,fxbuf[midline][254]*2
  195.             ,fxbuf[prev(midline)][255],fxbuf[succ(midline)][255]);
  196.      }
  197.      fwrite(lineg,1,256,fg);
  198.      fwrite(lined,1,256,fd);
  199.      fread(fxbuf[prev(midline)],1,1024,ffloat);
  200.      midline=succ(midline);
  201.      printf(" did %d \n",z);
  202.  }
  203.  /* do last line */
  204.  if(fxbuf[midline][0]==256.0)lineg[0]=lined[0]=0;
  205.  else{
  206.      lined[0]=256-fxbuf[midline][0];
  207.      lineg[0]=gradx(fxbuf[midline][0]*2,fxbuf[midline][1]*2,
  208.          fxbuf[midline][0]*2,fxbuf[prev(midline)][0]*2);
  209.  }
  210.  for(i=1;i<255;i++)
  211.   if(fxbuf[midline][i]==256.0) lineg[i]=lined[i]=0;
  212.   else{
  213.       lined[i]=256-fxbuf[0][i];
  214.       lineg[i]=gradx(fxbuf[midline][i-1],fxbuf[midline][i+1],
  215.          2*fxbuf[midline][i],2*fxbuf[prev(midline)][i]);
  216.  }
  217.  if(fxbuf[midline][255]==256.0)lineg[255]=lined[255]=0;
  218.  else{
  219.      lined[255]=256-fxbuf[midline][255];
  220.      lineg[255]=
  221.        gradx(fxbuf[midline][255]*2,fxbuf[midline][254]*2,
  222.          fxbuf[midline][255]*2,fxbuf[prev(midline)][255]*2);
  223.  }
  224.  fwrite(lineg,1,256,fg);
  225.  fwrite(lined,1,256,fd);
  226.  fclose(fg);
  227.  fclose(fd);
  228.  fclose(ffloat);
  229. }
  230. /**********************************************************/
  231. /**** MAIN ***** MAIN ***** MAIN ***** MAIN ***** MAIN ****/
  232. /**********************************************************/
  233. main()
  234. {
  235.  /* first get some parameters from user */
  236.  printf("Enter Zoom factor: ");
  237.  scanf("%f",&ZOOM);
  238.  printf("Enter Starting scan number: ");
  239.  scanf("%d",&FIRSTSLICE);
  240.  printf("Enter ending scan number: ");
  241.  scanf("%d",&LASTSLICE);
  242.  printf("Enter threshold number: ");
  243.  scanf("%d",&THRESHOLD);
  244.  THRESHOLD+=1024;
  245.  getdistances();
  246.  printf("doing left views");
  247.  doviews("xdis.dat","gll.out","dll.out",NLINES);
  248.  printf("doing bottom views");
  249.  doviews("zdis.dat","gto.out","dto.out",256);
  250.  printf("numberof lines %d\n",NLINES);
  251. }
  252.