00110   {
00111   *x = coords % xdim;
00112   coords /= xdim;
00113   *y = coords % ydim;
00114   *z = coords / ydim;
00115   }

00137   {
00138   register node *n;
00139   n = (node *)malloc(sizeof(node));
00140   n->value = v;
00141   n->ref_count = rc;
00142   n->subtree_ref_count = src;
00143   n->parent = p;
00144   n->left_child = l;
00145   n->right_child = r;
00146   return(n);
00147   }

00232   {
00233   register node *n;
00234   if((t->tail+1)%(THRESH_FILTER_LEN+1) == t->head)
00235     {
00236     fprintf(stderr, "delete_oldest: queue is empty!\n");
00237     return;
00238     }
00239   for(n = t->queue[t->head]->parent; n != NULLTREE; n = n->parent)
00240     n->subtree_ref_count--;
00241   n = t->queue[t->head];
00242   t->head = (t->head+1)%(THRESH_FILTER_LEN+1);
00243   if(n->ref_count == 1)
00244     prune(t, n);
00245   else
00246     {
00247     n->ref_count--;
00248     n->subtree_ref_count--;
00249     }
00250   }

00151   {
00152   free(n);
00153   }

00157   {
00158   register node *temp;
00159   while(n != NULLTREE)
00160     {
00161     destroy_tree(n->left_child);
00162     temp = n;
00163     n = n->right_child;
00164     destroy_node(temp);
00165     }
00166   }

00255   {
00256   register node *n;
00257   register int left_count, right_count, left_size, right_size, middle_position;
00258   if((t->tail+1)%(THRESH_FILTER_LEN+1) == t->head)
00259     {
00260     fprintf(stderr, "extract_median: queue is empty!\n");
00261     return(0);
00262     }
00263   n = t->tree;
00264   middle_position = n->subtree_ref_count/2+1;
00265   left_count = right_count = 0;
00266   left_size = (n->left_child==NULLTREE? 0: n->left_child->subtree_ref_count);
00267   right_size = (n->right_child==NULLTREE? 0: n->right_child->subtree_ref_count);
00268   while(abs(left_count+left_size - (right_count+right_size)) > n->ref_count)
00269   /*inv: 'left_count' is the number of values less than the median that have
00270     been excluded during traversal of the path from 't->tree' to 'n', and
00271     'left_size' is the size of the left subtree of the current node 'n'.
00272     'right_count' and 'right_size' are similar.*/
00273     {
00274     if(left_count+left_size+n->ref_count >= middle_position)
00275       {
00276       right_count += n->ref_count+right_size;
00277       n = n->left_child;
00278       }
00279     else
00280       {
00281       left_count += n->ref_count+left_size;
00282       n = n->right_child;
00283       }
00284     left_size = (n->left_child==NULLTREE? 0: n->left_child->subtree_ref_count);
00285     right_size = (n->right_child==NULLTREE? 0: n->right_child->subtree_ref_count);
00286     }
00287   return(n->value);
00288   }

00293   {
00294   register node *n, *p;
00295   if((t->tail+2)%(THRESH_FILTER_LEN+1) == t->head)
00296     {
00297     fprintf(stderr, "insert_newest: queue is full; deleting oldest to make room\n");
00298     delete_oldest(t);
00299     }
00300   t->tail = (t->tail+1)%(THRESH_FILTER_LEN+1);
00301   p = NULLTREE;
00302   n = t->tree;
00303   while((n != NULLTREE) && (n->value != v))
00304   /*inv: 'p' is the parent of 'n'.  All 'subtree_ref_count' fields on the path
00305     from 't->tree' to 'p' have been incremented.  The proper location for the
00306     new value 'v' lies somewhere in the subtree rooted at 'n'.*/
00307     {
00308     n->subtree_ref_count++;
00309     p = n;
00310     n = (v<n->value? n->left_child: n->right_child);
00311     }
00312   if(n == NULLTREE)
00313     {
00314     register node **graft_site;
00315     graft_site = (p==NULLTREE? &(t->tree):
00316       v<p->value? &(p->left_child): &(p->right_child)
00317       );
00318     *graft_site = create_node(v, 1, 1, p, NULLTREE, NULLTREE);
00319     t->queue[t->tail] = *graft_site;
00320     }
00321   else
00322     {
00323     n->ref_count++;
00324     n->subtree_ref_count++;
00325     t->queue[t->tail] = n;
00326     }
00327   }

