a_comp.c

关于生物光学程序的代码

#include<stdio.h>
#include<math.h>
#include <stdlib.h>
#include "A_comp.h"
#include "allocate.h"

#define LEN 512
#define EVEN 0
#define ODD 1

extern int VERBOSE;

/*#define DERENZO
 */
void A_comp_ij(
	    int im_row, int im_col,
	    struct sino_info *sino,
	    struct geom_info *geom,
	    float **pix_prof,
	    struct column *sp_col);



void pix_prof_comp(
 struct sino_info *sino,
 struct geom_info *geom,
 float ***pix_prof)
{

 int	i,j;
 double pi,ang,d1,d2,t,t_1,t_2,t_3,t_4,maxval,scale,rc;


printf("pix-prof: geom->n_x=%d geom->n_y=%d geom->delta=%f\n",geom->n_x,geom->n_y,geom->delta);


 *pix_prof = (float **)get_img(LEN,sino->n_theta,sizeof(float));

 pi = 4.0 * atan(1.0);

/* Fill profiles of pixels from all angles
 Each profile assumed 2 pixels wide
 This is used to speed calculation (too many to store) of the
        entries of the projection matrix.
 Values are scaled by "scale"
*/
/* WATCH THIS; ONLY FOR SQUARE PIXELS NOW 	*/

	scale = geom->delta; /* = pixel width, for k=0:scale=sino_delta_t */
                             /* k>0: sino_delta_t * 2^k                   */
                             /* (current values by calling program, not nec.*/
	/*
	printf("delta_x = %f, delta_theta = %f\n",scale,sino->delta_theta);*/

    rc = sin((double)(pi/4.0));
    for(i = 0; i < sino->n_theta; i++) {
	 ang = sino->theta_0+i*sino->delta_theta;

	 while(ang >= pi/2.0) /* subtract pi/2 until    ang < pi/2 */
		ang -= pi/2.0;

         if(ang <= pi/4.0)
                maxval = scale/cos((double) ang);
         else
                maxval = scale/cos((double) (pi/2.0-ang));

         d1 = rc * cos((double) (pi/4.0-ang));
         d2 = rc*fabs((double) (sin((double)(pi/4.0-ang))));


	 /* Comments added 3/3/98 by TF, so might be wrong...*/
	 /* t=1.0 is assumed to be the t-value of a projection which would */
	 /* intersect the center of the square pixel.*/
	 /* The profile (for each angle) is 2 pixels wide, that is, it goes from */
	 /* t=0 to t=+2 in steps 2/LEN */
	 /* maxval is the maximum  length of intersection of the pixel with the */
	 /* projection. this happens, when the projection enters and leaves the */
	 /* pixel at opposite sides . Before and after this t range of maximum */
	 /* intersection, the projection enters and leaves on sides which have a 90 degree*/
	 /* angle to each other. */
	 /* t_1 is the t-value where the projection starts intersecting the pixel at all (<1.0) */
	 /* t_2 is the value, where maxval is reached (<=1.0) */
	 /* t_3 is last value for which maxvalue holds (>=1.0) */
	 /* t_4 is t-value >1 where projection does not intersect pixel anymore */

         t_1 = 1.0-d1;
         t_2 = 1.0-d2;
         t_3 = 1.0+d2;
         t_4 = 1.0+d1;

	 /* Once t_i are known, linearly interpolate in range */
	 /* t1..t2 and in t3..t4. Between t2..t3 have maxval. */
	 /* Below t1 and above t4 have no intersection */
         for(j = 0; j < LEN; j++) {
                t = 2.0*j / (double)LEN;
                if(t<=t_1 || t>=t_4)
                        (*pix_prof)[i][j] = 0.0;
                else if(t>t_3)
                        (*pix_prof)[i][j] = maxval*(t_4-t)/(t_4-t_3);
                else if(t>=t_2)
                        (*pix_prof)[i][j] = maxval;
                else
                        (*pix_prof)[i][j] = maxval*(t-t_1)/(t_2-t_1);
         }
    }
	/*
	printf("Finished profiles computation\n");*/


} /* Close prof_comp()	*/


