ecatnorm.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/stat.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcimgio.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Normalization and optionally attenuation correction for CTI ECAT 931",
  "transmission or emission sinogram.",
  " ",
  "Usage: @P [Options] scnfile nrmfile [atnfile] scnfile2",
  " ",
  "Options:",
  " -rm",
  "     Normalization and attenuation correction are removed from a previously",
  "     corrected sinogram.",
  " -dtc=<y|N>",
  "     Dead-time correction is done (y), or not done (n, default);",
  "     if normalization is to be removed, then this option decides whether",
  "     dead-time correction is removed (y, default), or kept (n).",
  "     Reconstruction programs ecatfbp and ecatmrp correct the data for",
  "     dead-time, therefore it is not done by default here.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1: Normalization correction for transmission sinogram",
  "  @P a2345tr1.scn a2345sta.nrm a2345tr1_normalized.scn",
  " ",
  "Example 2: Normalization and attenuation correction for emission sinogram",
  "  @P a2345dy1.scn a2345sta.nrm a2345tr1.atn a2345dy1_corrected.scn",
  " ",
  "See also: ecatfbp, ecatmrp, atnmake, egetstrt, edecay, ecalibr, eframe",
  " ",
  "Keywords: ECAT, sinogram, reconstruction, normalization, 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 main(int argc, char **argv)
{
  int        ai, help=0, version=0, verbose=1;
  char       scnfile[FILENAME_MAX], nrmfile[FILENAME_MAX], 
             atnfile[FILENAME_MAX], scnfile2[FILENAME_MAX];
  int        dtc=0; // by default no, because done when scan is read to IMG
  int        normalization=1; // 0=remove corrections; 1=correct
  char      *cptr;
  int        ret;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  scnfile[0]=nrmfile[0]=atnfile[0]=scnfile2[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, "DTC=", 4)==0) {
      cptr+=4;
      if(strncasecmp(cptr, "YES", 1)==0 || strcasecmp(cptr, "ON")==0) {
        dtc=1; continue;
      }
      if(strncasecmp(cptr, "NO", 1)==0 || strcasecmp(cptr, "OFF")==0) {
        dtc=0; continue;
      }
    } else if(strcasecmp(cptr, "RM")==0) {
      normalization=0; 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(scnfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(nrmfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc-1) {strlcpy(atnfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(scnfile2, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: too many arguments.\n");
    return(1);
  }
 
  /* Is something missing? */
  if(!scnfile2[0]) {
    fprintf(stderr, "Error: missing command-line argument; use --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("scnfile := %s\n", scnfile);
    printf("nrmfile := %s\n", nrmfile);
    if(atnfile[0]) printf("atnfile := %s\n", atnfile);
    printf("scnfile2 := %s\n", scnfile2);
    printf("normalization := %d\n", normalization);
    printf("dtc := %d\n", dtc);
  }
 
  /* Check that input and output files do not have the same filename */
  if(strcasecmp(scnfile, scnfile2)==0) {
    fprintf(stderr, "Error: original sinogram must not be overwritten.\n");
    return(1);
  }
 
 
  /*
   *  Open CTI sinogram file
   */
  if(verbose==1) printf("reading %s\n", scnfile);
  FILE *fpscn=NULL;
  if(verbose>1) printf("opening %s\n", scnfile);
  if((fpscn=fopen(scnfile, "rb")) == NULL) {
    fprintf(stderr, "Error: cannot open file %s\n", scnfile); 
    return(2);
  }
 
  /*
   *  Read sinogram main header
   */
  if(verbose>1) printf("reading main header\n");
  ECAT63_mainheader scn_main_header;
  if((ret=ecat63ReadMainheader(fpscn, &scn_main_header))) {
    fprintf(stderr, "Error: cannot read main header.\n");
    fclose(fpscn); return(2);
  }
  if(verbose>10) ecat63PrintMainheader(&scn_main_header, stdout);
 
  /* Check the file_type */
  if(scn_main_header.file_type!=RAW_DATA) {
    fprintf(stderr, "Error: only ECAT 931 sinogram can be normalized.\n");
    fclose(fpscn); return(2);
  }
 
  /*
   *  Read sinogram matrix list
   */
  if(verbose>1) printf("reading matrix list\n");
  static MATRIXLIST scn_mlist; ecat63InitMatlist(&scn_mlist);
  ret=ecat63ReadMatlist(fpscn, &scn_mlist, verbose-1);
  if(ret) {
    fprintf(stderr, "Error: cannot read matrix list.\n");
    fclose(fpscn); return(2);
  }
  if(scn_mlist.matrixNr<=0) {
    fprintf(stderr, "Error: matrix list is empty.\n");
    fclose(fpscn); return(2);
  }
  if(verbose>1) printf("sinogram_matrixNr := %d\n", scn_mlist.matrixNr);
  if(ecat63CheckMatlist(&scn_mlist)) {
    fprintf(stderr, "Error: check the matrix list.\n");
    fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
    return(2);
  }
 
 
 
 
  /*
   *  Open CTI normalization file
   */
  if(verbose==1) printf("reading %s\n", nrmfile);
  FILE *fpnrm=NULL;
  if(verbose>1) printf("opening %s\n", nrmfile);
  if((fpnrm=fopen(nrmfile, "rb")) == NULL) {
    fprintf(stderr, "Error: cannot open file %s\n", nrmfile); 
    fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
    return(3);
  }
 
  /*
   *  Read normalization main header
   */
  if(verbose>1) printf("reading main header\n");
  ECAT63_mainheader nrm_main_header;
  if((ret=ecat63ReadMainheader(fpnrm, &nrm_main_header))) {
    fprintf(stderr, "Error: cannot read main header.\n");
    fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
    fclose(fpnrm); 
    return(3);
  }
  if(verbose>10) ecat63PrintMainheader(&nrm_main_header, stdout);
 
  /* Check the file_type */
  if(nrm_main_header.file_type!=RAW_DATA) { 
    // files are in sinogram format, not NORM_DATA 
    fprintf(stderr, "Error: check the normalization file.\n");
    fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
    fclose(fpnrm); 
    return(3);
  }
 
  /*
   *  Read normalization matrix list
   */
  if(verbose>1) printf("reading matrix list\n");
  static MATRIXLIST nrm_mlist; ecat63InitMatlist(&nrm_mlist);
  ret=ecat63ReadMatlist(fpnrm, &nrm_mlist, verbose-1);
  if(ret) {
    fprintf(stderr, "Error: cannot read matrix list.\n");
    fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
    fclose(fpnrm); 
    return(3);
  }
  if(nrm_mlist.matrixNr<=0) {
    fprintf(stderr, "Error: matrix list is empty.\n");
    fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
    fclose(fpnrm); 
    return(3);
  }
  if(verbose>1) printf("normalization_matrixNr := %d\n", nrm_mlist.matrixNr);
  if(ecat63CheckMatlist(&nrm_mlist)) {
    fprintf(stderr, "Error: check the matrix list.\n");
    fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
    fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
    return(3);
  }
 
 
  /*
   *  Check that sinogram and normalization file have the same nr of planes
   */
  int scn_plane_nr, scn_frame_nr, nrm_plane_nr, nrm_frame_nr;
  ret=ecat63GetPlaneAndFrameNr(&scn_mlist, &scn_main_header, 
                               &scn_plane_nr, &scn_frame_nr);
  if(ret==0) ret=ecat63GetPlaneAndFrameNr(&nrm_mlist, &nrm_main_header, 
                                          &nrm_plane_nr, &nrm_frame_nr);
  if(ret) {
    fprintf(stderr, "Error: invalid matrix list.\n");
    fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
    fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
    return(4);
  }
  if(verbose>1) {
    printf("sinogram_plane_nr := %d\n", scn_plane_nr);
    printf("normalization_plane_nr := %d\n", nrm_plane_nr);
    printf("sinogram_frame_nr := %d\n", scn_frame_nr);
  }
  if(scn_plane_nr!=nrm_plane_nr) {
    fprintf(stderr, "Error: plane numbers do not match.\n");
    fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
    fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
    return(4);
  }
  if(nrm_frame_nr>1) {
    fprintf(stderr, "Error: invalid normalization file.\n");
    fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
    fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
    return(4);
  }
 
 
  /*
   *  Optional attenuation data
   */
  FILE *fpatn=NULL;
  ECAT63_mainheader atn_main_header;
  static MATRIXLIST atn_mlist; ecat63InitMatlist(&atn_mlist);
 
  if(atnfile[0]) {
 
    /*
     *  Open CTI attenuation file
     */
    if(verbose==1) printf("reading %s\n", atnfile);
    if(verbose>1) printf("opening %s\n", atnfile);
    if((fpatn=fopen(atnfile, "rb")) == NULL) {
      fprintf(stderr, "Error: cannot open file %s\n", atnfile); 
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      return(5);
    }
 
    /*
     *  Read attenuation main header
     */
    if(verbose>1) printf("reading main header\n");
    if((ret=ecat63ReadMainheader(fpatn, &atn_main_header))) {
      fprintf(stderr, "Error: cannot read main header.\n");
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      if(fpatn!=NULL) fclose(fpatn); 
      return(5);
    }
    if(verbose>10) ecat63PrintMainheader(&atn_main_header, stdout);
 
    /* Check the file_type */
    if(atn_main_header.file_type!=ATTN_DATA) { 
      fprintf(stderr, "Error: check the attenuation file.\n");
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      if(fpatn!=NULL) fclose(fpatn); 
      return(5);
    }
 
    /*
     *  Read attenuation matrix list
     */
    if(verbose>1) printf("reading matrix list\n");
    ret=ecat63ReadMatlist(fpatn, &atn_mlist, verbose-1);
    if(ret) {
      fprintf(stderr, "Error: cannot read matrix list.\n");
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      if(fpatn!=NULL) fclose(fpatn); 
      return(5);
    }
    if(atn_mlist.matrixNr<=0) {
      fprintf(stderr, "Error: matrix list is empty.\n");
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);
      if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
      return(5);
    }
    if(verbose>1) printf("attenuation_matrixNr := %d\n", atn_mlist.matrixNr);
    if(ecat63CheckMatlist(&atn_mlist)) {
      fprintf(stderr, "Error: check the matrix list.\n");
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
      return(5);
    }
 
 
    /*
     *  Check that sinogram and attenuation file have the same nr of planes
     */
    int atn_plane_nr, atn_frame_nr;
    ret=ecat63GetPlaneAndFrameNr(&atn_mlist, &atn_main_header, 
                                 &atn_plane_nr, &atn_frame_nr);
    if(ret) {
      fprintf(stderr, "Error: invalid matrix list.\n");
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
      return(6);
    }
    if(verbose>1) {
      printf("attenuation_plane_nr := %d\n", atn_plane_nr);
    }
    if(scn_plane_nr!=atn_plane_nr) {
      fprintf(stderr, "Error: plane numbers do not match.\n");
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
      return(6);
    }
    if(atn_frame_nr>1) {
      fprintf(stderr, "Error: invalid attenuation file.\n");
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
      return(6);
    }
  } // attenuation file matrix list ready
 
 
 
 
  /*
   *  Create new sinogram data file and write the main header
   */
  if(verbose>0) printf("creating %s\n", scnfile2);
  FILE *fpout=NULL;
  fpout=ecat63Create(scnfile2, &scn_main_header);
  if(fpout==NULL) {
    fprintf(stderr, "Error: cannot write file %s\n", scnfile2); 
    fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
    fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
    if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
    return(11);
  }
 
 
  /*
   *  Process one matrix at a time
   */
  if(verbose>0) {
    char tmp[64];
    if(!atnfile[0]) strcpy(tmp, "");
    else strcpy(tmp, " and attenuation correction");
    if(normalization) printf("normalization%s of sinogram...\n", tmp);
    else printf("removing normalization%s from sinogram...\n", tmp);
    fflush(stdout);
  }
 
  ECAT63_scanheader scan_header, norm_header;
  ECAT63_attnheader attn_header;
  Matval matval;
  int nr=0, pxlNr=0, matnum=0, ni=0;
  float *scnmat=NULL, *nrmmat=NULL;
 
  for(int mi=nr=0; mi<scn_mlist.matrixNr; mi++) {
 
    /* Get plane and frame nr */
    mat_numdoc(scn_mlist.matdir[mi].matnum, &matval);
    /* Read sinogram subheader and scaled data */
    ret=ecat63ReadScanMatrix(fpscn, scn_mlist.matdir[mi].strtblk, 
          scn_mlist.matdir[mi].endblk, &scan_header, &scnmat);
    if(ret) {
      fprintf(stderr, "Error: cannot read sinogram matrix %u.\n",
         scn_mlist.matdir[mi].matnum);
      fflush(stderr);
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);
      if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
      fclose(fpout); remove(scnfile2);
      return(5);
    }
    if(verbose>8) {
      printf("Matrix: plane %d frame %d gate %d bed %d\n",
             matval.plane, matval.frame, matval.gate, matval.bed);
      fflush(stdout);
    }
    if(verbose>20) {
      ecat63PrintScanheader(&scan_header, stdout); fflush(stdout);
    } else if(verbose>1 && mi==0) {
      printf("Sinogram dimensions := %d x %d\n",
             scan_header.dimension_1, scan_header.dimension_2);
      fflush(stdout);
    }
    pxlNr=scan_header.dimension_1*scan_header.dimension_2;
    /* If dead-time correction was requested, then do it */
    if(dtc!=0 && scan_header.loss_correction_fctr>0.0) {
      if(normalization) {
        for(int i=0; i<pxlNr; i++) 
          scnmat[i]*=scan_header.loss_correction_fctr;
      } else {
        for(int i=0; i<pxlNr; i++) 
          scnmat[i]/=scan_header.loss_correction_fctr;
      }
    }
 
    /* Search corresponding matrix from the normalization file */
    matnum=mat_numcod(1, matval.plane, matval.gate, matval.data, matval.bed);
    for(ni=0; ni<nrm_mlist.matrixNr; ni++)
      if(nrm_mlist.matdir[ni].matnum==matnum)
        break;
    if(ni>=nrm_mlist.matrixNr) {
      fprintf(stderr, "Error: cannot find matching normalization matrix.\n");
      fflush(stderr);
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
      free(scnmat);
      fclose(fpout); remove(scnfile2);
      return(8);
    }
    /* Read normalization subheader and scaled data */
    ret=ecat63ReadScanMatrix(fpnrm, nrm_mlist.matdir[ni].strtblk, 
          nrm_mlist.matdir[ni].endblk, &norm_header, &nrmmat);
    if(ret) {
      fprintf(stderr, "Error: cannot read normalization matrix %u.\n",
         nrm_mlist.matdir[ni].matnum);
      fflush(stderr);
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
      free(scnmat);
      fclose(fpout); remove(scnfile2);
      return(8);
    }
    if(verbose>20) {
      ecat63PrintScanheader(&norm_header, stdout); fflush(stdout);
    } else if(verbose>1 && mi==0) {
      printf("Normalization dimensions := %d x %d\n",
             norm_header.dimension_1, norm_header.dimension_2);
      fflush(stdout);
    }
    if(norm_header.dimension_1!=scan_header.dimension_1
       || norm_header.dimension_2!=scan_header.dimension_2)
    {
      fprintf(stderr, "Error: incompatible matrix dimensions.\n");
      fflush(stderr);
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
      free(scnmat); free(nrmmat);
      fclose(fpout); remove(scnfile2);
      return(8);
    }
 
    /* Search corresponding matrix from the attenuation file */
    if(atnfile[0]) {
      matnum=mat_numcod(1, matval.plane, matval.gate, matval.data, matval.bed);
      int ai;
      for(ai=0; ai<atn_mlist.matrixNr; ai++)
        if(atn_mlist.matdir[ai].matnum==matnum)
          break;
      if(ai>=atn_mlist.matrixNr) {
        fprintf(stderr, "Error: cannot find matching attenuation matrix.\n");
        fflush(stderr);
        fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
        fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
        if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
        free(scnmat); free(nrmmat);
        fclose(fpout); remove(scnfile2);
        return(8);
      }
      /* Read attenuation subheader and scaled data */
      float *atnmat=NULL;
      ret=ecat63ReadAttnMatrix(fpatn, atn_mlist.matdir[ai].strtblk, 
            atn_mlist.matdir[ai].endblk, &attn_header, &atnmat);
      if(ret) {
        fprintf(stderr, "Error: cannot read attenuation matrix %u.\n",
                atn_mlist.matdir[ai].matnum);
        fflush(stderr);
        fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
        fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
        if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
        free(scnmat); free(nrmmat);
        fclose(fpout); remove(scnfile2); 
        return(8);
      }
      if(verbose>20) {
        ecat63PrintAttnheader(&attn_header, stdout); fflush(stdout);
      } else if(verbose>1 && mi==0) {
        printf("Attenuation dimensions := %d x %d\n",
               attn_header.dimension_1, attn_header.dimension_2);
        fflush(stdout);
      }
      /* Sometimes attenuation dimensions are just in wrong order */
      if((attn_header.dimension_1!=scan_header.dimension_1
          || attn_header.dimension_2!=scan_header.dimension_2) &&
         (attn_header.dimension_2!=scan_header.dimension_1
          || attn_header.dimension_1!=scan_header.dimension_2))
      {
        fprintf(stderr, "Error: incompatible matrix dimensions.\n");
        fflush(stderr);
        fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
        fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
        if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
        free(scnmat); free(nrmmat); free(atnmat);
        fclose(fpout); remove(scnfile2); 
        return(8);
      }
 
      /* If attenuation value is < 1.0E-06, then use that as a limit */
      for(int i=0; i<pxlNr; i++) if(atnmat[i]<1.0E-06) atnmat[i]=1.0E-06; 
 
      /* Multiply or divide sinogram matrix with attenuation data */
      if(normalization) {
        for(int i=0; i<pxlNr; i++) 
          scnmat[i]*=atnmat[i];
      } else {
        for(int i=0; i<pxlNr; i++) 
          scnmat[i]/=atnmat[i]; // attenuation data should be >= 1 
      }
 
      /* Attenuation data is no more needed for this matrix */
      free(atnmat); 
    }
 
    /* If normalization value is < 1.0E-06, then use that as a limit */
    for(int i=0; i<pxlNr; i++) if(nrmmat[i]<1.0E-06) nrmmat[i]=1.0E-06; 
 
    /* Multiply or divide sinogram data with normalization data */
    if(normalization) {
      for(int i=0; i<pxlNr; i++) 
        scnmat[i]*=nrmmat[i];
    } else {
      for(int i=0; i<pxlNr; i++)
        scnmat[i]/=nrmmat[i]; 
    }
    free(nrmmat);
 
 
    /* Write the corrected sinogram matrix */
    ret=ecat63WriteScanMatrix(fpout, scn_mlist.matdir[mi].matnum, &scan_header,
                              scnmat);
    if(ret) {
      fprintf(stderr, "Error: cannot write matrix.\n");
      fflush(stderr);
      if(verbose>1) printf("ret := %d\n", ret);
      fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
      fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
      if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
      free(scnmat); free(nrmmat);
      fclose(fpout); remove(scnfile2); 
      return(12);
    }
    free(scnmat);
    nr++; // counter
  } /* next matrix */
  if(verbose>0) {printf("%d matrices processed.\n", nr); fflush(stdout);}
 
  fclose(fpscn); ecat63EmptyMatlist(&scn_mlist); 
  fclose(fpnrm); ecat63EmptyMatlist(&nrm_mlist);  
  if(fpatn!=NULL) {fclose(fpatn); ecat63EmptyMatlist(&atn_mlist);} 
  fclose(fpout); 
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/