simelli.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"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Create a 3D PET image file in NIfTI 1S format, with pixels values of 1",
  "inside of a general ellipsoid and 0 outside of it,",
  "for simple simulations and SW testing.",
  "Matrix size and ellipsoid diameters must be given in pixels.",
  " ",
  "Usage: @P [Options] xdim ydim zdim xdiam ydiam zdiam imagefile",
  " ",
  "Options:",
  " -cx=<pixels> || -cy=<pixels> || -cz=<pixels>",
  "    Move the centre of ellipsoid from the image centre",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: dft2img, tac2nii, 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   dbname[FILENAME_MAX];
  int    dimx=0, dimy=0, dimz=0;
  double diamx=0.0, diamy=0.0, diamz=0.0;
  double cx=0.0, cy=0.0, cz=0.0;
  int    ret=0;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dbname[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') { /* options */
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(strncasecmp(cptr, "CX=", 3)==0) {
      cx=atofVerified(cptr+3); if(!isnan(cx)) continue;
    } else if(strncasecmp(cptr, "CY=", 3)==0) {
      cy=atofVerified(cptr+3); if(!isnan(cy)) continue;
    } else if(strncasecmp(cptr, "CZ=", 3)==0) {
      cz=atofVerified(cptr+3); if(!isnan(cz)) 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) {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) diamx=atofVerified(argv[ai++]);
  if(ai<argc) diamy=atofVerified(argv[ai++]);
  if(ai<argc) diamz=atofVerified(argv[ai++]);
  if(isnan(diamx) || isnan(diamy) || isnan(diamz)) {
    fprintf(stderr, "Error: invalid diameter.\n"); return(1);}
  if(ai<argc) {strlcpy(dbname, 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("dbname := %s\n", dbname);
    printf("dimx := %d\n", dimx);
    printf("dimy := %d\n", dimy);
    printf("dimz := %d\n", dimz);
    printf("diamx := %g\n", diamx);
    printf("diamy := %g\n", diamy);
    printf("diamz := %g\n", diamz);
    printf("cx := %g\n", cx);
    printf("cy := %g\n", cy);
    printf("cz := %g\n", cz);
    fflush(stdout);
  }
 
  /* Move centre to the middle of the matrix */
  cx+=0.5*(double)dimx+0.5;
  cy+=0.5*(double)dimy+0.5;
  cz+=0.5*(double)dimz+0.5;
 
 
  /*
   *  Make image data
   */
  if(verbose>1) printf("Allocate memory for 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"); return(4);}
  if(verbose>1) printf("Fill the template\n");
  for(int z=0; z<dimz; z++) for(int y=0; y<dimy; y++) for(int x=0; x<dimx; x++) {
    size_t pos=((size_t)dimy*dimx)*z + dimx*y + x;
    double dz=(2.0*((double)z-cz)/diamz);
    double dy=(2.0*((double)y-cy)/diamy);
    double dx=(2.0*((double)x-cx)/diamx);
    if(dz*dz+dy*dy+dx*dx<=1.0) idata[pos]=1; else idata[pos]=0;
  }
 
 
  /* 
   *  Set NIfTI header contents
   */
  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; // just one frame
  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;
  dsr.h1.cal_max=1.0; // This image has only zeroes and ones
  dsr.h1.cal_min=0.0; // This image has only zeroes and ones
  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, "simelli", 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 image.
   */
  if(verbose>1) printf("Make NIfTI file names\n");
  char hdrfile[FILENAME_MAX], imgfile[FILENAME_MAX];
  if(niftiCreateFNames(dbname, hdrfile, imgfile, NULL, IMG_FORMAT_NIFTI_1S)) {
    fprintf(stderr, "  Error: invalid NIfTI name %s\n", dbname);
    free(idata); return(11);
  }
 
  if(verbose>1) printf("Writing 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");
    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);
    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); free(idata);
    if(fileExist(hdrfile)) remove(hdrfile);
    if(fileExist(imgfile)) remove(imgfile);
    return(14);
  }
  /* Write pixel data */
  if(fwrite(idata, sizeof(int), pxlNr, fp) != pxlNr) {
    fprintf(stderr, "Error: cannot write template matrix.\n");
    fclose(fp); free(idata);
    if(fileExist(hdrfile)) remove(hdrfile);
    if(fileExist(imgfile)) remove(imgfile);
    return(15);
  }
  fclose(fp); free(idata);
  if(verbose>0) printf("written %s\n", imgfile);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/