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 Routine to extract a time-series (fixed index, variable ival)
12 from a dataset.
13 ixyz = spatial index of desired voxel
14 = ix + jy * n1 + kz * n1*n2
15 raw = 0 if you always want floats
16 = 1 if you want the truly stored data type
17
18 05 Nov 2001: Split into 2 functions -
19 THD_extract_series() produces a new image each time
20 THD_extract_array() copies data into a user-supplied array
21 -----------------------------------------------------------------*/
22
THD_extract_series(int ind,THD_3dim_dataset * dset,int raw)23 MRI_IMAGE * THD_extract_series( int ind , THD_3dim_dataset *dset , int raw )
24 {
25 int nv , typ , ii ;
26 MRI_IMAGE *im ;
27 void *imar ;
28
29 ENTRY("THD_extract_series") ;
30
31 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
32
33 nv = dset->dblk->nvals ;
34 if( raw ) typ = DSET_BRICK_TYPE(dset,0) ; /* type of output array */
35 else typ = MRI_float ;
36
37 im = mri_new( nv , 1 , typ ) ; /* output image */
38 imar = mri_data_pointer(im) ;
39
40 ii = THD_extract_array( ind , dset , raw , imar ) ; /* get data */
41
42 if( ii != 0 ){ mri_free(im) ; RETURN(NULL) ; } /* bad */
43
44 if( dset->taxis != NULL ){ /* 21 Oct 1996 */
45 float zz , tt ;
46 int kz = ind / ( dset->daxes->nxx * dset->daxes->nyy ) ;
47
48 zz = dset->daxes->zzorg + kz * dset->daxes->zzdel ;
49 tt = THD_timeof( 0 , zz , dset->taxis ) ;
50
51 im->xo = tt ; im->dx = dset->taxis->ttdel ; /* origin and delta */
52
53 if( dset->taxis->units_type == UNITS_MSEC_TYPE ){ /* convert to sec */
54 im->xo *= 0.001 ; im->dx *= 0.001 ;
55 }
56 } else {
57 im->xo = 0.0 ; im->dx = 1.0 ; /* 08 Nov 1996 */
58 }
59
60 MRI_floatscan(im) ; /* 10 Jun 2021 */
61 RETURN(im) ;
62 }
63
64 /*---------------------------------------------------------------------------
65 Return value is 0 for all is good, -1 for all is bad.
66 If raw == 0, uar is of type float.
67 If raw != 0, uar is of the same type as the dataset brick
68 (the user needs to know the type ahead of time if raw != 0).
69 Data goes into a user-supplied array.
70 -----------------------------------------------------------------------------*/
71
THD_extract_array(int ind,THD_3dim_dataset * dset,int raw,void * uar)72 int THD_extract_array( int ind, THD_3dim_dataset *dset, int raw, void *uar )
73 {
74 MRI_TYPE typ ;
75 int nv , ival , nb , nb1 ;
76 char *iar ; /* brick in the input */
77 float *far=NULL ; /* non-raw output */
78 void *tar=NULL ;
79
80 /** ENTRY("THD_extract_array") ; **/
81
82 if( ind < 0 || uar == NULL ||
83 !ISVALID_DSET(dset) || ind >= DSET_NVOX(dset) ) return(-1) ;
84
85 nv = dset->dblk->nvals ;
86 iar = DSET_ARRAY(dset,0) ;
87 if( iar == NULL ){ /* load data from disk? */
88 DSET_load(dset) ;
89 iar = DSET_ARRAY(dset,0); if( iar == NULL ) return(-1) ;
90 }
91 typ = DSET_BRICK_TYPE(dset,0) ; /* raw data type */
92
93 /* will extract nb bytes of raw data into array tar */
94
95 nb1 = mri_datum_size(typ); nb = nb1 * (nv+1); nb1 = nb1 * nv;
96 tar = (void *)calloc(1,nb) ;
97 NULL_CHECK(tar) ;
98
99 if( !raw ) far = (float *)uar ; /* non-raw output */
100
101 switch( typ ){
102
103 default: /* don't know what to do --> return nada */
104 free(tar); return(-1);
105 break ;
106
107 case MRI_byte:{
108 byte *ar = (byte *)tar , *bar ;
109 for( ival=0 ; ival < nv ; ival++ ){
110 bar = (byte *) DSET_ARRAY(dset,ival) ;
111 if( bar != NULL ) ar[ival] = bar[ind] ;
112 }
113 if( !raw ){
114 for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
115 }
116 }
117 break ;
118
119 case MRI_short:{
120 short *ar = (short *)tar , *bar ;
121 for( ival=0 ; ival < nv ; ival++ ){
122 bar = (short *) DSET_ARRAY(dset,ival) ;
123 if( bar != NULL ) ar[ival] = bar[ind] ;
124 }
125 if( !raw ){
126 for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
127 }
128 }
129 break ;
130
131 case MRI_float:{
132 float *ar = (float *)tar , *bar ;
133 for( ival=0 ; ival < nv ; ival++ ){
134 bar = (float *) DSET_ARRAY(dset,ival) ;
135 if( bar != NULL ) ar[ival] = bar[ind] ;
136 }
137 if( !raw ){
138 for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
139 }
140 }
141 break ;
142
143 #if 0
144 case MRI_int:{
145 int *ar = (int *)tar , *bar ;
146 for( ival=0 ; ival < nv ; ival++ ){
147 bar = (int *) DSET_ARRAY(dset,ival) ;
148 if( bar != NULL ) ar[ival] = bar[ind] ;
149 }
150 if( !raw ){
151 for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
152 }
153 }
154 break ;
155
156 case MRI_double:{
157 double *ar = (double *)tar , *bar ;
158 for( ival=0 ; ival < nv ; ival++ ){
159 bar = (double *) DSET_ARRAY(dset,ival) ;
160 if( bar != NULL ) ar[ival] = bar[ind] ;
161 }
162 if( !raw ){
163 for( ival=0 ; ival < nv ; ival++ ) far[ival] = ar[ival] ;
164 }
165 }
166 break ;
167 #endif
168
169 case MRI_complex:{
170 complex *ar = (complex *)tar , *bar ;
171 for( ival=0 ; ival < nv ; ival++ ){
172 bar = (complex *) DSET_ARRAY(dset,ival) ;
173 if( bar != NULL ) ar[ival] = bar[ind] ;
174 }
175 if( !raw ){
176 for( ival=0 ; ival < nv ; ival++ ) far[ival] = CABS(ar[ival]) ;
177 }
178 }
179 break ;
180
181 }
182
183 if( raw ){ memcpy(uar,tar,nb1); free(tar); return(0); }
184
185 thd_floatscan(nv,far) ; /* 10 Jun 2021 */
186
187 if( THD_need_brick_factor(dset) ){
188 for( ival=0 ; ival < nv ; ival++ )
189 if( DSET_BRICK_FACTOR(dset,ival) > 0.0 )
190 far[ival] *= DSET_BRICK_FACTOR(dset,ival) ;
191 }
192
193 free(tar); return(0);
194 }
195
196 /*---------------------------------------------------------------------------
197 Return value is 0 for all is good, -1 for all is bad.
198 Data goes into a user-supplied array.
199 -----------------------------------------------------------------------------*/
200
THD_extract_float_array(int ind,THD_3dim_dataset * dset,float * far)201 int THD_extract_float_array( int ind, THD_3dim_dataset *dset, float *far )
202 {
203 MRI_TYPE typ ;
204 int nv , ival , nb , nb1 ;
205 char *iar ; /* brick in the input */
206
207 if( ind < 0 || far == NULL ||
208 !ISVALID_DSET(dset) || ind >= DSET_NVOX(dset) ) return(-1) ;
209
210 nv = dset->dblk->nvals ;
211 typ = DSET_BRICK_TYPE(dset,0) ; /* raw data type */
212
213 switch( typ ){
214
215 default: /* don't know what to do --> return nada */
216 return(-1);
217 break ;
218
219 case MRI_byte:{
220 byte *bar ;
221 for( ival=0 ; ival < nv ; ival++ ){
222 bar = (byte *) DSET_ARRAY(dset,ival) ;
223 if( bar != NULL ) far[ival] = bar[ind] ;
224 }
225 }
226 break ;
227
228 case MRI_short:{
229 short *bar ;
230 for( ival=0 ; ival < nv ; ival++ ){
231 bar = (short *) DSET_ARRAY(dset,ival) ;
232 if( bar != NULL ) far[ival] = bar[ind] ;
233 }
234 }
235 break ;
236
237 case MRI_float:{
238 float *bar ;
239 for( ival=0 ; ival < nv ; ival++ ){
240 bar = (float *) DSET_ARRAY(dset,ival) ;
241 if( bar != NULL ) far[ival] = bar[ind] ;
242 }
243 }
244 break ;
245
246 case MRI_complex:{
247 complex *bar ;
248 for( ival=0 ; ival < nv ; ival++ ){
249 bar = (complex *) DSET_ARRAY(dset,ival) ;
250 if( bar != NULL ) far[ival] = CABS(bar[ind]) ;
251 }
252 }
253 break ;
254
255 }
256
257 thd_floatscan(nv,far) ; /* 10 Jun 2021 */
258
259 if( THD_need_brick_factor(dset) ){
260 for( ival=0 ; ival < nv ; ival++ )
261 if( DSET_BRICK_FACTOR(dset,ival) > 0.0 )
262 far[ival] *= DSET_BRICK_FACTOR(dset,ival) ;
263 }
264
265 return(0);
266 }
267
268 /*---------------------------------------------------------------------------*/
269
THD_get_float_value(int ind,int ival,THD_3dim_dataset * dset)270 float THD_get_float_value( int ind , int ival , THD_3dim_dataset *dset )
271 {
272 MRI_TYPE typ ; float val=0.0f ;
273
274 if( ind < 0 || ival < 0 || !ISVALID_DSET(dset) ||
275 ival >= DSET_NVALS(dset) || ind >= DSET_NVOX(dset) ) return val ;
276
277 typ = DSET_BRICK_TYPE(dset,ival) ; /* raw data type */
278
279 switch( typ ){
280
281 default: /* don't know what to do --> return nada */
282 return(0.0f);
283 break ;
284
285 case MRI_byte:{
286 byte *bar ;
287 bar = (byte *) DSET_ARRAY(dset,ival) ;
288 if( bar != NULL ) val = (float)bar[ind] ;
289 }
290 break ;
291
292 case MRI_short:{
293 short *bar ;
294 bar = (short *) DSET_ARRAY(dset,ival) ;
295 if( bar != NULL ) val = (float)bar[ind] ;
296 }
297 break ;
298
299 case MRI_float:{
300 float *bar ;
301 bar = (float *) DSET_ARRAY(dset,ival) ;
302 if( bar != NULL ) val = bar[ind] ;
303 }
304 break ;
305
306 case MRI_complex:{
307 complex *bar ;
308 bar = (complex *) DSET_ARRAY(dset,ival) ;
309 if( bar != NULL ) val = CABS(bar[ind]) ;
310 }
311 break ;
312
313 }
314
315 if( DSET_BRICK_FACTOR(dset,ival) > 0.0f )
316 val *= DSET_BRICK_FACTOR(dset,ival) ;
317
318 thd_floatscan(1,&val) ; /* 10 Jun 2021 */
319 return val ;
320 }
321
322 /*----------------------------------------------------------------------------*/
323
THD_voxel_is_constant(int ind,THD_3dim_dataset * dset)324 int THD_voxel_is_constant( int ind , THD_3dim_dataset *dset )
325 {
326 float *far ; int ii,nvox,nvals ;
327
328 if( !ISVALID_DSET(dset) ) return 1 ;
329 if( ind < 0 || ind >= DSET_NVOX(dset) ) return 1 ;
330
331 nvals = DSET_NVALS(dset) ; if( nvals == 1 ) return 1 ;
332 far = (float *)malloc(sizeof(float)*nvals) ; NULL_CHECK(far) ;
333 ii = THD_extract_array( ind , dset , 0 , far ) ;
334 if( ii < 0 ){ free(far); return 1; }
335 for( ii=1 ; ii < nvals && far[ii]==far[0]; ii++ ) ; /*nada*/
336 free(far) ; return (ii==nvals) ;
337 }
338
339 /*----------------------------------------------------------------------------
340 04 Feb 2000: do a bunch of timeseries at once (for efficiency)
341 27 Feb 2003: rearranged slightly for more efficiency
342 ------------------------------------------------------------------------------*/
343
THD_extract_many_series(int ns,int * ind,THD_3dim_dataset * dset)344 MRI_IMARR * THD_extract_many_series( int ns, int *ind, THD_3dim_dataset *dset )
345 {
346 MRI_IMARR *imar ; /* output */
347 MRI_IMAGE *im ;
348 int nv , ival , kk ;
349 char *iar ; /* brick in the input */
350 float **far ; /* 27 Feb 2003: ptrs to output */
351
352 ENTRY("THD_extract_many_series") ;
353
354 if( ns <= 0 || ind == NULL || dset == NULL ) RETURN( NULL );
355
356 /* try to load dataset */
357
358 nv = dset->dblk->nvals ;
359 iar = DSET_ARRAY(dset,0) ;
360 if( iar == NULL ){ /* if data needs to be loaded from disk */
361 (void) THD_load_datablock( dset->dblk ) ;
362 iar = DSET_ARRAY(dset,0) ;
363 if( iar == NULL ){
364 static int nerr=0 ;
365 if( nerr < 2 ){ ERROR_message("Can't load dataset %s",DSET_HEADNAME(dset)); nerr++; }
366 RETURN( NULL );
367 }
368 }
369
370 /* create output */
371
372 far = (float **) malloc(sizeof(float *)*ns) ; /* 27 Feb 2003 */
373 NULL_CHECK(far) ;
374 INIT_IMARR(imar) ;
375 for( kk=0 ; kk < ns ; kk++ ){
376 im = mri_new( nv , 1 , MRI_float ) ; /* N.B.: now does 0 fill */
377 far[kk] = MRI_FLOAT_PTR(im) ; /* ptr to kk-th output series */
378 ADDTO_IMARR(imar,im) ;
379 }
380
381 /* fill the output */
382
383 switch( DSET_BRICK_TYPE(dset,0) ){
384
385 default: /* don't know what to do --> return nada */
386 DESTROY_IMARR(imar) ; free(far) ;
387 RETURN( NULL );
388
389 case MRI_byte:{
390 byte * bar ;
391 for( ival=0 ; ival < nv ; ival++ ){
392 bar = (byte *) DSET_ARRAY(dset,ival) ;
393 if( bar != NULL ){
394 for( kk=0 ; kk < ns ; kk++ ){
395 far[kk][ival] = (float)bar[ind[kk]] ;
396 }
397 }
398 }
399 }
400 break ;
401
402 case MRI_short:{
403 short * bar ;
404 for( ival=0 ; ival < nv ; ival++ ){
405 bar = (short *) DSET_ARRAY(dset,ival) ;
406 if( bar != NULL ){
407 for( kk=0 ; kk < ns ; kk++ ){
408 far[kk][ival] = (float)bar[ind[kk]] ;
409 }
410 }
411 }
412 }
413 break ;
414
415 case MRI_float:{
416 float * bar ;
417 for( ival=0 ; ival < nv ; ival++ ){
418 bar = (float *) DSET_ARRAY(dset,ival) ;
419 if( bar != NULL ){
420 for( kk=0 ; kk < ns ; kk++ ){
421 far[kk][ival] = bar[ind[kk]] ;
422 }
423 }
424 }
425 }
426 break ;
427
428 #if 0
429 case MRI_int:{
430 int * bar ;
431 for( ival=0 ; ival < nv ; ival++ ){
432 bar = (int *) DSET_ARRAY(dset,ival) ;
433 if( bar != NULL ){
434 for( kk=0 ; kk < ns ; kk++ ){
435 far[kk][ival] = bar[ind[kk]] ;
436 }
437 }
438 }
439 }
440 break ;
441
442 case MRI_double:{
443 double * bar ;
444 for( ival=0 ; ival < nv ; ival++ ){
445 bar = (double *) DSET_ARRAY(dset,ival) ;
446 if( bar != NULL ){
447 for( kk=0 ; kk < ns ; kk++ ){
448 far[kk][ival] = (float)bar[ind[kk]] ;
449 }
450 }
451 }
452 }
453 break ;
454 #endif
455
456 case MRI_complex:{
457 complex * bar ;
458 for( ival=0 ; ival < nv ; ival++ ){
459 bar = (complex *) DSET_ARRAY(dset,ival) ;
460 if( bar != NULL ){
461 for( kk=0 ; kk < ns ; kk++ ){
462 far[kk][ival] = bar[ind[kk]].r ;
463 }
464 }
465 }
466 }
467 break ;
468
469 }
470
471 for( kk=0 ; kk < ns ; kk++ ){ /* 10 Jun 2021 */
472 im = IMARR_SUBIM(imar,kk) ;
473 MRI_floatscan(im) ;
474 }
475
476 /* scale outputs, if needed */
477
478 if( THD_need_brick_factor(dset) ){
479 MRI_IMAGE *qim ;
480 for( kk=0 ; kk < ns ; kk++ ){
481 im = IMARR_SUBIMAGE(imar,kk) ;
482 qim = mri_mult_to_float( dset->dblk->brick_fac , im ) ;
483 mri_free(im) ;
484 IMARR_SUBIMAGE(imar,kk) = qim ;
485 }
486 }
487
488 #if 0 /* 27 Feb 2003 */
489 /* convert to floats, if needed */
490
491 if( IMARR_SUBIMAGE(imar,0)->kind != MRI_float ){
492 MRI_IMAGE * qim ;
493 for( kk=0 ; kk < ns ; kk++ ){
494 im = IMARR_SUBIMAGE(imar,kk) ;
495 qim = mri_to_float( im ) ;
496 mri_free(im) ;
497 IMARR_SUBIMAGE(imar,kk) = qim ;
498 }
499 }
500 #endif
501
502 /* add time axis stuff to output images, if present */
503
504 if( dset->taxis != NULL ){
505 float zz , tt ;
506 int kz ;
507
508 for( kk=0 ; kk < ns ; kk++ ){
509 kz = ind[kk] / ( dset->daxes->nxx * dset->daxes->nyy ) ;
510 zz = dset->daxes->zzorg + kz * dset->daxes->zzdel ;
511 tt = THD_timeof( 0 , zz , dset->taxis ) ;
512 im = IMARR_SUBIMAGE(imar,kk) ;
513 im->xo = tt ; im->dx = dset->taxis->ttdel ; /* origin and delta */
514 if( dset->taxis->units_type == UNITS_MSEC_TYPE ){ /* convert to sec */
515 im->xo *= 0.001 ; im->dx *= 0.001 ;
516 }
517 }
518 } else {
519 for( kk=0 ; kk < ns ; kk++ ){
520 im = IMARR_SUBIMAGE(imar,kk) ;
521 im->xo = 0.0 ; im->dx = 1.0 ;
522 }
523 }
524
525 free(far) ; RETURN(imar);
526 }
527
528 /*---------------------------------------------------------------------------*/
529 /* Convert time series dataset to a 1D style MRI_IMAGE struct. [05 Mar 2008] */
530
THD_dset_to_1Dmri(THD_3dim_dataset * dset)531 MRI_IMAGE * THD_dset_to_1Dmri( THD_3dim_dataset *dset )
532 {
533 MRI_IMAGE *im ; float *far ;
534 int nx , ny , ii ;
535
536 ENTRY("THD_dset_to_1D") ;
537
538 if( !ISVALID_DSET(dset) ) RETURN(NULL) ;
539 DSET_load(dset) ;
540 if( !DSET_LOADED(dset) ) RETURN(NULL) ;
541
542 nx = DSET_NVALS(dset) ;
543 ny = DSET_NVOX(dset) ;
544 im = mri_new( nx , ny , MRI_float ) ; far = MRI_FLOAT_PTR(im) ;
545
546 for( ii=0 ; ii < ny ; ii++ )
547 THD_extract_array( ii , dset , 0 , far + ii*nx ) ;
548
549 RETURN(im) ;
550 }
551
552 /*----------------------------------------------------------------------------*/
553
THD_extract_many_arrays(int ns,int * ind,THD_3dim_dataset * dset,float * dsar)554 void THD_extract_many_arrays( int ns , int *ind ,
555 THD_3dim_dataset *dset , float *dsar )
556 {
557 int nv , ival , kk ;
558 char *iar ; /* brick in the input */
559 float **far ; /* ptrs to output */
560 float fac ;
561
562 ENTRY("THD_extract_many_arrays") ;
563
564 if( ns <= 0 || ind == NULL || dset == NULL || dsar == NULL ) EXRETURN ;
565
566 /* try to load dataset */
567
568 DSET_load(dset) ; if( !DSET_LOADED(dset) ) EXRETURN ;
569 nv = dset->dblk->nvals ;
570
571 far = (float **) malloc(sizeof(float *)*ns) ; /* 27 Feb 2003 */
572 NULL_CHECK(far) ;
573 for( kk=0 ; kk < ns ; kk++ ) far[kk] = dsar + (size_t)kk*(size_t)nv ;
574
575 /* fill the output */
576
577 switch( DSET_BRICK_TYPE(dset,0) ){
578
579 default: /* don't know what to do --> return nada */
580 free(far) ; EXRETURN ;
581
582 case MRI_byte:{
583 byte *bar ;
584 for( ival=0 ; ival < nv ; ival++ ){
585 bar = (byte *)DSET_ARRAY(dset,ival) ;
586 for( kk=0 ; kk < ns ; kk++ ) far[kk][ival] = (float)bar[ind[kk]] ;
587 }
588 }
589 break ;
590
591 case MRI_short:{
592 short *bar ;
593 for( ival=0 ; ival < nv ; ival++ ){
594 bar = (short *)DSET_ARRAY(dset,ival) ;
595 for( kk=0 ; kk < ns ; kk++ ) far[kk][ival] = (float)bar[ind[kk]] ;
596 }
597 }
598 break ;
599
600 case MRI_float:{
601 float *bar ;
602 for( ival=0 ; ival < nv ; ival++ ){
603 bar = (float *)DSET_ARRAY(dset,ival) ;
604 for( kk=0 ; kk < ns ; kk++ ) far[kk][ival] = bar[ind[kk]] ;
605 }
606 }
607 break ;
608
609 #if 0
610 case MRI_int:{
611 int *bar ;
612 for( ival=0 ; ival < nv ; ival++ ){
613 bar = (int *)DSET_ARRAY(dset,ival) ;
614 for( kk=0 ; kk < ns ; kk++ ) far[kk][ival] = (float)bar[ind[kk]] ;
615 }
616 }
617 break ;
618
619 case MRI_double:{
620 double *bar ;
621 for( ival=0 ; ival < nv ; ival++ ){
622 bar = (double *)DSET_ARRAY(dset,ival) ;
623 for( kk=0 ; kk < ns ; kk++ ) far[kk][ival] = (float)bar[ind[kk]] ;
624 }
625 }
626 break ;
627 #endif
628
629 case MRI_complex:{
630 complex *bar ;
631 for( ival=0 ; ival < nv ; ival++ ){
632 bar = (complex *)DSET_ARRAY(dset,ival) ;
633 for( kk=0 ; kk < ns ; kk++ ) far[kk][ival] = bar[ind[kk]].r ;
634 }
635 }
636 break ;
637
638 }
639
640 /* scale outputs, if needed */
641
642 for( ival=0 ; ival < nv ; ival++ ){
643 fac = DSET_BRICK_FACTOR(dset,ival) ;
644 if( fac > 0.0f && fac != 1.0f ){
645 for( kk=0 ; kk < ns ; kk++ ) far[kk][ival] *= fac ;
646 }
647 }
648
649 thd_floatscan( (size_t)(ns)*(size_t)nv , dsar ) ; /* 10 Jun 2021 */
650
651 free(far) ; EXRETURN ;
652 }
653