00618   {
00619   PLUGIN_interface *plint;
00620   if(ncall > 0)
00621     return (PLUGIN_interface *)0;       /*only one interface*/
00622 /*set titles and entry point*/
00623   plint = PLUTO_new_interface("Threshold", hint, help, PLUGIN_CALL_VIA_MENU, THRESH_main);
00624   PLUTO_add_hint(plint, hint);
00625 /*first line of dialogue box: input dataset*/
00626   PLUTO_add_option(plint, input_label, input_label, TRUE);
00627   PLUTO_add_dataset(plint, "Dataset", ANAT_SPGR_MASK | ANAT_EPI_MASK, 0, DIMEN_4D_MASK | BRICK_SHORT_MASK);
00628 /*second line of dialogue box: output dataset*/
00629   PLUTO_add_option(plint, output_label, output_label, TRUE);
00630   PLUTO_add_string(plint, "Prefix", 0, (char **)0, 19);
00631   return plint;
00632   }

00171   {
00172   register node **prune_site, *largest;
00173   register int ref_count_of_largest;
00174   prune_site = (n->parent==NULLTREE? &(t->tree): n==n->parent->left_child? &(n->parent->left_child): &(n->parent->right_child));
00175   if(n->left_child == NULLTREE)
00176     {
00177     *prune_site = n->right_child;
00178     if(*prune_site != NULLTREE)
00179       (*prune_site)->parent = n->parent;
00180     destroy_node(n);
00181     }
00182   else if(n->right_child == NULLTREE)
00183     {
00184     *prune_site = n->left_child;
00185     if(*prune_site != NULLTREE)
00186       (*prune_site)->parent = n->parent;
00187     destroy_node(n);
00188     }
00189   else
00190     {
00191   /*find the largest value in the left subtree of 'n'*/
00192     for(largest = n->left_child; largest->right_child != NULLTREE; largest = largest->right_child)
00193       ;
00194   /*adjust reference counts to reflect the pruning of this largest value*/
00195     ref_count_of_largest = largest->ref_count;
00196     for(largest = n->left_child; largest->right_child != NULLTREE; largest = largest->right_child)
00197       largest->subtree_ref_count -= ref_count_of_largest;
00198   /*prune the largest value by replacing it with its left subtree*/
00199     if(largest==largest->parent->left_child)
00200       {
00201       largest->parent->left_child = largest->left_child;
00202       if(largest->parent->left_child != NULLTREE)
00203         largest->parent->left_child->parent = largest->parent;
00204       }
00205     else
00206       {
00207       largest->parent->right_child = largest->left_child;
00208       if(largest->parent->right_child != NULLTREE)
00209         largest->parent->right_child->parent = largest->parent;
00210       }
00211   /*substitute this largest-valued node for node 'n'*/
00212     if(n->parent == NULLTREE)
00213       t->tree = largest;
00214     else if(n == n->parent->left_child)
00215       n->parent->left_child = largest;
00216     else
00217       n->parent->right_child = largest;
00218     largest->parent = n->parent;
00219     largest->left_child = n->left_child;
00220     largest->right_child = n->right_child;
00221     if(largest->left_child != NULLTREE)
00222       largest->left_child->parent = largest;
00223     if(largest->right_child != NULLTREE)
00224       largest->right_child->parent = largest;
00225     largest->subtree_ref_count = largest->ref_count + (largest->left_child==NULLTREE? 0: largest->left_child->subtree_ref_count) + (largest->right_child==NULLTREE? 0: largest->right_child->subtree_ref_count);
00226     destroy_node(n);
00227     }
00228   }

