00001 /*plug_threshold.c - AFNI plugin that creates a mask separating brain from
00002   non-brain, and saves this mask as a fim dataset.
00003   Copyright (c) 2000 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 /*Most scanners, including GE Signa and Siemens Vision systems, do not use the
00024   full 16-bit dynamic range.  THRESH_MAX_BRAIN can therefore be decreased to
00025   about 3000 to save memory, if desired.  It's set to 32766 here rather than to
00026   32767 since we want to be able to add one to it without causing a signed
00027   overflow.*/
00028 #define THRESH_MAX_BRAIN 0x7ffe
00029 
00030 /*A median-filter length of 9 seems to work for every data set that I've tested.
00031   If your scanner gives a banded distribution of voxel intensities instead of a
00032   uniform distribution, you may need to increase this value or apply some
00033   smoothing to the histogram.*/
00034 #define THRESH_FILTER_LEN 9
00035 
00036 #include <sys/types.h>
00037 #include <stdio.h>
00038 #include <stdlib.h>
00039 #if !(defined(DARWIN) || defined(FreeBSD))
00040 # include <malloc.h>
00041 #endif
00042 #include <strings.h>
00043 #include <math.h>
00044 
00045 #if !(defined(DARWIN) || defined(CYGWIN) || defined(FreeBSD)) /* RWCox patch */
00046 #include <values.h>
00047 #endif
00048 
00049 #ifndef MAXINT
00050 #define MAXINT (1<<30)
00051 #endif
00052 
00053 #include "afni.h"
00054 
00055 static char help[] =
00056   "This plugin generates a mask dataset that labels brain with the value 1\n"
00057   "and non-brain with the value 0.  The algorithm operates in two phases,\n"
00058   "first finding the best threshold intensity to separate brain from\n"
00059   "non-brain, then applying that threshold along with information on\n"
00060   "connected regions to generate the mask.\n\n"
00061 
00062   "The first phase of the algorithm begins by forming an image each of whose\n"
00063   "voxels is the mean of the time series for the corresponding voxel in the\n"
00064   "echo-planar dataset.  (The first two acquisitions are excluded from the\n"
00065   "computation of these means, in order to allow the tissue to reach a\n"
00066   "magnetic steady state.)  The program constructs a histogram of these\n"
00067   "means, and then computes the centroid of the log-transformed histogram.\n"
00068   "This centroid value is assumed to lie somewhere in the interval between\n"
00069   "the non-brain peak to the left and the brain peak to the right.  (This\n"
00070   "assumption holds as long as the field of view is appropriate for the size\n"
00071   "of the subject's head -- i.e., as long as the air peak doesn't dominate\n"
00072   "the log-transformed histogram.)  The maximum of the median-filtered\n"
00073   "histogram in the interval to the right of the centroid value is the brain\n"
00074   "peak.  The maximum to the left of the centroid value is the air peak.\n"
00075   "The minimum between these two peaks is the threshold between air, muscle,\n"
00076   "and bone to the left and brain to the right.\n\n"
00077 
00078   "The second phase of the algorithm region-grows from a corner of the image,\n"
00079   "stopping at voxels whose intensities exceed the threshold that was\n"
00080   "computed in the previous phase.  The set of voxels not touched by this\n"
00081   "region-growing step necessarily includes the brain, although it generally\n"
00082   "also includes other high-intensity structures such as the eyes and perhaps\n"
00083   "some smaller clumps of tissue elsewhere in the head.  To get rid of these\n"
00084   "smaller, non-brain islets, we region-grow from each untouched voxel, and\n"
00085   "identify the contiguous region whose volume is greatest.  This region is\n"
00086   "labelled as brain, and everything outside it is labelled as non-brain.\n\n"
00087 
00088 "AUTHOR\n\n"
00089 
00090   "This plugin was written by Matthew Belmonte <belmonte@mit.edu> of the\n"
00091   "MIT Student Information Processing Board, supported by a grant from the\n"
00092   "National Alliance for Autism Research.\n\n"
00093 
00094 "VERSION\n\n"
00095 
00096   "1.1   (14 June 2001)\n\n"
00097 
00098 "SEE ALSO\n\n"
00099 
00100   "Draw Dataset Plugin",
00101 
00102             hint[] = "mask out non-brain",
00103             input_label[] = "Input",
00104             output_label[] = "Output";
00105 
00106 /*encoding and decoding functions for mapping a triple to a single integer*/
00107 #define coord_encode(x, y, z, xdim, ydim) ((x) + (xdim)*((y) + (ydim)*(z)))
00108 
00109 static void coord_decode(int coords, int *x, int *y, int *z, int xdim, int ydim)
00110   {
00111   *x = coords % xdim;
00112   coords /= xdim;
00113   *y = coords % ydim;
00114   *z = coords / ydim;
00115   }
00116 
00117 /**********************************************************************
00118  * MEDIAN FILTERING ROUTINES                                          *
00119  * Adapted from median.c in the Gnuroscan package by Matthew Belmonte *
00120  **********************************************************************/
00121 
00122 typedef struct _node {
00123   int value;
00124   int ref_count, subtree_ref_count;
00125   struct _node *parent, *left_child, *right_child;
00126   } node;
00127 #define NULLTREE ((node *)0)
00128 
00129 typedef struct {
00130   node *tree;
00131   int head, tail;
00132   node *(queue[1+THRESH_FILTER_LEN]);
00133   } btree;
00134 
00135 /*Allocate a node whose fields contain the given values.*/
00136 static node *create_node(int v, int rc, int src, node *p, node *l, node *r)
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   }
00148 
00149 /*Deallocate a node.*/
00150 static void destroy_node(node *n)
00151   {
00152   free(n);
00153   }
00154 
00155 /*Deallocate all the nodes of the tree rooted at the given node 'n'.*/
00156 static void destroy_tree(node *n)
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   }
00167 
00168 /*Prune node 'n' from the tree rooted at 't->tree', rearranging the tree as
00169   necessary in order to preserve the descendants of 'n'.*/
00170 static void prune(btree *t, node *n)
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   }
00229 
00230 /*Delete from the given btree the node at the head of the associated queue.*/
00231 static void delete_oldest(btree *t)
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   }
00251 
00252 /*Return the median of all the values in the given btree.  Do not alter the
00253   btree.*/
00254 static int extract_median(btree *t)
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   }
00289 
00290 /*Insert the given value into the given btree and place it at the tail of the
00291   associated queue.*/
00292 static void insert_newest(int v, btree *t)
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   }
00328 
00329 /***************************
00330  * REGION-GROWING ROUTINES *
00331  ***************************/
00332 
00333 #define UNVISITED 0     /*not yet evaluated or queued*/
00334 #define INCLUDED 1      /*included in target region*/
00335 #define EXCLUDED 2      /*excluded from target region*/
00336 #define QUEUED 3        /*queued for evaluation*/
00337 /*Find the set of voxels in img whose intensities are less than the given
00338   threshold and that are contiguous with the given starting voxel (x,y,z).
00339   Label the corresponding voxels in mask_img with the value INCLUDED+region_num.
00340   (mask_img should have been initialised to UNVISITED by the calling routine.
00341   With UNVISITED defined as 0, calloc() or bzero() suffices to accomplish this
00342   initialisation.  Multiple calls can be made without intervening
00343   re-initialisations only in the case in which the regions being labelled are
00344   unconnected.)  This routine returns the size in voxels of the labelled
00345   region.*/
00346 static int THRESH_region_grow(
00347 short *img, short *mask_img, int *stack, int region_num,
00348 int x, int y, int z, int xdim, int ydim, int zdim, short threshold)
00349 #if 0
00350 short   *img,           /*source image*/
00351         *mask_img;      /*same dimensions as img, pre-initialised to UNVISITED*/
00352 int     *stack,         /*a block of xdim*ydim*zdim ints*/
00353         region_num,     /*unique number with which to label this region*/
00354         x, y, z,        /*starting location of this region*/
00355         xdim, ydim, zdim; /*dimensions*/
00356 short   threshold;      /*threshold below which pixels are included*/
00357 #endif
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   }
00399 
00400 /*On entry, img is the averaged echo-planar image.  On exit, each voxel in img
00401   has been replaced by a 1 if it's been identified as a brain voxel, or by a 0
00402   otherwise.  Returns 0 if OK, 1 if out of memory.*/
00403 static int THRESH_mask_brain(short *img, int xdim, int ydim, int zdim, short threshold)
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   }
00445 
00446 /*****************
00447  * MAIN ROUTINES *
00448  *****************/
00449 
00450 /*Scan from the right-hand end of the histogram until the median-filtered value
00451   of the histogram definitely exceeds zero (right-hand edge of the brain peak).
00452   Then continue scanning left, keeping track of the filtered histogram's
00453   minimum value, till its current value becomes very large (right-hand edge of
00454   the air peak).  The position of this minimum is the intensity threshold for
00455   distinguishing brain from non-brain.*/
00456 static short *THRESH_compute(THD_3dim_dataset *dset, int verbose)
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   }
00561 
00562 char *THRESH_main(PLUGIN_interface *plint)
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   }
00613 
00614 
00615 DEFINE_PLUGIN_PROTOTYPE
00616 
00617 PLUGIN_interface *PLUGIN_init(int ncall)
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   }