dcmxform.c

/*****************************************************************************/
#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 XFORM matrix from DICOM.",
  " ",
  "NOT for production use!",
  " ",
  "Usage: @P [-Options] dicomfile",
  " ",
  "Options:",
  " -xform=<Y|n>",
  "     List (y, default) or do not list (n) the xform matrix.",
  " -quatern=<y|N>",
  "     List (y) or do not list (n, default) the NIfTI quatern parameters.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: dcmlhdr, dcmframe, dcmmlist",
  " ",
  "Keywords: image, DICOM",
  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];
  int   listXform=1;   // 0=no, 1=yes
  int   listQuatern=0; // 0=no, 1=yes
  int   ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dcmfile[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;
    if(strncasecmp(cptr, "XFORM=", 6)==0) {
      if((listXform=tpcYesNo(cptr+6))>=0) continue;
    } else if(strncasecmp(cptr, "QUATERN=", 8)==0) {
      if((listQuatern=tpcYesNo(cptr+8))>=0) 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) {
    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);
    printf("listXform := %d\n", listXform);
    printf("listQuatern := %d\n", listQuatern);
  }
 
  /*
   *  Read file
   */
  DCMFILE dcm; dcmfileInit(&dcm);
  ret=dcmFileRead(dcmfile, &dcm, 1, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    dcmfileFree(&dcm); return(2);
  }
  if(verbose>3) {printf("max_tree_depth := %u\n", dcmfileMaxDepth(&dcm)); fflush(stdout);}
 
  /* ... and necessary DICOM parameters */
  if(verbose>1) {printf("reading Image Orientation (Patient)\n"); fflush(stdout);}
  double iop[6];
  ret=dcmImgOrient(&dcm, iop, verbose-3);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot read Image Orientation (Patient).\n");
    dcmfileFree(&dcm); return(3);
  }
  if(verbose>1) {printf("reading Image Position (Patient)\n"); fflush(stdout);}
  double ipp[3];
  ret=dcmImgPos(&dcm, ipp, verbose-3);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot read Image Position (Patient).\n");
    dcmfileFree(&dcm); return(4);
  }
  if(verbose>1) {printf("reading Voxel sizes\n"); fflush(stdout);}
  double pxlsize[3];
  ret=dcmImgPxlsize(&dcm, pxlsize, verbose-3);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot find voxel size.\n");
    dcmfileFree(&dcm); return(5);
  }
  dcmfileFree(&dcm);
  if(verbose>2) {
    printf("imgorient := %g", iop[0]); for(int i=1; i<6; i++) printf(", %g", iop[i]);
    printf("\n");
    printf("imgpos := %g", ipp[0]); for(int i=1; i<3; i++) printf(", %g", ipp[i]);
    printf("\n");
    printf("pxlsize := %g", pxlsize[0]); for(int i=1; i<3; i++) printf(", %g", pxlsize[i]);
    printf("\n");
  }
 
  /*
   *  Calculate XFORM matrix
   */
  double xform[16];
  ret=dcmImgXform(iop, pxlsize, ipp, xform, verbose-1);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot calculate XFORM.\n");
    if(verbose>1) printf("dcmXform() := %d\n", ret);
    return(7);
  }
  /* and print it, if requested */
  if(listXform!=0) {
    for(int i=0; i<4; i++) {
      printf("%g", xform[4*i]);
      for(int j=1; j<4; j++) printf("\t%g", xform[j+4*i]);
      printf("\n");
    }
  }
 
  /*
   *  Calculate Quatern parameters for NIfTI, if requested
   */
  if(listQuatern!=0) {
    double quatern[6];
    ret=dcmXformToQuatern(xform, quatern+0, quatern+3, verbose-3);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot calculate quatern parameters.\n");
      if(verbose>1) printf("dcmXformToQuatern() := %d\n", ret);
      return(8);
    }
 
    printf("quatern_b := %g\n", quatern[0]);
    printf("quatern_c := %g\n", quatern[1]);
    printf("quatern_d := %g\n", quatern[2]);
    printf("qoffset_x := %g\n", quatern[3]);
    printf("qoffset_y := %g\n", quatern[4]);
    printf("qoffset_z := %g\n", quatern[5]);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/