atnmake.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <unistd.h>
#include <ctype.h>
#include <dirent.h>
#include <sys/file.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
#include "libtpcrec.h"
/*****************************************************************************/
#ifdef HAVE_OMP_H
#include <omp.h>
#endif
/*****************************************************************************/
#define ATN_HI_MAX 100.0
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Computing PET attenuation correction data from blank and transmission",
  "sinograms (1) in ECAT 6.3 file format, to be applied to ECAT 931 and",
  "GE Advance 2D sinograms.",
  "Normalization file is required to reduce the noise in attenuation data,",
  "and emission sinogram still needs to be normalized.",
  "Transmission image can be reconstructed and saved optionally.",
  " ",
  "Usage: @P [Options] blankscan trscan normfile atnfile [trimage]",
  " ",
  "Options:",
  " -decay=<Y|n>",
  "     Decay correction for the time difference between blank and",
  "     transmission scans is done (y, default), or not done (n).",
  " -limit=<y|N>",
  "     Attenuation factors are set to values between 1.0 and 100.0 (y),",
  "     or not limited (n, default).",
  " -zoom=<value>",
  "     Set zoom factor for transmission image; by default 1.0 (no zoom).",
  " -dim=<value>",
  "     Set transmission image x and y dimensions; by default the ray number.",
  " -rot[ation]=<value>",
  "     Set transmission image rotation in degrees, -180 - 180; by default 0.",
  " -x=<value>",
  "     Set transmission image shift in x direction (in cm); by default 0.",
  " -y=<value>",
  "     Set transmission image shift in y direction (in cm); by default 0.",
  " -beta=<value>",
  "     Set beta value, usually 0.1 - 0.9; by default 0.9; affects also",
  "     attenuation factors.",
  " -mask=<3|5>",
  "     Set mask dimension for median filter, either 3 (3x3) or 5 (5x5 without",
  "     corner pixels); by default 3. Affects also attenuation factors.",
  " -iter=<value>",
  "     Set maximum number of iterations for transmission image reconstruction;",
  "     by default 45. For attenuation 120 iterations are always used.",
  " -skip=<value>",
  "     Set the number of iterations without prior in transmission image",
  "     reconstruction; by default 1.",
  " -os=<1|2|4|8|16>",
  "     Set the number of OS sets (acceleration) in transmission image",
  "     reconstruction; by default 1.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "  @P 28oct97bl1.scn s02892tr1.scn s02892sta.nrm s02892tr1.atn s02892tr1.img",
  " ",
  "References:",
  "1. Bettinardi V, Alenius S, Numminen P, Teras M, Gilardi MC, Fazio F,",
  "   Ruotsalainen U. Implementation and evaluation of an ordered subsets",
  "   reconstruction algorithm for transmission PET studies using median root",
  "   prior and inter-update median filtering.",
  "   Eur J Nucl Med. 2003; 30(2):222-231.",
  " ",
  "See also: ecatfbp, ecatmrp, ecatnorm, lmhdr, egetstrt, edecay",
  " ",
  "Keywords: ECAT, sinogram, reconstruction, attenuation",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