/* Comments added 3/3/98 TF */
/*   Acomp computes one column of the A-matrix, i.e. takes pixel location row,col as*/
/*   input and returns sp_col , where */
/*   int *spcol.row   is a list of the projection indices */ 
/*(proj. in 1D order (angle*n_t + t)) that intersect with the pixel */
/*   float  sp_col.val[i]  contains the length of the intersection with projection spcol.row[i] */
/*  spcol.n_row is number of projections that intersect with pixel */
void A_comp(
	    int im_row, int im_col,
	    struct sino_info *sino,
	    struct geom_info *geom,
	    float **pix_prof,
	    struct column *sp_col)
{

#ifdef STORE_A_MATRIX
  /* precompute entire A-matrix at first call, then 
     just return precomputed cols in subsequent calls */

  static int first=1;
  struct column sp_col_sgl;
  static struct column **sp_col_arr;
  int i,j,r;

  if (first==1)
    {
      first=0;fflush(stderr);
      printf("Precomputing A-matrix...");fflush(stdout);
      sp_col_sgl.val = (float *) get_spc(sino->n_t*sino->n_theta, sizeof(float));
      sp_col_sgl.row = (int *) get_spc(sino->n_t*sino->n_theta, sizeof(int));
      sp_col_sgl.multiplier = 1.0; 
      sp_col_sgl.A_cor = sp_col->A_cor;

      sp_col_arr=(struct column **)multialloc(sizeof(struct column),2,geom->n_y,geom->n_x);


      for (i=0;i<geom->n_y;i++)
	for (j=0;j<geom->n_x;j++)
	  {
	    A_comp_ij(i,j,sino,geom,pix_prof,&sp_col_sgl);
	    sp_col_arr[i][j].val=(float *) get_spc(sp_col_sgl.n_row, sizeof(float));
	    sp_col_arr[i][j].row = (int *) get_spc(sp_col_sgl.n_row, sizeof(int));
	    sp_col_arr[i][j].n_row = sp_col_sgl.n_row;
	    sp_col_arr[i][j].multiplier = sp_col_sgl.multiplier;
	    for (r=0;r<sp_col_sgl.n_row;r++)
	      {
		sp_col_arr[i][j].val[r] = sp_col_sgl.val[r];
		sp_col_arr[i][j].row[r] = sp_col_sgl.row[r];
	      }

	  }
      free((void*)sp_col_sgl.val);
      free((void*)sp_col_sgl.row);
      fflush(stderr);
      printf("done\n");fflush(stdout);
    }


  sp_col->n_row = sp_col_arr[im_row][im_col].n_row;
  sp_col->multiplier = sp_col_arr[im_row][im_col].multiplier;
  for (r=0;r<sp_col_arr[im_row][im_col].n_row;r++)
    {
      sp_col->val[r]=sp_col_arr[im_row][im_col].val[r];
      sp_col->row[r]=sp_col_arr[im_row][im_col].row[r];
    }

#else
  /* compute A_*j on the fly */
  A_comp_ij(im_row,im_col,sino,geom,pix_prof,sp_col);
#endif
}




void A_comp_ij(
	    int im_row, int im_col,
	    struct sino_info *sino,
	    struct geom_info *geom,
	    float **pix_prof,
	    struct column *sp_col)
{
  int	ind,ind_min,ind_max,pr,prof_ind,pind,i,proj_count;
  
  double Aval,t_min,t_max,ang,x,y;
#ifdef DERENZO 
  int stat;
#endif

  /* WATCH THIS; ONLY FOR SQUARE PIXELS NOW 	*/
  
  /*if (sp_col->A_cor.mode != 'e')
    {
    fprintf(stderr,"WARNING: Acomp not in emission mode\n");
    }*/

  
  y = geom->y_0 + im_row*geom->delta;
  x = geom->x_0 + im_col*geom->delta;
  
