00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #define THRESH_MAX_BRAIN 0x7ffe
00029
00030
00031
00032
00033
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))
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
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
00119
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
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
00150 static void destroy_node(node *n)
00151 {
00152 free(n);
00153 }
00154
00155
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
00169
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
00192 for(largest = n->left_child; largest->right_child != NULLTREE; largest = largest->right_child)
00193 ;
00194
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
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
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
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
00253
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
00270
00271
00272
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
00291
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
00305
00306
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
00331
00332
00333 #define UNVISITED 0
00334 #define INCLUDED 1
00335 #define EXCLUDED 2
00336 #define QUEUED 3
00337
00338
00339
00340
00341
00342
00343
00344
00345
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,
00351 *mask_img;
00352 int *stack,
00353 region_num,
00354 x, y, z,
00355 xdim, ydim, zdim;
00356 short threshold;
00357 #endif
00358 {
00359 register int dx, dy, dz;
00360 register int *sp,
00361 *sp2,
00362 region_size;
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
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
00394 while(sp2 != stack+xdim*ydim*zdim)
00395 mask_img[*(sp2++)] = UNVISITED;
00396
00397 return(region_size);
00398 }
00399
00400
00401
00402
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
00419
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));
00422
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
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
00448
00449
00450
00451
00452
00453
00454
00455
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
00476
00477
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++)
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
00509
00510
00511
00512
00513
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
00531
00532
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
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
00572
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
00582 PLUTO_next_option(plint);
00583 prefix = PLUTO_get_string(plint);
00584 if(!PLUTO_prefix_ok(prefix))
00585 return("bad prefix");
00586
00587
00588 DSET_load(dset);
00589 mask_img = THRESH_compute(dset, 1);
00590 if(mask_img == (short *)0)
00591 return("out of memory");
00592
00593
00594 mask = EDIT_empty_copy(dset);
00595 if(EDIT_dset_items(mask,
00596 ADN_prefix, prefix,
00597 ADN_malloc_type, DATABLOCK_MEM_MALLOC,
00598 ADN_datum_all, MRI_short,
00599 ADN_nvals, 1,
00600 ADN_ntt, 0,
00601 ADN_type, ISHEAD(dset)?
00602 HEAD_FUNC_TYPE: GEN_FUNC_TYPE,
00603 ADN_func_type, FUNC_FIM_TYPE,
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);
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;
00622
00623 plint = PLUTO_new_interface("Threshold", hint, help, PLUGIN_CALL_VIA_MENU, THRESH_main);
00624 PLUTO_add_hint(plint, hint);
00625
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
00629 PLUTO_add_option(plint, output_label, output_label, TRUE);
00630 PLUTO_add_string(plint, "Prefix", 0, (char **)0, 19);
00631 return plint;
00632 }