00137 {
00138    double sumxx , sumyy , sumxy , tx,ty , dxy ;
00139    void  *  xar , *  yar ;
00140    float * fxar , * fyar ;
00141    int ii , nxyz , ivx,ivy , itypx,itypy , fxar_new,fyar_new , nnn ;
00142 
00143    nxyz = DSET_NVOX(xset) ;
00144 
00145    /* load bricks */
00146 
00147    THD_load_datablock( xset->dblk ) ;
00148    ivx   = 0 ;
00149    itypx = DSET_BRICK_TYPE(xset,ivx) ;
00150    xar   = DSET_ARRAY(xset,ivx) ; if( xar == NULL ) return 0.0 ;
00151    if( itypx == MRI_float ){
00152       fxar = (float *) xar ; fxar_new = 0 ;
00153    } else {
00154       fxar = (float *) malloc( sizeof(float) * nxyz ) ; fxar_new = 1 ;
00155       EDIT_coerce_type( nxyz , itypx,xar , MRI_float,fxar ) ;
00156       PURGE_DSET( xset ) ;
00157    }
00158 
00159    THD_load_datablock( yset->dblk ) ;
00160    ivy   = 0 ;
00161    itypy = DSET_BRICK_TYPE(yset,ivy) ;
00162    yar   = DSET_ARRAY(yset,ivy) ; if( yar == NULL ) return 0.0 ;
00163    if( itypy == MRI_float ){
00164       fyar = (float *) yar ; fyar_new = 0 ;
00165    } else {
00166       fyar = (float *) malloc( sizeof(float) * nxyz ) ; fyar_new = 1 ;
00167       EDIT_coerce_type( nxyz , itypy,yar , MRI_float,fyar ) ;
00168       PURGE_DSET( yset ) ;
00169    }
00170 
00171    /* 29 Feb 2000: remove mean? */
00172 
00173    if( dm ){
00174       sumxx = sumyy = 0.0 ;
00175       for( nnn=ii=0 ; ii < nxyz ; ii++ ){
00176          if( mmm == NULL || mmm[ii] ){
00177             sumxx += fxar[ii] ; sumyy += fyar[ii] ; nnn++ ;
00178          }
00179       }
00180       if( nnn < 5 ) return 0.0 ;             /* ERROR */
00181       sumxx /= nnn ; sumyy /= nnn ;
00182       for( ii=0 ; ii < nxyz ; ii++ ){
00183          if( mmm == NULL || mmm[ii] ){
00184             fxar[ii] -= sumxx ; fyar[ii] -= sumyy ;
00185          }
00186       }
00187    }
00188 
00189    /* compute sums */
00190 
00191    sumxx = sumyy = sumxy = 0.0 ;
00192 
00193    for( ii=0 ; ii < nxyz ; ii++ ){
00194       if( mmm == NULL || mmm[ii] ){
00195          tx = fxar[ii] ; ty = fyar[ii] ;
00196          sumxx += tx * tx ;
00197          sumyy += ty * ty ;
00198          sumxy += tx * ty ;
00199       }
00200    }
00201 
00202    /* toss trash */
00203 
00204    if( fxar_new ) free( fxar ) ;
00205    if( fyar_new ) free( fyar ) ;
00206 
00207    /* compute result */
00208 
00209    dxy = sumxx * sumyy ; if( dxy <= 0.0 ) return 0.0 ;
00210    dxy = sumxy / sqrt(dxy) ;
00211    return dxy ;
00212 }

