1 #include "mrilib.h"
2
3 #ifdef USE_OMP
4 #include <omp.h>
5 #endif
6
7 #undef VECTIM_scan
8 #define VECTIM_scan(vv) \
9 do{ size_t nbad = thd_floatscan((size_t)(vv)->nvals*(size_t)(vv)->nvec,(vv)->fvec) ; \
10 if( nbad > 0 ) WARNING_message("found %lld bad values in vectim",(long long)nbad); \
11 } while(0)
12
13 /*--------------------------------------------------------------------------*/
14 /*! Convert a dataset to the MRI_vectim format, where each time
15 series is a contiguous set of values in an array, and the
16 voxel indexes whence came the values are also stored.
17 *//*------------------------------------------------------------------------*/
18
THD_dset_to_vectim(THD_3dim_dataset * dset,byte * mask,int ignore)19 MRI_vectim * THD_dset_to_vectim( THD_3dim_dataset *dset, byte *mask , int ignore )
20 {
21 byte *mmm=mask ;
22 MRI_vectim *mrv=NULL ;
23 int kk,iv , nvals , nvox , nmask ;
24 size_t ntot ;
25
26 ENTRY("THD_dset_to_vectim") ;
27
28 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
29 DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;
30
31 if( ignore < 0 ) ignore = 0 ;
32 nvals = DSET_NVALS(dset) - ignore ; if( nvals <= 0 ) RETURN(NULL) ;
33 nvox = DSET_NVOX(dset) ;
34
35 if( mmm != NULL ){
36 nmask = THD_countmask( nvox , mmm ) ; /* number to keep */
37 if( nmask <= 0 ) RETURN(NULL) ;
38 } else {
39 nmask = nvox ; /* keep them all */
40 #pragma omp critical (MALLOC)
41 mmm = (byte *)malloc(sizeof(byte)*nmask) ;
42 if( mmm == NULL ){
43 ERROR_message("THD_dset_to_vectim: out of memory A") ;
44 RETURN(NULL) ;
45 }
46 memset( mmm , 1 , sizeof(byte)*nmask ) ;
47 }
48
49 #pragma omp critical (MALLOC)
50 mrv = (MRI_vectim *)malloc(sizeof(MRI_vectim)) ;
51
52 mrv->nvec = nmask ;
53 mrv->nvals = nvals ;
54 mrv->ignore = ignore ;
55 #pragma omp critical (MALLOC)
56 mrv->ivec = (int *)malloc(sizeof(int)*nmask) ;
57 if( mrv->ivec == NULL ){
58 ERROR_message("THD_dset_to_vectim: out of memory") ;
59 free(mrv) ; if( mmm != mask ) free(mmm) ;
60 RETURN(NULL) ;
61 }
62 ntot = sizeof(float)*(size_t)nmask*(size_t)nvals ;
63 #pragma omp critical (MALLOC)
64 mrv->fvec = (float *)malloc(ntot) ;
65 if( mrv->fvec == NULL ){
66 ERROR_message("THD_dset_to_vectim: out of memory B -- tried to get %lld bytes",(long long)ntot) ;
67 free(mrv->ivec) ; free(mrv) ; if( mmm != mask ) free(mmm) ;
68 RETURN(NULL) ;
69 } else if( ntot > 1000000000 ){ /* one BILLION bytes */
70 INFO_message("THD_dset_to_vectim: allocated %lld bytes",(long long)ntot) ;
71 }
72
73 /* store desired voxel time series */
74
75 STATUS("create index list") ;
76 for( kk=iv=0 ; iv < nvox ; iv++ ){
77 if( mmm[iv] ) mrv->ivec[kk++] = iv ; /* build index list */
78 }
79
80 if( ignore > 0 ){ /* extract 1 at a time, save what we want */
81
82 float *var;
83 #pragma omp critical (MALLOC)
84 var = (float *)malloc(sizeof(float)*(nvals+ignore)) ;
85 STATUS("ignore > 0 --> extracting one at a time") ;
86 for( kk=iv=0 ; iv < nvox ; iv++ ){
87 if( mmm[iv] == 0 ) continue ;
88 (void)THD_extract_array( iv , dset , 0 , var ) ;
89 AAmemcpy( VECTIM_PTR(mrv,kk) , var+ignore , sizeof(float)*nvals ) ;
90 kk++ ;
91 }
92 free(var) ;
93
94 } else { /* do all at once: this way is a lot faster */
95
96 STATUS("ignore==0 --> extracting all at once") ;
97 THD_extract_many_arrays( nmask , mrv->ivec , dset , mrv->fvec ) ;
98
99 }
100
101 STATUS("setting parameters in vectim header") ;
102
103 mrv->nx = DSET_NX(dset) ; mrv->dx = fabs(DSET_DX(dset)) ;
104 mrv->ny = DSET_NY(dset) ; mrv->dy = fabs(DSET_DY(dset)) ;
105 mrv->nz = DSET_NZ(dset) ; mrv->dz = fabs(DSET_DZ(dset)) ;
106
107 DSET_UNMSEC(dset) ; mrv->dt = DSET_TR(dset) ;
108 if( mrv->dt <= 0.0f ) mrv->dt = 1.0f ;
109
110
111 if( mmm != mask ){
112 STATUS("free(mmm)") ;
113 free(mmm) ;
114 }
115
116 STATUS("VECTIM_scan()") ;
117 VECTIM_scan(mrv) ; /* 09 Nov 2010 */
118 RETURN(mrv) ;
119 }
120
121 /*--------------------------------------------------------------------------*/
122
THD_dset_censored_to_vectim(THD_3dim_dataset * dset,byte * mask,int nkeep,int * keep)123 MRI_vectim * THD_dset_censored_to_vectim( THD_3dim_dataset *dset,
124 byte *mask , int nkeep , int *keep )
125 {
126 byte *mmm=mask ;
127 MRI_vectim *mrv=NULL ;
128 int kk,iv,jj , nvals , nvox , nmask ;
129
130 ENTRY("THD_dset_censored_to_vectim") ;
131
132 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
133
134 if( nkeep <= 0 || keep == NULL ){
135 mrv = THD_dset_to_vectim( dset , mask , 0 ) ;
136 RETURN(mrv) ;
137 }
138
139 if( ! THD_subset_loaded(dset,nkeep,keep) ){
140 DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;
141 }
142
143 nvals = nkeep ;
144 nvox = DSET_NVOX(dset) ;
145
146 if( mmm != NULL ){
147 nmask = THD_countmask( nvox , mmm ) ; /* number to keep */
148 if( nmask <= 0 ) RETURN(NULL) ;
149 } else {
150 nmask = nvox ; /* keep them all */
151 #pragma omp critical (MALLOC)
152 mmm = (byte *)malloc(sizeof(byte)*nmask) ;
153 if( mmm == NULL ){
154 ERROR_message("THD_dset_to_vectim: out of memory C") ;
155 RETURN(NULL) ;
156 }
157 memset( mmm , 1 , sizeof(byte)*nmask ) ;
158 }
159
160 #pragma omp critical (MALLOC)
161 mrv = (MRI_vectim *)malloc(sizeof(MRI_vectim)) ;
162
163 mrv->nvec = nmask ;
164 mrv->nvals = nvals ;
165 mrv->ignore = 0 ;
166 #pragma omp critical (MALLOC)
167 mrv->ivec = (int *)malloc(sizeof(int)*nmask) ;
168 if( mrv->ivec == NULL ){
169 ERROR_message("THD_dset_to_vectim: out of memory D") ;
170 free(mrv) ; if( mmm != mask ) free(mmm) ;
171 RETURN(NULL) ;
172 }
173 #pragma omp critical (MALLOC)
174 mrv->fvec = (float *)malloc(sizeof(float)*(size_t)nmask*(size_t)nvals) ;
175 if( mrv->fvec == NULL ){
176 ERROR_message("THD_dset_to_vectim: out of memory E") ;
177 free(mrv->ivec) ; free(mrv) ; if( mmm != mask ) free(mmm) ;
178 RETURN(NULL) ;
179 }
180
181 /* store desired voxel time series */
182
183 for( kk=iv=0 ; iv < nvox ; iv++ ){
184 if( mmm[iv] ) mrv->ivec[kk++] = iv ; /* build index list */
185 }
186
187 #pragma omp critical (MALLOC)
188 { float *var = (float *)malloc(sizeof(float)*DSET_NVALS(dset)) ;
189 float *vpt ;
190 for( kk=iv=0 ; iv < nvox ; iv++ ){
191 if( mmm[iv] == 0 ) continue ;
192 vpt = VECTIM_PTR(mrv,kk) ; kk++ ;
193 if( nkeep > 1 ){
194 (void)THD_extract_array( iv , dset , 0 , var ) ;
195 for( jj=0 ; jj < nkeep ; jj++ ) vpt[jj] = var[keep[jj]] ;
196 } else {
197 vpt[0] = THD_get_float_value(iv,keep[0],dset) ;
198 }
199 }
200 free(var) ;
201 }
202
203 mrv->nx = DSET_NX(dset) ; mrv->dx = fabs(DSET_DX(dset)) ;
204 mrv->ny = DSET_NY(dset) ; mrv->dy = fabs(DSET_DY(dset)) ;
205 mrv->nz = DSET_NZ(dset) ; mrv->dz = fabs(DSET_DZ(dset)) ;
206
207 DSET_UNMSEC(dset) ; mrv->dt = DSET_TR(dset) ;
208 if( mrv->dt <= 0.0f ) mrv->dt = 1.0f ;
209
210 if( mmm != mask ) free(mmm) ;
211 VECTIM_scan(mrv) ; /* 09 Nov 2010 */
212 RETURN(mrv) ;
213 }
214
215 /*---------------------------------------------------------------------------*/
216
THD_vectim_data_tofile(MRI_vectim * mrv,char * fnam)217 int THD_vectim_data_tofile( MRI_vectim *mrv , char *fnam )
218 {
219 FILE *fp ; size_t nf , nw ;
220
221 if( mrv == NULL || fnam == NULL ) return 0 ;
222
223 fp = fopen( fnam , "w" ) ; if( fp == NULL ) return 0 ;
224 nf = ((size_t)mrv->nvec) * ((size_t)mrv->nvals) ;
225 nw = fwrite( mrv->fvec , sizeof(float) , nf , fp ) ;
226 fclose(fp) ; if( nw == nf ) return 1 ;
227 remove(fnam) ; return 0 ;
228 }
229
230 /*---------------------------------------------------------------------------*/
231
THD_vector_fromfile(int nvals,int iv,float * vv,FILE * fp)232 void THD_vector_fromfile( int nvals , int iv , float *vv , FILE *fp )
233 {
234 fseeko( fp , ((off_t)iv) * ((off_t)nvals) * ((off_t)sizeof(float)) , SEEK_SET ) ;
235 fread( vv , sizeof(float) , nvals , fp ) ;
236 return ;
237 }
238
239 /*---------------------------------------------------------------------------*/
240
THD_vectim_reload_fromfile(MRI_vectim * mrv,char * fnam)241 int THD_vectim_reload_fromfile( MRI_vectim *mrv , char *fnam )
242 {
243 FILE *fp ; size_t nf , nw ;
244
245 if( mrv == NULL || fnam == NULL ) return 0 ;
246
247 fp = fopen( fnam , "r" ) ; if( fp == NULL ) return 0 ;
248 nf = (size_t)mrv->nvec * (size_t)mrv->nvals ;
249 if( mrv->fvec == NULL ) mrv->fvec = (float *)malloc(sizeof(float)*nf) ;
250 nw = fread( mrv->fvec , sizeof(float) , nf , fp ) ;
251 fclose(fp) ; return (int)nw ;
252 }
253
254 /*---------------------------------------------------------------------------*/
255
THD_dset_to_vectim_stend(THD_3dim_dataset * dset,byte * mask,int start,int end)256 MRI_vectim * THD_dset_to_vectim_stend( THD_3dim_dataset *dset, byte *mask , int start, int end )
257 {
258 byte *mmm=mask ;
259 MRI_vectim *mrv=NULL ;
260 int kk,iv , nvals , nvox , nmask ;
261
262 ENTRY("THD_dset_to_vectim_stend") ;
263
264 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
265 DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;
266
267 if( start < 0 ) start = 0 ;
268 if( end < start ) end = DSET_NVALS(dset)-1 ;
269 nvals = end - start + 1 ; if( nvals <= 0 ) RETURN(NULL) ;
270 nvox = DSET_NVOX(dset) ;
271
272 if( mmm != NULL ){
273 nmask = THD_countmask( nvox , mmm ) ; /* number to keep */
274 if( nmask <= 0 ) RETURN(NULL) ;
275 } else {
276 nmask = nvox ; /* keep them all */
277 #pragma omp critical (MALLOC)
278 mmm = (byte *)malloc(sizeof(byte)*nmask) ;
279 if( mmm == NULL ){
280 ERROR_message("THD_dset_to_vectim: out of memory F") ;
281 RETURN(NULL) ;
282 }
283 memset( mmm , 1 , sizeof(byte)*nmask ) ;
284 }
285
286 #pragma omp critical (MALLOC)
287 mrv = (MRI_vectim *)malloc(sizeof(MRI_vectim)) ;
288
289 mrv->nvec = nmask ;
290 mrv->nvals = nvals ;
291 mrv->ignore = start ;
292 #pragma omp critical (MALLOC)
293 mrv->ivec = (int *)malloc(sizeof(int)*nmask) ;
294 if( mrv->ivec == NULL ){
295 ERROR_message("THD_dset_to_vectim: out of memory G") ;
296 free(mrv) ; if( mmm != mask ) free(mmm) ;
297 RETURN(NULL) ;
298 }
299 #pragma omp critical (MALLOC)
300 mrv->fvec = (float *)malloc(sizeof(float)*(size_t)nmask*(size_t)nvals) ;
301 if( mrv->fvec == NULL ){
302 ERROR_message("THD_dset_to_vectim: out of memory H") ;
303 free(mrv->ivec) ; free(mrv) ; if( mmm != mask ) free(mmm) ;
304 RETURN(NULL) ;
305 }
306
307 /* store desired voxel time series */
308
309 for( kk=iv=0 ; iv < nvox ; iv++ ){
310 if( mmm[iv] ) mrv->ivec[kk++] = iv ; /* build index list */
311 }
312
313 if( nvals < DSET_NVALS(dset) ){ /* extract 1 at a time, save what we want */
314
315 float *var;
316 #pragma omp critical (MALLOC)
317 var = (float *)malloc(sizeof(float)*(DSET_NVALS(dset))) ;
318 for( kk=iv=0 ; iv < nvox ; iv++ ){
319 if( mmm[iv] == 0 ) continue ;
320 (void)THD_extract_array( iv , dset , 0 , var ) ;
321 AAmemcpy( VECTIM_PTR(mrv,kk) , var+start , sizeof(float)*nvals ) ;
322 kk++ ;
323 }
324 free(var) ;
325
326 } else { /* do all at once: this way is a lot faster */
327
328 THD_extract_many_arrays( nmask , mrv->ivec , dset , mrv->fvec ) ;
329
330 }
331
332 mrv->nx = DSET_NX(dset) ; mrv->dx = fabs(DSET_DX(dset)) ;
333 mrv->ny = DSET_NY(dset) ; mrv->dy = fabs(DSET_DY(dset)) ;
334 mrv->nz = DSET_NZ(dset) ; mrv->dz = fabs(DSET_DZ(dset)) ;
335
336 DSET_UNMSEC(dset) ; mrv->dt = DSET_TR(dset) ;
337 if( mrv->dt <= 0.0f ) mrv->dt = 1.0f ;
338
339 if( mmm != mask ) free(mmm) ;
340 VECTIM_scan(mrv) ; /* 09 Nov 2010 */
341 RETURN(mrv) ;
342 }
343
344 /*--------------------------------------------------------------------------*/
345
THD_dset_to_vectim_byslice(THD_3dim_dataset * dset,byte * mask,int ignore,int kzbot,int kztop)346 MRI_vectim * THD_dset_to_vectim_byslice( THD_3dim_dataset *dset, byte *mask ,
347 int ignore , int kzbot , int kztop )
348 {
349 byte *mmm ;
350 MRI_vectim *mrv=NULL ;
351 int kk,iv , nvals , nvox , nmask , nxy , nz ;
352
353 ENTRY("THD_dset_to_vectim_byslice") ;
354
355 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
356 DSET_load(dset) ; if( !DSET_LOADED(dset) ) RETURN(NULL) ;
357
358 nvals = DSET_NVALS(dset) ; if( nvals <= 0 ) RETURN(NULL) ;
359 nvox = DSET_NVOX(dset) ;
360
361 nxy = DSET_NX(dset) * DSET_NY(dset) ; nz = DSET_NZ(dset) ;
362
363 if( kzbot < 0 ) kzbot = 0 ;
364 if( kztop >= nz ) kztop = nz-1 ;
365 if( kztop < kzbot ) RETURN(NULL) ;
366 if( kzbot == 0 && kztop == nz-1 ){
367 mrv = THD_dset_to_vectim( dset , mask, ignore ) ; RETURN(mrv) ;
368 }
369
370 /* make a mask that includes cutting out un-desirable slices */
371
372 { int ibot , itop , ii ;
373 #pragma omp critical (MALLOC)
374 mmm = (byte *)malloc(sizeof(byte)*nvox) ;
375 if( mask == NULL ) AAmemset( mmm , 1 , sizeof(byte)*nvox ) ;
376 else AAmemcpy( mmm , mask , sizeof(byte)*nvox ) ;
377 if( kzbot > 0 )
378 AAmemset( mmm , 0 , sizeof(byte)*kzbot *nxy ) ;
379 if( kztop < nz-1 )
380 AAmemset( mmm+(kztop+1)*nxy , 0 , sizeof(byte)*(nz-1-kztop)*nxy ) ;
381 }
382
383 /* and make the vectim using the standard function */
384
385 mrv = THD_dset_to_vectim( dset , mmm , ignore ) ;
386 free(mmm) ;
387 RETURN(mrv) ;
388 }
389
390 /*---------------------------------------------------------------------*/
391 /*-------- Catenates two datasets into vectim ZSS Jan 2010 --------*/
392
THD_2dset_to_vectim(THD_3dim_dataset * dset1,byte * mask1,THD_3dim_dataset * dset2,byte * mask2,int ignore)393 MRI_vectim * THD_2dset_to_vectim( THD_3dim_dataset *dset1, byte *mask1 ,
394 THD_3dim_dataset *dset2, byte *mask2 ,
395 int ignore )
396 {
397 byte *mmmv[2]={NULL, NULL}, *mmmt=NULL;
398 THD_3dim_dataset *dsetv[2]={NULL, NULL};
399 MRI_vectim *mrv=NULL ;
400 int kk2, kk,iv,id, nvals , nvoxv[2]={0,0} , nmaskv[2]={0,0} ;
401 int *ivvectmp=NULL;
402
403 ENTRY("THD_2dset_to_vectim") ;
404 mmmv[0] = mask1;
405 mmmv[1] = mask2;
406 dsetv[0] = dset1;
407 dsetv[1] = dset2;
408 for (id=0; id<2;++id) {
409 if( !ISVALID_DSET(dsetv[id]) ) RETURN(NULL) ;
410 DSET_load(dsetv[id]) ; if( !DSET_LOADED(dsetv[id]) ) RETURN(NULL) ;
411 nvoxv[id] = DSET_NVOX(dsetv[id]) ;
412 }
413 if (DSET_NVALS(dsetv[0]) != DSET_NVALS(dsetv[1])) {
414 RETURN(NULL) ;
415 }
416
417 if( ignore < 0 ) ignore = 0 ;
418 nvals = DSET_NVALS(dsetv[0]) - ignore ; if( nvals <= 0 ) RETURN(NULL) ;
419
420 for (id=0; id<2; ++id) {
421 if( mmmv[id] != NULL ){
422 nmaskv[id] = THD_countmask( nvoxv[id] , mmmv[id] ) ;/* number to keep */
423 if( nmaskv[id] <= 0 ) RETURN(NULL) ;
424 } else {
425 nmaskv[id] = nvoxv[id] ; /* keep them all */
426 #pragma omp critical (MALLOC)
427 mmmv[id] = (byte *)malloc(sizeof(byte)*nmaskv[id]) ;
428 if( mmmv[id] == NULL ){
429 ERROR_message("THD_2dset_to_vectim: out of memory I") ;
430 RETURN(NULL) ;
431 }
432 memset( mmmv[id] , 1 , sizeof(byte)*nmaskv[id] ) ;
433 }
434 }
435
436 #pragma omp critical (MALLOC)
437 mrv = (MRI_vectim *)malloc(sizeof(MRI_vectim)) ;
438
439 mrv->nvec = nmaskv[0]+nmaskv[1] ;
440 mrv->nvals = nvals ;
441 mrv->ignore = ignore ;
442 #pragma omp critical (MALLOC)
443 mrv->ivec = (int *)malloc(sizeof(int)*(nmaskv[0]+nmaskv[1])) ;
444 #pragma omp critical (MALLOC)
445 ivvectmp = (int *)malloc(sizeof(int)*(nmaskv[1])) ;
446 if( mrv->ivec == NULL || ivvectmp == NULL){
447 ERROR_message("THD_2dset_to_vectim: out of memory J") ;
448 if (mrv->ivec) free(mrv->ivec) ;
449 if (ivvectmp) free(ivvectmp) ;
450 free(mrv) ;
451 if( mmmv[0] != mask1 ) free(mmmv[0]) ;
452 if( mmmv[1] != mask2 ) free(mmmv[1]) ;
453 RETURN(NULL) ;
454 }
455 #pragma omp critical (MALLOC)
456 mrv->fvec = (float *)malloc(sizeof(float)*(nmaskv[0]+nmaskv[1])*(size_t)nvals) ;
457 if( mrv->fvec == NULL ){
458 ERROR_message("THD_2dset_to_vectim: out of memory K") ;
459 if (ivvectmp) free(ivvectmp) ;
460 free(mrv->ivec) ; free(mrv) ;
461 if( mmmv[0] != mask1 ) free(mmmv[0]) ;
462 if( mmmv[1] != mask2 ) free(mmmv[1]) ;
463 RETURN(NULL) ;
464 }
465
466 /* store desired voxel time series */
467
468 mmmt = mmmv[0];
469 for( kk=iv=0 ; iv < nvoxv[0] ; iv++ ){
470 if( mmmt[iv] ) mrv->ivec[kk++] = iv ; /* build index list to 1st dset */
471 }
472 mmmt = mmmv[1]; kk2 = 0;
473 for( iv=0 ; iv < nvoxv[1] ; iv++ ){
474 if( mmmt[iv] ) {
475 mrv->ivec[kk++] = iv + nvoxv[0] ;
476 /* build index list to 2nd dset*/
477 ivvectmp[kk2++] = iv;
478 }
479 }
480
481 if( ignore > 0 ){ /* extract 1 at a time, save what we want */
482
483 float *var;
484 #pragma omp critical (MALLOC)
485 var = (float *)malloc(sizeof(float)*(nvals+ignore)) ;
486 mmmt = mmmv[0];
487 for( kk=iv=0 ; iv < nvoxv[0] ; iv++ ){
488 if( mmmt[iv] == 0 ) continue ;
489 (void)THD_extract_array( iv , dsetv[0] , 0 , var ) ;
490 AAmemcpy( VECTIM_PTR(mrv,kk) , var+ignore , sizeof(float)*nvals ) ;
491 kk++ ;
492 }
493 mmmt = mmmv[1];
494 for( iv=0 ; iv < nvoxv[1] ; iv++ ){
495 if( mmmt[iv] == 0 ) continue ;
496 (void)THD_extract_array( iv , dsetv[1] , 0 , var ) ;
497 AAmemcpy( VECTIM_PTR(mrv,kk) , var+ignore , sizeof(float)*nvals ) ;
498 kk++ ;
499 }
500
501 free(var) ;
502
503 } else { /* do all at once: this way is a lot faster */
504
505 THD_extract_many_arrays( nmaskv[0] , mrv->ivec ,
506 dsetv[0] , mrv->fvec ) ;
507 THD_extract_many_arrays( nmaskv[1] , ivvectmp,
508 dsetv[1] , (mrv->fvec+(size_t)nmaskv[0]*(size_t)mrv->nvals) ) ;
509
510 }
511
512 mrv->nx = DSET_NX(dsetv[0]) + DSET_NX(dsetv[1]);
513 mrv->dx = fabs(DSET_DX(dsetv[0])) ;
514 mrv->ny = DSET_NY(dsetv[0]) ; mrv->dy = fabs(DSET_DY(dsetv[0])) ;
515 mrv->nz = DSET_NZ(dsetv[0]) ; mrv->dz = fabs(DSET_DZ(dsetv[0])) ;
516
517 DSET_UNMSEC(dsetv[0]) ; mrv->dt = DSET_TR(dsetv[0]) ;
518 if( mrv->dt <= 0.0f ) mrv->dt = 1.0f ;
519
520 if( mmmv[0] != mask1 ) free(mmmv[0]) ;
521 if( mmmv[1] != mask2 ) free(mmmv[1]) ;
522 if (ivvectmp) free(ivvectmp) ;
523
524 if (0) {
525 int ShowThisTs=38001;
526 float *fff=NULL;
527 fprintf(stderr,"++ ZSS mrv->nvec = %d, mrv->nvals = %d\n",
528 mrv->nvec, mrv->nvals);
529 for( kk=0 ; kk < mrv->nvals; ++kk) {
530 fff=mrv->fvec+((size_t)mrv->nvals*(size_t)ShowThisTs);
531 fprintf(stderr," %f \t", *(fff+kk));
532 }
533 }
534
535 VECTIM_scan(mrv) ; /* 09 Nov 2010 */
536 RETURN(mrv) ;
537 }
538
539 /*---------------------------------------------------------------------*/
540
THD_vectim_copy(MRI_vectim * mrv)541 MRI_vectim * THD_vectim_copy( MRI_vectim *mrv ) /* 08 Apr 2010 */
542 {
543 MRI_vectim *qrv ;
544
545 ENTRY("THD_vectim_copy") ;
546
547 if( mrv == NULL ) RETURN(NULL) ;
548
549 MAKE_VECTIM( qrv , mrv->nvec , mrv->nvals ) ;
550 qrv->ignore = mrv->ignore ;
551 AAmemcpy( qrv->ivec , mrv->ivec , sizeof(int)*mrv->nvec ) ;
552 AAmemcpy( qrv->fvec , mrv->fvec , sizeof(float)*(size_t)mrv->nvec*(size_t)mrv->nvals ) ;
553 qrv->nx = mrv->nx ; qrv->dx = mrv->dx ;
554 qrv->ny = mrv->ny ; qrv->dy = mrv->dy ;
555 qrv->nz = mrv->nz ; qrv->dz = mrv->dz ; qrv->dt = mrv->dt ;
556 RETURN(qrv) ;
557 }
558
559 /*---------------------------------------------------------------------*/
560
THD_vectim_copy_nonzero(MRI_vectim * mrv)561 MRI_vectim * THD_vectim_copy_nonzero( MRI_vectim *mrv ) /* 21 Sep 2010 */
562 {
563 MRI_vectim *qrv ;
564 int nvals , nvec , ii,jj , ngood ;
565 float *mar , *qar ;
566
567 ENTRY("THD_vectim_copy_nonzero") ;
568
569 if( mrv == NULL ) RETURN(NULL) ;
570 nvals = mrv->nvals ; nvec = mrv->nvec ;
571
572 for( ngood=ii=0 ; ii < nvec ; ii++ ){
573 mar = VECTIM_PTR(mrv,ii) ;
574 for( jj=0 ; jj < nvals && mar[jj] == 0.0f ; jj++ ) ; /*nada*/
575 if( jj < nvals ) ngood++ ;
576 }
577 if( ngood == 0 ) RETURN(NULL) ; /* nothing to do */
578 if( ngood == nvec ) RETURN( THD_vectim_copy(mrv) ) ; /* everything to do */
579
580 MAKE_VECTIM( qrv , ngood , nvals ) ;
581 qrv->ignore = mrv->ignore ;
582 for( ngood=ii=0 ; ii < nvec ; ii++ ){
583 mar = VECTIM_PTR(mrv,ii) ;
584 for( jj=0 ; jj < nvals && mar[jj] == 0.0f ; jj++ ) ; /*nada*/
585 if( jj < nvals ){
586 qrv->ivec[ngood] = mrv->ivec[ii] ;
587 qar = VECTIM_PTR(qrv,ngood) ;
588 AAmemcpy(qar,mar,sizeof(float)*nvals) ;
589 ngood++ ;
590 }
591 }
592
593 qrv->nx = mrv->nx ; qrv->dx = mrv->dx ;
594 qrv->ny = mrv->ny ; qrv->dy = mrv->dy ;
595 qrv->nz = mrv->nz ; qrv->dz = mrv->dz ; qrv->dt = mrv->dt ;
596 RETURN(qrv) ;
597 }
598
599 /*---------------------------------------------------------------------*/
600 /* Apply a function to each time series, presumably to edit
601 it in-place. For 3dGroupInCorr -- 10 May 2012 -- RWCox.
602 *//*-------------------------------------------------------------------*/
603
THD_vectim_applyfunc(MRI_vectim * mrv,void * vp)604 void THD_vectim_applyfunc( MRI_vectim *mrv , void *vp )
605 {
606 int iv , ii ;
607 void (*fp)(int,float *) = (void (*)(int,float *))(vp) ;
608
609 if( mrv == NULL || vp == NULL ) return ;
610
611 for( iv=0 ; iv < mrv->nvec ; iv++ ) fp( mrv->nvals, VECTIM_PTR(mrv,iv) ) ;
612 }
613
614 /*---------------------------------------------------------------------*/
615
THD_vectim_normalize(MRI_vectim * mrv)616 void THD_vectim_normalize( MRI_vectim *mrv )
617 {
618 int iv , ii ;
619
620 if( mrv == NULL ) return ;
621
622 for( iv=0 ; iv < mrv->nvec ; iv++ )
623 THD_normalize( mrv->nvals , VECTIM_PTR(mrv,iv) ) ;
624 }
625
626 /*---------------------------------------------------------------------*/
627
THD_vectim_dotprod(MRI_vectim * mrv,float * vec,float * dp,int ata)628 void THD_vectim_dotprod( MRI_vectim *mrv , float *vec , float *dp , int ata )
629 {
630 if( mrv == NULL || vec == NULL || dp == NULL ) return ;
631
632 AFNI_OMP_START ;
633 #pragma omp parallel if( mrv->nvec > 1 && mrv->nvec * mrv->nvals > 999999 )
634 { int nvec=mrv->nvec, nvals=mrv->nvals, nv1=nvals-1, iv, ii ; float sum, *fv ;
635 #pragma omp for
636 for( iv=0 ; iv < nvec ; iv++ ){
637 fv = VECTIM_PTR(mrv,iv) ;
638 for( sum=0.0f,ii=0 ; ii < nv1 ; ii+=2 )
639 sum += fv[ii]*vec[ii] + fv[ii+1]*vec[ii+1] ;
640 if( ii == nv1 ) sum += fv[ii]*vec[ii] ;
641 dp[iv] = (ata) ? logf((1.0001f+sum)/(1.0001f-sum)) : sum ;
642 }
643 } /* end OpenMP */
644 AFNI_OMP_END ;
645
646 thd_floatscan(mrv->nvec,dp) ; /* 16 May 2013 */
647
648 return ;
649 }
650
651 /*---------------------------------------------------------------------*/
652
THD_vectim_vectim_dot(MRI_vectim * arv,MRI_vectim * brv,float * dp)653 void THD_vectim_vectim_dot( MRI_vectim *arv, MRI_vectim *brv, float *dp )
654 {
655 int nvec , nvals , iv , ii ; float *av , *bv , sum ;
656
657 if( arv == NULL || brv == NULL || dp == NULL ) return ;
658 if( arv->nvec != brv->nvec || arv->nvals != brv->nvals ) return ;
659
660 nvec = arv->nvec ; nvals = arv->nvals ;
661 for( iv=0 ; iv < nvec ; iv++ ){
662 av = VECTIM_PTR(arv,iv) ; bv = VECTIM_PTR(brv,iv) ;
663 for( sum=0.0f,ii=0 ; ii < nvals ; ii++ ) sum += av[ii]*bv[ii] ;
664 dp[iv] = sum ;
665 }
666
667 thd_floatscan(nvec,dp) ; /* 16 May 2013 */
668
669 return ;
670 }
671
672 /*---------------------------------------------------------------------*/
673 /* 01 Mar 2010: Rank correlation. */
674
THD_vectim_spearman(MRI_vectim * mrv,float * vec,float * dp)675 void THD_vectim_spearman( MRI_vectim *mrv , float *vec , float *dp )
676 {
677 float *av , *bv , sav ;
678 int nvec, nvals, iv ;
679
680 if( mrv == NULL || vec == NULL || dp == NULL ) return ;
681
682 nvec = mrv->nvec ; nvals = mrv->nvals ;
683 #pragma omp critical (MALLOC)
684 av = (float *)malloc(sizeof(float)*nvals) ;
685 #pragma omp critical (MALLOC)
686 bv = (float *)malloc(sizeof(float)*nvals) ;
687
688 AAmemcpy( av , vec , sizeof(float)*nvals ) ;
689 sav = spearman_rank_prepare( nvals , av ) ; if( sav <= 0.0f ) sav = 1.e+9f ;
690
691 for( iv=0 ; iv < nvec ; iv++ ){
692 AAmemcpy( bv , VECTIM_PTR(mrv,iv) , sizeof(float)*nvals ) ;
693 dp[iv] = spearman_rank_corr( nvals , bv , sav , av ) ;
694 }
695
696 thd_floatscan(nvec,dp) ; /* 16 May 2013 */
697
698 free(bv) ; free(av) ; return ;
699 }
700
701 /*---------------------------------------------------------------------*/
702 /* 01 Mar 2010: Quadrant correlation. */
703
THD_vectim_quadrant(MRI_vectim * mrv,float * vec,float * dp)704 void THD_vectim_quadrant( MRI_vectim *mrv , float *vec , float *dp )
705 {
706 float *av , *bv , sav ;
707 int nvec, nvals, iv ;
708
709 if( mrv == NULL || vec == NULL || dp == NULL ) return ;
710
711 nvec = mrv->nvec ; nvals = mrv->nvals ;
712 #pragma omp critical (MALLOC)
713 av = (float *)malloc(sizeof(float)*nvals) ;
714 #pragma omp critical (MALLOC)
715 bv = (float *)malloc(sizeof(float)*nvals) ;
716
717 AAmemcpy( av , vec , sizeof(float)*nvals ) ;
718 sav = quadrant_corr_prepare( nvals , av ) ; if( sav <= 0.0f ) sav = 1.e+9f ;
719
720 for( iv=0 ; iv < nvec ; iv++ ){
721 AAmemcpy( bv , VECTIM_PTR(mrv,iv) , sizeof(float)*nvals ) ;
722 dp[iv] = quadrant_corr( nvals , bv , sav , av ) ;
723 }
724
725 thd_floatscan(nvec,dp) ; /* 16 May 2013 */
726
727 free(bv) ; free(av) ; return ;
728 }
729
730 /*---------------------------------------------------------------------*/
731 /* 30 Mar 2011: TicTacToe correlation. */
732
THD_vectim_tictactoe(MRI_vectim * mrv,float * vec,float * dp)733 void THD_vectim_tictactoe( MRI_vectim *mrv , float *vec , float *dp )
734 {
735 float *av , *bv , sav ;
736 int nvec, nvals, iv ;
737
738 if( mrv == NULL || vec == NULL || dp == NULL ) return ;
739
740 nvec = mrv->nvec ; nvals = mrv->nvals ;
741 #pragma omp critical (MALLOC)
742 av = (float *)malloc(sizeof(float)*nvals) ;
743 #pragma omp critical (MALLOC)
744 bv = (float *)malloc(sizeof(float)*nvals) ;
745
746 { float tbot , ttop ;
747 tbot = (float)AFNI_numenv("AFNI_TICTACTOE_BOT") ;
748 ttop = (float)AFNI_numenv("AFNI_TICTACTOE_TOP") ;
749 tictactoe_set_thresh( tbot , ttop ) ;
750 }
751
752 AAmemcpy( av , vec , sizeof(float)*nvals ) ;
753 sav = tictactoe_corr_prepare( nvals , av ) ; if( sav <= 0.0f ) sav = 1.e+9f ;
754
755 for( iv=0 ; iv < nvec ; iv++ ){
756 AAmemcpy( bv , VECTIM_PTR(mrv,iv) , sizeof(float)*nvals ) ;
757 dp[iv] = tictactoe_corr( nvals , bv , sav , av ) ;
758 }
759
760 thd_floatscan(nvec,dp) ; /* 16 May 2013 */
761
762 free(bv) ; free(av) ; return ;
763 }
764
765 /*----------------------------------------------------------------------------*/
766 /* 29 Apr 2010: Kendall Tau-b correlation. */
767
THD_vectim_ktaub(MRI_vectim * mrv,float * vec,float * dp)768 void THD_vectim_ktaub( MRI_vectim *mrv , float *vec , float *dp )
769 {
770 float *av , *aav , *bv , *dv ;
771 int nvec , nvals , iv , jv , *qv ;
772
773 ENTRY("THD_vectim_ktaub") ;
774
775 if( mrv == NULL || vec == NULL || dp == NULL ) EXRETURN ;
776
777 nvec = mrv->nvec ; nvals = mrv->nvals ;
778 #pragma omp critical (MALLOC)
779 av = (float *)malloc(sizeof(float)*nvals) ;
780 #pragma omp critical (MALLOC)
781 aav = (float *)malloc(sizeof(float)*nvals) ;
782 #pragma omp critical (MALLOC)
783 bv = (float *)malloc(sizeof(float)*nvals) ;
784 #pragma omp critical (MALLOC)
785 qv = (int *)malloc(sizeof(int )*nvals) ;
786
787 AAmemcpy( av , vec , sizeof(float)*nvals ) ;
788 for( jv=0 ; jv < nvals ; jv++ ) qv[jv] = jv ;
789 STATUS("qsort") ;
790 qsort_floatint( nvals , av , qv ) ;
791
792 STATUS("loop") ;
793 for( iv=0 ; iv < nvec ; iv++ ){
794 dv = VECTIM_PTR(mrv,iv) ;
795 for( jv=0 ; jv < nvals ; jv++ ) bv[jv] = dv[qv[jv]] ;
796 AAmemcpy( aav , av , sizeof(float)*nvals) ;
797 dp[iv] = kendallNlogN( aav , bv , nvals ) ;
798 }
799
800 thd_floatscan(nvec,dp) ; /* 16 May 2013 */
801
802 free(qv) ; free(bv) ; free(aav) ; free(av) ; EXRETURN ;
803 }
804
805 /*---------------------------------------------------------------------*/
806 /* 12 May 2012: Quantile correlation. */
807
THD_vectim_quantile(MRI_vectim * mrv,float * vec,float * dp)808 void THD_vectim_quantile( MRI_vectim *mrv , float *vec , float *dp )
809 {
810 float *av , *bv , sav ;
811 int nvec, nvals, iv ;
812
813 if( mrv == NULL || vec == NULL || dp == NULL ) return ;
814
815 nvec = mrv->nvec ; nvals = mrv->nvals ;
816 #pragma omp critical (MALLOC)
817 av = (float *)malloc(sizeof(float)*nvals) ;
818 #pragma omp critical (MALLOC)
819 bv = (float *)malloc(sizeof(float)*nvals) ;
820
821 AAmemcpy( av , vec , sizeof(float)*nvals ) ;
822 sav = quantile_prepare( nvals , av ) ; if( sav <= 0.0f ) sav = 1.e+9f ;
823
824 for( iv=0 ; iv < nvec ; iv++ ){
825 AAmemcpy( bv , VECTIM_PTR(mrv,iv) , sizeof(float)*nvals ) ;
826 dp[iv] = quantile_corr( nvals , bv , sav , av ) ;
827 }
828
829 thd_floatscan(nvec,dp) ; /* 16 May 2013 */
830
831 free(bv) ; free(av) ; return ;
832 }
833
834 /*---------------------------------------------------------------------*/
835 /* 04 May 2012: Distances. */
836 /* Special parameters:
837 abs: 0 --> Euclidian distance
838 1 --> City Block distance
839 xform: String flags for transforming distance.
840 If string contains "n_scale", scale distance by number
841 of values (dimensions) at each voxel
842 If string contains "inv", return the inverse of the distance
843 */
THD_vectim_distance(MRI_vectim * mrv,float * vec,float * dp,int abs,char * xform)844 void THD_vectim_distance( MRI_vectim *mrv , float *vec ,
845 float *dp, int abs, char *xform)
846 {
847
848 if( mrv == NULL || vec == NULL || dp == NULL ) return ;
849
850 AFNI_OMP_START ;
851 #pragma omp parallel if( mrv->nvec > 1 && mrv->nvec * mrv->nvals > 999999 )
852 { int nvec=mrv->nvec, nvals=mrv->nvals, nv1=nvals-1, iv, ii ;
853 float sum, *fv, a1, a2;
854 #pragma omp for
855 for( iv=0 ; iv < nvec ; iv++ ){
856 fv = VECTIM_PTR(mrv,iv) ;
857 for( sum=0.0f,ii=0 ; ii < nv1 ; ii+=2 ) {
858 a1 = fv[ii]-vec[ii]; a2 = fv[ii+1]-vec[ii+1];
859 if (!abs) sum += (a1*a1+a2*a2);
860 else sum += (ABS(a1)+ABS(a2));
861 }
862 if( ii == nv1 ) {
863 a1 = fv[ii]-vec[ii];
864 if (!abs) sum += a1*a1;
865 else sum += ABS(a1);
866 }
867 if (!abs) dp[iv] = sqrt(sum) ;
868 else dp[iv] = sum;
869 }
870 } /* end OpenMP */
871 AFNI_OMP_END ;
872
873 if (xform) {
874 float a1 = 1.0; int iv, nvec=mrv->nvec;
875 if (strstr(xform,"n_scale")) { a1 = (float)mrv->nvals; }
876 if (strstr(xform,"inv")) {
877 for( iv=0 ; iv < nvec ; iv++ ) {
878 if (dp[iv] != 0.0) {
879 dp[iv] = a1/dp[iv];
880 }
881 }
882 } else if (a1 != 1.0) {
883 for( iv=0 ; iv < nvec ; iv++ ) {
884 if (dp[iv] != 0.0) {
885 dp[iv] = dp[iv]/a1;
886 }
887 }
888 }
889 }
890
891 thd_floatscan(mrv->nvec,dp) ; /* 16 May 2013 */
892
893 return ;
894 }
895
896 /*----------------------------------------------------------------------------*/
897
THD_vectim_pearson_section(MRI_vectim * mrv,float * vec,float * dp,int ibot,int itop)898 void THD_vectim_pearson_section( MRI_vectim *mrv, float *vec, /* 07 Oct 2014 */
899 float *dp, int ibot, int itop )
900 {
901 if( mrv == NULL || vec == NULL || dp == NULL ) return ;
902 if( ibot < 0 ) ibot = 0 ;
903 if( itop >= mrv->nvals ) itop = mrv->nvals-1 ;
904
905 AFNI_OMP_START ;
906 #pragma omp parallel if( mrv->nvec > 1 && mrv->nvec * mrv->nvals > 999999 )
907 { int nvec=mrv->nvec, nvals=itop-ibot+1, iv ; float sum , *fv ;
908 #pragma omp for
909 for( iv=0 ; iv < nvec ; iv++ ){
910 fv = VECTIM_PTR(mrv,iv) ;
911 dp[iv] = THD_pearson_corr( nvals , vec+ibot , fv+ibot ) ;
912 }
913 } /* end OpenMP */
914 AFNI_OMP_END ;
915
916 thd_floatscan(mrv->nvec,dp) ; /* 16 May 2013 */
917
918 return ;
919 }
920
921
922
923 #define USE_VSTEP
924 #ifdef USE_VSTEP
925 /*---------------------------------------------------------------------------*/
926 /*! For voxel loop progress report. */
927
928 static int vstep = 0 ;
929
vstep_print(void)930 static void vstep_print(void)
931 {
932 static char xx[10] = "0123456789" ; static int vn=0 ;
933 fprintf(stderr , "%c" , xx[vn%10] ) ;
934 if( vn%10 == 9) fprintf(stderr,".") ;
935 vn++ ;
936 }
937 #endif
938
939 /*----------------------------------------------------------------------------*/
940
THD_vectim_pearsonBC(MRI_vectim * mrv,float srad,int sijk,int pv,float * par)941 void THD_vectim_pearsonBC( MRI_vectim *mrv, float srad, int sijk, int pv, float *par )
942 {
943 MCW_cluster *smask ;
944 int sqq,qq,pp,pijk,nsar,nyar,nlen,nx,ny,nz,nxy,ii,ik,ij,qi,qj,qk,qjk,mm,nmask ;
945 float **sar, **yar ;
946
947 ENTRY("THD_vectim_pearsonBC") ;
948
949 if( mrv == NULL || par == NULL ) EXRETURN ;
950 sqq = THD_vectim_ifind( sijk , mrv ) ; if( sqq < 0 ) EXRETURN ;
951 if( srad >= 0.0f ){
952 srad = MAX(srad,mrv->dx); srad = MAX(srad,mrv->dy); srad = MAX(srad,mrv->dz);
953 smask = MCW_spheremask(mrv->dx,mrv->dy,mrv->dz,1.001f*srad) ;
954 } else {
955 srad = MAX(-srad,1.01f) ;
956 smask = MCW_spheremask(1.0f,1.0f,1.0f,srad) ;
957 }
958 nmask = smask->num_pt ;
959 nlen = mrv->nvals ;
960
961 nx = mrv->nx ; ny = mrv->ny ; nz = mrv->nz ; nxy = nx*ny ;
962 ii = sijk % nx ; ik = sijk / nxy ; ij = (sijk-ik*nxy) / nx ;
963
964 #pragma omp critical (MALLOC)
965 sar = (float **)malloc(sizeof(float *)*nmask) ;
966 #pragma omp critical (MALLOC)
967 yar = (float **)malloc(sizeof(float *)*nmask) ;
968
969 for( nsar=mm=0 ; mm < nmask ; mm++ ){
970 qi = ii + smask->i[mm] ; if( qi < 0 || qi >= nx ) continue ;
971 qj = ij + smask->j[mm] ; if( qj < 0 || qj >= ny ) continue ;
972 qk = ik + smask->k[mm] ; if( qk < 0 || qk >= nz ) continue ;
973 qjk = qi + qj*nx + qk*nxy ;
974 qq = THD_vectim_ifind( qjk , mrv ) ;
975 if( qq >= 0 ) sar[nsar++] = VECTIM_PTR(mrv,qq) ;
976 }
977
978 #ifdef USE_VSTEP
979 vstep = (mrv->nvec > 999) ? mrv->nvec/50 : 0 ;
980 if( vstep ) fprintf(stderr," + Voxel loop [nmask=%d]: ",nmask) ;
981 #endif
982
983 for( pp=0 ; pp < mrv->nvec ; pp++ ){
984 if( pp == sqq ){ par[pp] = 1.0f ; continue ; }
985 #ifdef USE_VSTEP
986 if( vstep && pp%vstep==vstep-1 ) vstep_print() ;
987 #endif
988 pijk = mrv->ivec[pp] ;
989 ii = pijk % nx ; ik = pijk / nxy ; ij = (pijk-ik*nxy) / nx ;
990 for( nyar=mm=0 ; mm < nmask ; mm++ ){
991 qi = ii + smask->i[mm] ; if( qi < 0 || qi >= nx ) continue ;
992 qj = ij + smask->j[mm] ; if( qj < 0 || qj >= ny ) continue ;
993 qk = ik + smask->k[mm] ; if( qk < 0 || qk >= nz ) continue ;
994 qjk = qi + qj*nx + qk*nxy ;
995 qq = THD_vectim_ifind( qjk , mrv ) ;
996 if( qq >= 0 ) yar[nyar++] = VECTIM_PTR(mrv,qq) ;
997 }
998 par[pp] = THD_bootstrap_vectcorr( nlen , 50 , pv , 1 ,
999 nsar,sar , nyar,yar ) ;
1000 }
1001 #ifdef USE_VSTEP
1002 fprintf(stderr,"\n") ;
1003 #endif
1004
1005 EXRETURN ;
1006 }
1007
1008 /*----------------------------------------------------------------------------*/
1009
THD_vectim_subset_average(MRI_vectim * mrv,int nind,int * ind,float * ar)1010 int THD_vectim_subset_average( MRI_vectim *mrv, int nind, int *ind, float *ar )
1011 {
1012 int nvals , jj,kk,nkk=0 ; register int ii ; float *fv ;
1013
1014 ENTRY("THD_vectim_subset_average") ;
1015
1016 if( mrv == NULL || nind <= 0 || ind == NULL || ar == NULL ) RETURN(0) ;
1017
1018 nvals = mrv->nvals ;
1019
1020 for( ii=0 ; ii < nvals ; ii++ ) ar[ii] = 0.0f ;
1021
1022 for( jj=0 ; jj < nind ; jj++ ){
1023 kk = THD_vectim_ifind( ind[jj] , mrv ) ; if( kk < 0 ) continue ;
1024 fv = VECTIM_PTR(mrv,kk) ;
1025 for( ii=0 ; ii < nvals ; ii++ ) ar[ii] += fv[ii] ;
1026 nkk++ ;
1027 }
1028 if( nkk > 1 ){
1029 register float fac = 1.0f/nkk ;
1030 for( ii=0 ; ii < nvals ; ii++ ) ar[ii] *= fac ;
1031 }
1032
1033 RETURN(nkk) ;
1034 }
1035
1036 /*-----------------------------------------------------------*/
1037 /*! Determine size of a MRI_vectim struct from this dataset. */
1038
THD_vectim_size(THD_3dim_dataset * dset,byte * mask)1039 int64_t THD_vectim_size( THD_3dim_dataset *dset , byte *mask )
1040 {
1041 int nvals , nvox , nmask ;
1042 int64_t sz ;
1043
1044 ENTRY("THD_vectim_size") ;
1045
1046 if( !ISVALID_DSET(dset) ) RETURN(0) ;
1047
1048 nvals = DSET_NVALS(dset) ;
1049 nvox = DSET_NVOX(dset) ;
1050 if( mask != NULL ) nmask = THD_countmask( nvox , mask ) ;
1051 else nmask = DSET_NVOX(dset) ;
1052
1053 sz = ((int64_t)nmask) * ( ((int64_t)nvals) * sizeof(float) + sizeof(int) ) ;
1054 RETURN(sz) ;
1055 }
1056
1057 /*-----------------------------------------------------------------------*/
1058 /*! Binary search in a sorted integer array;
1059 return index of tt in ar if found, return -1 if not found.
1060 *//*---------------------------------------------------------------------*/
1061
bsearch_int(int tt,int nar,int * ar)1062 int bsearch_int( int tt , int nar , int *ar )
1063 {
1064 int targ , ii , jj , kk , nn ;
1065
1066 if( nar <= 0 || ar == NULL ) return -1 ; /* bad inputs */
1067
1068 targ = tt ; ii = 0 ; jj = nar-1 ; /* setup */
1069
1070 if( targ < ar[0] ) return -1 ; /* not found */
1071 else if( targ == ar[0] ) return 0 ; /* at start! */
1072
1073 if( targ > ar[jj] ) return -1 ; /* not found */
1074 else if( targ == ar[jj] ) return jj ; /* at end! */
1075
1076 /* at the start of this loop, we've already checked
1077 indexes ii and jj, so check the middle of them (kk),
1078 and if that doesn't work, make the middle the
1079 new ii or the new jj -- so again we will have
1080 checked both ii and jj when the loop iterates back. */
1081
1082 while( jj-ii > 1 ){
1083 kk = (ii+jj) / 2 ; /* midpoint */
1084 nn = ar[kk] - targ ; /* sign of difference */
1085 if( nn == 0 ) return kk ; /* found it! */
1086 if( nn < 0 ) ii = kk ; /* must be above kk */
1087 else jj = kk ; /* must be below kk */
1088 }
1089
1090 return -1 ;
1091 }
1092
1093 /*-----------------------------------------------------------------------*/
1094 /*! Find a voxel index in a MRI_vectim struct. If not found, return -1. */
1095
THD_vectim_ifind(int iv,MRI_vectim * mrv)1096 int THD_vectim_ifind( int iv , MRI_vectim *mrv )
1097 {
1098 if( mrv == NULL ) return -1 ; /* stoopid user */
1099 return bsearch_int( iv , mrv->nvec , mrv->ivec ) ;
1100 }
1101
1102 /*----------------------------------------------------------------------*/
1103 /*! Note that the dataset must already have bricks set up! Or else! */
1104
THD_vectim_to_dset(MRI_vectim * mrv,THD_3dim_dataset * dset)1105 void THD_vectim_to_dset( MRI_vectim *mrv , THD_3dim_dataset *dset )
1106 {
1107 int nvals , nvec , kk , ign ;
1108
1109 ENTRY("THD_vectim_to_dset") ;
1110
1111 if( mrv == NULL || !ISVALID_DSET(dset) ) EXRETURN ;
1112 if( mrv->nvals + mrv->ignore != DSET_NVALS(dset) ) EXRETURN ;
1113
1114 nvec = mrv->nvec ;
1115 nvals = mrv->nvals ;
1116 ign = mrv->ignore ;
1117
1118 if( ign == 0 ){
1119 for( kk=0 ; kk < nvec ; kk++ )
1120 THD_insert_series( mrv->ivec[kk] , dset ,
1121 nvals , MRI_float , VECTIM_PTR(mrv,kk) , 0 ) ;
1122 } else {
1123 float *var ;
1124 #pragma omp critical (MALLOC)
1125 var = (float *)malloc(sizeof(float)*(nvals+ign)) ;
1126 for( kk=0 ; kk < nvec ; kk++ ){
1127 (void)THD_extract_array( mrv->ivec[kk] , dset , 0 , var ) ;
1128 AAmemcpy( var+ign , VECTIM_PTR(mrv,kk) , sizeof(float)*nvals ) ;
1129 THD_insert_series( mrv->ivec[kk] , dset ,
1130 nvals , MRI_float , var , 0 ) ;
1131 }
1132 }
1133
1134 EXRETURN ;
1135 }
1136
1137 /*----------------------------------------------------------------------*/
1138 /* The jj-th time point in the vectim goes to the tlist[jj]-th time
1139 point in the output dataset [06 Aug 2013].
1140 *//*--------------------------------------------------------------------*/
1141
THD_vectim_to_dset_indexed(MRI_vectim * mrv,THD_3dim_dataset * dset,int * tlist)1142 void THD_vectim_to_dset_indexed( MRI_vectim *mrv ,
1143 THD_3dim_dataset *dset , int *tlist )
1144 {
1145 int nvals , nvec , jj,kk , tmax=0 ;
1146 float *tar , *var ;
1147
1148 ENTRY("THD_vectim_to_dset_indexed") ;
1149
1150 if( mrv == NULL || !ISVALID_DSET(dset) || tlist == NULL ) EXRETURN ;
1151
1152 nvec = mrv->nvec ;
1153 nvals = mrv->nvals ;
1154
1155 for( kk=0 ; kk < nvals ; kk++ ){
1156 if( tlist[kk] < 0 ) EXRETURN ;
1157 if( tlist[kk] > tmax ) tmax = tlist[kk] ;
1158 }
1159 tmax++ ; if( DSET_NVALS(dset) < tmax ) EXRETURN ;
1160
1161 tar = (float *)malloc(sizeof(float)*tmax) ;
1162
1163 for( kk=0 ; kk < nvec ; kk++ ){
1164 var = VECTIM_PTR(mrv,kk) ;
1165 for( jj=0 ; jj < tmax ; jj++ ) tar[jj] = 0.0f ;
1166 for( jj=0 ; jj < nvals ; jj++ ) tar[tlist[jj]] = var[jj] ;
1167 THD_insert_series( mrv->ivec[kk] , dset ,
1168 tmax , MRI_float , tar , 0 ) ;
1169 }
1170
1171 free(tar) ; EXRETURN ;
1172 }
1173
1174 /*----------------------------------------------------------------------*/
1175 /* The ilist[jj]-th point in the vectim goes into the jj-th index
1176 in the output dataset [06 Feb 2014].
1177 *//*--------------------------------------------------------------------*/
1178
THD_vectim_indexed_to_dset(MRI_vectim * mrv,int nlist,int * ilist,THD_3dim_dataset * dset)1179 void THD_vectim_indexed_to_dset( MRI_vectim *mrv, int nlist, int *ilist,
1180 THD_3dim_dataset *dset )
1181 {
1182 int nvals , nvec , jj,kk ;
1183 float *tar , *var ;
1184
1185 ENTRY("THD_vectim_indexed_to_dset") ;
1186
1187 if( mrv == NULL || !ISVALID_DSET(dset) ||
1188 nlist <= 0 || ilist == NULL || nlist > DSET_NVALS(dset) ){
1189 ERROR_message("THD_vectim_indexed_to_dset: illegal inputs (nlist=%d)",nlist) ;
1190 EXRETURN ;
1191 }
1192
1193 nvec = mrv->nvec ;
1194 nvals = mrv->nvals ;
1195
1196 for( kk=0 ; kk < nlist ; kk++ ){
1197 if( ilist[kk] < 0 || ilist[kk] >= nvals ){
1198 ERROR_message("THD_vectim_indexed_to_dset: illegal ilist[%d]=%d",kk,ilist[kk]) ;
1199 EXRETURN ;
1200 }
1201 }
1202
1203 tar = (float *)malloc(sizeof(float)*nlist) ;
1204
1205 for( kk=0 ; kk < nvec ; kk++ ){
1206 var = VECTIM_PTR(mrv,kk) ;
1207 for( jj=0 ; jj < nlist ; jj++ ) tar[jj] = var[ilist[jj]] ;
1208 THD_insert_series( mrv->ivec[kk] , dset ,
1209 nlist , MRI_float , tar , 0 ) ;
1210 }
1211
1212 free(tar) ; EXRETURN ;
1213 }
1214
1215 /*---------------------------------------------------------------------------*/
1216
THD_tcat_vectims(int nvim,MRI_vectim ** vim)1217 MRI_vectim * THD_tcat_vectims( int nvim , MRI_vectim **vim )
1218 {
1219 MRI_vectim *vout ;
1220 int iv , nvec , nvsum , vv , nvals ; size_t nvv ;
1221 float *vout_ptr , *vin_ptr ;
1222
1223 ENTRY("THD_tcat_vectims") ;
1224
1225 if( nvim <= 0 || vim == NULL ) RETURN(NULL) ;
1226
1227 if( nvim == 1 ){
1228 vout = THD_vectim_copy( vim[0] ) ; RETURN(vout) ;
1229 }
1230
1231 nvec = vim[0]->nvec ;
1232 nvsum = vim[0]->nvals ;
1233 for( iv=1 ; iv < nvim ; iv++ ){
1234 if( vim[iv]->nvec != nvec ) RETURN(NULL) ;
1235 nvsum += vim[iv]->nvals ;
1236 }
1237
1238 MAKE_VECTIM(vout,nvec,nvsum) ;
1239 vout->ignore = 0 ;
1240 vout->nx = vim[0]->nx ; vout->dx = vim[0]->dx ;
1241 vout->ny = vim[0]->ny ; vout->dy = vim[0]->dy ;
1242 vout->nz = vim[0]->nz ; vout->dz = vim[0]->dz ; vout->dt = vim[0]->dt ;
1243 AAmemcpy( vout->ivec , vim[0]->ivec , sizeof(int)*vim[0]->nvec ) ;
1244
1245 for( nvv=iv=0 ; iv < nvim ; iv++,nvv+=nvals ){
1246 nvals = vim[iv]->nvals ;
1247 for( vv=0 ; vv < nvec ; vv++ ){
1248 vout_ptr = VECTIM_PTR(vout,vv) + nvv ;
1249 vin_ptr = VECTIM_PTR(vim[iv],vv) ;
1250 AAmemcpy( vout_ptr , vin_ptr , sizeof(float)*nvals ) ;
1251 }
1252 }
1253
1254 RETURN(vout) ;
1255 }
1256
1257 /*---------------------------------------------------------------------------*/
1258
THD_dset_list_to_vectim(int nds,THD_3dim_dataset ** ds,byte * mask)1259 MRI_vectim * THD_dset_list_to_vectim( int nds, THD_3dim_dataset **ds, byte *mask )
1260 {
1261 MRI_vectim *vout , **vim ;
1262 int kk , jj ;
1263
1264 ENTRY("THD_dset_list_to_vectim") ;
1265
1266 if( nds < 1 || ds == NULL ) RETURN(NULL) ;
1267
1268 if( nds == 1 ) RETURN( THD_dset_to_vectim( ds[0] , mask , 0 ) ) ;
1269
1270 for( kk=0 ; kk < nds ; kk++ )
1271 if( !ISVALID_DSET(ds[kk]) ) RETURN(NULL) ;
1272
1273 #pragma omp critical (MALLOC)
1274 vim = (MRI_vectim **)malloc(sizeof(MRI_vectim *)*nds) ;
1275 for( kk=0 ; kk < nds ; kk++ ){
1276 vim[kk] = THD_dset_to_vectim( ds[kk] , mask , 0 ) ;
1277 /** DSET_unload( ds[kk] ) ; **/
1278 if( vim[kk] == NULL ){
1279 for( jj=0 ; jj < kk ; jj++ ) VECTIM_destroy(vim[jj]) ;
1280 free(vim) ; RETURN(NULL) ;
1281 }
1282 }
1283
1284 vout = THD_tcat_vectims( nds , vim ) ;
1285 for( jj=0 ; jj < nds ; jj++ ) VECTIM_destroy(vim[jj]) ;
1286 free(vim) ; RETURN(vout) ;
1287 }
1288
1289 /*---------------------------------------------------------------------------*/
1290
THD_dset_list_censored_to_vectim(int nds,THD_3dim_dataset ** ds,byte * mask,int nkeep,int * keep)1291 MRI_vectim * THD_dset_list_censored_to_vectim( int nds, THD_3dim_dataset **ds,
1292 byte *mask , int nkeep , int *keep )
1293 {
1294 MRI_vectim *vout , **vim ;
1295 int kk , jj ;
1296
1297 ENTRY("THD_dset_list_censored_to_vectim") ;
1298
1299 if( nds < 1 || ds == NULL ) RETURN(NULL) ;
1300
1301 if( nds == 1 ) /* trivial case */
1302 RETURN( THD_dset_censored_to_vectim( ds[0],mask,nkeep,keep ) ) ;
1303
1304 for( kk=0 ; kk < nds ; kk++ ){
1305 if( !ISVALID_DSET(ds[kk]) ) RETURN(NULL) ;
1306 if( DSET_NVALS(ds[kk]) != DSET_NVALS(ds[0]) ) RETURN(NULL) ;
1307 }
1308
1309 #pragma omp critical (MALLOC)
1310 vim = (MRI_vectim **)malloc(sizeof(MRI_vectim *)*nds) ;
1311 for( kk=0 ; kk < nds ; kk++ ){
1312 vim[kk] = THD_dset_censored_to_vectim( ds[kk] , mask , nkeep,keep ) ;
1313 /** DSET_unload( ds[kk] ) ; **/
1314 if( vim[kk] == NULL ){
1315 for( jj=0 ; jj < kk ; jj++ ) VECTIM_destroy(vim[jj]) ;
1316 free(vim) ; RETURN(NULL) ;
1317 }
1318 }
1319
1320 vout = THD_tcat_vectims( nds , vim ) ;
1321 for( jj=0 ; jj < nds ; jj++ ) VECTIM_destroy(vim[jj]) ;
1322 free(vim) ; RETURN(vout) ;
1323 }
1324
1325 /*---------------------------------------------------------------------------*/
1326
THD_check_vectim(MRI_vectim * mv,char * fname)1327 void THD_check_vectim( MRI_vectim *mv , char *fname )
1328 {
1329 int nvec , nvals ;
1330 float *vpt , vz ;
1331 int nbad , ii,jj ;
1332
1333 if( fname == NULL ) fname = "vectim check" ;
1334
1335 if( mv == NULL ){
1336 WARNING_message("%s :: bad input vector",fname); return;
1337 }
1338
1339 nvec = mv->nvec ;
1340 nvals = mv->nvals ;
1341
1342 /* scan each time series for constancy */
1343
1344 for( nbad=jj=0 ; jj < nvec ; jj++ ){
1345 vpt = VECTIM_PTR(mv,jj) ; vz = vpt[0] ;
1346 for( ii=1 ; ii < nvals && vpt[ii] == vz ; ii++ ) ; /*nada*/
1347 if( ii == nvals ) nbad++ ;
1348 }
1349 if( nbad > 0 && nvals > 1 )
1350 WARNING_message("%s :: %d vector%s constant",
1351 fname , nbad , (nbad==1) ? " is" : "s are" ) ;
1352
1353 /* scan each time point for constancy */
1354
1355 for( nbad=ii=0 ; ii < nvals ; ii++ ){
1356 vpt = VECTIM_PTR(mv,0) ; vz = vpt[ii] ;
1357 for( jj=1 ; jj < nvec ; jj++ ){
1358 vpt = VECTIM_PTR(mv,jj) ; if( vpt[ii] != vz ) break ;
1359 }
1360 if( jj == nvec ) nbad++ ;
1361 }
1362 if( nbad > 0 && nvec > 1 )
1363 WARNING_message("%s :: %d volume%s constant",
1364 fname , nbad , (nbad==1) ? " is" : "s are" ) ;
1365 return ;
1366 }
1367
1368 /*---------------------------------------------------------------------*/
1369
THD_xyzcat_vectims(int nvim,MRI_vectim ** vim)1370 MRI_vectim * THD_xyzcat_vectims( int nvim , MRI_vectim **vim )
1371 {
1372 MRI_vectim *vout ;
1373 int nx,ny,nz,nv , nvectot , iv,ii,jj ;
1374 float *vout_ptr , *vin_ptr ;
1375
1376 ENTRY("THD_xyzcat_vectims") ;
1377
1378 if( nvim <= 0 || vim == NULL ) RETURN(NULL) ;
1379
1380 if( nvim == 1 ){
1381 vout = THD_vectim_copy( vim[0] ) ; RETURN(vout) ;
1382 }
1383
1384 nx = vim[0]->nx ;
1385 ny = vim[0]->ny ;
1386 nz = vim[0]->nz ;
1387 nv = vim[0]->nvals ;
1388 nvectot = vim[0]->nvec ;
1389
1390 for( iv=1 ; iv < nvim ; iv++ ){
1391 if( vim[iv]->nx != nx ||
1392 vim[iv]->ny != ny ||
1393 vim[iv]->nz != nz ||
1394 vim[iv]->nvals != nv ) RETURN(NULL) ;
1395
1396 nvectot += vim[iv]->nvec ;
1397 }
1398
1399 MAKE_VECTIM( vout , nvectot , nv ) ;
1400 vout->nx = nx ; vout->dx = vim[0]->dx ;
1401 vout->ny = ny ; vout->dy = vim[0]->dy ;
1402 vout->nz = nz ; vout->dz = vim[0]->dz ; vout->dt = vim[0]->dt ;
1403
1404 for( jj=iv=0 ; iv < nvim ; iv++ ){
1405 for( ii=0 ; ii < vim[iv]->nvec ; ii++,jj++ ){
1406 vout->ivec[jj] = vim[iv]->ivec[ii] ;
1407 memcpy( VECTIM_PTR(vout,jj) , VECTIM_PTR(vim[iv],ii) , sizeof(float)*nv ) ;
1408 }
1409 }
1410
1411 RETURN(vout) ;
1412 }
1413
1414 /*---------------------------------------------------------------------*/
1415 /* To discard the output from this function, use mri_clear_and_free()
1416 NOT mri_free() -- the MRI_IMAGE data is just a pointer into the
1417 vectim, and trying to free it would be disastrous! [19 Nov 2021]
1418 To put ALL vectors into the temp image, use istart=numi=0.
1419 *//*-------------------------------------------------------------------*/
1420
THD_temp_subim_from_vectim(MRI_vectim * vim,int istart,int numi)1421 MRI_IMAGE * THD_temp_subim_from_vectim( MRI_vectim *vim ,
1422 int istart , int numi )
1423 {
1424 MRI_IMAGE *imout ;
1425
1426 ENTRY("THD_temp_subim_from_vectim") ;
1427
1428 if( vim == NULL ) RETURN(NULL) ;
1429 if( istart < 0 || istart >= vim->nvec ) RETURN(NULL) ;
1430 if( numi <= 0 ){
1431 numi = vim->nvec - istart ;
1432 } else {
1433 if( istart+numi-1 >= vim->nvec ) RETURN(NULL) ;
1434 }
1435
1436 imout = mri_new_vol_empty( vim->nvals , numi , 1 , MRI_float ) ;
1437 mri_fix_data_pointer( VECTIM_PTR(vim,istart) , imout ) ;
1438
1439 RETURN(imout) ;
1440 }
1441