atnmake.c File Reference

Compute PET attenuation correction data from blank and transmission data for ECAT 931 and GE Advance 2D sinograms in ECAT 6.3 file format. More...

Go to the source code of this file.

Functions
int	atnMake (char *trafile, char *blkfile, char *norfile, char *atnfile, char *imgfile, int decay, int limit, int keepNegat, int imgDim, float zoom, float shiftX, float shiftY, float rotation, int maxIterNr, int skipPriorNr, float beta, int maskDim, int osSetNr, char *status, int verbose)

Detailed Description

Compute PET attenuation correction data from blank and transmission data for ECAT 931 and GE Advance 2D sinograms in ECAT 6.3 file format.

Updated from program written by Jarkko Rintaluoma March 1995, and later modified by Sakari Alenius to apply MRP method, with further editions by Vesa Oikonen and Petri Numminen.

Copyright

(c) Turku PET Centre

Author

Vesa Oikonen

Definition in file atnmake.c.

Function Documentation

atnMake()

int atnMake	(	char *	trafile,
		char *	blkfile,
		char *	norfile,
		char *	atnfile,
		char *	imgfile,
		int	decay,
		int	limit,
		int	keepNegat,
		int	imgDim,
		float	zoom,
		float	shiftX,
		float	shiftY,
		float	rotation,
		int	maxIterNr,
		int	skipPriorNr,
		float	beta,
		int	maskDim,
		int	osSetNr,
		char *	status,
		int	verbose )

Read transmission scan file and blank scan file and calculate the attenuation correction data and the logarithmic correction data.

See also

trmrp

Returns

Returns 0 if ok.

Parameters

