home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource4 / 293_01 / angpp.c < prev    next >
C/C++ Source or Header  |  1989-08-25  |  9KB  |  321 lines
  1. /********************  angpp.c  ******************************
  2.  
  3.    Three dimensional surface reconstruction program
  4.  
  5.    This program assumes that the views from the main axis
  6.    directions were done thus an angle view between two main
  7.    axis views can be created from the depth_code viewa(floating
  8.    point). 
  9.  
  10.  Daniel Geist
  11.  Michael W. Vannier
  12.  
  13.  Mallinckrodt Institute of Radiology
  14.  Washington University School of Medicine
  15.  510 S. Kingshighway Blvd.
  16.  St. Louis, Mo. 63110
  17.  
  18.  Note:  This program is intended for the private non-commercial
  19.  use of interested individuals to provide a fuller explanation of
  20.  the algorithms described in "Three dimensional reconstruction in
  21.  Medical Imaging" by D. Geist and M. Vannier.
  22.  
  23.  1988
  24.  *************************************************************/
  25.  
  26. #include <stdio.h>
  27. #include <math.h>
  28. #define FLOAT_LINE 1024
  29. #define PI 3.141592653
  30.  
  31. /* the structure below is the data for a point on the surface projected
  32.    on the view plane                        */
  33. typedef struct DIS_REC {
  34.     float dist;            /*distance from view plane */
  35.     int indXY;             /* index of data on main axis view */
  36.     char XY;               /* Which main axis view obtained from (X or Y)*/
  37. };
  38.  
  39. struct DIS_REC distances[256]; /* one projected line */
  40.  
  41. int NLINES,IXmax,IYmax,xinc,yinc;
  42. int (*xdir)(),(*ydir)();
  43. float THETA,cosTheta,sinTheta,tgnTheta;
  44. float fxbuf[3][256],fybuf[3][256];/* input buffers */
  45.  
  46. /*  output files  */
  47. char fng[13],fnd[13],xfile[13],yfile[13],DR;
  48.  
  49. succ(i)
  50. int i;
  51. {
  52.     return(i==2?0:i+1);
  53. }
  54.  
  55. prev(i)
  56. int i;
  57. {
  58.     return(i==0?2:i-1);
  59. }
  60.  
  61. forward(i)
  62. int i;
  63. { 
  64.     return(i); 
  65. }
  66.  
  67. backward(i,start)
  68. int i;
  69. {
  70.     return(start-i); 
  71. }
  72.  
  73. /* PUTD - fill a DIS_REC with values */
  74. putd(pdis,D,i,xy_sym)
  75. struct DIS_REC *pdis;
  76. float D;
  77. int i,xy_sym;
  78.  
  79. { 
  80.     pdis->dist=D;
  81.     pdis->indXY=i;
  82.     pdis->XY=xy_sym;
  83. }
  84.  
  85. /* take one line from X and Y views and create a projrcted line */
  86. getdistances(linex,liney)
  87. float linex[],liney[];
  88.  
  89. { 
  90.     int i,IND,X,Y;
  91.     float D;
  92.     for(i=0;i<256;i++) distances[i].XY=0;
  93.     /* project Y-data onto image line */
  94.     for(i=0;i<256;i++) if(liney[i]!=256.0){
  95.         X=(*ydir)(i,255);
  96.         IND=IXmax -(X*sinTheta-liney[i]*cosTheta)+0.5;
  97.         D=liney[i]/sinTheta+(X-liney[i]/tgnTheta)*cosTheta;
  98.         if((IND>=0) && (IND<256)){
  99.             if(distances[IND].XY==0)putd(&distances[IND],D,i,1);
  100.             else if(distances[IND].dist>D)putd(&distances[IND],D,i,1);
  101.         }
  102.     }
  103.  
  104.     /*project X-data onto image plane */
  105.     for(i=0;i<256;i++) if(linex[i]!=256.0){
  106.         Y=(*xdir)(i,255);
  107.         IND=IXmax-(linex[i]*sinTheta-Y*cosTheta)+0.5;
  108.         D=Y/sinTheta+(linex[i]-Y/tgnTheta)*cosTheta;
  109.         if((IND>=0) && (IND<256)){
  110.             if(distances[IND].XY==0)putd(&distances[IND],D,i,2);
  111.             else if(distances[IND].dist>D)putd(&distances[IND],D,i,2);
  112.         }
  113.     }
  114.  
  115.     /* fill holes due to low resolution */
  116.     for(i=1;i<255;i++)if( (distances[i].XY==0) && (distances[i+1].XY!=0) &&
  117.         (distances[i-1].XY!=0))
  118.         putd(&distances[i],distances[i-1].dist,
  119.         distances[i-1].indXY,(int)distances[i-1].XY);
  120. }
  121.  
  122. /* returns gradient shade in point given the variations of the surface */
  123. unsigned char grad(x1,x2,y1,y2,z1,z2,x_fac,y_fac,z_fac)
  124. float x1,x2,y1,y2,z1,z2;
  125. int x_fac,y_fac,z_fac;
  126. {
  127.     float gx,gy,gz,G,nx,ny;
  128.     unsigned char gxint;
  129.     /* components of gradient */
  130.     gx=(x2-x1)/x_fac;
  131.     gy=(y2-y1)/y_fac;
  132.     gz=(z1-z2)/z_fac;
  133.     G=sqrt(gx*gx+gy*gy+gz*gz);
  134.     /*compute nx,ny normalized x,y component of gradient */
  135.     nx=gx/G;
  136.     ny=gy/G;
  137.     gxint=255*(nx*cosTheta+ny*sinTheta)+0.5; /*scale gradient shade by 256 */
  138.     return(gxint);
  139. }
  140.  
  141. doline(linex,linex1,linex2,liney,liney1,liney2,z_fac,fg,fd)
  142. float linex[],linex1[],linex2[],liney[],liney1[],liney2[];
  143. int z_fac;
  144. FILE *fg,*fd;
  145.  
  146. { 
  147.     int i;
  148.     unsigned char lined[256],lineg[256];
  149.     /* empty bit on image line */
  150.     for(i=0;i<256;i++)if(distances[i].XY==0)lineg[i]=lined[i]=0;
  151.     else{
  152.         lined[i]=(distances[i].dist<256)?255-distances[i].dist+0.5:0;
  153.         /* bit on image line projected from Y view */
  154.         if(distances[i].XY==1) switch(distances[i].indXY){
  155.         case 0:
  156.             lineg[i]=
  157.                 grad(liney[1]*yinc,liney[0]*yinc,(float)0,(float)2,
  158.             liney1[0],liney2[0],1,2,z_fac);
  159.             break;
  160.         case 255:
  161.             lineg[i]=
  162.                 grad(liney[255]*yinc,liney[254]*yinc,(float)0,(float)2,
  163.             liney1[255],liney2[255],1,2,z_fac);
  164.             break;
  165.         default:
  166.             lineg[i]=
  167.                 grad(liney[distances[i].indXY+yinc],
  168.             liney[distances[i].indXY-yinc],(float)0,(float)2,
  169.             liney1[distances[i].indXY],
  170.             liney2[distances[i].indXY],2,2,z_fac); 
  171.             break;
  172.         }
  173.         /* bit on image line projected from X view */
  174.         else switch(distances[i].indXY){
  175.         case 0:
  176.             lineg[i]=
  177.                 grad((float)0,(float)2,linex[1]*xinc,linex[0]*xinc,
  178.             linex1[0],linex2[0],2,1,z_fac);
  179.             break;
  180.         case 255:
  181.             lineg[i]=
  182.                 grad((float)0,(float)2,linex[255]*xinc,linex[254]*xinc,
  183.             linex1[255],linex2[255],2,1,z_fac);
  184.             break;
  185.         default:
  186.             lineg[i]=
  187.                 grad((float)0,(float)2,linex[distances[i].indXY+xinc],
  188.             linex[distances[i].indXY-xinc],
  189.             linex1[distances[i].indXY],
  190.             linex2[distances[i].indXY],2,2,z_fac);
  191.             break;
  192.         }
  193.     }
  194.  
  195.     fwrite(lineg,1,256,fg);
  196.     fwrite(lined,1,256,fd);
  197. }
  198.  
  199. doviews()
  200. {
  201.     FILE *fg,*fd,*fx,*fy;
  202.     int z,i,j,k,mid;
  203.     mid=1;
  204.     fd=fopen(fnd,"wb");
  205.     fg=fopen(fng,"wb");
  206.     fx=fopen(xfile,"rb");
  207.     fy=fopen(yfile,"rb");
  208.  
  209.     /*read first three lines */
  210.     for(i=0;i<3;i++){
  211.         fseek(fx,(long)FLOAT_LINE*(*ydir)(i,NLINES-1),SEEK_SET);
  212.         fread(fxbuf[i],1,FLOAT_LINE,fx);
  213.         fseek(fy,(long)FLOAT_LINE*(*xdir)(i,NLINES-1),SEEK_SET);
  214.         fread(fybuf[i],1,FLOAT_LINE,fy);
  215.     }
  216.  
  217.     /* do first line (forward differene)*/
  218.     getdistances(fxbuf[0],fybuf[0]);
  219.     doline(fxbuf[0],fxbuf[0],fxbuf[1],fybuf[0],fybuf[0],fybuf[1],1,fg,fd);
  220.  
  221.     /* do middle lines (central diffrence) */
  222.     printf("Begining computation of  views\n");
  223.     for(z=0;z<(NLINES-2);z++){      /*for each slice */
  224.         getdistances(fxbuf[mid],fybuf[mid]);
  225.         doline(fxbuf[mid],fxbuf[prev(mid)],fxbuf[succ(mid)],
  226.         fybuf[mid],fybuf[prev(mid)],fybuf[succ(mid)],2,fg,fd);
  227.         fseek(fx,(long)FLOAT_LINE*(*ydir)(i,NLINES-1),SEEK_SET);
  228.         fseek(fy,(long)FLOAT_LINE*(*xdir)(i,NLINES-1),SEEK_SET);
  229.         i++;
  230.         fread(fxbuf[prev(mid)],1,FLOAT_LINE,fx);
  231.         fread(fybuf[prev(mid)],1,FLOAT_LINE,fy);
  232.         mid=succ(mid);
  233.         printf(" did %d \n",z);
  234.     }
  235.  
  236.     /* do last line (backward difference)*/
  237.     getdistances(fxbuf[mid],fybuf[mid]);
  238.     doline(fxbuf[mid],fxbuf[prev(mid)],fxbuf[mid],
  239.     fybuf[mid],fybuf[prev(mid)],fybuf[mid],1,fg,fd);
  240.  
  241.     fclose(fg);
  242.     fclose(fd);
  243.     fclose(fx);
  244.     fclose(fy);
  245. }
  246.  
  247. getang(ang)
  248. int *ang;
  249. { 
  250.     scanf("%d",ang);
  251.     sprintf(fng,"gang%d.out",*ang);
  252.     sprintf(fnd,"dang%d.out",*ang);
  253.     if( ( *ang > 0) && ( *ang < 90) ){
  254.         xdir=forward;
  255.         ydir=forward;
  256.         sprintf(xfile,"%c:xdis1.dat",DR);
  257.         sprintf(yfile,"%c:ydis1.dat",DR);
  258.         xinc=1;
  259.         yinc=1;
  260.     }
  261.     else if( ( *ang > 90) && ( *ang < 180) ){
  262.         xdir=backward;
  263.         ydir=forward;
  264.         sprintf(xfile,"%c:ydis1.dat",DR);
  265.         sprintf(yfile,"%c:xdis2.dat",DR);
  266.         xinc=-1;
  267.         yinc=1;
  268.         *ang-=90;
  269.     }
  270.     else if( ( *ang > 180) && ( *ang < 270) ){
  271.         xdir=backward;
  272.         ydir=backward;
  273.         sprintf(xfile,"%c:xdis2.dat",DR);
  274.         sprintf(yfile,"%c:ydis2.dat",DR);
  275.         xinc=-1;
  276.         yinc=-1;
  277.         *ang-=180;
  278.     }
  279.     else if( ( *ang > 270) && ( *ang < 360) ){
  280.         xdir=forward;
  281.