Attachment 'mriio.c'

Download

   1 /*
   2  *       FILE NAME:   mriio.c
   3  *
   4  *       DESCRIPTION: utilities for reading/writing MRI data structure
   5  *
   6  *       AUTHOR:      Bruce Fischl
   7  *       DATE:        4/12/97
   8  *
   9  */
  10 
  11 /*-----------------------------------------------------
  12   INCLUDE FILES
  13   -------------------------------------------------------*/
  14 #define USE_ELECTRIC_FENCE 1
  15 #define _MRIIO_SRC
  16 
  17 #include <stdio.h>
  18 #include <stdlib.h>
  19 #include <string.h>
  20 #include <math.h>
  21 #include <string.h>
  22 #include <ctype.h>
  23 #include <unistd.h>
  24 #include <memory.h>
  25 #include <sys/stat.h>
  26 #include <sys/types.h>
  27 #include <ctype.h>
  28 #include <dirent.h>
  29 #include <time.h>
  30 #include <errno.h>
  31 #include <fcntl.h>
  32 
  33 #include "utils.h"
  34 #include "error.h"
  35 #include "proto.h"
  36 #include "mri.h"
  37 #include "macros.h"
  38 #include "diag.h"
  39 #include "volume_io.h"
  40 #include "region.h"
  41 #include "machine.h"
  42 #include "analyze.h"
  43 #include "fio.h"
  44 #include "mri_identify.h"
  45 #include "signa.h"
  46 #include "fio.h"
  47 #include "matfile.h"
  48 #include "math.h"
  49 #include "matrix.h"
  50 #include "diag.h"
  51 #include "chklc.h"
  52 #include "mriColorLookupTable.h"
  53 #include "DICOMRead.h"
  54 #include "imautils.h"
  55 #include "cma.h"
  56 #include "Bruker.h"
  57 #include "bfileio.h"
  58 #include "AFNI.h"
  59 #include "mghendian.h"
  60 #include "tags.h"
  61 #ifdef Darwin
  62 // /usr/include/zconf.h should typedef Byte, but doesnt on the Mac:
  63 typedef unsigned char  Byte;  /* 8 bits */
  64 #endif
  65 #include "nifti1.h"
  66 #include "znzlib.h"
  67 
  68 // unix director separator
  69 #define DIR_SEPARATOR '/'
  70 #define CURDIR "./"
  71 
  72 #define MM_PER_METER  1000.0f
  73 #define INFO_FNAME    "COR-.info"
  74 
  75 #ifdef Linux
  76 extern void swab(const void *from, void *to, size_t n);
  77 #endif
  78 
  79 #if 0
  80 static int NormalizeVector(float *v, int n);
  81 static MRI *mincRead2(char *fname, int read_volume);
  82 static int mincWrite2(MRI *mri, char *fname);
  83 static int GetMINCInfo(MRI *mri,
  84                        char *dim_names[4],
  85                        int   dim_sizes[4],
  86                        Real separations[4],
  87                        Real dircos[3][3],
  88                        Real VolCenterVox[3],
  89                        Real VolCenterWorld[3]);
  90 #endif
  91 
  92 static MRI *mri_read
  93 (char *fname, int type, int volume_flag, int start_frame, int end_frame);
  94 static MRI *corRead(char *fname, int read_volume);
  95 static int corWrite(MRI *mri, char *fname);
  96 static MRI *siemensRead(char *fname, int read_volume);
  97 static MRI *readGCA(char *fname) ;
  98 
  99 static MRI *mincRead(char *fname, int read_volume);
 100 static int mincWrite(MRI *mri, char *fname);
 101 static int bvolumeWrite(MRI *vol, char *fname_passed, int type);
 102 //static int bshortWrite(MRI *mri, char *fname_passed);
 103 //static int bfloatWrite(MRI *mri, char *stem);
 104 static int write_bhdr(MRI *mri, FILE *fp);
 105 static int read_bhdr(MRI *mri, FILE *fp);
 106 
 107 static MRI *bvolumeRead(char *fname_passed, int read_volume, int type);
 108 //static MRI *bshortRead(char *fname_passed, int read_volume);
 109 //static MRI *bfloatRead(char *fname_passed, int read_volume);
 110 static MRI *genesisRead(char *stem, int read_volume);
 111 static MRI *gelxRead(char *stem, int read_volume);
 112 
 113 static int CountAnalyzeFiles(char *analyzefname, int nzpad, char **ppstem);
 114 static MRI *analyzeRead(char *fname, int read_volume);
 115 static dsr *ReadAnalyzeHeader(char *hdrfile, int *swap,
 116                               int *mritype, int *bytes_per_voxel);
 117 static int DumpAnalyzeHeader(FILE *fp, dsr *hdr);
 118 
 119 
 120 static int analyzeWrite(MRI *mri, char *fname);
 121 static int analyzeWriteFrame(MRI *mri, char *fname, int frame);
 122 static int analyzeWriteSeries(MRI *mri, char *fname);
 123 static int analyzeWrite4D(MRI *mri, char *fname);
 124 
 125 static void swap_analyze_header(dsr *hdr);
 126 
 127 #if 0
 128 static int orient_with_register(MRI *mri);
 129 #endif
 130 static int nan_inf_check(MRI *mri);
 131 #ifdef VC_TO_CV
 132 static int voxel_center_to_center_voxel
 133 (MRI *mri, float *x, float *y, float *z);
 134 #endif
 135 static MRI *gdfRead(char *fname, int read_volume);
 136 static int gdfWrite(MRI *mri, char *fname);
 137 
 138 static MRI *ximgRead(char *fname, int read_volume);
 139 
 140 static MRI *nifti1Read(char *fname, int read_volume);
 141 static int nifti1Write(MRI *mri, char *fname);
 142 static MRI *niiRead(char *fname, int read_volume);
 143 static int niiWrite(MRI *mri, char *fname);
 144 static int niftiQformToMri(MRI *mri, struct nifti_1_header *hdr);
 145 static int mriToNiftiQform(MRI *mri, struct nifti_1_header *hdr);
 146 static int niftiSformToMri(MRI *mri, struct nifti_1_header *hdr);
 147 static void swap_nifti_1_header(struct nifti_1_header *hdr);
 148 static MRI *MRISreadCurvAsMRI(char *curvfile, int read_volume);
 149 
 150 /********************************************/
 151 
 152 static void short_local_buffer_to_image
 153 (short *buf, MRI *mri, int slice, int frame) ;
 154 
 155 static void int_local_buffer_to_image
 156 (int *buf, MRI *mri, int slice, int frame) ;
 157 
 158 static void long32_local_buffer_to_image
 159 (long32 *buf, MRI *mri, int slice, int frame) ;
 160 
 161 static void float_local_buffer_to_image
 162 (float *buf, MRI *mri, int slice, int frame) ;
 163 
 164 static void local_buffer_to_image
 165 (BUFTYPE *buf, MRI *mri, int slice, int frame) ;
 166 
 167 static MRI *sdtRead(char *fname, int read_volume);
 168 static MRI *mghRead(char *fname, int read_volume, int frame) ;
 169 static int mghWrite(MRI *mri, char *fname, int frame) ;
 170 static int mghAppend(MRI *mri, char *fname, int frame) ;
 171 
 172 /********************************************/
 173 
 174 extern int errno;
 175 
 176 extern char *Progname;
 177 /*char *Progname;*/
 178 
 179 static char *command_line;
 180 static char *subject_name;
 181 
 182 #define MAX_UNKNOWN_LABELS  100
 183 
 184 static short cma_field[512][512];
 185 static char unknown_labels[MAX_UNKNOWN_LABELS][STRLEN];
 186 static int n_unknown_labels;
 187 
 188 //////////////////////////////////////////////////////////////////////
 189 // this is a one way of setting direction cosine
 190 // when the direction cosine is not provided in the volume.
 191 // may not agree with the volume. what can we do?  Let them set by themselves.
 192 int setDirectionCosine(MRI *mri, int orientation)
 193 {
 194   switch(orientation)
 195     {
 196     case MRI_CORONAL: // x is from right to left.
 197                       // y is from top to neck, z is from back to front
 198       mri->x_r = -1; mri->y_r =  0; mri->z_r =  0; mri->c_r = 0;
 199       mri->x_a =  0; mri->y_a =  0; mri->z_a =  1; mri->c_a = 0;
 200       mri->x_s =  0; mri->y_s = -1; mri->z_s =  0; mri->c_s = 0;
 201       break;
 202     case MRI_SAGITTAL: // x is frp, back to front,
 203                        // y is from top to neck, z is from left to right
 204       mri->x_r =  0; mri->y_r =  0; mri->z_r =  1; mri->c_r = 0;
 205       mri->x_a =  1; mri->y_a =  0; mri->z_a =  0; mri->c_a = 0;
 206       mri->x_s =  0; mri->y_s = -1; mri->z_s =  0; mri->c_s = 0;
 207       break;
 208     case MRI_HORIZONTAL: // x is from right to left,
 209                          // y is from front to back, z is from neck to top
 210       mri->x_r = -1; mri->y_r =  0; mri->z_r =  0; mri->c_r = 0;
 211       mri->x_a =  0; mri->y_a = -1; mri->z_a =  0; mri->c_a = 0;
 212       mri->x_s =  0; mri->y_s =  0; mri->z_s =  1; mri->c_s = 0;
 213       break;
 214     default:
 215       ErrorReturn
 216         (ERROR_BADPARM,
 217          (ERROR_BADPARM, "setDirectionCosine():unknown slice direction"));
 218       break; // should not reach here (handled at the conversion)
 219     }
 220   return (NO_ERROR);
 221 }
 222 
 223 #define isOne(a)  (fabs(fabs(a)-1)<0.00001)
 224 
 225 // here I take the narrow view of slice_direction
 226 int getSliceDirection(MRI *mri)
 227 {
 228   int direction = MRI_UNDEFINED;
 229 
 230   if (isOne(mri->x_r) && isOne(mri->y_s) && isOne(mri->z_a))
 231     direction = MRI_CORONAL;
 232   else if (isOne(mri->x_a) && isOne(mri->y_s) && isOne(mri->z_r))
 233     direction = MRI_SAGITTAL;
 234   else if (isOne(mri->x_r) && isOne(mri->y_a) && isOne( mri->z_s))
 235     direction = MRI_HORIZONTAL;
 236   return direction;
 237 }
 238 
 239 // For surface, we currently cannot handle volumes with general slice direction
 240 // nor we cannot handle non-conformed volumes
 241 int mriOKforSurface(MRI *mri)
 242 {
 243   // first check slice direction
 244   if (getSliceDirection(mri) != MRI_CORONAL)
 245     return 0;
 246   // remove slice size limitation
 247   // else if (mri->width != 256 || mri->height != 256 || mri->depth != 256)
 248   //  return 0;
 249   // remove check for 1 mm size
 250   //  else if (mri->xsize != 1 || mri->ysize != 1 || mri->zsize != 1)
 251   //  return 0;
 252   else
 253     return 1;
 254 }
 255 
 256 int mriConformed(MRI *mri)
 257 {
 258   // first check slice direction
 259   if (getSliceDirection(mri) != MRI_CORONAL)
 260     return 0;
 261   else if (mri->width != 256 || mri->height != 256 || mri->depth != 256)
 262     return 0;
 263   else if (mri->xsize != 1 || mri->ysize != 1 || mri->zsize != 1)
 264     return 0;
 265   else
 266     return 1;
 267 }
 268 
 269 void setMRIforSurface(MRI *mri)
 270 {
 271   if (!mriOKforSurface(mri))
 272     ErrorExit(ERROR_BADPARM,
 273               "%s: the volume is not conformed, that is, "
 274               "the volume must be in CORONAL direction.\n", Progname) ;
 275 #if 0
 276   else
 277     {
 278       // we checked conformed in mriOKforSurface().
 279       // The only thing missing is c_(r,a,s) = 0
 280       // for surface creation assume that the
 281       // volume is conformed and c_(r,a,s) = 0
 282       mri->c_r=mri->c_a=mri->c_s = 0;
 283       if (mri->i_to_r__)
 284         MatrixFree(&mri->i_to_r__) ;
 285       if (mri->r_to_i__)
 286         MatrixFree(&mri->r_to_i__) ;
 287       mri->i_to_r__ = extract_i_to_r(mri);
 288       mri->r_to_i__ = extract_r_to_i(mri);
 289     }
 290 #endif
 291 }
 292 
 293 int mriio_command_line(int argc, char *argv[])
 294 {
 295 
 296   int i;
 297   int length;
 298   char *c;
 299 
 300   length = 0;
 301   for(i = 0;i < argc;i++)
 302     length += strlen(argv[i]);
 303 
 304   /* --- space for spaces and \0 --- */
 305   length += argc;
 306 
 307   command_line = (char *)malloc(length);
 308 
 309   c = command_line;
 310   for(i = 0;i < argc;i++)
 311     {
 312       strcpy(c, argv[i]);
 313       c += strlen(argv[i]);
 314       *c = (i == argc-1 ? '\0' : ' ');
 315       c++;
 316     }
 317 
 318   return(NO_ERROR);
 319 
 320 } /* end mriio_command_line() */
 321 
 322 int mriio_set_subject_name(char *name)
 323 {
 324 
 325   if(subject_name == NULL)
 326     subject_name = (char *)malloc(STRLEN);
 327 
 328   if(subject_name == NULL)
 329     {
 330       errno = 0;
 331       ErrorReturn(ERROR_NO_MEMORY, (ERROR_NO_MEMORY,
 332                                     "mriio_set_subject_name(): "
 333                                     "could't allocate %d bytes...!", STRLEN));
 334     }
 335 
 336   if(name == NULL)
 337     strcpy(subject_name, name);
 338   else
 339     {
 340       free(subject_name);
 341       subject_name = NULL;
 342     }
 343 
 344   return(NO_ERROR);
 345 
 346 } /* end mriio_set_subject_name() */
 347 
 348 int MRIgetVolumeName(char *string, char *name_only)
 349 {
 350 
 351   char *at, *pound;
 352 
 353   strcpy(name_only, string);
 354 
 355   if((at = strrchr(name_only, '@')) != NULL)
 356     *at = '\0';
 357 
 358   if((pound = strrchr(name_only, '#')) != NULL)
 359     *pound = '\0';
 360 
 361   return(NO_ERROR);
 362 
 363 } /* end MRIgetVolumeName() */
 364 
 365 static MRI *mri_read
 366 (char *fname, int type, int volume_flag, int start_frame, int end_frame)
 367 {
 368   MRI *mri, *mri2;
 369   IMAGE *I;
 370   char fname_copy[STRLEN];
 371   char *ptmpstr;
 372   char *at, *pound, *colon;
 373   char *ep;
 374   int i, j, k, t;
 375   int volume_frames;
 376 
 377   // sanity-checks
 378   if(fname == NULL)
 379     {
 380       errno = 0;
 381       ErrorReturn(NULL, (ERROR_BADPARM, "mri_read(): null fname!\n"));
 382     }
 383   if(fname[0] == 0)
 384     {
 385       errno = 0;
 386       ErrorReturn(NULL, (ERROR_BADPARM, "mri_read(): empty fname!\n"));
 387     }
 388 
 389   // if filename does not contain any directory separator, then add cwd
 390   if (!strchr(fname,DIR_SEPARATOR))
 391     {
 392       char *cwd = getcwd(NULL, 0); // posix 1 extension
 393                                    // (allocate as much space needed)
 394       if (cwd)
 395         {
 396           strcpy(fname_copy, cwd);
 397           strcat(fname_copy, "/");
 398           strcat(fname_copy, fname);
 399           free(cwd);
 400         }
 401       else // why fail?
 402         strcpy(fname_copy, fname);
 403     }
 404   else
 405     strcpy(fname_copy, fname);
 406 
 407   at = strrchr(fname_copy, '@');
 408   pound = strrchr(fname_copy, '#');
 409 
 410   if(at != NULL)
 411     {
 412       *at = '\0';
 413       at++;
 414     }
 415 
 416   if(pound != NULL)
 417     {
 418       *pound = '\0';
 419       pound++;
 420     }
 421 
 422   if(at != NULL)
 423     {
 424       type = string_to_type(at);
 425       if(type == MRI_VOLUME_TYPE_UNKNOWN)
 426         {
 427           errno = 0;
 428           ErrorReturn
 429             (NULL, (ERROR_BADPARM, "mri_read(): unknown type '%s'\n", at));
 430         }
 431     }
 432   else if(type == MRI_VOLUME_TYPE_UNKNOWN)
 433     {
 434       type = mri_identify(fname_copy);
 435       if(type == MRI_VOLUME_TYPE_UNKNOWN)
 436         {
 437           errno = 0;
 438           ErrorReturn
 439             (NULL,
 440              (ERROR_BADPARM,
 441               "mri_read(): couldn't determine type of file %s", fname_copy));
 442         }
 443     }
 444 
 445   if(pound != NULL)
 446     {
 447       colon = strchr(pound, ':');
 448       if(colon != NULL)
 449         {
 450           *colon = '\0';
 451           colon++;
 452           if(*colon == '\0')
 453             {
 454               errno = 0;
 455               ErrorReturn
 456                 (NULL,
 457                  (ERROR_BADPARM,
 458                   "mri_read(): bad frame specification ('%s:')\n", pound));
 459             }
 460 
 461           start_frame = strtol(pound, &ep, 10);
 462           if(*ep != '\0')
 463             {
 464               errno = 0;
 465               ErrorReturn
 466                 (NULL,
 467                  (ERROR_BADPARM,
 468                   "mri_read(): bad start frame ('%s')\n", pound));
 469             }
 470 
 471           end_frame = strtol(colon, &ep, 10);
 472           if(*ep != '\0')
 473             {
 474               errno = 0;
 475               ErrorReturn
 476                 (NULL,
 477                  (ERROR_BADPARM,
 478                   "mri_read(): bad end frame ('%s')\n", colon));
 479             }
 480 
 481         }
 482       else
 483         {
 484           start_frame = end_frame = strtol(pound, &ep, 10);
 485           if(*ep != '\0')
 486             {
 487               errno = 0;
 488               ErrorReturn
 489                 (NULL,
 490                  (ERROR_BADPARM,
 491                   "mri_read(): bad frame ('%s')\n", pound));
 492             }
 493         }
 494 
 495       if(start_frame < 0)
 496         {
 497           errno = 0;
 498           ErrorReturn
 499             (NULL,
 500              (ERROR_BADPARM,
 501               "mri_read(): frame (%d) is less than zero\n", start_frame));
 502         }
 503 
 504       if(end_frame < 0)
 505         {
 506           errno = 0;
 507           ErrorReturn
 508             (NULL,
 509              (ERROR_BADPARM,
 510               "mri_read(): frame (%d) is less than zero\n", end_frame));
 511         }
 512 
 513       if(start_frame > end_frame)
 514         {
 515           errno = 0;
 516           ErrorReturn(NULL, (ERROR_BADPARM,
 517                              "mri_read(): the start frame (%d) is "
 518                              "greater than the end frame (%d)\n",
 519                              start_frame, end_frame));
 520         }
 521 
 522     }
 523 
 524   if(type == MRI_CORONAL_SLICE_DIRECTORY)
 525     {
 526       mri = corRead(fname_copy, volume_flag);
 527     }
 528   else if(type == SIEMENS_FILE)
 529     {
 530       mri = siemensRead(fname_copy, volume_flag);
 531     }
 532   else if (type == MRI_GCA_FILE)
 533     {
 534       mri = readGCA(fname_copy) ;
 535     }
 536   else if(type == BHDR)
 537     {
 538       ptmpstr = bhdr_firstslicefname(fname_copy);
 539       t = bhdr_precision(fname_copy);
 540       mri = bvolumeRead(ptmpstr, volume_flag, t);
 541       free(ptmpstr);
 542     }
 543   else if(type == BSHORT_FILE)
 544     {
 545       //mri = bshortRead(fname_copy, volume_flag);
 546       mri = bvolumeRead(fname_copy, volume_flag, MRI_SHORT);
 547     }
 548   else if(type == BFLOAT_FILE)
 549     {
 550       //mri = bfloatRead(fname_copy, volume_flag);
 551       mri = bvolumeRead(fname_copy, volume_flag, MRI_FLOAT);
 552     }
 553   else if(type == GENESIS_FILE)
 554     {
 555       mri = genesisRead(fname_copy, volume_flag);
 556     }
 557   else if(type == SIGNA_FILE)
 558     {
 559       mri = signaRead(fname_copy, volume_flag);
 560     }
 561   else if(type == GE_LX_FILE)
 562     {
 563       mri = gelxRead(fname_copy, volume_flag);
 564     }
 565   else if(type == MRI_ANALYZE_FILE || type == MRI_ANALYZE4D_FILE)
 566     {
 567       mri = analyzeRead(fname_copy, volume_flag);
 568     }
 569   else if(type == BRIK_FILE)
 570     {
 571       mri = afniRead(fname_copy, volume_flag);
 572     }
 573   else if(type == MRI_MINC_FILE)
 574     {
 575       //mri = mincRead2(fname_copy, volume_flag);
 576       mri = mincRead(fname_copy, volume_flag);
 577     }
 578   else if(type == SDT_FILE)
 579     {
 580       mri = sdtRead(fname_copy, volume_flag);
 581     }
 582   else if(type == MRI_MGH_FILE)
 583     {
 584       mri = mghRead(fname_copy, volume_flag, -1);
 585     }
 586   else if(type == GDF_FILE)
 587     {
 588       mri = gdfRead(fname_copy, volume_flag);
 589     }
 590   else if(type == DICOM_FILE)
 591     {
 592       DICOMRead(fname_copy, &mri, volume_flag);
 593     }
 594   else if(type == SIEMENS_DICOM_FILE)
 595     {
 596       // mri_convert -nth option sets start_frame = nth.  otherwise -1
 597       mri = sdcmLoadVolume(fname_copy, volume_flag, start_frame);
 598       start_frame = -1;
 599       // in order to avoid the later processing on start_frame and end_frame
 600       // read the comment later on
 601       end_frame = 0;
 602     }
 603   else if (type == BRUKER_FILE)
 604     {
 605       mri = brukerRead(fname_copy, volume_flag);
 606     }
 607   else if(type == XIMG_FILE)
 608     {
 609       mri = ximgRead(fname_copy, volume_flag);
 610     }
 611   else if(type == NIFTI1_FILE)
 612     {
 613       mri = nifti1Read(fname_copy, volume_flag);
 614     }
 615   else if(type == NII_FILE)
 616     mri = niiRead(fname_copy, volume_flag);
 617   else if(type == MRI_CURV_FILE)
 618     mri = MRISreadCurvAsMRI(fname_copy, volume_flag);
 619   else if (type == IMAGE_FILE) {
 620     I = ImageRead(fname_copy);
 621     mri = ImageToMRI(I);
 622     ImageFree(&I);
 623   }
 624   else
 625     {
 626       fprintf(stderr,"mri_read(): type = %d\n",type);
 627       errno = 0;
 628       ErrorReturn(NULL, (ERROR_BADPARM, "mri_read(): code inconsistency "
 629                          "(file type recognized but not caught)"));
 630     }
 631 
 632   if (mri == NULL)  return(NULL) ;
 633 
 634   // update/cache the transform
 635   if (mri->i_to_r__)
 636     MatrixFree(&(mri->i_to_r__));
 637   mri->i_to_r__ = extract_i_to_r(mri);
 638 
 639   if (mri->r_to_i__)
 640     MatrixFree(&(mri->r_to_i__));
 641   mri->r_to_i__ = extract_r_to_i(mri);
 642 
 643   /* Compute the FOV from the vox2ras matrix (don't rely on what
 644      may or may not be in the file).*/
 645 
 646   if(start_frame == -1) return(mri);
 647 
 648   /* --- select frames --- */
 649 
 650   if(start_frame >= mri->nframes)
 651     {
 652       volume_frames = mri->nframes;
 653       MRIfree(&mri);
 654       errno = 0;
 655       ErrorReturn
 656         (NULL,
 657          (ERROR_BADPARM,
 658           "mri_read(): start frame (%d) is "
 659           "out of range (%d frames in volume)",
 660           start_frame, volume_frames));
 661     }
 662 
 663   if(end_frame >= mri->nframes)
 664     {
 665       volume_frames = mri->nframes;
 666       MRIfree(&mri);
 667       errno = 0;
 668       ErrorReturn
 669         (NULL,
 670          (ERROR_BADPARM,
 671           "mri_read(): end frame (%d) is out of range (%d frames in volume)",
 672           end_frame, volume_frames));
 673     }
 674 
 675   if(!volume_flag)
 676     {
 677 
 678       if(nan_inf_check(mri) != NO_ERROR)
 679         {
 680           MRIfree(&mri);
 681           return(NULL);
 682         }
 683       mri2 = MRIcopy(mri, NULL);
 684       MRIfree(&mri);
 685       mri2->nframes = (end_frame - start_frame + 1);
 686       return(mri2);
 687 
 688     }
 689 
 690   if(start_frame == 0 && end_frame == mri->nframes-1)
 691     {
 692       if(nan_inf_check(mri) != NO_ERROR)
 693         {
 694           MRIfree(&mri);
 695           return(NULL);
 696         }
 697       return(mri);
 698     }
 699 
 700   /* reading the whole volume and then copying only
 701      some frames is a bit inelegant,
 702      but it's easier for now
 703      to do this than to insert the appropriate
 704      code into the read function for each format... */
 705   /////////////////////////////////////////////////////////////////
 706   // This method does not work, because
 707   // 1. each frame has different acq parameters
 708   // 2. when making frames, it takes the info only from the first frame
 709   // Thus taking a frame out must be done at each frame extraction
 710   ////////////////////////////////////////////////////////////////////
 711   mri2 = MRIallocSequence(mri->width,
 712                           mri->height,
 713                           mri->depth,
 714                           mri->type,
 715                           (end_frame - start_frame + 1));
 716   MRIcopyHeader(mri, mri2);
 717   mri2->nframes = (end_frame - start_frame + 1);
 718   mri2->imnr0 = 1;
 719   mri2->imnr0 = mri2->nframes;
 720 
 721   if(mri2->type == MRI_UCHAR)
 722     for(t = 0;t < mri2->nframes;t++)
 723       for(i = 0;i < mri2->width;i++)
 724         for(j = 0;j < mri2->height;j++)
 725           for(k = 0;k < mri2->depth;k++)
 726             MRIseq_vox(mri2, i, j, k, t) =
 727               MRIseq_vox(mri, i, j, k, t + start_frame);
 728 
 729   if(mri2->type == MRI_SHORT)
 730     for(t = 0;t < mri2->nframes;t++)
 731       for(i = 0;i < mri2->width;i++)
 732         for(j = 0;j < mri2->height;j++)
 733           for(k = 0;k < mri2->depth;k++)
 734             MRISseq_vox(mri2, i, j, k, t) =
 735               MRISseq_vox(mri, i, j, k, t + start_frame);
 736 
 737   if(mri2->type == MRI_INT)
 738     for(t = 0;t < mri2->nframes;t++)
 739       for(i = 0;i < mri2->width;i++)
 740         for(j = 0;j < mri2->height;j++)
 741           for(k = 0;k < mri2->depth;k++)
 742             MRIIseq_vox(mri2, i, j, k, t) =
 743               MRIIseq_vox(mri, i, j, k, t + start_frame);
 744 
 745   if(mri2->type == MRI_LONG)
 746     for(t = 0;t < mri2->nframes;t++)
 747       for(i = 0;i < mri2->width;i++)
 748         for(j = 0;j < mri2->height;j++)
 749           for(k = 0;k < mri2->depth;k++)
 750             MRILseq_vox(mri2, i, j, k, t) =
 751               MRILseq_vox(mri, i, j, k, t + start_frame);
 752 
 753   if(mri2->type == MRI_FLOAT)
 754     for(t = 0;t < mri2->nframes;t++)
 755       for(i = 0;i < mri2->width;i++)
 756         for(j = 0;j < mri2->height;j++)
 757           for(k = 0;k < mri2->depth;k++)
 758             MRIFseq_vox(mri2, i, j, k, t) =
 759               MRIFseq_vox(mri, i, j, k, t + start_frame);
 760 
 761   if(nan_inf_check(mri) != NO_ERROR)
 762     {
 763       MRIfree(&mri);
 764       return(NULL);
 765     }
 766 
 767   MRIfree(&mri);
 768 
 769   return(mri2);
 770 
 771 } /* end mri_read() */
 772 
 773 static int nan_inf_check(MRI *mri)
 774 {
 775 
 776   int i, j, k, t;
 777 
 778   if(mri->type != MRI_FLOAT)
 779     return(NO_ERROR);
 780 
 781   for(i = 0;i < mri->width;i++)
 782     for(j = 0;j < mri->height;j++)
 783       for(k = 0;k < mri->depth;k++)
 784         for(t = 0;t < mri->nframes;t++)
 785           if(!devFinite((MRIFseq_vox(mri, i, j, k, t))))
 786             {
 787               if(devIsinf((MRIFseq_vox(mri, i, j, k, t))) != 0)
 788                 {
 789                   errno = 0;
 790                   ErrorReturn
 791                     (ERROR_BADPARM,
 792                      (ERROR_BADPARM,
 793                       "nan_inf_check(): Inf at voxel %d, %d, %d, %d",
 794                       i, j, k, t));
 795                 }
 796               else if(devIsnan((MRIFseq_vox(mri, i, j, k, t))))
 797                 {
 798                   errno = 0;
 799                   ErrorReturn
 800                     (ERROR_BADPARM,
 801                      (ERROR_BADPARM,
 802                       "nan_inf_check(): NaN at voxel %d, %d, %d, %d",
 803                       i, j, k, t));
 804                 }
 805               else
 806                 {
 807                   errno = 0;
 808                   ErrorReturn
 809                     (ERROR_BADPARM,
 810                      (ERROR_BADPARM,
 811                       "nan_inf_check(): bizarre value (not Inf, "
 812                       "not NaN, but not finite) at %d, %d, %d, %d",
 813                       i, j, k, t));
 814                 }
 815             }
 816 
 817   return(NO_ERROR);
 818 
 819 } /* end nan_inf_check() */
 820 
 821 MRI *MRIreadType(char *fname, int type)
 822 {
 823 
 824   MRI *mri;
 825 
 826   chklc();
 827 
 828   mri = mri_read(fname, type, TRUE, -1, -1);
 829 
 830   return(mri);
 831 
 832 } /* end MRIreadType() */
 833 
 834 MRI *MRIread(char *fname)
 835 {
 836   char  buf[STRLEN] ;
 837   MRI *mri = NULL;
 838 
 839   chklc() ;
 840 
 841   FileNameFromWildcard(fname, buf) ; fname = buf ;
 842   mri = mri_read(fname, MRI_VOLUME_TYPE_UNKNOWN, TRUE, -1, -1);
 843 
 844   /* some volume format needs to read many
 845      different files for slices (GE DICOM or COR).
 846      we make sure that mri_read() read the slices, not just one   */
 847   if (mri==NULL)
 848     return NULL;
 849 
 850   MRIremoveNaNs(mri, mri) ;
 851   return(mri);
 852 
 853 } /* end MRIread() */
 854 
 855 // allow picking one frame out of many frame
 856 // currently implemented only for Siemens dicom file
 857 MRI *MRIreadEx(char *fname, int nthframe)
 858 {
 859   char  buf[STRLEN] ;
 860   MRI *mri = NULL;
 861 
 862   chklc() ;
 863 
 864   FileNameFromWildcard(fname, buf) ; fname = buf ;
 865   mri = mri_read(fname, MRI_VOLUME_TYPE_UNKNOWN, TRUE, nthframe, nthframe);
 866 
 867   /* some volume format needs to read many
 868      different files for slices (GE DICOM or COR).
 869      we make sure that mri_read() read the slices, not just one   */
 870   if (mri==NULL)
 871     return NULL;
 872 
 873   MRIremoveNaNs(mri, mri) ;
 874   return(mri);
 875 
 876 } /* end MRIread() */
 877 
 878 MRI *MRIreadInfo(char *fname)
 879 {
 880 
 881   MRI *mri = NULL;
 882 
 883   mri = mri_read(fname, MRI_VOLUME_TYPE_UNKNOWN, FALSE, -1, -1);
 884 
 885   return(mri);
 886 
 887 } /* end MRIreadInfo() */
 888 
 889 /*---------------------------------------------------------------
 890   MRIreadHeader() - reads the MRI header of the given file name.
 891   If type is MRI_VOLUME_TYPE_UNKNOWN, then the type will be
 892   inferred from the file name.
 893   ---------------------------------------------------------------*/
 894 MRI *MRIreadHeader(char *fname, int type)
 895 {
 896   int usetype;
 897   MRI *mri = NULL;
 898   char modFname[STRLEN];
 899   struct stat stat_buf;
 900 
 901   usetype = type;
 902 
 903   if (!strchr(fname,DIR_SEPARATOR))
 904     {
 905       char *cwd = getcwd(NULL, 0); // posix 1 extension
 906                                    // (allocate as much space needed)
 907       if (cwd)
 908         {
 909           strcpy(modFname, cwd);
 910           strcat(modFname, "/");
 911           strcat(modFname, fname);
 912           free(cwd);
 913         }
 914       else // why fail?
 915         strcpy(modFname, fname);
 916     }
 917   else
 918     strcpy(modFname, fname);
 919 
 920   if(usetype == MRI_VOLUME_TYPE_UNKNOWN){
 921     usetype = mri_identify(modFname);
 922     if(usetype == MRI_VOLUME_TYPE_UNKNOWN)
 923       {
 924         // just check again
 925         if (stat(fname, &stat_buf) < 0)
 926           printf("ERROR: cound not find volume %s.  Does it exist?\n", fname);
 927         else
 928           printf("ERROR: could not determine type of %s\n",fname);
 929         return(NULL);
 930       }
 931   }
 932   mri = mri_read(modFname, usetype, FALSE, -1, -1);
 933 
 934   return(mri);
 935 
 936 } /* end MRIreadInfo() */
 937 
 938 int MRIwriteType(MRI *mri, char *fname, int type)
 939 {
 940   struct stat stat_buf;
 941   int error=0;
 942   char *fstem;
 943   char tmpstr[STRLEN];
 944 
 945   if(type == MRI_CORONAL_SLICE_DIRECTORY)
 946     {
 947       error = corWrite(mri, fname);
 948     }
 949   else
 950     {
 951       /* ----- all remaining types should write to
 952          a filename, not to within a directory,
 953          so check that it isn't an existing directory name we've been passed.
 954          failure to check will result in a segfault
 955          when attempting to write file data
 956          to a 'file' that is actually an existing directory ----- */
 957       if(stat(fname, &stat_buf) == 0){  // if can stat, then fname exists
 958         if(S_ISDIR(stat_buf.st_mode)){ // if is directory...
 959           errno = 0;
 960           ErrorReturn
 961             (ERROR_BADFILE,
 962              (ERROR_BADFILE,
 963               "MRIwriteType(): %s is an existing directory. (type=%d)\n"
 964               "                A filename should be specified instead.",
 965               fname, type));
 966         }
 967       }
 968     }
 969 
 970   /* ----- continue file writing... ----- */
 971 
 972   if(type == MRI_MINC_FILE)
 973     {
 974       error = mincWrite(mri, fname);
 975     }
 976   else if(type == BHDR)
 977     {
 978       fstem = bhdr_stem(fname);
 979       if(mri->type == MRI_SHORT){
 980         sprintf(tmpstr,"%s_000.bshort",fstem);
 981         error = bvolumeWrite(mri, tmpstr, MRI_SHORT);
 982       }
 983       else{
 984         sprintf(tmpstr,"%s_000.bfloat",fstem);
 985         error = bvolumeWrite(mri, tmpstr, MRI_FLOAT);
 986       }
 987       free(fstem);
 988     }
 989   else if(type == BSHORT_FILE)
 990     {
 991       error = bvolumeWrite(mri, fname, MRI_SHORT);
 992     }
 993   else if(type == BFLOAT_FILE)
 994     {
 995       error = bvolumeWrite(mri, fname, MRI_FLOAT);
 996     }
 997   else if(type == MRI_ANALYZE_FILE)
 998     {
 999       error = analyzeWrite(mri, fname);
1000     }
1001   else if(type == MRI_ANALYZE4D_FILE)
1002     {
1003       error = analyzeWrite4D(mri, fname);
1004     }
1005   else if(type == BRIK_FILE)
1006     {
1007       error = afniWrite(mri, fname);
1008     }
1009   else if(type == MRI_MGH_FILE)
1010     {
1011       error = mghWrite(mri, fname, -1);
1012     }
1013   else if(type == GDF_FILE)
1014     {
1015       error = gdfWrite(mri, fname);
1016     }
1017   else if(type == NIFTI1_FILE)
1018     {
1019       error = nifti1Write(mri, fname);
1020     }
1021   else if(type == NII_FILE) error = niiWrite(mri, fname);
1022   else if(type == GENESIS_FILE)
1023     {
1024       errno = 0;
1025       ErrorReturn
1026         (ERROR_BADPARM,
1027          (ERROR_BADPARM,
1028           "MRIwriteType(): writing of GENESIS file type not supported"));
1029     }
1030   else if(type == GE_LX_FILE)
1031     {
1032       errno = 0;
1033       ErrorReturn
1034         (ERROR_BADPARM,
1035          (ERROR_BADPARM,
1036           "MRIwriteType(): writing of GE LX file type not supported"));
1037     }
1038   else if(type == SIEMENS_FILE)
1039     {
1040       errno = 0;
1041       ErrorReturn
1042         (ERROR_BADPARM,
1043          (ERROR_BADPARM,
1044           "MRIwriteType(): writing of SIEMENS file type not supported"));
1045     }
1046   else if(type == SDT_FILE)
1047     {
1048       errno = 0;
1049       ErrorReturn
1050         (ERROR_BADPARM,
1051          (ERROR_BADPARM,
1052           "MRIwriteType(): writing of SDT file type not supported"));
1053     }
1054   else if(type == OTL_FILE)
1055     {
1056       errno = 0;
1057       ErrorReturn
1058         (ERROR_BADPARM,
1059          (ERROR_BADPARM,
1060           "MRIwriteType(): writing of OTL file type not supported"));
1061     }
1062   else if(type == RAW_FILE)
1063     {
1064       errno = 0;
1065       ErrorReturn
1066         (ERROR_BADPARM,
1067          (ERROR_BADPARM,
1068           "MRIwriteType(): writing of RAW file type not supported"));
1069     }
1070   else if(type == SIGNA_FILE)
1071     {
1072       errno = 0;
1073       ErrorReturn
1074         (ERROR_BADPARM,
1075          (ERROR_BADPARM,
1076           "MRIwriteType(): writing of SIGNA file type not supported"));
1077     }
1078   else if(type == DICOM_FILE)
1079     {
1080       errno = 0;
1081       ErrorReturn
1082         (ERROR_BADPARM,
1083          (ERROR_BADPARM,
1084           "MRIwriteType(): writing of DICOM file type not supported"));
1085     }
1086   else if(type == SIEMENS_DICOM_FILE)
1087     {
1088       errno = 0;
1089       ErrorReturn
1090         (ERROR_BADPARM,
1091          (ERROR_BADPARM,
1092           "MRIwriteType(): writing of SIEMENS DICOM file type not supported"));
1093     }
1094   else if(type == BRUKER_FILE)
1095     {
1096       errno = 0;
1097       ErrorReturn
1098         (ERROR_BADPARM,
1099          (ERROR_BADPARM,
1100           "MRIwriteType(): writing of BRUKER file type not supported"));
1101     }
1102   else if(type == XIMG_FILE)
1103     {
1104       errno = 0;
1105       ErrorReturn
1106         (ERROR_BADPARM,
1107          (ERROR_BADPARM,
1108           "MRIwriteType(): writing of XIMG file type not supported"));
1109     }
1110   else if(type == MRI_CORONAL_SLICE_DIRECTORY)
1111     {
1112       // already processed above
1113     }
1114   else
1115     {
1116       errno = 0;
1117       ErrorReturn(ERROR_BADPARM,
1118                   (ERROR_BADPARM,
1119                    "MRIwriteType(): code inconsistency "
1120                    "(file type recognized but not caught)"));
1121     }
1122 
1123   return(error);
1124 
1125 } /* end MRIwriteType() */
1126 
1127 int
1128 MRIwriteFrame(MRI *mri, char *fname, int frame)
1129 {
1130   MRI *mri_tmp ;
1131 
1132   mri_tmp =  MRIcopyFrame(mri, NULL, frame, 0) ;
1133   MRIwrite(mri_tmp, fname) ;
1134   MRIfree(&mri_tmp) ;
1135   return(NO_ERROR) ;
1136 }
1137 
1138 int MRIwrite(MRI *mri, char *fname)
1139 {
1140 
1141   int int_type = -1;
1142   int error;
1143 
1144   chklc() ;
1145   if((int_type = mri_identify(fname)) < 0)
1146     {
1147       errno = 0;
1148       ErrorReturn
1149         (ERROR_BADPARM,
1150          (ERROR_BADPARM, "unknown file type for file (%s)", fname));
1151     }
1152 
1153   error = MRIwriteType(mri, fname, int_type);
1154 
1155   return(error);
1156 
1157 } /* end MRIwrite() */
1158 
1159 
1160 /* ----- required header fields ----- */
1161 
1162 #define COR_ALL_REQUIRED 0x00001fff
1163 
1164 #define IMNR0_FLAG   0x00000001
1165 #define IMNR1_FLAG   0x00000002
1166 #define PTYPE_FLAG   0x00000004
1167 #define X_FLAG       0x00000008
1168 #define Y_FLAG       0x00000010
1169 #define THICK_FLAG   0x00000020
1170 #define PSIZ_FLAG    0x00000040
1171 #define STRTX_FLAG   0x00000080
1172 #define ENDX_FLAG    0x00000100
1173 #define STRTY_FLAG   0x00000200
1174 #define ENDY_FLAG    0x00000400
1175 #define STRTZ_FLAG   0x00000800
1176 #define ENDZ_FLAG    0x00001000
1177 
1178 /* trivially time course clean */
1179 static MRI *corRead(char *fname, int read_volume)
1180 {
1181 
1182   MRI *mri;
1183   struct stat stat_buf;
1184   char fname_use[STRLEN];
1185   char *fbase;
1186   FILE *fp;
1187   int i, j;
1188   char line[STRLEN];
1189   int imnr0, imnr1, x, y, ptype;
1190   double fov, thick, psiz, locatn; /* using floats to read creates problems
1191                                       when checking values
1192                                       (e.g. thick = 0.00100000005) */
1193   float strtx, endx, strty, endy, strtz, endz;
1194   float tr, te, ti, flip_angle ;
1195   int ras_good_flag;
1196   float x_r, x_a, x_s;
1197   float y_r, y_a, y_s;
1198   float z_r, z_a, z_s;
1199   float c_r, c_a, c_s;
1200   char xform[STRLEN];
1201   long gotten;
1202   char xform_use[STRLEN];
1203   char *cur_char;
1204   char *last_slash;
1205 
1206   /* ----- check that it is a directory we've been passed ----- */
1207   strcpy(fname_use, fname);
1208   if(stat(fname_use, &stat_buf) < 0)
1209     {
1210       errno = 0;
1211       ErrorReturn(NULL,
1212                   (ERROR_BADFILE, "corRead(): can't stat %s", fname_use));
1213     }
1214 
1215   if(!S_ISDIR(stat_buf.st_mode))
1216     {
1217       /* remove the last bit and try again */
1218       cur_char = fname_use;
1219       last_slash = cur_char;
1220       while( *(cur_char+1) != '\0' )
1221         {
1222           if(*cur_char == '/')
1223             last_slash = cur_char;
1224           cur_char++;
1225         }
1226       *last_slash = '\0';
1227       if(stat(fname_use, &stat_buf) < 0)
1228         {
1229           errno = 0;
1230           ErrorReturn
1231             (NULL,
1232              (ERROR_BADFILE, "corRead(): can't stat %s", fname_use));
1233         }
1234       if(!S_ISDIR(stat_buf.st_mode))
1235         {
1236           errno = 0;
1237           ErrorReturn
1238             (NULL,
1239              (ERROR_BADFILE, "corRead(): %s isn't a directory", fname_use));
1240         }
1241     }
1242 
1243   /* ----- copy the directory name and remove any trailing '/' ----- */
1244   fbase = &fname_use[strlen(fname_use)];
1245   if(*(fbase-1) != '/')
1246     {
1247       *fbase = '/';
1248       fbase++;
1249     }
1250 
1251   /* ----- read the header file ----- */
1252   sprintf(fbase, "COR-.info");
1253   if((fp = fopen(fname_use, "r")) == NULL)
1254     {
1255       errno = 0;
1256       ErrorReturn
1257         (NULL, (ERROR_BADFILE, "corRead(): can't open file %s", fname_use));
1258     }
1259 
1260   /* ----- defaults (a good idea for non-required values...) ----- */
1261   xform[0] = '\0';
1262   ras_good_flag = 0;
1263   x_r = x_a = x_s = 0.0;
1264   y_r = y_a = y_s = 0.0;
1265   z_r = z_a = z_s = 0.0;
1266   c_r = c_a = c_s = 0.0;
1267   flip_angle = tr = te = ti = 0.0;
1268   fov = 0.0;
1269   locatn = 0.0;
1270 
1271   gotten = 0x00;
1272 
1273   while(fgets(line, STRLEN, fp) != NULL)
1274     {
1275       if(strncmp(line, "imnr0 ", 6) == 0)
1276         {
1277           sscanf(line, "%*s %d", &imnr0);
1278           gotten = gotten | IMNR0_FLAG;
1279         }
1280       else if(strncmp(line, "imnr1 ", 6) == 0)
1281         {
1282           sscanf(line, "%*s %d", &imnr1);
1283           gotten = gotten | IMNR1_FLAG;
1284         }
1285       else if(strncmp(line, "ptype ", 6) == 0)
1286         {
1287           sscanf(line, "%*s %d", &ptype);
1288           gotten = gotten | PTYPE_FLAG;
1289         }
1290       else if(strncmp(line, "x ", 2) == 0)
1291         {
1292           sscanf(line, "%*s %d", &x);
1293           gotten = gotten | X_FLAG;
1294         }
1295       else if(strncmp(line, "y ", 2) == 0)
1296         {
1297           sscanf(line, "%*s %d", &y);
1298           gotten = gotten | Y_FLAG;
1299         }
1300       else if(strncmp(line, "fov ", 4) == 0)
1301         {
1302           sscanf(line, "%*s %lf", &fov);
1303         }
1304       else if(strncmp(line, "thick ", 6) == 0)
1305         {
1306           sscanf(line, "%*s %lf", &thick);
1307           gotten = gotten | THICK_FLAG;
1308         }
1309       else if(strncmp(line, "flip ", 5) == 0)
1310         {
1311           sscanf(line+11, "%f", &flip_angle);
1312           flip_angle = RADIANS(flip_angle) ;
1313         }
1314       else if(strncmp(line, "psiz ", 5) == 0)
1315         {
1316           sscanf(line, "%*s %lf", &psiz);
1317           gotten = gotten | PSIZ_FLAG;
1318         }
1319       else if(strncmp(line, "locatn ", 7) == 0)
1320         {
1321           sscanf(line, "%*s %lf", &locatn);
1322         }
1323       else if(strncmp(line, "strtx ", 6) == 0)
1324         {
1325           sscanf(line, "%*s %f", &strtx);
1326           gotten = gotten | STRTX_FLAG;
1327         }
1328       else if(strncmp(line, "endx ", 5) == 0)
1329         {
1330           sscanf(line, "%*s %f", &endx);
1331           gotten = gotten | ENDX_FLAG;
1332         }
1333       else if(strncmp(line, "strty ", 6) == 0)
1334         {
1335           sscanf(line, "%*s %f", &strty);
1336           gotten = gotten | STRTY_FLAG;
1337         }
1338       else if(strncmp(line, "endy ", 5) == 0)
1339         {
1340           sscanf(line, "%*s %f", &endy);
1341           gotten = gotten | ENDY_FLAG;
1342         }
1343       else if(strncmp(line, "strtz ", 6) == 0)
1344         {
1345           sscanf(line, "%*s %f", &strtz);
1346           gotten = gotten | STRTZ_FLAG;
1347         }
1348       else if(strncmp(line, "endz ", 5) == 0)
1349         {
1350           sscanf(line, "%*s %f", &endz);
1351           gotten = gotten | ENDZ_FLAG;
1352         }
1353       else if(strncmp(line, "tr ", 3) == 0)
1354         {
1355           sscanf(line, "%*s %f", &tr);
1356         }
1357       else if(strncmp(line, "te ", 3) == 0)
1358         {
1359           sscanf(line, "%*s %f", &te);
1360         }
1361       else if(strncmp(line, "ti ", 3) == 0)
1362         {
1363           sscanf(line, "%*s %f", &ti);
1364         }
1365       else if(strncmp(line, "ras_good_flag ", 14) == 0)
1366         {
1367           sscanf(line, "%*s %d", &ras_good_flag);
1368         }
1369       else if(strncmp(line, "x_ras ", 6) == 0)
1370         {
1371           sscanf(line, "%*s %f %f %f", &x_r, &x_a, &x_s);
1372         }
1373       else if(strncmp(line, "y_ras ", 6) == 0)
1374         {
1375           sscanf(line, "%*s %f %f %f", &y_r, &y_a, &y_s);
1376         }
1377       else if(strncmp(line, "z_ras ", 6) == 0)
1378         {
1379           sscanf(line, "%*s %f %f %f", &z_r, &z_a, &z_s);
1380         }
1381       else if(strncmp(line, "c_ras ", 6) == 0)
1382         {
1383           sscanf(line, "%*s %f %f %f", &c_r, &c_a, &c_s);
1384         }
1385       else if(strncmp(line, "xform", 5) == 0 ||
1386               strncmp(line, "transform", 9) == 0)
1387         {
1388           sscanf(line, "%*s %s", xform);
1389         }
1390     }
1391 
1392   fclose(fp);
1393 
1394   /* ----- check for required fields ----- */
1395   if((gotten & COR_ALL_REQUIRED) != COR_ALL_REQUIRED)
1396     {
1397       errno = 0;
1398       ErrorPrintf(ERROR_BADFILE, "missing fields in file %s:", fname_use);
1399       if(!(gotten & IMNR0_FLAG))
1400         ErrorPrintf(ERROR_BADFILE, "  imnr0 field missing");
1401       if(!(gotten & IMNR1_FLAG))
1402         ErrorPrintf(ERROR_BADFILE, "  imnr1 field missing");
1403       if(!(gotten & PTYPE_FLAG))
1404         ErrorPrintf(ERROR_BADFILE, "  ptype field missing");
1405       if(!(gotten & X_FLAG))
1406         ErrorPrintf(ERROR_BADFILE, "  x field missing");
1407       if(!(gotten & Y_FLAG))
1408         ErrorPrintf(ERROR_BADFILE, "  y field missing");
1409       if(!(gotten & THICK_FLAG))
1410         ErrorPrintf(ERROR_BADFILE, "  thick field missing");
1411       if(!(gotten & PSIZ_FLAG))
1412         ErrorPrintf(ERROR_BADFILE, "  psiz field missing");
1413       if(!(gotten & STRTX_FLAG))
1414         ErrorPrintf(ERROR_BADFILE, "  strtx field missing");
1415       if(!(gotten & ENDX_FLAG))
1416         ErrorPrintf(ERROR_BADFILE, "  endx field missing");
1417       if(!(gotten & STRTY_FLAG))
1418         ErrorPrintf(ERROR_BADFILE, "  strty field missing");
1419       if(!(gotten & ENDY_FLAG))
1420         ErrorPrintf(ERROR_BADFILE, "  endy field missing");
1421       if(!(gotten & STRTZ_FLAG))
1422         ErrorPrintf(ERROR_BADFILE, "  strtz field missing");
1423       if(!(gotten & ENDZ_FLAG))
1424         ErrorPrintf(ERROR_BADFILE, "  endz field missing");
1425       return(NULL);
1426     }
1427 
1428   /* ----- check for required but forced (constant) values ----- */
1429 
1430   if(imnr0 != 1)
1431     {
1432       printf
1433         ("warning: non-standard value for imnr0 "
1434          "(%d, usually 1) in file %s\n",
1435          imnr0, fname_use);
1436     }
1437 
1438   if(imnr1 != 256)
1439     {
1440       printf
1441         ("warning: non-standard value for imnr1 "
1442          "(%d, usually 256) in file %s\n",
1443          imnr1, fname_use);
1444     }
1445 
1446   if(ptype != 2)
1447     {
1448       printf("warning: non-standard value for ptype "
1449              "(%d, usually 2) in file %s\n",
1450              ptype, fname_use);
1451     }
1452 
1453   if(x != 256)
1454     {
1455       printf("warning: non-standard value for x "
1456              "(%d, usually 256) in file %s\n",
1457              x, fname_use);
1458     }
1459 
1460   if(y != 256)
1461     {
1462       printf("warning: non-standard value for y "
1463              "(%d, usually 256) in file %s\n",
1464              y, fname_use);
1465     }
1466 
1467   if(thick != 0.001)
1468     {
1469       printf("warning: non-standard value for thick "
1470              "(%g, usually 0.001) in file %s\n",
1471              thick, fname_use);
1472     }
1473 
1474   if(psiz != 0.001)
1475     {
1476       printf("warning: non-standard value for psiz "
1477              "(%g, usually 0.001) in file %s\n",
1478              psiz, fname_use);
1479     }
1480 
1481   /* ----- copy header information to an mri structure ----- */
1482 
1483   if(read_volume)
1484     mri = MRIalloc(x, y, imnr1 - imnr0 + 1, MRI_UCHAR);
1485   else
1486     mri = MRIallocHeader(x, y, imnr1 - imnr0 + 1, MRI_UCHAR);
1487 
1488   /* hack */
1489   /*
1490     printf("%g, %g, %g\n", x_r, x_a, x_s);
1491     printf("%g, %g, %g\n", y_r, y_a, y_s);
1492     printf("%g, %g, %g\n", z_r, z_a, z_s);
1493   */
1494   if(x_r == 0.0 && x_a == 0.0 && x_s == 0.0 && y_r == 0.0 \
1495      && y_a == 0.0 && y_s == 0.0 && z_r == 0.0 && z_a == 0.0 && z_s == 0.0)
1496     {
1497       fprintf(stderr,
1498               "----------------------------------------------------\n"
1499               "Could not find the direction cosine information.\n"
1500               "Will use the CORONAL orientation.\n"
1501               "If not suitable, please provide the information "
1502               "in COR-.info file\n"
1503               "----------------------------------------------------\n"
1504               );
1505       x_r = -1.0;
1506       y_s = -1.0;
1507       z_a = 1.0;
1508       ras_good_flag = 0;
1509     }
1510 
1511   mri->imnr0 = imnr0;
1512   mri->imnr1 = imnr1;
1513   mri->fov = (float)(fov * 1000);
1514   mri->thick = (float)(thick * 1000);
1515   mri->ps = (float)(psiz * 1000);
1516   mri->xsize = mri->ps;
1517   mri->ysize = mri->ps;
1518   mri->zsize = (float)(mri->thick);
1519   mri->xstart = strtx * 1000;
1520   mri->xend = endx * 1000;
1521   mri->ystart = strty * 1000;
1522   mri->yend = endy * 1000;
1523   mri->zstart = strtz * 1000;
1524   mri->zend = endz * 1000;
1525   strcpy(mri->fname, fname);
1526   mri->tr = tr;
1527   mri->te = te;
1528   mri->ti = ti;
1529   mri->flip_angle = flip_angle ;
1530   mri->ras_good_flag = ras_good_flag;
1531   mri->x_r = x_r;  mri->x_a = x_a;  mri->x_s = x_s;
1532   mri->y_r = y_r;  mri->y_a = y_a;  mri->y_s = y_s;
1533   mri->z_r = z_r;  mri->z_a = z_a;  mri->z_s = z_s;
1534   mri->c_r = c_r;  mri->c_a = c_a;  mri->c_s = c_s;
1535 
1536   if(strlen(xform) > 0)
1537     {
1538 
1539       if(xform[0] == '/')
1540         strcpy(mri->transform_fname, xform);
1541       else
1542         sprintf(mri->transform_fname, "%s/%s", fname, xform);
1543 
1544       strcpy(xform_use, mri->transform_fname);
1545 
1546       if(!FileExists(xform_use))
1547         {
1548 
1549           sprintf(xform_use, "%s/../transforms/talairach.xfm", fname);
1550 
1551           if(!FileExists(xform_use))
1552             printf("can't find talairach file '%s'\n", xform_use);
1553 
1554         }
1555 
1556       if(FileExists(xform_use))
1557         {
1558           if(input_transform_file(xform_use, &(mri->transform)) == NO_ERROR)
1559             {
1560               mri->linear_transform =
1561                 get_linear_transform_ptr(&mri->transform);
1562               mri->inverse_linear_transform =
1563                 get_inverse_linear_transform_ptr(&mri->transform);
1564               mri->free_transform = 1;
1565               strcpy(mri->transform_fname, xform_use);
1566               if (DIAG_VERBOSE_ON)
1567                 fprintf
1568                   (stderr,
1569                    "INFO: loaded talairach xform : %s\n",
1570                    mri->transform_fname);
1571             }
1572           else
1573             {
1574               errno = 0;
1575               ErrorPrintf
1576                 (ERROR_BAD_FILE, "error loading transform from %s", xform_use);
1577               mri->linear_transform = NULL;
1578               mri->inverse_linear_transform = NULL;
1579               mri->free_transform = 1;
1580               (mri->transform_fname)[0] = '\0';
1581             }
1582         }
1583 
1584     }
1585 
1586   if(!read_volume)
1587     return(mri);
1588 
1589   /* ----- read the data files ----- */
1590   for(i = mri->imnr0;i <= imnr1;i++)
1591     {
1592       sprintf(fbase, "COR-%03d", i);
1593       if((fp = fopen(fname_use, "r")) == NULL)
1594         {
1595           MRIfree(&mri);
1596           errno = 0;
1597           ErrorReturn
1598             (NULL,
1599              (ERROR_BADFILE, "corRead(): can't open file %s", fname_use));
1600         }
1601       for(j = 0;j < mri->height;j++)
1602         {
1603           if(fread(mri->slices[i-mri->imnr0][j], 1, mri->width, fp) <
1604              mri->width)
1605             {
1606               MRIfree(&mri);
1607               errno = 0;
1608               ErrorReturn
1609                 (NULL,
1610                  (ERROR_BADFILE,
1611                   "corRead(): error reading from file %s", fname_use));
1612             }
1613         }
1614       fclose(fp);
1615     }
1616 
1617   return(mri);
1618 
1619 } /* end corRead() */
1620 
1621 static int corWrite(MRI *mri, char *fname)
1622 {
1623 
1624   struct stat stat_buf;
1625   char fname_use[STRLEN];
1626   char *fbase;
1627   FILE *fp;
1628   int i, j;
1629   int rv;
1630 
1631   /* ----- check the mri structure for COR file compliance ----- */
1632 
1633   if(mri->slices == NULL)
1634     {
1635       errno = 0;
1636       ErrorReturn(ERROR_BADPARM, (ERROR_BADPARM,
1637                                   "corWrite(): mri structure to be "
1638                                   "written contains no voxel data"));
1639     }
1640 
1641   if(mri->imnr0 != 1)
1642     {
1643       printf
1644         ("non-standard value for imnr0 (%d, usually 1) in volume structure\n",
1645          mri->imnr0);
1646     }
1647 
1648   if(mri->imnr1 != 256)
1649     {
1650       printf
1651         ("non-standard value for imnr1 (%d, "
1652          "usually 256) in volume structure\n",
1653          mri->imnr1);
1654     }
1655 
1656   if(mri->type != MRI_UCHAR)
1657     {
1658       printf("non-standard value for type "
1659              "(%d, usually %d) in volume structure\n",
1660              mri->type, MRI_UCHAR);
1661     }
1662 
1663   if(mri->width != 256)
1664     {
1665       printf("non-standard value for width "
1666              "(%d, usually 256) in volume structure\n",
1667              mri->width);
1668     }
1669 
1670   if(mri->height != 256)
1671     {
1672       printf("non-standard value for height "
1673              "(%d, usually 256) in volume structure\n",
1674              mri->height);
1675     }
1676 
1677   if(mri->thick != 1)
1678     {
1679       printf("non-standard value for thick "
1680              "(%g, usually 1) in volume structure\n",
1681              mri->thick);
1682     }
1683 
1684   if(mri->ps != 1){
1685     printf("non-standard value for ps "
1686            "(%g, usually 1) in volume structure\n",
1687            mri->ps);
1688   }
1689 
1690   /* ----- copy the directory name and remove any trailing '/' ----- */
1691   strcpy(fname_use, fname);
1692   fbase = &fname_use[strlen(fname_use)];
1693   if(*(fbase-1) != '/')
1694     {
1695       *fbase = '/';
1696       fbase++;
1697       *fbase = '\0';
1698     }
1699 
1700   /* Create the directory */
1701   rv = mkdir(fname_use,(mode_t)-1);
1702   if(rv != 0 && errno != EEXIST){
1703     printf("ERROR: creating directory %s\n",fname_use);
1704     perror(NULL);
1705     return(1);
1706   }
1707 
1708   /* ----- check that it is a directory we've been passed ----- */
1709   if(stat(fname, &stat_buf) < 0){
1710     errno = 0;
1711     ErrorReturn
1712       (ERROR_BADFILE, (ERROR_BADFILE, "corWrite(): can't stat %s", fname));
1713   }
1714 
1715   if(!S_ISDIR(stat_buf.st_mode)){
1716     errno = 0;
1717     ErrorReturn
1718       (ERROR_BADFILE,
1719        (ERROR_BADFILE, "corWrite(): %s isn't a directory", fname));
1720   }
1721 
1722   sprintf(fbase, "COR-.info");
1723   if((fp = fopen(fname_use, "w")) == NULL)
1724     {
1725       errno = 0;
1726       ErrorReturn
1727         (ERROR_BADFILE,
1728          (ERROR_BADFILE,
1729           "corWrite(): can't open file %s for writing", fname_use));
1730     }
1731 
1732   fprintf(fp, "imnr0 %d\n", mri->imnr0);
1733   fprintf(fp, "imnr1 %d\n", mri->imnr1);
1734   fprintf(fp, "ptype %d\n", 2);
1735   fprintf(fp, "x %d\n", mri->width);
1736   fprintf(fp, "y %d\n", mri->height);
1737   fprintf(fp, "fov %g\n", mri->fov / 1000.0);
1738   fprintf(fp, "thick %g\n", mri->zsize / 1000.0); /* was xsize 11/1/01 */
1739   fprintf(fp, "psiz %g\n", mri->xsize / 1000.0);  /* was zsize 11/1/01 */
1740   fprintf(fp, "locatn %g\n", 0.0);
1741   fprintf(fp, "strtx %g\n", mri->xstart / 1000.0);
1742   fprintf(fp, "endx %g\n", mri->xend / 1000.0);
1743   fprintf(fp, "strty %g\n", mri->ystart / 1000.0);
1744   fprintf(fp, "endy %g\n", mri->yend / 1000.0);
1745   fprintf(fp, "strtz %g\n", mri->zstart / 1000.0);
1746   fprintf(fp, "endz %g\n", mri->zend / 1000.0);
1747   fprintf(fp, "tr %f\n", mri->tr);
1748   fprintf(fp, "te %f\n", mri->te);
1749   fprintf(fp, "ti %f\n", mri->ti);
1750   fprintf(fp, "flip angle %f\n", DEGREES(mri->flip_angle));
1751   fprintf(fp, "ras_good_flag %d\n", mri->ras_good_flag);
1752   fprintf(fp, "x_ras %f %f %f\n", mri->x_r, mri->x_a, mri->x_s);
1753   fprintf(fp, "y_ras %f %f %f\n", mri->y_r, mri->y_a, mri->y_s);
1754   fprintf(fp, "z_ras %f %f %f\n", mri->z_r, mri->z_a, mri->z_s);
1755   fprintf(fp, "c_ras %f %f %f\n", mri->c_r, mri->c_a, mri->c_s);
1756   if(strlen(mri->transform_fname) > 0)
1757     fprintf(fp, "xform %s\n", mri->transform_fname);
1758 
1759   fclose(fp);
1760 
1761   for(i = mri->imnr0;i <= mri->imnr1;i++)
1762     {
1763       sprintf(fbase, "COR-%03d", i);
1764       if((fp = fopen(fname_use, "w")) == NULL)
1765         {
1766           errno = 0;
1767           ErrorReturn
1768             (ERROR_BADFILE,
1769              (ERROR_BADFILE,
1770               "corWrite(): can't open file %s for writing", fname_use));
1771         }
1772       for(j = 0;j < mri->height;j++)
1773         {
1774           if(fwrite(mri->slices[i-mri->imnr0][j], 1, mri->width, fp) <
1775              mri->width)
1776             {
1777               fclose(fp);
1778               errno = 0;
1779               ErrorReturn
1780                 (ERROR_BADFILE,
1781                  (ERROR_BADFILE,
1782                   "corWrite(): error writing to file %s ", fname_use));
1783             }
1784         }
1785       fclose(fp);
1786     }
1787 
1788   return(NO_ERROR);
1789 
1790 } /* end corWrite() */
1791 
1792 static MRI *siemensRead(char *fname, int read_volume_flag)
1793 {
1794 
1795   int file_n, n_low, n_high;
1796   char fname_use[STRLEN];
1797   MRI *mri;
1798   FILE *fp;
1799   char *c, *c2;
1800   short rows, cols;
1801   int n_slices;
1802   double d, d2;
1803   double im_c_r, im_c_a, im_c_s;
1804   int i, j;
1805   int n_files, base_raw_matrix_size, number_of_averages;
1806   int mosaic_size;
1807   int mosaic;
1808   int mos_r, mos_c;
1809   char pulse_sequence_name[STRLEN], ps2[STRLEN];
1810   int n_t;
1811   int n_dangling_images, n_full_mosaics, mosaics_per_volume;
1812   int br, bc;
1813   MRI *mri_raw;
1814   int t, s;
1815   int slice_in_mosaic;
1816   int file;
1817   char ima[4];
1818   IMAFILEINFO *ifi;
1819 
1820   /* ----- stop compiler complaints ----- */
1821   mri = NULL;
1822   mosaic_size = 0;
1823 
1824 
1825   ifi = imaLoadFileInfo(fname);
1826   if(ifi == NULL){
1827     printf("ERROR: siemensRead(): %s\n",fname);
1828     return(NULL);
1829   }
1830 
1831   strcpy(fname_use, fname);
1832 
1833   /* ----- point to ".ima" ----- */
1834   c = fname_use + (strlen(fname_use) - 4);
1835 
1836   if(c < fname_use)
1837     {
1838       errno = 0;
1839       ErrorReturn
1840         (NULL,
1841          (ERROR_BADPARM,
1842           "siemensRead(): bad file name %s (must end in '.ima' or '.IMA')",
1843           fname_use));
1844     }
1845   if(strcmp(".ima", c) == 0)
1846     sprintf(ima, "ima");
1847   else if(strcmp(".IMA", c) == 0)
1848     sprintf(ima, "IMA");
1849   else
1850     {
1851       errno = 0;
1852       ErrorReturn
1853         (NULL,
1854          (ERROR_BADPARM,
1855           "siemensRead(): bad file name %s (must end in '.ima' or '.IMA')",
1856           fname_use));
1857     }
1858 
1859   c2 = c;
1860   for(c--;isdigit((int)(*c));c--);
1861   c++;
1862 
1863   /* ----- c now points to the first digit in the last number set
1864      (e.g. to the "5" in 123-4-567.ima) ----- */
1865 
1866   /* ----- count down and up from this file --
1867      max and min image number within the series ----- */
1868   *c2 = '\0';
1869   file_n = atol(c);
1870   *c2 = '.';
1871 
1872   if(!FileExists(fname_use))
1873     {
1874       errno = 0;
1875       ErrorReturn
1876         (NULL,
1877          (ERROR_BADFILE, "siemensRead(): file %s doesn't exist", fname_use));
1878     }
1879 
1880   /* --- get the low image number --- */
1881   n_low = file_n - 1;
1882   sprintf(c, "%d.%s", n_low, ima);
1883   while(FileExists(fname_use))
1884     {
1885       n_low--;
1886       sprintf(c, "%d.%s", n_low, ima);
1887     }
1888   n_low++;
1889 
1890   /* --- get the high image number --- */
1891   n_high = file_n + 1;
1892   sprintf(c, "%d.%s", n_high, ima);
1893   while(FileExists(fname_use))
1894     {
1895       n_high++;
1896       sprintf(c, "%d.%s", n_high, ima);
1897     }
1898   n_high--;
1899 
1900   n_files = n_high - n_low + 1;
1901 
1902   sprintf(c, "%d.%s", n_low, ima);
1903   if((fp = fopen(fname_use, "r")) == NULL)
1904     {
1905       errno = 0;
1906       ErrorReturn
1907         (NULL,
1908          (ERROR_BADFILE,
1909           "siemensRead(): can't open file %s (low = %d, this = %d, high = %d)",
1910           fname_use, n_low, file_n, n_high));
1911     }
1912 
1913   fseek(fp, 4994, SEEK_SET);
1914   fread(&rows, 2, 1, fp);
1915   rows = orderShortBytes(rows);
1916   fseek(fp, 4996, SEEK_SET);
1917   fread(&cols, 2, 1, fp);
1918   cols = orderShortBytes(cols);
1919   fseek(fp, 4004, SEEK_SET);
1920   fread(&n_slices, 4, 1, fp);
1921   n_slices = orderIntBytes(n_slices);
1922   if (n_slices==0)
1923     {
1924       errno = 0;
1925       ErrorReturn
1926         (NULL,
1927          (ERROR_BADFILE,
1928           "\n\nPlease try with the option '-it siemens_dicom'.\n "
1929           "The handling failed assuming the old siemens format.\n"))
1930         }
1931   fseek(fp, 2864, SEEK_SET);
1932   fread(&base_raw_matrix_size, 4, 1, fp);
1933   base_raw_matrix_size = orderIntBytes(base_raw_matrix_size);
1934   fseek(fp, 1584, SEEK_SET);
1935   fread(&number_of_averages, 4, 1, fp);
1936   number_of_averages = orderIntBytes(number_of_averages);
1937   memset(pulse_sequence_name, 0x00, STRLEN);
1938   fseek(fp, 3009, SEEK_SET);
1939   fread(&pulse_sequence_name, 1, 65, fp);
1940 
1941   /* --- scout --- */
1942   strcpy(ps2, pulse_sequence_name);
1943   StrLower(ps2);
1944   if(strstr(ps2, "scout") != NULL)
1945     {
1946       errno = 0;
1947       ErrorReturn
1948         (NULL,
1949          (ERROR_BADPARM,
1950           "siemensRead(): series appears to be a scout "
1951           "(sequence file name is %s)", pulse_sequence_name));
1952     }
1953 
1954   /* --- structural --- */
1955   if(n_slices == 1 && ! ifi->IsMosaic)
1956     {
1957       n_slices = n_files;
1958       n_t = 1;
1959       if(base_raw_matrix_size != rows || base_raw_matrix_size != cols)
1960         {
1961           errno = 0;
1962           ErrorReturn
1963             (NULL,
1964              (ERROR_BADPARM, "siemensRead(): bad file/base matrix sizes"));
1965         }
1966       mos_r = mos_c = 1;
1967       mosaic_size = 1;
1968     }
1969   else
1970     {
1971 
1972       if(rows % base_raw_matrix_size != 0)
1973         {
1974           errno = 0;
1975           ErrorReturn
1976             (NULL,
1977              (ERROR_BADPARM, "siemensRead(): file rows (%hd) not"
1978               " divisible by image rows (%d)", rows, base_raw_matrix_size));
1979         }
1980       if(cols % base_raw_matrix_size != 0)
1981         {
1982           errno = 0;
1983           ErrorReturn
1984             (NULL,
1985              (ERROR_BADPARM,
1986               "siemensRead(): file cols (%hd)"
1987               " not divisible by image cols (%d)",
1988               cols, base_raw_matrix_size));
1989         }
1990 
1991       mos_r = rows / base_raw_matrix_size;
1992       mos_c = cols / base_raw_matrix_size;
1993       mosaic_size = mos_r * mos_c;
1994 
1995       n_dangling_images = n_slices % mosaic_size;
1996       n_full_mosaics = (n_slices - n_dangling_images) / mosaic_size;
1997 
1998       mosaics_per_volume = n_full_mosaics + (n_dangling_images == 0 ? 0 : 1);
1999 
2000       if(n_files % mosaics_per_volume != 0)
2001         {
2002           errno = 0;
2003           ErrorReturn(NULL, (ERROR_BADPARM, "siemensRead(): files in volume "
2004                              "(%d) not divisible by mosaics per volume (%d)",
2005                              n_files, mosaics_per_volume));
2006         }
2007 
2008       n_t = n_files / mosaics_per_volume;
2009 
2010     }
2011 
2012   if(read_volume_flag)
2013     mri = MRIallocSequence
2014       (base_raw_matrix_size, base_raw_matrix_size, n_slices, MRI_SHORT, n_t);
2015   else
2016     {
2017       mri = MRIallocHeader
2018         (base_raw_matrix_size, base_raw_matrix_size, n_slices, MRI_SHORT);
2019       mri->nframes = n_t;
2020     }
2021 
2022   /* --- pixel sizes --- */
2023   /* --- mos_r and mos_c factors are strange, but they're there... --- */
2024   fseek(fp, 5000, SEEK_SET);
2025   fread(&d, 8, 1, fp);
2026   mri->xsize = mos_r * orderDoubleBytes(d);
2027   fseek(fp, 5008, SEEK_SET);
2028   fread(&d, 8, 1, fp);
2029   mri->ysize = mos_c * orderDoubleBytes(d);
2030 
2031   /* --- slice distance factor --- */
2032   fseek(fp, 4136, SEEK_SET);
2033   fread(&d, 8, 1, fp);
2034   d = orderDoubleBytes(d);
2035   if(d == -19222) /* undefined (Siemens code) -- I assume this to mean 0 */
2036     d = 0.0;
2037   /* --- slice thickness --- */
2038   fseek(fp, 1544, SEEK_SET);
2039   fread(&d2, 8, 1, fp);
2040   d2 = orderDoubleBytes(d2);
2041   /* --- distance between slices --- */
2042   mri->zsize = (1.0 + d) * d2;
2043 
2044   /* --- field of view (larger of height, width fov) --- */
2045   fseek(fp, 3744, SEEK_SET);
2046   fread(&d, 8, 1, fp);
2047   d = orderDoubleBytes(d);
2048   fseek(fp, 3752, SEEK_SET);
2049   fread(&d2, 8, 1, fp);
2050   d2 = orderDoubleBytes(d);
2051   mri->fov = (d > d2 ? d : d2);
2052 
2053   mri->thick = mri->zsize;
2054   mri->ps = mri->xsize;
2055 
2056   strcpy(mri->fname, fname);
2057 
2058   mri->location = 0.0;
2059 
2060   fseek(fp, 2112, SEEK_SET) ;
2061   mri->flip_angle = RADIANS(freadDouble(fp)) ;  /* was in degrees */
2062 
2063   fseek(fp, 1560, SEEK_SET);
2064   fread(&d, 8, 1, fp);
2065   mri->tr = orderDoubleBytes(d);
2066   fseek(fp, 1568, SEEK_SET);
2067   fread(&d, 8, 1, fp);
2068   mri->te = orderDoubleBytes(d);
2069   fseek(fp, 1576, SEEK_SET);
2070   fread(&d, 8, 1, fp);
2071   mri->ti = orderDoubleBytes(d);
2072 
2073   fseek(fp, 3792, SEEK_SET);
2074   fread(&d, 8, 1, fp);  mri->z_r = -orderDoubleBytes(d);
2075   fread(&d, 8, 1, fp);  mri->z_a =  orderDoubleBytes(d);
2076   fread(&d, 8, 1, fp);  mri->z_s = -orderDoubleBytes(d);
2077 
2078   fseek(fp, 3832, SEEK_SET);
2079   fread(&d, 8, 1, fp);  mri->x_r = -orderDoubleBytes(d);
2080   fread(&d, 8, 1, fp);  mri->x_a =  orderDoubleBytes(d);
2081   fread(&d, 8, 1, fp);  mri->x_s = -orderDoubleBytes(d);
2082 
2083   fseek(fp, 3856, SEEK_SET);
2084   fread(&d, 8, 1, fp);  mri->y_r = -orderDoubleBytes(d);
2085   fread(&d, 8, 1, fp);  mri->y_a =  orderDoubleBytes(d);
2086   fread(&d, 8, 1, fp);  mri->y_s = -orderDoubleBytes(d);
2087 
2088   fseek(fp, 3768, SEEK_SET);
2089   fread(&im_c_r, 8, 1, fp);  im_c_r = -orderDoubleBytes(im_c_r);
2090   fread(&im_c_a, 8, 1, fp);  im_c_a =  orderDoubleBytes(im_c_a);
2091   fread(&im_c_s, 8, 1, fp);  im_c_s = -orderDoubleBytes(im_c_s);
2092 
2093   mri->c_r = im_c_r - (mosaic_size - 1) * mri->z_r * mri->zsize + \
2094     ((mri->depth - 1.0) / 2.0) * mri->z_r * mri->zsize;
2095   mri->c_a = im_c_a - (mosaic_size - 1) * mri->z_a * mri->zsize + \
2096     ((mri->depth - 1.0) / 2.0) * mri->z_a * mri->zsize;
2097   mri->c_s = im_c_s - (mosaic_size - 1) * mri->z_s * mri->zsize + \
2098     ((mri->depth - 1.0) / 2.0) * mri->z_s * mri->zsize;
2099 
2100   mri->ras_good_flag = 1;
2101 
2102   fseek(fp, 3760, SEEK_SET);
2103   fread(&i, 4, 1, fp);
2104   i = orderIntBytes(i);
2105 
2106 #if 0
2107   if(i == 1 || i == 2)
2108     mri->slice_direction = MRI_HORIZONTAL;
2109   else if(i == 3 || i == 5)
2110     mri->slice_direction = MRI_CORONAL;
2111   else if(i == 4 || i == 6)
2112     mri->slice_direction = MRI_SAGITTAL;
2113   else
2114 #endif
2115     if (i < 1 || i > 6)
2116       {
2117         errno = 0;
2118         ErrorReturn
2119           (NULL,
2120            (ERROR_BADFILE,
2121             "siemensRead(): bad slice direction (%d) in file %s",
2122             i, fname_use));
2123       }
2124 
2125   mri->xstart = -mri->width * mri->xsize / 2.0;
2126   mri->xend = -mri->xstart;
2127   mri->ystart = -mri->height * mri->ysize / 2.0;
2128   mri->yend = -mri->ystart;
2129   mri->zstart = -mri->depth * mri->zsize / 2.0;
2130   mri->zend = -mri->zstart;
2131 
2132   /*
2133     printf("%d, %d; %d, %hd, %hd; %d\n", n_files, number_of_averages,
2134     base_raw_matrix_size, rows, cols,
2135     slices);
2136   */
2137   /*
2138     rows, cols, brms, mosaic i, j, mosaic size, n slices, n files, n_t
2139   */
2140   fclose(fp);
2141 
2142   mri->imnr0 = 1;
2143   mri->imnr1 = mri->depth;
2144   /*
2145     printf("%d, %d, %d, %d\n",
2146     mri->width, mri->height, mri->depth, mri->nframes);
2147   */
2148   if(read_volume_flag)
2149     {
2150 
2151       mri_raw = MRIalloc(rows, cols, n_files, MRI_SHORT);
2152 
2153       for(file_n = n_low;file_n <= n_high;file_n++)
2154         {
2155 
2156           sprintf(c, "%d.%s", file_n, ima);
2157           if((fp = fopen(fname_use, "r")) == NULL)
2158             {
2159               MRIfree(&mri);
2160               errno = 0;
2161               ErrorReturn(NULL, (ERROR_BADFILE,
2162                                  "siemensRead(): can't open file %s "
2163                                  "(low = %d, this = %d, high = %d)",
2164                                  fname_use, n_low, file_n, n_high));
2165             }
2166 
2167           fseek(fp, 6144, SEEK_SET);
2168 
2169           for(i = 0;i < rows;i++)
2170             {
2171               fread(&MRISvox(mri_raw, 0, i, file_n - n_low),
2172                     sizeof(short), cols, fp);
2173 #if (BYTE_ORDER == LITTLE_ENDIAN)
2174               swab(&MRISvox(mri_raw, 0, i, file_n - n_low), \
2175                    &MRISvox(mri_raw, 0, i, file_n - n_low),
2176                    sizeof(short) * cols);
2177 #endif
2178             }
2179 
2180           fclose(fp);
2181 
2182         }
2183 
2184       for(t = 0;t < mri->nframes;t++)
2185         {
2186           for(s = 0;s < mri->depth;s++)
2187             {
2188               slice_in_mosaic = s % mosaic_size;
2189               mosaic = (s - slice_in_mosaic) / mosaic_size;
2190               file = mri->nframes * mosaic + t;
2191               /*
2192                 printf("s, t = %d, %d; f, sm = %d, %d\n",
2193                 s, t, file, slice_in_mosaic);
2194               */
2195               bc = slice_in_mosaic % mos_r;
2196               br = (slice_in_mosaic - bc) / mos_r;
2197 
2198               for(i = 0;i < mri->width;i++)
2199                 {
2200                   for(j = 0;j < mri->height;j++)
2201                     {
2202                       MRISseq_vox(mri, i, j, s, t) = \
2203                         MRISvox
2204                         (mri_raw, mri->width *
2205                          bc + i, mri->height * br + j, file);
2206                     }
2207                 }
2208 
2209             }
2210         }
2211 
2212       MRIfree(&mri_raw);
2213 
2214     }
2215 
2216   return(mri);
2217 
2218 } /* end siemensRead() */
2219 /*-----------------------------------------------------------*/
2220 static MRI *mincRead(char *fname, int read_volume)
2221 {
2222   // Real wx, wy, wz;
2223   MRI *mri;
2224   Volume vol;
2225   VIO_Status status;
2226   char *dim_names[4];
2227   int dim_sizes[4];
2228   int ndims;
2229   int dtype;
2230   volume_input_struct input_info;
2231   Real separations[4];
2232   Real voxel[4];
2233   Real worldr, worlda, worlds;
2234   Real val;
2235   int i, j, k, t;
2236   float xfov, yfov, zfov;
2237   Real f;
2238   BOOLEAN sflag = TRUE ;
2239   Transform *pVox2WorldLin;
2240   General_transform *pVox2WorldGen;
2241 
2242   /* ----- read the volume ----- */
2243 
2244   /*   we no longer much around axes and thus */
2245   /*   it is safe to read in the standard way    */
2246   dim_names[0] = MIxspace;
2247   dim_names[1] = MIyspace;
2248   dim_names[2] = MIzspace;
2249   dim_names[3] = MItime;
2250 
2251 #if 0
2252   dim_names[0] = MIzspace;
2253   dim_names[1] = MIyspace;
2254   dim_names[2] = MIxspace;
2255   dim_names[3] = MItime;
2256 #endif
2257 
2258   if(!FileExists(fname))
2259     {
2260       errno = 0;
2261       ErrorReturn(NULL, (ERROR_BADFILE,
2262                          "mincRead(): can't find file %s", fname));
2263     }
2264 
2265   status = start_volume_input(fname, 0, dim_names, NC_UNSPECIFIED, 0, 0, 0,
2266                               TRUE, &vol, NULL, &input_info);
2267 
2268   if (Gdiag & DIAG_VERBOSE_ON && DIAG_SHOW){
2269     printf("status = %d\n", status);
2270     printf("n_dimensions = %d\n", get_volume_n_dimensions(vol));
2271     printf("nc_data_type = %d\n", vol->nc_data_type);
2272   }
2273 
2274   if(status != OK){
2275     errno = 0;
2276     ErrorReturn(NULL, (ERROR_BADPARM, "mincRead(): error reading "
2277                        "volume from file %s", fname));
2278   }
2279 
2280   /* ----- check the number of dimensions ----- */
2281   ndims = get_volume_n_dimensions(vol);
2282   if(ndims != 3 && ndims != 4){
2283     errno = 0;
2284     ErrorReturn(NULL, (ERROR_BADPARM, "mincRead(): %d dimensions "
2285                        "in file; expecting 3 or 4", ndims));
2286   }
2287 
2288   /* ----- get the dimension sizes ----- */
2289   get_volume_sizes(vol, dim_sizes);
2290 
2291   /* --- one time point if there are only three dimensions in the file --- */
2292   if(ndims == 3) dim_sizes[3] = 1;
2293 
2294   if ((Gdiag & DIAG_SHOW) && DIAG_VERBOSE_ON)
2295     printf("DimSizes: %d, %d, %d, %d, %d\n", ndims, dim_sizes[0],
2296            dim_sizes[1], dim_sizes[2], dim_sizes[3]);
2297   if ((Gdiag & DIAG_SHOW) && DIAG_VERBOSE_ON)
2298     printf("DataType: %d\n", vol->nc_data_type);
2299 
2300   dtype = get_volume_nc_data_type(vol, &sflag) ;
2301   dtype = orderIntBytes(vol->nc_data_type) ;
2302 
2303   /* ----- get the data type ----- */
2304   if(vol->nc_data_type == NC_BYTE || vol->nc_data_type == NC_CHAR)
2305     dtype = MRI_UCHAR;
2306   else if(vol->nc_data_type == NC_SHORT)
2307     dtype = MRI_SHORT;
2308   else if(vol->nc_data_type == NC_LONG)
2309     dtype = MRI_LONG;
2310   else if(vol->nc_data_type == NC_FLOAT || vol->nc_data_type == NC_DOUBLE)
2311     dtype = MRI_FLOAT;
2312   else
2313     {
2314       errno = 0;
2315       ErrorReturn(NULL, (ERROR_BADPARM, "mincRead(): bad data type "
2316                          "(%d) in input file %s", vol->nc_data_type, fname));
2317     }
2318 
2319   /* ----- allocate the mri structure ----- */
2320   if(read_volume)
2321     mri = MRIallocSequence(dim_sizes[0], dim_sizes[1], dim_sizes[2],
2322                            dtype, dim_sizes[3]);
2323   else{
2324     mri = MRIallocHeader(dim_sizes[0], dim_sizes[1], dim_sizes[2], dtype);
2325     mri->nframes = dim_sizes[3];
2326   }
2327 
2328   /* ----- set up the mri structure ----- */
2329   get_volume_separations(vol, separations);
2330   mri->xsize = fabs(separations[0]);
2331   mri->ysize = fabs(separations[1]);
2332   mri->zsize = fabs(separations[2]);
2333   mri->ps    = mri->xsize;
2334   mri->thick = mri->zsize;
2335 
2336   mri->x_r = vol->direction_cosines[0][0];
2337   mri->x_a = vol->direction_cosines[0][1];
2338   mri->x_s = vol->direction_cosines[0][2];
2339 
2340   mri->y_r = vol->direction_cosines[1][0];
2341   mri->y_a = vol->direction_cosines[1][1];
2342   mri->y_s = vol->direction_cosines[1][2];
2343 
2344   mri->z_r = vol->direction_cosines[2][0];
2345   mri->z_a = vol->direction_cosines[2][1];
2346   mri->z_s = vol->direction_cosines[2][2];
2347 
2348   if(separations[0] < 0){
2349     mri->x_r = -mri->x_r;  mri->x_a = -mri->x_a;  mri->x_s = -mri->x_s;
2350   }
2351   if(separations[1] < 0){
2352     mri->y_r = -mri->y_r;  mri->y_a = -mri->y_a;  mri->y_s = -mri->y_s;
2353   }
2354   if(separations[2] < 0){
2355     mri->z_r = -mri->z_r;  mri->z_a = -mri->z_a;  mri->z_s = -mri->z_s;
2356   }
2357 
2358   voxel[0] = (mri->width) / 2.0;
2359   voxel[1] = (mri->height) / 2.0;
2360   voxel[2] = (mri->depth) / 2.0;
2361   voxel[3] = 0.0;
2362   convert_voxel_to_world(vol, voxel, &worldr, &worlda, &worlds);
2363   mri->c_r = worldr;
2364   mri->c_a = worlda;
2365   mri->c_s = worlds;
2366 
2367   mri->ras_good_flag = 1;
2368 
2369   mri->xend = mri->xsize * mri->width / 2.0;   mri->xstart = -mri->xend;
2370   mri->yend = mri->ysize * mri->height / 2.0;  mri->ystart = -mri->yend;
2371   mri->zend = mri->zsize * mri->depth / 2.0;   mri->zstart = -mri->zend;
2372 
2373   xfov = mri->xend - mri->xstart;
2374   yfov = mri->yend - mri->ystart;
2375   zfov = mri->zend - mri->zstart;
2376 
2377   mri->fov = ( xfov > yfov ? (xfov > zfov ? xfov : zfov ) :
2378                (yfov > zfov ? yfov : zfov ) );
2379 
2380   strcpy(mri->fname, fname);
2381 
2382   //////////////////////////////////////////////////////////////////////////
2383   // test transform their way and our way:
2384   // MRIvoxelToWorld(mri, voxel[0], voxel[1], voxel[2], &wx, &wy, &wz);
2385   // printf("MNI calculated %.2f, %.2f, %.2f
2386   //     vs. MRIvoxelToWorld: %.2f, %.2f, %.2f\n",
2387   //     worldr, worlda, worlds, wx, wy, wz);
2388 
2389   /* ----- copy the data from the file to the mri structure ----- */
2390   if(read_volume){
2391 
2392     while(input_more_of_volume(vol, &input_info, &f));
2393 
2394     for(i = 0;i < mri->width;i++) {
2395       for(j = 0;j < mri->height;j++) {
2396         for(k = 0;k < mri->depth;k++) {
2397           for(t = 0;t < mri->nframes;t++) {
2398             val = get_volume_voxel_value(vol, i, j, k, t, 0);
2399             if(mri->type == MRI_UCHAR)
2400               MRIseq_vox(mri, i, j, k, t) = (unsigned char)val;
2401             if(mri->type == MRI_SHORT)
2402               MRISseq_vox(mri, i, j, k, t) = (short)val;
2403             if(mri->type == MRI_LONG)
2404               MRILseq_vox(mri, i, j, k, t) = (long)val;
2405             if(mri->type == MRI_FLOAT)
2406               MRIFseq_vox(mri, i, j, k, t) = (float)val;
2407           }
2408         }
2409       }
2410     }
2411   }
2412 
2413   pVox2WorldGen = get_voxel_to_world_transform(vol);
2414   pVox2WorldLin = get_linear_transform_ptr(pVox2WorldGen);
2415   if ((Gdiag & DIAG_SHOW) && DIAG_VERBOSE_ON)
2416     {
2417       printf("MINC Linear Transform\n");
2418       for(i=0;i<4;i++){
2419         for(j=0;j<4;j++) printf("%7.4f ",pVox2WorldLin->m[j][i]);
2420         printf("\n");
2421       }
2422     }
2423 
2424   delete_volume_input(&input_info);
2425   delete_volume(vol);
2426 
2427 
2428   return(mri);
2429 
2430 } /* end mincRead() */
2431 #if 0
2432 /*-----------------------------------------------------------*/
2433 static MRI *mincRead2(char *fname, int read_volume)
2434 {
2435 
2436   MRI *mri;
2437   Volume vol;
2438   VIO_Status status;
2439   char *dim_names[4];
2440   int dim_sizes[4];
2441   int ndims;
2442   int dtype;
2443   volume_input_struct input_info;
2444   Real separations[4];
2445   Real voxel[4];
2446   Real worldr, worlda, worlds;
2447   Real val;
2448   int i, j, k, t;
2449   float xfov, yfov, zfov;
2450   Real f;
2451   BOOLEAN sflag = TRUE ;
2452   MATRIX *T;
2453   Transform *pVox2WorldLin;
2454   General_transform *pVox2WorldGen;
2455 
2456   /* Make sure file exists */
2457   if(!FileExists(fname)){
2458     errno = 0;
2459     ErrorReturn(NULL, (ERROR_BADFILE,
2460                        "mincRead(): can't find file %s", fname));
2461   }
2462 
2463   /* Specify read-in order */
2464   dim_names[0] = MIxspace; /* cols */
2465   dim_names[1] = MIzspace; /* rows */
2466   dim_names[2] = MIyspace; /* slices */
2467   dim_names[3] = MItime;   /* time */
2468 
2469   /* Read the header info into vol. input_info needed for further reads */
2470   status = start_volume_input(fname, 0, dim_names, NC_UNSPECIFIED, 0, 0, 0,
2471                               TRUE, &vol, NULL, &input_info);
2472 
2473   if (Gdiag & DIAG_VERBOSE_ON && DIAG_SHOW){
2474     printf("status = %d\n", status);
2475     printf("n_dimensions = %d\n", get_volume_n_dimensions(vol));
2476     printf("nc_data_type = %d\n", vol->nc_data_type);
2477   }
2478 
2479   if(status != OK){
2480     errno = 0;
2481     ErrorReturn(NULL, (ERROR_BADPARM, "mincRead(): error reading "
2482                        "volume from file %s", fname));
2483   }
2484 
2485   /* ----- check the number of dimensions ----- */
2486   ndims = get_volume_n_dimensions(vol);
2487   if(ndims != 3 && ndims != 4){
2488     errno = 0;
2489     ErrorReturn(NULL, (ERROR_BADPARM, "mincRead(): %d dimensions "
2490                        "in file; expecting 3 or 4", ndims));
2491   }
2492 
2493   /* ----- get the dimension sizes ----- */
2494   get_volume_sizes(vol, dim_sizes);
2495 
2496   /* --- one time point if there are only three dimensions in the file --- */
2497   if(ndims == 3) dim_sizes[3] = 1;
2498 
2499   dtype = get_volume_nc_data_type(vol, &sflag) ;
2500   dtype = orderIntBytes(vol->nc_data_type) ;
2501   get_volume_separations(vol, separations);
2502 
2503   for(i=0;i<4;i++)
2504     printf("%d %s %3d  %7.3f %6.4f %6.4f %6.4f \n",
2505            i,dim_names[i],dim_sizes[i],separations[i],
2506            vol->direction_cosines[i][0],vol->direction_cosines[i][1],
2507            vol->direction_cosines[i][2]);
2508   if ((Gdiag & DIAG_SHOW) && DIAG_VERBOSE_ON)
2509     printf("DataType: %d\n", vol->nc_data_type);
2510 
2511   /* Translate data type to that of mri structure */
2512   switch(vol->nc_data_type){
2513   case NC_BYTE:   dtype = MRI_UCHAR; break;
2514   case NC_CHAR:   dtype = MRI_UCHAR; break;
2515   case NC_SHORT:  dtype = MRI_SHORT; break;
2516   case NC_LONG:   dtype = MRI_LONG;  break;
2517   case NC_FLOAT:  dtype = MRI_FLOAT; break;
2518   case NC_DOUBLE: dtype = MRI_FLOAT; break;
2519   default:
2520     errno = 0;
2521     ErrorReturn(NULL, (ERROR_BADPARM, "mincRead(): bad data type "
2522                        "(%d) in input file %s", vol->nc_data_type, fname));
2523     break;
2524   }
2525 
2526   /* ----- allocate the mri structure ----- */
2527   if(read_volume)
2528     mri = MRIallocSequence(dim_sizes[0], dim_sizes[1], dim_sizes[2],
2529                            dtype, dim_sizes[3]);
2530   else{
2531     mri = MRIallocHeader(dim_sizes[0], dim_sizes[1], dim_sizes[2], dtype);
2532     mri->nframes = dim_sizes[3];
2533   }
2534 
2535   /* ----- set up the mri structure ----- */
2536   get_volume_separations(vol, separations);
2537   mri->xsize = fabs(separations[0]); /* xsize = col   resolution */
2538   mri->ysize = fabs(separations[1]); /* ysize = row   resolution */
2539   mri->zsize = fabs(separations[2]); /* zsize = slice resolution */
2540   mri->ps    = mri->xsize;
2541   mri->thick = mri->zsize;
2542 
2543   /* column direction cosines */
2544   mri->x_r = vol->direction_cosines[0][0];
2545   mri->x_a = vol->direction_cosines[0][1];
2546   mri->x_s = vol->direction_cosines[0][2];
2547 
2548   /* row direction cosines */
2549   mri->y_r = vol->direction_cosines[1][0];
2550   mri->y_a = vol->direction_cosines[1][1];
2551   mri->y_s = vol->direction_cosines[1][2];
2552 
2553   /* slice direction cosines */
2554   mri->z_r = vol->direction_cosines[2][0];
2555   mri->z_a = vol->direction_cosines[2][1];
2556   mri->z_s = vol->direction_cosines[2][2];
2557 
2558   if(separations[0] < 0){
2559     mri->x_r = -mri->x_r;  mri->x_a = -mri->x_a;  mri->x_s = -mri->x_s;
2560   }
2561   if(separations[1] < 0){
2562     mri->y_r = -mri->y_r;  mri->y_a = -mri->y_a;  mri->y_s = -mri->y_s;
2563   }
2564   if(separations[2] < 0){
2565     mri->z_r = -mri->z_r;  mri->z_a = -mri->z_a;  mri->z_s = -mri->z_s;
2566   }
2567 
2568   /* Get center point */       // don't.  our convention is different
2569   voxel[0] = (mri->width)/2.;  // (mri->width  - 1) / 2.0;
2570   voxel[1] = (mri->height)/2.; // (mri->height - 1) / 2.0;
2571   voxel[2] = (mri->depth)/2.;  //(mri->depth  - 1) / 2.0;
2572   voxel[3] = 0.0;
2573   convert_voxel_to_world(vol, voxel, &worldr, &worlda, &worlds);
2574   mri->c_r = worldr;
2575   mri->c_a = worlda;
2576   mri->c_s = worlds;
2577   printf("Center Voxel: %7.3f %7.3f %7.3f\n",voxel[0],voxel[1],voxel[2]);
2578   printf("Center World: %7.3f %7.3f %7.3f\n",mri->c_r,mri->c_a,mri->c_s);
2579 
2580   mri->ras_good_flag = 1;
2581 
2582   mri->xend = mri->xsize * mri->width / 2.0; mri->xstart = -mri->xend;
2583   mri->yend = mri->ysize * mri->height/ 2.0; mri->ystart = -mri->yend;
2584   mri->zend = mri->zsize * mri->depth / 2.0; mri->zstart = -mri->zend;
2585 
2586   xfov = mri->xend - mri->xstart;
2587   yfov = mri->yend - mri->ystart;
2588   zfov = mri->zend - mri->zstart;
2589 
2590   mri->fov =
2591     (xfov > yfov ? (xfov > zfov ? xfov : zfov) : (yfov > zfov ? yfov : zfov));
2592 
2593   strcpy(mri->fname, fname);
2594 
2595   T = MRIxfmCRS2XYZ(mri,0);
2596   printf("Input Coordinate Transform (CRS2XYZ)-------\n");
2597   MatrixPrint(stdout,T);
2598   MatrixFree(&T);
2599   pVox2WorldGen = get_voxel_to_world_transform(vol);
2600   pVox2WorldLin = get_linear_transform_ptr(pVox2WorldGen);
2601   if ((Gdiag & DIAG_SHOW) && DIAG_VERBOSE_ON)
2602     {
2603       printf("MINC Linear Transform ----------------------\n");
2604       for(i=0;i<4;i++){
2605         for(j=0;j<4;j++) printf("%7.4f ",pVox2WorldLin->m[j][i]);
2606         printf("\n");
2607       }
2608       printf("-------------------------------------------\n");
2609     }
2610 
2611   /* ----- copy the data from the file to the mri structure ----- */
2612   if(read_volume){
2613 
2614     while(input_more_of_volume(vol, &input_info, &f));
2615 
2616     for(i = 0;i < mri->width;i++) {
2617       for(j = 0;j < mri->height;j++) {
2618         for(k = 0;k < mri->depth;k++) {
2619           for(t = 0;t < mri->nframes;t++) {
2620             val = get_volume_voxel_value(vol, i, j, k, t, 0);
2621             switch(mri->type){
2622             case MRI_UCHAR: MRIseq_vox(mri,i,j,k,t) = (unsigned char)val;break;
2623             case MRI_SHORT: MRISseq_vox(mri,i,j,k,t) = (short)val;break;
2624             case MRI_LONG:  MRILseq_vox(mri,i,j,k,t) = (long)val;break;
2625             case MRI_FLOAT: MRIFseq_vox(mri,i,j,k,t) = (float)val;break;
2626             }
2627           }
2628         }
2629       }
2630     }
2631   }
2632 
2633   delete_volume_input(&input_info);
2634   delete_volume(vol);
2635 
2636   return(mri);
2637 
2638 } /* end mincRead2() */
2639 /*-----------------------------------------------------------*/
2640 static int GetMINCInfo(MRI *mri,
2641                        char *dim_names[4],
2642                        int   dim_sizes[4],
2643                        Real separations[4],
2644                        Real dircos[3][3],
2645                        Real VolCenterVox[3],
2646                        Real VolCenterWorld[3])
2647 {
2648   int xspacehit, yspacehit, zspacehit;
2649   float col_dc_x, col_dc_y, col_dc_z;
2650   float row_dc_x, row_dc_y, row_dc_z;
2651   float slc_dc_x, slc_dc_y, slc_dc_z;
2652   Real col_dc_sign, row_dc_sign, slc_dc_sign;
2653   int err,i,j;
2654 
2655   col_dc_x = fabs(mri->x_r);
2656   col_dc_y = fabs(mri->x_a);
2657   col_dc_z = fabs(mri->x_s);
2658   col_dc_sign = 1;
2659 
2660   row_dc_x = fabs(mri->y_r);
2661   row_dc_y = fabs(mri->y_a);
2662   row_dc_z = fabs(mri->y_s);
2663   row_dc_sign = 1;
2664 
2665   slc_dc_x = fabs(mri->z_r);
2666   slc_dc_y = fabs(mri->z_a);
2667   slc_dc_z = fabs(mri->z_s);
2668   slc_dc_sign = 1;
2669 
2670   xspacehit = 0;
2671   yspacehit = 0;
2672   zspacehit = 0;
2673 
2674   /* Name the Column Axis */
2675   if(col_dc_x >= row_dc_x && col_dc_x >= slc_dc_x){
2676     dim_names[0] = MIxspace;
2677     VolCenterVox[0] = (mri->width-1)/2.0;
2678     if(mri->x_r < 0) col_dc_sign = -1;
2679     xspacehit = 1;
2680   }
2681   else if(col_dc_y >= row_dc_y && col_dc_y >= slc_dc_y){
2682     dim_names[0] = MIyspace;
2683     VolCenterVox[0] = (mri->height-1)/2.0;
2684     if(mri->x_a < 0) col_dc_sign = -1;
2685     yspacehit = 1;
2686   }
2687   else if(col_dc_z >= row_dc_z && col_dc_z >= slc_dc_z){
2688     dim_names[0] = MIzspace;
2689     VolCenterVox[0] = (mri->depth-1)/2.0;
2690     if(mri->x_s < 0) col_dc_sign = -1;
2691     zspacehit = 1;
2692   }
2693 
2694   /* Name the Row Axis */
2695   if(!xspacehit && row_dc_x >= slc_dc_x){
2696     dim_names[1] = MIxspace;
2697     VolCenterVox[1] = (mri->width-1)/2.0;
2698     if(mri->y_r < 0) row_dc_sign = -1;
2699     xspacehit = 1;
2700   }
2701   else if(!yspacehit && row_dc_y >= slc_dc_y){
2702     dim_names[1] = MIyspace;
2703     VolCenterVox[1] = (mri->height-1)/2.0;
2704     if(mri->y_a < 0) row_dc_sign = -1;
2705     yspacehit = 1;
2706   }
2707   else if(!zspacehit && row_dc_z >= slc_dc_z){
2708     dim_names[1] = MIzspace;
2709     VolCenterVox[1] = (mri->depth-1)/2.0;
2710     if(mri->y_s < 0) row_dc_sign = -1;
2711     zspacehit = 1;
2712   }
2713 
2714   /* Name the Slice Axis */
2715   if(!xspacehit){
2716     dim_names[2] = MIxspace;
2717     VolCenterVox[2] = (mri->width-1)/2.0;
2718     if(mri->z_r < 0) slc_dc_sign = -1;
2719     xspacehit = 1;
2720   }
2721   else if(!yspacehit){
2722     dim_names[2] = MIyspace;
2723     VolCenterVox[2] = (mri->height-1)/2.0;
2724     if(mri->z_a < 0) slc_dc_sign = -1;
2725     yspacehit = 1;
2726   }
2727   if(!zspacehit){
2728     dim_names[2] = MIzspace;
2729     VolCenterVox[2] = (mri->depth-1)/2.0;
2730     if(mri->z_s < 0) slc_dc_sign = -1;
2731     zspacehit = 1;
2732   }
2733 
2734   /* Check for errors in the Axis Naming*/
2735   err = 0;
2736   if(!xspacehit){
2737     printf("ERROR: could not assign xspace\n");
2738     err = 1;
2739   }
2740   if(!yspacehit){
2741     printf("ERROR: could not assign yspace\n");
2742     err = 1;
2743   }
2744   if(!zspacehit){
2745     printf("ERROR: could not assign zspace\n");
2746     err = 1;
2747   }
2748   if(err) return(1);
2749 
2750   /* Name the Frame Axis */
2751   dim_names[3] = MItime;
2752 
2753   /* Set world center */
2754   VolCenterWorld[0] = mri->c_r;
2755   VolCenterWorld[1] = mri->c_a;
2756   VolCenterWorld[2] = mri->c_s;
2757 
2758   /* Set dimension lengths */
2759   dim_sizes[0] = mri->width;
2760   dim_sizes[1] = mri->height;
2761   dim_sizes[2] = mri->depth;
2762   dim_sizes[3] = mri->nframes;
2763 
2764   /* Set separations */
2765   separations[0] = col_dc_sign * mri->xsize;
2766   separations[1] = row_dc_sign * mri->ysize;
2767   separations[2] = slc_dc_sign * mri->zsize;
2768   separations[3] = mri->tr;
2769 
2770   /* Set direction Cosines */
2771   dircos[0][0] = col_dc_sign * mri->x_r;
2772   dircos[0][1] = col_dc_sign * mri->x_a;
2773   dircos[0][2] = col_dc_sign * mri->x_s;
2774 
2775   dircos[1][0] = row_dc_sign * mri->y_r;
2776   dircos[1][1] = row_dc_sign * mri->y_a;
2777   dircos[1][2] = row_dc_sign * mri->y_s;
2778 
2779   dircos[2][0] = slc_dc_sign * mri->z_r;
2780   dircos[2][1] = slc_dc_sign * mri->z_a;
2781   dircos[2][2] = slc_dc_sign * mri->z_s;
2782 
2783   /* This is a hack for the case where the dircos are the default.
2784      The MINC routines will not save the direction cosines as
2785      NetCDF attriubtes if they correspond to the default. */
2786   for(i=0;i<3;i++)
2787     for(j=0;j<3;j++)
2788       if(fabs(dircos[i][j] < .00000000001)) dircos[i][j] = .0000000000000001;
2789 
2790   return(0);
2791 }
2792 
2793 /*-----------------------------------------------------------*/
2794 static int mincWrite2(MRI *mri, char *fname)
2795 {
2796   Volume minc_volume;
2797   nc_type nc_data_type = NC_BYTE;
2798   int ndim,i,j;
2799   Real min, max;
2800   float fmin, fmax;
2801   char *dim_names[4];
2802   int   dim_sizes[4];
2803   Real separations[4];
2804   Real dircos[3][3];
2805   Real VolCenterVox[3];
2806   Real VolCenterWorld[3];
2807   int signed_flag = 0;
2808   int r,c,s,f;
2809   Real VoxVal = 0.0;
2810   int return_value;
2811   VIO_Status status;
2812   MATRIX *T;
2813   Transform *pVox2WorldLin;
2814   General_transform *pVox2WorldGen;
2815 
2816   /* Get the min and max for the volume */
2817   if((return_value = MRIlimits(mri, &fmin, &fmax)) != NO_ERROR)
2818     return(return_value);
2819   min = (Real)fmin;
2820   max = (Real)fmax;
2821 
2822   /* Translate mri type to NetCDF type */
2823   switch(mri->type){
2824   case MRI_UCHAR: nc_data_type = NC_BYTE;  signed_flag = 0; break;
2825   case MRI_SHORT: nc_data_type = NC_SHORT; signed_flag = 1; break;
2826   case MRI_INT:   nc_data_type = NC_LONG;  signed_flag = 1; break;
2827   case MRI_LONG:  nc_data_type = NC_LONG;  signed_flag = 1; break;
2828   case MRI_FLOAT: nc_data_type = NC_FLOAT; signed_flag = 1; break;
2829   }
2830 
2831   /* Assign default direction cosines, if needed */
2832   if(mri->ras_good_flag == 0){
2833     setDirectionCosine(mri, MRI_CORONAL);
2834   }
2835 
2836   GetMINCInfo(mri, dim_names, dim_sizes, separations, dircos,
2837               VolCenterVox,VolCenterWorld);
2838 
2839   if(mri->nframes > 1) ndim = 4;
2840   else                 ndim = 3;
2841 
2842   minc_volume = create_volume(ndim, dim_names, nc_data_type,
2843                               signed_flag, min, max);
2844   set_volume_sizes(minc_volume, dim_sizes);
2845   alloc_volume_data(minc_volume);
2846 
2847   for(i=0;i<3;i++){
2848     printf("%d %s %4d %7.3f  %7.3f %7.3f %7.3f \n",
2849            i,dim_names[i],dim_sizes[i],separations[i],
2850            dircos[i][0],dircos[i][1],dircos[i][2]);
2851     set_volume_direction_cosine(minc_volume, i, dircos[i]);
2852   }
2853   printf("Center Voxel: %7.3f %7.3f %7.3f\n",
2854          VolCenterVox[0],VolCenterVox[1],VolCenterVox[2]);
2855   printf("Center World: %7.3f %7.3f %7.3f\n",
2856          VolCenterWorld[0],VolCenterWorld[1],VolCenterWorld[2]);
2857 
2858   set_volume_separations(minc_volume, separations);
2859   set_volume_translation(minc_volume, VolCenterVox, VolCenterWorld);
2860 
2861   T = MRIxfmCRS2XYZ(mri,0);
2862   printf("MRI Transform -----------------------------\n");
2863   MatrixPrint(stdout,T);
2864   pVox2WorldGen = get_voxel_to_world_transform(minc_volume);
2865   pVox2WorldLin = get_linear_transform_ptr(pVox2WorldGen);
2866   if ((Gdiag & DIAG_SHOW) && DIAG_VERBOSE_ON)
2867     {
2868       printf("MINC Linear Transform ----------------------\n");
2869       for(i=0;i<4;i++){
2870         for(j=0;j<4;j++) printf("%7.4f ",pVox2WorldLin->m[j][i]);
2871         printf("\n");
2872       }
2873       printf("--------------------------------------------\n");
2874     }
2875   MatrixFree(&T);
2876 
2877   printf("Setting Volume Values\n");
2878   for(f = 0; f < mri->nframes; f++) {     /* frames */
2879     for(c = 0; c < mri->width;   c++) {     /* columns */
2880       for(r = 0; r < mri->height;  r++) {     /* rows */
2881         for(s = 0; s < mri->depth;   s++) {     /* slices */
2882 
2883           switch(mri->type){
2884           case MRI_UCHAR: VoxVal = (Real)MRIseq_vox(mri,  c, r, s, f);break;
2885           case MRI_SHORT: VoxVal = (Real)MRISseq_vox(mri, c, r, s, f);break;
2886           case MRI_INT:   VoxVal = (Real)MRIIseq_vox(mri, c, r, s, f);break;
2887           case MRI_LONG:  VoxVal = (Real)MRILseq_vox(mri, c, r, s, f);break;
2888           case MRI_FLOAT: VoxVal = (Real)MRIFseq_vox(mri, c, r, s, f);break;
2889           }
2890           set_volume_voxel_value(minc_volume, c, r, s, f, 0, VoxVal);
2891         }
2892       }
2893     }
2894   }
2895 
2896   printf("Writing Volume\n");
2897   status = output_volume((STRING)fname, nc_data_type, signed_flag, min, max,
2898                          minc_volume, (STRING)"", NULL);
2899   printf("Cleaning Up\n");
2900   delete_volume(minc_volume);
2901 
2902   if(status) {
2903     printf("ERROR: mincWrite: output_volume exited with %d\n",status);
2904     return(1);
2905   }
2906 
2907   printf("mincWrite2: done\n");
2908   return(NO_ERROR);
2909 }
2910 
2911 
2912 /*-------------------------------------------------------
2913   NormalizeVector() - in-place vector normalization
2914   -------------------------------------------------------*/
2915 static int NormalizeVector(float *v, int n)
2916 {
2917   float sum2;
2918   int i;
2919   sum2 = 0;
2920   for(i=0;i<n;i++) sum2 += v[i];
2921   sum2 = sqrt(sum2);
2922   for(i=0;i<n;i++) v[i] /= sum2;
2923   return(0);
2924 }
2925 #endif
2926 /*----------------------------------------------------------*/
2927 /* time course clean */
2928 static int mincWrite(MRI *mri, char *fname)
2929 {
2930 
2931   Volume minc_volume;
2932   STRING dimension_names[4] = { "xspace", "yspace", "zspace", "time" };
2933   nc_type nc_data_type;
2934   Real min, max;
2935   float fmin, fmax;
2936   int dimension_sizes[4];
2937   Real separations[4];
2938   Real dir_cos[4];
2939   int return_value;
2940   Real voxel[4], world[4];
2941   int signed_flag;
2942   int di_x, di_y, di_z;
2943   int vi[4];
2944   /*   int r, a, s; */
2945   /*   float r_max; */
2946   VIO_Status status;
2947 
2948   /* di gives the bogus minc index
2949      di[0] is width, 1 is height, 2 is depth, 3 is time if
2950      there is a time dimension of length > 1
2951      minc wants these to be ras */
2952 
2953   /* here: minc volume index 0 is r, 1 is a, 2 is s */
2954   /* mri is lia *//* (not true for some volumes *)*/
2955 
2956   /* ----- get the orientation of the volume ----- */
2957   if(mri->ras_good_flag == 0)
2958     {
2959       setDirectionCosine(mri, MRI_CORONAL);
2960     }
2961 
2962   /* we don't muck around axes anymore */
2963   di_x = 0;
2964   di_y = 1;
2965   di_z = 2;
2966 
2967 #if 0
2968   /*  The following remapping is only valid for COR files */
2969   /*  In order to handle arbitrary volumes, we don't muck around */
2970   /*  the axes anymore  ... tosa */
2971 
2972   /* orig->minc map ***********************************/
2973   /* r axis is mapped to (x_r, y_r, z_r) voxel coords */
2974   /* thus find the biggest value among x_r, y_r, z_r  */
2975   /* that one corresponds to the r-axis               */
2976   r = 0;
2977   r_max = fabs(mri->x_r);
2978   if(fabs(mri->y_r) > r_max)
2979     {
2980       r_max = fabs(mri->y_r);
2981       r = 1;
2982     }
2983   if(fabs(mri->z_r) > r_max)
2984     r = 2;
2985 
2986   /* a axis is mapped to (x_a, y_a, z_a) voxel coords */
2987   if(r == 0)
2988     a = (fabs(mri->y_a) > fabs(mri->z_a) ? 1 : 2);
2989   else if(r == 1)
2990     a = (fabs(mri->x_a) > fabs(mri->z_a) ? 0 : 2);
2991   else
2992     a = (fabs(mri->x_a) > fabs(mri->y_a) ? 0 : 1);
2993 
2994   /* s axis is mapped to (x_s, y_s, z_s) voxel coords */
2995   /* use the rest to figure                           */
2996   s = 3 - r - a;
2997 
2998   /* ----- set the appropriate minc axes to this orientation ----- */
2999   /* r gives the mri structure axis of the lr coordinate */
3000   /* lr = minc xspace = 0 */
3001   /* a ..... of the pa coordinate *//* pa = minc yspace = 1 */
3002   /* s ..... of the is coordinate *//* is = minc zspace = 2 */
3003 
3004   /* di of this axis must be set to 2 */
3005   /* ... and so on */
3006 
3007   /* minc->orig map **************************************/
3008   /* you don't need the routine above but do a similar thing */
3009   /* to get exactly the same, i.e.                           */
3010   /* x-axis corresponds to (x_r, x_a, x_s) minc coords */
3011   /* y-axis                (y_r, y_a, y_s) minc coords */
3012   /* z-axis                (z_r, z_a, z_s) minc coords */
3013   /* thus we are doing too much work                   */
3014 
3015   if(r == 0)
3016     {
3017       if(a == 1)
3018         {
3019           di_x = 0;
3020           di_y = 1;
3021           di_z = 2;
3022         }
3023       else
3024         {
3025           di_x = 0;
3026           di_y = 2;
3027           di_z = 1;
3028         }
3029     }
3030   else if(r == 1)
3031     {
3032       if(a == 0)
3033         {
3034           di_x = 1;
3035           di_y = 0;
3036           di_z = 2;
3037         }
3038       else
3039         {
3040           di_x = 2;
3041           di_y = 0;
3042           di_z = 1;
3043         }
3044     }
3045   else
3046     {
3047       if(a == 0)
3048         {
3049           di_x = 1;
3050           di_y = 2;
3051           di_z = 0;
3052         }
3053       else
3054         {
3055           di_x = 2;
3056           di_y = 1;
3057           di_z = 0;
3058         }
3059     }
3060 #endif
3061 
3062   /* ----- set the data type ----- */
3063   if(mri->type == MRI_UCHAR)
3064     {
3065       nc_data_type = NC_BYTE;
3066       signed_flag = 0;
3067     }
3068   else if(mri->type == MRI_SHORT)
3069     {
3070       nc_data_type = NC_SHORT;
3071       signed_flag = 1;
3072     }
3073   else if(mri->type == MRI_INT)
3074     {
3075       nc_data_type = NC_LONG;
3076       signed_flag = 1;
3077     }
3078   else if(mri->type == MRI_LONG)
3079     {
3080       nc_data_type = NC_LONG;
3081       signed_flag = 1;
3082     }
3083   else if(mri->type == MRI_FLOAT)
3084     {
3085       nc_data_type = NC_FLOAT;
3086       signed_flag = 1;
3087     }
3088   else
3089     {
3090       errno = 0;
3091       ErrorReturn
3092         (ERROR_BADPARM,
3093          (ERROR_BADPARM,
3094           "mincWrite(): bad data type (%d) in mri structure",
3095           mri->type));
3096     }
3097 
3098   if((return_value = MRIlimits(mri, &fmin, &fmax)) != NO_ERROR)
3099     return(return_value);
3100 
3101   min = (Real)fmin;
3102   max = (Real)fmax;
3103 
3104   if(mri->nframes == 1)
3105     minc_volume = create_volume
3106       (3, dimension_names, nc_data_type, signed_flag, min, max);
3107   else
3108     minc_volume = create_volume
3109       (4, dimension_names, nc_data_type, signed_flag, min, max);
3110 
3111   /* di_(x,y,z) is the map from minc to orig */
3112   /* minc dimension size is that of di_x, etc. */
3113   dimension_sizes[di_x] = mri->width;
3114   dimension_sizes[di_y] = mri->height;
3115   dimension_sizes[di_z] = mri->depth;
3116   dimension_sizes[3] = mri->nframes;
3117 
3118   set_volume_sizes(minc_volume, dimension_sizes);
3119 
3120   alloc_volume_data(minc_volume);
3121 
3122   separations[di_x] = (Real)(mri->xsize);
3123   separations[di_y] = (Real)(mri->ysize);
3124   separations[di_z] = (Real)(mri->zsize);
3125   separations[3] = 1.0; // appears to do nothing
3126   set_volume_separations(minc_volume, separations);
3127   /* has side effect to change transform and thus must be set first */
3128 
3129   dir_cos[0] = (Real)mri->x_r;
3130   dir_cos[1] = (Real)mri->x_a;
3131   dir_cos[2] = (Real)mri->x_s;
3132   set_volume_direction_cosine(minc_volume, di_x, dir_cos);
3133 
3134   dir_cos[0] = (Real)mri->y_r;
3135   dir_cos[1] = (Real)mri->y_a;
3136   dir_cos[2] = (Real)mri->y_s;
3137   set_volume_direction_cosine(minc_volume, di_y, dir_cos);
3138 
3139   dir_cos[0] = (Real)mri->z_r;
3140   dir_cos[1] = (Real)mri->z_a;
3141   dir_cos[2] = (Real)mri->z_s;
3142   set_volume_direction_cosine(minc_volume, di_z, dir_cos);
3143 
3144   voxel[di_x] = mri->width / 2.0; // promoted to double
3145   voxel[di_y] = mri->height/ 2.0;
3146   voxel[di_z] = mri->depth / 2.0;
3147   voxel[3] = 0.0;
3148   world[0] = (Real)(mri->c_r);
3149   world[1] = (Real)(mri->c_a);
3150   world[2] = (Real)(mri->c_s);
3151   world[3] = 0.0;
3152   set_volume_translation(minc_volume, voxel, world);
3153 
3154   /* get the position from (vi[di_x], vi[di_y], vi[di_z]) orig position     */
3155   /*      put the value to (vi[0], vi[1], vi[2]) minc volume                */
3156   /* this will keep the orientation of axes the same as the original        */
3157 
3158   /* vi[n] gives the index of the variable along minc axis x */
3159   /* vi[di_x] gives the index of the variable
3160      along minc axis di_x, or along mri axis x */
3161   for(vi[3] = 0;vi[3] < mri->nframes;vi[3]++) {           /* frames */
3162     for(vi[di_x] = 0;vi[di_x] < mri->width;vi[di_x]++) {   /* columns */
3163       for(vi[di_y] = 0;vi[di_y] < mri->height;vi[di_y]++) { /* rows */
3164         for(vi[di_z] = 0;vi[di_z] < mri->depth;vi[di_z]++) { /* slices */
3165 
3166           if(mri->type == MRI_UCHAR)
3167             set_volume_voxel_value
3168               (minc_volume, vi[0], vi[1], vi[2], vi[3], 0,
3169                (Real)MRIseq_vox(mri, vi[di_x], vi[di_y], vi[di_z], vi[3]));
3170           if(mri->type == MRI_SHORT)
3171             set_volume_voxel_value
3172               (minc_volume, vi[0], vi[1], vi[2], vi[3], 0,
3173                (Real)MRISseq_vox(mri, vi[di_x], vi[di_y], vi[di_z], vi[3]));
3174           if(mri->type == MRI_INT)
3175             set_volume_voxel_value
3176               (minc_volume, vi[0], vi[1], vi[2], vi[3], 0,
3177                (Real)MRIIseq_vox(mri, vi[di_x], vi[di_y], vi[di_z], vi[3]));
3178           if(mri->type == MRI_LONG)
3179             set_volume_voxel_value
3180               (minc_volume, vi[0], vi[1], vi[2], vi[3], 0,
3181                (Real)MRILseq_vox(mri, vi[di_x], vi[di_y], vi[di_z], vi[3]));
3182           if(mri->type == MRI_FLOAT)
3183             set_volume_voxel_value
3184               (minc_volume, vi[0], vi[1], vi[2], vi[3], 0,
3185                (Real)MRIFseq_vox(mri, vi[di_x], vi[di_y], vi[di_z], vi[3]));
3186         }
3187       }
3188     }
3189   }
3190 
3191   status = output_volume((STRING)fname, nc_data_type, signed_flag, min, max,
3192                          minc_volume, (STRING)"", NULL);
3193   delete_volume(minc_volume);
3194 
3195   if(status) {
3196     printf("ERROR: mincWrite: output_volume exited with %d\n",status);
3197     return(1);
3198   }
3199 
3200   return(NO_ERROR);
3201 
3202 } /* end mincWrite() */
3203 
3204 
3205 /*----------------------------------------------------------------
3206   bvolumeWrite() - replaces bshortWrite and bfloatWrite.
3207   Bug: fname_passed must the the stem, not the full file name.
3208   -----------------------------------------------------------------*/
3209 static int bvolumeWrite(MRI *vol, char *fname_passed, int type)
3210 {
3211 
3212   int i, j, t;
3213   char fname[STRLEN];
3214   short *bufshort;
3215   float *buffloat;
3216   FILE *fp;
3217   int result;
3218   MRI *subject_info = NULL;
3219   char subject_volume_dir[STRLEN];
3220   char *subjects_dir;
3221   char *sn;
3222   char analyse_fname[STRLEN], register_fname[STRLEN];
3223   char output_dir[STRLEN];
3224   char *c;
3225   int od_length;
3226   char t1_path[STRLEN];
3227   MATRIX *cf, *bf, *ibf, *af, *iaf, *as, *bs, *cs, *ics, *r;
3228   MATRIX *r1, *r2, *r3, *r4;
3229   float det;
3230   int bad_flag;
3231   int l;
3232   char stem[STRLEN];
3233   char *c1, *c2, *c3;
3234   struct stat stat_buf;
3235   char subject_dir[STRLEN];
3236   int dealloc, nslices, nframes;
3237   MRI *mri;
3238   float min,max;
3239   int swap_bytes_flag, size, bufsize;
3240   char *ext;
3241   void *buf;
3242 
3243   /* check the type and set the extension and size*/
3244   switch(type) {
3245   case MRI_SHORT: ext = "bshort"; size = sizeof(short); break;
3246   case MRI_FLOAT: ext = "bfloat"; size = sizeof(float); break;
3247   default:
3248     fprintf
3249       (stderr,"ERROR: bvolumeWrite: type (%d) is not short or float\n", type);
3250     return(1);
3251   }
3252 
3253   if(vol->type != type){
3254     if (DIAG_VERBOSE_ON)
3255       printf("INFO: bvolumeWrite: changing type\n");
3256     nslices = vol->depth;
3257     nframes = vol->nframes;
3258     vol->depth = nslices*nframes;
3259     vol->nframes = 1;
3260     MRIlimits(vol,&min,&max);
3261     if (DIAG_VERBOSE_ON)
3262       printf("INFO: bvolumeWrite: range %g %g\n",min,max);
3263     mri = MRIchangeType(vol,type,min,max,1);
3264     if(mri == NULL) {
3265       fprintf(stderr,"ERROR: bvolumeWrite: MRIchangeType\n");
3266       return(1);
3267     }
3268     vol->depth = nslices;
3269     vol->nframes = nframes;
3270     mri->depth = nslices;
3271     mri->nframes = nframes;
3272     dealloc = 1;
3273   }
3274   else{
3275     mri = vol;
3276     dealloc = 0;
3277   }
3278 
3279   /* ----- get the stem from the passed file name ----- */
3280   /*
3281     four options:
3282     1. stem_xxx.bshort
3283     2. stem.bshort
3284     3. stem_xxx
3285     4. stem
3286     other possibles:
3287     stem_.bshort
3288   */
3289   l = strlen(fname_passed);
3290   c1 = fname_passed + l - 11;
3291   c2 = fname_passed + l - 7;
3292   c3 = fname_passed + l - 4;
3293   strcpy(stem, fname_passed);
3294   if(c1 > fname_passed)
3295     {
3296       if( (*c1 == '_' && strcmp(c1+4, ".bshort") == 0) ||
3297           (*c1 == '_' && strcmp(c1+4, ".bfloat") == 0) )
3298         stem[(int)(c1-fname_passed)] = '\0';
3299     }
3300   if(c2 > fname_passed)
3301     {
3302       if( (*c2 == '_' && strcmp(c2+4, ".bshort") == 0) ||
3303           (*c2 == '_' && strcmp(c2+4, ".bfloat") == 0) )
3304         stem[(int)(c2-fname_passed)] = '\0';
3305 
3306     }
3307   if(c3 > fname_passed)
3308     {
3309       if(*c3 == '_')
3310         stem[(int)(c3-fname_passed)] = '\0';
3311     }
3312   printf("INFO: bvolumeWrite: stem = %s\n",stem);
3313 
3314   c = strrchr(stem, '/');
3315   if(c == NULL)
3316     output_dir[0] = '\0';
3317   else
3318     {
3319       od_length = (int)(c - stem);
3320       strncpy(output_dir, stem, od_length);
3321       /* -- leaving the trailing '/' on a directory is not my
3322          usual convention, but here it's a load easier if there's
3323          no directory in stem... -ch -- */
3324       output_dir[od_length] = '/';
3325       output_dir[od_length+1] = '\0';
3326     }
3327 
3328   sprintf(analyse_fname, "%s%s", output_dir, "analyse.dat");
3329   sprintf(register_fname, "%s%s", output_dir, "register.dat");
3330 
3331   bufsize  = mri->width * size;
3332   bufshort = (short *)malloc(bufsize);
3333   buffloat = (float *)malloc(bufsize);
3334 
3335   buf = NULL; /* shuts up compiler */
3336   if(type == MRI_SHORT)  buf = bufshort;
3337   else                   buf = buffloat;
3338 
3339   swap_bytes_flag = 0;
3340 #if (BYTE_ORDER==LITTLE_ENDIAN)
3341   swap_bytes_flag = 1;
3342 #endif
3343 
3344   //printf("--------------------------------\n");
3345   //MRIdump(mri,stdout);
3346   //MRIdumpBuffer(mri,stdout);
3347   //printf("--------------------------------\n");
3348 
3349   for(i = 0;i < mri->depth;i++)
3350     {
3351 
3352       /* ----- write the header file ----- */
3353       sprintf(fname, "%s_%03d.hdr", stem, i);
3354       if((fp = fopen(fname, "w")) == NULL)
3355         {
3356           if(dealloc) MRIfree(&mri);
3357           free(bufshort);
3358           free(buffloat);
3359           errno = 0;
3360           ErrorReturn
3361             (ERROR_BADFILE, (ERROR_BADFILE,
3362                              "bvolumeWrite(): can't open file %s", fname));
3363         }
3364       fprintf(fp, "%d %d %d %d\n", mri->height, mri->width, mri->nframes, 0);
3365       fclose(fp);
3366 
3367       /* ----- write the data file ----- */
3368       sprintf(fname, "%s_%03d.%s", stem, i, ext);
3369       if((fp = fopen(fname, "w")) == NULL)
3370         {
3371           if(dealloc) MRIfree(&mri);
3372           free(bufshort);
3373           free(buffloat);
3374           errno = 0;
3375           ErrorReturn
3376             (ERROR_BADFILE, (ERROR_BADFILE,
3377                              "bvolumeWrite(): can't open file %s", fname));
3378         }
3379 
3380       for(t = 0;t < mri->nframes;t++) {
3381         for(j = 0;j < mri->height;j++) {
3382           memcpy(buf,mri->slices[t*mri->depth + i][j], bufsize);
3383 
3384           if(swap_bytes_flag){
3385             if(type == MRI_SHORT) byteswapbufshort((void*)buf, bufsize);
3386             else                  byteswapbuffloat((void*)buf, bufsize);
3387           }
3388           fwrite(buf, size, mri->width, fp);
3389         }
3390       }
3391 
3392       fclose(fp);
3393     }
3394 
3395   free(bufshort);
3396   free(buffloat);
3397 
3398   sn = subject_name;
3399   if(mri->subject_name[0] != '\0')
3400     sn = mri->subject_name;
3401 
3402   if(sn != NULL)
3403     {
3404       if((subjects_dir = getenv("SUBJECTS_DIR")) == NULL)
3405         {
3406           errno = 0;
3407           ErrorPrintf
3408             (ERROR_BADPARM,
3409              "bvolumeWrite(): environment variable SUBJECTS_DIR unset");
3410           if(dealloc) MRIfree(&mri);
3411         }
3412       else
3413         {
3414 
3415           sprintf(subject_dir, "%s/%s", subjects_dir, sn);
3416           if(stat(subject_dir, &stat_buf) < 0)
3417             {
3418               fprintf
3419                 (stderr,
3420                  "can't stat %s; writing to bhdr instead\n", subject_dir);
3421             }
3422           else
3423             {
3424               if(!S_ISDIR(stat_buf.st_mode))
3425                 {
3426                   fprintf
3427                     (stderr,
3428                      "%s is not a directory; writing to bhdr instead\n",
3429                      subject_dir);
3430                 }
3431               else
3432                 {
3433                   sprintf(subject_volume_dir, "%s/mri/T1", subject_dir);
3434                   subject_info = MRIreadInfo(subject_volume_dir);
3435                   if(subject_info == NULL)
3436                     {
3437                       sprintf(subject_volume_dir, "%s/mri/orig", subject_dir);
3438                       subject_info = MRIreadInfo(subject_volume_dir);
3439                       if(subject_info == NULL)
3440                         fprintf(stderr,
3441                                 "can't read the subject's orig or T1 volumes; "
3442                                 "writing to bhdr instead\n");
3443                     }
3444                 }
3445             }
3446 
3447         }
3448 
3449     }
3450 
3451   if(subject_info != NULL)
3452     {
3453       if(subject_info->ras_good_flag == 0)
3454         {
3455           setDirectionCosine(subject_info, MRI_CORONAL);
3456         }
3457     }
3458 
3459   cf = bf = ibf = af = iaf = as = bs = cs = ics = r = NULL;
3460   r1 = r2 = r3 = r4 = NULL;
3461 
3462   /* ----- write the register.dat and analyse.dat  or bhdr files ----- */
3463   if(subject_info != NULL)
3464     {
3465 
3466       bad_flag = FALSE;
3467 
3468       if((as = MatrixAlloc(4, 4, MATRIX_REAL)) == NULL)
3469         {
3470           errno = 0;
3471           ErrorPrintf(ERROR_BADPARM, "bvolumeWrite(): error creating matrix");
3472           bad_flag = TRUE;
3473         }
3474       stuff_four_by_four
3475         (as, subject_info->x_r, subject_info->y_r, subject_info->z_r, subject_info->c_r,
3476          subject_info->y_r, subject_info->y_r, subject_info->y_r, subject_info->c_r,
3477          subject_info->z_r, subject_info->z_r, subject_info->z_r, subject_info->c_r,
3478          0.0,               0.0,               0.0,               1.0);
3479 
3480       if((af = MatrixAlloc(4, 4, MATRIX_REAL)) == NULL)
3481         {
3482           errno = 0;
3483           ErrorPrintf(ERROR_BADPARM, "bvolumeWrite(): error creating matrix");
3484           bad_flag = TRUE;
3485         }
3486       stuff_four_by_four(af, mri->x_r, mri->y_r, mri->z_r, mri->c_r,
3487                          mri->y_r, mri->y_r, mri->y_r, mri->c_r,
3488                          mri->z_r, mri->z_r, mri->z_r, mri->c_r,
3489                          0.0,      0.0,      0.0,      1.0);
3490 
3491       if((bs = MatrixAlloc(4, 4, MATRIX_REAL)) == NULL)
3492         {
3493           errno = 0;
3494           ErrorPrintf(ERROR_BADPARM, "bvolumeWrite(): error creating matrix");
3495           bad_flag = TRUE;
3496         }
3497       stuff_four_by_four(bs, 1, 0, 0, (subject_info->width  - 1) / 2.0,
3498                          0, 1, 0, (subject_info->height - 1) / 2.0,
3499                          0, 0, 1, (subject_info->depth  - 1) / 2.0,
3500                          0, 0, 0,                              1.0);
3501 
3502       if((bf = MatrixAlloc(4, 4, MATRIX_REAL)) == NULL)
3503         {
3504           errno = 0;
3505           ErrorPrintf(ERROR_BADPARM, "bvolumeWrite(): error creating matrix");
3506           bad_flag = TRUE;
3507         }
3508       stuff_four_by_four(bf, 1, 0, 0, (mri->width  - 1) / 2.0,
3509                          0, 1, 0, (mri->height - 1) / 2.0,
3510                          0, 0, 1, (mri->depth  - 1) / 2.0,
3511                          0, 0, 0,                     1.0);
3512 
3513       if((cs = MatrixAlloc(4, 4, MATRIX_REAL)) == NULL)
3514         {
3515           errno = 0;
3516           ErrorPrintf(ERROR_BADPARM, "bvolumeWrite(): error creating matrix");
3517           bad_flag = TRUE;
3518         }
3519       stuff_four_by_four
3520         (cs,
3521          -subject_info->xsize, 0, 0,(subject_info->width  * mri->xsize) / 2.0,\
3522          0, 0, subject_info->zsize, -(subject_info->depth  * mri->zsize) / 2.0, \
3523          0, -subject_info->ysize, 0,  (subject_info->height * mri->ysize) / 2.0, \
3524          0, 0, 0, 1);
3525 
3526       if((cf = MatrixAlloc(4, 4, MATRIX_REAL)) == NULL)
3527         {
3528           errno = 0;
3529           ErrorPrintf(ERROR_BADPARM, "bvolumeWrite(): error creating matrix");
3530           bad_flag = TRUE;
3531         }
3532       stuff_four_by_four
3533         (cf,
3534          -mri->xsize, 0,          0,  (mri->width  * mri->xsize) / 2.0,
3535          0,           0, mri->zsize, -(mri->depth  * mri->zsize) / 2.0,
3536          0, -mri->ysize,          0,  (mri->height * mri->ysize) / 2.0,
3537          0,           0,          0,                                 1);
3538 
3539       if(bad_flag)
3540         {
3541           errno = 0;
3542           ErrorPrintf
3543             (ERROR_BADPARM,
3544              "bvolumeWrite(): error creating one "
3545              "or more matrices; aborting register.dat "
3546              "write and writing bhdr instead");
3547           MRIfree(&subject_info);
3548         }
3549 
3550     }
3551 
3552   if(subject_info != NULL)
3553     {
3554 
3555       bad_flag = FALSE;
3556 
3557       if((det = MatrixDeterminant(as)) == 0.0)
3558         {
3559           errno = 0;
3560           ErrorPrintf
3561             (ERROR_BADPARM,
3562              "bvolumeWrite(): bad determinant in "
3563              "matrix (check structural volume)");
3564           bad_flag = TRUE;
3565         }
3566       if((det = MatrixDeterminant(bs)) == 0.0)
3567         {
3568           errno = 0;
3569           ErrorPrintf
3570             (ERROR_BADPARM,
3571              "bvolumeWrite(): bad determinant in "
3572              "matrix (check structural volume)");
3573           bad_flag = TRUE;
3574         }
3575       if((det = MatrixDeterminant(cs)) == 0.0)
3576         {
3577           errno = 0;
3578           ErrorPrintf
3579             (ERROR_BADPARM,
3580              "bvolumeWrite(): bad determinant in "
3581              "matrix (check structural volume)");
3582           bad_flag = TRUE;
3583         }
3584 
3585       if((det = MatrixDeterminant(af)) == 0.0)
3586         {
3587           errno = 0;
3588           ErrorPrintf
3589             (ERROR_BADPARM,
3590              "bvolumeWrite(): bad determinant in "
3591              "matrix (check functional volume)");
3592           bad_flag = TRUE;
3593         }
3594       if((det = MatrixDeterminant(bf)) == 0.0)
3595         {
3596           errno = 0;
3597           ErrorPrintf
3598             (ERROR_BADPARM,
3599              "bvolumeWrite(): bad determinant in "
3600              "matrix (check functional volume)");
3601           bad_flag = TRUE;
3602         }
3603       if((det = MatrixDeterminant(cf)) == 0.0)
3604         {
3605           errno = 0;
3606           ErrorPrintf
3607             (ERROR_BADPARM,
3608              "bvolumeWrite(): bad determinant in "
3609              "matrix (check functional volume)");
3610           bad_flag = TRUE;
3611         }
3612 
3613       if(bad_flag)
3614         {
3615           errno = 0;
3616           ErrorPrintf
3617             (ERROR_BADPARM, "bvolumeWrite(): one or more zero determinants;"
3618              " aborting register.dat write and writing bhdr instead");
3619           MRIfree(&subject_info);
3620         }
3621 
3622     }
3623 
3624   if(subject_info != NULL)
3625     {
3626 
3627       bad_flag = FALSE;
3628 
3629       if((iaf = MatrixInverse(af, NULL)) == NULL)
3630         {
3631           errno = 0;
3632           ErrorPrintf(ERROR_BADPARM, "bvolumeWrite(): error inverting matrix");
3633           bad_flag = TRUE;
3634         }
3635       if((ibf = MatrixInverse(bf, NULL)) == NULL)
3636         {
3637           errno = 0;
3638           ErrorPrintf(ERROR_BADPARM, "bvolumeWrite(): error inverting matrix");
3639           bad_flag = TRUE;
3640         }
3641       if((ics = MatrixInverse(cs, NULL)) == NULL)
3642         {
3643           errno = 0;
3644           ErrorPrintf(ERROR_BADPARM, "bvolumeWrite(): error inverting matrix");
3645           bad_flag = TRUE;
3646         }
3647 
3648       if(bad_flag)
3649         {
3650           errno = 0;
3651           ErrorPrintf
3652             (ERROR_BADPARM, "bvolumeWrite(): one or more zero determinants; "
3653              "aborting register.dat write and writing bhdr instead");
3654           MRIfree(&subject_info);
3655         }
3656     }
3657 
3658   bad_flag = FALSE;
3659 
3660   if(subject_info != NULL)
3661     {
3662 
3663       if((r1 = MatrixMultiply(bs, ics, NULL)) == NULL)
3664         {
3665           bad_flag = TRUE;
3666           MRIfree(&subject_info);
3667         }
3668 
3669     }
3670 
3671   if(subject_info != NULL)
3672     {
3673 
3674       if((r2 = MatrixMultiply(as, r1, NULL)) == NULL)
3675         {
3676           bad_flag = TRUE;
3677           MRIfree(&subject_info);
3678         }
3679 
3680     }
3681 
3682   if(subject_info != NULL)
3683     {
3684 
3685       if((r3 = MatrixMultiply(iaf, r2, NULL)) == NULL)
3686         {
3687           bad_flag = TRUE;
3688           MRIfree(&subject_info);
3689         }
3690 
3691     }
3692 
3693   if(subject_info != NULL)
3694     {
3695 
3696       if((r4 = MatrixMultiply(ibf, r3, NULL)) == NULL)
3697         {
3698           bad_flag = TRUE;
3699           MRIfree(&subject_info);
3700         }
3701 
3702     }
3703 
3704   if(subject_info != NULL)
3705     {
3706 
3707       if((r = MatrixMultiply(cf, r4, NULL)) == NULL)
3708         {
3709           bad_flag = TRUE;
3710           MRIfree(&subject_info);
3711         }
3712 
3713     }
3714 
3715   if(bad_flag)
3716     {
3717       errno = 0;
3718       ErrorPrintf
3719         (ERROR_BADPARM, "bvolumeWrite(): error during matrix multiplications; "
3720          "aborting register.dat write and writing bhdr instead");
3721     }
3722 
3723   if( as != NULL)  MatrixFree( &as);
3724   if( bs != NULL)  MatrixFree( &bs);
3725   if( cs != NULL)  MatrixFree( &cs);
3726   if( af != NULL)  MatrixFree( &af);
3727   if( bf != NULL)  MatrixFree( &bf);
3728   if( cf != NULL)  MatrixFree( &cf);
3729   if(iaf != NULL)  MatrixFree(&iaf);
3730   if(ibf != NULL)  MatrixFree(&ibf);
3731   if(ics != NULL)  MatrixFree(&ics);
3732   if( r1 != NULL)  MatrixFree( &r1);
3733   if( r2 != NULL)  MatrixFree( &r2);
3734   if( r3 != NULL)  MatrixFree( &r3);
3735   if( r4 != NULL)  MatrixFree( &r4);
3736 
3737   if(subject_info != NULL)
3738     {
3739 
3740       if(mri->path_to_t1[0] == '\0')
3741         sprintf(t1_path, ".");
3742       else
3743         strcpy(t1_path, mri->path_to_t1);
3744 
3745       if(FileExists(analyse_fname))
3746         fprintf(stderr, "warning: overwriting file %s\n", analyse_fname);
3747 
3748       if((fp = fopen(analyse_fname, "w")) == NULL)
3749         {
3750           MRIfree(&subject_info);
3751           errno = 0;
3752           ErrorReturn
3753             (ERROR_BADFILE,
3754              (ERROR_BADFILE,
3755               "bvolumeWrite(): couldn't open file %s for writing",
3756               analyse_fname));
3757         }
3758 
3759       fprintf(fp, "%s\n", t1_path);
3760       fprintf(fp, "%s_%%03d.bshort\n", stem);
3761       fprintf(fp, "%d %d\n", mri->depth, mri->nframes);
3762       fprintf(fp, "%d %d\n", mri->width, mri->height);
3763 
3764       fclose(fp);
3765 
3766       if(FileExists(analyse_fname))
3767         fprintf(stderr, "warning: overwriting file %s\n", register_fname);
3768 
3769       if((fp = fopen(register_fname, "w")) == NULL)
3770         {
3771           MRIfree(&subject_info);
3772           errno = 0;
3773           ErrorReturn
3774             (ERROR_BADFILE,
3775              (ERROR_BADFILE,
3776               "bvolumeWrite(): couldn't open file %s for writing",
3777               register_fname));
3778         }
3779 
3780       fprintf(fp, "%s\n", sn);
3781       fprintf(fp, "%g\n", mri->xsize);
3782       fprintf(fp, "%g\n", mri->zsize);
3783       fprintf(fp, "%g\n", 1.0);
3784       fprintf(fp, "%g %g %g %g\n", \
3785               *MATRIX_RELT(r, 1, 1),
3786               *MATRIX_RELT(r, 1, 2),
3787               *MATRIX_RELT(r, 1, 3),
3788               *MATRIX_RELT(r, 1, 4));
3789       fprintf(fp, "%g %g %g %g\n", \
3790               *MATRIX_RELT(r, 2, 1),
3791               *MATRIX_RELT(r, 2, 2),
3792               *MATRIX_RELT(r, 2, 3),
3793               *MATRIX_RELT(r, 2, 4));
3794       fprintf(fp, "%g %g %g %g\n", \
3795               *MATRIX_RELT(r, 3, 1),
3796               *MATRIX_RELT(r, 3, 2),
3797               *MATRIX_RELT(r, 3, 3),
3798               *MATRIX_RELT(r, 3, 4));
3799       fprintf(fp, "%g %g %g %g\n", \
3800               *MATRIX_RELT(r, 4, 1),
3801               *MATRIX_RELT(r, 4, 2),
3802               *MATRIX_RELT(r, 4, 3),
3803               *MATRIX_RELT(r, 4, 4));
3804 
3805       fclose(fp);
3806 
3807       MatrixFree(&r);
3808 
3809     }
3810 
3811   if(subject_info == NULL)
3812     {
3813       sprintf(fname, "%s.bhdr", stem);
3814       if((fp = fopen(fname, "w")) == NULL)
3815         {
3816           if(dealloc) MRIfree(&mri);
3817           errno = 0;
3818           ErrorReturn
3819             (ERROR_BADFILE,
3820              (ERROR_BADFILE, "bvolumeWrite(): can't open file %s", fname));
3821         }
3822 
3823       result = write_bhdr(mri, fp);
3824 
3825       fclose(fp);
3826 
3827       if(result != NO_ERROR)
3828         return(result);
3829 
3830     }
3831   else
3832     MRIfree(&subject_info);
3833 
3834   if(dealloc) MRIfree(&mri);
3835 
3836   return(NO_ERROR);
3837 
3838 } /* end bvolumeWrite() */
3839 
3840 static MRI *get_b_info
3841 (char *fname_passed, int read_volume, char *directory, char *stem, int type)
3842 {
3843 
3844   MRI *mri, *mri2;
3845   FILE *fp;
3846   int nslices=0, nt;
3847   int nx, ny;
3848   int result;
3849   char fname[STRLEN];
3850   char extension[STRLEN];
3851   char bhdr_name[STRLEN];
3852 
3853   if(type == MRI_SHORT)      sprintf(extension, "bshort");
3854   else if(type == MRI_FLOAT) sprintf(extension, "bfloat");
3855   else{
3856     errno = 0;
3857     ErrorReturn(NULL, (ERROR_UNSUPPORTED,
3858                        "internal error: get_b_info() passed type %d", type));
3859   }
3860 
3861   result = decompose_b_fname(fname_passed, directory, stem);
3862   if(result != NO_ERROR) return(NULL);
3863 
3864   if(directory[0] == '\0') sprintf(directory, ".");
3865 
3866 #if 0
3867 
3868   char sn[STRLEN];
3869   char fname_descrip[STRLEN];
3870   float fov_x, fov_y, fov_z;
3871   MATRIX *m;
3872   int res;
3873   char register_fname[STRLEN], analyse_fname[STRLEN];
3874   float ipr, st, intensity;
3875   float m11, m12, m13, m14;
3876   float m21, m22, m23, m24;
3877   float m31, m32, m33, m34;
3878   float m41, m42, m43, m44;
3879   char t1_path[STRLEN];
3880 
3881   /* ----- try register.dat and analyse.dat, then bhdr, then defaults ----- */
3882   sprintf(register_fname, "%s/register.dat", directory);
3883   sprintf(analyse_fname, "%s/analyse.dat", directory);
3884 
3885   if((fp = fopen(register_fname, "r")) != NULL)
3886     {
3887 
3888       fscanf(fp, "%s", sn);
3889       fscanf(fp, "%f", &ipr);
3890       fscanf(fp, "%f", &st);
3891       fscanf(fp, "%f", &intensity);
3892       fscanf(fp, "%f %f %f %f", &m11, &m12, &m13, &m14);
3893       fscanf(fp, "%f %f %f %f", &m21, &m22, &m23, &m24);
3894       fscanf(fp, "%f %f %f %f", &m31, &m32, &m33, &m34);
3895       fscanf(fp, "%f %f %f %f", &m41, &m42, &m43, &m44);
3896       fclose(fp);
3897 
3898       if((fp = fopen(analyse_fname, "r")) != NULL)
3899         {
3900 
3901           fscanf(fp, "%s", t1_path);
3902           fscanf(fp, "%s", fname_descrip);
3903           fscanf(fp, "%d %d", &nslices, &nt);
3904           fscanf(fp, "%d %d", &nx, &ny);
3905           fclose(fp);
3906 
3907           if(read_volume)
3908             {
3909               mri = MRIallocSequence(nx, ny, nslices, MRI_SHORT, nt);
3910             }
3911           else
3912             {
3913               mri = MRIallocHeader(nx, ny, nslices, MRI_SHORT);
3914               mri->nframes = nt;
3915             }
3916 
3917           strcpy(mri->fname, fname_passed);
3918           mri->imnr0 = 1;
3919           mri->imnr1 = nslices;
3920           mri->xsize = mri->ysize = mri->ps = ipr;
3921           mri->zsize = mri->thick = st;
3922 
3923           fov_x = mri->xsize * mri->width;
3924           fov_y = mri->ysize * mri->height;
3925           fov_z = mri->zsize * mri->depth;
3926 
3927           mri->fov = (fov_x > fov_y ? (fov_x > fov_z ? fov_x : fov_z) :
3928                       (fov_y > fov_z ? fov_y : fov_z));
3929 
3930           mri->xend = fov_x / 2.0;
3931           mri->xstart = -mri->xend;
3932           mri->yend = fov_y / 2.0;
3933           mri->ystart = -mri->yend;
3934           mri->zend = fov_z / 2.0;
3935           mri->zstart = -mri->zend;
3936 
3937           mri->brightness = intensity;
3938           strcpy(mri->subject_name, sn);
3939           strcpy(mri->path_to_t1, t1_path);
3940           strcpy(mri->fname_format, fname_descrip);
3941 
3942           m = MatrixAlloc(4, 4, MATRIX_REAL);
3943           if(m == NULL)
3944             {
3945               MRIfree(&mri);
3946               errno = 0;
3947               ErrorReturn(NULL, (ERROR_NO_MEMORY, "error allocating "
3948                                  "matrix in %s read", extension));
3949             }
3950           stuff_four_by_four(m, m11, m12, m13, m14,
3951                              m21, m22, m23, m24,
3952                              m31, m32, m33, m34,
3953                              m41, m42, m43, m44);
3954           mri->register_mat = m;
3955           res = orient_with_register(mri);
3956           if(res != NO_ERROR){
3957             MRIfree(&mri);
3958             return(NULL);
3959           }
3960           return(mri);
3961         }
3962     }
3963 #endif
3964 
3965   mri = MRIallocHeader(1, 1, 1, type);
3966 
3967   /* ----- try to read the stem.bhdr ----- */
3968   sprintf(bhdr_name, "%s/%s.bhdr", directory, stem);
3969   if((fp = fopen(bhdr_name, "r")) != NULL)
3970     {
3971       read_bhdr(mri, fp);
3972       sprintf(fname, "%s/%s_000.hdr", directory, stem);
3973       if((fp = fopen(fname, "r")) == NULL){
3974         MRIfree(&mri);
3975         errno = 0;
3976         ErrorReturn(NULL, (ERROR_BADFILE, "cannot open %s", fname));
3977       }
3978       fscanf(fp, "%d %d %d %*d", &ny, &nx, &nt);
3979       mri->nframes = nt;
3980       fclose(fp);
3981     }
3982   else
3983     { /* ----- get defaults ----- */
3984       fprintf
3985         (stderr,
3986          "-----------------------------------------------------------------\n"
3987          "Could not find the direction cosine information.\n"
3988          "Will use the CORONAL orientation.\n"
3989          "If not suitable, please provide the information in %s file\n"
3990          "-----------------------------------------------------------------\n",
3991          bhdr_name);
3992       sprintf(fname, "%s/%s_000.hdr", directory, stem);
3993       if((fp = fopen(fname, "r")) == NULL)
3994         {
3995           MRIfree(&mri);
3996           errno = 0;
3997           ErrorReturn(NULL, (ERROR_BADFILE, "can't find file %s (last resort);"
3998                              "bailing out on read", fname));
3999         }
4000       fscanf(fp, "%d %d %d %*d", &ny, &nx, &nt);
4001       fclose(fp);
4002 
4003       /* --- get the number of slices --- */
4004       sprintf(fname, "%s/%s_000.%s", directory, stem, extension);
4005       if(!FileExists(fname)){
4006         MRIfree(&mri);
4007         errno = 0;
4008         ErrorReturn(NULL, (ERROR_BADFILE, "can't find file %s; "
4009                            "bailing out on read", fname));
4010       }
4011       for(nslices = 0;FileExists(fname);nslices++)
4012         sprintf(fname, "%s/%s_%03d.%s", directory, stem, nslices, extension);
4013       nslices--;
4014 
4015       mri->width   = nx;
4016       mri->height  = ny;
4017       mri->depth   = nslices;
4018       mri->nframes = nt;
4019 
4020       mri->imnr0 = 1;
4021       mri->imnr1 = nslices;
4022 
4023       mri->thick = mri->ps = 1.0;
4024       mri->xsize = mri->ysize = mri->zsize = 1.0;
4025 
4026       setDirectionCosine(mri, MRI_CORONAL);
4027 
4028       mri->ras_good_flag = 0;
4029 
4030       strcpy(mri->fname, fname_passed);
4031 
4032     }
4033 
4034   mri->imnr0 = 1;
4035   mri->imnr1 = nslices;
4036 
4037   mri->xend = mri->width  * mri->xsize / 2.0;
4038   mri->xstart = -mri->xend;
4039   mri->yend = mri->height * mri->ysize / 2.0;
4040   mri->ystart = -mri->yend;
4041   mri->zend = mri->depth  * mri->zsize / 2.0;
4042   mri->zstart = -mri->zend;
4043 
4044   mri->fov = ((mri->xend - mri->xstart) > (mri->yend - mri->ystart) ?
4045               (mri->xend - mri->xstart) : (mri->yend - mri->ystart));
4046 
4047 
4048   if(read_volume)
4049     {
4050       mri2 = MRIallocSequence(mri->width, mri->height, mri->depth,
4051                               mri->type, mri->nframes);
4052       MRIcopyHeader(mri, mri2);
4053       MRIfree(&mri);
4054       mri = mri2;
4055     }
4056 
4057   return(mri);
4058 
4059 } /* end get_b_info() */
4060 
4061 /*-------------------------------------------------------------------
4062   bvolumeRead() - this replaces bshortRead and bfloatRead.
4063   -------------------------------------------------------------------*/
4064 static MRI *bvolumeRead(char *fname_passed, int read_volume, int type)
4065 {
4066 
4067   MRI *mri;
4068   FILE *fp;
4069   char fname[STRLEN];
4070   char directory[STRLEN];
4071   char stem[STRLEN];
4072   int swap_bytes_flag;
4073   int slice, frame, row, k;
4074   int nread;
4075   char *ext;
4076   int size;
4077   float min, max;
4078 
4079   /* check the type and set the extension and size*/
4080   switch(type) {
4081   case MRI_SHORT: ext = "bshort"; size = sizeof(short); break;
4082   case MRI_FLOAT: ext = "bfloat"; size = sizeof(float); break;
4083   default:
4084     fprintf(stderr,"ERROR: bvolumeRead: type (%d) is not "
4085             "short or float\n", type);
4086     return(NULL);
4087   }
4088 
4089   /* Get the header info (also allocs if needed) */
4090   mri = get_b_info(fname_passed, read_volume, directory, stem, type);
4091   if(mri == NULL) return(NULL);
4092 
4093   /* If not reading the volume, return now */
4094   if(! read_volume) return(mri);
4095 
4096   /* Read in the header of the first slice to get the endianness */
4097   sprintf(fname, "%s/%s_%03d.hdr", directory, stem, 0);
4098   if((fp = fopen(fname, "r")) == NULL){
4099     fprintf(stderr, "ERROR: can't open file %s; assuming big-endian bvolume\n",
4100             fname);
4101     swap_bytes_flag = 0;
4102   }
4103   else{
4104     fscanf(fp, "%*d %*d %*d %d", &swap_bytes_flag);
4105 #if (BYTE_ORDER == LITTLE_ENDIAN)
4106     swap_bytes_flag = !swap_bytes_flag;
4107 #endif
4108     fclose(fp);
4109   }
4110 
4111   /* Go through each slice */
4112   for(slice = 0;slice < mri->depth; slice++){
4113 
4114     /* Open the file for this slice */
4115     sprintf(fname, "%s/%s_%03d.%s", directory, stem, slice, ext);
4116     if((fp = fopen(fname, "r")) == NULL){
4117       MRIfree(&mri);
4118       errno = 0;
4119       ErrorReturn(NULL, (ERROR_BADFILE,
4120                          "bvolumeRead(): error opening file %s", fname));
4121     }
4122     //fprintf(stderr, "Reading %s ... \n", fname);
4123     /* Loop through the frames */
4124     for(frame = 0; frame < mri->nframes; frame ++){
4125       k = slice + mri->depth*frame;
4126 
4127       /* Loop through the rows */
4128       for(row = 0;row < mri->height; row++){
4129 
4130         /* read in all the columns for a row */
4131         nread = fread(mri->slices[k][row], size, mri->width, fp);
4132         if( nread != mri->width){
4133           fclose(fp);
4134           MRIfree(&mri);
4135           errno = 0;
4136           ErrorReturn(NULL, (ERROR_BADFILE, "bvolumeRead(): "
4137                              "error reading from file %s", fname));
4138         }
4139 
4140         if(swap_bytes_flag){
4141           if(type == MRI_SHORT)
4142             swab(mri->slices[k][row], mri->slices[k][row], mri->width * size);
4143           else
4144             byteswapbuffloat((void*)mri->slices[k][row],size*mri->width);
4145         }
4146       } /* row loop */
4147     } /* frame loop */
4148 
4149     fclose(fp);
4150   } /* slice loop */
4151 
4152   MRIlimits(mri,&min,&max);
4153   printf("INFO: bvolumeRead: min = %g, max = %g\n",min,max);
4154 
4155   mri->imnr0 = 1;
4156   mri->imnr1 = mri->depth;
4157   mri->thick = mri->zsize;
4158 
4159   return(mri);
4160 
4161 } /* end bvolumeRead() */
4162 
4163 
4164 #if 0
4165 static int orient_with_register(MRI *mri)
4166 {
4167 
4168   MRI *subject_mri;
4169   char *subjects_dir;
4170   char subject_directory[STRLEN];
4171   MATRIX *sa, *fa;
4172   MATRIX *sr, *fr;
4173   MATRIX *rinv, *sainv;
4174   MATRIX *r1, *r2;
4175   int res;
4176   float det;
4177 
4178   subject_mri = NULL;
4179   if((subjects_dir = getenv("SUBJECTS_DIR")) != NULL)
4180     {
4181       sprintf
4182         (subject_directory, "%s/%s/mri/T1", subjects_dir, mri->subject_name);
4183       if((subject_mri = MRIreadInfo(subject_directory)) == NULL)
4184         {
4185           errno = 0;
4186           ErrorPrintf
4187             (ERROR_BADPARM,
4188              "can't get get information from %s; ", subject_directory);
4189           sprintf
4190             (subject_directory,
4191              "%s/%s/mri/T1", subjects_dir, mri->subject_name);
4192           ErrorPrintf
4193             (ERROR_BADPARM, "trying %s instead...\n", subject_directory);
4194           if((subject_mri = MRIreadInfo(subject_directory)) == NULL)
4195             ErrorPrintf
4196               (ERROR_BADPARM,
4197                "can't get information from %s\n", subject_directory);
4198         }
4199     }
4200   else
4201     {
4202       errno = 0;
4203       ErrorPrintf
4204         (ERROR_BADPARM, "can't get environment variable SUBJECTS_DIR");
4205     }
4206 
4207   if(subject_mri == NULL)
4208     {
4209 
4210       errno = 0;
4211       ErrorPrintf(ERROR_BADPARM, "guessing at COR- orientation...\n");
4212 
4213       subject_mri = MRIallocHeader(256, 256, 256, MRI_UCHAR);
4214       subject_mri->fov = 256.0;
4215       subject_mri->thick = subject_mri->ps = 1.0;
4216       subject_mri->xsize = subject_mri->ysize = subject_mri->zsize = 1.0;
4217       setDirectionCosine(subject_mri, MRI_CORONAL);
4218       subject_mri->ras_good_flag = 1;
4219     }
4220 
4221   det = MatrixDeterminant(mri->register_mat);
4222   if(det == 0.0)
4223     {
4224       MRIfree(&subject_mri);
4225       errno = 0;
4226       ErrorPrintf
4227         (ERROR_BADPARM,
4228          "orient_with_register(): registration matrix has zero determinant");
4229       ErrorPrintf(ERROR_BADPARM, "matrix is:");
4230       MatrixPrint(stderr, mri->register_mat);
4231       return(ERROR_BADPARM);
4232     }
4233 
4234   rinv = MatrixInverse(mri->register_mat, NULL);
4235 
4236   if(rinv == NULL)
4237     {
4238       MRIfree(&subject_mri);
4239       errno = 0;
4240       ErrorPrintf
4241         (ERROR_BADPARM,
4242          "orient_with_register(): error inverting registration matrix");
4243       ErrorPrintf(ERROR_BADPARM, "matrix is:");
4244       MatrixPrint(stderr, mri->register_mat);
4245       return(ERROR_BADPARM);
4246     }
4247 
4248   sr = extract_i_to_r(subject_mri);
4249 
4250   sa = MatrixAlloc(4, 4, MATRIX_REAL);
4251   fa = MatrixAlloc(4, 4, MATRIX_REAL);
4252 
4253   if(sr == NULL || sa == NULL || fa == NULL)
4254     {
4255       if(sr != NULL)
4256         MatrixFree(&sr);
4257       if(sa != NULL)
4258         MatrixFree(&sa);
4259       if(fa != NULL)
4260         MatrixFree(&fa);
4261       MatrixFree(&rinv);
4262       MRIfree(&subject_mri);
4263       errno = 0;
4264       ErrorReturn
4265