1 /*****************************************************************************
2    Major portions of this software are copyrighted by the Medical College
3    of Wisconsin, 1994-2000, and are released under the Gnu General Public
4    License, Version 2.  See the file README.Copyright for details.
5 ******************************************************************************/
6 
7 #include "mrilib.h"
8 #include "thd.h"
9 
10 /*-------------------------------------------------------------------------*/
11 /*! Load the statistics of a dataset [modified Nov 15 1995]
12 ---------------------------------------------------------------------------*/
13 
THD_load_statistics(THD_3dim_dataset * dset)14 void THD_load_statistics( THD_3dim_dataset *dset )
15 {
16    int ii , mmin , mmax , ibr , nbr ;
17    short *brkk ;
18    THD_brick_stats *bsold ;
19 
20    /*-- sanity checks --*/
21 
22    if( ! ISVALID_3DIM_DATASET(dset) ) return ;
23 
24    nbr = THD_count_databricks( dset->dblk ) ;
25 
26    /*-- 3/24/95: if don't get data, try for the warp parent --*/
27 
28    if( nbr == 0 ){
29 
30      if( ! ISVALID_3DIM_DATASET(dset->warp_parent) ) return ;     /* nothing */
31      if( dset->warp_parent == dset )                 return ;
32 
33      RELOAD_STATS( dset->warp_parent ) ;                       /* recursion! */
34      if( ! ISVALID_STATISTIC(dset->warp_parent->stats) ) return ; /* nothing */
35 
36      if( dset->stats == NULL ){                     /* create if not present */
37        dset->stats = myRwcNew( THD_statistics ) ;
38        ADDTO_KILL( dset->kl , dset->stats ) ;
39        dset->stats->type   = STATISTICS_TYPE ;
40        dset->stats->parent = (RwcPointer) dset ;
41        dset->stats->bstat  = NULL ;
42      }
43 
44      bsold               = dset->stats->bstat ;
45      dset->stats->nbstat = dset->dblk->nvals ;
46      dset->stats->bstat  = (THD_brick_stats *)
47                             RwcRealloc( (char *) bsold ,
48                                 sizeof(THD_brick_stats) * dset->dblk->nvals ) ;
49      if( bsold != dset->stats->bstat )
50         REPLACE_KILL( dset->kl , bsold , dset->stats->bstat ) ;
51 
52      /** copy stats from warp parent for each brick **/
53 
54      for( ibr=0 ; ibr < dset->dblk->nvals ; ibr++ ){
55        if( ibr < dset->warp_parent->stats->nbstat )
56          dset->stats->bstat[ibr] = dset->warp_parent->stats->bstat[ibr] ;
57        else
58          INVALIDATE_BSTAT( dset->stats->bstat[ibr] ) ;
59      }
60 
61      return ;
62    }
63 
64    /*-- if here, have good data in this dataset --*/
65 
66    if( dset->stats == NULL ){                  /* create if not present */
67      dset->stats = myRwcNew( THD_statistics ) ;
68      ADDTO_KILL( dset->kl , dset->stats ) ;
69      dset->stats->type   = STATISTICS_TYPE ;
70      dset->stats->parent = (RwcPointer) dset ;
71      dset->stats->bstat  = NULL ;
72    }
73 
74    bsold               = dset->stats->bstat ;
75    dset->stats->nbstat = dset->dblk->nvals ;
76    dset->stats->bstat  = (THD_brick_stats *)
77                           RwcRealloc( (char *) bsold ,
78                               sizeof(THD_brick_stats) * dset->dblk->nvals ) ;
79    if( bsold != dset->stats->bstat )
80      REPLACE_KILL( dset->kl , bsold , dset->stats->bstat ) ;
81 
82    /* 3/24/95: load stats for each sub-brick, not just the first */
83 
84    for( ibr=0 ; ibr < dset->dblk->nvals ; ibr++ ){  /* for each sub-brick */
85      dset->stats->bstat[ibr] = THD_get_brick_stats( DSET_BRICK(dset,ibr) ) ;
86 
87      /* 11/21/95: allow for scaling factor that may be applied to data */
88 
89      if( DSET_BRICK_FACTOR(dset,ibr) > 0.0 ){
90        dset->stats->bstat[ibr].min *= DSET_BRICK_FACTOR(dset,ibr) ;
91        dset->stats->bstat[ibr].max *= DSET_BRICK_FACTOR(dset,ibr) ;
92      }
93    }
94    return ;
95 }
96 
97 /*--------------------------------------------------------------*/
98 /*! Compute statistics for a 3D brick of varying data type
99 ----------------------------------------------------------------*/
100 
THD_get_brick_stats(MRI_IMAGE * im)101 THD_brick_stats THD_get_brick_stats( MRI_IMAGE *im )
102 {
103    int64_t ii , nvox ;  /* allow for big volumes    26 Mar 2016 [rickr] */
104    register float bot , top ;
105    void *br ;
106    THD_brick_stats bst ;
107 
108    bst.min = bst.max = 0 ;  /* got to give them something */
109 
110    if( im == NULL ) return bst ;
111    br   = mri_data_pointer( im ) ; if( br == NULL ) return bst ;
112    nvox = im->nvox ;
113 
114    switch( im->kind ){
115 
116      default:
117         INVALIDATE_BSTAT(bst) ; bot = top = 0 ;
118      break ;
119 
120      case MRI_rgb:{
121        register byte *ar = (byte *) br ; register float val ;
122        bot = top = 0 ;
123        for( ii=0 ; ii < nvox ; ii++ ){      /* scale to brightness */
124          val =  0.299 * ar[3*ii]           /* between 0 and 255   */
125               + 0.587 * ar[3*ii+1]
126               + 0.114 * ar[3*ii+2] ;
127               if( bot > val ) bot = val ;
128          else if( top < val ) top = val ;
129        }
130      }
131      break ;
132 
133      case MRI_byte:{
134        register byte *ar = (byte *) br ;
135        bot = top = ar[0] ;
136        for( ii=1 ; ii < nvox ; ii++ ){
137               if( bot > ar[ii] ) bot = ar[ii] ;
138          else if( top < ar[ii] ) top = ar[ii] ;
139        }
140      }
141      break ;
142 
143      case MRI_short:{
144        register short *ar = (short *) br ;
145        bot = top = ar[0] ;
146        for( ii=1 ; ii < nvox ; ii++ ){
147               if( bot > ar[ii] ) bot = ar[ii] ;
148          else if( top < ar[ii] ) top = ar[ii] ;
149        }
150      }
151      break ;
152 
153 #if 0
154      case MRI_int:{
155        register int *ar = (int *) br ;
156        bot = top = ar[0] ;
157        for( ii=1 ; ii < nvox ; ii++ ){
158               if( bot > ar[ii] ) bot = ar[ii] ;
159          else if( top < ar[ii] ) top = ar[ii] ;
160        }
161      }
162      break ;
163 #endif
164 
165      case MRI_float:{
166        register float *ar = (float *) br ;
167        bot = top = ar[0] ;
168        for( ii=1 ; ii < nvox ; ii++ ){
169               if( bot > ar[ii] ) bot = ar[ii] ;
170          else if( top < ar[ii] ) top = ar[ii] ;
171        }
172      }
173      break ;
174 
175 #if 0
176      case MRI_double:{
177        register double *ar = (double *) br ;
178        bot = top = ar[0] ;
179        for( ii=1 ; ii < nvox ; ii++ ){
180               if( bot > ar[ii] ) bot = ar[ii] ;
181          else if( top < ar[ii] ) top = ar[ii] ;
182        }
183      }
184      break ;
185 #endif
186 
187      case MRI_complex:{
188        register complex *ar = (complex *) br ;
189        register float zz ;
190        bot = top = CABS(ar[0]) ;
191        for( ii=1 ; ii < nvox ; ii++ ){
192           zz = CABS(ar[ii]) ;
193                if( bot > zz ) bot = zz ;
194           else if( top < zz ) top = zz ;
195        }
196      }
197      break ;
198 
199    } /* end of switch on sub-brick type */
200 
201    bst.min = bot ; bst.max = top ;
202    return bst ;
203 }
204 
205 /*----------------------------------------------------------------------*/
206 /*! Update the statistics of a dataset
207    (use only if new bricks are added -- see EDIT_add_bricklist)
208 ------------------------------------------------------------------------*/
209 
THD_update_statistics(THD_3dim_dataset * dset)210 void THD_update_statistics( THD_3dim_dataset *dset )
211 {
212    RwcBoolean good ;
213    int ii , mmin , mmax , ibr , nbsold , nbr ;
214    THD_brick_stats *bsold ;
215    short *brkk ;
216 
217    /*-- sanity checks --*/
218 
219    if( ! ISVALID_3DIM_DATASET(dset) ) return ;
220 
221    nbr = THD_count_databricks( dset->dblk ) ;
222 
223    if( nbr == 0 ) return ;
224 
225    /*-- if here, have good data in this dataset --*/
226 
227    if( dset->stats == NULL ){                  /* create if not present */
228      dset->stats = myRwcNew( THD_statistics ) ;
229      ADDTO_KILL( dset->kl , dset->stats ) ;
230      dset->stats->type   = STATISTICS_TYPE ;
231      dset->stats->parent = (RwcPointer) dset ;
232      dset->stats->bstat  = NULL ;
233      dset->stats->nbstat = 0 ;
234      nbsold              = 0 ;
235    } else {
236      nbsold              = dset->stats->nbstat ;
237    }
238 
239    if( dset->dblk->nvals > nbsold ){
240      bsold               = dset->stats->bstat ;
241      dset->stats->nbstat = dset->dblk->nvals ;
242      dset->stats->bstat  = (THD_brick_stats *)
243                            RwcRealloc( (char *) bsold ,
244                                 sizeof(THD_brick_stats) * dset->dblk->nvals ) ;
245      if( bsold != dset->stats->bstat )
246        REPLACE_KILL( dset->kl , bsold , dset->stats->bstat ) ;
247 
248      for( ibr=nbsold ; ibr < dset->dblk->nvals ; ibr++ )  /* 11 Mar 2005 */
249        INVALIDATE_BSTAT( dset->stats->bstat[ibr] ) ;
250    }
251 
252    /* 28 Apr 1997: load stats for new sub-bricks, not all */
253 
254    for( ibr=0 ; ibr < dset->dblk->nvals ; ibr++ ){
255 
256      if( ibr >= nbsold || ! ISVALID_BSTAT(dset->stats->bstat[ibr]) ){
257         dset->stats->bstat[ibr] = THD_get_brick_stats( DSET_BRICK(dset,ibr) ) ;
258 
259         if( DSET_BRICK_FACTOR(dset,ibr) > 0.0 ){
260           dset->stats->bstat[ibr].min *= DSET_BRICK_FACTOR(dset,ibr) ;
261           dset->stats->bstat[ibr].max *= DSET_BRICK_FACTOR(dset,ibr) ;
262         }
263      }
264    }
265    return ;
266 }
267 
268 /*----------------------------------------------------------------------*/
269 /*! Update the statistics of one sub-brick in a dataset. [29 Mar 2005]
270 ------------------------------------------------------------------------*/
271 
THD_update_one_bstat(THD_3dim_dataset * dset,int iv)272 void THD_update_one_bstat( THD_3dim_dataset *dset , int iv )
273 {
274    RwcBoolean good ;
275    int ii , mmin , mmax , ibr , nbsold , nbr ;
276    THD_brick_stats *bsold ;
277    short *brkk ;
278 
279    /*-- sanity checks --*/
280 
281    if( ! ISVALID_3DIM_DATASET(dset)     ) return ;
282    if( iv < 0 || iv >= DSET_NVALS(dset) ) return ;
283 
284    /*-- if here, have good data in this dataset --*/
285 
286    if( dset->stats == NULL ){                  /* create if not present */
287      dset->stats = myRwcNew( THD_statistics ) ;
288      ADDTO_KILL( dset->kl , dset->stats ) ;
289      dset->stats->type   = STATISTICS_TYPE ;
290      dset->stats->parent = (RwcPointer) dset ;
291      dset->stats->bstat  = NULL ;
292      dset->stats->nbstat = 0 ;
293      nbsold              = 0 ;
294    } else {
295      nbsold              = dset->stats->nbstat ;
296    }
297 
298    if( dset->dblk->nvals > nbsold ){
299      bsold               = dset->stats->bstat ;
300      dset->stats->nbstat = dset->dblk->nvals ;
301      dset->stats->bstat  = (THD_brick_stats *)
302                            RwcRealloc( (char *) bsold ,
303                                 sizeof(THD_brick_stats) * dset->dblk->nvals ) ;
304      if( bsold != dset->stats->bstat )
305        REPLACE_KILL( dset->kl , bsold , dset->stats->bstat ) ;
306 
307      for( ibr=nbsold ; ibr < dset->dblk->nvals ; ibr++ )  /* 11 Mar 2005 */
308        INVALIDATE_BSTAT( dset->stats->bstat[ibr] ) ;
309    }
310 
311    if( iv >= nbsold || ! ISVALID_BSTAT(dset->stats->bstat[iv]) ){
312      dset->stats->bstat[iv] = THD_get_brick_stats( DSET_BRICK(dset,iv) ) ;
313 
314      if( DSET_BRICK_FACTOR(dset,iv) > 0.0 ){
315        dset->stats->bstat[iv].min *= DSET_BRICK_FACTOR(dset,iv) ;
316        dset->stats->bstat[iv].max *= DSET_BRICK_FACTOR(dset,iv) ;
317      }
318    }
319 
320    return ;
321 }
322 
323 /*
324    Multiply values in dset by fac
325    Return number of sub-bricks for which scaling could not be done.
326 */
THD_dset_scale(THD_3dim_dataset * aset,float fac)327 int THD_dset_scale(THD_3dim_dataset *aset, float fac)
328 {
329    int ii, jj, err=0;
330    float fac0 = 1.0f, *fv=NULL;
331 
332 ENTRY("THD_dset_scale");
333 
334    if( !ISVALID_DSET(aset) ){
335      ERROR_message("THD_dset_scale: invalid dataset input :(") ;
336      RETURN(666) ;
337    }
338    if( fac == 0.0f ){
339      ERROR_message("THD_dset_scale: can't scale by 0 :(") ;
340      RETURN(DSET_NVALS(aset)) ;
341    }
342 
343    DSET_load(aset) ;
344    if( !DSET_LOADED(aset) ){
345      ERROR_message("THD_dset_scale: can't load dataset :(") ;
346      RETURN(DSET_NVALS(aset)) ;
347    }
348 
349    for (ii=0; ii<DSET_NVALS(aset); ++ii) {
350       switch (DSET_BRICK_TYPE(aset,ii)) {
351          case MRI_short:
352          case MRI_byte:
353             fac0 = DSET_BRICK_FACTOR(aset,ii);
354             if (fac0 == 0.0f) fac0 = 1.0f ;
355             EDIT_BRICK_FACTOR( aset,ii,fac0*fac ) ;
356             break;
357          case MRI_float:
358             fv = (float *)DSET_ARRAY(aset,ii);
359             if( fv == NULL ){
360               ERROR_message("THD_dset_scale: bad float array at ii=%d",ii) ;
361               err++ ;
362             } else {
363               for (jj=0; jj<DSET_NVOX(aset); ++jj) {
364                 fv[jj] *= fac;
365               }
366             }
367             break;
368          default:
369             if(!err)
370                ERROR_message("THD_dset_scale: not ready for data type %d [%s]\n"
371                              "   --> Sub-bricks of such types are untouched :(",
372                              DSET_BRICK_TYPE(aset,ii) ,
373                              MRI_type_string(DSET_BRICK_TYPE(aset,ii)) );
374             ++err;
375       }
376    }
377    DSET_KILL_STATS(aset); THD_load_statistics(aset);
378    if (err > 0) {
379       ERROR_message("THD_dset_scale: A total of %d sub-bricks were not scaled", err);
380    }
381 
382    RETURN(err);
383 }
384 
385 /*--------------------------------------------------------------*/
386 /* Number of sub-bricks with nonzero data somewhere inside.
387    If this returns zero, the dataset is all zero. [17 Jan 2017]
388 *//*------------------------------------------------------------*/
389 
THD_count_nonzero_bricks(THD_3dim_dataset * dset)390 int THD_count_nonzero_bricks( THD_3dim_dataset *dset )
391 {
392    int nzc=0 , iv , nvals ;
393 
394    if( !ISVALID_DSET(dset) ) return 0 ;
395 
396    DSET_load(dset) ;
397    iv = THD_count_databricks( dset->dblk ) ;
398    if( iv == 0 ) return 0 ;
399 
400    nvals = DSET_NVALS(dset) ;
401    for( iv=0 ; iv < nvals ; iv++ ){
402      nzc += ( mri_allzero(DSET_BRICK(dset,iv)) == 0 ) ;
403    }
404 
405    return nzc ;
406 }
407