00008 {
00009    THD_3dim_dataset *xset , *cset ;
00010    int nopt=1 , method=PEARSON , do_autoclip=0 ;
00011    int nvox , nvals , ii,jj,kout , polort=1 , ix,jy,kz ;
00012    MRI_IMAGE *xsim , *ysim ;
00013    float     *xsar , *ysar ;
00014    short     *car ;
00015    char * prefix = "ATcorr" ;
00016    byte * mmm=NULL ;
00017    int   nmask , abuc=1 ;
00018    char str[32] ;
00019 
00020    /*----*/
00021 
00022    if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00023       printf("Usage: 3dAutoTcorrelate [options] dset\n"
00024              "Computes the correlation coefficient between each pair of\n"
00025              "voxels in the input dataset, and stores the output into\n"
00026              "a new anatomical bucket dataset.\n"
00027              "\n"
00028              "Options:\n"
00029              "  -pearson  = Correlation is the normal Pearson (product moment)\n"
00030              "                correlation coefficient [default].\n"
00031              "  -spearman = Correlation is the Spearman (rank) correlation\n"
00032              "                coefficient.\n"
00033              "  -quadrant = Correlation is the quadrant correlation coefficient.\n"
00034              "\n"
00035              "  -polort m = Remove polynomical trend of order 'm', for m=-1..3.\n"
00036              "                [default is m=1; removal is by least squares].\n"
00037              "                Using m=-1 means no detrending; this is only useful\n"
00038              "                for data/information that has been pre-processed.\n"
00039              "\n"
00040              "  -autoclip = Clip off low-intensity regions in the dataset,\n"
00041              "  -automask =  so that the correlation is only computed between\n"
00042              "               high-intensity (presumably brain) voxels.  The\n"
00043              "               intensity level is determined the same way that\n"
00044              "               3dClipLevel works.\n"
00045              "\n"
00046              "  -prefix p = Save output into dataset with prefix 'p'\n"
00047              "               [default prefix is 'ATcorr'].\n"
00048              "\n"
00049              "  -time     = Save output as a 3D+time dataset instead\n"
00050              "               of a anat bucket.\n"
00051              "\n"
00052              "Notes:\n"
00053              " * The output dataset is anatomical bucket type of shorts.\n"
00054              " * The output file might be gigantic and you might run out\n"
00055              "    of memory running this program.  Use at your own risk!\n"
00056              " * The program prints out an estimate of its memory usage\n"
00057              "    when it starts.  It also prints out a progress 'meter'\n"
00058              "    of 1 dot per 10 output sub-bricks.\n"
00059              " * This is a quick hack for Peter Bandettini. Now pay up.\n"
00060              "\n"
00061              "-- RWCox - Jan 31 2002\n"
00062             ) ;
00063       exit(0) ;
00064    }
00065 
00066    mainENTRY("3dAutoTcorrelate main"); machdep(); PRINT_VERSION("3dAutoTcorrelate");
00067    AFNI_logger("3dAutoTcorrelate",argc,argv);
00068 
00069    /*-- option processing --*/
00070 
00071    while( nopt < argc && argv[nopt][0] == '-' ){
00072 
00073       if( strcmp(argv[nopt],"-time") == 0 ){
00074          abuc = 0 ; nopt++ ; continue ;
00075       }
00076 
00077       if( strcmp(argv[nopt],"-autoclip") == 0 ||
00078           strcmp(argv[nopt],"-automask") == 0   ){
00079 
00080          do_autoclip = 1 ; nopt++ ; continue ;
00081       }
00082 
00083       if( strcmp(argv[nopt],"-pearson") == 0 ){
00084          method = PEARSON ; nopt++ ; continue ;
00085       }
00086 
00087       if( strcmp(argv[nopt],"-spearman") == 0 ){
00088          method = SPEARMAN ; nopt++ ; continue ;
00089       }
00090 
00091       if( strcmp(argv[nopt],"-quadrant") == 0 ){
00092          method = QUADRANT ; nopt++ ; continue ;
00093       }
00094 
00095       if( strcmp(argv[nopt],"-prefix") == 0 ){
00096          prefix = argv[++nopt] ;
00097          if( !THD_filename_ok(prefix) ){
00098             fprintf(stderr,"** Illegal value after -prefix!\n");exit(1);
00099          }
00100          nopt++ ; continue ;
00101       }
00102 
00103       if( strcmp(argv[nopt],"-polort") == 0 ){
00104          char *cpt ;
00105          int val = strtod(argv[++nopt],&cpt) ;
00106          if( *cpt != '\0' || val < -1 || val > 3 ){
00107             fprintf(stderr,"** Illegal value after -polort!\n");exit(1);
00108          }
00109          polort = val ; nopt++ ; continue ;
00110       }
00111 
00112       fprintf(stderr,"** Illegal option: %s\n",argv[nopt]) ; exit(1) ;
00113    }
00114 
00115    /*-- open dataset, check for legality --*/
00116 
00117    if( nopt >= argc ){
00118       fprintf(stderr,"*** Need a dataset on command line!?\n"); exit(1);
00119    }
00120 
00121    xset = THD_open_dataset( argv[nopt] ) ;
00122    if( xset == NULL ){
00123       fprintf(stderr,"** Can't open dataset %s\n",argv[nopt]); exit(1);
00124    }
00125    if( DSET_NUM_TIMES(xset) < 2 ){
00126       fprintf(stderr,"** Input dataset %s is not 3D+time\n",argv[nopt]); exit(1);
00127    }
00128    DSET_load(xset) ;
00129    if( !DSET_LOADED(xset) ){
00130       fprintf(stderr,"*** Can't read dataset bricks from %s\n",argv[nopt-1]);
00131       exit(1);
00132    }
00133 
00134    /*-- compute mask array, if desired --*/
00135 
00136    nvox = DSET_NVOX(xset) ; nvals = DSET_NVALS(xset) ;
00137 
00138    if( do_autoclip ){
00139       mmm   = THD_automask( xset ) ;
00140       nmask = THD_countmask( nvox , mmm ) ;
00141       fprintf(stderr,"++ %d voxels survive -autoclip\n",nmask) ;
00142       if( nmask < 2 ) exit(1) ;
00143    } else {
00144       nmask = nvox ;
00145       fprintf(stderr,"++ computing for all %d voxels\n",nmask) ;
00146    }
00147 
00148    /*-- create output dataset --*/
00149 
00150    cset = EDIT_empty_copy( xset ) ;
00151 
00152    if( abuc ){
00153      EDIT_dset_items( cset ,
00154                         ADN_prefix    , prefix         ,
00155                         ADN_nvals     , nmask          ,
00156                         ADN_ntt       , 0              , /* no time axis */
00157                         ADN_type      , HEAD_ANAT_TYPE ,
00158                         ADN_func_type , ANAT_BUCK_TYPE ,
00159                       ADN_none ) ;
00160    } else {
00161      EDIT_dset_items( cset ,
00162                         ADN_prefix    , prefix         ,
00163                         ADN_nvals     , nmask          ,
00164                         ADN_ntt       , nmask          ,  /* num times */
00165                         ADN_ttdel     , 1.0            ,  /* fake TR */
00166                         ADN_nsl       , 0              ,  /* no slice offsets */
00167                         ADN_type      , HEAD_ANAT_TYPE ,
00168                         ADN_func_type , ANAT_EPI_TYPE  ,
00169                       ADN_none ) ;
00170    }
00171 
00172    if( THD_is_file(DSET_HEADNAME(cset)) ){
00173       fprintf(stderr,"** Output dataset %s already exists!\n",
00174               DSET_HEADNAME(cset)) ;
00175       exit(1) ;
00176    }
00177 
00178    { double nb = (double)(xset->dblk->total_bytes) ;
00179      nb += (double)(nmask) * (double)(nvox) * sizeof(short) ;
00180      nb /= (1024.0*1024.0) ;
00181      fprintf(stderr,
00182        "++ Memory required = %.1f Mbytes for %d output sub-bricks\n",nb,nmask);
00183      if( nb > 2000.0 ){
00184        fprintf(stderr,"** ERROR: can't use more than 2000 Mbytes of memory!\n");
00185        exit(1) ;
00186      }
00187    }
00188 
00189    tross_Make_History( "3dAutoTcorrelate" , argc,argv , cset ) ;
00190 
00191    /* loop over voxels, correlate */
00192 
00193    kout = 0 ;
00194    for( ii=0 ; ii < nvox ; ii++ ){  /* voxel to correlate with */
00195 
00196       if( mmm != NULL && mmm[ii] == 0 ) continue ; /* skip it */
00197 
00198       if( kout > 0 && kout%10 == 0 )
00199          fprintf(stderr, (kout%1000==0) ? "!"
00200                         :(kout%100 ==0) ? ":" : "." ) ;  /* progress meter */
00201 
00202       /* create and modify output brick */
00203 
00204       EDIT_substitute_brick( cset,kout, MRI_short, NULL ) ; /* make array */
00205       car = DSET_ARRAY(cset,kout) ;                         /* get array  */
00206 
00207       EDIT_BRICK_TO_FICO(cset,kout,nvals,1,polort+1) ;  /* stat params  */
00208       EDIT_BRICK_FACTOR(cset,kout,0.0001) ;             /* scale factor */
00209 
00210       ix = DSET_index_to_ix(cset,ii) ;                  /* brick label  */
00211       jy = DSET_index_to_jy(cset,ii) ;
00212       kz = DSET_index_to_kz(cset,ii) ;
00213       sprintf(str,"%03d_%03d_%03d",ix,jy,kz) ;
00214       EDIT_BRICK_LABEL(cset,kout,str) ;
00215 
00216       /* get ref time series from this voxel */
00217 
00218       xsim = THD_extract_series(ii,xset,0) ; xsar = MRI_FLOAT_PTR(xsim) ;
00219 
00220       DETREND_polort(polort,nvals,xsar) ;  /* remove polynomial trend */
00221 
00222       for( jj=0 ; jj < nvox ; jj++ ){  /* loop over voxels, correlate w/ref */
00223 
00224          if( mmm != NULL && mmm[jj] == 0 ){  /* the easy case */
00225             car[jj] = 0 ; continue ;
00226          }
00227 
00228          ysim = THD_extract_series(jj,xset,0) ; ysar = MRI_FLOAT_PTR(ysim) ;
00229          DETREND_polort(polort,nvals,ysar) ;
00230 
00231          switch( method ){                    /* correlate */
00232             default:
00233             case PEARSON:
00234               car[jj] = rint(10000.0*THD_pearson_corr (nvals,xsar,ysar));
00235             break;
00236             case SPEARMAN:
00237               car[jj] = rint(10000.0*THD_spearman_corr(nvals,xsar,ysar));
00238             break;
00239             case QUADRANT:
00240               car[jj] = rint(10000.0*THD_quadrant_corr(nvals,xsar,ysar));
00241             break;
00242          }
00243 
00244          mri_free(ysim) ;
00245 
00246       } /* end of loop over voxels in this sub-brick */
00247 
00248       mri_free(xsim) ;
00249       kout++ ;         /* move to next output brick */
00250 
00251    } /* end of loop over ref voxels */
00252 
00253    /* toss the other trash */
00254 
00255    DSET_unload(xset); if( mmm != NULL ) free(mmm);
00256 
00257    /* finito */
00258 
00259 
00260    DSET_write(cset) ;
00261    fprintf(stderr,"++ Wrote output: %s\n",DSET_BRIKNAME(cset)) ;
00262    exit(0) ;
00263 }