00590 {
00591    PLUGIN_interface * plint ;
00592 
00593    /*-- set titles and call point --*/
00594 
00595    plint = PLUTO_new_interface( "Histogram: CC" ,
00596                                 "Histogram of Correlation Coefficient" ,
00597                                 c_helpstring ,
00598                                 PLUGIN_CALL_VIA_MENU , CORREL_main  ) ;
00599 
00600    PLUTO_add_hint( plint , "Histogram: Correlation Coefficient" ) ;
00601 
00602    PLUTO_set_sequence( plint , "A:afniinfo:dset" ) ;
00603 
00604    /*-- first line of input --*/
00605 
00606    PLUTO_add_option( plint , "Source" , "Source" , TRUE ) ;
00607 
00608    PLUTO_add_dataset(  plint ,
00609                        "Dataset" ,        /* label next to button   */
00610                        ANAT_ALL_MASK ,    /* take any anat datasets */
00611                        FUNC_FIM_MASK ,    /* only allow fim funcs   */
00612                        DIMEN_4D_MASK |    /* need 3D+time datasets  */
00613                        BRICK_ALLREAL_MASK /* need real-valued datasets */
00614                     ) ;
00615 
00616    PLUTO_add_number( plint ,
00617                      "Ignore" ,   /* label next to chooser */
00618                      0 ,          /* smallest possible value */
00619                      999 ,        /* largest possible value */
00620                      0 ,          /* decimal shift (none in this case) */
00621                      4 ,          /* default value */
00622                      TRUE         /* allow user to edit value? */
00623                    ) ;
00624 
00625    /*-- second line of input --*/
00626 
00627    PLUTO_add_option( plint , "Vector" , "Vector" , TRUE ) ;
00628 
00629    PLUTO_add_timeseries( plint , "1D File" ) ;
00630 
00631    PLUTO_add_number( plint , "Polort" , 0,MAX_POLORT,0,1,FALSE ) ;
00632 
00633    /*-- (optional) third line of input --*/
00634 
00635    PLUTO_add_option( plint , "Mask" , "Mask" , FALSE ) ;
00636    PLUTO_add_dataset( plint , "Dataset" ,
00637                                     ANAT_ALL_MASK , FUNC_ALL_MASK ,
00638                                     DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00639    PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
00640 
00641    /*-- (optional) fourth line of input --*/
00642 
00643    PLUTO_add_option( plint , "Range"  , "Range" , FALSE ) ;
00644    PLUTO_add_number( plint , "Bottom" , -99999,99999, 1, 0,1 ) ;
00645    PLUTO_add_number( plint , "Top"    , -99999,99999,-1, 0,1 ) ;
00646 
00647 
00648    return plint ;
00649 }

00656 {
00657    MCW_idcode * idc ;
00658    THD_3dim_dataset * input_dset , * mask_dset = NULL ;
00659    MRI_IMAGE * tsim , * flim ;
00660    int ignore , nvox , ntim , polort , miv , it , ip , nupdt , nbin ;
00661    int mcount , ii , jj ;
00662    float * tsar ;
00663    float mask_bot=666.0 , mask_top=-666.0 , hbot,htop ;
00664    byte * mmm ;
00665    char buf[THD_MAX_NAME+16] , * tag ;
00666 
00667    PCOR_references * pc_ref ;
00668    PCOR_voxel_corr * pc_vc ;
00669    int fim_nref ;
00670    float * ref_vec , * vval , * zval , * aval ;
00671    int   * hbin , * jbin , * kbin , *jist[2] ;
00672    float sum , sumq , dbin , gval,rval,gg , sqp , zmid,zmed,zsig ;
00673    float pstar,zstar,zplus,zminus,psum,msum ;
00674 
00675    /*--------------------------------------------------------------------*/
00676    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00677 
00678    if( plint == NULL )
00679       return "************************\n"
00680              "CORREL_main:  NULL input\n"
00681              "************************"  ;
00682 
00683    /*-- read 1st line --*/
00684 
00685    PLUTO_next_option(plint) ;
00686    idc        = PLUTO_get_idcode(plint) ;
00687    input_dset = PLUTO_find_dset(idc) ;
00688    if( input_dset == NULL )
00689       return "*******************************\n"
00690              "CORREL_main:  bad input dataset\n"
00691              "*******************************"  ;
00692 
00693    nvox = DSET_NVOX(input_dset) ;
00694    ntim = DSET_NUM_TIMES(input_dset) ;
00695 
00696    ignore = (int) PLUTO_get_number(plint) ;
00697    if( ignore >= ntim-5 || ignore < 0 )
00698       return "******************************\n"
00699              "CORREL_main:  bad ignore count\n"
00700              "******************************" ;
00701 
00702    DSET_load(input_dset) ;
00703    if( DSET_ARRAY(input_dset,0) == NULL )
00704       return "********************************\n"
00705              "CORREL_main:  can't load dataset\n"
00706              "********************************"  ;
00707 
00708    /*-- read 2nd line --*/
00709 
00710    PLUTO_next_option(plint) ;
00711    tsim = PLUTO_get_timeseries(plint) ;
00712    if( tsim == NULL || tsim->nx < ntim )
00713       return "*****************************\n"
00714              "CORREL_main: bad input vector\n"
00715              "*****************************"  ;
00716 
00717    flim = mri_to_float(tsim) ; tsar = MRI_FLOAT_PTR(flim) ;
00718 
00719    polort = (int) PLUTO_get_number(plint) ;
00720 
00721    /*-- read optional lines --*/
00722 
00723    while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
00724 
00725       /*-- Mask range of values --*/
00726 
00727       if( strcmp(tag,"Range") == 0 ){
00728          if( mask_dset == NULL ){
00729             mri_free(flim) ;
00730             return "******************************************\n"
00731                    "CORREL_main:  Can't use Range without Mask\n"
00732                    "******************************************"  ;
00733          }
00734 
00735          mask_bot = PLUTO_get_number(plint) ;
00736          mask_top = PLUTO_get_number(plint) ;
00737          continue ;
00738       }
00739 
00740       /*-- Mask itself --*/
00741 
00742       if( strcmp(tag,"Mask") == 0 ){
00743 
00744          idc       = PLUTO_get_idcode(plint) ;
00745          mask_dset = PLUTO_find_dset(idc) ;
00746 
00747          if( mask_dset == NULL ){
00748             mri_free(flim) ;
00749             return "******************************\n"
00750                    "CORREL_main:  bad mask dataset\n"
00751                    "******************************"  ;
00752          }
00753 
00754          if( DSET_NVOX(mask_dset) != nvox ){
00755            mri_free(flim) ;
00756            return "************************************************************\n"
00757                   "CORREL_main: mask input dataset doesn't match source dataset\n"
00758                   "************************************************************" ;
00759          }
00760 
00761          miv = (int) PLUTO_get_number(plint) ;
00762          if( miv >= DSET_NVALS(mask_dset) || miv < 0 ){
00763             mri_free(flim) ;
00764             return "****************************************************\n"
00765                    "CORREL_main: mask dataset sub-brick index is illegal\n"
00766                    "****************************************************"  ;
00767          }
00768 
00769          DSET_load(mask_dset) ;
00770          if( DSET_ARRAY(mask_dset,miv) == NULL ){
00771             mri_free(flim) ;
00772             return "*************************************\n"
00773                    "CORREL_main:  can't load mask dataset\n"
00774                    "*************************************"  ;
00775          }
00776          continue ;
00777       }
00778    }
00779 
00780    /*------------------------------------------------------*/
00781    /*---------- At this point, the inputs are OK ----------*/
00782 
00783    /*-- build a byte mask array --*/
00784 
00785    if( mask_dset == NULL ){
00786       mmm = (byte *) malloc( sizeof(byte) * nvox ) ;
00787       if( mmm == NULL )
00788          return " \n*** Can't malloc workspace! ***\n" ;
00789       memset( mmm , 1, nvox ) ; mcount = nvox ;
00790    } else {
00791 
00792       mmm = THD_makemask( mask_dset , miv , mask_bot , mask_top ) ;
00793       if( mmm == NULL )
00794          return " \n*** Can't make mask for some reason! ***\n" ;
00795       mcount = THD_countmask( nvox , mmm ) ;
00796 
00797       if( !EQUIV_DSETS(mask_dset,input_dset) ) DSET_unload(mask_dset) ;
00798       if( mcount < 3 ){
00799          free(mmm) ;
00800          return " \n*** Less than 3 voxels survive the mask! ***\n" ;
00801       }
00802       sprintf(buf," \n"
00803                   " %d voxels in the mask\n"
00804                   " out of %d dataset voxels\n ",mcount,nvox) ;
00805       PLUTO_popup_transient(plint,buf) ;
00806    }
00807 
00808    /*-- setup to do the FIM calculation --*/
00809 
00810    fim_nref = polort+2 ;
00811    ref_vec = (float *) malloc( sizeof(float) * fim_nref ) ;
00812    vval    = (float *) malloc( sizeof(float) * mcount ) ;
00813 
00814    pc_ref = new_PCOR_references( fim_nref ) ;
00815    pc_vc  = new_PCOR_voxel_corr( mcount , fim_nref ) ;
00816 
00817    /*-- do FIM --*/
00818 
00819    for( nupdt=0,it=ignore ; it < ntim ; it++ ){
00820       if( tsar[it] >= WAY_BIG ) continue ;           /* skip this */
00821 
00822       ref_vec[0] = 1.0 ;
00823       for( ip=1 ; ip <= polort ; ip++ )
00824          ref_vec[ip] = ref_vec[ip-1] * ((float)it) ;
00825 
00826       ref_vec[ip] = tsar[it] ; /* vector value */
00827 
00828       update_PCOR_references( ref_vec , pc_ref ) ;
00829 
00830       /*-- load values into vval --*/
00831 
00832       switch( DSET_BRICK_TYPE(input_dset,it) ){
00833 
00834          case MRI_short:{
00835             short * bar = (short *) DSET_ARRAY(input_dset,it) ;
00836             float mfac = DSET_BRICK_FACTOR(input_dset,it) ;
00837             if( mfac == 0.0 ) mfac = 1.0 ;
00838             for( ii=jj=0 ; ii < nvox ; ii++ )
00839                if( mmm[ii] ) vval[jj++] = mfac*bar[ii] ;
00840          }
00841          break ;
00842 
00843          case MRI_byte:{
00844             byte * bar = (byte *) DSET_ARRAY(input_dset,it) ;
00845             float mfac = DSET_BRICK_FACTOR(input_dset,it) ;
00846             if( mfac == 0.0 ) mfac = 1.0 ;
00847             for( ii=jj=0 ; ii < nvox ; ii++ )
00848                if( mmm[ii] ) vval[jj++] = mfac*bar[ii] ;
00849          }
00850          break ;
00851 
00852          case MRI_float:{
00853             float * bar = (float *) DSET_ARRAY(input_dset,it) ;
00854             float mfac = DSET_BRICK_FACTOR(input_dset,it) ;
00855             if( mfac == 0.0 ) mfac = 1.0 ;
00856                for( ii=jj=0 ; ii < nvox ; ii++ )
00857                   if( mmm[ii] ) vval[jj++] = mfac*bar[ii] ;
00858          }
00859          break ;
00860       }
00861 
00862       PCOR_update_float( vval , pc_ref , pc_vc ) ;
00863       nupdt++ ;
00864    }
00865 
00866    free(ref_vec) ; mri_free(flim) ; free(mmm) ;
00867 
00868    /*-- get correlation coefficient --*/
00869 
00870    PCOR_get_pcor( pc_ref , pc_vc , vval ) ;
00871 
00872    free_PCOR_references(pc_ref) ; free_PCOR_voxel_corr(pc_vc) ;
00873 
00874    /*-- compute statistics --*/
00875 
00876    sum = 0.0 ;
00877    for( ii=0 ; ii < mcount ; ii++ ) sum += vval[ii] ;
00878    sum /= mcount ; sumq = 0.0 ;
00879    for( ii=0 ; ii < mcount ; ii++ )
00880       sumq += (vval[ii]-sum)*(vval[ii]-sum) ;
00881    sumq = sqrt(sumq/mcount) ;
00882 
00883    /*-- get robust statistics of Fisher z-transform --*/
00884 
00885    zval = (float *) malloc( sizeof(float) * mcount ) ;
00886    aval = (float *) malloc( sizeof(float) * mcount ) ;
00887    for( ii=0 ; ii < mcount ; ii++ ) zval[ii] = atanh(vval[ii]) ;
00888    qsort_float( mcount , zval ) ;
00889    if( mcount%2 == 1 )              /* median */
00890       zmid = zval[mcount/2] ;
00891    else
00892       zmid = 0.5 * ( zval[mcount/2] + zval[mcount/2-1] ) ;
00893 
00894    for( ii=0 ; ii < mcount ; ii++ ) aval[ii] = fabs(zval[ii]-zmid) ;
00895    qsort_float( mcount , aval ) ;
00896    if( mcount%2 == 1 )              /* MAD = median absolute deviation */
00897       zmed = aval[mcount/2] ;
00898    else
00899       zmed = 0.5 * ( aval[mcount/2] + aval[mcount/2-1] ) ;
00900    zsig = 1.4826 * zmed ;           /* estimate st.dev. */
00901                                     /* 1/1.4826 = sqrt(2)*erfinv(0.5) */
00902    free(aval) ;
00903    free(zval) ;
00904 
00905    /*-- do histogram --*/
00906 
00907    hbot = -1.0 ; htop = 1.0 ; nbin = 100 ; dbin = (htop-hbot)/nbin ;
00908 
00909    hbin = (int *) calloc((nbin+1),sizeof(int)) ;
00910    jbin = (int *) calloc((nbin+1),sizeof(int)) ;  /* 04 Oct 1999 */
00911 
00912 #ifndef FIT_FISHER
00913    sqp = 1.0/(sqrt(2.0*PI)*sumq) ;                /* Gaussian fit */
00914    for( ii=0 ; ii < nbin ; ii++ ){                /* to rho data  */
00915       gval = hbot + (ii+0.5)*dbin - sum ;
00916       gval = sqp * exp( -0.5*gval*gval/(sumq*sumq) ) ;
00917       jbin[ii] = (int)( mcount * dbin * gval + 0.5 ) ;
00918    }
00919 #else
00920    sqp = 1.0/(sqrt(2.0*PI)*zsig) ;                /* Gaussian fit */
00921    for( ii=0 ; ii < nbin ; ii++ ){                /* to z(rho)    */
00922       rval = hbot + (ii+0.5)*dbin ;
00923       gval = atanh(rval) - zmid ;
00924       gval = sqp * exp( -0.5*gval*gval/(zsig*zsig) )/sqrt(1.0-rval*rval) ;
00925       jbin[ii] = (int)( mcount * dbin * gval + 0.5 ) ;
00926    }
00927    sum  = tanh(zmid) ;
00928    sumq = 0.5*( tanh(zmid+zmed) - tanh(zmid-zmed) ) ;
00929 #endif
00930 
00931    kbin = (int *) calloc((nbin+1),sizeof(int)) ;
00932    for( ii=0 ; ii < nbin ; ii++ ){
00933       gval = correl_t2p( fabs(hbot+ii*dbin) ,
00934                          (double)nupdt , (double)1 , (double)(polort+1) ) ;
00935       sqp = correl_t2p( fabs(hbot+(ii+1)*dbin) ,
00936                         (double)nupdt , (double)1 , (double)(polort+1) ) ;
00937       kbin[ii] = (int)( 0.5*mcount*fabs(gval-sqp) ) ;
00938    }
00939    jist[0] = jbin ; jist[1] = kbin ;
00940 
00941    flim = mri_new_vol_empty( mcount,1,1 , MRI_float ) ;
00942    mri_fix_data_pointer( vval , flim ) ;
00943    mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
00944 
00945    { char * ps = my_getenv("PTAIL") ;
00946      float pp=0.0 ;
00947      if( ps != NULL ) pp = strtod(ps,NULL) ;
00948      set_unusuality_tail(pp) ;
00949    }
00950 
00951    for( ii=0 ; ii < mcount ; ii++ ) vval[ii] = -vval[ii] ;
00952    msum = unusuality( mcount , vval ) ;
00953    for( ii=0 ; ii < mcount ; ii++ ) vval[ii] = -vval[ii] ;
00954    psum = unusuality( mcount , vval ) ;
00955 
00956    sprintf(buf,"%s^.%s:\\rho_{mid}=%.2f\\pm%.2f,u=%.1f,%.1f",
00957                DSET_FILECODE(input_dset),
00958                (tsim->name != NULL) ? THD_trailname(tsim->name,0) : " " ,
00959                sum,sumq , psum,msum ) ;
00960 #ifdef DO_GREEN
00961    PLUTO_histoplot( nbin,hbot,htop,hbin ,
00962                     "Correlation Coefficient",NULL,buf , 2,jist ) ;
00963 #else
00964    PLUTO_histoplot( nbin,hbot,htop,hbin ,
00965                     "Correlation Coefficient",NULL,buf , 1,jist ) ;
00966 #endif
00967 
00968    mri_clear_data_pointer(flim) ; mri_free(flim) ;
00969    free(vval) ; free(hbin) ; free(jbin) ; free(kbin) ;
00970    return NULL ;
00971 }

