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