dcmmatrix.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 "tpcdcm.h"
/*****************************************************************************/
 
/*****************************************************************************/
void dcmmatrixInit(
  DCMMATRIX *m
) {
  if(m==NULL) return;
  m->filename=NULL;
  m->acqDate[0]=m->acqTime[0]=(char)0;
  m->frame=m->plane=0;
  m->frameStart=m->frameDur=0.0;
}
/*****************************************************************************/
 
/*****************************************************************************/
void dcmmatrixFree(
  DCMMATRIX *m
) {
  if(m==NULL) return;
  if(m->filename!=NULL) free(m->filename);
  dcmmatrixInit(m);
}
/*****************************************************************************/
 
/*****************************************************************************/
void dcmmlInit(
  DCMML *d
) {
  if(d==NULL) return;
  d->nr=d->anr=0;
  d->m=(DCMMATRIX*)NULL;
}
/*****************************************************************************/
 
/*****************************************************************************/
void dcmmlFree(
  DCMML *d
) {
  if(d==NULL) return;
  if(d->m!=NULL) {
    for(unsigned int i=0; i<d->anr; i++) dcmmatrixFree(d->m+i);
    free(d->m);
  }
  dcmmlInit(d);
}
/*****************************************************************************/
 
/*****************************************************************************/
int dcmmlAllocate(
  DCMML *d,
  int mNr
) {
  if(d==NULL) return TPCERROR_FAIL;
  /* Check if there is enough space already */
  unsigned int newNr, addNr;
  newNr=d->nr+mNr; 
  if(d->anr>=newNr) return TPCERROR_OK;
  addNr=newNr-d->anr;
 
  /* (Re)allocate memory. */
  DCMMATRIX *mptr;
  mptr=(DCMMATRIX*)realloc(d->m, sizeof(DCMMATRIX)*newNr);
  if(mptr==NULL) return TPCERROR_OUT_OF_MEMORY;
  d->m=mptr;
  mptr=d->m+d->anr; for(unsigned int i=0; i<addNr; i++) dcmmatrixInit(mptr++);
  d->anr=newNr;
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int dcmMListRead(
  IFT *ift,
  DCMML *ml,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>1) printf("%s()\n", __func__);
  if(ift==NULL || ml==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(ift->keyNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Delete any previous data */
  dcmmlFree(ml);
 
  /* 
   *  Read the first file in list
   */
  int fi=0, ret=TPCERROR_OK;
  DCMFILE dcm; dcmfileInit(&dcm);
  ret=dcmFileRead(ift->item[fi].value, &dcm, 1, status);
  if(ret!=TPCERROR_OK) {dcmfileFree(&dcm); return(ret);}
 
  /* Read some information common to the whole image */
  DCMTAG tag;
  DCMITEM *iptr;
 
  /* 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) {
    ret=TPCERROR_NO_KEY; statusSet(status, __func__, __FILE__, __LINE__, ret);
    dcmfileFree(&dcm); return(ret);
  }
  /* 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>2) 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;
  }
 
 
  if(multiframe) {  /* Multi-frame file */
 
    /* Get the number of frames and planes */
    if(verbose>1) {printf("reading frame and slice number\n"); fflush(stdout);}
    unsigned short int frameNr=0;
    /* Get frame number 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;
    if(verbose>1) printf("frameNr := %u\n", frameNr); 
    unsigned short int sliceNr=0;
    /* Get plane number from 'In Stack Position Number' */
    tag.group=0x0020; tag.element=0x9057;
    if(dcmTagIntRange(dcm.item, &tag, &a, &b, verbose-100)==0) sliceNr=b;
    else {
      /* Get plane number from 'Dimension Index Values' */
      tag.group=0x0020; tag.element=0x9157;
      if(dcmTagIntRange(dcm.item, &tag, &a, &b, verbose-1)==0) sliceNr=b;
    }
    if(verbose>1) printf("sliceNr := %u\n", sliceNr); 
    if(frameNr<=0 || sliceNr<=0) {
      if(verbose>0) fprintf(stderr, "Error: cannot find frame/slice number.\n");
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
      dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
    }
 
    /* Get Decay Correction */
    short unsigned int decayCorrection=0; // 0=NO, 1=YES
    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, read time of decay correction; otherwise read start time */
    char zeroDateTime[32]; zeroDateTime[0]=(char)0;
    if(decayCorrection==0) {
      /* Get Acquisition DateTime */
      tag.group=0x0008; tag.element=0x002A;
      iptr=dcmFindTag(dcm.item, 0, &tag, 0);
      if(iptr!=NULL) dcmDT2intl(iptr->rd, zeroDateTime);
    } else {
      tag.group=0x0018; tag.element=0x9701;
      iptr=dcmFindTag(dcm.item, 0, &tag, 0);
      if(iptr!=NULL) dcmDT2intl(iptr->rd, zeroDateTime);
    }
    if(!zeroDateTime[0]) {
      if(verbose>0) fprintf(stderr, "Error: missing Acquisition Date and Time.\n");
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
      dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
    }
    if(verbose>4) printf("zeroDateTime := %s\n", zeroDateTime);
 
    /* Initially allocate for frameNr x sliceNr matrices */
    ret=dcmmlAllocate(ml, sliceNr*frameNr);
    if(ret!=TPCERROR_OK) {
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      dcmfileFree(&dcm); return ret;
    }
 
    /*
     *  Get the information from one matrix at a time
     */
    /* Find Per Frame Functional Groups Sequence */
    tag.group=0x5200; tag.element=0x9230;
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr==NULL) {
      if(verbose>0) {
        fprintf(stderr, "Error: cannot find Per Frame Functional Groups Sequence.\n"); 
        fflush(stderr);
      }
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
      dcmmlFree(ml); dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
    }
    /* 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; // not an error!
      if(verbose>10) printf(" found Frame Content Sequence\n");
      /* Check that we have room */
      if(ml->nr>=ml->anr) {
        ret=dcmmlAllocate(ml, 100);
        if(ret!=TPCERROR_OK) break;
      }
      /* Set file name */
      ml->m[ml->nr].filename=strdup(ift->item[0].value);
      /* Find Temporal Position Index */
      DCMITEM *jptr;
      jptr=dcmFindDownTag(iptr->child_item, 0, &tagfr, 0);
      if(jptr==NULL) {
        if(verbose>0) {printf("cannot find Temporal Position Index.\n"); fflush(stdout);}
        ret=TPCERROR_NO_KEY; break;
      }
      if(verbose>20) dcmitemPrint(jptr);
      ml->m[ml->nr].frame=(unsigned int)dcmitemGetInt(jptr);
      //unsigned short int frameIndex=(unsigned short int)(dcmitemGetInt(jptr)-1);
      /* Find slice index */
      DCMTAG tagpl; tagpl.group=0x0020; tagpl.element=0x9057;
      jptr=dcmFindDownTag(iptr->child_item, 0, &tagpl, 0);
      if(jptr==NULL) {
        tagpl.group=0x0020; tagpl.element=0x9157;
        jptr=dcmFindDownTag(iptr->child_item, 0, &tagpl, 0);
      }
      if(jptr==NULL) {
        if(verbose>0) {printf("cannot find Slice Index.\n"); fflush(stdout);}
        ret=TPCERROR_NO_KEY; break;
      }
      if(verbose>20) dcmitemPrint(jptr);
      ml->m[ml->nr].plane=(unsigned int)dcmitemGetInt(jptr);
      /* 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]) {
        if(verbose>0) {printf("cannot find Frame Acquisition DateTime.\n"); fflush(stdout);}
        ret=TPCERROR_NO_KEY; break;
      }
      if(verbose>20) printf("facqDateTime := %s\n", facqDateTime);
      strlcpy(ml->m[ml->nr].acqDate, facqDateTime, 11);
      if(strlen(facqDateTime)>17) strlcpy(ml->m[ml->nr].acqTime, facqDateTime+11, 9);
      /* Calculate the Frame start time */
      ml->m[ml->nr].frameStart=strDateTimeDifference(facqDateTime, zeroDateTime);
      /* Find Frame Acquisition Duration */
      jptr=dcmFindDownTag(iptr->child_item, 0, &tagdur, 0);
      if(jptr!=NULL) ml->m[ml->nr].frameDur=0.001*dcmitemGetReal(jptr);
 
      /* prepare for the next matrix */
      fptr=iptr->next_item;
      ml->nr++;
    }
 
    dcmfileFree(&dcm);
    if(ret!=TPCERROR_OK) {
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      dcmmlFree(ml); 
      return(ret);
    }
 
  } else {   /* Single-frame files */
 
    /* Get the number of frames (time slices) */
    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 {
      if(dynamic!=0) { // error, if dynamic
        ret=TPCERROR_NO_KEY; statusSet(status, __func__, __FILE__, __LINE__, ret);
        dcmfileFree(&dcm); return(ret);
      } 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) */
    unsigned short int sliceNr=0;
    tag.group=0x0054; tag.element=0x0081;
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr==NULL) {
      ret=TPCERROR_NO_KEY; statusSet(status, __func__, __FILE__, __LINE__, ret);
      dcmfileFree(&dcm); return(ret);
    }
    sliceNr=(unsigned short int)dcmitemGetInt(iptr);
    if(verbose>1) printf("sliceNr := %u\n", sliceNr); 
 
    /* Get Decay correction type */
    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(strcasecmp(buf, "NONE")==0) decayCorrection=0;
      else if(strcasecmp(buf, "START")==0) decayCorrection=1;
      else if(strcasecmp(buf, "ADMIN")==0) decayCorrection=2;
      free(buf);
    }
 
    /* 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;
      }
      /* 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);
          }
        }
      }
      if(injDateTime[0] && verbose>3) 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] && verbose>3) printf("seriesDateTime := %s\n", seriesDateTime); 
 
 
    /* Initially allocate for frameNr x sliceNr matrices */
    ret=dcmmlAllocate(ml, sliceNr*frameNr);
    if(ret!=TPCERROR_OK) {
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      dcmfileFree(&dcm); return ret;
    }
 
 
    /*
     *  Get the information from one file at a time
     */
    //char startDateTime[32]; startDateTime[0]=(char)0;
    //int fi=0, beds=0;
    do {
 
      /* 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) {ret=TPCERROR_NO_KEY; break;}
      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; 
      }
      /* Set file name, plane, and frame into matrix list */
      if(ml->nr>=ml->anr) {
        ret=dcmmlAllocate(ml, 100);
        if(ret!=TPCERROR_OK) break;
      }
      ml->m[ml->nr].filename=strdup(ift->item[fi].value);
      ml->m[ml->nr].plane=1+sliceIndex;
      ml->m[ml->nr].frame=1+frameIndex;
 
      /* 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]) {
        strlcpy(ml->m[ml->nr].acqDate, acqDateTime, 11);
        if(strlen(acqDateTime)>17) strlcpy(ml->m[ml->nr].acqTime, acqDateTime+11, 9);
      }
 
      /* Calculate the Frame start time, if possible */
      if(acqDateTime[0]) {
        double t1=nan("");
        if(decayCorrection==2 && injDateTime[0]) {
          /* Injection time is the reference time */
          t1=strDateTimeDifference(acqDateTime, injDateTime);
          ml->m[ml->nr].frameStart=t1;
        } else if(seriesDateTime[0]) {
          /* Use series time, for now */
          t1=strDateTimeDifference(acqDateTime, seriesDateTime);
          ml->m[ml->nr].frameStart=t1;
        }
      }
 
      /* Get Actual Frame Duration */
      tag.group=0x0018; tag.element=0x1242;
      iptr=dcmFindTag(dcm.item, 0, &tag, 0);
      if(iptr!=NULL) ml->m[ml->nr].frameDur=0.001*dcmitemGetReal(iptr);
 
 
      /* Done with this file */
      ml->nr++;
      dcmfileFree(&dcm);
 
      /* Read the next file in list */
      fi++; if(fi>=ift->keyNr) {ret=TPCERROR_OK; break;}
      ret=dcmFileRead(ift->item[fi].value, &dcm, 1, status);
      if(ret!=TPCERROR_OK) break;
 
    } while(1);
    dcmfileFree(&dcm);
    if(ret!=TPCERROR_OK) {
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      dcmmlFree(ml); 
      return(ret);
    }
  }
 
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
 
static int dcmmlQSortPlane(const void *m1, const void *m2)
{
  if(((DCMMATRIX*)m1)->plane < ((DCMMATRIX*)m2)->plane) return(-1);
  else if(((DCMMATRIX*)m1)->plane > ((DCMMATRIX*)m2)->plane) return(+1);
  else if(((DCMMATRIX*)m1)->frame < ((DCMMATRIX*)m2)->frame) return(-1);
  else if(((DCMMATRIX*)m1)->frame > ((DCMMATRIX*)m2)->frame) return(+1);
  else return(0);
}
 
int dcmmlSortByPlane(
  DCMML *d,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  /* Check that required data exists */
  if(d==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  if(d->nr<2) return TPCERROR_OK;
  qsort(d->m, d->nr, sizeof(DCMMATRIX), dcmmlQSortPlane);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/