tac2nii.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpccsv.h"
#include "tpctac.h"
#include "tpcdcm.h"
#include "tpcecat.h"
#include "tpcnifti.h"
#include "tpcimage.h"
/*****************************************************************************/
 
/*****************************************************************************/
int simRCDims(
  const long long N,
  const int mxdim,
  const int mydim,
  const int mzdim,
  int *rnx,
  int *rny,
  int *rnz,
  int *rcxdim,
  int *rcydim,
  int *rczdim
);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Create a 4D PET image file in NIfTI 1S format, with contents from",
  "the user-specified TAC file, for software testing.",
  "Image volume is divided into rectangular cuboids, containing as pixel",
  "values the regional TAC values.",
  "Frame times are written in SIF as specified in the TAC file.",
  " ",
  "Usage: @P [Options] tacfile xdim ydim zdim imagefile [template]",
  " ",
  "If file name for template NIfTI is given, the rectangular cuboids are",
  "are saved in it, with integer pixel values starting from 1.",
  " ",
  "Options:",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: dft2img, flat2img, img2tif, img2dft, simboxes, pxl2mask",
  " ",
  "Keywords: image, NIfTI, simulation, software testing",
  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  tacfile[FILENAME_MAX], dbname[FILENAME_MAX], tname[FILENAME_MAX];
  int   dimx=0, dimy=0, dimz=0;
  int   ret=0;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=dbname[0]=tname[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* options */
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
  
  /* 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 */
  ret=0;
  if(ai<argc) {strlcpy(tacfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {if(atoiCheck(argv[ai++], &dimx)) ret++;}
  if(ai<argc) {if(atoiCheck(argv[ai++], &dimy)) ret++;}
  if(ai<argc) {if(atoiCheck(argv[ai++], &dimz)) ret++;}
  if(ret || dimx<1 || dimy<1 || dimz<1) {
    fprintf(stderr, "Error: invalid dimension.\n");
    return(1);  
  }
  if(ai<argc) {strlcpy(dbname, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(tname, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing or wrong? */
  if(!dbname[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("tacfile := %s\n", tacfile);
    printf("dbname := %s\n", dbname);
    if(tname[0]) printf("tname := %s\n", tname);
    printf("dimx := %d\n", dimx);
    printf("dimy := %d\n", dimy);
    printf("dimz := %d\n", dimz);
    fflush(stdout);
  }
 
  /*
   *  Read regional data
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  TAC tac;
  tacInit(&tac);
  ret=tacRead(&tac, tacfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    printf("isframe := %d\n", tac.isframe);
    fflush(stdout);
  }
  /* Get the min and max pixel values for the header and scaling */
  double min, max;
  if(tacYRange(&tac, -1, &min, &max, NULL, NULL, NULL, NULL)) {
    fprintf(stderr, "Error: invalid TAC contents.\n");
    tacFree(&tac);
    return(2);
  }
  if(verbose>2) printf("min := %g\nmax := %g\n", min, max);
 
 
  /*
   *  Determine dimensions of rectangular cuboids
   */
  if(verbose>1) printf("Determine dimensions of rectangular cuboids\n");
  int rdimx=0, rdimy=0, rdimz=0;
  int rnx=0, rny=0, rnz=0;
  if(simRCDims(tac.tacNr, dimx, dimy, dimz, &rnx, &rny, &rnz, &rdimx, &rdimy, &rdimz)) {
    fprintf(stderr, "Error: incompatible input data.\n");
    tacFree(&tac);
    return(3);
  }
  if(verbose>3)
    printf("%d into %d x %d x %d -> %d x %d x %d (%d x %d x %d)\n", 
           tac.tacNr, dimx, dimy, dimz, rdimx, rdimy, rdimz, rnx, rny, rnz);
 
 
  /*
   *  Make template data
   */
  if(verbose>1) printf("Allocate memory for template xyz matrix\n");
  size_t pxlNr=(size_t)dimz*dimy*dimx;
  int *idata;
  idata=(int*)calloc(pxlNr, sizeof(int));
  if(idata==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    tacFree(&tac);
    return(4);
  }
  /* Fill the template with TAC numbers */
  if(verbose>1) printf("Fill the template\n");
  {
    int xi, yi, zi;
    int x[2], y[2], z[2];
    xi=yi=zi=0;
    for(int ri=0; ri<tac.tacNr; ri++) {
      x[0]=xi*rdimx; x[1]=x[0]+rdimx-1;
      y[0]=yi*rdimy; y[1]=y[0]+rdimy-1;
      z[0]=zi*rdimz; z[1]=z[0]+rdimz-1;
      if(verbose>4)
        printf("  %s : %d,%d,%d - %d,%d,%d\n", tac.c[ri].name, x[0],y[0],z[0], x[1],y[1],z[1]);
      for(int cz=z[0]; cz<=z[1]; cz++)
        for(int cy=y[0]; cy<=y[1]; cy++)
          for(int cx=x[0]; cx<=x[1]; cx++) {
            size_t pos=((size_t)dimy*dimx)*cz + dimx*cy + cx;
            idata[pos]=1+ri;
          }
      xi++; if(xi==rnx) {xi=0; yi++; if(yi==rny) {yi=0; zi++; if(zi==rnz) break;}}
    }
  }
 
 
  /* 
   *  Set NIfTI header contents, first for the template
   */
  if(verbose>1) printf("Fill NIfTI header\n");
  NIFTI_DSR dsr;
  dsr.n=1;
  /* Set NIfTI byte order to current machines byte order */
  dsr.byte_order=endianLittle();
  /* Initiate header structures with zeroes */
  memset(&dsr.h1, 0, sizeof(NIFTI_1_HEADER));
  memset(&dsr.e, 0, sizeof(NIFTI_EXTENDER));
  /* Set header */
  dsr.h1.sizeof_hdr=NIFTI1_HEADER_SIZE;
  strcpy(dsr.h1.data_type, "");
  strlcpy(dsr.h1.db_name, dbname, 17);
  dsr.h1.extents=16384; // not used in NIfTI, but required for Analyze compatibility
  dsr.h1.regular='r'; // not used in NIfTI, but required for Analyze compatibility
  dsr.h1.dim_info='\0'; // MRI slice ordering
  /* Image dimension */
  for(int i=0; i<8; i++) dsr.h1.dim[i]=1;
  dsr.h1.dim[0]=3;
  dsr.h1.dim[1]=dimx;
  dsr.h1.dim[2]=dimy;
  dsr.h1.dim[3]=dimz;
  dsr.h1.dim[4]=1;
  dsr.h1.intent_p1=0.0;
  dsr.h1.intent_p2=0.0;
  dsr.h1.intent_p3=0.0;
  dsr.h1.intent_code=NIFTI_INTENT_NONE;
  dsr.h1.datatype=NIFTI_DT_SIGNED_INT; // For the template
  dsr.h1.bitpix=32;
  dsr.h1.slice_start=0;
  for(int i=0; i<8; i++) dsr.h1.pixdim[i]=0.0;
  // https://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/qsform.html
  dsr.h1.pixdim[0]=1.0; // Set to either 1.0 or -1.0
  dsr.h1.pixdim[1]=1.0; // pixel size in x dimension
  dsr.h1.pixdim[2]=1.0; // pixel size in y dimension
  dsr.h1.pixdim[3]=1.0; // pixel size in z dimension
  dsr.h1.vox_offset=352; // Would be 0 for 1D format
  dsr.h1.scl_slope=1.0; // no need to scale pixel values
  dsr.h1.scl_inter=0.0; // no need to scale pixel values
  dsr.h1.slice_end=0;
  dsr.h1.slice_code=0;
  dsr.h1.xyzt_units=NIFTI_UNITS_MM+NIFTI_UNITS_SEC;
  dsr.h1.cal_max=tac.tacNr; // For the template
  dsr.h1.cal_min=0.0; // there may be voxels that have no TAC defined
  dsr.h1.slice_duration=0.0;
  dsr.h1.toffset=0.0;
  dsr.h1.glmax=dsr.h1.cal_max; // unused in NIfTI
  dsr.h1.glmin=0; // unused in NIfTI
  strlcpy(dsr.h1.descrip, "tac2nii", 80);
  strcpy(dsr.h1.aux_file, "");
  dsr.h1.qform_code=0;
  dsr.h1.sform_code=0;
  dsr.h1.quatern_b=0;
  dsr.h1.quatern_c=0;
  dsr.h1.quatern_d=0;
  dsr.h1.qoffset_x=0;
  dsr.h1.qoffset_y=0;
  dsr.h1.qoffset_z=0;
  for(int i=0; i<4; i++) dsr.h1.srow_x[i]=0;
  for(int i=0; i<4; i++) dsr.h1.srow_y[i]=0;
  for(int i=0; i<4; i++) dsr.h1.srow_z[i]=0;
  strcpy(dsr.h1.intent_name, "");
  strcpy(dsr.h1.magic, "n+1");  // Would be "ni1" for 1D format
  /* Extension is left as 0 0 0 0 */
 
 
  /*
   *  Write template, if required.
   */
  if(tname[0]) {
 
    if(verbose>1) printf("Make NIfTI file names for the template\n");
    char hdrfile[FILENAME_MAX], imgfile[FILENAME_MAX];
    if(niftiCreateFNames(tname, hdrfile, imgfile, NULL, IMG_FORMAT_NIFTI_1S)) {
      fprintf(stderr, "  Error: invalid NIfTI name %s\n", tname);
      tacFree(&tac); free(idata);
      return(11);
    }
 
    if(verbose>1) printf("Writing template NIfTI header\n");
    /* Delete previous NIfTI */
    if(fileExist(hdrfile)) remove(hdrfile);
    if(fileExist(imgfile)) remove(imgfile);
    /* Write NIfTI header */
    if(niftiWriteHeader(hdrfile, &dsr, verbose-1)) {
      fprintf(stderr, "Error: cannot write template header.\n");
      tacFree(&tac); free(idata);
      return(12);
    }
 
    if(verbose>1) printf("Writing NIfTI image data\n");
    FILE *fp=fopen(imgfile, "r+b");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open %s for write.\n", imgfile);
      tacFree(&tac); free(idata);
      if(fileExist(hdrfile)) remove(hdrfile);
      if(fileExist(imgfile)) remove(imgfile);
      return(13);
    }
    /* Move file pointer to the place of matrix data start */
    if(fseeko(fp, (size_t)dsr.h1.vox_offset, SEEK_SET)!=0) {
      fprintf(stderr, "Error: invalid file write position.\n");
      fclose(fp); tacFree(&tac); free(idata);
      if(fileExist(hdrfile)) remove(hdrfile);
      if(fileExist(imgfile)) remove(imgfile);
      return(14);
    }
    /* Write template data */
    if(fwrite(idata, sizeof(int), pxlNr, fp) != pxlNr) {
      fprintf(stderr, "Error: cannot write template matrix.\n");
      fclose(fp); tacFree(&tac); free(idata);
      if(fileExist(hdrfile)) remove(hdrfile);
      if(fileExist(imgfile)) remove(imgfile);
      return(15);
    }
    fclose(fp);
    if(verbose>0) printf("written %s\n", imgfile);
  }
 
  /*
   *  Edit NIfTI header for the 4D image
   */
  if(verbose>1) printf("Set NIfTI header\n");
  dsr.h1.dim[0]=4;
  dsr.h1.dim[4]=tac.sampleNr;
  dsr.h1.datatype=NIFTI_DT_FLOAT; // data as floats, so no need to scale
  dsr.h1.bitpix=32;
  dsr.h1.cal_max=max;
  if(min>0.0) dsr.h1.cal_min=0.0; // For voxels that do not represent any TAC
  else dsr.h1.cal_min=min;
  dsr.h1.glmax=max; // unused in NIfTI
  dsr.h1.glmin=dsr.h1.cal_min; // unused in NIfTI
  strlcpy(dsr.h1.descrip, "tac2nii", 80);
 
  /* Make NIfTI filenames */
  if(verbose>1) printf("Make NIfTI file names\n");
  char hdrfile[FILENAME_MAX], imgfile[FILENAME_MAX], siffile[FILENAME_MAX];
  if(niftiCreateFNames(dbname, hdrfile, imgfile, siffile, IMG_FORMAT_NIFTI_1S)) {
    fprintf(stderr, "  Error: invalid NIfTI name %s\n", dbname);
    tacFree(&tac); free(idata);
    return(21);
  }
 
  /* Allocate memory for float data for one frame */
  if(verbose>1) printf("Allocate memory for one xyz matrix\n");
  float *fdata;
  fdata=(float*)calloc(pxlNr, sizeof(float));
  if(fdata==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    tacFree(&tac); free(idata);
    return(6);
  }
 
  /*
   *  Write NIfTI header
   */
  if(verbose>1) printf("Writing NIfTI header\n");
  /* Delete previous NIfTI */
  /* It does not need to be valid NIfTI format, just that the filenames match */
  if(fileExist(hdrfile)) remove(hdrfile);
  if(fileExist(imgfile)) remove(imgfile);
  if(fileExist(siffile)) remove(siffile);
  /* Write NIfTI header */
  {
    if(niftiWriteHeader(hdrfile, &dsr, verbose-1)) {
      fprintf(stderr, "Error: cannot write header.\n");
      tacFree(&tac); free(fdata); free(idata);
      return(22);
    }
  }
 
  /*
   *  Write matrix data
   */
  if(verbose>1) printf("Writing NIfTI image data\n");
  FILE *fp=fopen(imgfile, "r+b");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open %s for write.\n", imgfile);
    tacFree(&tac); free(fdata); free(idata);
    if(fileExist(hdrfile)) remove(hdrfile);
    if(fileExist(imgfile)) remove(imgfile);
    return(23);
  }
  /* Move file pointer to the place of matrix data start */
  if(fseeko(fp, (size_t)dsr.h1.vox_offset, SEEK_SET)!=0) {
    fprintf(stderr, "Error: invalid file write position.\n");
    fclose(fp); tacFree(&tac); free(fdata); free(idata);
    if(fileExist(hdrfile)) remove(hdrfile);
    if(fileExist(imgfile)) remove(imgfile);
    return(24);
  }
  /* Write frame data */
  for(int fi=0; fi<tac.sampleNr; fi++) {
    if(verbose>8) printf("Writing frame %d\n", 1+fi);
    /* Fill matrix data */
    memset(fdata, 0, pxlNr*sizeof(float));
    for(size_t i=0; i<pxlNr; i++) {
      if(idata[i]==0) fdata[i]=0.0;
      else fdata[i]=tac.c[idata[i]-1].y[fi];
    }
    if(fwrite(fdata, sizeof(float), pxlNr, fp) != pxlNr) {
      fprintf(stderr, "Error: cannot write image matrix.\n");
      fclose(fp); tacFree(&tac); free(fdata); free(idata);
      if(fileExist(hdrfile)) remove(hdrfile);
      if(fileExist(imgfile)) remove(imgfile);
      return(25);
    }
  }
  fclose(fp); free(fdata); free(idata);
  if(verbose>0) printf("written %s\n", imgfile);
 
  /*
   *  Write SIF
   */
  if(verbose>1) printf("Making SIF\n");
  fp=fopen(siffile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot write %s\n", siffile);
    tacFree(&tac);
    return(31);
  }
  tac.tacNr=0;
  if(tacWriteSIF(&tac, fp, 0, NULL)) { 
    fprintf(stderr, "Error: cannot write %s\n", siffile);
    tacFree(&tac); fclose(fp);
    return(31);
  }
  fclose(fp);
  if(verbose>0) printf("written %s\n", siffile);
 
  tacFree(&tac);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
/*****************************************************************************/
 
/*****************************************************************************/
int simRCDims(
  const long long N,
  const int mxdim,
  const int mydim,
  const int mzdim,
  int *rnx,
  int *rny,
  int *rnz,
  int *rcxdim,
  int *rcydim,
  int *rczdim
) {
  if(rnx!=NULL) *rnx=0;
  if(rny!=NULL) *rny=0;
  if(rnz!=NULL) *rnz=0;
  if(rcxdim!=NULL) *rcxdim=0;
  if(rcydim!=NULL) *rcydim=0;
  if(rczdim!=NULL) *rczdim=0;
  if(N<1 || mxdim<1 || mydim<1 || mzdim<1) return(1);
 
  int xd, yd, zd;
  xd=mxdim; yd=mydim; zd=mzdim;
 
  long long nx=1, ny=1, nz=1;
  long long maxn=nx*ny*nz;
  long long n=N-maxn;
  while(n>0 && (nz<mzdim || ny<mydim || nx<mxdim)) {
//printf("  n=%lld nx=%lld ny=%lld nz=%lld\n", n, nx, ny, nz);
    if(n>0 && nx<mxdim) {
      nx++; maxn=nx*ny*nz; n=N-maxn;
    }
    if(n>0 && ny<mydim) {
      ny++; maxn=nx*ny*nz; n=N-maxn;
    }
    if(n>0 && nz<mzdim) {
      nz++; maxn=nx*ny*nz; n=N-maxn;
    }
  }
  if(n>0) return(2);
 
  xd/=nx; yd/=ny; zd/=nz;
 
  if(rnx!=NULL) *rnx=nx;
  if(rny!=NULL) *rny=ny;
  if(rnz!=NULL) *rnz=nz;
  if(rcxdim!=NULL) *rcxdim=xd;
  if(rcydim!=NULL) *rcydim=yd;
  if(rczdim!=NULL) *rczdim=zd;
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/