imagedcm.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include "tpcimage.h"
#include "tpcfileutil.h"
/*****************************************************************************/
//#include <sys/stat.h>
//#include <unistd.h>
/*****************************************************************************/
 
/*****************************************************************************/
int imgReadDICOM(
  IMG *img,
  const char *fname,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(..., %s, ...)\n", __func__, fname); fflush(stdout);}
 
  if(img==NULL || strnlen(fname, 2)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
 
  /* Delete any old contents in data structure */
  imgFree(img);
 
  /* Read DICOM file */
  DCMFILE dcm; dcmfileInit(&dcm);
  int ret=dcmFileRead(fname, &dcm, 0, status);
  if(ret!=TPCERROR_OK) {dcmfileFree(&dcm); return(ret);}
 
//verbose=20;
 
  /* Get modality */
  imgmodality modality=IMG_MODALITY_UNKNOWN;
  {
    DCMTAG tag; dcmTagSet(&tag, 0x0008, 0x0060); DCMITEM *iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr==NULL) {
      if(verbose>0) {fprintf(stderr, "Error: cannot find modality.\n"); fflush(stderr);}
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
      dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
    }
    char *buf=dcmValueString(iptr); 
    if(verbose>12) printf("Modality := '%s'\n", buf);
    if(strncmp(buf, "PT", 2)==0) modality=IMG_MODALITY_PET;
    else if(strncmp(buf, "CT", 2)==0) modality=IMG_MODALITY_CT;
    else if(strncmp(buf, "MR", 2)==0) modality=IMG_MODALITY_MRI;
    else if(strncmp(buf, "NM", 2)==0) modality=IMG_MODALITY_SPECT;
    else if(verbose>0) fprintf(stderr, "Warning: modality '%s' not supported.\n", buf);
    free(buf);
  }
 
  /* Get the number of frames and planes */
  if(verbose>10) {printf("reading matrix dimensions\n"); fflush(stdout);}
  unsigned short int imgdim[4];
  ret=dcmImgDim(&dcm, imgdim, verbose-2);
  if(ret!=0 || imgdim[0]<=0 || imgdim[1]<=0 || imgdim[2]<=0 || imgdim[3]<=0) {
    if(verbose>0) {fprintf(stderr, "Error: cannot find matrix dimensions.\n"); fflush(stderr);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
    dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
  }
  unsigned short int dimx, dimy, sliceNr, frameNr;
  dimx=imgdim[0]; dimy=imgdim[1]; sliceNr=imgdim[2]; frameNr=imgdim[3];
  if(verbose>11) 
    printf("dimx := %u\ndimy := %u\ndimz := %u\ndimt := %u\n", dimx, dimy, sliceNr, frameNr); 
 
  /* 
   *  Allocate memory for IMG data
   */
  ret=imgAllocate(img, sliceNr, dimy, dimx, frameNr, status);
  if(ret!=TPCERROR_OK) {dcmfileFree(&dcm); return(ret);}
 
 
  /* Get the status of decay correction and isotope */
  ret=dcmImgIsotope(&dcm, &img->isot, &img->decayCorrection, verbose-2);
  if(ret!=0 && verbose>0) {
    fprintf(stderr, "Warning: decay correction and isotope unknown.\n"); fflush(stderr);
  }
  if(verbose>11) {
    printf("decayCorrection := %s\n", decayDescr(img->decayCorrection));
    printf("isotope := %s\n", isotopeName(img->isot));
    fflush(stdout);
  }
 
  /* If decay is corrected, read time of decay correction; otherwise read start time */
  char zeroDateTime[32]; zeroDateTime[0]=(char)0;
  if(img->decayCorrection>=DECAY_CORRECTED) {
    DCMTAG tag; dcmTagSet(&tag, 0x0018, 0x9701);
    DCMITEM *iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr!=NULL) dcmDT2intl(iptr->rd, zeroDateTime);
  } else {
    /* Get Acquisition DateTime */
    DCMTAG tag; dcmTagSet(&tag, 0x0008, 0x002A);
    DCMITEM *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);
    imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
  }
  if(verbose>14) printf("zeroDateTime := %s\n", zeroDateTime);
 
  /* Get frame times */
  if(modality==IMG_MODALITY_PET || modality==IMG_MODALITY_SPECT) {
    /* Initiate frame times to NaN */
    for(unsigned short int i=0; i<img->dimt; i++) img->x1[i]=img->x2[i]=img->x[i]=nanf("");
 
    /* Find Per Frame Functional Groups Sequence */
    DCMTAG tag; dcmTagSet(&tag, 0x5200, 0x9230);
    DCMITEM *iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr==NULL) {
      if(verbose>0) fprintf(stderr, "Error: Per Frame Functional Groups Sequence not found.\n");
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
      imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
    }
    /* Loop through all Frame Content Sequences under it */
    dcmTagSet(&tag, 0x0020, 0x9111);
    DCMTAG tagstart; dcmTagSet(&tagstart, 0x0018, 0x9074);
    DCMTAG tagdur; dcmTagSet(&tagdur, 0x0018, 0x9220);
    DCMTAG tagfr; dcmTagSet(&tagfr, 0x0020, 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>20) 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>30) 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]) {
        if(verbose>0) fprintf(stderr, "Error: missing Frame Acquisition DateTime.\n");
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
        imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
      }
      if(verbose>20) printf("facqDateTime := %s\n", facqDateTime); 
      /* Calculate the Frame start time */
      float t1=(float)strDateTimeDifference(facqDateTime, zeroDateTime);
      /* Set frame start time */
      if(frameIndex<img->dimt) {
        if(isnan(img->x1[frameIndex])) {
          /* Save frame start time, since there is no previous value for this frame */
          img->x1[frameIndex]=t1;
          if(verbose>29) printf("t1[%u]=%g\n", frameIndex, t1);
        } else {
          /* Check whether previous frame start time is the same as current */
          if(!floatMatch(img->x1[frameIndex], t1, 0.001)) {
            // not the same, this probably is whole-body study with more than one bed pos
            if(verbose>1) printf("  different 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);
      if(frameIndex<img->dimt) {
        img->x2[frameIndex]=img->x1[frameIndex]+(float)0.001*dcmitemGetReal(jptr);
        img->x[frameIndex]=0.5*(img->x1[frameIndex]+img->x2[frameIndex]);
      }
      /* prepare for the next frame */
      fptr=iptr->next_item;
    }
  }
 
  /* Get pixel sizes and image position */
  {
    double pxlsize[3];
    ret=dcmImgPxlsize(&dcm, pxlsize, verbose-2);
    if(ret!=0) {
      if(verbose>0) fprintf(stderr, "Error: missing Pixel size.\n");
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
      imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
    }
    img->sizex=(float)pxlsize[0];
    img->sizey=(float)pxlsize[1];
    img->sizez=(float)pxlsize[2];
    if(verbose>13) printf("  pixel_size := %g x %g x %g\n", img->sizex, img->sizey, img->sizez);
 
    double ipp[3];
    int ret1=dcmImgPos(&dcm, ipp, verbose-2);
    if(ret1==0) {
      for(int i=0; i<3; i++) img->ipp[i]=ipp[i];
    } else {
      if(verbose>0) fprintf(stderr, "Warning: missing Image Position (Patient).\n");
    }
    double iop[6];
    int ret2=dcmImgOrient(&dcm, iop, verbose-2);
    if(ret2==0) {
      for(int i=0; i<6; i++) img->iop[i]=iop[i];
    } else {
      if(verbose>0) fprintf(stderr, "Warning: missing Image Orientation (Patient).\n");
    }
    if(ret1==0 && ret2==0) { // calculate xform into IMG
      double xform[16];
      ret=dcmImgXform(iop, pxlsize, ipp, xform, verbose-2);
      if(ret==0) {
        img->xform[0]=1;
        for(int i=0; i<12; i++) img->srow[i]=(float)xform[i];
        double quatern[3], qoffset[3];
        if(dcmXformToQuatern(xform, quatern, qoffset, verbose-2)==0) {
          for(int i=0; i<3; i++) img->quatern[i]=(float)quatern[i];
          for(int i=0; i<3; i++) img->quatern[3+i]=(float)qoffset[i];
        }
      }
    }
  }
 
 
  /* Get Accession Number and write it in IMG as studyNr */
  {
    img->studyNr[0]=(char)0;
    DCMITEM *iptr; DCMTAG tag; dcmTagSet(&tag, 0x0008, 0x0050);
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr!=NULL) {
      char *buf=dcmValueString(iptr); if(verbose>12) printf("Accession Number := '%s'\n", buf);
      if(buf!=NULL && strlen(buf)>1) strncpyClean(img->studyNr, buf, MAX_STUDYNR_LEN+1);
      free(buf);
    }
    if(!strcasecmp(img->studyNr, "empty")) img->studyNr[0]=(char)0;
  }
  /* If not successful, then try Patient ID */
  if(!img->studyNr[0]) {
    DCMITEM *iptr; DCMTAG tag; dcmTagSet(&tag, 0x0010, 0x0020);
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr!=NULL) {
      char *buf=dcmValueString(iptr); if(verbose>12) printf("Patient ID := '%s'\n", buf);
      if(buf!=NULL && strlen(buf)>1) strncpyClean(img->studyNr, buf, MAX_STUDYNR_LEN+1);
      free(buf);
    }
  } 
 
 
  /* Get pixel data scaling factors */
  double rescale_intercept=0.0, rescale_slope=0.0;
  {
    DCMTAG tag; DCMITEM *iptr, *jptr;
    dcmTagSet(&tag, 0x0040, 0x9096);
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr==NULL) {
      if(verbose>0) fprintf(stderr, "Error: missing Real World Value Mapping Sequence.\n");
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
      imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
    }
    dcmTagSet(&tag, 0x0040, 0x9224);
    jptr=dcmFindDownTag(iptr->child_item, 0, &tag, 0);
    if(jptr!=NULL) rescale_intercept=dcmitemGetReal(jptr);
    dcmTagSet(&tag, 0x0040, 0x9225);
    jptr=dcmFindDownTag(iptr->child_item, 0, &tag, 0);
    if(jptr!=NULL) rescale_slope=dcmitemGetReal(jptr);
    if(rescale_slope<=0.0 && verbose>0) {
      fprintf(stderr, "Warning: rescale_slope := %g\n", rescale_slope); fflush(stderr);
    } else if(verbose>14) {
      printf("rescale_slope := %g\n", rescale_slope);
      printf("rescale_intercept := %g\n", rescale_intercept); // always 0 for PET
    }
    /* Unit */
    dcmTagSet(&tag, 0x0040, 0x9210);
    jptr=dcmFindDownTag(iptr->child_item, 0, &tag, 0);
    if(jptr!=NULL) {
      char *buf=dcmValueString(jptr); if(verbose>14) printf("  buf := '%s'\n", buf);
      strClean(buf); img->cunit=unitIdentify(buf);
      free(buf);
    }
    if(verbose>14) printf("image_unit := %s\n", unitName(img->cunit));
  }
 
  /* Fill IMG pixel values */
  {
    DCMTAG tag; dcmTagSet(&tag, 0x7fe0, 0x0010);
    DCMITEM *iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr==NULL || iptr->rd==NULL) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
      imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_DATA);
    }
    unsigned long long bi=0, byteNr=iptr->vl;
    if(verbose>13) printf("  pixel data byteNr := %llu\n", byteNr);
 
    char *cptr=iptr->rd;
    short int s;
    for(unsigned int ti=0; ti<img->dimt; ti++) {
      for(unsigned int zi=0; zi<img->dimz; zi++)
        for(unsigned int yi=0; yi<img->dimy; yi++)
          for(unsigned int xi=0; xi<img->dimx; xi++) {
            if(bi>=byteNr) {printf("missing pixel data\n"); fflush(stdout); break;}
            memcpy(&s, cptr, 2); cptr+=2; bi+=2;
            img->m[zi][yi][xi][ti]=(float)(rescale_slope*(double)s + rescale_intercept);
          }
    }
  }
 
 
  dcmfileFree(&dcm);
 
  /* Set IMG header fields */
  img->modality=modality;
  img->format=IMG_FORMAT_DICOM;
  img->content=IMG_CONTENT_IMAGE;
  strcpy(img->scanStart, zeroDateTime);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int imgWriteDICOM(
  IMG *img,
  const char *fname,
  TPCSTATUS *status
) {
  if(status==NULL) return(TPCERROR_FAIL);
  if(status->verbose>1) {printf("%s(..., %s, ...)\n", __func__, fname); fflush(stdout);}
  if(img==NULL || img->dimt<1 || img->dimz<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA); return(status->error);}
  if(fname==NULL || strnlen(fname, 1)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME); return(status->error);}
 
  /*
   *  Make DCMFILE data structure
   */
  DCMFILE dcm; dcmfileInit(&dcm);
  /* Set Transfer UID */
  dcm.truid=DCM_TRUID_LEE;
  {
    DCMTAG tag;
    unsigned int ul;
    char buf[256];
 
    /* Meta Element Group */
    dcmTagSet(&tag, 0x0002, 0x0000); ul=0;
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UL, 0xFFFFFFFF, (char*)&ul, 0);
    dcmTagSet(&tag, 0x0002, 0x0001); // File Meta Information Version
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_OB, 2, "\0\1", 0);
    dcmTagSet(&tag, 0x0002, 0x0002); // Media Stored SOP Class UID
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, "1.2.840.10008.5.1.4.1.1.128", 0);
    dcmTagSet(&tag, 0x0002, 0x0003); // Media Stored SOP Instance UID
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, "1.2.826.0.1.3417726.3.0.20180616", 0);
    dcmTagSet(&tag, 0x0002, 0x0010); // Transfer Syntax UID
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, dcmTrUIDString(dcm.truid), 0);
    dcmTagSet(&tag, 0x0002, 0x0012); // Implementation Class UID
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, "1.3.6.1.4.1.25403.1.1.1", 0);
    dcmTagSet(&tag, 0x0002, 0x0013); // Implementation Version Name
    snprintf(buf, 256, "TPCCLIB_%d.%d.%d", tpcclib_VERSION_MAJOR, tpcclib_VERSION_MINOR, tpcclib_VERSION_PATCH);
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_SH, 0xFFFFFFFF, buf, 0);
 
    /* Identifying Group */
    dcmTagSet(&tag, 0x0008, 0x0000); ul=0;
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UL, 0xFFFFFFFF, (char*)&ul, 0);
    dcmTagSet(&tag, 0x0008, 0x0008); // Image Type
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_CS, 17, "ORIGINAL/PRIMARY", 0);
    {
      time_t t=time(NULL);
      struct tm *local_tm; local_tm=localtime(&t);
      dcmTagSet(&tag, 0x0008, 0x0012); // Instance Creation Date
      strftime(buf, 256, "%Y%m%d", local_tm);
      dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_DA, 0xFFFFFFFF, buf, 0);
      dcmTagSet(&tag, 0x0008, 0x0013); // Instance Creation Time
      strftime(buf, 256, "%H%M%S", local_tm);
      dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_DA, 0xFFFFFFFF, buf, 0);
    }
    dcmTagSet(&tag, 0x0008, 0x0016); // SOP Class UID; see dcmsop.c
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, "1.2.840.10008.5.1.4.1.1.128", 0);
    dcmTagSet(&tag, 0x0008, 0x0018); // SOP Instance UID
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, "1.2.826.0.1.3417726.3.0.20180616", 0);
    dcmTagSet(&tag, 0x0008, 0x0060); // Modality
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_CS, 3, "PT", 0);
 
    /* Acquisition Group */
    dcmTagSet(&tag, 0x0018, 0x0000); ul=0;
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UL, 0xFFFFFFFF, (char*)&ul, 0);
    dcmTagSet(&tag, 0x0018, 0x0050);  // Slice thickness (mm)
    snprintf(buf, 256, "%g", img->sizez);
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_DS, 0xFFFFFFFF, buf, 0);
    dcmTagSet(&tag, 0x0018, 0x1075); // half-life
    if(img->isot!=ISOTOPE_UNKNOWN) {
      snprintf(buf, 256, "%g", 60.0*isotopeHalflife(img->isot));
      dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_DS, 0xFFFFFFFF, buf, 0);
    }
 
    /* Relationship Group */
    dcmTagSet(&tag, 0x0020, 0x0000); ul=0;
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UL, 0xFFFFFFFF, (char*)&ul, 0);
 
 
 
  }
 
 
  /* Create folder, or if it already exists check that it is empty */
  if(pathExist(fname)) {
    if(status->verbose>2) {
      fprintf(stdout, "  checking that folder '%s' is empty\n", fname); fflush(stdout);}
    if(pathFileNr(fname)>0) {
      dcmfileFree(&dcm);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FILE_EXISTS);
      return(status->error);
    }
  } else {
    if(status->verbose>2) {
      fprintf(stdout, "  creating subdirectory '%s'\n", fname); fflush(stdout);}
    if(pathCreate(fname)) {
      dcmfileFree(&dcm);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
      return(status->error);
    }
  }
 
  /*
   *  Write each image plane in its own DICOM file
   */
  for(int zi=0; zi<img->dimz; zi++) {
    char dname[strlen(fname)+25];
    snprintf(dname, strlen(fname)+25, "%s/%0*d.dcm", fname, 1+(int)floor(log10(img->dimz)), 1+zi);
    if(status->verbose>2) {printf("  writing '%s'\n", dname); fflush(stdout);}
 
    /* Write */
    dcmFileWrite(dname, &dcm, status);
    if(status->error!=TPCERROR_OK) {
      dcmfileFree(&dcm);
      dcmFileRemove(dname, NULL);
      pathRemove(fname);
      return(status->error);
    }
  } // next plane
  dcmfileFree(&dcm);
 
  if(0) { // Remove files and empty folder in case of error
    dcmFileRemove(fname, NULL);
    pathRemove(fname);
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/