00001 /*plug_permtest.c - AFNI plugin that applies a permutation test to a 3D+time
00002   dataset to create a fizt dataset.
00003   Copyright (c) 2000 - 2005 Matthew Belmonte
00004 
00005   This program is free software; you can redistribute it and/or
00006   modify it under the terms of the GNU General Public License
00007   as published by the Free Software Foundation; either version 2
00008   of the License, or (at your option) any later version.
00009 
00010   This program is distributed in the hope that it will be useful,
00011   but WITHOUT ANY WARRANTY; without even the implied warranty of
00012   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013   GNU General Public License for more details.
00014 
00015   You should have received a copy of the GNU General Public License
00016   along with this program; if not, write to the Free Software
00017   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
00018 
00019   If you find this program useful, please send mail to Matthew Belmonte
00020   <belmonte@mit.edu>.
00021 */
00022 
00023 #define NUM_ITERS 10000
00024 #define NUM_COORDS 2
00025 
00026 #define BLAST 33333.0
00027 #include <sys/types.h>
00028 #include <stdio.h>
00029 #include <stdlib.h>
00030 #include <math.h>
00031 #include <string.h>
00032 #include "afni.h"
00033 
00034 /*a macro that wraps THD_pval_to_stat() (which in turn evaluates to a call to
00035   the one-tailed probability-to-z-score function normal_p2t()) so that it
00036   behaves like Gary Perlman's two-tailed zcrit() function in Netlib*/
00037 #define critz(p) \
00038    (((p) >= 0.5)? THD_pval_to_stat(2.0*(1.0-(p)), FUNC_ZT_TYPE, (float *)0) \
00039                 : -THD_pval_to_stat(2.0*(p), FUNC_ZT_TYPE, (float *)0))
00040 
00041 static char help[] =
00042   "This plugin implements a permutation test, a nonparametric statistical\n"
00043   "method that avoids the pitfall of over-correcting for multiple comparisons\n"
00044   "since it implicitly takes into account spatial correlations in the data.\n"
00045   "If you use this software, please take a moment to send mail to the author,\n"
00046   "belmonte@mit.edu, and cite the following paper in your report:\n\n"
00047 
00048   "Matthew K Belmonte and Deborah A Yurgelun-Todd, `Permutation Testing Made\n"
00049   "Practical for Functional Magnetic Resonance Image Analysis',\n"
00050   "IEEE Transactions on Medical Imaging 20(3):243-248 (March 2001).\n\n"
00051 
00052 "USAGE\n\n"
00053 
00054   "`Input' is the echo-planar dataset for which the permutation test is to be\n"
00055   "computed.\n\n"
00056 
00057   "`Ideal' is the ideal time series against which the echo-planar data are to\n"
00058   "be correlated.\n\n"
00059 
00060   "`Ort', an optional parameter, is a time series to which the echo-planar\n"
00061   "data are to be orthogonalised.  This feature is useful when one wants to\n"
00062   "analyse responses to two superposed conditions independently.  (Note that\n"
00063   "the data are always detrended, i.e. orthogonalised to a linear series --\n"
00064   "this optional orthogonalisation is in addition to the detrending step.)\n\n"
00065 
00066   "`Output' is the output dataset.  The activation data saved in this dataset\n"
00067   "consist of a sub-brick of regression coefficients and a sub-brick of\n"
00068   "probability values.  As AFNI provides no standard format for storing\n"
00069   "probabilities, these are mapped to z-scores and the entire dataset is\n"
00070   "saved as a FIZT dataset.\n\n"
00071 
00072   "`alpha level' is the tail probability at which to stop marking voxels as\n"
00073   "activated.  (Note that the slider in the `Define Function' panel will\n"
00074   "have no effect below this value.)  You should avoid setting this value\n"
00075   "very high, lest the Phase 3 algorithm run out of substitute points in too\n"
00076   "many cases (for an explanation of substitutes, see the article cited above).\n\n"
00077 
00078   "`two-tailed', `one-tailed positive', and `one-tailed negative' select the\n"
00079   "tail(s) of interest.  Exactly one of these options must be selected.\n\n"
00080 
00081   "`Mask', an optional parameter, is a FIM dataset whose three spatial\n"
00082   "dimensions match those of the echo-planar dataset.  Only those\n"
00083   "coordinates whose mask values are between 1 and the largest positive\n"
00084   "short integer are included in the permutation test.  These inclusion\n"
00085   "limits can be overridden using the optional `least mask value' and\n"
00086   "`greatest mask value' parameters.  Mask datasets can be generated\n"
00087   "automatically using the Threshold plugin (q.v.), or using 3dAutomask\n\n"
00088 
00089 "AUTHOR\n\n"
00090 
00091   "This plugin was written by Matthew Belmonte <belmonte@mit.edu> of the\n"
00092   "MIT Student Information Processing Board, supported by a grant from the\n"
00093   "National Alliance for Autism Research.\n\n"
00094 
00095 "REVISION HISTORY\n\n"
00096 
00097   "1.3  20 January 2005   added support for 3dAutomask and other byte-valued masks\n"
00098   "1.2  19 December 2002  fixed a dtree bug that caused a rare crash in phase 3\n"
00099   "1.1  14 June 2001      cosmetic changes only\n"
00100   "1.0  4 January 2001    initial release\n\n"
00101 
00102 "SEE ALSO\n\n"
00103 
00104   "Threshold Plugin  or  3dAutomask program\n"
00105   "Draw Dataset Plugin\n"
00106   "(These normally are applied before the permutation test.)",
00107 
00108             hint[] = "compute FIM with permutation test",
00109             input_label[] = "Input",
00110             ts_label[] = "Ideal",
00111             ort_label[] = "Ort",
00112             output_label[] = "Output",
00113             alpha_label[] = "alpha level (0,1]",
00114             tails2[] = "two-tailed",
00115             tails1pos[] = "one-tailed positive",
00116             tails1neg[] = "one-tailed negative",
00117             mask_label[] = "Mask",
00118             masklo_label[] = "least mask value",
00119             maskhi_label[] = "greatest mask value";
00120 
00121 /*data structures for trees*/
00122 
00123 /*singly-indexed tree: ordered on correlations*/
00124 typedef struct _snode {
00125         double corr;                    /*correlation (key)*/
00126         int coords;                     /*coordinates: x+xdim*(y+ydim*z) */
00127         struct _snode *left, *right;    /*children*/
00128         } SNODE;
00129 
00130 typedef struct {
00131         SNODE *mem,                     /*base address of tree node storage*/
00132               *next,                    /*next free tree storage address*/
00133               *root;                    /*==tree_mem if tree is non-null*/
00134         int xdim, ydim;                 /*maximum coordinate values*/
00135         } STREE;
00136 
00137 #define stree_null(t) ((t)->root == (SNODE *)0)
00138 
00139 typedef struct {                        /*pairs (r, (x,y,z))*/
00140         double corr;
00141         int coords;
00142         } CLIST;
00143 
00144 /*doubly-indexed tree: ordered on correlations and on coordinates*/
00145 typedef struct _dnode {
00146         CLIST *data;                    /*correlations and coordinates*/
00147         int dlen;                       /*# of valid entries in data[]*/
00148         int size;                       /*size of this correlation subtree*/
00149         struct _dnode *corr_lchild,     /*children in correlation tree*/
00150                       *corr_rchild,
00151                       *corr_parent,     /*parent in correlation tree*/
00152                       *coord_lchild,    /*children in coordinate tree*/
00153                       *coord_rchild;
00154         } DNODE;
00155 
00156 typedef struct {
00157         DNODE *mem,                     /*base address of tree node storage*/
00158               *next,                    /*next free tree storage address*/
00159               *corr_root, *coord_root;  /*==mem if tree is non-null*/
00160         int xdim, ydim;                 /*maximum coordinate values*/
00161         } DTREE;
00162 
00163 #define dtree_size(t) (((t)->corr_root == (DNODE *)0)? 0: ((t)->corr_root->size))
00164 
00165 /*encoding and decoding functions for mapping a triple to a single integer*/
00166 #define coord_encode(x, y, z, xdim, ydim) ((x) + (xdim)*((y) + (ydim)*(z)))
00167 
00168 static void coord_decode(int coords, int *x, int *y, int *z, int xdim, int ydim)
00169   {
00170   *x = coords % xdim;
00171   coords /= xdim;
00172   *y = coords % ydim;
00173   *z = coords / ydim;
00174   }
00175 
00176 #define sign(x) (((x)<0)? -1: 1)
00177 
00178 #ifdef PERMTEST_DEBUG
00179 
00180 /*flush output buffers on SIGUSR1 - useful for batch debugging*/
00181 #include <signal.h>
00182 static void flush(int sig)
00183   {
00184   fflush(stdout);
00185   fflush(stderr);
00186   signal(SIGUSR1, flush);
00187   }
00188 
00189 /*slave routine to stree_traverse()*/
00190 void rcsv_stree_traverse(SNODE *t, int indent)
00191   {
00192   register int i;
00193   while(t)
00194     {
00195     for(i = 0; i != indent; i++)
00196       putchar(' ');
00197     indent += 2;
00198     printf("%lx: corr=%f coords=%d left=%lx right=%lx\n", (long)t, t->corr, t->coords, (long)(t->left), (long)(t->right));
00199     rcsv_stree_traverse(t->left, indent);
00200     t = t->right;
00201     }
00202   }
00203 
00204 /*Print a traversal of the given stree.*/
00205 void stree_traverse(STREE *t)
00206   {
00207   printf("t=%lx\n  mem=%lx  next=%lx = mem+%ld  root=%lx xdim=%d ydim=%d\n",
00208         (long)t,
00209         (long)(t->mem),
00210         (long)(t->next),
00211         (((long)(t->next))-(long)(t->mem))/sizeof(SNODE),
00212         (long)(t->root),
00213         t->xdim, t->ydim);
00214   rcsv_stree_traverse(t->root, 0);
00215   putchar('\n');
00216   }
00217 
00218 /*slave routine to dtree_traverse()*/
00219 void rcsv_dtree_traverse(DNODE *t, int indent, int order)
00220   {
00221   register int i;
00222   while(t)
00223     {
00224     for(i = 0; i != indent; i++)
00225       putchar(' ');
00226     indent += 2;
00227     printf("%lx: data=", (long)t);
00228     for(i = 0; i != t->dlen; i++)
00229       printf("%f %d ", t->data[i].corr, t->data[i].coords);
00230     printf("size=%d corr_lchild=%lx corr_rchild=%lx corr_parent=%lx coord_lchild=%lx coord_rchild=%lx\n", t->size, (long)(t->corr_lchild), (long)(t->corr_rchild), (long)(t->corr_parent), (long)(t->coord_lchild), (long)(t->coord_rchild));
00231     if(order)
00232       {
00233       rcsv_dtree_traverse(t->coord_lchild, indent, order);
00234       t = t->coord_rchild;
00235       }
00236     else
00237       {
00238       rcsv_dtree_traverse(t->corr_lchild, indent);
00239       t = t->corr_rchild;
00240       }
00241     }
00242   }
00243 
00244 /*Print a traversal of a dtree.  Use corr links if order=0, coord links if
00245   order=1.*/
00246 void dtree_traverse(DTREE *t, int order)
00247   {
00248   printf("t=%lx\n  mem=%lx  next=%lx = mem+%ld  coord_root=%lx corr_root=%lx xdim=%d ydim=%d\n",
00249         (long)t,
00250         (long)(t->mem),
00251         (long)(t->next),
00252         (((long)(t->next))-(long)(t->mem))/sizeof(DNODE),
00253         (long)(t->coord_root),
00254         (long)(t->corr_root),
00255         t->xdim, t->ydim);
00256   rcsv_dtree_traverse((order? t->coord_root: t->corr_root), 0, order);
00257   putchar('\n');
00258   }
00259 
00260 /*slave routine to dtree_check_sizes()*/
00261 int rcsv_dtree_check_sizes(DNODE *t)
00262   {
00263   int s;
00264   if(t)
00265     {
00266     s = 1+rcsv_dtree_check_sizes(t->corr_lchild)+rcsv_dtree_check_sizes(t->corr_rchild);
00267     if(s != t->size)
00268       {
00269       printf("%lx->size = %d should be %d\n", (long)t, t->size, s);
00270       fflush(stdout);
00271       }
00272     return(t->size);
00273     }
00274   return 0;
00275   }
00276 
00277 /*Check that the size fields of a dtree are consistent with each other and with
00278   the structure of the tree.*/
00279 int dtree_check_sizes(DTREE *t)
00280   {
00281   return(rcsv_dtree_check_sizes(t->corr_root));
00282   }
00283 
00284 /*slave routine to dtree_circtest()*/
00285 void rcsv_dtree_circtest(DNODE *t, DNODE *mem, char *mark, int n, int order)
00286   {
00287   register int m;
00288   while(t)
00289     {
00290     m = (int)((((long)t)-(long)mem)/(((long)(mem+1))-((long)mem)));
00291     if((m >= n) || (m < 0))
00292       printf("out-of-bounds reference %lx (mem+%d)\n", (long)t, m);
00293     else if(mark[m])
00294       printf("duplicate reference %lx (mem+%d)\n", (long)t, m);
00295     else
00296       {
00297       mark[m] = 1;
00298       if(order)
00299         {
00300         rcsv_dtree_circtest(t->coord_lchild, mem, mark, n, order);
00301         t = t->coord_rchild;
00302         }
00303       else
00304         {
00305         rcsv_dtree_circtest(t->corr_lchild, mem, mark, n, order);
00306         t = t->corr_rchild;
00307         }
00308       }
00309     }
00310   }
00311 
00312 /*Test the given dtree for circular or out-of-bounds references.  Use corr links
00313   if order=0, coord links if order=1.*/
00314 void dtree_circtest(DTREE *t, int order)
00315   {
00316   char *mark;
00317   int n;
00318   n = (int)((((long)(t->next))-(long)(t->mem))/(((long)(t->mem+1))-((long)(t->mem))));
00319   mark = calloc(n, 1);
00320   if(mark)
00321     {
00322     rcsv_dtree_circtest((order? t->coord_root: t->corr_root), t->mem, mark, n, order);
00323     free(mark);
00324     }
00325   else
00326     fprintf(stderr, "dtree_circtest: couldn't calloc mark\n");
00327   }
00328 
00329 /*slave routine to dtree_check_parent_pointers()*/
00330 int rcsv_dtree_check_parent_pointers(DNODE *t)
00331   {
00332   int bad;
00333   bad = 0;
00334   while(t != (DNODE *)0)
00335     {
00336     if(t->corr_lchild)
00337       {
00338       if(t->corr_lchild->corr_parent != t)
00339         {
00340         printf("%lx->corr_lchild = %lx, but %lx->corr_parent = %lx\n",
00341           (long)t, (long)(t->corr_lchild), (long)(t->corr_lchild),
00342           (long)(t->corr_lchild->corr_parent));
00343         bad = 1;
00344         }
00345       if(rcsv_dtree_check_parent_pointers(t->corr_lchild)) bad = 1;
00346       }
00347     if((t->corr_rchild) && (t->corr_rchild->corr_parent != t))
00348       {
00349       printf("%lx->corr_rchild = %lx, but %lx->corr_parent = %lx\n",
00350         (long)t, (long)(t->corr_rchild), (long)(t->corr_rchild),
00351         (long)(t->corr_rchild->corr_parent));
00352       bad = 1;
00353       }
00354     t = t->corr_rchild;
00355     }
00356   return bad;
00357   }
00358 
00359 /*Verify that the corr_parent fields inside t are consistent with the
00360   corr_lchild and corr_rchild fields.*/
00361 int dtree_check_parent_pointers(DTREE *t)
00362   {
00363   return(rcsv_dtree_check_parent_pointers(t->corr_root));
00364   }
00365 
00366 /*Print hex dump of an IEEE double-precision floating-point number, MSB first.*/
00367 void fdump(double *x)
00368   {
00369   double test = 2.0;
00370   int i;
00371   register char *c;
00372   c = (char *)x;
00373   if(*(char *)&test)
00374     for(i = 7; i >= 0; i--)
00375       printf("%02x", c[i]&0xff);
00376   else
00377     for(i = 0; i != 8; i++)
00378       printf("%02x", c[i]&0xff);
00379   }
00380 
00381 #endif
00382 
00383 /*Which metric is used for comparison of correlation values depends on which
00384   tail(s) of the distribution interest us.  For the two-tailed case, our
00385   `maximum' should be the most extreme value from either tail, so we use the
00386   absolute value function fabs() from libm.  For the negative one-tailed case,
00387   our `maximum' should be the most negative value, so we use the negation
00388   function neg().  For the positive one-tailed case, our maximum should be the
00389   most positive value, so we use the identity function id().*/
00390 
00391 static double id(double r)
00392   {
00393   return(r);
00394   }
00395 
00396 static double neg(double r)
00397   {
00398   return(-r);
00399   }
00400 
00401 static double (*magnitude)(double);
00402 
00403 /*Allocate storage and initialise a tree that can contain up to nelem nodes.*/
00404 static STREE *stree_create(int xdim, int ydim, int nelem)
00405   {
00406   STREE *t;
00407   t = (STREE *)malloc(sizeof(*t));
00408   if(t == (STREE *)0)
00409     return((STREE *)0);
00410   t->mem = (SNODE *)calloc(nelem, sizeof(*(t->mem)));
00411   if(t->mem == (SNODE *)0)
00412     {
00413     free(t);
00414     return((STREE *)0);
00415     }
00416   t->next = t->mem;
00417   t->root = (SNODE *)0;
00418   t->xdim = xdim;
00419   t->ydim = ydim;
00420   return(t);
00421   }
00422 
00423 /*Deallocate a tree.*/
00424 static void stree_destroy(STREE *t)
00425   {
00426   free(t->mem);
00427   free(t);
00428   }
00429 
00430 /*Insert the given value into the tree.*/
00431 static void stree_insert(STREE *t, double r, int x, int y, int z)
00432   {
00433   register SNODE **p;
00434   double rmag;
00435   p = &(t->root);
00436   rmag = (*magnitude)(r);
00437 /*inv: The path from t->root to p is labelled with values whose magnitudes
00438   bound or equal that of r.*/
00439   while(*p != (SNODE *)0)
00440     p = (rmag < (*magnitude)((*p)->corr))? &((*p)->left): &((*p)->right);
00441   t->next->corr = r;
00442   t->next->coords = coord_encode(x, y, z, t->xdim, t->ydim);
00443   t->next->left = (SNODE *)0;
00444   t->next->right = (SNODE *)0;
00445   *p = t->next++;
00446   }
00447 
00448 /*Delete the maximum correlation value from the given tree, and return that
00449   value and its associated coordinates.*/
00450 static double stree_extract_max(STREE *t, int *x, int *y, int *z)
00451   {
00452   register SNODE **p, **pp;
00453   SNODE *target;
00454   if(t->root == (SNODE *)0)
00455     {
00456     fprintf(stderr, "Error: STREE underflow\n");
00457     return(0.0);
00458     }
00459   if((t->root->left == (SNODE *)0) && (t->root->right == (SNODE *)0))
00460     {
00461     target = t->root;
00462     t->root = (SNODE *)0;
00463     t->next = t->mem;
00464     }
00465   else
00466     {
00467     p = &(t->root);
00468   /*inv: **pp and **p are successive links in a chain of right children
00469     beginning at t->root.*/
00470     do
00471       {
00472       pp = p;
00473       p = &((*p)->right);
00474       } while((*p) != (SNODE *)0);
00475     target = *pp;
00476     *pp = (*pp)->left;
00477     }
00478   coord_decode(target->coords, x, y, z, t->xdim, t->ydim);
00479   return(target->corr);
00480   }
00481 
00482 /*Allocate storage and initialise a dtree that can contain up to nelem nodes.*/
00483 static DTREE *dtree_create(int xdim, int ydim, int nelem)
00484   {
00485   DTREE *t;
00486   t = (DTREE *)malloc(sizeof(*t));
00487   if(t == (DTREE *)0)
00488     return((DTREE *)0);
00489   t->mem = (DNODE *)calloc(nelem, sizeof(*(t->mem)));
00490   if(t->mem == (DNODE *)0)
00491     {
00492     free(t);
00493     return((DTREE *)0);
00494     }
00495   t->next = t->mem;
00496   t->corr_root = (DNODE *)0;
00497   t->coord_root = (DNODE *)0;
00498   t->xdim = xdim;
00499   t->ydim = ydim;
00500   return(t);
00501   }
00502 
00503 /*Deallocate a dtree.*/
00504 static void dtree_destroy(DTREE *t)
00505   {
00506   free(t->mem);
00507   free(t);
00508   }
00509 
00510 /*Insert the given value into the dtree, using the given storage location.*/
00511 static void dtree_insert_at_node(DTREE *t, DNODE *node)
00512   {
00513   register DNODE **p,                   /*pointer into the coord tree*/
00514                  *q, *qparent;          /*pointers to adjacent corr tree nodes*/
00515   p = &(t->coord_root);
00516 /*inv: The path from t->coord_root to *p is labelled with values that bound or
00517   equal node->data->coords.*/
00518   while(*p != (DNODE *)0)
00519     p = ((node->data->coords > (*p)->data->coords)? &((*p)->coord_rchild): &((*p)->coord_lchild));
00520 /*link the node into the coord tree*/
00521   node->coord_lchild = (DNODE *)0;
00522   node->coord_rchild = (DNODE *)0;
00523   *p = node;
00524 /*link the node into the corr tree*/
00525   node->size = 1;
00526   node->corr_lchild = (DNODE *)0;
00527   node->corr_rchild = (DNODE *)0;
00528   if(t->corr_root == (DNODE *)0)
00529     {
00530   /*insertion into an empty tree is a special case*/
00531     node->corr_parent = (DNODE *)0;
00532     t->corr_root = node;
00533     }
00534   else
00535     {
00536     q = t->corr_root;
00537   /*inv: The path from t->corr_root to q is labelled with values that bound or
00538     equal r.*/
00539     do
00540       {
00541       q->size++;
00542       qparent = q;
00543       q = ((node->data->corr > q->data->corr)? q->corr_rchild: q->corr_lchild);
00544       } while(q != (DNODE *)0);
00545   /*link the node into the correlation tree*/
00546     node->corr_parent = qparent;
00547     if(node->data->corr > qparent->data->corr)
00548       qparent->corr_rchild = node;
00549     else
00550       qparent->corr_lchild = node;
00551     }
00552   }
00553 
00554 static int num_coords;
00555 
00556 /*Insert the given value into the dtree.*/
00557 static void dtree_insert(DTREE *t, CLIST *clist)
00558   {
00559   t->next->data = clist;
00560   t->next->dlen = num_coords;
00561   dtree_insert_at_node(t, t->next++);
00562   }
00563 
00564 static int num_coords_exhausted;
00565 
00566 /*p points to a coord_lchild, coord_rchild, or coord_root field that is to be
00567   unlinked from its coord children and its corr children and corr parent within
00568   the tree t.*/
00569 void dtree_unlink_node(DTREE *t, DNODE **p)
00570   {
00571   register DNODE **q;
00572   DNODE *node, *parent, *temp;
00573   node = *p;
00574   if(node->coord_lchild == (DNODE *)0)
00575     *p = node->coord_rchild;
00576   else if(node->coord_rchild == (DNODE *)0)
00577     *p = node->coord_lchild;
00578   else
00579     {
00580   /*inv: *q is the rightmost node at the current level in the left subtree of
00581     (*p).*/
00582     for(q = &(node->coord_lchild); (*q)->coord_rchild != (DNODE *)0; q = &((*q)->coord_rchild))
00583       ;
00584     (*q)->coord_rchild = node->coord_rchild;
00585     *p = *q;
00586     if(q != &(node->coord_lchild))
00587       {
00588       DNODE *temp;
00589       temp = (*q)->coord_lchild;
00590       (*q)->coord_lchild = node->coord_lchild;
00591       *q = temp;
00592       }
00593     }
00594 /*The node has been unlinked from the coord tree.  Now unlink it from the corr
00595   tree, using a procedure analogous to that which was implemented above for the
00596   coord tree.  First, decrement the size counts from the current position back
00597   to the corr root:*/
00598   for(parent = node->corr_parent; parent != (DNODE *)0; parent = parent->corr_parent)
00599     parent->size--;
00600 /*Make p point to the link that leads into the subtree rooted at node.*/
00601   p = ((node->corr_parent == (DNODE *)0)? &(t->corr_root):
00602         (node == node->corr_parent->corr_lchild)?
00603                 &(node->corr_parent->corr_lchild):
00604                 &(node->corr_parent->corr_rchild));
00605 /*The cases of fewer than two children are easy to handle:*/
00606   if(node->corr_lchild == (DNODE *)0)
00607     {
00608     if(node->corr_rchild)
00609       node->corr_rchild->corr_parent = node->corr_parent;
00610     *p = node->corr_rchild;
00611     }
00612   else if(node->corr_rchild == (DNODE *)0)
00613     {
00614     if(node->corr_lchild)
00615       node->corr_lchild->corr_parent = node->corr_parent;
00616     *p = node->corr_lchild;
00617     }
00618 /*In the two-child case, we find the rightmost node of the left subtree and
00619   promote it.*/
00620   else
00621     {
00622   /*inv: *q is the rightmost node at the current level in the left subtree of
00623     (*p).*/
00624     for(q = &(node->corr_lchild); (*q)->corr_rchild != (DNODE *)0; q = &((*q)->corr_rchild))
00625       (*q)->size--;
00626     (*q)->corr_rchild = node->corr_rchild;
00627     if(node->corr_rchild)
00628       node->corr_rchild->corr_parent = *q;
00629     *p = *q;
00630     parent = (*q)->corr_parent;
00631     (*q)->corr_parent = node->corr_parent;
00632     if(q != &(node->corr_lchild))
00633       {
00634   /*If the left subtree of the rightmost node is non-null, promote it to fill
00635     the place that has just been vacated by the promoted rightmost node.*/
00636       temp = (*q)->corr_lchild;
00637       (*q)->corr_lchild = node->corr_lchild;
00638       node->corr_lchild->corr_parent = *q;
00639       *q = temp;
00640       if(temp)
00641         temp->corr_parent = parent;
00642       }
00643     (*p)->size = 1
00644         + ((*p)->corr_lchild? (*p)->corr_lchild->size: 0)
00645         + ((*p)->corr_rchild? (*p)->corr_rchild->size: 0);
00646     }
00647   }
00648 
00649 /*Delete from the dtree all nodes that are labelled with the given coordinates.
00650   For each of the deleted nodes, if any substitute coordinates remain in the
00651   coordinate list, use those substitute coordinates to re-insert the node into
00652   the dtree.  Otherwise increment num_coords_exhausted.*/
00653 static void dtree_delete_coords(DTREE *t, int x, int y, int z, short *deleted)
00654 /*deleted[x+xdim*(y+ydim*z)] is nonzero iff (x,y,z) has not yet been deleted*/
00655   {
00656   register DNODE **p,                   /*pointer into coord tree*/
00657                  *del;                  /*pointer to deleted node*/
00658   int target_coords_not_found,          /*flag indicates coords not yet found*/
00659       coords;                           /*encoded form of target coordinates*/
00660   coords = coord_encode(x, y, z, t->xdim, t->ydim);
00661   p = &(t->coord_root);
00662   target_coords_not_found = 1;
00663 /*inv: All remaining nodes whose coordinates match the target coords are in the
00664   subtree rooted at *p.  All matching nodes outside this subtree have been
00665   replaced by nodes with substitute coordinates and correlations.*/
00666   while(target_coords_not_found && (*p != (DNODE *)0))
00667     {
00668     if(coords > (*p)->data->coords)
00669       p = &((*p)->coord_rchild);
00670     else if(coords < (*p)->data->coords)
00671       p = &((*p)->coord_lchild);
00672     else
00673       {
00674     /*found a matching node*/
00675       target_coords_not_found = 0;
00676       do
00677         {
00678         del = *p;
00679         dtree_unlink_node(t, p); /*note: *p != del after this call!*/
00680       /*re-insert the node with substitute data, if any substitutes remain*/
00681         do {del->dlen--; del->data++;
00682       /*check that substitutes haven't been exhausted (dlen != 0) and the
00683         coordinates associated with the substitute at the head of the list are
00684         still not valid (deleted[del->data->coords] != 0)*/
00685            } while((del->dlen != 0) && deleted[del->data->coords]);
00686       /*If valid substitute data are available, re-insert this node using those
00687         substitutes.  Otherwise, increment the count of occurrences of failed
00688         substitutions.*/
00689         if(del->dlen != 0)
00690           dtree_insert_at_node(t, del);
00691         else
00692           num_coords_exhausted++;
00693         } while(*p && (coords == (*p)->data->coords));
00694 /*By clearing the flag target_coords_not_found on entry to the inner loop, we
00695   ensure that the outer loop will terminate once the inner loop's guard has
00696   become false.  The correctness of this strategy depends on the structure of
00697   the auxiliary routine dtree_unlink_node() and on the following line in
00698   dtree_insert_at_node() above:
00699     p = ((node->data->coords > (*p)->data->coords)? &((*p)->coord_rchild):
00700                                                     &((*p)->coord_lchild));
00701   This line makes sets of equal coordinates into *left*-recursive rather than
00702   right-recursive trees, and accordingly, dtree_unlink_node() looks in the
00703   *left* subtree for a node to promote into the place of the deleted node as
00704   (*p).*/
00705       }
00706     }
00707   }
00708 
00709 /*Return a fraction in the interval (0,1) that describes the ordinal position
00710   that the given correlation value would occupy within the ordered contents of
00711   the dtree, were it to be added to the tree.*/
00712 static double dtree_position(DTREE *t, double r)
00713   {
00714   register DNODE *p;
00715   register int rank;
00716   p = t->corr_root;
00717   if((p == (DNODE *)0) || (p->size == 1))
00718     return(0.5);
00719   rank = 0;
00720 /*inv: rank/2 is the number of correlation values in the tree that are less
00721   than r and are not within the subtree rooted at p, plus half of the number
00722   of correlation values in the tree that equal r and are not within the subtree
00723   rooted at p.*/
00724   while(p != (DNODE *)0)
00725     if(r > p->data->corr)
00726       {
00727       rank += 2*(1+(p->corr_lchild? p->corr_lchild->size: 0));
00728       p = p->corr_rchild;
00729       }
00730     else
00731       {
00732       if(r == p->data->corr)
00733         rank++; /*values equal to the target count only half*/
00734       p = p->corr_lchild;
00735       }
00736   if(rank == 2*t->corr_root->size)
00737     rank--;
00738   else if(rank == 0)
00739     rank++;
00740   return(rank / (2.0*t->corr_root->size));
00741   }
00742 
00743 /*Randomly permute the integer array a[0..n-1].
00744   Richard Durstenfeld, `Algorithm 235: Random Permutation [G6]',
00745   Communications of the Association for Computing Machinery 7:7:420 (1964).*/
00746 static void shuffle(int *a, int n)
00747   {
00748   register int i, j, b;
00749   for(i = n-1; i != 0; i--)
00750     {
00751     j = (int)(random()%(i+1));
00752     b = a[i]; a[i] = a[j]; a[j] = b;
00753     }
00754   }
00755 
00756 /*Prepare static sums for correlate().*/
00757 static void correlate_prep(float *x, float *y, int nxy, double *sy, double *syy)
00758   {
00759   register int n;       /*length of time series, excluding ignored points*/
00760   *sy = *syy = 0.0;
00761   n = 0;
00762   while(nxy--)
00763     {
00764     if(*x < BLAST)
00765       {
00766       *syy += *y**y;
00767       *sy += *y;
00768       n++;
00769       }
00770     x++; y++;
00771     }
00772   *syy -= *sy**sy/n;
00773   }
00774 
00775 /*Correlate the ideal time series x with the actual time series y, where both
00776   series are of length nxy.  If the pointers r, m, and b are non-null, they
00777   receive, respectively, the correlation coefficient and the slope and
00778   intercept of the regression line.  A call to correlate() must be preceded by
00779   calls to correlate_prep() of the following forms:
00780   correlate_prep(x, y, nxy, &sx, &sxx);
00781   correlate_prep(x, y, nxy, &sy, &syy);*/
00782 static void correlate(
00783         float *x,       /*ideal time series (points >= BLAST are ignored)*/
00784         float *y,       /*actual time series*/
00785         int nxy,        /*length of time series (including any ignored points)*/
00786         double sx,      /*sum of x[i] | 0<=i<nxy and x[i]<BLAST*/
00787         double sxx,     /*sum of x[i]^2 | 0<=i<nxy and x[i]<BLAST*/
00788         double sy,      /*sum of y[i] | 0<=i<nxy and x[i]<BLAST*/
00789         double syy,     /*sum of y[i]^2 | 0<=i<nxy and x[i]<BLAST*/
00790         double *r,      /*correlation coefficient*/
00791         double *m,      /*slope*/
00792         double *b)      /*intercept*/
00793   {
00794   register double sxy;
00795   register int n;       /*length of time series, excluding ignored points*/
00796   double slope;
00797   if(sxx == 0.0)
00798     {
00799     if(r) *r = 0.0;
00800     if(m) *m = 0.0;
00801     if(b)
00802       {
00803       n = 0;
00804       while(nxy--)
00805         if(x[nxy] < BLAST)
00806           n++;
00807       *b = sy/n;
00808       }
00809     }
00810   else
00811     {
00812     sxy = 0.0;
00813     n = 0;
00814     while(nxy--)
00815       {
00816       if(*x < BLAST)
00817         {
00818         sxy += *x**y;
00819         n++;
00820         }
00821       x++; y++;
00822       }
00823     sxy -= sx*sy/n;
00824     if(r)
00825       *r = ((syy == 0.0)? 0.0: (sxy/sqrt(sxx*syy)));
00826     slope = sxy/sxx;
00827     if(m) *m = slope;
00828     if(b) *b = (sy - slope*sx)/n;
00829     }
00830   }
00831 
00832 /*Orthogonalise y to x, where m and b have been set by a call to correlate().
00833   Store the result in new_y (which may point to the same location as y, in the
00834   case in which the original value of y need not be saved).  For linear
00835   detrending, input x[i] = i, 0<=i<n.*/
00836 static void orthogonalise(float *x, float *y, float *new_y, int n, double m, double b)
00837   {
00838   while(n--)
00839     if(x[n] < BLAST)
00840       new_y[n] = y[n]-(m*x[n]+b);
00841   }
00842 
00843 /*Compute the permutation test.  Return 0 if successful, 1 if error.*/
00844 static int PERMTEST_compute(
00845   PLUGIN_interface *plint,      /*AFNI plugin interface*/
00846   THD_3dim_dataset *dset,       /*AFNI 3D+time dataset*/
00847   MRI_IMAGE *ref_ts,            /*reference time series*/
00848   MRI_IMAGE *ort_ts,            /*time series against which to orthogonalise*/
00849   double pcrit,                 /*critical probability*/
00850   int one_tailed,       /*-1=one-tailed (-), 0=two-tailed, 1=one-tailed (+)*/
00851         /*returned parameters:*/
00852   short **intensities,          /*address of pointer to intensities*/
00853   short **zvals,                /*address of pointer to z-scores*/
00854   double *fim_scale,            /*scaling factor for intensities*/
00855   THD_3dim_dataset *mask,       /*mask dataset, or NULL for no mask*/
00856   float masklo,                 /*lower and upper bounds on masked range*/
00857   float maskhi,
00858   int verbose)          /*1 for verbose info on coordinates, 0 otherwise*/
00859   {
00860   register int t, iter, xyz;    /*indices and counters*/
00861   int x, y, z,
00862       xdim, ydim, zdim,         /*spatial dimensions*/
00863       tdim,                     /*temporal dimension of dataset*/
00864       tsize,                    /*# of included (< BLAST) points, tsize<=tdim*/
00865       ignore,                   /*# of ignored (>=BLAST) pts at head of series*/
00866       t2,
00867       *tindex,                  /*temporal sequence of all included images*/
00868       *sequence;                /*permuted temporal sequence*/
00869   float *vox_xyz,               /*time series for one voxel*/
00870         *vox,                   /*3D+time data indexed by z,y,x,t*/
00871         *ts,                  /*storage for linear series (ts[t]=t, 0<=t<tsize),
00872                                  and later for randomised time series*/
00873         *fit_coeff,             /*fit coefficients*/
00874         mask_val;               /*value of mask at current voxel*/
00875   STREE *actual_corr;
00876   DTREE *randomised_corr;
00877   CLIST *coord_mem,                     /*data area for coord-corr pair lists*/
00878         *coord_next;                    /*next free location in coord_mem*/
00879   double sx_trend, sx_ref, sx_ort, *sy, /*correlation parameters*/
00880          sxx_trend, sxx_ref, sxx_ort, *syy,
00881          corr, slope, intercept,        /*regression coefficients*/
00882          max_corr;                      /*correlation w/ max absolute value*/
00883   int percent_done, prev_percent_done,  /*statistics for progress meter*/
00884       not_done;
00885 /*set dimensions*/
00886   xdim = dset->daxes->nxx;
00887   ydim = dset->daxes->nyy;
00888   zdim = dset->daxes->nzz;
00889   tdim = DSET_NUM_TIMES(dset);
00890   switch(one_tailed)
00891     {
00892     case -1: magnitude = neg; break;
00893     case 0: magnitude = fabs; break;
00894     case 1: magnitude = id;
00895     }
00896   num_coords = NUM_COORDS;
00897   num_coords_exhausted = 0;
00898   if(ort_ts != (MRI_IMAGE *)0) /*need to save extra points if orthogonalising*/
00899     num_coords *= 2;
00900 /*allocate storage for functional intensities*/
00901   *intensities = (short *)calloc(xdim*ydim*zdim, sizeof(**intensities));
00902   if(*intensities == (short *)0)
00903     return 1;
00904 /*allocate storage for z-scores*/
00905   *zvals = (short *)calloc(xdim*ydim*zdim, sizeof(**zvals));
00906   if(*zvals == (short *)0)
00907     {
00908     free(*intensities);
00909     return 1;
00910     }
00911   /*no need to initialise (*zvals)[] to zeroes, since calloc() does it for us*/
00912   /*Initialise temporal sequences:  `tindex' is a mapping from the time index
00913     used in statistical tests to the time index used for storage in the AFNI
00914     dataset; it excludes ignored time points.  `sequence' is initially the same
00915     as tindex but will be permuted during each iteration of the resampling
00916     procedure.*/
00917   tindex = calloc(tdim, sizeof(*tindex));
00918   if(tindex == (int *)0)
00919     {
00920     free(*zvals);
00921     free(*intensities);
00922     return 1;
00923     }
00924   sequence = calloc(tdim, sizeof(*sequence));
00925   if(sequence == (int *)0)
00926     {
00927     free(tindex);
00928     free(*zvals);
00929     free(*intensities);
00930     return 1;
00931     }
00932   for((t = 0), (tsize = 0); t != tdim; t++)
00933     if(MRI_FLOAT_PTR(ref_ts)[t] < BLAST)
00934       {
00935       tindex[tsize] = t;
00936       sequence[tsize++] = t;
00937       }
00938   if(tsize < 4)         /*make sure we have enough points for regression*/
00939     {
00940     free(sequence);
00941     free(tindex);
00942     free(*zvals);
00943     free(*intensities);
00944     return 1;
00945     }
00946 /*initialise data*/
00947   vox = calloc(xdim*ydim*zdim*tdim, sizeof(*vox));
00948   if(vox == (float *)0)
00949     {
00950     free(sequence);
00951     free(tindex);
00952     free(*zvals);
00953     free(*intensities);
00954     return 1;
00955     }
00956   sy = calloc(xdim*ydim*zdim, sizeof(*sy));
00957   if(sy == (double *)0)
00958     {
00959     free(vox);
00960     free(sequence);
00961     free(tindex);
00962     free(*zvals);
00963     free(*intensities);
00964     return 1;
00965     }
00966   syy = calloc(xdim*ydim*zdim, sizeof(*syy));
00967   if(syy == (double *)0)
00968     {
00969     free(sy);
00970     free(vox);
00971     free(sequence);
00972     free(tindex);
00973     free(*zvals);
00974     free(*intensities);
00975     return 1;
00976     }
00977   ts = calloc(tdim, sizeof(*ts));
00978   if(ts == (float *)0)
00979     {
00980     free(syy);
00981     free(sy);
00982     free(vox);
00983     free(sequence);
00984     free(tindex);
00985     free(*zvals);
00986     free(*intensities);
00987     return 1;
00988     }
00989   actual_corr = stree_create(xdim, ydim, xdim*ydim*zdim);
00990   if(actual_corr == (STREE *)0)
00991     {
00992     free(syy);
00993     free(sy);
00994     free(ts);
00995     free(vox);
00996     free(sequence);
00997     free(tindex);
00998     free(*zvals);
00999     free(*intensities);
01000     return 1;
01001     }
01002   fit_coeff = (float *)calloc(xdim*ydim*zdim, sizeof(*fit_coeff));
01003   if(fit_coeff == (float *)0)
01004     {
01005     stree_destroy(actual_corr);
01006     free(syy);
01007     free(sy);
01008     free(ts);
01009     free(vox);
01010     free(sequence);
01011     free(tindex);
01012     free(*zvals);
01013     free(*intensities);
01014     return 1;
01015     }
01016 
01017 /******************************************************************************
01018  * PHASE 1: Copy data, remove linear trend, optionally orthogonalise against  *
01019  *          a user-supplied series, and correlate with the ideal series.      *
01020  ******************************************************************************/
01021   for(ignore = 0; MRI_FLOAT_PTR(ref_ts)[ignore] >= BLAST; ignore++)
01022     ;
01023 /*linear series to use in detrending*/
01024   for(t = 0; t != tdim; t++)
01025     ts[t] = (float)t;
01026 /*sum of the linear series, and sum of its squares*/
01027   sx_trend = (double)(((tdim-1)*tdim - (ignore-1)*ignore)/2);
01028   sxx_trend = ((double)((tdim*(tdim*(2*tdim-3)+1) - ignore*(ignore*(2*ignore-3)+1))/6)) - sx_trend*sx_trend/(tdim-ignore);
01029 /*compute the sum and sum of squares of the ideal time series*/
01030   correlate_prep(MRI_FLOAT_PTR(ref_ts)+ignore, MRI_FLOAT_PTR(ref_ts)+ignore, tdim-ignore, &sx_ref, &sxx_ref);
01031 /*compute the sum and sum of squares of the orthogonal series, if specified*/
01032   if(ort_ts != (MRI_IMAGE *)0)
01033     correlate_prep(MRI_FLOAT_PTR(ort_ts), MRI_FLOAT_PTR(ort_ts), tdim, &sx_ort, &sxx_ort);
01034   *fim_scale = 0.0;
01035   mask_val = masklo;
01036   xyz = 0;
01037   for(z = 0; z !=  zdim; z++)
01038     for(y = 0; y != ydim; y++)
01039       for(x = 0; x != xdim; x++)
01040         {
01041         if(mask != NULL)
01042           {
01043           mask_val = DSET_BRICK_FACTOR(mask, 0);
01044           if(mask_val == 0.0)
01045             mask_val = 1.0;
01046           if(DSET_BRICK_TYPE(mask, 0) == MRI_short)
01047             mask_val *= ((short *)DSET_ARRAY(mask, 0))[xyz];
01048           else
01049             mask_val *= ((char *)DSET_ARRAY(mask, 0))[xyz];
01050           }
01051         if((mask_val >= masklo) && (mask_val <= maskhi))
01052           {
01053           vox_xyz = vox + tdim*xyz;
01054           for(t = 0; t != tdim; t++)
01055             {
01056             vox_xyz[t] = DSET_BRICK_FACTOR(dset, t);
01057             if(vox_xyz[t] == 0.0)
01058               vox_xyz[t] = 1.0;
01059             vox_xyz[t] *= ((short *)DSET_ARRAY(dset, t))[xyz];
01060             }
01061         /*detrend*/
01062           correlate_prep(ts+ignore, vox_xyz+ignore, tdim-ignore, sy+xyz, syy+xyz);
01063           correlate(ts+ignore, vox_xyz+ignore, tdim-ignore, sx_trend, sxx_trend, sy[xyz], syy[xyz], (double *)0, &slope, &intercept);
01064           orthogonalise(ts, vox_xyz, vox_xyz, tdim, slope, intercept+slope*ignore);
01065         /*remove the orthogonalisation time series, if one has been supplied*/
01066           if(ort_ts != (MRI_IMAGE *)0)
01067             {
01068             correlate_prep(MRI_FLOAT_PTR(ort_ts), vox_xyz, tdim, sy+xyz, syy+xyz);
01069             correlate(MRI_FLOAT_PTR(ort_ts), vox_xyz, tdim, sx_ort, sxx_ort, sy[xyz], syy[xyz], (double *)0, &slope, &intercept);
01070             orthogonalise(MRI_FLOAT_PTR(ort_ts), vox_xyz, vox_xyz, tdim, slope, intercept);
01071             }
01072         /*correlate*/
01073           correlate_prep(MRI_FLOAT_PTR(ref_ts)+ignore, vox_xyz+ignore, tdim-ignore, sy+xyz, syy+xyz);
01074           correlate(MRI_FLOAT_PTR(ref_ts)+ignore, vox_xyz+ignore, tdim-ignore, sx_ref, sxx_ref, sy[xyz], syy[xyz], &corr, &slope, (double *)0);
01075         /*save regression coefficient*/
01076           fit_coeff[xyz] = slope;
01077         /*retain greatest slope for use in computing scaling factor*/
01078           slope = fabs(slope);
01079           if(slope > *fim_scale)
01080             *fim_scale = slope;
01081         /*save correlation coefficient*/
01082           stree_insert(actual_corr, corr, x, y, z);
01083           }
01084         xyz++;
01085         }
01086   /*make sure that the mask left us *some* voxels*/
01087   if(stree_null(actual_corr))
01088     {
01089     fprintf(stderr, "PERMTEST: All voxels are masked!\n");
01090     return 1;
01091     }
01092   *fim_scale /= MRI_TYPE_maxval[MRI_short];
01093 /*copy fit coefficients for all voxels, scaled for storage in shorts*/
01094   for(t = 0; t != xdim*ydim*zdim; t++)
01095     (*intensities)[t] = (short)(fit_coeff[t] / *fim_scale);
01096   free(fit_coeff);
01097   randomised_corr = dtree_create(xdim, ydim, NUM_ITERS);
01098   if(randomised_corr == (DTREE *)0)
01099     {
01100     stree_destroy(actual_corr);
01101     free(syy);
01102     free(sy);
01103     free(ts);
01104     free(vox);
01105     free(sequence);
01106     free(tindex);
01107     free(*zvals);
01108     free(*intensities);
01109     return 1;
01110     }
01111   coord_mem = (CLIST *)calloc(NUM_ITERS*num_coords, sizeof(*coord_mem));
01112   if(coord_mem == (CLIST *)0)
01113     {
01114     dtree_destroy(randomised_corr);
01115     stree_destroy(actual_corr);
01116     free(syy);
01117     free(sy);
01118     free(ts);
01119     free(vox);
01120     free(sequence);
01121     free(tindex);
01122     free(*zvals);
01123     free(*intensities);
01124     return 1;
01125     }
01126   coord_next = coord_mem;
01127 /*initialise progress meter*/
01128   prev_percent_done = -1;
01129   PLUTO_popup_meter(plint);
01130 
01131 /******************************************************************************
01132  * PHASE 2: Compute the empirical distribution.                               *
01133  ******************************************************************************/
01134   for(iter = 0; iter != NUM_ITERS; iter++)
01135     {
01136     percent_done = (100*iter)/NUM_ITERS;
01137     if(percent_done > prev_percent_done)
01138       {
01139       prev_percent_done = percent_done;
01140       PLUTO_set_meter(plint, percent_done);
01141       }
01142     shuffle(sequence, tsize);
01143     for(t = 0; t != num_coords; t++)
01144       coord_next[t].corr = 0.0;
01145     xyz = 0;
01146     for(z = 0; z !=  zdim; z++)
01147       for(y = 0; y != ydim; y++)
01148         for(x = 0; x != xdim; x++)
01149           {
01150           if(mask != NULL)
01151             {
01152             mask_val = DSET_BRICK_FACTOR(mask, 0);
01153             if(mask_val == 0.0)
01154               mask_val = 1.0;
01155             if(DSET_BRICK_TYPE(mask, 0) == MRI_short)
01156               mask_val *= ((short *)DSET_ARRAY(mask, 0))[xyz];
01157             else
01158               mask_val *= ((char *)DSET_ARRAY(mask, 0))[xyz];
01159             }
01160           if((mask_val >= masklo) && (mask_val <= maskhi))
01161             {
01162             vox_xyz = vox + tdim*xyz;
01163             for((t = ignore), (t2 = 0); t != tdim; t++)
01164               if(MRI_FLOAT_PTR(ref_ts)[t] < BLAST)
01165                 ts[t] = vox_xyz[sequence[t2++]];
01166             correlate(MRI_FLOAT_PTR(ref_ts)+ignore, ts+ignore, tdim-ignore, sx_ref, sxx_ref, sy[xyz], syy[xyz], &corr, (double *)0, (double *)0);
01167           /*Since num_coords is tiny, any sorting algorithm faster
01168             than linear insertion wouldn't be worth the trouble.*/
01169             for(t = 0; t < num_coords; t++)
01170               if((*magnitude)(corr) > (*magnitude)(coord_next[t].corr))
01171                 {
01172                 t2 = t;
01173                 while(++t != num_coords)
01174                   {
01175                   coord_next[t].corr = coord_next[t-1].corr;
01176                   coord_next[t].coords = coord_next[t-1].coords;
01177                   }
01178                 coord_next[t2].corr = corr;
01179                 coord_next[t2].coords = coord_encode(x, y, z, xdim, ydim);
01180                 }
01181             }
01182           xyz++;
01183           }
01184     dtree_insert(randomised_corr, coord_next);
01185     if(verbose)
01186       printf("iter=%d max_corr=%f\n", iter, coord_next->corr);
01187     coord_next += num_coords;
01188     }
01189 #ifdef PERMTEST_DEBUG
01190   dtree_traverse(randomised_corr, 0);   /*dump distribution to stdout*/
01191 #endif
01192   PLUTO_set_meter(plint, 100);
01193 
01194 /******************************************************************************
01195  * PHASE 3: Rank correlations from the actual data set within the empirical   *
01196  *          distribution.                                                     *
01197  ******************************************************************************/
01198   do
01199     {
01200     not_done = 0;
01201   /*extract the most significant *correlation* (in [-1,1])*/
01202     max_corr = stree_extract_max(actual_corr, &x, &y, &z);
01203   /*map it to an adjusted *probability* (in [0,1])*/
01204     corr = dtree_position(randomised_corr, max_corr);
01205     if(verbose)
01206       printf("(%2d,%2d,%2d) raw r=%f adjusted p=%f", x, y, z, max_corr, corr);
01207     if(((one_tailed == 0) && (fabs(corr-0.5) >= 0.5-pcrit/2.0))
01208      ||((one_tailed == -1) && (corr <= pcrit))
01209      ||((one_tailed == 1) && (corr >= 1.0-pcrit)))
01210         {
01211         slope = critz(corr);
01212         if(verbose)
01213           printf(" zcrit=%f\n", slope);
01214         (*zvals)[x + xdim*(y + ydim*z)] = (short)(slope*FUNC_ZT_SCALE_SHORT);
01215         not_done = 1;
01216 /******************************************************************************
01217  *percent completion of this phase of the algorithm can be expressed as follows:
01218  *      100.0 * ((one_tailed == 1)? 1.0-corr:                                 *
01219  *                ((one_tailed == 0) && (corr > 0.5))?                        *
01220  *                  2.0*(1.0-corr):                                           *
01221  *                    corr) / pcrit;                                          *
01222  ******************************************************************************/
01223       /*since this point is now defined as activated, exclude its data from the
01224         empirical distribution*/
01225         dtree_delete_coords(randomised_corr, x, y, z, *zvals);
01226         }
01227     } while(not_done && (dtree_size(randomised_corr) > NUM_ITERS/2));
01228                         /*^^^ stop if we've deleted too many coordinates ^^^*/
01229   if(verbose)
01230     putchar('\n');
01231   free(coord_mem);
01232   dtree_destroy(randomised_corr);
01233   stree_destroy(actual_corr);
01234   free(syy);
01235   free(sy);
01236   free(ts);
01237   free(vox);
01238   free(sequence);
01239   free(tindex);
01240   return 0;
01241   }
01242 
01243 /*In the event of a conflict with some future revision of the AFNI widgets
01244   code, PERMTEST_set_logo() and PERMTEST_reset_logo() can be redefined as null
01245   functions and everything will still work right (though it would be a shame to
01246   lose the free advertising!).*/
01247 
01248 static Pixmap sipb_pixmap;              /*pixmap of SIPB logo*/
01249 
01250 static void PERMTEST_set_logo(PLUGIN_interface *plint)
01251   {
01252   Pixel bg_pix, fg_pix;                 /*colours from control window*/
01253 #define sipb_width 90
01254 #define sipb_height 90
01255   static char sipb_bits[] = {           /*bitmap of SIPB logo*/
01256   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01257   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01258   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01259   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01260   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01261   0x1f, 0x00, 0x00, 0x30, 0x00, 0x00, 0x00, 0x18, 0x00, 0x00, 0xff, 0x03, 
01262   0x1f, 0x00, 0x00, 0x30, 0x00, 0x00, 0x00, 0x18, 0x00, 0x00, 0xf8, 0x03, 
01263   0x1f, 0x00, 0x00, 0x30, 0x00, 0x00, 0x00, 0x18, 0x00, 0x00, 0xf0, 0x03, 
01264   0x1f, 0x00, 0x00, 0x30, 0x00, 0x00, 0x00, 0x18, 0x00, 0x00, 0xf0, 0x03, 
01265   0x3f, 0x00, 0x0c, 0xf0, 0x00, 0x7c, 0x00, 0x7e, 0x00, 0x03, 0xe0, 0x03, 
01266   0x7f, 0x00, 0x18, 0xf0, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x0f, 0xe0, 0x03, 
01267   0x7f, 0x00, 0x30, 0xf0, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x1f, 0xe0, 0x03, 
01268   0xff, 0x00, 0x60, 0xf0, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x1f, 0xe0, 0x03, 
01269   0xff, 0x01, 0x60, 0xf0, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x0f, 0xe0, 0x03, 
01270   0xff, 0x03, 0xc0, 0xff, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x0f, 0xe0, 0x03, 
01271   0xff, 0x07, 0x80, 0xff, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x03, 0xe0, 0x03, 
01272   0xff, 0x0f, 0x00, 0xff, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x00, 0xf0, 0x03, 
01273   0xff, 0x1f, 0x00, 0xff, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x00, 0xf8, 0x03, 
01274   0xff, 0x3f, 0x00, 0xff, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x00, 0xe0, 0x03, 
01275   0xff, 0x3f, 0x00, 0xff, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x00, 0xe0, 0x03, 
01276   0xff, 0x1f, 0x80, 0xff, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x03, 0xc0, 0x03, 
01277   0xff, 0x0f, 0xc0, 0xff, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x0f, 0xc0, 0x03, 
01278   0xff, 0x07, 0xe0, 0xff, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x1f, 0xc0, 0x03, 
01279   0xff, 0x03, 0x70, 0xf0, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x3f, 0xc0, 0x03, 
01280   0xff, 0x01, 0x78, 0xf0, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x3f, 0xc0, 0x03, 
01281   0xff, 0x00, 0x3c, 0xf0, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x1f, 0xc0, 0x03, 
01282   0x7f, 0x00, 0x00, 0xf0, 0x01, 0xfe, 0x00, 0xff, 0x00, 0x0f, 0xc0, 0x03, 
01283   0x3f, 0x00, 0x00, 0xf0, 0x00, 0x7c, 0x00, 0x7e, 0x00, 0x07, 0xc0, 0x03, 
01284   0x1f, 0x00, 0x00, 0x30, 0x00, 0x10, 0x00, 0x18, 0x00, 0x00, 0xe0, 0x03, 
01285   0x1f, 0x00, 0x00, 0x30, 0x00, 0x10, 0x00, 0x18, 0x00, 0x00, 0xf0, 0x03, 
01286   0x1f, 0x00, 0x00, 0x30, 0x00, 0x10, 0x00, 0x18, 0x00, 0x00, 0xf8, 0x03, 
01287   0x1f, 0x00, 0x00, 0x30, 0x00, 0x10, 0x00, 0x18, 0x00, 0x00, 0xff, 0x03, 
01288   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01289   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01290   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01291   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01292   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01293   0xff, 0xff, 0x40, 0xfe, 0x3f, 0x9f, 0x04, 0xfe, 0xff, 0xff, 0xff, 0x03, 
01294   0xff, 0xff, 0x73, 0xfe, 0x3f, 0x8e, 0x9c, 0xff, 0xff, 0xff, 0xff, 0x03, 
01295   0xff, 0xff, 0x73, 0x38, 0x3e, 0x8e, 0x9c, 0xff, 0xff, 0xff, 0xff, 0x03, 
01296   0xff, 0xff, 0x73, 0x92, 0x3c, 0x84, 0x9c, 0xff, 0xff, 0xff, 0xff, 0x03, 
01297   0xff, 0xff, 0x73, 0x12, 0x3c, 0x95, 0x9c, 0xff, 0xff, 0xff, 0xff, 0x03, 
01298   0xff, 0xff, 0x73, 0x92, 0x3f, 0x91, 0x9c, 0xff, 0xff, 0xff, 0xff, 0x03, 
01299   0xff, 0xff, 0x73, 0x92, 0x3c, 0x9b, 0x9c, 0xff, 0xff, 0xff, 0xff, 0x03, 
01300   0xff, 0xff, 0x73, 0x32, 0x3e, 0x9b, 0x9c, 0xff, 0xff, 0xff, 0xff, 0x03, 
01301   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01302   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01303   0xff, 0xff, 0xc3, 0xfc, 0x9f, 0xff, 0xcf, 0xff, 0xff, 0xff, 0xff, 0x03, 
01304   0xff, 0xff, 0x99, 0xfc, 0x9f, 0xff, 0xcf, 0xff, 0xff, 0xff, 0xff, 0x03, 
01305   0xff, 0xff, 0xf1, 0x48, 0x86, 0xb1, 0x8c, 0xff, 0xff, 0xff, 0xff, 0x03, 
01306   0xff, 0xff, 0xc3, 0x4c, 0x92, 0x24, 0xc9, 0xff, 0xff, 0xff, 0xff, 0x03, 
01307   0xff, 0xff, 0x8f, 0x4c, 0x9a, 0x20, 0xc9, 0xff, 0xff, 0xff, 0xff, 0x03, 
01308   0xff, 0xff, 0x9d, 0x4c, 0x9a, 0x3c, 0xc9, 0xff, 0xff, 0xff, 0xff, 0x03, 
01309   0xff, 0xff, 0x99, 0x4c, 0x92, 0x24, 0xc9, 0xff, 0xff, 0xff, 0xff, 0x03, 
01310   0xff, 0xff, 0xc3, 0x99, 0x86, 0x31, 0x99, 0xff, 0xff, 0xff, 0xff, 0x03, 
01311   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01312   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01313   0xff, 0xff, 0xe1, 0x1f, 0xff, 0xff, 0xff, 0x99, 0xff, 0xff, 0xff, 0x03, 
01314   0xff, 0xff, 0xf3, 0xcf, 0xff, 0xff, 0xff, 0xf9, 0xff, 0xff, 0xff, 0x03, 
01315   0xff, 0xff, 0xb3, 0x8c, 0xb1, 0x48, 0x8e, 0x91, 0xb1, 0xfc, 0xff, 0x03, 
01316   0xff, 0xff, 0x33, 0xc9, 0x24, 0x92, 0x34, 0x99, 0x24, 0xf9, 0xff, 0x03, 
01317   0xff, 0xff, 0x33, 0xc9, 0x24, 0x93, 0x0c, 0x99, 0x24, 0xf9, 0xff, 0x03, 
01318   0xff, 0xff, 0x33, 0xc9, 0x24, 0x93, 0x24, 0x99, 0x24, 0xf9, 0xff, 0x03, 
01319   0xff, 0xff, 0x33, 0xc9, 0x24, 0x93, 0x24, 0x99, 0x24, 0xf9, 0xff, 0x03, 
01320   0xff, 0xff, 0x21, 0xc9, 0x31, 0x93, 0x4c, 0x92, 0x31, 0xf9, 0xff, 0x03, 
01321   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01322   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01323   0xff, 0xff, 0xc1, 0xff, 0xff, 0xff, 0xff, 0xf9, 0xff, 0xff, 0xff, 0x03, 
01324   0xff, 0xff, 0x99, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01325   0xff, 0xff, 0x99, 0x12, 0xc7, 0x38, 0x8e, 0x29, 0xa7, 0xff, 0xff, 0x03, 
01326   0xff, 0xff, 0x99, 0x48, 0x52, 0x92, 0x24, 0x49, 0x92, 0xff, 0xff, 0x03, 
01327   0xff, 0xff, 0xc1, 0x4c, 0x72, 0x30, 0x8e, 0x49, 0x9a, 0xff, 0xff, 0x03, 
01328   0xff, 0xff, 0xf9, 0x4c, 0x72, 0xfe, 0x3c, 0x49, 0x9a, 0xff, 0xff, 0x03, 
01329   0xff, 0xff, 0xf9, 0x4c, 0x52, 0x92, 0x24, 0x49, 0x92, 0xff, 0xff, 0x03, 
01330   0xff, 0xff, 0xf9, 0x1c, 0xc7, 0x38, 0x8e, 0x49, 0x86, 0xff, 0xff, 0x03, 
01331   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x9f, 0xff, 0xff, 0x03, 
01332   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xc7, 0xff, 0xff, 0x03, 
01333   0xff, 0xff, 0xc1, 0xff, 0xff, 0xe7, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01334   0xff, 0xff, 0x99, 0xff, 0xff, 0xe7, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01335   0xff, 0xff, 0x99, 0x71, 0x2c, 0xe1, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01336   0xff, 0xff, 0xc1, 0xa4, 0x89, 0xe4, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01337   0xff, 0xff, 0x99, 0x64, 0xc8, 0xe6, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01338   0xff, 0xff, 0x99, 0x24, 0xc9, 0xe6, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01339   0xff, 0xff, 0x99, 0x24, 0xc9, 0xe4, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01340   0xff, 0xff, 0xc1, 0x71, 0xc2, 0xe1, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01341   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01342   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01343   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01344   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01345   0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0x03, 
01346   };
01347   /*Initialise SIPB logo, in the current controller window only.  Leave all
01348     other controller windows with the default MCW logo.  (Adapted from
01349     afni_widg:AFNI_make_widgets).*/
01350   XtVaGetValues(plint->im3d->vwid->top_form,
01351                 XmNforeground, &bg_pix,     /*note reversal of fg & bg colours*/
01352                 XmNbackground, &fg_pix, NULL);
01353   sipb_pixmap = XCreatePixmapFromBitmapData(
01354                 XtDisplay(plint->im3d->vwid->top_shell),
01355                 RootWindowOfScreen(XtScreen(plint->im3d->vwid->top_shell)),
01356                 sipb_bits, sipb_width, sipb_height,
01357                 fg_pix, bg_pix,
01358                 DefaultDepthOfScreen(XtScreen(plint->im3d->vwid->top_shell)));
01359 /*shameless self-aggrandisement (adapted from afni:AFNI_set_cursor)*/
01360   XtVaSetValues(plint->im3d->vwid->picture, XmNlabelPixmap, sipb_pixmap, NULL);
01361   }
01362 
01363 static void PERMTEST_reset_logo(PLUGIN_interface *plint)
01364   {
01365   XtVaSetValues(plint->im3d->vwid->picture, XmNlabelPixmap, XmUNSPECIFIED_PIXMAP, NULL);
01366   XFreePixmap(XtDisplay(plint->im3d->vwid->top_shell), sipb_pixmap);
01367   }
01368 
01369 char *PERMTEST_main(PLUGIN_interface *plint)
01370   {
01371   int t;
01372   char *prefix;
01373   MRI_IMAGE *ref_time_series, *orthogonal_time_series;
01374   THD_3dim_dataset *dset, *mask, *new_dset;
01375   short *intensities, *zvals;
01376   double fim_scale,                     /*scaling factor for short image data*/
01377          pcrit;                         /*alpha level*/
01378   char *optiontag;                      /*one of tails2, tails1pos, tails1neg*/
01379   int one_tailed;                       /*0, 1, or -1*/
01380   float masklo, maskhi;                 /*bounds on mask dataset*/
01381   static char tails_err[] = "exactly one of two-tailed, one-tailed positive, or one-tailed negative must be chosen";
01382 
01383 #ifdef PERMTEST_DEBUG
01384   void (*handler)(int);
01385 #endif
01386 
01387   if(plint == (PLUGIN_interface *)0)
01388     return "PERMTEST_main: null input";
01389 
01390 /*Make sure that the input dataset exists and is stored in a format that we can
01391   understand (MRI_short)*/
01392   PLUTO_next_option(plint);
01393   dset = PLUTO_find_dset(PLUTO_get_idcode(plint));
01394   if(dset == (THD_3dim_dataset *)0)
01395     return "bad dataset";
01396   for(t = 0; t != dset->dblk->nvals; t++)
01397 /*to do: make copies of this routine that work on MRI_byte and on MRI_float*/
01398     if(DSET_BRICK_TYPE(dset, t) != MRI_short)
01399       return("permutation test on non-short values is not implemented");
01400 
01401 /*Make sure that the time series exists and is stored in a format that we can
01402   understand (MRI_float), and contains enough points to cover the time
01403   dimension of the specified dataset*/
01404   PLUTO_next_option(plint);
01405   ref_time_series = PLUTO_get_timeseries(plint);
01406   if((ref_time_series == (MRI_IMAGE *)0)
01407    ||(ref_time_series->kind != MRI_float)
01408    ||(ref_time_series->nx < DSET_NUM_TIMES(dset)))
01409     return("bad time series");
01410 
01411 /*Read the orthogonalisation time series, if it exists*/
01412   optiontag = PLUTO_get_optiontag(plint);
01413   if(strcmp(optiontag, ort_label))
01414     orthogonal_time_series = (MRI_IMAGE *)0;
01415   else
01416     {
01417     orthogonal_time_series = PLUTO_get_timeseries(plint);
01418     if((orthogonal_time_series == (MRI_IMAGE *)0)
01419      ||(orthogonal_time_series->kind != MRI_float)
01420      ||(orthogonal_time_series->nx < DSET_NUM_TIMES(dset)))
01421       return("bad ort");
01422     PLUTO_next_option(plint);
01423     }
01424 
01425 /*Make sure that the prefix specified for the new dataset is a valid prefix*/
01426   prefix = PLUTO_get_string(plint);
01427   if(!PLUTO_prefix_ok(prefix))
01428     return("bad prefix");
01429 
01430 /*Read the alpha level*/
01431   PLUTO_next_option(plint);
01432   pcrit = PLUTO_get_number(plint);
01433 
01434 /*Exactly one of two-tailed, one-tailed positive, one-tailed negative*/
01435   optiontag = PLUTO_get_optiontag(plint);
01436   if(optiontag == (char *)0)
01437     return(tails_err);
01438   if(strcmp(optiontag, tails2))
01439     {
01440     if(strcmp(optiontag, tails1pos))
01441       {
01442       if(strcmp(optiontag, tails1neg))
01443         return(tails_err);
01444       else
01445         one_tailed = -1;
01446       }
01447     else
01448       one_tailed = 1;
01449     }
01450   else
01451     one_tailed = 0;
01452   optiontag = PLUTO_get_optiontag(plint);
01453   if((optiontag != (char *)0) && strcmp(optiontag, mask_label))
01454     return(tails_err);
01455 
01456 /*Optional mask dataset*/
01457   masklo = 1.0;
01458   maskhi = 32767.0;
01459   if(optiontag == (char *)0)
01460     mask = (THD_3dim_dataset *)0;
01461   else
01462     {
01463     mask = PLUTO_find_dset(PLUTO_get_idcode(plint));
01464     if(mask == (THD_3dim_dataset *)0)
01465       return("bad mask");
01466     if((DSET_BRICK_TYPE(mask, 0) != MRI_short)
01467     && (DSET_BRICK_TYPE(mask, 0) != MRI_byte))
01468       return("mask brick type must be byte or short integer");
01469     optiontag = PLUTO_get_optiontag(plint);
01470     if(optiontag != (char *)0)
01471       {
01472       if(strcmp(optiontag, masklo_label))
01473         maskhi = PLUTO_get_number(plint);
01474       else
01475         {
01476         masklo = PLUTO_get_number(plint);
01477         if(PLUTO_get_optiontag(plint) != (char *)0)
01478           maskhi = PLUTO_get_number(plint);
01479         }
01480       }
01481     DSET_load(mask);
01482     }
01483 
01484 #ifdef PERMTEST_DEBUG
01485   handler = signal(SIGUSR1, flush);
01486 #endif
01487 
01488 /*toot our own horn :-)*/
01489   PERMTEST_set_logo(plint);
01490 /*Make sure source dataset is in memory*/
01491   DSET_load(dset);
01492   if(PERMTEST_compute(plint, dset, ref_time_series, orthogonal_time_series, pcrit, one_tailed, &intensities, &zvals, &fim_scale, mask, masklo, maskhi, 0))
01493     {
01494     PERMTEST_reset_logo(plint);
01495     return("out of memory");
01496     }
01497   if(num_coords_exhausted)
01498     printf("%d of %d points (%d%%) were deleted from the distribution due to\nexhausted coordinates.  If this fraction is larger than about 10%%, consider\nrecompiling plug_permtest.so with a NUM_COORDS larger than the current value of\n%d.\n", num_coords_exhausted, NUM_ITERS, (100*num_coords_exhausted+NUM_ITERS/2)/NUM_ITERS, NUM_COORDS);
01499 
01500 /*create the output dataset*/
01501   new_dset = EDIT_empty_copy(dset);
01502   if(EDIT_dset_items(new_dset,
01503         ADN_prefix, prefix,
01504         ADN_malloc_type, DATABLOCK_MEM_MALLOC, /*hold in r/w memory*/
01505         ADN_datum_all, MRI_short,       /*store as (scaled) short ints*/
01506         ADN_nvals, 2,                   /*2 sub-bricks: intensity + z-score*/
01507         ADN_ntt, 0,                     /*no time dimension*/
01508         ADN_type, ISHEAD(dset)?
01509           HEAD_FUNC_TYPE: GEN_FUNC_TYPE, /*functional image*/
01510         ADN_func_type, FUNC_ZT_TYPE,    /*intensity + z-score*/
01511         ADN_none))
01512     {
01513     PERMTEST_reset_logo(plint);
01514     return("EDIT_dset_items error");
01515     }
01516   EDIT_BRICK_LABEL(new_dset, 0, "Fit Coef");
01517   mri_fix_data_pointer(intensities, DSET_BRICK(new_dset, 0));
01518   EDIT_BRICK_FACTOR(new_dset, 0, (float)fim_scale);
01519   EDIT_BRICK_LABEL(new_dset, 1, "z-score");
01520   EDIT_BRICK_TO_FIZT(new_dset, 1);
01521   mri_fix_data_pointer(zvals, DSET_BRICK(new_dset, 1));
01522   EDIT_BRICK_FACTOR(new_dset, 1, (float)(1.0/FUNC_ZT_SCALE_SHORT));
01523   DSET_write(new_dset);
01524   PLUTO_add_dset(plint, new_dset, DSET_ACTION_MAKE_CURRENT);
01525 #ifdef PERMTEST_DEBUG
01526   signal(SIGUSR1, handler);
01527 #endif
01528   PERMTEST_reset_logo(plint);
01529   return (char *)0;
01530   }
01531 
01532 
01533 DEFINE_PLUGIN_PROTOTYPE
01534 
01535 PLUGIN_interface *PLUGIN_init(int ncall)
01536   {
01537   PLUGIN_interface *plint;
01538   if(ncall > 0)
01539     return (PLUGIN_interface *)0;       /*only one interface*/
01540 /*set titles and entry point*/
01541   plint = PLUTO_new_interface("Permutation Test", hint, help, PLUGIN_CALL_VIA_MENU, PERMTEST_main);
01542   PLUTO_add_hint(plint, hint);
01543 /*first line of dialogue box: input dataset*/
01544   PLUTO_add_option(plint, input_label, input_label, TRUE);
01545   PLUTO_add_dataset(plint, "Dataset", ANAT_SPGR_MASK | ANAT_EPI_MASK, 0, DIMEN_4D_MASK | BRICK_SHORT_MASK);
01546 /*second line of dialogue box: input time series*/
01547   PLUTO_add_option(plint, ts_label, ts_label, TRUE);
01548   PLUTO_add_timeseries(plint, "Reference Time Series");
01549 /*third line of dialogue box: time series against which to orthogonalise*/
01550   PLUTO_add_option(plint, ort_label, ort_label, FALSE);
01551   PLUTO_add_timeseries(plint, "Orthogonalisation Time Series");
01552 /*fourth line of dialogue box: output dataset*/
01553   PLUTO_add_option(plint, output_label, output_label, TRUE);
01554   PLUTO_add_string(plint, "Prefix", 0, (char **)0, 19);
01555 /*fifth line of dialogue box: alpha level (range 10^-4..1, default 0.05)*/
01556   PLUTO_add_option(plint, alpha_label, alpha_label, TRUE);
01557   PLUTO_add_number(plint, "alpha level", 1, 10000, 4, 500, 1);
01558 /*penultimate lines of dialogue box: tail options*/
01559   PLUTO_add_option(plint, tails2, tails2, FALSE);
01560   PLUTO_add_option(plint, tails1pos, tails1pos, FALSE);
01561   PLUTO_add_option(plint, tails1neg, tails1neg, FALSE);
01562 /*last lines of dialogue box: mask dataset and bounds*/
01563   PLUTO_add_option(plint, mask_label, mask_label, FALSE);
01564   PLUTO_add_dataset(plint, "mask dataset", 0, FUNC_FIM_MASK, DIMEN_3D_MASK | BRICK_SHORT_MASK | BRICK_BYTE_MASK);
01565   PLUTO_add_option(plint, masklo_label, masklo_label, FALSE);
01566   PLUTO_add_number(plint, "voxel is masked if >=", 0, 0x7fff, 0, 1, 1);
01567   PLUTO_add_option(plint, maskhi_label, maskhi_label, FALSE);
01568   PLUTO_add_number(plint, "voxel is masked if <=", 0, 0x7fff, 0, 1, 1);
01569   return plint;
01570   }