dcmframe.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpctac.h"
#include "tpcdcm.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "List the time frame information of a PET image in DICOM format.",
  "The frame start times and lengths can be saved in SIF.",
  " ",
  "NOT for production use!",
  " ",
  "Usage: @P [-Options] dicomfile [SIF]",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: dcmlhdr, dcmmlist, eframe, tacframe",
  " ",
  "Keywords: image, DICOM, time, SIF",
  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  dcmfile[FILENAME_MAX], siffile[FILENAME_MAX];
  int   ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dcmfile[0]=siffile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    //char *cptr; cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-3;
  
  /* 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(dcmfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(siffile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(!dcmfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dcmfile := %s\n", dcmfile);
    if(siffile[0]) printf("siffile := %s\n", siffile);
  }
 
  /*
   *  Get list of DICOM files belonging to the same image
   */
  IFT fl; iftInit(&fl);
  ret=dcmFileList(dcmfile, &fl, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    return(2);
  }
  if(verbose>1) {printf("fileNr := %d\n", fl.keyNr); fflush(stdout);}
  if(verbose>6) iftWrite(&fl, stdout, NULL);
 
  /* 
   *  Read the first file in list
   */
  DCMFILE dcm; dcmfileInit(&dcm);
  ret=dcmFileRead(fl.item[0].value, &dcm, 1, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    iftFree(&fl); dcmfileFree(&dcm); return(3);
  }
 
  /* Read some information common to the whole image */
  DCMTAG tag;
  DCMITEM *iptr;
 
  /* Get Decay Correction */
  short unsigned int decayCorrection=0; // 0=NONE, 1=START, 2=ADMIN
  tag.group=0x0054; tag.element=0x1102;
  iptr=dcmFindTag(dcm.item, 0, &tag, 0);
  if(iptr!=NULL) {
    char *buf=dcmValueString(iptr); 
    if(verbose>2) printf("Decay Correction := %s\n", buf);
    if(strcasestr(buf, "NONE")!=NULL) decayCorrection=0;
    else if(strcasestr(buf, "START")!=NULL) decayCorrection=1;
    else if(strcasestr(buf, "ADMIN")!=NULL) decayCorrection=2;
    free(buf);
  } else {
    tag.group=0x0018; tag.element=0x9758;
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr!=NULL) {
      char *buf=dcmValueString(iptr); 
      if(verbose>2) printf("Decay Correction := %s\n", buf);
      if(strcasestr(buf, "YES")!=NULL) decayCorrection=1;
      free(buf);
    }
  }
  if(verbose>1) printf("decayCorrection := %d\n", decayCorrection);
 
  /* If decay is corrected to tracer administration time, the we must read it */
  char injDateTime[32]; injDateTime[0]=(char)0;
  if(decayCorrection==2) {
    /* Get the Tracer Start DateTime */
    tag.group=0x0018; tag.element=0x1078;
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr!=NULL) {
      if(dcmDT2intl(iptr->rd, injDateTime)==NULL) iptr=NULL;
      else if(verbose>3) printf("0018,1078 -> %s\n", injDateTime);
    }
    /* If not successful, the try other tags */
    if(iptr==NULL) {
      char s1[16], s2[16];
      /* Get the Tracer Start Time */
      tag.group=0x0018; tag.element=0x1072;
      iptr=dcmFindTag(dcm.item, 0, &tag, 0);
      if(iptr!=NULL) {
        /* Get the Series Date */
        tag.group=0x0008; tag.element=0x0021;
        DCMITEM *jptr=dcmFindTag(dcm.item, 0, &tag, 0);
        if(jptr!=NULL) {
          if(dcmDA2intl(jptr->rd, s1)!=NULL && dcmTM2intl(iptr->rd, s2)!=NULL)
            sprintf(injDateTime, "%s %s", s1, s2);
        }
      }
    }
    /* Check that we got injection date and time, since we will need it */
    if(!injDateTime[0]) {
      fprintf(stderr, "Error: missing Tracer Administration Time.\n");
      iftFree(&fl); dcmfileFree(&dcm);
      return(4);
    }
    if(verbose>1) printf("injection_time := %s\n", injDateTime); 
  }
 
  /* Get the Series DateTime
     https://dicom.nema.org/MEDICAL/dicom/current/output/chtml/part03/sect_C.8.9.html#figure_C.8.9.1.1.11-1a
   */
  char seriesDateTime[32]; seriesDateTime[0]=(char)0;
  tag.group=0x0008; tag.element=0x0021;
  iptr=dcmFindTag(dcm.item, 0, &tag, 0);
  if(iptr!=NULL) {
    tag.group=0x0008; tag.element=0x0031;
    DCMITEM *jptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(jptr!=NULL) {
      char s1[16], s2[16];
      if(dcmDA2intl(iptr->rd, s1)!=NULL && dcmTM2intl(jptr->rd, s2)!=NULL)
        sprintf(seriesDateTime, "%s %s", s1, s2);
    }
  }
  if(!seriesDateTime[0]) {
    fprintf(stderr, "Error: missing Series Date and Time.\n");
    iftFree(&fl); dcmfileFree(&dcm); return(5);
  } else {
    if(verbose>2) printf("seriesDateTime := %s\n", seriesDateTime); 
  }
 
 
  /* Get and check the series type */
  tag.group=0x0054; tag.element=0x1000;
  iptr=dcmFindTag(dcm.item, 0, &tag, 0);
  if(iptr==NULL) {
    if(verbose>1) fprintf(stderr, "Warning: cannot find series type.\n");
  } else {
    char *buf=dcmValueString(iptr); 
    if(verbose>1) printf("seriesType := %s\n", buf); 
    /* Check that data is either dynamic or static */
    if(strncasecmp(buf, "DYNAMIC", 7) && strncasecmp(buf, "STATIC", 6)) {
      fprintf(stderr, "Error: data series type '%s' not supported.\n", buf);
      free(buf);
      iftFree(&fl); dcmfileFree(&dcm); return(5);
    }
    free(buf);
  }
 
 
  /* Get and check the Image Type */
  int dynamic=0; // 0=no, 1=yes
  tag.group=0x0008; tag.element=0x0008;
  iptr=dcmFindTag(dcm.item, 0, &tag, 0);
  if(iptr==NULL) {
    fprintf(stderr, "Error: cannot find image type.\n");
    iftFree(&fl); dcmfileFree(&dcm); return(5);
  }
  /* Check if dynamic; see DICOM 3.3 C.7.6.16.2.2.6
     https://dicom.nema.org/MEDICAL/dicom/current/output/chtml/part03/sect_C.7.6.16.2.html
   */
  {
    char *buf=dcmValueString(iptr); 
    if(verbose>1) printf("imageType := %s\n", buf); 
    if(strcasestr(buf, "DYNAMIC")!=NULL) dynamic=1; else dynamic=0;
    free(buf);
  }
  if(verbose>1) {printf("dynamic := %d\n", dynamic); fflush(stdout);}
 
 
  /* Check whether Per-frame Functional Groups Sequence is set */
  int multiframe=0; // 0=no, 1=yes
  tag.group=0x5200; tag.element=0x9230;
  iptr=dcmFindTag(dcm.item, 0, &tag, 0);
  if(iptr!=NULL) {
    if(verbose>1) fprintf(stderr, "Notice: Per-frame Functional Groups Sequence is available.\n");
    multiframe=1;
  }
 
 
  /* Get the number of frames (time slices) */
  if(verbose>1) {printf("reading frame number\n"); fflush(stdout);}
  unsigned short int frameNr=0;
  tag.group=0x0054; tag.element=0x0101;
  iptr=dcmFindTag(dcm.item, 0, &tag, 0);
  if(iptr!=NULL) { 
    frameNr=(unsigned short int)dcmitemGetInt(iptr);
  } else {
    /* Get it from 'Temporal position index' */
    tag.group=0x0020; tag.element=0x9128;
    int a, b;
    if(dcmTagIntRange(dcm.item, &tag, &a, &b, verbose-100)==0) frameNr=b;
  }
  /* Check that frame number > 0 if dynamic */
  if(frameNr==0) {
    if(dynamic!=0) { // error, if dynamic
      fprintf(stderr, "Error: cannot find frame number.\n");
      iftFree(&fl); dcmfileFree(&dcm); return(4);
    } else { // assume 1, if static or whole-body
      frameNr=1;
    }
  }
  if(verbose>1) printf("frameNr := %u\n", frameNr); 
 
 
 
  /* Get the number of slices (needed to understand image index) */
  if(verbose>1) {printf("reading plane number\n"); fflush(stdout);}
  unsigned short int sliceNr=0;
  tag.group=0x0054; tag.element=0x0081;
  iptr=dcmFindTag(dcm.item, 0, &tag, 0);
  if(iptr!=NULL) {
    sliceNr=(unsigned short int)dcmitemGetInt(iptr);
  } else {
    /* Get it from 'In Stack Position Number' */
    tag.group=0x0020; tag.element=0x9057;
    int a, b;
    if(dcmTagIntRange(dcm.item, &tag, &a, &b, verbose-100)==0) sliceNr=b;
    else {
      /* Get it from 'Dimension Index Values' */
      tag.group=0x0020; tag.element=0x9157;
      if(dcmTagIntRange(dcm.item, &tag, &a, &b, verbose-1)==0) sliceNr=b;
    }
  }
  if(sliceNr<=0) {
    fprintf(stderr, "Error: cannot find slice number.\n");
    iftFree(&fl); dcmfileFree(&dcm); return(4);
  }
  if(verbose>2) printf("sliceNr := %u\n", sliceNr); 
 
 
  /* Get the isotope */
  if(verbose>1) {printf("reading isotope\n"); fflush(stdout);}
  int isotope_code=ISOTOPE_UNKNOWN;
  tag.group=0x0018; tag.element=0x1075;
  iptr=dcmFindTag(dcm.item, 0, &tag, 0);
  if(iptr!=NULL) {
    double halflife=dcmitemGetReal(iptr);
    if(verbose>2) printf("isotope_halflife := %g\n", halflife);
    isotope_code=isotopeIdentifyHalflife(halflife/60.0);
  }
  if(isotope_code==ISOTOPE_UNKNOWN && verbose>0) {
    fprintf(stderr, "Warning: isotope not identified.\n");
  }
 
 
  /* Get accession (study) number */
  if(verbose>1) {printf("reading accession number\n"); fflush(stdout);}
  char studyNr[17]; studyNr[0]=(char)0;
  tag.group=0x0008; tag.element=0x0050;
  iptr=dcmFindTag(dcm.item, 0, &tag, 0);
  if(iptr!=NULL) {
    strlcpy(studyNr, iptr->rd, 17);
    strClean(studyNr);
    strReplaceChar(studyNr, ' ', '_');
  }
  if(strlen(studyNr)==0) {
    if(verbose>1) {printf("reading series description\n"); fflush(stdout);}
    tag.group=0x0008; tag.element=0x0103e;
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr!=NULL) {
      strlcpy(studyNr, iptr->rd, 17);
      strClean(studyNr);
      char *cptr=strrchr(studyNr, ' '); if(cptr!=NULL) *cptr=(char)0;
      cptr=strrchr(studyNr, '/'); if(cptr!=NULL) *cptr=(char)0;
    }
  }
  if(verbose>1) printf("studyNr := %s\n", studyNr);
 
 
 
  /*
   *  Allocate memory for frame time data
   */
  if(verbose>1) {printf("allocating memory for frame times\n"); fflush(stdout);}
  TAC tac; tacInit(&tac);
  ret=tacAllocate(&tac, frameNr, 3);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(ret));
    iftFree(&fl); dcmfileFree(&dcm);
    return(5);
  }
  tac.format=TAC_FORMAT_SIF;
  tac.sampleNr=frameNr;
  tac.isframe=1;
  tac.tunit=UNIT_SEC;
  tac.cunit=UNIT_COUNTS;
  tac.tacNr=2;
  /* Set SIF compatible titles into TAC struct */
  strcpy(tac.c[0].name, "Prompts");
  strcpy(tac.c[1].name, "Randoms");
  strcpy(tac.c[2].name, "Trues"); // computed, not saved
  /* Isotope */
  tacSetHeaderIsotope(&tac.h, isotopeName(isotope_code));
  /* Studynumber */
  if(studyNr[0]) tacSetHeaderStudynr(&tac.h, studyNr);
  /* Set frame times to NaN */
  for(int fi=0; fi<frameNr; fi++) tac.x1[fi]=tac.x2[fi]=tac.x[fi]=nan("");
 
 
  /*
   *  Get frame information from one file at a time
   */
  if(verbose>1) {printf("reading frame information\n"); fflush(stdout);}
  char startDateTime[32]; startDateTime[0]=(char)0;
  int fi=0, beds=0;
 
  if(multiframe==0) {
    /*
     *  Get frame information from one file at a time
     */
    if(verbose>1) {printf("frames is separate files\n"); fflush(stdout);}
    while(1) {
 
      /* Get the Image Index; DICOM PS3.3 C.8.9.4.1.9 */
      unsigned short int imageIndex=0;
      tag.group=0x0054; tag.element=0x1330;
      iptr=dcmFindTag(dcm.item, 0, &tag, 0);
      if(iptr==NULL) {
        fprintf(stderr, "Error: cannot find Image Index.\n");
        tacFree(&tac); iftFree(&fl); dcmfileFree(&dcm);
        return(6);
      }
      imageIndex=(unsigned short int)dcmitemGetInt(iptr);
      /* Calculate the current frame and slice (index starting from 0) */
      unsigned short int frameIndex, sliceIndex; 
      {
        div_t meh=div(imageIndex-1, sliceNr);
        frameIndex=meh.quot; 
        sliceIndex=meh.rem; 
      }
      if(verbose>3) {
        printf("imageIndex := %u\n", imageIndex);
        printf("frameIndex := %u\n", frameIndex);
        printf("sliceIndex := %u\n", sliceIndex);
      }
 
      /* If dynamic image, then read certain tags, just for fun */
      if(dynamic && verbose>5) {
        tag.group=0x0020; tag.element=0x9128;
        iptr=dcmFindTag(dcm.item, 0, &tag, 0);
        if(iptr!=0) {
          char *buf=dcmValueString(iptr); 
          /*if(verbose>6)*/ printf("Temporal Position Index := %s\n", buf); 
          free(buf);
        }
        tag.group=0x0020; tag.element=0x9056;
        iptr=dcmFindTag(dcm.item, 0, &tag, 0);
        if(iptr!=0) {
          char *buf=dcmValueString(iptr); 
          /*if(verbose>6)*/ printf("Stack ID := %s\n", buf); 
          free(buf);
        }
        tag.group=0x0020; tag.element=0x9057;
        iptr=dcmFindTag(dcm.item, 0, &tag, 0);
        if(iptr!=0) {
          char *buf=dcmValueString(iptr); 
          /*if(verbose>6)*/ printf("In-Stack Position Number := %s\n", buf); 
          free(buf);
        }
      }
 
      /* Get Frame Reference Time, currently just for fun */
      tag.group=0x0054; tag.element=0x1300;
      iptr=dcmFindTag(dcm.item, 0, &tag, 0);
      if(iptr!=NULL) {
        char *buf=dcmValueString(iptr); 
        if(verbose>5) printf("Frame Reference Time := %s\n", buf); 
        free(buf);
      }
 
 
      /* Get acquisition date and time */
      char acqDateTime[32]; acqDateTime[0]=(char)0;
      /* Try first Acquisition Datetime */
      tag.group=0x0008; tag.element=0x002A;
      iptr=dcmFindTag(dcm.item, 0, &tag, 0);
      if(iptr!=NULL) dcmDT2intl(iptr->rd, acqDateTime);
      /* If that did not work, try Acquisition Date and Time */
      if(!acqDateTime[0]) {
        tag.group=0x0008; tag.element=0x0022;
        iptr=dcmFindTag(dcm.item, 0, &tag, 0);
        if(iptr!=NULL) {
          tag.group=0x0008; tag.element=0x0032;
          DCMITEM *jptr=dcmFindTag(dcm.item, 0, &tag, 0);
          if(jptr!=NULL) {
            char s1[16], s2[16];
            if(dcmDA2intl(iptr->rd, s1)!=NULL && dcmTM2intl(jptr->rd, s2)!=NULL)
              sprintf(acqDateTime, "%s %s", s1, s2);
          }
        }
      }
      if(!acqDateTime[0]) {
        fprintf(stderr, "Error: missing Acquisition Date and Time.\n");
        tacFree(&tac); iftFree(&fl); dcmfileFree(&dcm);
        return(6);
      }
      if(verbose>4) printf("acqDateTime := %s\n", acqDateTime); 
      /* If this is the first frame, then set the start time of the whole image */
      if(frameIndex==0) strcpy(startDateTime, acqDateTime);
 
      /* Calculate the Frame start time */
      double t1;
      if(decayCorrection==2) {
        /* Injection time is the reference time */
        t1=strDateTimeDifference(acqDateTime, injDateTime);
      } else {
        /* Use series time, for now */
        t1=strDateTimeDifference(acqDateTime, seriesDateTime);
      }
      /* Set frame start time */
      if(frameIndex<frameNr) {
        if(isnan(tac.x1[frameIndex])) {
          /* Save frame start time, since there is no previous value for this frame */
          tac.x1[frameIndex]=t1;
          if(verbose>3) printf("t1[%d]=%g at sliceIndex := %u\n", frameIndex, t1, sliceIndex);
        } else {
          /* Check whether previous frame start time is the same as current */
          if(!doubleMatch(tac.x1[frameIndex], t1, 0.001)) {
            beds=1; // not the same, this probably is whole-body study with more than one bed pos
            if(verbose>3) printf("t1=%g at sliceIndex := %u\n", t1, sliceIndex);
          }
        }
      }
 
      /* Get Actual Frame Duration */
      double fdur=0.0;
      tag.group=0x0018; tag.element=0x1242;
      iptr=dcmFindTag(dcm.item, 0, &tag, 0);
      if(iptr!=NULL) fdur=dcmitemGetReal(iptr);
      if(verbose>3) printf("actualFrameDuration := %g\n", fdur); 
      /* For now, put frame duration in place of frame end time */
      if(frameIndex<frameNr) tac.x2[frameIndex]=0.001*fdur;
 
      /* Try to find prompts and Randoms */
      double prompts=0.0, randoms=0.0;
      /* total prompts */
      tag.group=0x0054; tag.element=0x1310;
      iptr=dcmFindTag(dcm.item, 0, &tag, 0);
      if(iptr!=NULL) prompts=dcmitemGetReal(iptr);
      if(prompts<1.0E-03) {
        tag.group=0x0009; tag.element=0x1071;
        iptr=dcmFindTag(dcm.item, 0, &tag, 0);
        if(iptr!=NULL) prompts=dcmitemGetReal(iptr);
      }
      if(prompts<1.0E-03) {
        tag.group=0x0009; tag.element=0x10A7;
        iptr=dcmFindTag(dcm.item, 0, &tag, 0);
        if(iptr!=NULL) prompts=dcmitemGetReal(iptr);
      }
      tac.c[0].y[frameIndex]=prompts;
      /* randoms, delayed */
      tag.group=0x0054; tag.element=0x1311;
      iptr=dcmFindTag(dcm.item, 0, &tag, 0);
      if(iptr!=NULL) randoms=dcmitemGetReal(iptr);
      if(randoms<1.0E-03) {
        tag.group=0x0009; tag.element=0x1072;
        iptr=dcmFindTag(dcm.item, 0, &tag, 0);
        if(iptr!=NULL) randoms=dcmitemGetReal(iptr);
      }
      if(randoms<1.0E-03) {
        tag.group=0x0009; tag.element=0x1022;
        iptr=dcmFindTag(dcm.item, 0, &tag, 0);
        if(iptr!=NULL) randoms=dcmitemGetReal(iptr);
      }
      tac.c[1].y[frameIndex]=randoms;
 
      /* Load the next file */
      dcmfileFree(&dcm);
      fi++; if(fi==fl.keyNr) break;
      ret=dcmFileRead(fl.item[fi].value, &dcm, 1, &status);
      if(ret!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&tac); iftFree(&fl); dcmfileFree(&dcm); return(7);
      }
 
    } // process next file
    iftFree(&fl);
    dcmfileFree(&dcm);
 
    if(verbose>0 && beds!=0) {
      fprintf(stdout, "Note: this seems to be whole-body study\n");
    }
 
  } else {
    /*
     *  Get frame information from the single multi-frame file
     */
    if(verbose>1) {printf("multi-frame file\n"); fflush(stdout);}
 
    /* Get Acquisition DateTime */
    char acqDateTime[32]; acqDateTime[0]=(char)0;
    tag.group=0x0008; tag.element=0x002A;
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr!=NULL) dcmDT2intl(iptr->rd, acqDateTime);
    if(!acqDateTime[0]) {
      fprintf(stderr, "Error: missing Acquisition Date and Time.\n");
      tacFree(&tac); iftFree(&fl); dcmfileFree(&dcm);
      return(6);
    }
    if(verbose>4) printf("acqDateTime := %s\n", acqDateTime); 
    /* Set the start time of the whole image */
    strcpy(startDateTime, acqDateTime);
 
    /* Get Decay Correction DateTime */
    char dcDateTime[32]; dcDateTime[0]=(char)0;
    tag.group=0x0018; tag.element=0x9701;
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    /* Required, if data is decay corrected */
    if(iptr!=NULL) dcmDT2intl(iptr->rd, dcDateTime);
    if(decayCorrection>0 && !dcDateTime[0]) {
      fprintf(stderr, "Error: missing Decay Correction DateTime.\n");
      tacFree(&tac); iftFree(&fl); dcmfileFree(&dcm);
      return(6);
    }
    if(dcDateTime[0] && verbose>4) printf("dcDateTime := %s\n", dcDateTime);
 
    /* Find Per Frame Functional Groups Sequence */
    tag.group=0x5200; tag.element=0x9230;
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr==NULL) {
      fprintf(stderr, "Error: Per Frame Functional Groups Sequence not found.\n");
      tacFree(&tac); iftFree(&fl); dcmfileFree(&dcm); return(7);
    }
    /* Loop through all Frame Content Sequences under it */
    tag.group=0x0020; tag.element=0x9111;
    DCMTAG tagstart; tagstart.group=0x0018; tagstart.element=0x9074;
    DCMTAG tagdur; tagdur.group=0x0018; tagdur.element=0x9220;
    DCMTAG tagfr; tagfr.group=0x0020; tagfr.element=0x9128;
    DCMITEM *fptr=iptr->child_item;
    while(fptr!=NULL) {
      /* Find Frame Content Sequence */
      iptr=dcmFindDownTag(fptr, 0, &tag, 0); if(iptr==NULL) break;
      if(verbose>10) printf(" found Frame Content Sequence\n");
      DCMITEM *jptr;
      /* Find Temporal Position Index */
      jptr=dcmFindDownTag(iptr->child_item, 0, &tagfr, 0);
      if(jptr==NULL) {fptr=fptr->next_item; continue;}
      if(verbose>20) dcmitemPrint(jptr);
      unsigned short int frameIndex=(unsigned short int)(dcmitemGetInt(jptr)-1);
      /* Find Frame Acquisition DateTime */
      jptr=dcmFindDownTag(iptr->child_item, 0, &tagstart, 0);
      char facqDateTime[32]; facqDateTime[0]=(char)0;
      if(jptr!=NULL) dcmDT2intl(jptr->rd, facqDateTime);
      if(!facqDateTime[0]) {
        fprintf(stderr, "Error: missing Frame Acquisition DateTime.\n");
        tacFree(&tac); iftFree(&fl); dcmfileFree(&dcm);
        return(6);
      }
      if(verbose>20) printf("facqDateTime := %s\n", facqDateTime); 
      /* Calculate the Frame start time */
      double t1;
      if(decayCorrection>0) {
        /* Decay Correction DatTime is the reference time */
        t1=strDateTimeDifference(facqDateTime, dcDateTime);
      } else {
        /* Use series time, for now */
        t1=strDateTimeDifference(facqDateTime, acqDateTime);
      }
      /* Set frame start time */
      if(frameIndex<frameNr) {
        if(isnan(tac.x1[frameIndex])) {
          /* Save frame start time, since there is no previous value for this frame */
          tac.x1[frameIndex]=t1;
          if(verbose>29) printf("t1[%d]=%g\n", frameIndex, t1);
        } else {
          /* Check whether previous frame start time is the same as current */
          if(!doubleMatch(tac.x1[frameIndex], t1, 0.001)) {
            beds=1; // not the same, this probably is whole-body study with more than one bed pos
            if(verbose>20) printf("t1=%g\n", t1);
          }
        }
      }
 
      /* Find Frame Acquisition Duration */
      jptr=dcmFindDownTag(iptr->child_item, 0, &tagdur, 0);
      if(jptr==NULL) {fptr=fptr->next_item; continue;}
      if(verbose>20) dcmitemPrint(jptr);
      /* For now, put frame duration in place of frame end time */
      if(frameIndex<frameNr) tac.x2[frameIndex]=0.001*dcmitemGetReal(jptr);
      /* prepare for the next frame */
      fptr=iptr->next_item;
    }
  }
 
 
 
  /* Set scan start time */
  if(decayCorrection==2) {
    /* Decay is corrected to injection time, therefore put it here too */
    tacSetHeaderScanstarttime(&tac.h, injDateTime);
  } else {
    /* Decay is not corrected, or is corrected to acquisition start time */
    tacSetHeaderScanstarttime(&tac.h, startDateTime);
    /* Correct frame start times, too */
    for(int i=0; i<tac.sampleNr; i++)
      tac.x1[i]=tac.x1[i]-strDateTimeDifference(startDateTime, seriesDateTime);
  }
 
  /* Calculate frame end time */
  for(int i=0; i<tac.sampleNr; i++) tac.x2[i]+=tac.x1[i];
 
 
  /* 
   *  Write SIF data stored in TAC struct in SIF format
   */
  FILE *fp=stdout;
  if(siffile[0]) {
    if(verbose>1) printf("writing %s\n", siffile);
    fp=fopen(siffile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&tac); return(11);
    }
  }
  ret=tacWriteSIF(&tac, fp, 0, &status);
  tacFree(&tac); if(siffile[0]) fclose(fp);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    return(11);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/