(((p) >= 0.5)? THD_pval_to_stat(2.0*(1.0-(p)), FUNC_ZT_TYPE, (float *)0) \
                : -THD_pval_to_stat(2.0*(p), FUNC_ZT_TYPE, (float *)0))

00169   {
00170   *x = coords % xdim;
00171   coords /= xdim;
00172   *y = coords % ydim;
00173   *z = coords / ydim;
00174   }

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   }

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   }

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   }

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   }

00505   {
00506   free(t->mem);
00507   free(t);
00508   }

00558   {
00559   t->next->data = clist;
00560   t->next->dlen = num_coords;
00561   dtree_insert_at_node(t, t->next++);
00562   }

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   }

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   }

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   }

00392   {
00393   return(r);
00394   }

00397   {
00398   return(-r);
00399   }

00837   {
00838   while(n--)
00839     if(x[n] < BLAST)
00840       new_y[n] = y[n]-(m*x[n]+b);
00841   }

00851                              :*/
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   }

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   }

01364   {
01365   XtVaSetValues(plint->im3d->vwid->picture, XmNlabelPixmap, XmUNSPECIFIED_PIXMAP, NULL);
01366   XFreePixmap(XtDisplay(plint->im3d->vwid->top_shell), sipb_pixmap);
01367   }

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   }

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   }

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   }

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   }

00425   {
00426   free(t->mem);
00427   free(t);
00428   }

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   }

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   }