00457   {
00458   register int t, x, y, z;
00459   int xdim, ydim, zdim, nvox, tdim, histo_min, histo_max, centroid;
00460   double centroid_num, centroid_denom, ln;
00461   long sum;
00462   short *img;
00463   short air_peak, brain_peak, cutoff;
00464   btree filter;
00465   int histogram[1+THRESH_MAX_BRAIN];
00466   xdim = dset->daxes->nxx;
00467   ydim = dset->daxes->nyy;
00468   zdim = dset->daxes->nzz;
00469   nvox = xdim*ydim*zdim;
00470   tdim = DSET_NUM_TIMES(dset);
00471   bzero((char *)histogram, (1+THRESH_MAX_BRAIN)*sizeof(*histogram));
00472   img = (short *)calloc(nvox, sizeof(short));
00473   if(img == (short *)0)
00474     return((short *)0);
00475   /*Form an image each of whose voxels is the mean of the time series for the
00476     corresponding voxel in the original image.  Construct a histogram of these
00477     mean intensities.*/
00478   for(z = 0; z != zdim; z++)
00479     for(y = 0; y != ydim; y++)
00480       for(x = 0; x != xdim; x++)
00481         {
00482         sum = 0L;
00483         for (t = 2; t < tdim; t++) /*t=2 to stabilise transverse magnetisation*/
00484           sum += ((short *)DSET_ARRAY(dset, t))[x + xdim*(y + ydim*z)];
00485         sum = (sum+(tdim-2)/2)/(tdim-2);
00486         img[coord_encode(x, y, z, xdim, ydim)] = (short)sum;
00487         if((sum >= 0) && (sum <= THRESH_MAX_BRAIN))
00488           histogram[sum]++;
00489         }
00490   centroid_num = centroid_denom = 0.0;
00491   for(x = THRESH_MAX_BRAIN; x != 0; x--)
00492     {
00493     ln = log((double)(1+histogram[x]));
00494     centroid_num += x*ln;
00495     centroid_denom += ln;
00496     }
00497   centroid = (int)(centroid_num/centroid_denom);
00498 
00499   filter.tree = NULLTREE;
00500   filter.head = 0;
00501   filter.tail = -1;
00502   x = THRESH_MAX_BRAIN;
00503   while(x > THRESH_MAX_BRAIN-THRESH_FILTER_LEN)
00504     insert_newest(histogram[--x], &filter);
00505   histo_max = -1;
00506   histo_min = MAXINT;
00507   cutoff = brain_peak = THRESH_MAX_BRAIN;
00508   /*inv: filter contains histogram[x..x+THRESH_FILTER_LEN-1], histo_max is the
00509     maximum median-filtered value of histogram[x+1+THRESH_FILTER_LEN/2..
00510     THRESH_MAX_BRAIN-(THRESH_FILTER_LEN+1)/2, brain_peak is the index at which
00511     histo_max occurs, histo_min is the minimum median-filtered value of
00512     histogram[x+1+THRESH_FILTER_LEN/2..brain_peak], and cutoff is the index at
00513     which histo_min occurs.*/
00514   while(x > centroid-THRESH_FILTER_LEN/2)
00515     {
00516     y = extract_median(&filter);
00517     if(y > histo_max)
00518       {
00519       cutoff = brain_peak = x+THRESH_FILTER_LEN/2;
00520       histo_min = histo_max = y;
00521       }
00522     else if(y < histo_min)
00523       {
00524       cutoff = x+THRESH_FILTER_LEN/2;
00525       histo_min = y;
00526       }
00527     delete_oldest(&filter);
00528     insert_newest(histogram[--x], &filter);
00529     }
00530   /*inv: filter contains histogram[x..x+THRESH_FILTER_LEN-1], histo_min is the
00531     minimum median-filtered value of histogram[x+1+THRESH_FILTER_LEN/2..
00532     brain_peak], and cutoff is the index at which histo_min occurs.*/
00533   while((x >= 0) && ((y = extract_median(&filter)) <= brain_peak))
00534     {
00535     if(y < histo_min)
00536       {
00537       histo_min = y;
00538       cutoff = x+THRESH_FILTER_LEN/2;
00539       }
00540     delete_oldest(&filter);
00541     insert_newest(histogram[--x], &filter);
00542     }
00543   for(z = cutoff-THRESH_FILTER_LEN/2, y = z+THRESH_FILTER_LEN; z < y; z++)
00544     if(histogram[z] < histo_min)
00545       {
00546       histo_min = histogram[z];
00547       cutoff = z;
00548       }
00549   if(verbose)
00550     printf("centroid %d, brain peak %d, air peak edge %d, threshold %d\n", centroid, brain_peak, x+THRESH_FILTER_LEN/2, cutoff);
00551   destroy_tree(filter.tree);
00552   free(histogram);
00553   /*region-grow from the edge of the image inward, filling with zeroes*/
00554   if(THRESH_mask_brain(img, xdim, ydim, zdim, cutoff))
00555     {
00556     free(img);
00557     return((short *)0);
00558     }
00559   return(img);
00560   }

