00001 #include "mrilib.h"
00002 #include <signal.h>
00003 
00004 /*******************************************************************/
00005 /*!  Write an AFNI dataset as a MINC file.
00006      - fname = MINC filename
00007      - dset = AFNI dataset
00008      - flags = logical OR of various masks:
00009      - MINC_FLOATIZE_MASK => force output in float format
00010 
00011    Return value is 1 if went OK, 0 if not.
00012 -------------------------------------------------------------------*/
00013 
00014 int THD_write_minc( char *fname, THD_3dim_dataset *dset, int flags )
00015 {
00016    static char *pg=NULL ; static int first=1 ;
00017 
00018    int nvals,nvox , iv , floatize , ii , datum , good=1 ;
00019    THD_dataxes *dax ;
00020    char  axcode[3] ;
00021    float axstep[3] , axstart[3] ;
00022    int   axnum[3] ;
00023    char *cmd ;
00024    FILE *fp ;
00025    void *bar ;
00026    float *far=NULL ;
00027 
00028 ENTRY("THD_write_minc") ;
00029 
00030    /*-- check inputs for goodness --*/
00031 
00032    if( !THD_filename_ok(fname) || fname[0] == '-' ){
00033       fprintf(stderr,"** ERROR: Illegal filename for MINC output: %s\n",
00034               (fname != NULL) ? fname : "(null)" ) ;
00035       RETURN(0) ;
00036    }
00037 
00038    if( !ISVALID_DSET(dset) ){
00039       fprintf(stderr,
00040               "** ERROR: Illegal input dataset for MINC output: %s\n",
00041               fname ) ;
00042       RETURN(0) ;
00043    }
00044 
00045    /*-- find rawtominc program to do the real work --*/
00046 
00047    if( first ){
00048      pg = THD_find_executable("rawtominc") ; first = 0 ;
00049    }
00050 
00051    if( pg == NULL ){
00052      fprintf(stderr,
00053              "** ERROR: Can't write MINC file without program rawtominc: %s\n",
00054              fname);
00055      RETURN(0) ;
00056    }
00057 
00058    /*-- load dataset from disk, if need be --*/
00059 
00060    DSET_load(dset) ;
00061    if( !DSET_LOADED(dset) ){
00062      fprintf(stderr,
00063              "** ERROR: Can't write MINC file since dataset isn't loaded: %s\n",
00064              fname) ;
00065      RETURN(0) ;
00066    }
00067 
00068    /*-- check sub-bricks for uniformity in type, legal type, etc. --*/
00069 
00070    nvals    = DSET_NVALS(dset) ;
00071    nvox     = DSET_NVOX(dset) ;
00072    floatize = (flags & MINC_FLOATIZE_MASK) != 0 ;
00073    ii       = DSET_BRICK_TYPE(dset,0) ;
00074 
00075    for( iv=0 ; iv < nvals ; iv++ ){
00076      if( DSET_BRICK_TYPE(dset,iv) == MRI_complex ){
00077        fprintf(stderr,
00078                "** ERROR: Can't write MINC file since dataset is complex: %s\n",
00079                fname) ;
00080        RETURN(0) ;
00081      }
00082      if( DSET_BRICK_TYPE(dset,iv) == MRI_rgb ){
00083        fprintf(stderr,
00084                "** ERROR: Can't write MINC file since dataset is RGB: %s\n",
00085                fname) ;
00086        RETURN(0) ;
00087      }
00088 
00089      /* check if must convert to floats */
00090 
00091      if( !floatize ){
00092         if( DSET_BRICK_TYPE(dset,iv)   == MRI_float ||
00093             DSET_BRICK_FACTOR(dset,iv) != 0.0       ||
00094             DSET_BRICK_TYPE(dset,iv)   != ii          ){
00095 
00096           floatize = 1 ;
00097           fprintf(stderr,
00098                   "++ NOTE: Writing MINC file in float format: %s\n",
00099                   fname) ;
00100         }
00101      }
00102    }
00103 
00104    if( floatize ) datum = MRI_float ;
00105    else           datum = ii ;
00106 
00107    /*-- load geometry information --*/
00108 
00109    dax = dset->daxes ;
00110 
00111    axcode[0] = ORIENT_xyz[ dax->xxorient ] ; axnum[0] = dax->nxx ;
00112    axcode[1] = ORIENT_xyz[ dax->yyorient ] ; axnum[1] = dax->nyy ;
00113    axcode[2] = ORIENT_xyz[ dax->zzorient ] ; axnum[2] = dax->nzz ;
00114 
00115    axstep[0] = dax->xxdel ; axstart[0] = dax->xxorg ;
00116    axstep[1] = dax->yydel ; axstart[1] = dax->yyorg ;
00117    axstep[2] = dax->zzdel ; axstart[2] = dax->zzorg ;
00118 
00119    if( axcode[0] == 'x' || axcode[0] == 'y' ){
00120      axstep[0] = -axstep[0] ; axstart[0] = -axstart[0] ;
00121    }
00122 
00123    if( axcode[1] == 'x' || axcode[1] == 'y' ){
00124      axstep[1] = -axstep[1] ; axstart[1] = -axstart[1] ;
00125    }
00126 
00127    if( axcode[2] == 'x' || axcode[2] == 'y' ){
00128      axstep[2] = -axstep[2] ; axstart[2] = -axstart[2] ;
00129    }
00130 
00131    /*-- start to create command --*/
00132 
00133    cmd = AFMALL(char, 65500) ; /* long enough?  */
00134    strcpy(cmd,pg) ;      /* basic command */
00135 
00136    /* axes orientation */
00137 
00138    sprintf(cmd+strlen(cmd)," -%c%c%c",axcode[2],axcode[1],axcode[0]) ;
00139 
00140    /* input and output data type */
00141 
00142    sprintf(cmd+strlen(cmd)," -%s -o%s",
00143            MRI_TYPE_name[datum],MRI_TYPE_name[datum]) ;
00144 
00145    /* axis stuff */
00146 
00147    sprintf(cmd+strlen(cmd),
00148            " -%cstep %.3f -%cstep %.3f -%cstep %.3f"
00149            " -%cstart %.3f -%cstart %.3f -%cstart %.3f"
00150            " -xdircos 1 0 0 -ydircos 0 1 0 -zdircos 0 0 1" ,
00151            axcode[0],axstep[0] ,axcode[1],axstep[1] ,axcode[2],axstep[2]  ,
00152            axcode[0],axstart[0],axcode[1],axstart[1],axcode[2],axstart[2]  ) ;
00153 
00154    /*-- do we create a time step attribute? --*/
00155 
00156    if( nvals > 1 & DSET_NUM_TIMES(dset) > 1 && DSET_TIMESTEP(dset) > 0.0 ){
00157      float tr = DSET_TIMESTEP(dset) ;
00158      sprintf(cmd+strlen(cmd)," -dattribute time:step=%.3f",tr) ;
00159      sprintf(cmd+strlen(cmd)," -frame_times '") ;
00160      for( iv=0 ; iv < nvals ; iv++ )
00161        sprintf(cmd+strlen(cmd),"%.3f ",iv*tr) ;
00162      strcat(cmd,"'") ;
00163    }
00164 
00165    /*-- Is this Talaraich? --*/
00166 
00167    if( dset->view_type == VIEW_TALAIRACH_TYPE ){
00168      sprintf(cmd+strlen(cmd)," -sattribute xspace:spacetype=talairach_") ;
00169      sprintf(cmd+strlen(cmd)," -sattribute yspace:spacetype=talairach_") ;
00170      sprintf(cmd+strlen(cmd)," -sattribute zspace:spacetype=talairach_") ;
00171    }
00172 
00173    /*-- integer datasets need to be scanned    --*/
00174    /*-- and have valid_range set (stupid MINC) --*/
00175 
00176    if( !floatize ){
00177      sprintf(cmd+strlen(cmd)," -scan_range") ;
00178      switch( datum ){
00179         case MRI_short:
00180           sprintf(cmd+strlen(cmd)," -range 0 32767") ; break ;
00181 
00182         case MRI_byte:
00183           sprintf(cmd+strlen(cmd)," -range 0 255")   ; break ;
00184      }
00185    }
00186 
00187    /*-- add output file name --*/
00188 
00189    sprintf(cmd+strlen(cmd)," %s",fname) ;
00190 
00191    /*-- add # sub-bricks, if > 1 --*/
00192 
00193    if( nvals > 1 ) sprintf(cmd+strlen(cmd)," %d",nvals) ;
00194 
00195    /*-- add number of points along each axis --*/
00196 
00197    sprintf(cmd+strlen(cmd)," %d %d %d" , axnum[2],axnum[1],axnum[0] ) ;
00198 
00199    /*-- execute the command as a filter --*/
00200 
00201    signal( SIGPIPE , SIG_IGN ) ; errno = 0 ;
00202    fp = popen( cmd , "w" ) ;
00203    if( fp == NULL ){
00204      fprintf(stderr,"** ERROR: Can't open MINC output filter: %s\a\n",cmd) ;
00205      if( errno != 0 ) perror("** Unix error message") ;
00206      free(cmd) ; RETURN(0) ;
00207    }
00208 
00209    /*-- allocate space for floatizing --*/
00210 
00211    if( floatize ){
00212      far = (float *) malloc(sizeof(float)*nvox) ;
00213      if( far == NULL ){
00214        fprintf(stderr,
00215                "** ERROR: Can't write MINC file due to lack of memory: %s\n",
00216                fname) ;
00217        free(cmd) ; RETURN(0) ;
00218      }
00219    }
00220 
00221    /*-- loop through sub-bricks, convert to floats if needed, write to pipe --*/
00222 
00223    for( iv=0 ; iv < nvals ; iv++ ){
00224       bar = DSET_ARRAY(dset,iv) ; errno = 0 ;
00225       if( floatize ){
00226         EDIT_coerce_scale_type( nvox ,
00227                                 DSET_BRICK_FACTOR(dset,iv) ,
00228                                 DSET_BRICK_TYPE(dset,iv)   , bar ,
00229                                 MRI_float                  , far  ) ;
00230         ii = fwrite( far , sizeof(float)         , nvox , fp ) ;
00231       } else {
00232         ii = fwrite( bar , mri_datum_size(datum) , nvox , fp ) ;
00233       }
00234 
00235       DSET_unload_one(dset,iv) ;
00236 
00237       if( ii < nvox ){
00238         fprintf(stderr,
00239                 "** ERROR: fwrite to MINC failed at iv=%d: %s\n",
00240                 iv , fname ) ;
00241         if( errno ) perror("fwrite") ;
00242         if( ii == 0 ){ good=0; break; }       /* don't try to continue */
00243       }
00244    }
00245 
00246    if( floatize ) free(far) ;  /* no longer needed */
00247 
00248    /*-- close pipe --*/
00249 
00250    ii = pclose(fp) ;
00251    if( ii == -1 ){
00252      perror("** ERROR in MINC write filter") ;
00253      fprintf(stderr,"** MINC filter command was %s\n\a",cmd) ;
00254      free(cmd) ; RETURN(0) ;
00255    }
00256 
00257    /*-- and we are done --*/
00258 
00259    free(cmd) ; RETURN(good) ;
00260 }