trafile	Pointer to transmission sinogram filename.
blkfile	Pointer to blank sinogram filename.
norfile	Pointer to normalization sinogram filename.
atnfile	Pointer to filename for the attenuation data to be written.
imgfile	Pointer to filename for the transmission image to be written; enter NULL or empty string if not needed.
decay	Decay correction for scan time difference: 0=no, 1=yes.
limit	Limit: 0=no; 1=attenuation factors are set between 1.0 and ATN_HI_MAX.
keepNegat	Keep negative sinogram values (1), or set negative sinogram values to zero (0), like old program scnnorm did.
imgDim	Image dimension (size, usually 128 or 256); must be an even number.
zoom	Zoom factor (for example 2.45 for the brain); 1=no zooming.
shiftX	Possible shifting in x dimension (mm).
shiftY	Possible shifting in y dimension (mm).
rotation	Possible image rotation, -180 - +180 (in degrees).
maxIterNr	Nr of iterations, for example 150.
skipPriorNr	Number of iterations to skip before prior; usually 1.
beta	Beta, 0.01 - 0.9; usually 0.3 for emission, 0.9 for transmission.
maskDim	Median filter mask dimension; 3 or 5 (9 or 21 pixels).
osSetNr	Number of Ordered Subset sets; 1, 2, 4, ... 128.
status	Pointer to a string (allocated for at least 64 chars) where error message or other execution status will be written; enter NULL, if not needed
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 115 of file atnmake.c.

  {
  int ret;
 
  /* Check the input */
  if(status!=NULL) sprintf(status, "invalid filename");
  if(trafile==NULL || !trafile[0]) return(1);
  if(blkfile==NULL || !blkfile[0]) return(1);
  if(norfile==NULL || !blkfile[0]) return(1);
  if(atnfile==NULL || !atnfile[0]) return(1);
 
 
  /* 
   *  Open the transmission, blank, and normalization sinograms for reading
   */
  FILE *fptra=NULL, *fpblk=NULL, *fpnor=NULL;
  if(verbose>0) printf("opening %s\n", trafile);
  if((fptra=fopen(trafile, "rb")) == NULL) {
    if(status!=NULL) sprintf(status, "cannot open transmission sinogram"); 
    return(2);
  }
  if(verbose>0) printf("opening %s\n", blkfile);
  if((fpblk=fopen(blkfile, "rb")) == NULL) {
    if(status!=NULL) sprintf(status, "cannot open blank sinogram"); 
    fclose(fptra);
    return(2);
  }
  if(verbose>0) printf("opening %s\n", norfile);
  if((fpnor=fopen(norfile, "rb")) == NULL) {
    if(status!=NULL) sprintf(status, "cannot open normalization file"); 
    fclose(fptra); fclose(fpblk);
    return(2);
  }
 
  /* 
   *  Read main headers
   */
  if(verbose>1) printf("reading transmission main header\n");
  ECAT63_mainheader tra_main_header;
  if((ret=ecat63ReadMainheader(fptra, &tra_main_header))) {
    if(status!=NULL) sprintf(status, "cannot read transmission main header");
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
  if(verbose>100) ecat63PrintMainheader(&tra_main_header, stdout);
  if(verbose>4) printf("file_type := %d\n", tra_main_header.file_type);
 
  if(verbose>1) printf("reading blank main header\n");
  ECAT63_mainheader blk_main_header;
  if((ret=ecat63ReadMainheader(fpblk, &blk_main_header))) {
    if(status!=NULL) sprintf(status, "cannot read blank main header");
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
  if(verbose>100) ecat63PrintMainheader(&blk_main_header, stdout);
  if(verbose>4) printf("file_type := %d\n", blk_main_header.file_type);
 
  if(verbose>1) printf("reading normalization main header\n");
  ECAT63_mainheader nor_main_header;
  if((ret=ecat63ReadMainheader(fpnor, &nor_main_header))) {
    if(status!=NULL) sprintf(status, "cannot read blank main header");
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
  if(verbose>100) ecat63PrintMainheader(&nor_main_header, stdout);
  if(verbose>4) printf("file_type := %d\n", nor_main_header.file_type);
 
  /* Check the file_type */
  if(tra_main_header.file_type!=RAW_DATA || blk_main_header.file_type!=RAW_DATA
     || nor_main_header.file_type!=RAW_DATA) 
  {
    if(status!=NULL) sprintf(status, "file is not a sinogram");
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
 
  /* Read scan start times */
  time_t tra_start, blk_start;
  tra_start=ecat63Scanstarttime(&tra_main_header);
  blk_start=ecat63Scanstarttime(&blk_main_header);
  if(verbose>1) {
    char buf[32];
    printf("Blank scan start time := %s\n", ctime_r_int(&blk_start, buf));
    printf("Transmission scan start time := %s\n", ctime_r_int(&tra_start, buf));
  }
  if(decay!=0 && (tra_start==0 || blk_start==0)) {
    if(status!=NULL) sprintf(status, "missing scan start time");
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
  double start_time_diff=difftime(tra_start, blk_start);
  if(verbose>1) {printf("scan start time difference := %g\n", start_time_diff);}
 
  /* Get isotope halflife */
  double halflife=-1.0;
  if(tra_main_header.isotope_halflife>0.0) halflife=tra_main_header.isotope_halflife;
  else if(blk_main_header.isotope_halflife>0.0) halflife=blk_main_header.isotope_halflife;
  if(verbose>1 && halflife>0.0) {printf("isotope halflife := %g\n", halflife);}
  if(decay!=0 && halflife<=0.0) {
    if(status!=NULL) sprintf(status, "missing isotope halflife");
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
 
 
 
  /*
   *  Read the matrix lists
   */
  if(verbose>1) printf("reading transmission matrix list\n");
  static MATRIXLIST tra_mlist; ecat63InitMatlist(&tra_mlist);
  ret=ecat63ReadMatlist(fptra, &tra_mlist, verbose-2);
  if(ret) {
    if(status!=NULL) sprintf(status, "cannot read sinogram matrix list");
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
 
  if(verbose>1) printf("reading blank matrix list\n");
  static MATRIXLIST blk_mlist; ecat63InitMatlist(&blk_mlist);
  ret=ecat63ReadMatlist(fpblk, &blk_mlist, verbose-2);
  if(ret) {
    if(status!=NULL) sprintf(status, "cannot read blank matrix list");
    ecat63EmptyMatlist(&tra_mlist);
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
 
  if(verbose>1) printf("reading normalization matrix list\n");
  static MATRIXLIST nor_mlist; ecat63InitMatlist(&nor_mlist);
  ret=ecat63ReadMatlist(fpnor, &nor_mlist, verbose-2);
  if(ret) {
    if(status!=NULL) sprintf(status, "cannot read normalization matrix list");
    ecat63EmptyMatlist(&tra_mlist); ecat63EmptyMatlist(&blk_mlist);
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
 
  /* Check the contents of matrix lists */
  if(tra_mlist.matrixNr<=0 || blk_mlist.matrixNr<=0 || nor_mlist.matrixNr<=0) {
    if(status!=NULL) sprintf(status, "empty matrix list");
    ecat63EmptyMatlist(&tra_mlist); ecat63EmptyMatlist(&blk_mlist); 
    ecat63EmptyMatlist(&nor_mlist);
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
  if(ecat63CheckMatlist(&tra_mlist) || ecat63CheckMatlist(&blk_mlist) || 
     ecat63CheckMatlist(&nor_mlist)) 
  {
    if(status!=NULL) sprintf(status, "invalid matrix list");
    ecat63EmptyMatlist(&tra_mlist); ecat63EmptyMatlist(&blk_mlist);
    ecat63EmptyMatlist(&nor_mlist);
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
  if(verbose>3) {
    printf("transmission_matrixNr := %d\n", tra_mlist.matrixNr);
    printf("blank_matrixNr := %d\n", blk_mlist.matrixNr);
    printf("normalization_matrixNr := %d\n", nor_mlist.matrixNr);
  }
  if(tra_mlist.matrixNr!=blk_mlist.matrixNr ||
     tra_mlist.matrixNr!=nor_mlist.matrixNr ||
     tra_main_header.num_planes != blk_main_header.num_planes ||
     tra_main_header.num_planes != nor_main_header.num_planes)
  {
    if(status!=NULL) sprintf(status, "different plane numbers");
    ecat63EmptyMatlist(&tra_mlist); ecat63EmptyMatlist(&blk_mlist);
    ecat63EmptyMatlist(&nor_mlist);
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(3);
  }
 
  /*
   *  Open output file(s) and write main header(s)
   */
  if(verbose>0) printf("writing main header in %s\n", atnfile);
  ECAT63_mainheader atn_main_header;
  ecat63CopyMainheader(&tra_main_header, &atn_main_header);
  atn_main_header.data_type = IEEE_R4;    /* float */
  atn_main_header.file_type = ATTN_DATA;
  atn_main_header.calibration_units = 2;
  strcpy(atn_main_header.user_process_code, "tMRP");
  FILE *fpatn=NULL;
  fpatn=ecat63Create(atnfile, &atn_main_header);
  if(fpatn==NULL) {
    if(status!=NULL) sprintf(status, "cannot write attenuation file\n"); 
    ecat63EmptyMatlist(&tra_mlist); ecat63EmptyMatlist(&blk_mlist);
    ecat63EmptyMatlist(&nor_mlist);
    fclose(fptra); fclose(fpblk); fclose(fpnor);
    return(11);
  }
 
  FILE *fpimg=NULL;
  if(imgfile!=NULL && imgfile[0]) {
    if(verbose>0) printf("writing main header in %s\n", imgfile);
    ECAT63_mainheader img_main_header;
    ecat63CopyMainheader(&tra_main_header, &img_main_header);
    img_main_header.file_type = IMAGE_DATA;
    img_main_header.data_type = SUN_I2;    /* short int */
    img_main_header.calibration_units = 2;
    img_main_header.num_frames=1;
    strcpy(img_main_header.user_process_code, "tMRP");
    fpimg=ecat63Create(imgfile, &img_main_header);
    if(fpimg==NULL) {
      if(status!=NULL) sprintf(status, "cannot write transmission image\n"); 
      ecat63EmptyMatlist(&tra_mlist); ecat63EmptyMatlist(&blk_mlist);
      ecat63EmptyMatlist(&nor_mlist);
      fclose(fptra); fclose(fpblk); fclose(fpnor);
      remove(atnfile);
      return(12);
    }
  }
 
 
  /*
   *  Process one matrix at a time
   */
  if(verbose>0) printf("processing matrices...\n");
  int cret=0;
  for(int mi=0; mi<tra_mlist.matrixNr; mi++) {
    if(verbose==1) {fprintf(stdout, "."); fflush(stdout);}
 
    FILE *lfptra=fptra;
    FILE *lfpblk=fpblk;
    FILE *lfpnor=fpnor;
    FILE *lfpatn=fpatn;
    FILE *lfpimg=fpimg;
 
    ECAT63_scanheader tra_header, blk_header, nor_header;
    ECAT63_attnheader atn_header; 
    ECAT63_imageheader img_header; 
    Matval matval;
    int dimx, dimy, scnpxlNr=0, imgpxlNr=0, matnum=0, ni=0, bi=0;
    float *tramat, *blkmat, *normat, f;
    int ret=0;
 
    /* Get plane and frame nr */
    mat_numdoc(tra_mlist.matdir[mi].matnum, &matval);
 
    /* Read transmission sinogram subheader and scaled data */
    if(verbose>1) printf("reading transmission matrix %d\n", 1+mi);
    ret=ecat63ReadScanMatrix(lfptra, tra_mlist.matdir[mi].strtblk, 
          tra_mlist.matdir[mi].endblk, &tra_header, &tramat);
    if(ret) {
      if(status!=NULL) 
        sprintf(status, "cannot read transmission sinogram matrix %u", tra_mlist.matdir[mi].matnum);
      cret=ret;
      continue;
    }
    if(verbose>2)
      printf("Matrix: plane %d frame %d gate %d bed %d\n",
             matval.plane, matval.frame, matval.gate, matval.bed);
    if(verbose>200) {
      ecat63PrintScanheader(&tra_header, stdout);
    } else if(verbose>1 && mi==0) {
      printf("Transmission sinogram dimensions := %d x %d\n",
             tra_header.dimension_1, tra_header.dimension_2);
      printf("Transmission sinogram frame duration := %g\n",
             0.001*(double)tra_header.frame_duration);
      printf("Transmission sinogram dead-time correction := %g\n",
             tra_header.loss_correction_fctr);
    }
    dimx=tra_header.dimension_1;
    dimy=tra_header.dimension_2;
    scnpxlNr=dimx*dimy;
    if(imgDim<2) imgDim=dimx; // transmission dim to number of rays by default
    /* Dead-time correction */
    if(tra_header.loss_correction_fctr>0.0) {
      for(int i=0; i<scnpxlNr; i++) 
        tramat[i]*=tra_header.loss_correction_fctr;
    }
    /* Divide counts by frame duration */
    if(tra_header.frame_duration<1) {
      if(status!=NULL) sprintf(status, "missing frame duration in transmission sinogram");
      free(tramat);
      cret=ret;
      continue;
    }
    f=1000.0/tra_header.frame_duration;
    for(int i=0; i<scnpxlNr; i++) tramat[i]*=f; 
 
    /* Search corresponding matrix from the normalization file */
    matnum=mat_numcod(matval.frame, matval.plane, matval.gate, matval.data, matval.bed);
    for(ni=0; ni<nor_mlist.matrixNr; ni++)
      if(nor_mlist.matdir[ni].matnum==matnum)
        break;
    if(ni>=nor_mlist.matrixNr) {
      if(status!=NULL) 
        sprintf(status, "cannot find matrix %u in normalization sinogram",
                nor_mlist.matdir[ni].matnum);
      free(tramat);
      cret=ret;
      continue;
    }
    /* Read normalization subheader and scaled data */
    if(verbose>1) printf("reading normalization matrix %d\n", 1+ni);
    ret=ecat63ReadScanMatrix(lfpnor, nor_mlist.matdir[ni].strtblk, 
          nor_mlist.matdir[ni].endblk, &nor_header, &normat);
    if(ret) {
      if(status!=NULL) 
        sprintf(status, "cannot read blank sinogram matrix %u", nor_mlist.matdir[ni].matnum);
      free(tramat);
      cret=ret;
      continue;
    }
    if(verbose>200) {
      ecat63PrintScanheader(&nor_header, stdout);
    } else if(verbose>1 && mi==0) {
      printf("Normalization sinogram dimensions := %d x %d\n",
             nor_header.dimension_1, nor_header.dimension_2);
      printf("Normalization sinogram frame duration := %g\n",
             0.001*(double)nor_header.frame_duration);
      printf("Normalization sinogram dead-time correction := %g\n",
             nor_header.loss_correction_fctr);
    }
    if(tra_header.dimension_1!=nor_header.dimension_1
       || tra_header.dimension_2!=nor_header.dimension_2)
    {
      if(status!=NULL) sprintf(status, "incompatible matrix dimensions");
      free(tramat); free(normat);
      cret=ret;
      continue;
    }
 
    /* Normalization correction for transmission sinogram */
    for(int i=0; i<scnpxlNr; i++) tramat[i]*=normat[i];
    /* Set negatives to zero */
    if(keepNegat==0) 
      for(int i=0; i<scnpxlNr; i++) if(tramat[i]<0.0) tramat[i]=0.0;
 
    /* Search corresponding matrix from the blank file */
    matnum=mat_numcod(matval.frame, matval.plane, matval.gate, matval.data, 
                      matval.bed);
    for(bi=0; bi<blk_mlist.matrixNr; bi++)
      if(blk_mlist.matdir[bi].matnum==matnum)
        break;
    if(bi>=blk_mlist.matrixNr) {
      if(status!=NULL) 
        sprintf(status, "cannot find matrix %u in blank sinogram", blk_mlist.matdir[bi].matnum);
      free(tramat); free(normat);
      cret=ret;
      continue;
    }
    /* Read blank subheader and scaled data */
    if(verbose>1) printf("reading blank matrix %d\n", 1+bi);
    ret=ecat63ReadScanMatrix(lfpblk, blk_mlist.matdir[bi].strtblk, 
          blk_mlist.matdir[bi].endblk, &blk_header, &blkmat);
    if(ret) {
      if(status!=NULL) 
        sprintf(status, "cannot read blank sinogram matrix %u", blk_mlist.matdir[bi].matnum);
      free(tramat); free(normat);
      cret=ret;
      continue;
    }
    if(verbose>200) {
      ecat63PrintScanheader(&blk_header, stdout);
    } else if(verbose>1 && mi==0) {
      printf("Blank sinogram dimensions := %d x %d\n",
             blk_header.dimension_1, blk_header.dimension_2);
      printf("Blank sinogram frame duration := %g\n",
             0.001*(double)blk_header.frame_duration);
      printf("Blank sinogram dead-time correction := %g\n",
             blk_header.loss_correction_fctr);
    }
    if(tra_header.dimension_1!=blk_header.dimension_1
       || tra_header.dimension_2!=blk_header.dimension_2)
    {
      if(status!=NULL) sprintf(status, "incompatible matrix dimensions");
      free(tramat); free(normat); free(blkmat);
      cret=ret;
      continue;
    }
    /* Dead-time correction */
    if(blk_header.loss_correction_fctr>0.0) {
      for(int i=0; i<scnpxlNr; i++) 
        blkmat[i]*=blk_header.loss_correction_fctr;
    }
    /* Divide counts by frame duration */
    if(blk_header.frame_duration<1) {
      if(status!=NULL) sprintf(status, "missing frame duration in blank sinogram");
      free(tramat); free(normat); free(blkmat);
      cret=ret;
      continue;
    }
    f=1000.0/blk_header.frame_duration;
    for(int i=0; i<scnpxlNr; i++) blkmat[i]*=f; 
 
    /* Normalization correction for blank sinogram */
    for(int i=0; i<scnpxlNr; i++) blkmat[i]*=normat[i];
    /* Set negatives to zero */
    if(keepNegat==0) for(int i=0; i<scnpxlNr; i++) if(blkmat[i]<0.0) blkmat[i]=0.0;
 
 
    /* Correct transmission scan for the isotope decay in radiation source
       that may have happened between blank and transmission scans;
       it is not important in patient scans, but may be very important for
       phantom studies and calibration measurements.
       Frame scan durations are considered to be insignificant for decay.
     */
    double dcf=1.0;
    if(decay!=0) {
      double lambda=hl2lambda(halflife);
      double t=start_time_diff;
      t+=0.001*(double)tra_header.frame_start_time;
      t-=0.001*(double)blk_header.frame_start_time;
      dcf=hlLambda2factor(lambda, t, -1.0);
      if(verbose>1 && mi==0) {
        printf("frame start time difference := %g\n", t);
        printf("decay correction factor := %g\n", dcf);
      }
      for(int i=0; i<scnpxlNr; i++) tramat[i]*=dcf; 
    }
 
 
    /* Set up the attenuation sub header information */
    if(verbose>3) printf("setting attenuation sub-header\n");
    if(mi==0) {
      /* Set fields that are common for all matrices */
      memset(&atn_header, 0, sizeof(ECAT63_attnheader));
      atn_header.data_type=IEEE_R4;   /* float */
      atn_header.attenuation_type=1;  /* 1 for measured attenuation */
      atn_header.scale_factor=1.0; 
      atn_header.x_origin=0.0;
      atn_header.y_origin=0.0;
      atn_header.x_radius=0.0;
      atn_header.y_radius=0.0;
      atn_header.tilt_angle=tra_main_header.gantry_tilt; 
      atn_header.attenuation_coeff=1.0;
      atn_header.sample_distance=tra_header.sample_distance;
    }
    atn_header.dimension_1 = tra_header.dimension_1; 
    atn_header.dimension_2 = tra_header.dimension_2; 
 
    /* Set up the log attenuation sub header information, if needed */
    if(lfpimg!=NULL && mi==0) {
      if(verbose>6) printf("setting transmission image sub-header\n");
      /* Set fields that are common for all matrices */
      img_header.data_type=SUN_I2;   /* short int */
      img_header.num_dimensions=2;
      img_header.dimension_1=imgDim;
      img_header.dimension_2=imgDim;
      img_header.x_origin=10.0*shiftX;
      img_header.y_origin=10.0*shiftY;
      img_header.recon_scale=zoom;
      img_header.quant_scale=1.0;
      img_header.pixel_size=tra_header.sample_distance*(float)dimx/(zoom*(float)imgDim);
      img_header.slice_width=tra_main_header.plane_separation;
      img_header.frame_start_time=0;
      img_header.frame_duration=0; // to prevent accidental decay correction
      img_header.frame_duration=0;
      img_header.filter_code=3;
      img_header.image_rotation=rotation;
      img_header.decay_corr_fctr=dcf;
      img_header.loss_corr_fctr=1.0;
      img_header.quant_units=2;
      img_header.ecat_calibration_fctr=1.0;
      img_header.filter_params[0]=beta;
      img_header.filter_params[1]=(float)maskDim;
      img_header.filter_params[2]=(float)maxIterNr;
      img_header.filter_params[3]=(float)1;
      img_header.filter_params[4]=(float)osSetNr;
      img_header.filter_params[5]=(float)skipPriorNr;
      strcpy(img_header.annotation, "tMRP");
    }
    if(lfpimg!=NULL) { // separately for each matrix
      img_header.plane_eff_corr_fctr=nor_header.scale_factor;
      img_header.scan_matrix_num=tra_mlist.matdir[mi].matnum;
      img_header.norm_matrix_num=nor_mlist.matdir[ni].matnum;
      img_header.atten_cor_mat_num=tra_mlist.matdir[mi].matnum;
    }
 
 
    /*
     *  Reconstruct image containing log-transformed attenuation correction factors.
     *  It must be log-transformed, because actual attenuation factors are >1 for the whole sinogram.
     *  At this step there must be no zoom, shifts, or rotation, and
     *  image dimension must be the same as sinogram width (nr or rays or bins).
     */
    if(verbose>1) {printf("matrix reconstruction\n"); fflush(stdout);}
    imgpxlNr=dimx*dimx;
    float imgmat[imgpxlNr];
    ret=trmrp(blkmat, tramat, dimx, imgmat, 120, 1,  
              dimx, dimy, maskDim, 1.0, beta, 
              10.0*tra_main_header.axial_fov, 10.0*tra_header.sample_distance, 
              1, 1, 0.0, 0.0, 0.0,
              verbose-6);
    if(ret) {
      if(status!=NULL) sprintf(status, "cannot calculate attenuation correction factors");
      if(verbose>1) printf("trmrp_return_value := %d\n", ret);
      free(tramat); free(normat); free(blkmat);
      cret=ret;
      continue;
    }
    if(verbose>1) {
      for(int i=0; i<imgpxlNr; i++) 
        if(!isfinite(imgmat[i])) {
          printf("  inf in the image!\n");
          break;
        }
    }
 
 
    /*
     *  If requested transmission image dimension is the same as scan ray number, and
     *  zoom, shifts or rotation were not set, then simply save the reconstructed image;
     *  otherwise we have to reconstruct it later again.
     */
    if(lfpimg!=NULL && imgDim==dimx && zoom==1.0 && shiftX==0.0 && shiftY==0.0 && rotation==0.0) {
      if(verbose>1) printf("writing transmission image matrix\n");
      matnum=mat_numcod(1, matval.plane, matval.gate, matval.data, matval.bed);
      ret=ecat63WriteImageMatrix(lfpimg, matnum, &img_header, imgmat);
      if(ret) {
        if(status!=NULL) sprintf(status, "cannot write transmission image matrix");
        if(verbose>1) printf("ret := %d\n", ret);
        free(tramat); free(normat); free(blkmat);
        cret=ret;
        continue;
      }
    }
 
    /*
     *  Reprojection to sinogram space
     */
    if(verbose>1) {printf("matrix reprojection\n"); fflush(stdout);}
    float atnmat[scnpxlNr];
    ret=reprojection(imgmat, dimx, dimx, dimy, 1.0, atnmat, verbose-10);
    if(ret) {
      if(status!=NULL) sprintf(status, "cannot reproject attenuation correction data");
      if(verbose>1) printf("trmrp_return_value := %d\n", ret);
      free(tramat); free(normat); free(blkmat);
      cret=ret;
      continue;
    }
    if(verbose>1) {
      for(int i=0; i<scnpxlNr; i++) 
        if(!isfinite(atnmat[i])) {
          printf("  inf in the log attenuation data!\n");
          break;
        }
    }
 
    /*
     *  Convert from log values to correction factors for attenuation file
     */
    if(verbose>2) {printf("conversion from logarithms\n"); fflush(stdout);}
#if(0)
    if(tra_header.sample_distance>0.001) f=10.0*tra_header.sample_distance;
    else f=1.0;
    for(int i=0; i<scnpxlNr; i++) atnmat[i]=expf(atnmat[i]*f);
#endif
    for(int i=0; i<scnpxlNr; i++) atnmat[i]=expf(atnmat[i]);
    if(verbose>1) {
      for(int i=0; i<scnpxlNr; i++) 
        if(!isfinite(atnmat[i])) {
          printf("  inf in the attenuation data!\n");
          break;
        }
    }
    for(int i=0; i<scnpxlNr; i++) if(!isfinite(atnmat[i])) atnmat[i]=1.0;
 
    /* Fix too low and high values, if requested */
    if(limit>0) {
      for(int i=0; i<scnpxlNr; i++) {
        if(!(atnmat[i]>=1.0)) atnmat[i]=1.0;
        else if(atnmat[i]>ATN_HI_MAX) atnmat[i]=ATN_HI_MAX;
      }
    }
 
 
    /* Write the attenuation matrix */
    if(verbose>1) printf("writing attenuation matrix\n");
    matnum=mat_numcod(1, matval.plane, matval.gate, matval.data, matval.bed);
    ret=ecat63WriteAttn(lfpatn, matnum, &atn_header, atnmat);
    if(ret) {
      if(status!=NULL) sprintf(status, "cannot write attenuation matrix");
      if(verbose>1) printf("ret := %d\n", ret);
      free(tramat); free(normat); free(blkmat);
      cret=ret;
      continue;
    }
 
 
    /*
     *  Reconstruct and save transmission image, in case user requested it, and
     *  user gave zoom, shifts, rotation, or image dimensions that is different
     *  from the sinogram width (nr or rays or bins).
     */
    if(lfpimg!=NULL && (imgDim!=dimx || zoom!=1.0 || shiftX!=0.0 || shiftY!=0.0 || rotation!=0.0)) {
      if(verbose>1) {printf("transmission matrix reconstruction\n"); fflush(stdout);}
      imgpxlNr=imgDim*imgDim;
      float imgmat2[imgpxlNr];
      ret=trmrp(blkmat, tramat, imgDim, imgmat2, maxIterNr, osSetNr,  
                dimx, dimy, maskDim, zoom, beta, 
                10.0*tra_main_header.axial_fov, 10.0*tra_header.sample_distance, 
                skipPriorNr, 1, shiftX, shiftY, rotation,
                verbose-7);
      if(ret) {
        if(status!=NULL) sprintf(status, "cannot reconstruct transmission image");
        if(verbose>1) printf("trmrp_return_value := %d\n", ret);
      } else {
        f=10.0*tra_header.sample_distance;
        for(int i=0; i<imgpxlNr; i++) imgmat2[i]*=f;
        if(verbose>1) printf("writing transmission image matrix\n");
        matnum=mat_numcod(1, matval.plane, matval.gate, matval.data, matval.bed);
        ret=ecat63WriteImageMatrix(lfpimg, matnum, &img_header, imgmat2);
        if(ret) {
          if(status!=NULL) sprintf(status, "cannot write transmission image matrix");
          if(verbose>1) printf("ret := %d\n", ret);
        }
      }
      if(ret) {
        free(tramat); free(normat); free(blkmat);
        cret=ret;
        continue;
      }
    }
 
    free(tramat); free(normat); free(blkmat);
  } /* next matrix */
  if(verbose==1) {fprintf(stdout, "\n"); fflush(stdout);}
  if(verbose>0) {fprintf(stdout, "done.\n"); fflush(stdout);}
 
 
  fclose(fpatn); if(fpimg!=NULL) fclose(fpimg);
  ecat63EmptyMatlist(&tra_mlist); ecat63EmptyMatlist(&blk_mlist);
  ecat63EmptyMatlist(&nor_mlist);
  fclose(fptra); fclose(fpblk); fclose(fpnor);
  if(cret) {remove(atnfile); if(fpimg!=NULL) remove(imgfile);}
  return(cret);
}