00563   {
00564   int t;
00565   char *prefix;
00566   THD_3dim_dataset *dset, *mask;
00567   short *mask_img;
00568   if(plint == (PLUGIN_interface *)0)
00569     return "THRESH_main: null input";
00570 
00571 /*Make sure that the input dataset exists and is stored in a format that we can
00572   understand (MRI_short)*/
00573   PLUTO_next_option(plint);
00574   dset = PLUTO_find_dset(PLUTO_get_idcode(plint));
00575   if(dset == (THD_3dim_dataset *)0)
00576     return "bad dataset";
00577   for(t = 0; t != dset->dblk->nvals; t++)
00578     if(DSET_BRICK_TYPE(dset, t) != MRI_short)
00579       return("thresholding on non-short values is not implemented");
00580 
00581 /*Make sure that the prefix specified for the new dataset is a valid prefix*/
00582   PLUTO_next_option(plint);
00583   prefix = PLUTO_get_string(plint);
00584   if(!PLUTO_prefix_ok(prefix))
00585     return("bad prefix");
00586 
00587 /*Make sure source dataset is in memory*/
00588   DSET_load(dset);
00589   mask_img = THRESH_compute(dset, 1);
00590   if(mask_img == (short *)0)
00591     return("out of memory");
00592 
00593 /*create the output dataset*/
00594   mask = EDIT_empty_copy(dset);
00595   if(EDIT_dset_items(mask,
00596         ADN_prefix, prefix,
00597         ADN_malloc_type, DATABLOCK_MEM_MALLOC, /*hold in r/w memory*/
00598         ADN_datum_all, MRI_short,       /*store as (scaled) short ints*/
00599         ADN_nvals, 1,                   /*single sub-brick*/
00600         ADN_ntt, 0,                     /*no time dimension*/
00601         ADN_type, ISHEAD(dset)?
00602           HEAD_FUNC_TYPE: GEN_FUNC_TYPE, /*functional image*/
00603         ADN_func_type, FUNC_FIM_TYPE,   /*intensity + z-score*/
00604         ADN_none))
00605     return("EDIT_dset_items error");
00606   EDIT_BRICK_LABEL(mask, 0, "Mask");
00607   mri_fix_data_pointer(mask_img, DSET_BRICK(mask, 0));
00608   EDIT_BRICK_FACTOR(mask, 0, 0.0);      /*no scaling*/
00609   DSET_write(mask);
00610   PLUTO_add_dset(plint, mask, DSET_ACTION_NONE);
00611   return (char *)0;
00612   }

