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 "afni.h"
8
9 #ifndef ALLOW_PLUGINS
10 # error "Plugins not properly set up -- see machdep.h"
11 #endif
12
13 /***********************************************************************
14 Plugin to compute power spectrum of a 3D+time dataset.
15 This is a moderately complex example, showing how to deal
16 with different data types on input and output.
17 ************************************************************************/
18
19 /*------------- string to 'help' the user -------------*/
20
21 static char helpstring[] =
22 " Purpose: Compute 'Power Spectrum' of a 3D+time dataset.\n"
23 " Input items are:\n"
24 " Input = 3D+time dataset to analyze\n"
25 "\n"
26 " Output: Prefix = Filename prefix for new dataset\n"
27 " Datum = How to store results\n"
28 "\n"
29 " Ignore Count = How many points to ignore at start\n"
30 " Taper Percent = Amount of data to taper (Hamming)\n"
31 " FFT Length = Fourier transform size to use [N]\n"
32 " (If N > size of data, data will be zero)\n"
33 " (padded. 'shortest' means to use N just)\n"
34 " (above the length of the time series. )\n"
35 "\n"
36 " The output dataset will be stored in the 3D+time format, with\n"
37 " the 'time' index actually being frequency. The frequency grid\n"
38 " spacing will be 1/(N*dt), where N=FFT length and dt = input\n"
39 " dataset time spacing.\n"
40 "\n"
41 " The method used is the simplest known: squared periodogram.\n"
42 " A single FFT is done (i.e., each point has DOF=2.)\n"
43 ;
44
45 /*------------- strings for output format -------------*/
46
47 static char * type_strings[]
48 = { "as Input" , "Byte" , "Short" , "Float" } ;
49
50 #define NUM_TYPE_STRINGS (sizeof(type_strings)/sizeof(char *))
51
52 /*------------- strings for FFT length -------------*/
53
54 static char * fft_strings[] =
55 #if 0
56 { "shortest", "32", "64", "128", "256", "512", "1024", "2048", "4096" } ;
57 #else
58 /* 3* 15* 2** 5* */
59 { "shortest", "32" , "48" , "60" , "64" , "80" ,
60 "96" , "120" , "128" , "160" ,
61 "192" , "240" , "256" , "320" ,
62 "384" , "480" , "512" , "640" ,
63 "768" , "960" , "1024" , "1280" ,
64 "1536" , "1920" , "2048" } ;
65 #endif
66
67 #define NUM_FFT_STRINGS (sizeof(fft_strings)/sizeof(char *))
68
69 /*--------------- prototypes for internal routines ---------------*/
70
71 char * POWER_main( PLUGIN_interface * ) ; /* the entry point */
72
73 #undef ALLOW_TESTING
74 #ifdef ALLOW_TESTING
75 PLUGIN_interface * TEST_init(void) ;
76 char * TEST_main( PLUGIN_interface * ) ; /* the entry point */
77 #endif
78
79 /***********************************************************************
80 Set up the interface to the user:
81 1) Create a new interface using "PLUTO_new_interface";
82
83 2) For each line of inputs, create the line with "PLUTO_add_option"
84 (this line of inputs can be optional or mandatory);
85
86 3) For each item on the line, create the item with
87 "PLUTO_add_dataset" for a dataset chooser,
88 "PLUTO_add_string" for a string chooser,
89 "PLUTO_add_number" for a number chooser.
90 ************************************************************************/
91
92
93 DEFINE_PLUGIN_PROTOTYPE
94
PLUGIN_init(int ncall)95 PLUGIN_interface * PLUGIN_init( int ncall )
96 {
97 PLUGIN_interface * plint ; /* will be the output of this routine */
98
99 if( ncall > 1 ) return NULL ; /* two interfaces */
100
101 #ifdef ALLOW_TESTING
102 if( ncall == 1 ) return TEST_init() ;
103 #else
104 if( ncall == 1 ) return NULL ;
105 #endif
106
107 /*---------------- set titles and call point ----------------*/
108
109 CHECK_IF_ALLOWED("POWERSPECTRUM","Power Spectrum") ; /* 30 Sep 2016 */
110
111 plint = PLUTO_new_interface( "Power Spectrum" ,
112 "Power Spectrum of a 3D+time Dataset" ,
113 helpstring ,
114 PLUGIN_CALL_VIA_MENU , POWER_main ) ;
115
116 PLUTO_add_hint( plint , "Power Spectrum of a 3D+time Dataset" ) ;
117
118 PLUTO_set_sequence( plint , "A:newdset:statistics" ) ;
119
120 /*--------- 1st line: Input dataset ---------*/
121
122 PLUTO_add_option( plint ,
123 "Input" , /* label at left of input line */
124 "Input" , /* tag to return to plugin */
125 TRUE /* is this mandatory? */
126 ) ;
127
128 PLUTO_add_dataset( plint ,
129 "---->>" , /* label next to button */
130 ANAT_ALL_MASK , /* take any anat datasets */
131 FUNC_FIM_MASK , /* only allow fim funcs */
132 DIMEN_4D_MASK | /* need 3D+time datasets */
133 BRICK_ALLREAL_MASK /* need real-valued datasets */
134 ) ;
135
136 /*---------- 2nd line: Output dataset ----------*/
137
138 PLUTO_add_option( plint ,
139 "Output" , /* label at left of input line */
140 "Output" , /* tag to return to plugin */
141 TRUE /* is this mandatory? */
142 ) ;
143
144 PLUTO_add_string( plint ,
145 "Prefix" , /* label next to textfield */
146 0,NULL , /* no fixed strings to choose among */
147 19 /* 19 spaces for typing in value */
148 ) ;
149
150 PLUTO_add_string( plint ,
151 "Datum" , /* label next to chooser button */
152 NUM_TYPE_STRINGS , /* number of strings to choose among */
153 type_strings , /* list of strings to choose among */
154 0 /* index of default string */
155 ) ;
156
157 /*--------- Other lines: Parameters ---------*/
158
159 PLUTO_add_option( plint , "Ignore" , "Ignore" , TRUE ) ;
160
161 PLUTO_add_number( plint ,
162 "Count" , /* label next to chooser */
163 0 , /* smallest possible value */
164 999 , /* largest possible value */
165 0 , /* decimal shift (none in this case) */
166 4 , /* default value */
167 TRUE /* allow user to edit value? */
168 ) ;
169
170 PLUTO_add_option( plint , "Taper" , "Taper" , TRUE ) ;
171
172 PLUTO_add_number( plint ,
173 "Percent" , /* label next to chooser */
174 0 , /* smallest possible value */
175 10 , /* largest possible value */
176 -1 , /* decimal shift (1 right == 0 to 100) */
177 0 , /* default value (with shift == 0) */
178 FALSE /* allow user to edit value? */
179 ) ;
180
181 PLUTO_add_option( plint , "FFT" , "FFT" , TRUE ) ;
182
183 PLUTO_add_string( plint ,
184 "Length" , /* label next to chooser */
185 NUM_FFT_STRINGS , /* number of strings to choose among */
186 fft_strings , /* list of strings to choose among */
187 0 /* index of default string */
188 ) ;
189
190 /*--------- done with interface setup ---------*/
191
192 return plint ;
193 }
194
195 /***************************************************************************
196 Main routine for this plugin (will be called from AFNI).
197 If the return string is not NULL, some error transpired, and
198 AFNI will popup the return string in a message box.
199 ****************************************************************************/
200
201 /*------------------ macros to return workspace at exit -------------------*/
202
203 #undef FREEUP
204 #define FREEUP(x) do{ if((x) != NULL){free((x)); (x)=NULL;} } while(0)
205
206 #define FREE_WORKSPACE \
207 do{ FREEUP(bptr) ; FREEUP(sptr) ; FREEUP(fptr) ; \
208 FREEUP(fout) ; FREEUP(cxar) ; FREEUP(tar) ; \
209 FREEUP(fxar) ; FREEUP(fyar) ; FREEUP(dtr) ; \
210 } while(0)
211
212 /*-------------------------------------------------------------------------*/
213
POWER_main(PLUGIN_interface * plint)214 char * POWER_main( PLUGIN_interface * plint )
215 {
216 MCW_idcode * idc ; /* input dataset idcode */
217 THD_3dim_dataset * old_dset , * new_dset ; /* input and output datasets */
218 char * new_prefix , * str ; /* strings from user */
219 int new_datum , ignore , nfft , ninp , /* control parameters */
220 old_datum , nuse , ntaper , ktbot ;
221 float taper ;
222
223 byte ** bptr = NULL ; /* one of these will be the array of */
224 short ** sptr = NULL ; /* pointers to input dataset sub-bricks */
225 float ** fptr = NULL ; /* (depending on input datum type) */
226
227 complex * cxar = NULL ; /* will be array of data to FFT */
228 float * fxar = NULL ; /* array loaded from input dataset */
229 float * fyar = NULL ; /* array loaded from input dataset */
230 float ** fout = NULL ; /* will be array of output floats */
231
232 float * tar = NULL ; /* will be array of taper coefficients */
233 float * dtr = NULL ; /* will be array of detrending coeff */
234
235 float dfreq , pfact , phi , xr,xi , yr,yi ;
236 float x0,x1 , y0,y1 , d0fac,d1fac ;
237 int nfreq , nvox , perc , new_units ;
238 int istr , ii,iip , ibot,itop , kk , icx ; /* temp variables */
239
240 /*--------------------------------------------------------------------*/
241 /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
242
243 /*--------- go to first input line ---------*/
244
245 PLUTO_next_option(plint) ;
246
247 idc = PLUTO_get_idcode(plint) ; /* get dataset item */
248 old_dset = PLUTO_find_dset(idc) ; /* get ptr to dataset */
249 if( old_dset == NULL )
250 return "*************************\n"
251 "Cannot find Input Dataset\n"
252 "*************************" ;
253
254 /*--------- go to second input line ---------*/
255
256 PLUTO_next_option(plint) ;
257
258 new_prefix = PLUTO_get_string(plint) ; /* get string item (the output prefix) */
259 if( ! PLUTO_prefix_ok(new_prefix) ) /* check if it is OK */
260 return "************************\n"
261 "Output Prefix is illegal\n"
262 "************************" ;
263
264 str = PLUTO_get_string(plint) ; /* get string item (the datum type) */
265 istr = PLUTO_string_index( str , /* find it in the list it came from */
266 NUM_TYPE_STRINGS ,
267 type_strings ) ;
268 switch( istr ){
269 default:
270 case 0:
271 new_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* use old dataset type */
272 break ;
273
274 case 1: new_datum = MRI_byte ; break ; /* assign type of user's choice */
275 case 2: new_datum = MRI_short ; break ;
276 case 3: new_datum = MRI_float ; break ;
277 }
278
279 /*--------- go to next input lines ---------*/
280
281 PLUTO_next_option(plint) ; /* skip to next line */
282 ignore = PLUTO_get_number(plint) ; /* get number item (ignore) */
283
284 PLUTO_next_option(plint) ; /* skip to next line */
285 taper = PLUTO_get_number(plint) * 0.01 ; /* get number item (taper %) */
286
287 /* compute FFT length to use */
288
289 PLUTO_next_option(plint) ; /* skip to next line */
290
291 str = PLUTO_get_string(plint) ; /* get string item for FFT count */
292 ninp = DSET_NUM_TIMES(old_dset) ; /* number of values in input */
293 nuse = ninp - ignore ; /* number of values to actually use */
294
295 if( nuse < 4 )
296 return "*****************************\n"
297 "Not enough time points to FFT\n"
298 "*****************************" ;
299
300 if( strcmp(str,fft_strings[0]) == 0 ){
301
302 /*-- get next larger power-of-2 --*/
303 #if 0
304 for( nfft=32 ; nfft < nuse ; nfft *= 2 ) ; /* loop until nfft >= nuse */
305 #else
306 nfft = csfft_nextup_even(nuse) ;
307 #endif
308
309 } else {
310 nfft = strtol( str , NULL , 10 ) ; /* just convert string to integer */
311 }
312
313 /* if the input FFT length is less than the data length,
314 tell the user and truncate the amount of data to use */
315
316 if( nfft < nuse ){
317 char str[256] ;
318
319 sprintf( str , "******************************\n"
320 "Warning:\n"
321 " Number of points in FFT =%4d\n"
322 " is less than available data\n"
323 " in time series = %d\n"
324 "******************************" ,
325 nfft , nuse ) ;
326
327 PLUTO_popup_transient( plint , str ) ;
328
329 nuse = nfft ; /* can't use more data than the FFT length */
330 }
331
332 /* compute the number of output points and the output grid spacing */
333
334 nfreq = nfft / 2 ; /* # frequencies */
335 dfreq = 1.0 / (nfft * DSET_TIMESTEP(old_dset) ) ; /* frequency grid */
336
337 switch( DSET_TIMEUNITS(old_dset) ){
338 case UNITS_MSEC_TYPE: dfreq *= 1000.0 ; new_units = UNITS_HZ_TYPE ; break;
339 case UNITS_SEC_TYPE: new_units = UNITS_HZ_TYPE ; break;
340 case UNITS_HZ_TYPE: new_units = UNITS_SEC_TYPE; break;
341
342 default: new_units = DSET_TIMEUNITS(old_dset) ; break ; /* shouldn't happen */
343 }
344
345 /*------------------------------------------------------*/
346 /*---------- At this point, the inputs are OK ----------*/
347
348 PLUTO_popup_meter( plint ) ; /* popup a progress meter */
349
350 /*--------- set up pointers to each sub-brick in the input dataset ---------*/
351
352 DSET_load( old_dset ) ; /* must be in memory before we get pointers to it */
353
354 old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum type */
355
356 switch( old_datum ){ /* pointer type depends on input datum type */
357
358 default:
359 return "******************************\n"
360 "Illegal datum in Input Dataset\n"
361 "******************************" ;
362
363 /** create array of pointers into old dataset sub-bricks **/
364 /** Note that we skip the first 'ignore' sub-bricks here **/
365
366 /*--------- input is bytes ----------*/
367 /* voxel #i at time #k is bptr[k][i] */
368 /* for i=0..nvox-1 and k=0..nuse-1. */
369
370 case MRI_byte:
371 bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;
372 if( bptr == NULL ) return "Malloc\nFailure!\n [bptr]" ;
373 for( kk=0 ; kk < nuse ; kk++ )
374 bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk+ignore) ;
375 break ;
376
377 /*--------- input is shorts ---------*/
378 /* voxel #i at time #k is sptr[k][i] */
379 /* for i=0..nvox-1 and k=0..nuse-1. */
380
381 case MRI_short:
382 sptr = (short **) malloc( sizeof(short *) * nuse ) ;
383 if( sptr == NULL ) return "Malloc\nFailure!\n [sptr]" ;
384 for( kk=0 ; kk < nuse ; kk++ )
385 sptr[kk] = (short *) DSET_ARRAY(old_dset,kk+ignore) ;
386 break ;
387
388 /*--------- input is floats ---------*/
389 /* voxel #i at time #k is fptr[k][i] */
390 /* for i=0..nvox-1 and k=0..nuse-1. */
391
392 case MRI_float:
393 fptr = (float **) malloc( sizeof(float *) * nuse ) ;
394 if( fptr == NULL ) return "Malloc\nFailure!\n [fptr]" ;
395 for( kk=0 ; kk < nuse ; kk++ )
396 fptr[kk] = (float *) DSET_ARRAY(old_dset,kk+ignore) ;
397 break ;
398
399 } /* end of switch on input type */
400
401 /*---- allocate space for 2 voxel timeseries and 1 FFT ----*/
402
403 cxar = (complex *) malloc( sizeof(complex) * nfft ) ; /* FFT */
404 fxar = (float *) malloc( sizeof(float) * nuse ) ; /* input */
405 fyar = (float *) malloc( sizeof(float) * nuse ) ; /* input */
406 if( cxar == NULL || fxar == NULL || fyar == NULL ){
407 FREE_WORKSPACE ;
408 return "Malloc\nFailure!\n [cxar]" ;
409 }
410
411 /*--------- make space for taper coefficient array ---------*/
412
413 tar = (float *) malloc( sizeof(float) * MAX(nuse,nfreq) ) ;
414 dtr = (float *) malloc( sizeof(float) * nuse ) ;
415
416 if( tar == NULL || dtr == NULL ){
417 FREE_WORKSPACE ;
418 return "Malloc\nFailure!\n [tar]" ;
419 }
420
421 ntaper = (int)(0.5 * taper * nuse + 0.49) ; /* will taper data over */
422 phi = PI / MAX(ntaper,1) ; /* kk=0..ntaper-1 on left */
423 ktbot = nuse - ntaper ; /* kk=ktbot..nuse-1 on right */
424 pfact = 0.0 ; /* sum of taper**2 */
425
426 for( kk=0 ; kk < nuse ; kk++ ){ /* Hamming-ize */
427 if( kk < ntaper )
428 tar[kk] = 0.54 - 0.46 * cos(kk*phi) ; /* ramp up */
429 else if( kk >= ktbot )
430 tar[kk] = 0.54 + 0.46 * cos((kk-ktbot+1)*phi) ; /* ramp down */
431 else
432 tar[kk] = 1.0 ; /* in the middle */
433
434 pfact += tar[kk] * tar[kk] ;
435
436 dtr[kk] = kk - 0.5 * (nuse-1) ; /* factors for linear detrending */
437 }
438
439 d0fac = 1.0 / nuse ;
440 d1fac = 12.0 / nuse / (nuse*nuse - 1.0) ;
441
442 /*--- compute factor to go from |FFT|**2 to PSD;
443 includes the scaling needed for loss of energy with tapering ---*/
444
445 pfact = DSET_TIMESTEP(old_dset) / pfact ;
446
447 /*--- include scaling factors for sub-bricks, if any ---*/
448
449 for( kk=0 ; kk < nuse ; kk++ )
450 if( DSET_BRICK_FACTOR(old_dset,kk+ignore) > 0.0 )
451 tar[kk] *= DSET_BRICK_FACTOR(old_dset,kk+ignore) ;
452
453 /*---------------------- make a new dataset ----------------------*/
454
455 new_dset = EDIT_empty_copy( old_dset ) ; /* start with copy of old one */
456
457 { char * his = PLUTO_commandstring(plint) ;
458 tross_Copy_History( old_dset , new_dset ) ;
459 tross_Append_History( new_dset , his ) ; free(his) ;
460 }
461
462 /*-- edit some of its internal parameters --*/
463
464 ii = EDIT_dset_items(
465 new_dset ,
466 ADN_prefix , new_prefix , /* filename prefix */
467 ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
468 ADN_datum_all , new_datum , /* atomic datum */
469 ADN_nvals , nfreq , /* # sub-bricks */
470 ADN_ntt , nfreq , /* # time points */
471 ADN_ttorg , dfreq , /* time origin */
472 ADN_ttdel , dfreq , /* time step */
473 ADN_ttdur , dfreq , /* time duration */
474 ADN_nsl , 0 , /* z-axis time slicing */
475 ADN_tunits , new_units , /* time units */
476 ADN_none ) ;
477
478 if( ii != 0 ){
479 THD_delete_3dim_dataset( new_dset , False ) ;
480 FREE_WORKSPACE ;
481 return "***********************************\n"
482 "Error while creating output dataset\n"
483 "***********************************" ;
484 }
485
486 /*------ make floating point output sub-bricks
487 (only at the end will scale to byte or shorts)
488
489 Output #ii at freq #kk will go into fout[kk][ii],
490 for kk=0..nfreq-1, and for ii=0..nvox-1. ------*/
491
492 nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
493
494 fout = (float **) malloc( sizeof(float *) * nfreq ) ; /* ptrs to sub-bricks */
495
496 if( fout == NULL ){
497 THD_delete_3dim_dataset( new_dset , False ) ;
498 FREE_WORKSPACE ;
499 return "Malloc\nFailure!\n [fout]" ;
500 }
501
502 for( kk=0 ; kk < nfreq ; kk++ ){
503 fout[kk] = (float *) malloc( sizeof(float) * nvox ) ; /* sub-brick # kk */
504 if( fout[kk] == NULL ) break ;
505 }
506
507 if( kk < nfreq ){
508 for( ; kk >= 0 ; kk-- ) FREEUP(fout[kk]) ; /* free all we did get */
509 THD_delete_3dim_dataset( new_dset , False ) ;
510 FREE_WORKSPACE ;
511 return "Malloc\nFailure!\n [arrays]" ;
512 }
513
514 { char buf[128] ;
515 ii = (nfreq * nvox * sizeof(float)) / (1024*1024) ;
516 sprintf( buf , " \n"
517 "*** 3D+time Power Spectral Density:\n"
518 "*** Using %d MBytes of workspace,\n "
519 "*** with FFT length = %d\n" , ii,nfft ) ;
520 PLUTO_popup_transient( plint , buf ) ;
521 }
522
523 /*----------------------------------------------------*/
524 /*----- Setup has ended. Now do some real work. -----*/
525
526 /***** loop over voxels *****/
527
528 for( ii=0 ; ii < nvox ; ii += 2 ){ /* 2 time series at a time */
529
530 iip = (ii+1) % nvox ; /* voxel after ii */
531
532 /*** load data from input dataset, depending on type ***/
533
534 switch( old_datum ){
535
536 /*** input = bytes ***/
537
538 case MRI_byte:
539 for( kk=0 ; kk < nuse ; kk++ ){
540 fxar[kk] = bptr[kk][ii] ;
541 fyar[kk] = bptr[kk][iip] ;
542 }
543 break ;
544
545 /*** input = shorts ***/
546
547 case MRI_short:
548 for( kk=0 ; kk < nuse ; kk++ ){
549 fxar[kk] = sptr[kk][ii] ;
550 fyar[kk] = sptr[kk][iip] ;
551 }
552 break ;
553
554 /*** input = floats ***/
555
556 case MRI_float:
557 for( kk=0 ; kk < nuse ; kk++ ){
558 fxar[kk] = fptr[kk][ii] ;
559 fyar[kk] = fptr[kk][iip] ;
560 }
561 break ;
562
563 } /* end of switch over input type */
564
565 /*** detrend:
566 x0 = sum( fxar[kk] )
567 x1 = sum( fxar[kk] * (kk-0.5*(N-1)) )
568 x0 is used to remove the mean of fxar
569 x1 is used to remove the linear trend of fxar ***/
570
571 x0 = x1 = y0 = y1 = 0.0 ;
572 for( kk=0 ; kk < nuse ; kk++ ){
573 x0 += fxar[kk] ; x1 += fxar[kk] * dtr[kk] ;
574 y0 += fyar[kk] ; y1 += fyar[kk] * dtr[kk] ;
575 }
576
577 x0 *= d0fac ; x1 *= d1fac ; /* factors to remove mean and trend */
578 y0 *= d0fac ; y1 *= d1fac ;
579
580 for( kk=0 ; kk < nuse ; kk++ ){
581 fxar[kk] -= (x0 + x1 * dtr[kk]) ; /* remove mean and trend here! */
582 fyar[kk] -= (y0 + y1 * dtr[kk]) ;
583 }
584
585 /*** taper, scale, and put into cxar array ***/
586
587 for( kk=0 ; kk < nuse ; kk++ ){
588 cxar[kk].r = fxar[kk] * tar[kk] ;
589 cxar[kk].i = fyar[kk] * tar[kk] ;
590 }
591
592 /*** load zeros after where data was put ***/
593
594 for( kk=nuse ; kk < nfft ; kk++ ) cxar[kk].r = cxar[kk].i = 0.0 ;
595
596 /***** do the FFT (at long last) *****/
597
598 csfft_cox( -1 , nfft , cxar ) ;
599
600 /***** now compute output into corresponding voxels in fout *****/
601
602 /*--- Let x = fxar (1st real time series)
603 y = fyar (2nd real time series)
604 z = cxar (complex time series) = x + i y
605 N = nfft (length of FFT)
606
607 Then after FFT, since x and y are real, we have
608 zhat[k] = xhat[k] + i yhat[k] > for k=1..N/2
609 zhat[N-k]* = xhat[k] - i yhat[k]
610
611 so we can untangle the FFTs of x and y by
612 xhat[k] = 0.5 ( zhat[k] + zhat[N-k]* )
613 yhat[k] = 0.5 ( zhat[k] - zhat[N-k]* ) / i
614
615 This is the basis for doing 2 time series at once. ---*/
616
617 for( kk=1 ; kk <= nfreq ; kk++ ){
618 xr = 0.5 * ( cxar[kk].r + cxar[nfft-kk].r ) ; /* Re xhat[kk] */
619 xi = 0.5 * ( cxar[kk].i - cxar[nfft-kk].i ) ; /* Im xhat[kk] */
620 yr = 0.5 * ( cxar[kk].i + cxar[nfft-kk].i ) ; /* Re yhat[kk] */
621 yi = 0.5 * ( cxar[kk].r - cxar[nfft-kk].r ) ; /*-Im yhat[kk] */
622
623 fout[kk-1][ii] = pfact * (xr*xr + xi*xi) ;
624 fout[kk-1][iip] = pfact * (yr*yr + yi*yi) ;
625 }
626
627 perc = (100 * ii) / nvox ; /* display percentage done */
628 PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
629
630 } /* end of outer loop over 2 voxels at a time */
631
632 DSET_unload( old_dset ) ; /* don't need this no more */
633
634 /*------------------------------------------------------------*/
635 /*------- The output is now in fout[kk][ii],
636 for kk=0..nfreq-1 , ii=0..nvox-1.
637 We must now put this into the output dataset -------*/
638
639 switch( new_datum ){
640
641 /*** output is floats is the simplest:
642 we just have to attach the fout bricks to the dataset ***/
643
644 case MRI_float:
645 for( kk=0 ; kk < nfreq ; kk++ )
646 EDIT_substitute_brick( new_dset , kk , MRI_float , fout[kk] ) ;
647 break ;
648
649 /*** output is shorts:
650 we have to create a scaled sub-brick from fout ***/
651
652 case MRI_short:{
653 short * bout ;
654 float fac ;
655
656 for( kk=0 ; kk < nfreq ; kk++ ){ /* loop over sub-bricks */
657
658 /*-- get output sub-brick --*/
659
660 bout = (short *) malloc( sizeof(short) * nvox ) ;
661 if( bout == NULL ){
662 fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
663 return("\nFinal malloc error in plug_power!\n") ;
664 /* EXIT(1) ;*/
665 }
666
667 /*-- find scaling and then scale --*/
668
669 fac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[kk] ) ;
670 if( fac > 0.0 ){
671 fac = 32767.0 / fac ;
672 EDIT_coerce_scale_type( nvox,fac ,
673 MRI_float,fout[kk] , MRI_short,bout ) ;
674 fac = 1.0 / fac ;
675 }
676
677 free( fout[kk] ) ; /* don't need this anymore */
678
679 /*-- put output brick into dataset, and store scale factor --*/
680
681 EDIT_substitute_brick( new_dset , kk , MRI_short , bout ) ;
682 tar[kk] = fac ;
683
684 perc = (100 * kk) / nfreq ;
685 PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
686 }
687
688 /*-- save scale factor array into dataset --*/
689
690 EDIT_dset_items( new_dset , ADN_brick_fac , tar , ADN_none ) ;
691
692 }
693 break ;
694
695 /*** output is bytes (byte = unsigned char)
696 we have to create a scaled sub-brick from fout ***/
697
698 case MRI_byte:{
699 byte * bout ;
700 float fac ;
701
702 for( kk=0 ; kk < nfreq ; kk++ ){ /* loop over sub-bricks */
703
704 /*-- get output sub-brick --*/
705
706 bout = (byte *) malloc( sizeof(byte) * nvox ) ;
707 if( bout == NULL ){
708 fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
709 return("\nFinal malloc error in plug_power!\n\a") ;
710 /* EXIT(1) ;*/
711 }
712
713 /*-- find scaling and then scale --*/
714
715 fac = MCW_vol_amax( nvox,1,1 , MRI_float , fout[kk] ) ;
716 if( fac > 0.0 ){
717 fac = 255.0 / fac ;
718 EDIT_coerce_scale_type( nvox,fac ,
719 MRI_float,fout[kk] , MRI_byte,bout ) ;
720 fac = 1.0 / fac ;
721 }
722
723 free( fout[kk] ) ; /* don't need this anymore */
724
725 /*-- put output brick into dataset, and store scale factor --*/
726
727 EDIT_substitute_brick( new_dset , kk , MRI_byte , bout ) ;
728 tar[kk] = fac ;
729
730 perc = (100 * kk) / nfreq ;
731 PLUTO_set_meter( plint , perc ) ; /* on the progress meter */
732 }
733
734 /*-- save scale factor array into dataset --*/
735
736 EDIT_dset_items( new_dset , ADN_brick_fac , tar , ADN_none ) ;
737 }
738 break ;
739
740 } /* end of switch on output data type */
741
742 /*-------------- Cleanup and go home ----------------*/
743
744 PLUTO_set_meter( plint , 100 ) ; /* set progress meter to 100% */
745
746 PLUTO_add_dset( plint , new_dset , DSET_ACTION_MAKE_CURRENT ) ;
747
748 FREE_WORKSPACE ;
749 return NULL ; /* null string returned means all was OK */
750 }
751
752 #ifdef ALLOW_TESTING
753 /*****************************************************************************
754 -----------------------------------------------------------------------------
755 Create the second interface within this plugin.
756 -----------------------------------------------------------------------------*/
757
TEST_init(void)758 PLUGIN_interface * TEST_init( void )
759 {
760 PLUGIN_interface * plint ; /* will be the output of this routine */
761
762 /*---------------- set titles and call point ----------------*/
763
764 plint = PLUTO_new_interface( "Testing" ,
765 "Testing, Testing, 1-2-3 ..." ,
766 NULL ,
767 PLUGIN_CALL_VIA_MENU , TEST_main ) ;
768
769 PLUTO_add_hint( plint , "1-2-3, 1-2-3, ..." ) ;
770
771 PLUTO_add_option( plint ,
772 "Input" , /* label at left of input line */
773 "Input" , /* tag to return to plugin */
774 TRUE /* is this mandatory? */
775 ) ;
776
777 PLUTO_add_dataset_list( plint ,
778 "Datasets" , /* label next to button */
779 ANAT_ALL_MASK , /* take any anat datasets */
780 FUNC_FIM_MASK , /* only allow fim funcs */
781 DIMEN_4D_MASK | /* need 3D+time datasets */
782 BRICK_ALLREAL_MASK /* need real-valued datasets */
783 ) ;
784 return plint ;
785 }
786
TEST_main(PLUGIN_interface * plint)787 char * TEST_main( PLUGIN_interface * plint )
788 {
789 MRI_IMAGE * tsim ;
790 MCW_idclist * idclist ;
791 MCW_idcode * idc ;
792 THD_3dim_dataset * dset ;
793 char str[256] ;
794 int id ;
795
796 /*--------- go to first input line ---------*/
797
798 PLUTO_next_option(plint) ;
799
800 idclist = PLUTO_get_idclist(plint) ;
801 if( PLUTO_idclist_count(idclist) == 0 )
802 return " \nNo input dataset list!\n " ;
803
804 id = 0 ;
805 do {
806 idc = PLUTO_idclist_next(idclist) ;
807 dset = PLUTO_find_dset(idc) ;
808 if( dset == NULL ) return NULL ;
809 id++ ;
810 sprintf(str, " \nDataset %d = %s\n nx = %d\n ny = %d\n nz = %d\n " ,
811 id , DSET_FILECODE(dset) , dset->daxes->nxx,dset->daxes->nyy,dset->daxes->nzz ) ;
812
813 PLUTO_popup_transient( plint , str ) ;
814 } while(1) ;
815 return NULL ;
816 }
817 #endif /* ALLOW_TESTING */
818