  proj_count = 0;
  for(pr=0; pr < sino->n_theta; pr++)
    {

#ifdef DERENZO
      if (2*(int)(pr/2)==pr)
	{
	  stat=EVEN;
	  sino->t_0 = -300.4147;
	}
      else
	{
	  stat=ODD;
	  sino->t_0 = -297.2854;
	}
#endif
      
      
      ang = sino->theta_0 + pr*(sino->delta_theta);
      pind = pr*sino->n_t;
      
      /* Need profile to contain 2 pixel widths       */
      /* pixel center is    y*cos((double) ang) - x*sin((double) ang) */
      /* so go down to tmin=center-1  and upto  t_max-center+1 */
      /* consistent with profiles computation t=0..2 */
      t_min = y*cos((double) ang) - x*sin((double) ang) - geom->delta;
      t_max = t_min + 2.0*geom->delta;
      
      /* Check whether pixel is hit at all or t-range below even the first parallel proj. */
      if(t_max < sino->t_0) /* This also prevents over-reach (with rounding of
			       negative numbers)	*/
	continue;
      
      /* From physical t_min and t_max, compute the index range of parallel proj. */
      /* A more general implementation: T.F. May 7, 99 */
      ind_min = ceil((t_min-sino->t_0)/sino->delta_t); 
      ind_max = (t_max-sino->t_0)/sino->delta_t;

      /* Old implementation: Comments by T.F. May 7, 99
	 The following assumes that the reconstruction
	 origin (x,y)=(0,0) lies at 
               0.5*sino->n_t-0.5, 
	 i.e. in the
	 center of the parallel projection array. 
	 Then for general t, we have
              ind= t/sino->delta_t + 0.5*sino->n_t-0.5
	 Now for ind_min, we want the ceil-operator, so add 1,
	 which results in the +0.5 in the end.
	 For ind_max want floor operator, which is implicit
	 in the conversion to int.
	 
	 Notice that this assumes that t_0 is such
	 that the center of the par. proj is
	 at index (n_t-1)/2.

	 It seems to me that Suhails implementation 
	 below assumes that
	     geom->dia=sino->n_t*sino->delta_t;
	 otherwise, I can't make sense out of it, i.e. 
	 dividing t by geom.n_x * geom.delta in general 
	 doesn't make sense to me since t is already in
	 physical coordinates, not in units of geom.delta.
      */
      /* ind_min = (t_min/geom->dia+0.5)*sino->n_t+0.5;
	 ind_max = (t_max/geom->dia+0.5)*sino->n_t-0.5;
      */




      /*printf("x=%f y=%f t_min=%f tmax=%f ind_min=%d ind_max=%d ind_min*delta_t+t_0=%f ind_max*delta_t+t_0=%f\n",x,y,t_min,t_max,ind_min,ind_max,ind_min*sino->delta_t+sino->t_0,ind_max*sino->delta_t+sino->t_0);*/
   

      /* Fix this 4/91 to prevent over-reach at ends  */
      /* clip ends of t-range */
      ind_min = (ind_min<0) ? 0:ind_min;
      ind_max = (ind_max>=sino->n_t) ? sino->n_t-1:ind_max;
      for(i = ind_min; i <=ind_max; i++) 
	{
	  
	  ind = pind+i; /* linear index in 1D array of all proj. (angle*n_t + t)*/
	  
	  /* get profile_index that is closest to this projection */
	  /* profiles are in discrete t-steps of 2pixels/LEN = 2.0*geom->delta / LEN */
	  prof_ind = LEN*(sino->t_0+i*sino->delta_t-t_min)/(2.0*geom->delta);
	  
	  if(prof_ind >= LEN || prof_ind <0)
	    {
	      /*
		printf("Error in index computation for profiles %d\n",prof_ind);
		*/
	      if(prof_ind == LEN)
		prof_ind = LEN-1;
	      else if (prof_ind == -1)
		prof_ind = 0;
	      else
		{
		  printf("Exiting w/ error in Acomp line 127.\nprof_ind=%d out of range 0..%d",prof_ind,LEN-1);
		  exit(-1);
		}
	    }
	  /* Assign output */
	  if((Aval=pix_prof[pr][prof_ind]) > 0.0) 
            {
              if (sp_col->A_cor.mode=='t')
	        sp_col->val[proj_count] = Aval;
              else
                if (sp_col->A_cor.mode=='e') 
		  /* in emis-mode, divide by forward proj. trn-sinogram */
		  /* to account for attenuation. */            
		  {
		    sp_col->val[proj_count] = Aval*(sp_col->A_cor.Y_t[pr][i]);
		  }		
		else
		  {
		    fprintf(stderr,"A_comp Error: Unknown mode %c\n\n",sp_col->A_cor.mode);
		    exit(-1);
		  }
	      sp_col->row[proj_count] = ind;/* Index for ray    */
	      proj_count++;
	    }
	}
    } 
  sp_col->n_row = proj_count;
  /*
    printf("proj_count = %d\n",proj_count);
    */
  
  
}/* Close A_comp()	*/