/*****************************************************************************/
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
) {
  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);
}
/*****************************************************************************/
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  char    blkfile[FILENAME_MAX], trafile[FILENAME_MAX], norfile[FILENAME_MAX],
          atnfile[FILENAME_MAX], imgfile[FILENAME_MAX];
  int     decayCorrection=1;
  int     limit=0; /* 0=no; 1=limits applied */
  int     maxIterNr=45;
  int     skipPriorNr=1;
  float   beta=0.9;
  float   rotate=0.0;
  float   zoom=1.0;
  float   shiftX=0.0, shiftY=0.0;
  int     imgDim=0;
  int     maskDim=3;
  int     osSetNr=1;
  char   *cptr;
  int     ret;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  blkfile[0]=trafile[0]=norfile[0]=atnfile[0]=imgfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* options */
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "DECAY=", 6)==0) {
      cptr+=6;
      if(strncasecmp(cptr, "YES", 1)==0 || strcasecmp(cptr, "ON")==0) {
        decayCorrection=1; continue;
      }
      if(strncasecmp(cptr, "NO", 1)==0 || strcasecmp(cptr, "OFF")==0) {
        decayCorrection=0; continue;
      }
    } else if(strncasecmp(cptr, "LIMIT=", 6)==0) {
      cptr+=6;
      if(strncasecmp(cptr, "YES", 1)==0 || strcasecmp(cptr, "ON")==0) {
        limit=1; continue;
      }
      if(strncasecmp(cptr, "NO", 1)==0 || strcasecmp(cptr, "OFF")==0) {
        limit=0; continue;
      }
    } else if(strncasecmp(cptr, "ITER=", 5)==0) {
      maxIterNr=atoi(cptr+5); if(maxIterNr>1) continue;
    } else if(strncasecmp(cptr, "SKIP=", 5)==0) {
      skipPriorNr=atoi(cptr+5); if(skipPriorNr>0) continue;
    } else if(strncasecmp(cptr, "OS=", 3)==0) {
      osSetNr=atoi(cptr+3); 
      if(osSetNr==0 || osSetNr==1 || osSetNr==2) continue;
      if(osSetNr==4 || osSetNr==8 || osSetNr==16) continue;
      if(osSetNr==32 || osSetNr==64 || osSetNr==128) continue;
    } else if(strncasecmp(cptr, "DIM=", 4)==0) {
      imgDim=atoi(cptr+4); if(imgDim>1 && imgDim<=2048) continue;
    } else if(strncasecmp(cptr, "MASK=", 5)==0) {
      maskDim=atoi(cptr+5); if(maskDim==3 || maskDim==5) continue;
    } else if(strncasecmp(cptr, "ZOOM=", 5)==0) {
      zoom=atof_dpi(cptr+5); if(zoom>0.1 && zoom<100.0) continue;
    } else if(strncasecmp(cptr, "BETA=", 5)==0) {
      beta=atof_dpi(cptr+5); if(beta>0.05 && beta<100.0) continue;
    } else if(strncasecmp(cptr, "ROTATION=", 9)==0) {
      double v;
      ret=atof_with_check(cptr+9, &v);
      if(ret==0 && v>=-180.0 && v<=180.0) {rotate=v; continue;}
    } else if(strncasecmp(cptr, "ROT=", 4)==0) {
      double v;
      ret=atof_with_check(cptr+4, &v);
      if(ret==0 &&  v>-360.0 && v<360.0) {rotate=v; continue;}
    } else if(strncasecmp(cptr, "X=", 2)==0) {
      double v;
      ret=atof_with_check(cptr+2, &v);
      if(ret==0 && v>-50.0 && v<50.0) {shiftX=10.0*v; continue;}
    } else if(strncasecmp(cptr, "Y=", 2)==0) {
      double v;
      ret=atof_with_check(cptr+2, &v);
      if(ret==0 && v>-50.0 && v<50.0) {shiftY=10.0*v; continue;}
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) {strlcpy(blkfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(trafile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(norfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(atnfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(imgfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: too many arguments.\n");
    return(1);
  }
 
  /* Is something missing? */
  if(!atnfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("blkfile := %s\n", blkfile);
    printf("trafile := %s\n", trafile);
    printf("norfile := %s\n", norfile);
    printf("atnfile := %s\n", atnfile);
    if(imgfile[0]) printf("imgfile := %s\n", imgfile);
    printf("decayCorrection := %d\n", decayCorrection);
    printf("limit := %d\n", limit);
    printf("beta := %g\n", beta);
    printf("zoom := %g\n", zoom);
    if(imgDim>0) printf("dimension := %d\n", imgDim);
    printf("mask := %d\n", maskDim);
    printf("rotation := %g\n", rotate);
    printf("x_shift_mm := %g\n", shiftX);
    printf("y_shift_mm := %g\n", shiftY);
    printf("maxIterNr := %d\n", maxIterNr);
    printf("skipPriorNr := %d\n", skipPriorNr);
    printf("osSetNr := %d\n", osSetNr);
  }
 
  /* Check that input and output files do not have the same filename */
  if(!strcasecmp(atnfile, blkfile) || !strcasecmp(atnfile, trafile)) {
    fprintf(stderr, "Error: original data must not be overwritten.\n");
    return(1);
  }
  if(!strcasecmp(imgfile, blkfile) || !strcasecmp(imgfile, trafile)) {
    fprintf(stderr, "Error: original data must not be overwritten.\n");
    return(1);
  }
  /* Output filenames must have correct extension */
  if(strcasestr(atnfile, ".atn")==NULL) {
    fprintf(stderr, "Error: attenuation file must have extension .atn\n");
    return(1);
  }
  if(imgfile[0] && strcasestr(imgfile, ".img")==NULL) {
    fprintf(stderr, "Error: transmission image must have extension .img\n");
    return(1);
  }
 
 
  /* Call the function that actually does all the work */
  char buf[128];
  ret=atnMake(
    trafile, blkfile, norfile, atnfile, imgfile, decayCorrection, limit, 1,
    imgDim, zoom, shiftX, shiftY, rotate,
    maxIterNr, skipPriorNr, beta, maskDim, osSetNr,
    buf, verbose
  );
  if(ret!=0) {
    fprintf(stderr, "Error: %s.\n", buf);
    if(verbose>3) printf("ret := %d\n", ret);
    return(1);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/