00404   {
00405   register int x, y, z, region;
00406   int region_size, max_region_size, max_region;
00407   short *mask_img;
00408   int *stack;
00409   mask_img = calloc(xdim*ydim*zdim, sizeof(*mask_img));
00410   if(mask_img == (short *)0)
00411     return 1;
00412   stack = (int *)calloc(xdim*ydim*zdim, sizeof(*stack));
00413   if(stack == (int *)0)
00414     {
00415     free(mask_img);
00416     return 1;
00417     }
00418 /*make all non-brain regions in mask_img INCLUDED, and leave all brain regions
00419   UNVISITED*/
00420   THRESH_region_grow(img, mask_img, stack, 0, 0, 0, 0, xdim, ydim, zdim, threshold);
00421   bzero((char *)img, xdim*ydim*zdim*sizeof(*img)); /*fill img with UNVISITED*/
00422 /*now label all the volumes that were not INCLUDED in mask_img*/
00423   region = 1;
00424   max_region_size = 0;
00425   for(z = 0; z != zdim; z++)
00426     for(y = 0; y != ydim; y++)
00427       for(x = 0; x != xdim; x++)
00428         if(img[coord_encode(x, y, z, xdim, ydim)] == UNVISITED)
00429           {
00430           region_size = THRESH_region_grow(mask_img, img, stack, region, x, y, z, xdim, ydim, zdim, INCLUDED);
00431           if(region_size > max_region_size)
00432             {
00433             max_region_size = region_size;
00434             max_region = region;
00435             }
00436           region++;
00437           }
00438 /*zero all but the largest region, which is assumed to be the brain*/
00439   for(z = 0; z != xdim*ydim*zdim; z++)
00440     img[z] = (img[z] == INCLUDED+max_region);
00441   free(stack);
00442   free(mask_img);
00443   return 0;
00444   }

00358   {
00359   register int dx, dy, dz;      /*increments to current voxel index*/
00360   register int  *sp,            /*stack of voxels scheduled for visiting*/
00361                 *sp2,           /*stack of excluded voxels to unlabel*/
00362                 region_size;    /*count of included voxels*/
00363   region_size = 0;
00364   sp = stack;
00365   sp2 = stack+xdim*ydim*zdim;
00366   *sp = coord_encode(x, y, z, xdim, ydim);
00367   mask_img[*(sp++)] = QUEUED;
00368 /*inv: the stack contains all reachable points that still need to be visited*/
00369   while(sp != stack)
00370     if(mask_img[*(--sp)] == QUEUED)
00371       {
00372       if(img[*sp] < threshold)
00373         {
00374         mask_img[*sp] = INCLUDED+region_num;
00375         region_size++;
00376         coord_decode(*sp, &x, &y, &z, xdim, ydim);
00377         for(dz = -1; dz <= 1; dz++)
00378           for(dy = -1; dy <= 1; dy++)
00379             for(dx = -1; dx <= 1; dx++)
00380               if((x+dx >= 0) && (x+dx < xdim) && (y+dy >= 0) && (y+dy < ydim) && (z+dz >= 0) && (z+dz < zdim) && (dx || dy || dz))
00381                 {
00382                 *sp = coord_encode(x+dx, y+dy, z+dz, xdim, ydim);
00383                 if(mask_img[*sp] == UNVISITED)
00384                   mask_img[*(sp++)] = QUEUED;
00385                 }
00386         }
00387       else
00388         {
00389         mask_img[*sp] = EXCLUDED;
00390         *(--sp2) = *sp;
00391         }
00392       }
00393 /*remove all EXCLUDED labels and return these voxels to the UNVISITED state*/
00394   while(sp2 != stack+xdim*ydim*zdim)
00395     mask_img[*(sp2++)] = UNVISITED;
00396 /*the visited areas of mask_img now consist entirely of UNVISITED or INCLUDED*/
00397   return(region_size);
00398   }