#undef LEN 

#include <mat.h>


#ifdef MATLAB5

#include <matrix.h>

void read_mat_sino_info(
  MATFile *fin,       /* pointer to file */
  struct sino_info *sino
  )
{
 mxArray *pm;

 pm = matGetArray(fin,"n_theta");
 if(pm != NULL){
  sino->n_theta = (int) mxGetScalar(pm);
  }
 mxDestroyArray(pm);

 pm = matGetArray(fin,"n_t");
 if(pm != NULL){
  sino->n_t = (int) mxGetScalar(pm);
  }
 mxDestroyArray(pm);

 pm = matGetArray(fin,"delta_theta");
 if(pm != NULL){
  sino->delta_theta = mxGetScalar(pm);
  }
 mxDestroyArray(pm);

 pm = matGetArray(fin,"delta_t");
 if(pm != NULL){
  sino->delta_t = mxGetScalar(pm);
  }
 mxDestroyArray(pm);

 pm = matGetArray(fin,"theta_0");
 if(pm != NULL){
  sino->theta_0 = mxGetScalar(pm);
  }
 mxDestroyArray(pm);

 pm = matGetArray(fin,"t_0");
 if(pm != NULL){
  sino->t_0 = mxGetScalar(pm);
  }
 mxDestroyArray(pm);

 sino->dec_a = 1;
 sino->dec_t = 1;

 if(VERBOSE)
 {
 printf("delta_t = %f\n",sino->delta_t);
 printf("delta_theta = %f\n",sino->delta_theta);
 printf("t_0 = %f\n",sino->t_0);
 printf("theta_0 = %f\n",sino->theta_0);
 printf("n_t = %d\n",sino->n_t);
 printf("n_theta = %d\n",sino->n_theta);
 }

}

#else

#include <mat.h>

void read_mat_sino_info(
  MATFile *fin,       /* pointer to file */
  struct sino_info *sino
  )
{
 Matrix *pm;

 pm = matGetMatrix(fin,"n_theta");
 if(pm != NULL){
  sino->n_theta = (int) mxGetScalar(pm);
  }
 pm = matGetMatrix(fin,"n_t");
 if(pm != NULL){
  sino->n_t = (int) mxGetScalar(pm);
  }
 pm = matGetMatrix(fin,"delta_theta");
 if(pm != NULL){
  sino->delta_theta = mxGetScalar(pm);
  }
 pm = matGetMatrix(fin,"delta_t");
 if(pm != NULL){
  sino->delta_t = mxGetScalar(pm);
  }
 pm = matGetMatrix(fin,"theta_0");
 if(pm != NULL){
  sino->theta_0 = mxGetScalar(pm);
  }
 pm = matGetMatrix(fin,"t_0");
 if(pm != NULL){
  sino->t_0 = mxGetScalar(pm);
  }

 sino->dec_a = 1;
 sino->dec_t = 1;

 if(VERBOSE)
 {
 printf("delta_t = %f\n",sino->delta_t);
 printf("delta_theta = %f\n",sino->delta_theta);
 printf("t_0 = %f\n",sino->t_0);
 printf("theta_0 = %f\n",sino->theta_0);
 printf("n_t = %d\n",sino->n_t);
 printf("n_theta = %d\n",sino->n_theta);
 }

}

#endif