00133 {
00134    MCW_idcode * idc ;
00135    THD_3dim_dataset * input_dset , * mask_dset=NULL ;
00136    int iv , mcount , nvox , ii,jj , nbin=-1 , do_mval=0,mval ;
00137    float mask_bot=666.0 , mask_top=-666.0 , hbot,htop ;
00138    float val_bot=666.0  , val_top=-666.0 ;
00139    char * tag , * str , buf[THD_MAX_NAME+16] ;
00140    byte * mmm ;
00141    MRI_IMAGE * flim ;
00142    float     * flar ;
00143    int       * hbin ;
00144    int smooth=0 ;      /* 03 Dec 2004 */
00145 
00146    int miv=0 ;
00147 
00148    int maxcount=0 ; /* 01 Mar 2001 */
00149    float hrad=0.0 ; /* 20 Mar 2001 */
00150 
00151    char * histout=NULL ; /* 05 Feb 2002 - VR */
00152    FILE * HISTOUT=NULL ; /* 05 Feb 2002 - VR */
00153    int writehist=0     ; /* 05 Feb 2002 - VR */
00154    float dx            ; /* 05 Feb 2002 - VR */
00155 
00156    /*--------------------------------------------------------------------*/
00157    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
00158 
00159    if( plint == NULL )
00160       return "***********************\n"
00161              "HISTO_main:  NULL input\n"
00162              "***********************"  ;
00163 
00164    /*-- read 1st line --*/
00165 
00166    PLUTO_next_option(plint) ;
00167    idc        = PLUTO_get_idcode(plint) ;
00168    input_dset = PLUTO_find_dset(idc) ;
00169    if( input_dset == NULL )
00170       return "******************************\n"
00171              "HISTO_main:  bad input dataset\n"
00172              "******************************"  ;
00173 
00174    iv = (int) PLUTO_get_number(plint) ;
00175    if( iv >= DSET_NVALS(input_dset) || iv < 0 )
00176       return "********************************\n"
00177              "HISTO_main:  bad input sub-brick\n"
00178              "********************************" ;
00179 
00180    DSET_load(input_dset) ;
00181    if( DSET_ARRAY(input_dset,iv) == NULL )
00182       return "*******************************\n"
00183              "HISTO_main:  can't load dataset\n"
00184              "*******************************"  ;
00185    nvox = DSET_NVOX(input_dset) ;
00186 
00187    /*-- read optional lines --*/
00188 
00189    while( (tag=PLUTO_get_optiontag(plint)) != NULL ){
00190 
00191       /*-- Dataset range of values --*/
00192 
00193       if( strcmp(tag,"Values") == 0 ){
00194          val_bot = PLUTO_get_number(plint) ;
00195          val_top = PLUTO_get_number(plint) ;
00196          do_mval = (val_bot < val_top) ;
00197          continue ;
00198       }
00199 
00200       /*-- Number of bins --*/
00201 
00202       if( strcmp(tag,"Bins") == 0 ){
00203          nbin     = PLUTO_get_number(plint) ;
00204          maxcount = PLUTO_get_number(plint) ;
00205          smooth   = PLUTO_get_number(plint) ;  /* 03 Dec 2004 */
00206          continue ;
00207       }
00208 
00209       /*-- Mask range of values --*/
00210 
00211       if( strcmp(tag,"Range") == 0 ){
00212          if( mask_dset == NULL )
00213             return "*****************************************\n"
00214                    "HISTO_main:  Can't use Range without Mask\n"
00215                    "*****************************************"  ;
00216 
00217          mask_bot = PLUTO_get_number(plint) ;
00218          mask_top = PLUTO_get_number(plint) ;
00219          continue ;
00220       }
00221 
00222       /*-- Mask itself --*/
00223 
00224       if( strcmp(tag,"Mask") == 0 ){
00225 
00226          idc       = PLUTO_get_idcode(plint) ;
00227          mask_dset = PLUTO_find_dset(idc) ;
00228 
00229          if( mask_dset == NULL )
00230             return "*****************************\n"
00231                    "HISTO_main:  bad mask dataset\n"
00232                    "*****************************"  ;
00233 
00234          if( DSET_NVOX(mask_dset) != nvox )
00235            return "***********************************************************\n"
00236                   "HISTO_main: mask input dataset doesn't match source dataset\n"
00237                   "***********************************************************" ;
00238 
00239          miv = (int) PLUTO_get_number(plint) ;  /* 06 Aug 1998 */
00240          if( miv >= DSET_NVALS(mask_dset) || miv < 0 )
00241             return "***************************************************\n"
00242                    "HISTO_main: mask dataset sub-brick index is illegal\n"
00243                    "***************************************************"  ;
00244 
00245          DSET_load(mask_dset) ;
00246          if( DSET_ARRAY(mask_dset,miv) == NULL )
00247             return "************************************\n"
00248                    "HISTO_main:  can't load mask dataset\n"
00249                    "************************************"  ;
00250          continue ;
00251       }
00252 
00253       /*-- 20 Mar 2001: Aboot --*/
00254 
00255       if( strcmp(tag,"Aboot") == 0 ){
00256          hrad = PLUTO_get_number(plint) ;
00257          continue ;
00258       }
00259 
00260       /*-- 05 Feb 2002: Output - VR --*/
00261 
00262       if( strcmp(tag,"Output") == 0 ){
00263          histout = PLUTO_get_string(plint) ;
00264          writehist = 1 ;
00265          continue ;
00266       }
00267    }
00268 
00269    /*------------------------------------------------------*/
00270    /*---------- At this point, the inputs are OK ----------*/
00271 
00272    /*-- build a byte mask array --*/
00273 
00274    if( mask_dset == NULL ){
00275       mmm = (byte *) malloc( sizeof(byte) * nvox ) ;
00276       if( mmm == NULL )
00277          return " \n*** Can't malloc workspace! ***\n" ;
00278       memset( mmm , 1, nvox ) ; mcount = nvox ;
00279    } else {
00280 
00281       mmm = THD_makemask( mask_dset , miv , mask_bot , mask_top ) ;
00282       if( mmm == NULL )
00283          return " \n*** Can't make mask for some reason! ***\n" ;
00284       mcount = THD_countmask( nvox , mmm ) ;
00285 
00286       if( !EQUIV_DSETS(mask_dset,input_dset) ) DSET_unload(mask_dset) ;
00287       if( mcount < 3 ){
00288          free(mmm) ;
00289          return " \n*** Less than 3 voxels survive the mask! ***\n" ;
00290       }
00291       sprintf(buf," \n"
00292                   " %d voxels in the mask\n"
00293                   " out of %d dataset voxels\n ",mcount,nvox) ;
00294       PLUTO_popup_transient(plint,buf) ;
00295    }
00296 
00297    /*-- 20 Mar 2001: modify mask via Aboot Radius, if present --*/
00298 
00299    if( hrad > 0.0 ){
00300       MCW_cluster * cl ;
00301       short *di,*dj,*dk ;
00302       int nd , xx,yy,zz , dd , nx,ny,nz,nxy, nx1,ny1,nz1 , ip,jp,kp ;
00303 
00304       cl = MCW_build_mask( 0,0,0 , fabs(DSET_DX(input_dset)) ,
00305                                    fabs(DSET_DY(input_dset)) ,
00306                                    fabs(DSET_DZ(input_dset)) , hrad ) ;
00307 
00308       if( cl == NULL || cl->num_pt < 6 ){
00309          KILL_CLUSTER(cl);
00310          PLUTO_popup_transient(plint, " \n"
00311                                       " Aboot Radius too small\n"
00312                                       " for this dataset!\n"     ) ;
00313       } else {
00314          ADDTO_CLUSTER(cl,0,0,0,0) ;
00315          di = cl->i ; dj = cl->j ; dk = cl->k ; nd = cl->num_pt ;
00316          nx = DSET_NX(input_dset) ; nx1 = nx-1 ;
00317          ny = DSET_NY(input_dset) ; ny1 = ny-1 ; nxy  = nx*ny ;
00318          nz = DSET_NZ(input_dset) ; nz1 = nz-1 ;
00319          xx = plint->im3d->vinfo->i1 ;
00320          yy = plint->im3d->vinfo->j2 ;
00321          zz = plint->im3d->vinfo->k3 ;
00322          for( dd=0 ; dd < nd ; dd++ ){
00323             ip = xx+di[dd] ; if( ip < 0 || ip > nx1 ) continue ;
00324             jp = yy+dj[dd] ; if( jp < 0 || jp > ny1 ) continue ;
00325             kp = zz+dk[dd] ; if( kp < 0 || kp > nz1 ) continue ;
00326             mmm[ip+jp*nx+kp*nxy]++ ;
00327          }
00328          KILL_CLUSTER(cl) ;
00329          for( dd=0 ; dd < nvox ; dd++ ) if( mmm[dd] == 1 ) mmm[dd] = 0 ;
00330          mcount = THD_countmask( nvox , mmm ) ;
00331 
00332          if( mcount < 3 ){
00333             free(mmm) ;
00334             return " \n*** Less than 3 voxels survive the mask+radius! ***\n" ;
00335          }
00336          sprintf(buf," \n"
00337                      " %d voxels in the mask+radius\n"
00338                      " out of %d dataset voxels\n ",mcount,nvox) ;
00339          PLUTO_popup_transient(plint,buf) ;
00340       }
00341    }
00342 
00343    /*-- check for text output of histogram - 05 Feb 2002 - VR --*/
00344 
00345    if ( writehist )
00346    {
00347       static char hbuf[1024] ;
00348       if ( ( histout == NULL ) || ( strlen (histout) == 0 ) ){
00349          sprintf( hbuf , "%s.histog" , DSET_PREFIX(input_dset) ) ;
00350       } else {
00351          strcpy( hbuf , histout ) ;
00352          if( strstr(hbuf,".hist") == NULL ) strcat( hbuf , ".histog" ) ;
00353       }
00354       histout = hbuf ;
00355 
00356       if (THD_is_file(histout))
00357       {
00358          free(mmm) ;
00359 
00360          return "*******************************\n"
00361                 "Outfile exists, won't overwrite\n"
00362                 "*******************************\n" ;
00363       }
00364       else {
00365          HISTOUT = fopen (histout, "w") ;
00366          if( HISTOUT == NULL ){
00367             free(mmm) ;
00368             return "**********************************\n"
00369                    "Can't open Outfile for some reason\n"
00370                    "**********************************\n" ;
00371          }
00372       }
00373    }
00374 
00375    /*-- allocate an array to histogrammatize --*/
00376 
00377    flim = mri_new( mcount , 1 , MRI_float ) ;
00378    flar = MRI_FLOAT_PTR(flim) ;
00379 
00380    /*-- load values into this array --*/
00381 
00382    switch( DSET_BRICK_TYPE(input_dset,iv) ){
00383       default:
00384          free(mmm) ; mri_free(flim) ;
00385          return "*** Can't use source dataset -- illegal data type! ***" ;
00386 
00387       case MRI_short:{
00388          short * bar = (short *) DSET_ARRAY(input_dset,iv) ;
00389          float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
00390          if( mfac == 0.0 ) mfac = 1.0 ;
00391          if( do_mval ){
00392             float val ;
00393             for( ii=jj=0 ; ii < nvox ; ii++ ){
00394                if( mmm[ii] ){
00395                   val = mfac*bar[ii] ;
00396                   if( val >= val_bot && val <= val_top ) flar[jj++] = val ;
00397                }
00398             }
00399             mval = jj ;
00400          } else {
00401             for( ii=jj=0 ; ii < nvox ; ii++ )
00402                if( mmm[ii] ) flar[jj++] = mfac*bar[ii] ;
00403          }
00404       }
00405       break ;
00406 
00407       case MRI_byte:{
00408          byte * bar = (byte *) DSET_ARRAY(input_dset,iv) ;
00409          float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
00410          if( mfac == 0.0 ) mfac = 1.0 ;
00411          if( do_mval ){
00412             float val ;
00413             for( ii=jj=0 ; ii < nvox ; ii++ ){
00414                if( mmm[ii] ){
00415                   val = mfac*bar[ii] ;
00416                   if( val >= val_bot && val <= val_top ) flar[jj++] = val ;
00417                }
00418             }
00419             mval = jj ;
00420          } else {
00421             for( ii=jj=0 ; ii < nvox ; ii++ )
00422                if( mmm[ii] ) flar[jj++] = mfac*bar[ii] ;
00423          }
00424       }
00425       break ;
00426 
00427       case MRI_float:{
00428          float * bar = (float *) DSET_ARRAY(input_dset,iv) ;
00429          float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
00430          if( mfac == 0.0 ) mfac = 1.0 ;
00431          if( do_mval ){
00432             float val ;
00433             for( ii=jj=0 ; ii < nvox ; ii++ ){
00434                if( mmm[ii] ){
00435                   val = mfac*bar[ii] ;
00436                   if( val >= val_bot && val <= val_top ) flar[jj++] = val ;
00437                }
00438             }
00439             mval = jj ;
00440          } else {
00441             for( ii=jj=0 ; ii < nvox ; ii++ )
00442                if( mmm[ii] ) flar[jj++] = mfac*bar[ii] ;
00443          }
00444       }
00445       break ;
00446    }
00447 
00448    if( do_mval ){
00449       if( mval == 0 ){
00450          free(mmm) ; mri_free(flim) ;
00451          return "*** Can't use source dataset -- no data in Values range ***" ;
00452       }
00453       flim->nx = flim->nvox = mcount = mval ;
00454    }
00455 
00456    /*-- set range and size of histogram --*/
00457 
00458    if( val_bot > val_top ){
00459       hbot = mri_min(flim) ; htop = mri_max(flim) ;
00460       if( hbot >= htop ){
00461          free(mmm) ; mri_free(flim) ;
00462          return "***********************************\n"
00463                 "Selected voxels have no data range!\n"
00464                 "***********************************"  ;
00465       }
00466    } else {
00467       hbot = val_bot ; htop = val_top ;
00468    }
00469 
00470    if( nbin < 10 || nbin > 1000 ){
00471       switch( DSET_BRICK_TYPE(input_dset,iv) ){
00472          case MRI_float:
00473             nbin = (int) sqrt((double)mcount) ;
00474          break ;
00475 
00476          case MRI_short:
00477          case MRI_byte:{
00478             float mfac = DSET_BRICK_FACTOR(input_dset,iv) ;
00479             if( mfac == 0.0 || mfac == 1.0 )
00480                nbin = (int)( htop - hbot ) ;
00481             else
00482                nbin = (int) sqrt((double)mcount) ;
00483          }
00484          break ;
00485 
00486       }
00487       if( nbin < 10 ) nbin = 10 ; else if( nbin > 1000 ) nbin = 1000 ;
00488    }
00489 
00490    /*-- actually compute and plot histogram --*/
00491 
00492    hbin = (int *) calloc((nbin+1),sizeof(int)) ;
00493 
00494    mri_histogram( flim , hbot,htop , TRUE , nbin,hbin ) ;
00495 
00496    if( smooth > 0 ){  /* 03 Dec 2004 */
00497      int nwid=smooth , *gbin=(int *)calloc((nbin+1),sizeof(int)) , ibot,itop ;
00498      float ws,wss , *wt ;
00499 
00500      ws = 0.0 ;
00501      wt = (float *)malloc(sizeof(float)*(2*nwid+1)) ;
00502      for( ii=0 ; ii <= 2*nwid ; ii++ ){
00503        wt[ii] = nwid-abs(nwid-ii) + 0.5f ;
00504        ws += wt[ii] ;
00505      }
00506      for( ii=0 ; ii <= 2*nwid ; ii++ ) wt[ii] /= ws ;
00507 
00508      for( jj=0 ; jj <= nbin ; jj++ ){
00509        ibot = jj-nwid ; if( ibot < 0    ) ibot = 0 ;
00510        itop = jj+nwid ; if( itop > nbin ) itop = nbin ;
00511        ws = wss = 0.0 ;
00512        for( ii=ibot ; ii <= itop ; ii++ ){
00513          ws += wt[nwid-jj+ii] * hbin[ii] ; wss += wt[nwid-jj+ii] ;
00514        }
00515        gbin[jj] = rint(ws/wss) ;
00516      }
00517      memcpy(hbin,gbin,sizeof(int)*(nbin+1)) ;
00518      free((void *)wt) ; free((void *)gbin) ;
00519    }
00520 
00521    if( maxcount > 0 ){
00522       for( ii=0 ; ii <= nbin ; ii++ ) hbin[ii] = MIN( hbin[ii] , maxcount ) ;
00523    }
00524    sprintf(buf,"\\noesc %s[%d] %d voxels",DSET_FILECODE(input_dset),iv,mcount);
00525    PLUTO_histoplot( nbin,hbot,htop,hbin , NULL , NULL ,  buf , 0,NULL ) ;
00526 
00527    /*-- 05 Feb 2002: Output - VR --*/
00528 
00529    if ( HISTOUT != NULL )
00530    {
00531       if( hbot >= htop ){ hbot = 0.0 ; htop = nbin ;}
00532 
00533       dx = (htop-hbot)/nbin ;
00534 
00535       for( ii=0 ; ii <= nbin ; ii++ )
00536          fprintf (HISTOUT, "%12.6f %13d \n", hbot+ii*dx, hbin[ii]) ;
00537 
00538       fclose (HISTOUT) ;
00539 
00540       fprintf (stderr, "%s written to disk \n", histout) ;
00541    }
00542 
00543    /*-- go home to mama --*/
00544 
00545    free(hbin) ; free(mmm) ; mri_free(flim) ; return NULL ;
00546 }