00013 {
00014    double dxy ;
00015    int narg , ndset , nvox , demean=0 ;
00016    THD_3dim_dataset * xset , * yset , * mask_dset=NULL ;
00017    float mask_bot=666.0 , mask_top=-666.0 ;
00018    byte * mmm=NULL ;
00019 
00020    /*-- read command line arguments --*/
00021 
00022    if( argc < 3 || strncmp(argv[1],"-help",5) == 0 ){
00023       printf("Usage: 3ddot [options] dset1 dset2\n"
00024              "Output = correlation coefficient between 2 dataset bricks\n"
00025              "         - you can use sub-brick selectors on the dsets\n"
00026              "         - the result is a number printed to stdout"
00027              "Options:\n"
00028              "  -mask mset   Means to use the dataset 'mset' as a mask:\n"
00029              "                 Only voxels with nonzero values in 'mset'\n"
00030              "                 will be averaged from 'dataset'.  Note\n"
00031              "                 that the mask dataset and the input dataset\n"
00032              "                 must have the same number of voxels.\n"
00033              "  -mrange a b  Means to further restrict the voxels from\n"
00034              "                 'mset' so that only those mask values\n"
00035              "                 between 'a' and 'b' (inclusive) will\n"
00036              "                 be used.  If this option is not given,\n"
00037              "                 all nonzero values from 'mset' are used.\n"
00038              "                 Note that if a voxel is zero in 'mset', then\n"
00039              "                 it won't be included, even if a < 0 < b.\n"
00040              "  -demean      Means to remove the mean from each volume\n"
00041              "                 prior to computing the correlation.\n"
00042             ) ;
00043 
00044       printf("\n" MASTER_SHORTHELP_STRING ) ;
00045 
00046       exit(0) ;
00047    }
00048 
00049    narg = 1 ;
00050    while( narg < argc && argv[narg][0] == '-' ){
00051 
00052       if( strncmp(argv[narg],"-demean",5) == 0 ){
00053          demean++ ;
00054          narg++ ; continue ;
00055       }
00056 
00057       if( strncmp(argv[narg],"-mask",5) == 0 ){
00058          if( mask_dset != NULL ){
00059             fprintf(stderr,"*** Cannot have two -mask options!\n") ; exit(1) ;
00060          }
00061          if( narg+1 >= argc ){
00062             fprintf(stderr,"*** -mask option requires a following argument!\n");
00063              exit(1) ;
00064          }
00065          mask_dset = THD_open_dataset( argv[++narg] ) ;
00066          if( mask_dset == NULL ){
00067             fprintf(stderr,"*** Cannot open mask dataset!\n") ; exit(1) ;
00068          }
00069          if( DSET_BRICK_TYPE(mask_dset,0) == MRI_complex ){
00070             fprintf(stderr,"*** Cannot deal with complex-valued mask dataset!\n");
00071             exit(1) ;
00072          }
00073          narg++ ; continue ;
00074       }
00075 
00076       if( strncmp(argv[narg],"-mrange",5) == 0 ){
00077          if( narg+2 >= argc ){
00078             fprintf(stderr,"*** -mrange option requires 2 following arguments!\n")
00079 ;
00080              exit(1) ;
00081          }
00082          mask_bot = strtod( argv[++narg] , NULL ) ;
00083          mask_top = strtod( argv[++narg] , NULL ) ;
00084          if( mask_top < mask_top ){
00085             fprintf(stderr,"*** -mrange inputs are illegal!\n") ; exit(1) ;
00086          }
00087          narg++ ; continue ;
00088       }
00089 
00090       fprintf(stderr,"*** Unknown option: %s\n",argv[narg]) ; exit(1) ;
00091    }
00092 
00093    /* should have at least 2 more arguments */
00094 
00095    ndset = argc - narg ;
00096    if( ndset <= 1 ){
00097       fprintf(stderr,"*** No input datasets!?\n") ; exit(1) ;
00098    }
00099 
00100    xset = THD_open_dataset( argv[narg++] ) ;
00101    yset = THD_open_dataset( argv[narg++] ) ;
00102    if( xset == NULL || yset == NULL ){
00103       fprintf(stderr,"*** cannot open both input datasets!\n") ; exit(1) ;
00104    }
00105    if( DSET_NUM_TIMES(xset) > 1 || DSET_NUM_TIMES(yset) > 1 ){
00106       fprintf(stderr,"*** cannot use time-dependent datasets!\n") ; exit(1) ;
00107    }
00108    nvox = DSET_NVOX(xset) ;
00109    if( nvox != DSET_NVOX(yset) ){
00110       fprintf(stderr,"*** input datasets dimensions don't match!\n");exit(1);
00111    }
00112 
00113    /* make a byte mask from mask dataset */
00114 
00115    if( mask_dset != NULL ){
00116       int mcount ;
00117       if( DSET_NVOX(mask_dset) != nvox ){
00118          fprintf(stderr,"*** Input and mask datasets are not same dimensions!\n");
00119          exit(1) ;
00120       }
00121       mmm = THD_makemask( mask_dset , 0 , mask_bot,mask_top ) ;
00122       mcount = THD_countmask( nvox , mmm ) ;
00123       fprintf(stderr,"+++ %d voxels in the mask\n",mcount) ;
00124       if( mcount <= 5 ){
00125          fprintf(stderr,"*** Mask is too small!\n");exit(1);
00126       }
00127       DSET_delete(mask_dset) ;
00128    }
00129 
00130    dxy = DSET_cor( xset , yset , mmm , demean ) ;
00131    printf("%g\n",dxy) ;
00132    exit(0) ;
00133 }