00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 #include "thd.h"
00009 
00010 /*----------------------------------------------------------------------*/
00011 /*! Load the statistics of a dataset (modified Nov 15 1995)
00012 ------------------------------------------------------------------------*/
00013 
00014 void THD_load_statistics( THD_3dim_dataset *dset )
00015 {
00016    int ii , mmin , mmax , ibr , nbr ;
00017    short *brkk ;
00018    THD_brick_stats *bsold ;
00019 
00020    /*-- sanity checks --*/
00021 
00022    if( ! ISVALID_3DIM_DATASET(dset) ) return ;
00023 
00024    nbr = THD_count_databricks( dset->dblk ) ;
00025 
00026    /*-- 3/24/95: if don't get data, try for the warp parent --*/
00027 
00028    if( nbr == 0 ){
00029 
00030      if( ! ISVALID_3DIM_DATASET(dset->warp_parent) ) return ;     /* nothing */
00031      if( dset->warp_parent == dset )                 return ;
00032 
00033      RELOAD_STATS( dset->warp_parent ) ;                       /* recursion! */
00034      if( ! ISVALID_STATISTIC(dset->warp_parent->stats) ) return ; /* nothing */
00035 
00036      if( dset->stats == NULL ){                     /* create if not present */
00037        dset->stats = myXtNew( THD_statistics ) ;
00038        ADDTO_KILL( dset->kl , dset->stats ) ;
00039        dset->stats->type   = STATISTICS_TYPE ;
00040        dset->stats->parent = (XtPointer) dset ;
00041        dset->stats->bstat  = NULL ;
00042      }
00043 
00044      bsold               = dset->stats->bstat ;
00045      dset->stats->nbstat = dset->dblk->nvals ;
00046      dset->stats->bstat  = (THD_brick_stats *)
00047                             XtRealloc( (char *) bsold ,
00048                                 sizeof(THD_brick_stats) * dset->dblk->nvals ) ;
00049      if( bsold != dset->stats->bstat )
00050         REPLACE_KILL( dset->kl , bsold , dset->stats->bstat ) ;
00051 
00052      /** copy stats from warp parent for each brick **/
00053 
00054      for( ibr=0 ; ibr < dset->dblk->nvals ; ibr++ ){
00055        if( ibr < dset->warp_parent->stats->nbstat )
00056          dset->stats->bstat[ibr] = dset->warp_parent->stats->bstat[ibr] ;
00057        else
00058          INVALIDATE_BSTAT( dset->stats->bstat[ibr] ) ;
00059      }
00060 
00061      return ;
00062    }
00063 
00064    /*-- if here, have good data in this dataset --*/
00065 
00066    if( dset->stats == NULL ){                  /* create if not present */
00067      dset->stats = myXtNew( THD_statistics ) ;
00068      ADDTO_KILL( dset->kl , dset->stats ) ;
00069      dset->stats->type   = STATISTICS_TYPE ;
00070      dset->stats->parent = (XtPointer) dset ;
00071      dset->stats->bstat  = NULL ;
00072    }
00073 
00074    bsold               = dset->stats->bstat ;
00075    dset->stats->nbstat = dset->dblk->nvals ;
00076    dset->stats->bstat  = (THD_brick_stats *)
00077                           XtRealloc( (char *) bsold ,
00078                               sizeof(THD_brick_stats) * dset->dblk->nvals ) ;
00079    if( bsold != dset->stats->bstat )
00080      REPLACE_KILL( dset->kl , bsold , dset->stats->bstat ) ;
00081 
00082    /* 3/24/95: load stats for each sub-brick, not just the first */
00083 
00084    for( ibr=0 ; ibr < dset->dblk->nvals ; ibr++ ){  /* for each sub-brick */
00085      dset->stats->bstat[ibr] = THD_get_brick_stats( DSET_BRICK(dset,ibr) ) ;
00086 
00087      /* 11/21/95: allow for scaling factor that may be applied to data */
00088 
00089      if( DSET_BRICK_FACTOR(dset,ibr) > 0.0 ){
00090        dset->stats->bstat[ibr].min *= DSET_BRICK_FACTOR(dset,ibr) ;
00091        dset->stats->bstat[ibr].max *= DSET_BRICK_FACTOR(dset,ibr) ;
00092      }
00093    }
00094    return ;
00095 }
00096 
00097 /*--------------------------------------------------------------*/
00098 /*! Compute statistics for a 3D brick of varying data type
00099 ----------------------------------------------------------------*/
00100 
00101 THD_brick_stats THD_get_brick_stats( MRI_IMAGE *im )
00102 {
00103    register int ii , nvox ;
00104    register float bot , top ;
00105    void *br ;
00106    THD_brick_stats bst ;
00107 
00108    bst.min = bst.max = 0 ;  /* got to give them something */
00109 
00110    if( im == NULL ) return bst ;
00111    br   = mri_data_pointer( im ) ; if( br == NULL ) return bst ;
00112    nvox = im->nvox ;
00113 
00114    switch( im->kind ){
00115 
00116      default:
00117         INVALIDATE_BSTAT(bst) ;
00118      break ;
00119 
00120      case MRI_rgb:{
00121        register byte *ar = (byte *) br ; register float val ;
00122        bot = top = 0 ;
00123        for( ii=0 ; ii < nvox ; ii++ ){      /* scale to brightness */
00124          val =  0.299 * ar[3*ii]           /* between 0 and 255   */
00125               + 0.587 * ar[3*ii+1]
00126               + 0.114 * ar[3*ii+2] ;
00127               if( bot > val ) bot = val ;
00128          else if( top < val ) top = val ;
00129        }
00130      }
00131      break ;
00132 
00133      case MRI_byte:{
00134        register byte *ar = (byte *) br ;
00135        bot = top = ar[0] ;
00136        for( ii=1 ; ii < nvox ; ii++ ){
00137               if( bot > ar[ii] ) bot = ar[ii] ;
00138          else if( top < ar[ii] ) top = ar[ii] ;
00139        }
00140      }
00141      break ;
00142 
00143      case MRI_short:{
00144        register short *ar = (short *) br ;
00145        bot = top = ar[0] ;
00146        for( ii=1 ; ii < nvox ; ii++ ){
00147               if( bot > ar[ii] ) bot = ar[ii] ;
00148          else if( top < ar[ii] ) top = ar[ii] ;
00149        }
00150      }
00151      break ;
00152 
00153 #if 0
00154      case MRI_int:{
00155        register int *ar = (int *) br ;
00156        bot = top = ar[0] ;
00157        for( ii=1 ; ii < nvox ; ii++ ){
00158               if( bot > ar[ii] ) bot = ar[ii] ;
00159          else if( top < ar[ii] ) top = ar[ii] ;
00160        }
00161      }
00162      break ;
00163 #endif
00164 
00165      case MRI_float:{
00166        register float *ar = (float *) br ;
00167        bot = top = ar[0] ;
00168        for( ii=1 ; ii < nvox ; ii++ ){
00169               if( bot > ar[ii] ) bot = ar[ii] ;
00170          else if( top < ar[ii] ) top = ar[ii] ;
00171        }
00172      }
00173      break ;
00174 
00175 #if 0
00176      case MRI_double:{
00177        register double *ar = (double *) br ;
00178        bot = top = ar[0] ;
00179        for( ii=1 ; ii < nvox ; ii++ ){
00180               if( bot > ar[ii] ) bot = ar[ii] ;
00181          else if( top < ar[ii] ) top = ar[ii] ;
00182        }
00183      }
00184      break ;
00185 #endif
00186 
00187      case MRI_complex:{
00188        register complex *ar = (complex *) br ;
00189        register float zz ;
00190        bot = top = CSQR(ar[0]) ;
00191        for( ii=1 ; ii < nvox ; ii++ ){
00192           zz = CSQR(ar[ii]) ;
00193                if( bot > zz ) bot = zz ;
00194           else if( top < zz ) top = zz ;
00195        }
00196        bot = sqrt(bot) ; top = sqrt(top) ;
00197      }
00198      break ;
00199 
00200    } /* end of switch on sub-brick type */
00201 
00202    bst.min = bot ; bst.max = top ;
00203    return bst ;
00204 }
00205 
00206 /*----------------------------------------------------------------------*/
00207 /*! Update the statistics of a dataset
00208    (use only if new bricks are added -- see EDIT_add_bricklist)
00209 ------------------------------------------------------------------------*/
00210 
00211 void THD_update_statistics( THD_3dim_dataset *dset )
00212 {
00213    Boolean good ;
00214    int ii , mmin , mmax , ibr , nbsold , nbr ;
00215    THD_brick_stats *bsold ;
00216    short *brkk ;
00217 
00218    /*-- sanity checks --*/
00219 
00220    if( ! ISVALID_3DIM_DATASET(dset) ) return ;
00221 
00222    nbr = THD_count_databricks( dset->dblk ) ;
00223 
00224    if( nbr == 0 ) return ;
00225 
00226    /*-- if here, have good data in this dataset --*/
00227 
00228    if( dset->stats == NULL ){                  /* create if not present */
00229      dset->stats = myXtNew( THD_statistics ) ;
00230      ADDTO_KILL( dset->kl , dset->stats ) ;
00231      dset->stats->type   = STATISTICS_TYPE ;
00232      dset->stats->parent = (XtPointer) dset ;
00233      dset->stats->bstat  = NULL ;
00234      dset->stats->nbstat = 0 ;
00235      nbsold              = 0 ;
00236    } else {
00237      nbsold              = dset->stats->nbstat ;
00238    }
00239 
00240    if( dset->dblk->nvals > nbsold ){
00241      bsold               = dset->stats->bstat ;
00242      dset->stats->nbstat = dset->dblk->nvals ;
00243      dset->stats->bstat  = (THD_brick_stats *)
00244                            XtRealloc( (char *) bsold ,
00245                                 sizeof(THD_brick_stats) * dset->dblk->nvals ) ;
00246      if( bsold != dset->stats->bstat )
00247        REPLACE_KILL( dset->kl , bsold , dset->stats->bstat ) ;
00248 
00249      for( ibr=nbsold ; ibr < dset->dblk->nvals ; ibr++ )  /* 11 Mar 2005 */
00250        INVALIDATE_BSTAT( dset->stats->bstat[ibr] ) ;
00251    }
00252 
00253    /* 28 Apr 1997: load stats for new sub-bricks, not all */
00254 
00255    for( ibr=0 ; ibr < dset->dblk->nvals ; ibr++ ){
00256 
00257      if( ibr >= nbsold || ! ISVALID_BSTAT(dset->stats->bstat[ibr]) ){
00258         dset->stats->bstat[ibr] = THD_get_brick_stats( DSET_BRICK(dset,ibr) ) ;
00259 
00260         if( DSET_BRICK_FACTOR(dset,ibr) > 0.0 ){
00261           dset->stats->bstat[ibr].min *= DSET_BRICK_FACTOR(dset,ibr) ;
00262           dset->stats->bstat[ibr].max *= DSET_BRICK_FACTOR(dset,ibr) ;
00263         }
00264      }
00265    }
00266    return ;
00267 }
00268 
00269 /*----------------------------------------------------------------------*/
00270 /*! Update the statistics of one sub-brick in a dataset. [29 Mar 2005]
00271 ------------------------------------------------------------------------*/
00272 
00273 void THD_update_one_bstat( THD_3dim_dataset *dset , int iv )
00274 {
00275    Boolean good ;
00276    int ii , mmin , mmax , ibr , nbsold , nbr ;
00277    THD_brick_stats *bsold ;
00278    short *brkk ;
00279 
00280    /*-- sanity checks --*/
00281 
00282    if( ! ISVALID_3DIM_DATASET(dset)     ) return ;
00283    if( iv < 0 || iv >= DSET_NVALS(dset) ) return ;
00284 
00285    /*-- if here, have good data in this dataset --*/
00286 
00287    if( dset->stats == NULL ){                  /* create if not present */
00288      dset->stats = myXtNew( THD_statistics ) ;
00289      ADDTO_KILL( dset->kl , dset->stats ) ;
00290      dset->stats->type   = STATISTICS_TYPE ;
00291      dset->stats->parent = (XtPointer) dset ;
00292      dset->stats->bstat  = NULL ;
00293      dset->stats->nbstat = 0 ;
00294      nbsold              = 0 ;
00295    } else {
00296      nbsold              = dset->stats->nbstat ;
00297    }
00298 
00299    if( dset->dblk->nvals > nbsold ){
00300      bsold               = dset->stats->bstat ;
00301      dset->stats->nbstat = dset->dblk->nvals ;
00302      dset->stats->bstat  = (THD_brick_stats *)
00303                            XtRealloc( (char *) bsold ,
00304                                 sizeof(THD_brick_stats) * dset->dblk->nvals ) ;
00305      if( bsold != dset->stats->bstat )
00306        REPLACE_KILL( dset->kl , bsold , dset->stats->bstat ) ;
00307 
00308      for( ibr=nbsold ; ibr < dset->dblk->nvals ; ibr++ )  /* 11 Mar 2005 */
00309        INVALIDATE_BSTAT( dset->stats->bstat[ibr] ) ;
00310    }
00311 
00312    if( iv >= nbsold || ! ISVALID_BSTAT(dset->stats->bstat[iv]) ){
00313      dset->stats->bstat[iv] = THD_get_brick_stats( DSET_BRICK(dset,iv) ) ;
00314 
00315      if( DSET_BRICK_FACTOR(dset,iv) > 0.0 ){
00316        dset->stats->bstat[iv].min *= DSET_BRICK_FACTOR(dset,iv) ;
00317        dset->stats->bstat[iv].max *= DSET_BRICK_FACTOR(dset,iv) ;
00318      }
00319    }
00320 
00321    return ;
00322 }