00056 {
00057    PLUGIN_interface * plint ;
00058 
00059 #if 0
00060    if( ncall > 0 ) return NULL ;  /* only one interface */
00061 #else
00062    if( ncall == 1 ) return CORREL_init() ;
00063    if( ncall >  1 ) return NULL ;
00064 #endif
00065 
00066    /*-- set titles and call point --*/
00067 
00068    plint = PLUTO_new_interface( "Histogram" ,
00069                                 "Histogram of Dataset Brick" ,
00070                                 helpstring ,
00071                                 PLUGIN_CALL_VIA_MENU , HISTO_main  ) ;
00072 
00073    PLUTO_add_hint( plint , "Histogram of Dataset Brick" ) ;
00074 
00075    PLUTO_set_sequence( plint , "A:afniinfo:dsethistog" ) ;
00076 
00077    PLUTO_set_runlabels( plint , "Plot+Keep" , "Plot+Close" ) ;  /* 04 Nov 2003 */
00078 
00079    /*-- first line of input --*/
00080 
00081    PLUTO_add_option( plint , "Source" , "Source" , TRUE ) ;
00082    PLUTO_add_dataset(plint , "Dataset" ,
00083                                     ANAT_ALL_MASK , FUNC_ALL_MASK ,
00084                                     DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00085 
00086    PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
00087 
00088    /*-- second line of input --*/
00089 
00090    PLUTO_add_option( plint , "Values" , "Values" , FALSE ) ;
00091    PLUTO_add_number( plint , "Bottom" , -99999,99999, 0,  1,1 ) ;
00092    PLUTO_add_number( plint , "Top"    , -99999,99999, 0, -1,1 ) ;
00093 
00094    /*-- third line of input --*/
00095 
00096    PLUTO_add_option( plint , "Bins" , "Bins" , FALSE ) ;
00097    PLUTO_add_number( plint , "Number" , 10,1000,0, 100,1 ) ;
00098    PLUTO_add_number( plint , "Max Count" , 0,999999999,0 , 0,1 ) ;
00099    PLUTO_add_number( plint , "Smooth"    , 0,99,0 , 0,1 ) ;  /* 03 Dec 2004 */
00100 
00101    /*-- fourth line of input --*/
00102 
00103    PLUTO_add_option( plint , "Mask" , "Mask" , FALSE ) ;
00104    PLUTO_add_dataset( plint , "Dataset" ,
00105                                     ANAT_ALL_MASK , FUNC_ALL_MASK ,
00106                                     DIMEN_ALL_MASK | BRICK_ALLREAL_MASK ) ;
00107    PLUTO_add_number( plint , "Sub-brick" , 0,9999,0 , 0,1 ) ;
00108 
00109    /*-- fifth line of input --*/
00110 
00111    PLUTO_add_option( plint , "Range"  , "Range" , FALSE ) ;
00112    PLUTO_add_number( plint , "Bottom" , -99999,99999, 0,  1,1 ) ;
00113    PLUTO_add_number( plint , "Top"    , -99999,99999, 0, -1,1 ) ;
00114 
00115    /*-- sixth line of input [20 Mar 2001] --*/
00116 
00117    PLUTO_add_option( plint , "Aboot" , "Aboot" , FALSE ) ;
00118    PLUTO_add_number( plint , "Radius" , 2,100,0,10,1 ) ;
00119 
00120    /*-- seventh line of input [05 Feb 2002 - VR] --*/
00121 
00122    PLUTO_add_option( plint , "Output", "Output", FALSE ) ;
00123    PLUTO_add_string( plint , "Filename", 0 , NULL , 20 ) ;
00124 
00125    return plint ;
00126 }

   " Purpose: Plot a histogram of data from a dataset brick.\n"
   "\n"
   " Source:  Dataset   = data to be processed\n"
   "          Sub-brick = which one to use\n\n"
   " Values:  Bottom    = minimum value from dataset to include\n"
   "          Top       = maximum value from dataset to include\n\n"
   " Bins:    Number    = number of bins to use\n"
   "          Max Count = maximum count per bin\n\n"
   "          Smooth    = +/- bin range to smooth histogram\n"
   " Mask:    Dataset   = masking dataset\n"
   "          Sub-brick = which one to use\n\n"
   " Range:   Bottom    = min value from mask dataset to use\n"
   "          Top       = max value from mask dataset to use\n"
   "          [if Bottom >  Top, then all nonzero mask voxels will be used; ]\n"
   "          [if Bottom <= Top, then only nonzero mask voxels in this range]\n"
   "          [                  will be used in computing the statistics.  ]\n"
   " Aboot:   If activated, then only voxels within a distance of Radius mm\n"
   "          of the current crosshairs will be used in the histogram.\n"
   " Output:  Name of the ascii file to which histogram values are written.\n"
   "\n"
   " Author -- RW Cox - 30 September 1999\n"
   " Output feature added by V Roopchansingh, Feb 2002\n"