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 Functional Imageof 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 'ASL a3/d3' 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 time points to ignore at start\n"
30   " \n"
31   "   Scale:  Scale the a3, d3 calculations\n"
32   " \n"
33   "   Notes: \n"
34   "    a3- three point moving average of the time series (BOLD) \n"
35   "                    a3(i)=image(i-1)+2*image(i)+image(i+1) \n"
36   "    d3- three point moving difference of the time series (Perfusion) \n"
37   "         d3(i)=(-1*i)*(image(i-1)-2*image(i)+1*image(i+1)) \n"
38   "\n"
39   "    Datum Type- be aware that rounding errors associated with scaling to \n"
40   "                byte/short can occur \n"
41   "   \n"
42   "\n"
43   " ASL A3/D3 Version 1.3 Written by Y.Behzadi 8/18/02\n"
44   " fixed problem with odd # of time points and sign\n"
45   "\n"
46 
47 ;
48 
49 /*------------- strings for output format -------------*/
50 
51 static char * type_strings[]
52   = { "Float" ,  "Byte" , "Short" , "as Input"  } ;
53 
54 static char * type_stringsx[]
55   = { "a3,d3,avg(d3)" } ;
56 
57 #define NUM_TYPE_STRINGS (sizeof(type_strings)/sizeof(char *))
58 #define NUM_TYPE_STRINGSX (sizeof(type_stringsx)/sizeof(char *))
59 
60 /*--------------- prototypes for internal routines ---------------*/
61 
62 char * POWER_main( PLUGIN_interface * ) ;  /* the entry point */
63 
64 #undef ALLOW_TESTING
65 #ifdef ALLOW_TESTING
66 PLUGIN_interface * TEST_init(void) ;
67 char * TEST_main( PLUGIN_interface * ) ;  /* the entry point */
68 #endif
69 
70 /***********************************************************************
71    Set up the interface to the user:
72     1) Create a new interface using "PLUTO_new_interface";
73 
74     2) For each line of inputs, create the line with "PLUTO_add_option"
75          (this line of inputs can be optional or mandatory);
76 
77     3) For each item on the line, create the item with
78         "PLUTO_add_dataset" for a dataset chooser,
79         "PLUTO_add_string"  for a string chooser,
80         "PLUTO_add_number"  for a number chooser.
81 ************************************************************************/
82 
PLUGIN_init(int ncall)83 PLUGIN_interface * PLUGIN_init( int ncall )
84 {
85    PLUGIN_interface * plint ;     /* will be the output of this routine */
86 
87    if( ncall > 1 ) return NULL ;  /* two interfaces */
88    CHECK_IF_ALLOWED("ASL","ASL a3/d3") ;  /* 30 Sep 2016 */
89 
90 #ifdef ALLOW_TESTING
91    if( ncall == 1 ) return TEST_init() ;
92 #else
93    if( ncall == 1 ) return NULL ;
94 #endif
95 
96    /*---------------- set titles and call point ----------------*/
97 
98    plint = PLUTO_new_interface( "ASL a3/d3" ,
99                                 "ASL a3/d3 of a 3D+time Dataset" ,
100                                 helpstring ,
101                                 PLUGIN_CALL_VIA_MENU , POWER_main  ) ;
102 
103    PLUTO_add_hint( plint , "ASL a3/d3 of a 3D+time Dataset" ) ;
104 
105    PLUTO_set_sequence( plint , "A:newdset:statistics" ) ;
106 
107    /*--------- 1st line: Input dataset ---------*/
108 
109    PLUTO_add_option( plint ,
110                      "Input" ,  /* label at left of input line */
111                      "Input" ,  /* tag to return to plugin */
112                      TRUE       /* is this mandatory? */
113                    ) ;
114 
115    PLUTO_add_dataset(  plint ,
116                        "---->>" ,         /* label next to button   */
117                        ANAT_ALL_MASK ,    /* take any anat datasets */
118                        FUNC_FIM_MASK ,    /* only allow fim funcs   */
119                        DIMEN_4D_MASK |    /* need 3D+time datasets  */
120                        BRICK_ALLREAL_MASK /* need real-valued datasets */
121                     ) ;
122 
123    /*---------- 2nd line: Output dataset ----------*/
124 
125    PLUTO_add_option( plint ,
126                      "Output" ,  /* label at left of input line */
127                      "Output" ,  /* tag to return to plugin */
128                      TRUE        /* is this mandatory? */
129                    ) ;
130 
131    PLUTO_add_string(   plint ,
132                        "Prefix" ,  /* label next to textfield */
133                        0,NULL ,    /* no fixed strings to choose among */
134                        19          /* 19 spaces for typing in value */
135                    ) ;
136 
137    PLUTO_add_string(   plint ,
138                        "Datum" ,          /* label next to chooser button */
139                        NUM_TYPE_STRINGS , /* number of strings to choose among */
140                        type_strings ,     /* list of strings to choose among */
141                        0                  /* index of default string */
142                    ) ;
143 
144    /*--------- Other lines: Parameters ---------*/
145 
146    PLUTO_add_option( plint , "Ignore" , "Ignore" , TRUE ) ;
147 
148    PLUTO_add_number( plint ,
149                      "Count" ,   /* label next to chooser */
150                      0 ,         /* smallest possible value */
151                      999 ,       /* largest possible value */
152                      0 ,         /* decimal shift (none in this case) */
153                      4 ,         /* default value */
154                      TRUE        /* allow user to edit value? */
155                    ) ;
156    PLUTO_add_string(   plint ,
157                        "Image Output" ,     /* label next to chooser button */
158                        NUM_TYPE_STRINGSX , /* number of strings to choose among */
159                        type_stringsx ,     /* list of strings to choose among */
160                        0                  /* index of default string */
161                    ) ;
162 
163    PLUTO_add_option( plint , "Scale" , "Scale" , TRUE ) ;
164 
165    PLUTO_add_number( plint ,
166                      "Scale Factor" ,   /* label next to chooser */
167                      0 ,         /* smallest possible value */
168                      999 ,       /* largest possible value */
169                      0 ,         /* decimal shift (none in this case) */
170                      1 ,         /* default value */
171                      TRUE        /* allow user to edit value? */
172                    ) ;
173 
174    return plint ;
175 }
176 
177 /***************************************************************************
178   Main routine for this plugin (will be called from AFNI).
179   If the return string is not NULL, some error transpired, and
180   AFNI will popup the return string in a message box.
181 ****************************************************************************/
182 
183 /*------------------ macros to return workspace at exit -------------------*/
184 
185 #undef  FREEUP
186 #define FREEUP(x) do{ if((x) != NULL){free((x)); (x)=NULL;} } while(0)
187 
188 #define FREE_WORKSPACE  do{ FREEUP(bptr) ;      \
189                             FREEUP(sptr) ;      \
190                             FREEUP(fptr) ;      \
191                             FREEUP(foutD3) ;    \
192                             FREEUP(this) ;      \
193                             FREEUP(tarD3) ;     \
194                             FREEUP(tarA3) ;     \
195                             FREEUP(taravgD3) ;  \
196                             FREEUP(foutA3) ;    \
197                             FREEUP(foutavgD3) ; \
198                         } while(0)
199 
200 /*-------------------------------------------------------------------------*/
201 
POWER_main(PLUGIN_interface * plint)202 char * POWER_main( PLUGIN_interface * plint )
203 {
204    MCW_idcode * idc ;                          /* input dataset idcode */
205    THD_3dim_dataset * old_dset , * new_dsetD3 , * new_dsetA3, * new_dsetavgD3 ; /* input and output datasets */
206    char * new_prefix , * str , * namestr, * filename;                 /* strings from user */
207    int   new_datum , ignore , nfft , ninp ,    /* control parameters */
208          old_datum , nuse , ntaper , ktbot,
209          image_type, scale,OutputFlag ,numT,flip;
210   float avFac;
211 
212    byte   ** bptr  = NULL ;  /* one of these will be the array of */
213    short  ** sptr  = NULL ;  /* pointers to input dataset sub-bricks */
214    float  ** fptr  = NULL ;  /* (depending on input datum type) */
215 
216 
217 
218    float   * this  = NULL ;  /* array loaded from input dataset */
219 
220 
221    float  ** foutD3  = NULL ;  /* will be array of output floats */
222    float  ** foutA3  = NULL ;  /* will be array of output floats */
223    float  ** foutavgD3  = NULL ;  /* will be array of output floats */
224 
225    float   * tarD3   = NULL ;  /* will be array of taper coefficients */
226    float   * tarA3   = NULL ;  /* will be array of taper coefficients */
227    float   * taravgD3   = NULL ;  /* will be array of taper coefficients */
228 
229 
230    /*float   * flip;*/
231    float   * numAv;
232    float dfreq , pfact , phi , xr,xi , yr,yi ;
233    float x0,x1 , y0,y1 , d0fac,d1fac ;
234    int   nfreq , nvox , perc , new_units ;
235    int   istr , ii,iip , ibot,itop , kk , icx ;       /* temp variables */
236 
237    new_prefix = (char *)calloc(100, sizeof(char));
238    filename = (char *)calloc(100, sizeof(char));
239    str = (char *)calloc(100, sizeof(char));
240    namestr = (char *)calloc(100, sizeof(char));
241    OutputFlag=0;
242    /*--------------------------------------------------------------------*/
243    /*----- Check inputs from AFNI to see if they are reasonable-ish -----*/
244 
245    /*--------- go to first input line ---------*/
246 
247    PLUTO_next_option(plint) ;
248 
249    idc      = PLUTO_get_idcode(plint) ;   /* get dataset item */
250    old_dset = PLUTO_find_dset(idc) ;      /* get ptr to dataset */
251    namestr  = DSET_PREFIX(old_dset) ;
252 
253 
254    if( old_dset == NULL )
255       return "*************************\n"
256              "Cannot find Input Dataset\n"
257              "*************************"  ;
258 
259    /*--------- go to second input line ---------*/
260 
261    PLUTO_next_option(plint) ;
262 
263   filename = PLUTO_get_string(plint) ;   /* get string item (the output prefix) */
264 
265   sprintf(new_prefix,"%s%s",filename,"_D3");
266 
267   if (strcmp(new_prefix,"_D3")==0){
268      OutputFlag=1;
269      sprintf(new_prefix,"%s%s",namestr,"_D3");
270   }
271 
272 
273    if (! PLUTO_prefix_ok(new_prefix) ){
274      PLUTO_popup_transient(plint,new_prefix);
275      return "*************************\n"
276              "Output filename already exists\n"
277              "*************************"  ;
278      }
279 
280 
281    PLUTO_popup_transient(plint,"Output file tags set automatically");
282 
283 
284    str  = PLUTO_get_string(plint) ;              /* get string item (the datum type) */
285    istr = PLUTO_string_index( str ,              /* find it in the list it came from */
286                               NUM_TYPE_STRINGS ,
287                               type_strings ) ;
288    switch( istr ){
289       default:
290       case 0:
291          new_datum = MRI_float ; break ;
292 	 break ;
293 
294       case 1: new_datum = MRI_byte  ; break ;  /* assign type of user's choice */
295       case 2: new_datum = MRI_short ; break ;
296       case 3: new_datum = DSET_BRICK_TYPE( old_dset , 0 ) ;  /* use old dataset type */
297    }
298 
299   /*--------- go to next input lines ---------*/
300 
301    PLUTO_next_option(plint) ;                 /* skip to next line */
302    ignore = PLUTO_get_number(plint) ;         /* get number item (ignore) */
303 
304 
305 
306 
307    ninp = DSET_NUM_TIMES(old_dset) ;   /* number of values in input */
308    nuse = ninp;              /* number of values to actually use */
309    nfreq=nuse;
310    nfft=nuse;
311 
312 
313    str  = PLUTO_get_string(plint) ;              /* get string item (the datum type) */
314    istr = PLUTO_string_index( str ,              /* find it in the list it came from */
315                               NUM_TYPE_STRINGSX ,
316                               type_stringsx ) ;
317    switch( istr ){
318       default:
319       case 0: image_type = 0; break;
320            }
321 
322   PLUTO_next_option(plint) ;                 /* skip to next line */
323   scale = PLUTO_get_number(plint) ;         /* get number item (scale) */
324 
325 
326    /*------------------------------------------------------*/
327    /*---------- At this point, the inputs are OK ----------*/
328 
329    PLUTO_popup_meter( plint ) ;  /* popup a progress meter */
330 
331    /*--------- set up pointers to each sub-brick in the input dataset ---------*/
332 
333    DSET_load( old_dset ) ;  /* must be in memory before we get pointers to it */
334 
335    old_datum = DSET_BRICK_TYPE( old_dset , 0 ) ; /* get old dataset datum type */
336 
337    switch( old_datum ){  /* pointer type depends on input datum type */
338 
339       default:
340          return "******************************\n"
341                 "Illegal datum in Input Dataset\n"
342                 "******************************"  ;
343 
344       /** create array of pointers into old dataset sub-bricks **/
345       /** Note that we skip the first 'ignore' sub-bricks here **/
346 
347       /*--------- input is bytes ----------*/
348       /* voxel #i at time #k is bptr[k][i] */
349       /* for i=0..nvox-1 and k=0..nuse-1.  */
350 
351       case MRI_byte:
352          bptr = (byte **) malloc( sizeof(byte *) * nuse ) ;
353          if( bptr == NULL ) return "Malloc\nFailure!\n [bptr]" ;
354          for( kk=0 ; kk < nuse ; kk++ )
355             bptr[kk] = (byte *) DSET_ARRAY(old_dset,kk) ;
356       break ;
357 
358       /*--------- input is shorts ---------*/
359       /* voxel #i at time #k is sptr[k][i] */
360       /* for i=0..nvox-1 and k=0..nuse-1.  */
361 
362       case MRI_short:
363          sptr = (short **) malloc( sizeof(short *) * nuse ) ;
364          if( sptr == NULL ) return "Malloc\nFailure!\n [sptr]" ;
365          for( kk=0 ; kk < nuse ; kk++ )
366             sptr[kk] = (short *) DSET_ARRAY(old_dset,kk) ;
367       break ;
368 
369       /*--------- input is floats ---------*/
370       /* voxel #i at time #k is fptr[k][i] */
371       /* for i=0..nvox-1 and k=0..nuse-1.  */
372 
373       case MRI_float:
374          fptr = (float **) malloc( sizeof(float *) * nuse ) ;
375          if( fptr == NULL ) return "Malloc\nFailure!\n [fptr]" ;
376          for( kk=0 ; kk < nuse ; kk++ )
377             fptr[kk] = (float *) DSET_ARRAY(old_dset,kk) ;
378       break ;
379 
380    } /* end of switch on input type */
381 
382    /*---- allocate space for 2 voxel timeseries and 1 FFT ----*/
383 
384 
385 
386    this = (float *)   malloc( sizeof(float) * nuse ) ;   /* input */
387    tarD3 = (float *) malloc( sizeof(float) * MAX(nuse,nfreq) ) ;
388    tarA3 = (float *) malloc( sizeof(float) * MAX(nuse,nfreq) ) ;
389    taravgD3 = (float *) malloc( sizeof(float) * MAX(nuse,nfreq) ) ;
390    /*flip = (float *)malloc( sizeof(float) * 1);*/
391    numAv = (float *)malloc( sizeof(float) * 1);
392 
393 
394   numT=nuse-ignore;
395 
396   if (OutputFlag==1)
397   sprintf(new_prefix,"%s%s",namestr,"_D3");
398   else
399   sprintf(new_prefix,"%s%s",filename,"_D3");
400 
401   new_dsetD3 = EDIT_empty_copy( old_dset );
402 
403   { char * his = PLUTO_commandstring(plint) ;
404   tross_Copy_History( old_dset , new_dsetD3 ) ;
405   tross_Append_History( new_dsetD3 , his ) ; free(his) ;
406   }
407 
408    	/*-- edit some of its internal parameters --*/
409 
410   ii = EDIT_dset_items(
411        new_dsetD3 ,
412          ADN_prefix      , new_prefix ,           /* filename prefix */
413          ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
414          ADN_datum_all   , new_datum ,            /* atomic datum */
415 	 ADN_nvals	      , numT ,
416 	 ADN_ntt	,numT,
417          ADN_none ) ;
418 
419 
420 
421   if (OutputFlag==1)
422   sprintf(new_prefix,"%s%s",namestr,"_A3");
423   else
424   sprintf(new_prefix,"%s%s",filename,"_A3");
425 
426   numT=nuse-ignore;
427   new_dsetA3 = EDIT_empty_copy( old_dset );
428 
429   { char * his = PLUTO_commandstring(plint) ;
430   tross_Copy_History( old_dset , new_dsetA3 ) ;
431   tross_Append_History( new_dsetA3 , his ) ; free(his) ;
432   }
433 
434    	/*-- edit some of its internal parameters --*/
435 
436   ii = EDIT_dset_items(
437        new_dsetA3 ,
438          ADN_prefix      , new_prefix ,           /* filename prefix */
439          ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
440          ADN_datum_all   , new_datum ,            /* atomic datum */
441 	 ADN_nvals	      , numT,
442 	 ADN_ntt	,numT,
443          ADN_none ) ;
444 
445 
446 
447   if (OutputFlag==1)
448   sprintf(new_prefix,"%s%s",namestr,"_avgD3");
449   else
450   sprintf(new_prefix,"%s%s",filename,"_avgD3");
451 
452   new_dsetavgD3 = EDIT_empty_copy( old_dset );
453 
454   { char * his = PLUTO_commandstring(plint) ;
455   tross_Copy_History( old_dset , new_dsetavgD3 ) ;
456   tross_Append_History( new_dsetavgD3 , his ) ; free(his) ;
457   }
458 
459    	/*-- edit some of its internal parameters --*/
460 
461   ii = EDIT_dset_items(
462         new_dsetavgD3 ,
463           ADN_prefix      , new_prefix ,           /* filename prefix */
464           ADN_malloc_type , DATABLOCK_MEM_MALLOC , /* store in memory */
465           ADN_datum_all   , new_datum ,            /* atomic datum */
466 	  ADN_nvals	      , 1,
467 	  ADN_ntt	,1,
468           ADN_none ) ;
469 
470 
471 
472 
473 
474    /*---------------------- make a new dataset ----------------------*/
475 
476 /*-------------------making a new dataset------------------------------------*/
477 
478 
479 
480 
481 
482    /*------ make floating point output sub-bricks
483             (only at the end will scale to byte or shorts)
484 
485             Output #ii at freq #kk will go into fout[kk][ii],
486             for kk=0..nfreq-1, and for ii=0..nvox-1.          ------*/
487 
488    nvox = old_dset->daxes->nxx * old_dset->daxes->nyy * old_dset->daxes->nzz ;
489 
490    foutD3 = (float **) malloc( sizeof(float *) * nuse ) ;  /* ptrs to sub-bricks */
491    foutA3 = (float **) malloc( sizeof(float *) * nuse ) ;  /* ptrs to sub-bricks */
492    foutavgD3 = (float **) malloc( sizeof(float *) * 1 ) ;  /* ptrs to sub-bricks */
493 
494 
495    if( foutD3 == NULL | foutA3 == NULL | foutavgD3 == NULL){
496       THD_delete_3dim_dataset( new_dsetD3 , False ) ;
497       THD_delete_3dim_dataset( new_dsetA3 , False ) ;
498       THD_delete_3dim_dataset( new_dsetavgD3 , False ) ;
499       FREE_WORKSPACE ;
500       return "Malloc\nFailure!\n [fout]" ;
501    }
502 
503    for( kk=0 ; kk < nfreq ; kk++ ){
504       foutD3[kk] = (float *) malloc( sizeof(float) * nvox ) ; /* sub-brick # kk */
505       foutA3[kk] = (float *) malloc( sizeof(float) * nvox ) ; /* sub-brick # kk */
506       foutavgD3[0] = (float *) malloc( sizeof(float) * nvox ) ; /* sub-brick # kk */
507       if( foutD3[kk] == NULL ) break ;
508       if( foutA3[kk] == NULL ) break ;
509       if( foutavgD3[0] == NULL ) break ;
510    }
511 
512    if( kk < nfreq ){
513       for( ; kk >= 0 ; kk-- ){
514        FREEUP(foutD3[kk]) ;
515        FREEUP(foutA3[kk]) ;
516        FREEUP(foutavgD3[0]) ;
517        }/* free all we did get */
518       THD_delete_3dim_dataset( new_dsetD3 , False ) ;
519       THD_delete_3dim_dataset( new_dsetA3 , False ) ;
520       THD_delete_3dim_dataset( new_dsetavgD3 , False ) ;
521       FREE_WORKSPACE ;
522       return "Malloc\nFailure!\n [arrays]" ;
523    }
524 
525    { char buf[128] ;
526      ii = (nfreq * nvox * sizeof(float)) / (1024*1024) ;
527      sprintf( buf , "  \n"
528                     "*** 3D+time ASL a3/d3:\n"
529                     "*** Using %d MBytes of workspace,\n "
530                     "*** with # time points = %d\n" , ii,numT ) ;
531      PLUTO_popup_transient( plint , buf ) ;
532    }
533 
534    /*----------------------------------------------------*/
535    /*----- Setup has ended.  Now do some real work. -----*/
536 
537    /***** loop over voxels *****/
538 
539 /* *(flip)=scale; */
540 
541 *(numAv)= nuse-ignore;
542 
543    for( ii=0 ; ii < nvox ; ii ++ ){  /* time series */
544 
545       switch( old_datum ){
546 
547 	case MRI_byte:
548             for( kk=0 ; kk < nuse ; kk++ ){
549             	this[kk] =  bptr[kk][ii] ;
550              }
551 
552          break ;
553 
554          case MRI_short:
555             for( kk=0 ; kk < nuse ; kk++ ){
556              this[kk] =  sptr[kk][ii] ;
557 
558             }
559          break ;
560 
561          case MRI_float:
562             for( kk=0 ; kk < nuse ; kk++ ){
563              this[kk] =  fptr[kk][ii] ;
564 
565             }
566 
567          break ;
568       }
569 
570       flip=scale*pow(-1,ignore+1);
571 
572       for( kk=0 ; kk < nuse-ignore ; kk++ ){
573 
574       		if (kk==nuse-1-ignore){
575         		*(*(foutD3+kk)+ii)=
576 			flip*( *(this+kk+ignore-1)-*(this+kk+ignore) );
577 
578 			*(*(foutA3+kk)+ii)=
579 			2*(*(this+kk+ignore-1)+*(this+kk+ignore));
580 
581 
582 			}
583 		else if (kk==0){
584 						/*D3 tag - control*/
585         		*(*(foutD3+kk)+ii)=
586 			flip*( *(this+kk+ignore)-*(this+kk+ignore+1) );
587 
588 			*(*(foutA3+kk)+ii)=
589 			2*(*(this+kk+ignore)+*(this+kk+ignore+1));
590 
591 			}
592 
593 		else{
594 			*(*(foutD3+kk)+ii)=
595 			flip*( 1*(*(this+kk+ignore-1))+-2*(*(this+kk+ignore))+1*(*(this+kk+ignore+1)) );
596 
597 			*(*(foutA3+kk)+ii)=
598 			((*(this+kk+ignore-1))+2*(*(this+kk+ignore))+(*(this+kk+ignore+1)));
599 
600 			flip=-1*flip;
601 
602 
603 			}
604 
605 
606 	}
607 
608 
609 
610       for( kk=0 ; kk < nuse-ignore ; kk++ )
611      *(*(foutavgD3)+ii)= *(*(foutavgD3)+ii)+(*(*(foutD3+kk)+ii));
612 
613      *(*(foutavgD3)+ii)=*(*(foutavgD3)+ii) / (*(numAv));
614 
615 
616       }
617 
618    DSET_unload( old_dset ) ;  /* don't need this no more */
619 
620    switch( new_datum ){
621 
622       /*** output is floats is the simplest:
623            we just have to attach the fout bricks to the dataset ***/
624 
625       case MRI_float:
626          for( kk=0 ; kk < nuse-ignore ; kk++ )
627             EDIT_substitute_brick( new_dsetD3 , kk , MRI_float , foutD3[kk] ) ;
628       break ;
629 
630       /*** output is shorts:
631            we have to create a scaled sub-brick from fout ***/
632 
633       case MRI_short:{
634          short * boutD3 ;
635          float facD3 ;
636 
637          for( kk=0 ; kk < nuse-ignore ; kk++ ){  /* loop over sub-bricks */
638 
639             /*-- get output sub-brick --*/
640 
641             boutD3 = (short *) malloc( sizeof(short) * nvox ) ;
642             if( boutD3 == NULL ){
643                fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
644                EXIT(1) ;
645             }
646 
647             /*-- find scaling and then scale --*/
648 
649             facD3  = MCW_vol_amax( nvox,1,1 , MRI_float , foutD3[kk] ) ;
650             if( facD3  > 0.0 ){
651                facD3  = 32767.0 / facD3  ;
652                EDIT_coerce_scale_type( nvox,facD3  ,
653                                        MRI_float,foutD3[kk] , MRI_short,boutD3  ) ;
654                facD3  = 1.0 / facD3  ;
655             }
656 
657             free( foutD3[kk] ) ;  /* don't need this anymore */
658 
659             /*-- put output brick into dataset, and store scale factor --*/
660 
661             EDIT_substitute_brick( new_dsetD3 , kk , MRI_short , boutD3  ) ;
662             tarD3 [kk] = facD3  ;
663 
664 
665          }
666 
667          /*-- save scale factor array into dataset --*/
668 
669          EDIT_dset_items( new_dsetD3 , ADN_brick_fac , tarD3  , ADN_none ) ;
670 
671       }
672       break ;
673 
674       /*** output is bytes (byte = unsigned char)
675            we have to create a scaled sub-brick from fout ***/
676 
677       case MRI_byte:{
678          byte * boutD3  ;
679          float facD3  ;
680 
681          for( kk=0 ; kk < nuse-ignore ; kk++ ){  /* loop over sub-bricks */
682 
683             /*-- get output sub-brick --*/
684 
685             boutD3  = (byte *) malloc( sizeof(byte) * nvox ) ;
686             if( boutD3  == NULL ){
687                fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
688                EXIT(1) ;
689             }
690 
691             /*-- find scaling and then scale --*/
692 
693             facD3  = MCW_vol_amax( nvox,1,1 , MRI_float , foutD3[kk] ) ;
694             if( facD3  > 0.0 ){
695                facD3  = 255.0 / facD3  ;
696                EDIT_coerce_scale_type( nvox,facD3  ,
697                                        MRI_float,foutD3[kk] , MRI_byte,boutD3 ) ;
698                facD3 = 1.0 / facD3  ;
699             }
700 
701             free( foutD3[kk] ) ;  /* don't need this anymore */
702 
703             /*-- put output brick into dataset, and store scale factor --*/
704 
705             EDIT_substitute_brick( new_dsetD3 , kk , MRI_byte , boutD3  ) ;
706             tarD3 [kk] = facD3  ;
707 
708 
709          }
710 
711          /*-- save scale factor array into dataset --*/
712 
713          EDIT_dset_items( new_dsetD3 , ADN_brick_fac , tarD3  , ADN_none ) ;
714       }
715       break ;
716 
717    } /* end of switch on output data type */
718 
719 
720    switch( new_datum ){
721 
722       /*** output is floats is the simplest:
723            we just have to attach the fout bricks to the dataset ***/
724 
725       case MRI_float:
726          for( kk=0 ; kk < nuse-ignore ; kk++ )
727             EDIT_substitute_brick( new_dsetA3 , kk , MRI_float , foutA3[kk] ) ;
728       break ;
729 
730       /*** output is shorts:
731            we have to create a scaled sub-brick from fout ***/
732 
733       case MRI_short:{
734          short * boutA3 ;
735          float facA3 ;
736 
737          for( kk=0 ; kk < nuse-ignore ; kk++ ){  /* loop over sub-bricks */
738 
739             /*-- get output sub-brick --*/
740 
741             boutA3 = (short *) malloc( sizeof(short) * nvox ) ;
742             if( boutA3 == NULL ){
743                fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
744                EXIT(1) ;
745             }
746 
747             /*-- find scaling and then scale --*/
748 
749             facA3 = MCW_vol_amax( nvox,1,1 , MRI_float , foutA3[kk] ) ;
750             if( facA3 > 0.0 ){
751                facA3 = 32767.0 / facA3 ;
752                EDIT_coerce_scale_type( nvox,facA3 ,
753                                        MRI_float,foutA3[kk] , MRI_short,boutA3 ) ;
754                facA3 = 1.0 / facA3 ;
755             }
756 
757             free( foutA3[kk] ) ;  /* don't need this anymore */
758 
759             /*-- put output brick into dataset, and store scale factor --*/
760 
761             EDIT_substitute_brick( new_dsetA3 , kk , MRI_short , boutA3 ) ;
762             tarA3[kk] = facA3 ;
763 
764 
765          }
766 
767          /*-- save scale factor array into dataset --*/
768 
769          EDIT_dset_items( new_dsetA3 , ADN_brick_fac , tarA3 , ADN_none ) ;
770 
771       }
772       break ;
773 
774       /*** output is bytes (byte = unsigned char)
775            we have to create a scaled sub-brick from fout ***/
776 
777       case MRI_byte:{
778          byte * boutA3 ;
779          float facA3 ;
780 
781          for( kk=0 ; kk < nuse-ignore ; kk++ ){  /* loop over sub-bricks */
782 
783             /*-- get output sub-brick --*/
784 
785             boutA3 = (byte *) malloc( sizeof(byte) * nvox ) ;
786             if( boutA3 == NULL ){
787                fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
788                EXIT(1) ;
789             }
790 
791             /*-- find scaling and then scale --*/
792 
793             facA3 = MCW_vol_amax( nvox,1,1 , MRI_float , foutA3[kk] ) ;
794             if( facA3 > 0.0 ){
795                facA3 = 255.0 / facA3 ;
796                EDIT_coerce_scale_type( nvox,facA3 ,
797                                        MRI_float,foutA3[kk] , MRI_byte,boutA3 ) ;
798                facA3 = 1.0 / facA3 ;
799             }
800 
801             free( foutA3[kk] ) ;  /* don't need this anymore */
802 
803             /*-- put output brick into dataset, and store scale factor --*/
804 
805             EDIT_substitute_brick( new_dsetA3 , kk , MRI_byte , boutA3 ) ;
806             tarA3[kk]= facA3 ;
807 
808 
809          }
810 
811          /*-- save scale factor array into dataset --*/
812 
813          EDIT_dset_items( new_dsetA3 , ADN_brick_fac , tarA3 , ADN_none ) ;
814       }
815       break ;
816 
817    } /* end of switch on output data type */
818 
819 
820      switch( new_datum ){
821 
822       case MRI_float:{
823 
824             EDIT_substitute_brick( new_dsetavgD3 , 0 , MRI_float , foutavgD3[0] ) ;
825 
826 
827     }
828       break ;
829 
830       case MRI_short:{
831          short * boutavgD3 ;
832          float facavgD3 ;
833 
834             boutavgD3 = (short *) malloc( sizeof(short) * nvox ) ;
835             if( boutavgD3 == NULL ){
836                fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
837                EXIT(1) ;
838             }
839 
840             facavgD3 = MCW_vol_amax( nvox,1,1 , MRI_float , foutavgD3[0] ) ;
841             if( facavgD3 > 0.0 ){
842                facavgD3 = 32767.0 / facavgD3 ;
843                EDIT_coerce_scale_type( nvox,facavgD3 ,
844                                        MRI_float,foutavgD3[0] , MRI_short,boutavgD3 ) ;
845                facavgD3 = 1.0 / facavgD3 ;
846             }
847 
848 
849 
850             EDIT_substitute_brick( new_dsetavgD3 , 0 , MRI_short , boutavgD3 ) ;
851             taravgD3[0] = facavgD3 ;
852 
853              EDIT_dset_items( new_dsetavgD3 , ADN_brick_fac , taravgD3 , ADN_none ) ;
854 
855 
856 
857       }
858       break ;
859 
860       case MRI_byte:{
861          byte * boutavgD3 ;
862          float facavgD3 ;
863 
864 
865             boutavgD3 = (byte *) malloc( sizeof(byte) * nvox ) ;
866             if( boutavgD3 == NULL ){
867                fprintf(stderr,"\nFinal malloc error in plug_power!\n\a") ;
868                EXIT(1) ;
869             }
870 
871             facavgD3 = MCW_vol_amax( nvox,1,1 , MRI_float , foutavgD3[0] ) ;
872             if( facavgD3 > 0.0 ){
873                facavgD3 = 255.0 / facavgD3 ;
874                EDIT_coerce_scale_type( nvox,facavgD3 ,
875                                        MRI_float,foutavgD3[0] , MRI_byte,boutavgD3 ) ;
876                facavgD3 = 1.0 / facavgD3 ;
877             }
878 
879 
880 
881             EDIT_substitute_brick( new_dsetavgD3 , 0 , MRI_byte , boutavgD3 ) ;
882             taravgD3[0]= facavgD3 ;
883 
884             EDIT_dset_items( new_dsetavgD3 , ADN_brick_fac , taravgD3 , ADN_none ) ;
885 
886 
887 
888 
889       }
890       break ;
891 
892    } /* endasda of switch on output data type */
893 
894 
895 
896 
897    /*-------------- Cleanup and go home ----------------*/
898 
899 
900 
901    PLUTO_add_dset( plint , new_dsetD3 , DSET_ACTION_NONE ) ;
902   PLUTO_add_dset( plint , new_dsetA3 , DSET_ACTION_NONE ) ;
903   PLUTO_add_dset( plint , new_dsetavgD3 , DSET_ACTION_NONE ) ;
904 
905 
906 
907    FREE_WORKSPACE ;
908    free(numAv);
909 
910 
911    return NULL ;  /* null string returned means all was OK */
912 }
913 
914 #ifdef ALLOW_TESTING
915 /*****************************************************************************
916  -----------------------------------------------------------------------------
917            Create the second interface within this plugin.
918  -----------------------------------------------------------------------------*/
919 
TEST_init(void)920 PLUGIN_interface * TEST_init( void )
921 {
922    PLUGIN_interface * plint ;     /* will be the output of this routine */
923 
924    /*---------------- set titles and call point ----------------*/
925 
926    plint = PLUTO_new_interface( "Testing" ,
927                                 "Testing, Testing, 1-2-3 ..." ,
928                                 NULL ,
929                                 PLUGIN_CALL_VIA_MENU , TEST_main  ) ;
930 
931    PLUTO_add_hint( plint , "1-2-3, 1-2-3, ..." ) ;
932 
933    PLUTO_add_option( plint ,
934                      "Input" ,  /* label at left of input line */
935                      "Input" ,  /* tag to return to plugin */
936                      TRUE       /* is this mandatory? */
937                    ) ;
938 
939    PLUTO_add_dataset_list(  plint ,
940                             "Datasets" ,       /* label next to button   */
941                             ANAT_ALL_MASK ,    /* take any anat datasets */
942                             FUNC_FIM_MASK ,    /* only allow fim funcs   */
943                             DIMEN_4D_MASK |    /* need 3D+time datasets  */
944                             BRICK_ALLREAL_MASK /* need real-valued datasets */
945                          ) ;
946    return plint ;
947 }
948 
TEST_main(PLUGIN_interface * plint)949 char * TEST_main( PLUGIN_interface * plint )
950 {
951    MRI_IMAGE * tsim ;
952    MCW_idclist * idclist ;
953    MCW_idcode * idc ;
954    THD_3dim_dataset * dset ;
955    char str[256] ;
956    int id ;
957 
958    /*--------- go to first input line ---------*/
959 
960    PLUTO_next_option(plint) ;
961 
962    idclist = PLUTO_get_idclist(plint) ;
963    if( PLUTO_idclist_count(idclist) == 0 )
964       return " \nNo input dataset list!\n " ;
965 
966    id = 0 ;
967    do {
968       idc  = PLUTO_idclist_next(idclist) ;
969       dset = PLUTO_find_dset(idc) ;
970       if( dset == NULL ) return NULL ;
971       id++ ;
972       sprintf(str, " \nDataset %d = %s\n nx = %d\n ny = %d\n nz = %d\n " ,
973               id , DSET_FILECODE(dset) , dset->daxes->nxx,dset->daxes->nyy,dset->daxes->nzz ) ;
974 
975       PLUTO_popup_transient( plint , str ) ;
976    } while(1) ;
977    return NULL ;
978 }
979 #endif  /* ALLOW_TESTING */
980