1 #include "mrilib.h"
2
3 /*-----------------------------------------------------------------*/
4 /*! Compute median brick of a dataset - 12 Aug 2001
5 - 05 Nov 2001: Modified to use THD_extract_array() instead
6 of THD_extract_series()
7 -------------------------------------------------------------------*/
8
THD_median_brick(THD_3dim_dataset * dset)9 MRI_IMAGE * THD_median_brick( THD_3dim_dataset *dset )
10 {
11 int nvox , nvals , ii ;
12 MRI_IMAGE *tsim , *medim ;
13 float *medar ;
14 float *tsar ; /* 05 Nov 2001 */
15
16 ENTRY("THD_median_brick") ;
17
18 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
19 DSET_load(dset) ;
20 if( !DSET_LOADED(dset) ) RETURN(NULL) ;
21
22 nvals = DSET_NVALS(dset) ;
23 tsim = DSET_BRICK(dset,0) ;
24
25 if( nvals == 1 ){
26 medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;
27 RETURN(medim) ;
28 }
29
30 medim = mri_new_conforming( tsim , MRI_float ) ;
31 medar = MRI_FLOAT_PTR(medim) ;
32 nvox = DSET_NVOX(dset) ;
33
34 tsar = (float *) calloc( sizeof(float),nvals+1 ) ; /* 05 Nov 2001 */
35 for( ii=0 ; ii < nvox ; ii++ ){
36 THD_extract_array( ii , dset , 0 , tsar ) ; /* 05 Nov 2001 */
37 medar[ii] = qmed_float( nvals , tsar ) ;
38 }
39
40 free(tsar) ; RETURN(medim) ;
41 }
42
43 /*-----------------------------------------------------------------*/
44 /*! Compute MAD brick of a dataset - 07 Dec 2006
45 -------------------------------------------------------------------*/
46
THD_mad_brick(THD_3dim_dataset * dset)47 MRI_IMAGE * THD_mad_brick( THD_3dim_dataset *dset )
48 {
49 int nvox , nvals , ii ;
50 MRI_IMAGE *tsim , *madim ;
51 float *madar ;
52 float *tsar ;
53
54 ENTRY("THD_mad_brick") ;
55
56 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
57
58 nvals = DSET_NVALS(dset) ; if( nvals == 1 ) RETURN(NULL) ;
59
60 DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;
61 tsim = DSET_BRICK(dset,0) ;
62
63 madim = mri_new_conforming( tsim , MRI_float ) ;
64 madar = MRI_FLOAT_PTR(madim) ;
65 nvox = DSET_NVOX(dset) ;
66
67 tsar = (float *) calloc( sizeof(float),nvals+1 ) ;
68 for( ii=0 ; ii < nvox ; ii++ ){
69 THD_extract_array( ii , dset , 0 , tsar ) ;
70 qmedmad_float( nvals , tsar , NULL , madar+ii ) ;
71 }
72
73 free(tsar) ; RETURN(madim) ;
74 }
75
76 /*-----------------------------------------------------------------*/
77 /*! Compute median and MAD bricks of a dataset - 07 Dec 2006
78 -------------------------------------------------------------------*/
79
THD_medmad_bricks(THD_3dim_dataset * dset)80 MRI_IMARR * THD_medmad_bricks( THD_3dim_dataset *dset )
81 {
82 int nvox , nvals , ii ;
83 MRI_IMAGE *tsim , *madim, *medim ;
84 float *madar, *medar ;
85 MRI_IMARR *imar ;
86 float *tsar ;
87
88 ENTRY("THD_medmad_bricks") ;
89
90 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
91
92 nvals = DSET_NVALS(dset) ; if( nvals == 1 ) RETURN(NULL) ;
93
94 DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;
95 tsim = DSET_BRICK(dset,0) ;
96
97 madim = mri_new_conforming( tsim , MRI_float ) ;
98 madar = MRI_FLOAT_PTR(madim) ;
99 medim = mri_new_conforming( tsim , MRI_float ) ;
100 medar = MRI_FLOAT_PTR(medim) ;
101 nvox = DSET_NVOX(dset) ;
102
103 tsar = (float *) calloc( sizeof(float),nvals+1 ) ;
104 for( ii=0 ; ii < nvox ; ii++ ){
105 THD_extract_array( ii , dset , 0 , tsar ) ;
106 qmedmad_float( nvals , tsar , medar+ii , madar+ii ) ;
107 }
108
109 free(tsar) ;
110 INIT_IMARR(imar) ; ADDTO_IMARR(imar,medim) ; ADDTO_IMARR(imar,madim) ;
111 RETURN(imar) ;
112 }
113
114 /*-----------------------------------------------------------------*/
115 /*! Compute median and MAD bricks of an image array
116 -------------------------------------------------------------------*/
117
IMARR_medmad_bricks(MRI_IMARR * dmar)118 MRI_IMARR * IMARR_medmad_bricks( MRI_IMARR *dmar )
119 {
120 int nvox , nvals , ii , kk ;
121 MRI_IMAGE *tsim , *madim, *medim ;
122 float *madar, *medar ;
123 MRI_IMARR *imar ;
124 float *tsar , **dar ;
125
126 ENTRY("IMARR_medmad_bricks") ;
127
128 if( dmar == NULL || IMARR_COUNT(dmar) < 2 ) RETURN(NULL) ;
129
130 nvals = IMARR_COUNT(dmar) ;
131 tsim = IMARR_SUBIM(dmar,0) ;
132 madim = mri_new_conforming( tsim , MRI_float ) ;
133 madar = MRI_FLOAT_PTR(madim) ;
134 medim = mri_new_conforming( tsim , MRI_float ) ;
135 medar = MRI_FLOAT_PTR(medim) ;
136 nvox = tsim->nvox ;
137
138 dar = (float **)malloc(sizeof(float *)*nvals) ;
139 for( kk=0 ; kk < nvals ; kk++ )
140 dar[kk] = MRI_FLOAT_PTR( IMARR_SUBIM(dmar,kk) ) ;
141
142 tsar = (float *) calloc( sizeof(float),nvals+1 ) ;
143 for( ii=0 ; ii < nvox ; ii++ ){
144 for( kk=0 ; kk < nvals ; kk++ ) tsar[kk] = dar[kk][ii] ;
145 qmedmad_float( nvals , tsar , medar+ii , madar+ii ) ;
146 }
147
148 free(tsar) ; free(dar) ;
149 INIT_IMARR(imar) ; ADDTO_IMARR(imar,medim) ; ADDTO_IMARR(imar,madim) ;
150 RETURN(imar) ;
151 }
152
153 /*-----------------------------------------------------------------*/
154 /*! Compute mean and sigma bricks of a dataset - 07 Dec 2006
155 -------------------------------------------------------------------*/
156
THD_meansigma_bricks(THD_3dim_dataset * dset)157 MRI_IMARR * THD_meansigma_bricks( THD_3dim_dataset *dset )
158 {
159 int nvox , nvals , ii ;
160 MRI_IMAGE *tsim , *sigim, *mnnim ;
161 float *sigar, *mnnar ;
162 MRI_IMARR *imar ;
163 float *tsar ;
164
165 ENTRY("THD_meansigma_bricks") ;
166
167 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
168
169 nvals = DSET_NVALS(dset) ; if( nvals == 1 ) RETURN(NULL) ;
170
171 DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;
172 tsim = DSET_BRICK(dset,0) ;
173
174 sigim = mri_new_conforming( tsim , MRI_float ) ;
175 sigar = MRI_FLOAT_PTR(sigim) ;
176 mnnim = mri_new_conforming( tsim , MRI_float ) ;
177 mnnar = MRI_FLOAT_PTR(mnnim) ;
178 nvox = DSET_NVOX(dset) ;
179
180 tsar = (float *) calloc( sizeof(float),nvals+1 ) ;
181 for( ii=0 ; ii < nvox ; ii++ ){
182 THD_extract_array( ii , dset , 0 , tsar ) ;
183 meansigma_float( nvals , tsar , mnnar+ii , sigar+ii ) ;
184 }
185
186 free(tsar) ;
187 INIT_IMARR(imar) ; ADDTO_IMARR(imar,mnnim) ; ADDTO_IMARR(imar,sigim) ;
188 RETURN(imar) ;
189 }
190
191 /*-----------------------------------------------------------------*/
192 /*! Compute mean brick of a dataset. [15 Apr 2005 - RWCox]
193 -------------------------------------------------------------------*/
194
THD_mean_brick(THD_3dim_dataset * dset)195 MRI_IMAGE * THD_mean_brick( THD_3dim_dataset *dset )
196 {
197 int nvox , nvals , ii , jj ;
198 MRI_IMAGE *tsim , *medim ;
199 float *medar , sum,fac ;
200 float *tsar ;
201
202 ENTRY("THD_mean_brick") ;
203
204 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
205 DSET_load(dset) ;
206 if( !DSET_LOADED(dset) ) RETURN(NULL) ;
207
208 nvals = DSET_NVALS(dset) ; fac = 1.0 / nvals ;
209 tsim = DSET_BRICK(dset,0) ;
210
211 if( nvals == 1 ){
212 medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;
213 RETURN(medim) ;
214 }
215
216 medim = mri_new_conforming( tsim , MRI_float ) ;
217 medar = MRI_FLOAT_PTR(medim) ;
218 nvox = DSET_NVOX(dset) ;
219
220 tsar = (float *) calloc( sizeof(float),nvals+1 ) ;
221 for( ii=0 ; ii < nvox ; ii++ ){
222 THD_extract_array( ii , dset , 0 , tsar ) ;
223 for( sum=0.0,jj=0 ; jj < nvals ; jj++ ) sum += tsar[jj] ;
224 medar[ii] = fac * sum ;
225 }
226
227 free(tsar) ; RETURN(medim) ;
228 }
229
230 /*-----------------------------------------------------------------*/
231 /*! Compute RMS brick of a dataset. [15 Apr 2005 - RWCox]
232 -------------------------------------------------------------------*/
233
THD_rms_brick(THD_3dim_dataset * dset)234 MRI_IMAGE * THD_rms_brick( THD_3dim_dataset *dset )
235 {
236 int nvox , nvals , ii , jj ;
237 MRI_IMAGE *tsim , *medim ;
238 float *medar , sum,fac ;
239 float *tsar ;
240
241 ENTRY("THD_rms_brick") ;
242
243 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
244 DSET_load(dset) ;
245 if( !DSET_LOADED(dset) ) RETURN(NULL) ;
246
247 nvals = DSET_NVALS(dset) ; fac = 1.0 / nvals ;
248 tsim = DSET_BRICK(dset,0) ;
249
250 if( nvals == 1 ){
251 medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;
252 RETURN(medim) ;
253 }
254
255 medim = mri_new_conforming( tsim , MRI_float ) ;
256 medar = MRI_FLOAT_PTR(medim) ;
257 nvox = DSET_NVOX(dset) ;
258
259 tsar = (float *) calloc( sizeof(float),nvals+1 ) ;
260 for( ii=0 ; ii < nvox ; ii++ ){
261 THD_extract_array( ii , dset , 0 , tsar ) ;
262 for( sum=0.0,jj=0 ; jj < nvals ; jj++ ) sum += tsar[jj]*tsar[jj] ;
263 medar[ii] = sqrtf(fac * sum) ;
264 }
265
266 free(tsar) ; RETURN(medim) ;
267 }
268
269 /*-----------------------------------------------------------------*/
270 /*! Compute average abs brick of a dataset. [11 May 2009 - RWCox]
271 -------------------------------------------------------------------*/
272
THD_aveabs_brick(THD_3dim_dataset * dset)273 MRI_IMAGE * THD_aveabs_brick( THD_3dim_dataset *dset )
274 {
275 int nvox , nvals , ii , jj ;
276 MRI_IMAGE *tsim , *medim ;
277 float *medar , sum,fac ;
278 float *tsar ;
279
280 ENTRY("THD_aveabs_brick") ;
281
282 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
283 DSET_load(dset) ;
284 if( !DSET_LOADED(dset) ) RETURN(NULL) ;
285
286 nvals = DSET_NVALS(dset) ; fac = 1.0 / nvals ;
287 tsim = DSET_BRICK(dset,0) ;
288 nvox = DSET_NVOX(dset) ;
289
290 if( nvals == 1 ){
291 medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;
292 medar = MRI_FLOAT_PTR(medim) ;
293 for( ii=0 ; ii < nvox ; ii++ ) medar[ii] = fabsf(medar[ii]) ;
294 RETURN(medim) ;
295 }
296
297 medim = mri_new_conforming( tsim , MRI_float ) ;
298 medar = MRI_FLOAT_PTR(medim) ;
299 tsar = (float *) calloc( sizeof(float),nvals+1 ) ;
300
301 for( ii=0 ; ii < nvox ; ii++ ){
302 THD_extract_array( ii , dset , 0 , tsar ) ;
303 for( sum=0.0,jj=0 ; jj < nvals ; jj++ ) sum += fabsf(tsar[jj]) ;
304 medar[ii] = fac * sum ;
305 }
306
307 free(tsar) ; RETURN(medim) ;
308 }
309
310 /*-----------------------------------------------------------------*/
311 /*! Compute max abs brick of a dataset. [08 Jan 2019 - RWCox]
312 -------------------------------------------------------------------*/
313
THD_maxabs_brick(THD_3dim_dataset * dset)314 MRI_IMAGE * THD_maxabs_brick( THD_3dim_dataset *dset )
315 {
316 int nvox , nvals , ii , jj ;
317 MRI_IMAGE *tsim , *medim ;
318 float *medar , sum,val ;
319 float *tsar ;
320
321 ENTRY("THD_maxabs_brick") ;
322
323 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
324 DSET_load(dset) ;
325 if( !DSET_LOADED(dset) ) RETURN(NULL) ;
326
327 nvals = DSET_NVALS(dset) ;
328 tsim = DSET_BRICK(dset,0) ;
329 nvox = DSET_NVOX(dset) ;
330
331 if( nvals == 1 ){
332 medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;
333 medar = MRI_FLOAT_PTR(medim) ;
334 for( ii=0 ; ii < nvox ; ii++ ) medar[ii] = fabsf(medar[ii]) ;
335 RETURN(medim) ;
336 }
337
338 medim = mri_new_conforming( tsim , MRI_float ) ;
339 medar = MRI_FLOAT_PTR(medim) ;
340 tsar = (float *) calloc( sizeof(float),nvals+1 ) ;
341
342 for( ii=0 ; ii < nvox ; ii++ ){
343 THD_extract_array( ii , dset , 0 , tsar ) ;
344 for( sum=0.0,jj=0 ; jj < nvals ; jj++ ){
345 val = fabsf(tsar[jj]) ; if( sum < val ) sum = val ;
346 }
347 medar[ii] = sum ;
348 }
349
350 free(tsar) ; RETURN(medim) ;
351 }
352
353 /*-------------------------------------------------------------------*/
354 /*! Compute average positive brick of a dataset. [09 Nov 2020 - RWCox]
355 ---------------------------------------------------------------------*/
356
THD_avepos_brick(THD_3dim_dataset * dset)357 MRI_IMAGE * THD_avepos_brick( THD_3dim_dataset *dset )
358 {
359 int nvox , nvals , ii , jj ;
360 MRI_IMAGE *tsim , *medim ;
361 float *medar , sum,fac ;
362 float *tsar ;
363
364 ENTRY("THD_avepos_brick") ;
365
366 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
367 DSET_load(dset) ;
368 if( !DSET_LOADED(dset) ) RETURN(NULL) ;
369
370 nvals = DSET_NVALS(dset) ; fac = 1.0 / nvals ;
371 tsim = DSET_BRICK(dset,0) ;
372 nvox = DSET_NVOX(dset) ;
373
374 if( nvals == 1 ){
375 medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;
376 medar = MRI_FLOAT_PTR(medim) ;
377 for( ii=0 ; ii < nvox ; ii++ ){
378 if( medar[ii] < 0.0f ) medar[ii] = 0.0f ;
379 }
380 RETURN(medim) ;
381 }
382
383 medim = mri_new_conforming( tsim , MRI_float ) ;
384 medar = MRI_FLOAT_PTR(medim) ;
385 tsar = (float *) calloc( sizeof(float),nvals+1 ) ;
386
387 for( ii=0 ; ii < nvox ; ii++ ){
388 THD_extract_array( ii , dset , 0 , tsar ) ;
389 for( sum=0.0,jj=0 ; jj < nvals ; jj++ ){
390 if( tsar[jj] > 0.0f ) sum += tsar[jj] ;
391 }
392 medar[ii] = fac * sum ;
393 }
394
395 free(tsar) ; RETURN(medim) ;
396 }
397