00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 2001, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   This program finds local extrema (minima or maxima) of the input
00010   dataset values for each sub-brick of the input dataset.  The extrema
00011   may be determined either for each volume, or for each individual slice.  
00012   Only those voxels whose corresponding intensity value is greater than 
00013   the user specified data threshold will be considered.
00014 
00015   File:    3dExtrema.c
00016   Author:  B. Douglas Ward
00017   Date:    12 April 2001
00018 
00019   Mod:     Added call to AFNI_logger.
00020   Date:    15 August 2001
00021 
00022 */
00023 /*---------------------------------------------------------------------------*/
00024 
00025 #define PROGRAM_NAME "3dExtrema"                     /* name of this program */
00026 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
00027 #define PROGRAM_INITIAL "12 April 2001" /* date of initial program release */
00028 #define PROGRAM_LATEST "15 August 2001"     /* date of last program revision */
00029 
00030 
00031 /*---------------------------------------------------------------------------*/
00032 /*
00033   Include header files and source code files.
00034 */
00035 
00036 #include "mrilib.h"
00037 
00038 
00039 /*---------------------------------------------------------------------------*/
00040 /*
00041   Structure declarations 
00042 */
00043 
00044 struct extrema;
00045 
00046 typedef struct extrema
00047 {
00048   float intensity;                     /* (averaged) extrema intensity */
00049   float * centroid;                    /* (averaged) extrema location */
00050   int count;                           /* count of extrema in average */
00051   float sum;                           /* sum of extrema intensities */
00052   float nearest_dist;                  /* distance to nearest extrema */
00053   struct extrema * nearest_extrema;    /* pointer to nearest extrema */
00054   struct extrema * next_extrema;       /* pointer to next extrema in list */
00055 } extrema;
00056 
00057 
00058 /*-------------------------- global data ------------------------------------*/
00059 
00060 #define  EX_MAX_LL 5000      /* maximum number of linked lists of extrema */
00061 #define  EX_MAX_NV 26        /* maximum number of neighboring voxels */
00062 #define  EX_DIMENSION 3      /* 3 dimensional space */
00063 
00064 #define  EX_GT 1             /* indices for binary relations */
00065 #define  EX_GE 2
00066 #define  EX_LT 3
00067 #define  EX_LE 4
00068 
00069 static THD_coorder EX_cord;  /* get orientation from AFNI_ORIENT environment */
00070 
00071 static int EX_nx, EX_ny, EX_nz,           /* dataset dimensions in voxels */
00072            EX_nxy, EX_nxyz;               
00073 
00074 static int        EX_quiet      = 0;      /* flag for suppress screen output */
00075 static int        EX_relation   = 1;      /* flag for binary relation */
00076 static int        EX_maxima     = 1;      /* flag for maxima or minima */
00077 static int        EX_strict     = 1;      /* flag for strict or parial ineq. */
00078 static int        EX_interior   = 1;      /* flag for interior or closure */
00079 static int        EX_slice      = 1;      /* flag for slice or volume */
00080 static int        EX_sort       = 1;      /* flag for sort extrema */
00081 static int        EX_merge      = 1;      /* flag for remove or average */
00082 static float      EX_sep_dist   = 0.0;    /* minimum separation distance */
00083 
00084 static float      EX_mask_thr   = 1.0;    /* mask threshold */
00085 static float      EX_data_thr   = 0.0;    /* data threshold */
00086 static int        EX_nbricks    = 0;      /* number of output subbricks */
00087 
00088 static extrema *  EX_head_extrema[EX_MAX_LL];    /* linked lists of extrema */
00089 static int        EX_num_ll = 0;                 /* number of linked lists */
00090 
00091 static int EX_offset[EX_MAX_NV];   /* relative indices of neighboring voxels */
00092 
00093 static char * EX_mask_filename    = NULL;
00094 static char * EX_output_prefix    = NULL;
00095 static char * EX_session          = "./";
00096 
00097 static byte * EX_mask = NULL;             /* mask for voxels above thr. */
00098 
00099 static char * commandline = NULL;          /* command line for history notes */
00100 
00101 static THD_3dim_dataset * EX_dset = NULL;  /* first input dataset */  
00102 
00103 /*---------------------------------------------------------------------------*/
00104 
00105 /** macro to test a malloc-ed pointer for validity **/
00106 
00107 #define MTEST(ptr) \
00108 if((ptr)==NULL) \
00109 ( printf ("Cannot allocate memory \n"),  exit(1) )
00110      
00111 
00112 /*---------------------------------------------------------------------------*/
00113 
00114 /** macro to open a dataset and make it ready for processing **/
00115 
00116 #define DOPEN(ds,name)                                                        \
00117 do{ int pv ; (ds) = THD_open_dataset((name)) ;                                \
00118        if( !ISVALID_3DIM_DATASET((ds)) ){                                     \
00119           fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \
00120        THD_load_datablock( (ds)->dblk ) ;                                     \
00121        pv = DSET_PRINCIPAL_VALUE((ds)) ;                                      \
00122        if( DSET_ARRAY((ds),pv) == NULL ){                                     \
00123           fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); }   \
00124        if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){                         \
00125           fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1);        }     \
00126        break ; } while (0)
00127 
00128 
00129 /*---------------------------------------------------------------------------*/
00130 
00131 /** macro to return pointer to correct location in brick for current processing **/
00132 
00133 #define SUB_POINTER(ds,vv,ind,ptr)                                            \
00134    do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){                                  \
00135          default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1);       \
00136             case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ;  \
00137                             (ptr) = (void *)( fim + (ind) ) ;                 \
00138             } break ;                                                         \
00139             case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ;     \
00140                             (ptr) = (void *)( fim + (ind) ) ;                 \
00141             } break ;                                                         \
00142             case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ;  \
00143                              (ptr) = (void *)( fim + (ind) ) ;                \
00144             } break ; } break ; } while(0)
00145 
00146 
00147 /*---------------------------------------------------------------------------*/
00148 /*
00149    Print error message and stop.
00150 */
00151 
00152 void EX_error (char * message)
00153 {
00154   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00155   exit(1);
00156 }
00157 
00158 
00159 /*---------------------------------------------------------------------------*/
00160 /*
00161   Create an empty extrema.
00162 */
00163   
00164 extrema * initialize_extrema ()
00165 {
00166   extrema * extrema_ptr = NULL;
00167 
00168   extrema_ptr = (extrema *) malloc (sizeof(extrema));
00169   MTEST (extrema_ptr);
00170   
00171   extrema_ptr->intensity = 0.0;
00172   extrema_ptr->centroid = NULL;
00173   extrema_ptr->count = 0;
00174   extrema_ptr->sum = 0.0;
00175   extrema_ptr->nearest_dist = 0.0;
00176   extrema_ptr->nearest_extrema = NULL;
00177   extrema_ptr->next_extrema = NULL;
00178 
00179   return (extrema_ptr);
00180   
00181 }
00182 
00183 
00184 /*---------------------------------------------------------------------------*/
00185 /*
00186   Print the contents of one extrema.
00187 */
00188 
00189 void print_extrema (int index, extrema * extrema_ptr)
00190 {
00191   int j;
00192   char str[30];
00193 
00194   sprintf (str, "%5d", index);
00195   printf ("%5s", str);
00196 
00197   printf ("%15.3f", extrema_ptr->intensity);
00198 
00199   for (j = 0;  j < EX_DIMENSION;  j++)
00200     printf ("%10.2f", extrema_ptr->centroid[j]);
00201 
00202 
00203   printf("%10d", extrema_ptr->count);
00204 
00205   if (extrema_ptr->nearest_dist < 1.0e+20)
00206     printf ("%10.3f", extrema_ptr->nearest_dist);
00207 
00208   printf ("\n");
00209 }
00210 
00211 
00212 /*---------------------------------------------------------------------------*/
00213 /*
00214   Print the contents of all extrema in the linked list.
00215 */
00216 
00217 void print_all_extrema (int ivolume, int islice, extrema * extrema_ptr)
00218 {
00219   int index;
00220 
00221   printf ("\n");
00222   if (EX_maxima)
00223     printf ("Maxima ");
00224   else
00225     printf ("Minima ");
00226 
00227   if (EX_slice)
00228     printf ("for Volume #%d and Slice #%d: \n", ivolume, islice);
00229   else
00230     printf ("for Volume #%d: \n", ivolume);
00231 
00232 
00233   if (extrema_ptr != NULL)
00234     {
00235       printf ("%5s",     "Index");
00236       printf ("%15s",    "Intensity");
00237       printf ("%6s[mm]", ORIENT_tinystr[EX_cord.xxor]);
00238       printf ("%6s[mm]", ORIENT_tinystr[EX_cord.yyor]);
00239       printf ("%6s[mm]", ORIENT_tinystr[EX_cord.zzor]);
00240       printf ("%10s",    "Count");
00241       printf ("%6s[mm]", "Dist");
00242       printf ("\n");
00243       
00244       printf ( "%5s", "-----");
00245       printf ("%15s", "---------");
00246       printf ("%10s", "------");
00247       printf ("%10s", "------");
00248       printf ("%10s", "------");
00249       printf ("%10s", "-----");
00250       printf ("%10s", "--------");
00251       printf ("\n");
00252     }
00253 
00254   index = 0;
00255   while (extrema_ptr != NULL)
00256     {
00257       index++;
00258       print_extrema (index, extrema_ptr);
00259       extrema_ptr = extrema_ptr->next_extrema;
00260     }
00261 }
00262 
00263 
00264 /*---------------------------------------------------------------------------*/
00265 /*
00266   Calculate the distance between two extrema as the Euclidean distance
00267   between their centroids.
00268 */
00269 
00270 float extrema_distance (extrema * aextrema, extrema * bextrema)
00271 {
00272   float sumsqr;
00273   float delta;
00274   int i;
00275 
00276   sumsqr = 0.0;
00277   for (i = 0;  i < EX_DIMENSION;  i++)
00278     {
00279       delta = aextrema->centroid[i] - bextrema->centroid[i];
00280       sumsqr += delta * delta;
00281     }
00282 
00283   return (sqrt(sumsqr));
00284   
00285 }
00286 
00287 
00288 /*---------------------------------------------------------------------------*/
00289 /*
00290   Find the extrema which is nearest to new_extrema.
00291   Set the nearest_dist and nearest_extrema structure elements accordingly.
00292 */
00293   
00294 void find_nearest_extrema (extrema * new_extrema, extrema * head_extrema)
00295 {
00296   const float MAX_DIST = 1.0e+30;
00297 
00298   extrema * extrema_ptr = NULL;
00299   float dist;
00300 
00301 
00302   /*----- Initialize nearest extrema elements -----*/
00303   new_extrema->nearest_dist = MAX_DIST;
00304   new_extrema->nearest_extrema = NULL;
00305 
00306 
00307   extrema_ptr = head_extrema;
00308 
00309   while (extrema_ptr != NULL)
00310     {
00311       if (extrema_ptr != new_extrema)
00312         {
00313           dist = extrema_distance (new_extrema, extrema_ptr);
00314 
00315           if (dist < new_extrema->nearest_dist)
00316             {
00317               new_extrema->nearest_dist = dist;
00318               new_extrema->nearest_extrema = extrema_ptr;
00319             }
00320 
00321           if (dist < extrema_ptr->nearest_dist)
00322             {
00323               extrema_ptr->nearest_dist = dist;
00324               extrema_ptr->nearest_extrema = new_extrema;
00325             }
00326         }
00327 
00328       extrema_ptr = extrema_ptr->next_extrema;
00329     }
00330 }
00331 
00332 
00333 /*---------------------------------------------------------------------------*/
00334 /*
00335   Add a new extrema to the linked list of extrema.
00336 */
00337 
00338 extrema * add_extrema (extrema * new_extrema, extrema * head_extrema)
00339 {
00340 
00341   new_extrema->next_extrema = head_extrema;
00342 
00343   find_nearest_extrema (new_extrema, head_extrema);
00344 
00345   return (new_extrema);
00346 }
00347 
00348 
00349 /*---------------------------------------------------------------------------*/
00350 /*
00351   Create a new extrema containing a single voxel, and add to list of extrema.
00352 */
00353 
00354 extrema * new_extrema (float * centroid, float intensity,
00355                        extrema * head_extrema)
00356 {
00357   extrema * extrema_ptr = NULL;
00358 
00359   extrema_ptr = initialize_extrema ();
00360 
00361   extrema_ptr->intensity = intensity;
00362   extrema_ptr->centroid  = centroid;
00363   extrema_ptr->count     = 1;
00364   extrema_ptr->sum       = intensity;
00365 
00366   head_extrema = add_extrema (extrema_ptr, head_extrema);
00367 
00368   return (head_extrema);
00369   
00370 }
00371 
00372 
00373 /*---------------------------------------------------------------------------*/
00374 /*
00375   Deallocate memory for this extrema.
00376 */
00377 
00378 void delete_extrema (extrema * extrema_ptr)
00379 {
00380   if (extrema_ptr != NULL)
00381     {
00382       if (extrema_ptr->centroid != NULL)
00383         {
00384           free (extrema_ptr->centroid);
00385           extrema_ptr->centroid = NULL;
00386         }
00387 
00388       free (extrema_ptr);
00389     }
00390 }
00391 
00392 
00393 /*---------------------------------------------------------------------------*/
00394 /*
00395   Remove one extrema from linked list of extrema.
00396   Reset extrema pointers, and recalculate nearest extrema distances 
00397   where needed.  Finally, delete the extrema from memory.
00398 */
00399 
00400 extrema * remove_extrema (extrema * aextrema, extrema * head_extrema)
00401 {
00402   extrema * extrema_ptr = NULL;
00403   extrema * next_extrema = NULL;
00404   
00405 
00406   if (head_extrema == NULL)  return;
00407 
00408 
00409   /*----- Remove aextrema from list; reset pointers -----*/
00410   if (head_extrema == aextrema)
00411     head_extrema = head_extrema->next_extrema;
00412   else
00413     {
00414       extrema_ptr = head_extrema;
00415       next_extrema = extrema_ptr->next_extrema;
00416       while (next_extrema != NULL)
00417         {
00418           if (next_extrema == aextrema)
00419             extrema_ptr->next_extrema = next_extrema->next_extrema;
00420           else
00421             extrema_ptr = next_extrema;
00422           
00423           next_extrema = extrema_ptr->next_extrema;
00424         }
00425     }
00426 
00427 
00428   /*----- Recalculate distances to nearest extrema -----*/
00429   if (head_extrema != NULL)
00430     {
00431       extrema_ptr = head_extrema;
00432       while (extrema_ptr != NULL)
00433         {
00434           if (extrema_ptr->nearest_extrema == aextrema) 
00435             {
00436               find_nearest_extrema (extrema_ptr, head_extrema);
00437             }
00438           extrema_ptr = extrema_ptr->next_extrema;
00439         }
00440     }
00441 
00442 
00443   /*----- Delete aextrema -----*/
00444   delete_extrema (aextrema);
00445 
00446   return (head_extrema);
00447 }
00448 
00449 
00450 /*---------------------------------------------------------------------------*/
00451 /*
00452   Average two extrema into one new extrema.
00453 */
00454 
00455 extrema * average_extrema (extrema * aextrema, extrema * bextrema)
00456 {
00457   extrema * abextrema = NULL;
00458   int na, nb;
00459   int i;
00460 
00461 
00462   /*----- Create new extrema -----*/
00463   abextrema = initialize_extrema ();
00464 
00465 
00466   /*----- Sum counts of extrema -----*/
00467   na = aextrema->count;
00468   nb = bextrema->count;
00469   abextrema->count = na + nb;
00470 
00471 
00472   /*----- Average of extrema intensities -----*/
00473   abextrema->intensity = 
00474     (na*aextrema->intensity + nb*bextrema->intensity) / (na+nb);
00475 
00476 
00477   abextrema->centroid = (float *) malloc (sizeof(float) * EX_DIMENSION);
00478   MTEST (abextrema->centroid);
00479 
00480 
00481   /*----- Average of extrema locations -----*/
00482   for (i = 0;  i < EX_DIMENSION;  i++)
00483     abextrema->centroid[i] = 
00484       (na*aextrema->centroid[i] + nb*bextrema->centroid[i]) / (na + nb);
00485   
00486 
00487   return (abextrema);
00488 }
00489 
00490 
00491 /*---------------------------------------------------------------------------*/
00492 /*
00493   Weighted average of two extrema into one new extrema.
00494 */
00495 
00496 extrema * weight_extrema (extrema * aextrema, extrema * bextrema)
00497 {
00498   extrema * abextrema = NULL;
00499   float suma, sumb;
00500   int i;
00501 
00502 
00503   /*----- Create new extrema -----*/
00504   abextrema = initialize_extrema ();
00505 
00506 
00507   /*----- Sum counts of extrema -----*/
00508   abextrema->count = aextrema->count + bextrema->count;
00509 
00510 
00511   /*----- Sum intensities of extrema -----*/
00512   suma = aextrema->sum;
00513   sumb = bextrema->sum;
00514   abextrema->sum = suma + sumb;
00515 
00516 
00517   /*----- Weighted average of extrema intensities -----*/
00518   abextrema->intensity = 
00519     (suma*aextrema->intensity + sumb*bextrema->intensity) / (suma + sumb);
00520 
00521 
00522   abextrema->centroid = (float *) malloc (sizeof(float) * EX_DIMENSION);
00523   MTEST (abextrema->centroid);
00524 
00525 
00526   /*----- Weighted average of extrema locations -----*/
00527   for (i = 0;  i < EX_DIMENSION;  i++)
00528     abextrema->centroid[i] = 
00529       (suma*aextrema->centroid[i] + sumb*bextrema->centroid[i]) / (suma+sumb);
00530   
00531 
00532   return (abextrema);
00533 }
00534 
00535 
00536 /*---------------------------------------------------------------------------*/
00537 /*
00538   Merge two extrema.
00539 */
00540   
00541 extrema * merge_extrema (extrema * aextrema, extrema * bextrema, 
00542                          extrema * head_extrema)
00543 {
00544   extrema * abextrema = NULL;
00545 
00546 
00547   switch (EX_merge)
00548     {
00549     case 1:   /*----- Merge extrema by removing the weakest -----*/
00550 
00551       if ((EX_maxima && ((aextrema->intensity) > (bextrema->intensity)))
00552           || (!EX_maxima && ((aextrema->intensity) < (bextrema->intensity))))
00553         {
00554           aextrema->count++;
00555           head_extrema = remove_extrema (bextrema, head_extrema);
00556         }
00557       else
00558         {
00559           bextrema->count++;      
00560           head_extrema = remove_extrema (aextrema, head_extrema);
00561         }
00562       break;
00563 
00564 
00565     case 2:   /*----- Merge extrema using average -----*/
00566 
00567       /*----- Merge two extrema into one new extrema -----*/
00568       abextrema = average_extrema (aextrema, bextrema);
00569 
00570       /*----- Remove the two original extrema from the linked list -----*/
00571       head_extrema = remove_extrema (aextrema, head_extrema);
00572       head_extrema = remove_extrema (bextrema, head_extrema);
00573 
00574 
00575       /*----- Add the merged extrema to the linked list -----*/
00576       head_extrema = add_extrema (abextrema, head_extrema);
00577   
00578       break;
00579 
00580 
00581     case 3:   /*----- Merge extrema using weighted average -----*/
00582  
00583       /*----- Merge two extrema into one new extrema -----*/
00584       abextrema = weight_extrema (aextrema, bextrema);
00585 
00586       /*----- Remove the two original extrema from the linked list -----*/
00587       head_extrema = remove_extrema (aextrema, head_extrema);
00588       head_extrema = remove_extrema (bextrema, head_extrema);
00589 
00590 
00591       /*----- Add the merged extrema to the linked list -----*/
00592       head_extrema = add_extrema (abextrema, head_extrema);
00593 
00594       break;
00595 
00596     }
00597 
00598   return (head_extrema);
00599 
00600 }
00601 
00602 
00603 /*---------------------------------------------------------------------------*/
00604 /*
00605   Sort extrema in order of intensity.
00606 */
00607 
00608 extrema * sort_extrema (extrema * head_extrema)
00609 {
00610   extrema * i  = NULL; 
00611   extrema * ip = NULL; 
00612   extrema * m  = NULL;
00613   extrema * mp = NULL;
00614   extrema * j  = NULL;
00615   extrema * jp = NULL;
00616   extrema * guard = NULL;
00617 
00618 
00619   /*----- Create guard extrema in case head extrema must be replaced -----*/
00620   guard = initialize_extrema();
00621   guard->next_extrema = head_extrema;
00622   ip = guard;
00623 
00624   while (ip->next_extrema != NULL)
00625     {
00626       /*----- Initialize search for next largest extrema -----*/
00627       i = ip->next_extrema;  /* current top of list */
00628       mp = ip;               /* extrema pointing to next largest extrema */
00629       m = i;                 /* next largest extrema */
00630       jp = i;
00631 
00632       /*----- Search through list for next largest extrema -----*/
00633       while (jp->next_extrema != NULL)
00634         {
00635           j = jp->next_extrema;
00636           if ( ((EX_maxima) && ((j->intensity) > (m->intensity)))
00637             || ((!EX_maxima) && ((j->intensity) < (m->intensity))) )
00638             {
00639               mp = jp;
00640               m = j;
00641             }
00642           jp = j;
00643         }
00644 
00645       /*----- Now move next largest extrema to top of list -----*/
00646       if (m != i)
00647         {
00648           ip->next_extrema = m;
00649           mp->next_extrema = m->next_extrema;
00650           m->next_extrema = i;
00651           i = m;
00652         }
00653 
00654       /*----- Move down the list -----*/
00655       ip = i;
00656         
00657     }
00658 
00659   
00660   /*----- Replace head extrema -----*/
00661   head_extrema = guard->next_extrema;
00662   delete_extrema (guard);
00663 
00664   return (head_extrema);
00665 }
00666 
00667 
00668 /*---------------------------------------------------------------------------*/
00669 /*
00670   Agglomerate extrema by merging the two extrema which are closest together.
00671 */
00672 
00673 extrema * agglomerate_extrema (extrema * head_extrema, float sep_dist)
00674 {
00675   const float MAX_DIST = 1.0e+30;
00676 
00677   extrema * extrema_ptr = NULL;
00678   extrema * aextrema    = NULL;
00679   extrema * bextrema    = NULL;
00680   float min_dist;
00681 
00682 
00683   min_dist = 0.0;
00684 
00685   while (min_dist < sep_dist)
00686     {
00687       /*----- Sort list of extrema -----*/
00688       if (EX_sort)  head_extrema = sort_extrema (head_extrema);
00689 
00690       /*----- Find the two extrema which are closest together -----*/
00691       min_dist = MAX_DIST;
00692       extrema_ptr = head_extrema;
00693       while (extrema_ptr != NULL)
00694         {
00695           if (extrema_ptr->nearest_dist < min_dist)
00696             {
00697               min_dist = extrema_ptr->nearest_dist;
00698               aextrema = extrema_ptr;
00699               bextrema = extrema_ptr->nearest_extrema;
00700             } 
00701           extrema_ptr = extrema_ptr->next_extrema;
00702         }
00703       
00704       /*
00705         printf ("min_dist = %f \n", min_dist);
00706       */
00707 
00708 
00709       /*----- Merge these two extrema -----*/
00710       if (min_dist < sep_dist)
00711         head_extrema = merge_extrema (aextrema, bextrema, head_extrema);
00712     }
00713 
00714   return (head_extrema);
00715 }
00716 
00717 
00718 /*---------------------------------------------------------------------------*/
00719 /*
00720   Convert from integer voxel coordinates to mm.
00721 */
00722 
00723 float * to_3dmm (int ixyz)
00724 {
00725   float * location = NULL;
00726   float x, y, z;
00727   int ix, jy, kz;
00728   THD_fvec3 fv;
00729   THD_ivec3 iv;
00730   
00731 
00732   location = (float *) malloc (sizeof(float) * EX_DIMENSION);
00733 
00734   IJK_TO_THREE (ixyz, ix, jy, kz, EX_nx, EX_nxy);
00735 
00736   iv.ijk[0] = ix;  iv.ijk[1] = jy;  iv.ijk[2] = kz;
00737 
00738   fv = THD_3dind_to_3dmm (EX_dset, iv);
00739   
00740   fv = THD_3dmm_to_dicomm (EX_dset, fv);
00741   
00742   x = fv.xyz[0];  y = fv.xyz[1];  z = fv.xyz[2];
00743   
00744   THD_dicom_to_coorder (&EX_cord, &x, &y, &z);
00745 
00746   location[0] = x;  location[1] = y;  location[2] = z;
00747 
00748   return (location);
00749 }
00750 
00751 
00752 /*---------------------------------------------------------------------------*/
00753 /*
00754   Find all extrema within one volume (or slice).
00755 */
00756 
00757 extrema * find_extrema 
00758 (
00759   float * fvol,                  /* volume or slice of floating point data */ 
00760   int num_nv,                    /* number of neighboring voxels */
00761   int nfirst,                    /* index of first voxel in volume or slice */
00762   int nlast                      /* index of last voxel in volume or slice */
00763 )
00764 
00765 {
00766   int ix, jy, kz, ixyz;          /* current voxel indices */
00767   int it, ijk;                   /* neighboring voxel indices   */
00768   int nx, ny, nz, nxy, nxyz;     /* numbers of voxels */
00769   float fval;                    /* intensity of current voxel */
00770   float * location = NULL;       /* current voxel coordinates */
00771   extrema * head_extrema = NULL; /* linked list of extrema */
00772   int passed_test;               /* flag for voxel is valid extema */
00773 
00774   
00775   /*----- Initialize local variables -----*/
00776   nx = EX_nx;      ny = EX_ny;      nz = EX_nz;  
00777   nxy = nx * ny;   nxyz = nx*ny*nz;
00778   
00779   
00780   /*----- Loop over voxels -----*/
00781   for (ixyz = nfirst;  ixyz < nlast;  ixyz++)
00782     {
00783 
00784       /*----- First, check if voxel is inside the mask -----*/
00785       if (!EX_mask[ixyz]) continue;
00786 
00787       /*----- Does voxel satisfy data threshold criterion? -----*/
00788       fval = fvol[ixyz];
00789       if (fabs(fval) < EX_data_thr)  continue; 
00790         
00791       /*----- Begin loop over neighboring voxels -----*/
00792       it = 0;
00793       passed_test = 1;
00794       while ((it < num_nv) && passed_test)
00795         {
00796           ijk = ixyz + EX_offset[it];
00797 
00798           /*----- Check for valid neighbor index -----*/
00799           if ((ijk < nfirst) || (ijk >= nlast))
00800             {  
00801               if (EX_interior)  passed_test = 0;
00802             }
00803   
00804           /*----- Early exit if neighbor is not inside the mask -----*/
00805           else if (!EX_mask[ijk])  
00806             {
00807               if (EX_interior)  passed_test = 0;
00808             }
00809             
00810           /*----- Test binary relation -----*/
00811           else
00812             switch (EX_relation)
00813               {
00814               case EX_GT: 
00815                 if (fval <= fvol[ijk])
00816                   passed_test = 0;  break;
00817               case EX_GE: 
00818                 if (fval <  fvol[ijk])
00819                   passed_test = 0;  break;
00820               case EX_LT: 
00821                 if (fval >= fvol[ijk])
00822                   passed_test = 0;  break;
00823               case EX_LE: 
00824                 if (fval >  fvol[ijk])
00825                   passed_test = 0;  break;
00826               }
00827 
00828           it++;
00829         }
00830 
00831 
00832       /*----- If extrema are found, save relevant information -----*/ 
00833       if (passed_test) 
00834         {
00835           location = to_3dmm (ixyz);
00836           head_extrema = new_extrema (location, fval, head_extrema);
00837         }
00838     }
00839   
00840 
00841   /*----- Agglomerate extrema -----*/
00842   if (EX_sep_dist > 0.0)  head_extrema = agglomerate_extrema (head_extrema,
00843                                                               EX_sep_dist);
00844 
00845 
00846   /*----- Sort list of extrema -----*/
00847   if (EX_sort)  head_extrema = sort_extrema (head_extrema);
00848 
00849 
00850   return (head_extrema);
00851 }
00852 
00853 
00854 /*---------------------------------------------------------------------------*/
00855 /*
00856   Convert from mm to integer voxel coordinates.
00857 */
00858 
00859 int from_3dmm (float * location)
00860 {
00861   float x, y, z;
00862   THD_fvec3 fv;
00863   int ix, jy, kz;
00864   THD_ivec3 iv;
00865   int ixyz;
00866 
00867   x = location[0];  y = location[1];  z = location[2];
00868 
00869   THD_coorder_to_dicom (&EX_cord, &x, &y, &z);
00870 
00871   fv.xyz[0] = x;  fv.xyz[1] = y;  fv.xyz[2] = z;
00872 
00873   fv = THD_dicomm_to_3dmm (EX_dset, fv);
00874  
00875   iv = THD_3dmm_to_3dind (EX_dset, fv);
00876     
00877   ix = iv.ijk[0];  jy = iv.ijk[1];  kz = iv.ijk[2];
00878  
00879   if (ix < 0) ix = 0;   if (ix > EX_nx-1) ix = EX_nx-1;
00880   if (jy < 0) jy = 0;   if (jy > EX_ny-1) jy = EX_ny-1;
00881   if (kz < 0) kz = 0;   if (kz > EX_nz-1) kz = EX_nz-1;
00882 
00883   ixyz = THREE_TO_IJK (ix, jy, kz, EX_nx, EX_nxy);
00884 
00885   return (ixyz);
00886 
00887 }
00888 
00889 
00890 /*---------------------------------------------------------------------------*/
00891 /*
00892   Save voxel contents of this extrema into byte array (sub-brick).
00893 */
00894 
00895 void save_extrema (extrema * extrema_ptr, byte iextrema, byte * bar)
00896 {
00897   int ixyz;
00898 
00899   ixyz = from_3dmm (extrema_ptr->centroid);
00900 
00901   bar[ixyz] = 1;
00902 
00903 }
00904 
00905 
00906 /*---------------------------------------------------------------------------*/
00907 /*
00908   Save voxel contents of all extrema into byte array (sub-brick).
00909 */
00910 
00911 void save_all_extrema (extrema * extrema_ptr, byte * bar)
00912 {
00913   byte iextrema = 0;
00914 
00915   while (extrema_ptr != NULL)
00916     {
00917       iextrema++;
00918       save_extrema (extrema_ptr, iextrema, bar);
00919       extrema_ptr = extrema_ptr->next_extrema;
00920     }
00921 
00922 }
00923 
00924 
00925 /*---------------------------------------------------------------------------*/
00926 /*
00927   Display help file.
00928 */
00929 
00930 void EX_Syntax(void)
00931 {
00932    printf(
00933     "This program finds local extrema (minima or maxima) of the input       \n"
00934     "dataset values for each sub-brick of the input dataset.  The extrema   \n"
00935     "may be determined either for each volume, or for each individual slice.\n"
00936     "Only those voxels whose corresponding intensity value is greater than  \n"
00937     "the user specified data threshold will be considered.                  \n"
00938     "\nUsage: 3dExtrema  options  datasets                                  \n"
00939     "where the options are:                                                 \n"
00940    ) ;
00941 
00942    printf(
00943     "-prefix pname    = Use 'pname' for the output dataset prefix name.     \n"
00944     "  OR                 [default = NONE; only screen output]              \n"
00945     "-output pname                                                          \n"
00946     "                                                                       \n"
00947     "-session dir     = Use 'dir' for the output dataset session directory. \n"
00948     "                     [default='./'=current working directory]          \n"
00949     "                                                                       \n"
00950     "-quiet           = Flag to suppress screen output                      \n"
00951     "                                                                       \n"
00952     "-mask_file mname = Use mask statistic from file mname.                 \n"
00953     "                   Note: If file mname contains more than 1 sub-brick, \n"
00954     "                   the mask sub-brick must be specified!               \n"
00955     "-mask_thr m        Only voxels whose mask statistic is greater         \n"
00956     "                   than m in abolute value will be considered.         \n"
00957     "                                                                       \n"
00958     "-data_thr d        Only voxels whose value (intensity) is greater      \n"
00959     "                   than d in abolute value will be considered.         \n"
00960     "                                                                       \n"
00961     "-sep_dist d        Min. separation distance [mm] for distinct extrema  \n"
00962     "                                                                       \n"
00963     "Choose type of extrema (one and only one choice):                      \n"
00964     "-minima            Find local minima.                                  \n"
00965     "-maxima            Find local maxima.                                  \n"
00966     "                                                                       \n"
00967     "Choose form of binary relation (one and only one choice):              \n"
00968     "-strict            >  for maxima,  <  for minima                       \n"
00969     "-partial           >= for maxima,  <= for minima                       \n"
00970     "                                                                       \n"
00971     "Choose boundary criteria (one and only one choice):                    \n"
00972     "-interior          Extrema must be interior points (not on boundary)   \n"
00973     "-closure           Extrema may be boudary points                       \n"
00974     "                                                                       \n"
00975     "Choose domain for finding extrema (one and only one choice):           \n"
00976     "-slice             Each slice is considered separately                 \n"
00977     "-volume            The volume is considered as a whole                 \n"
00978     "                                                                       \n"
00979     "Choose option for merging of extrema (one and only one choice):        \n"
00980     "-remove            Remove all but strongest of neighboring extrema     \n"
00981     "-average           Replace neighboring extrema by average              \n"
00982     "-weight            Replace neighboring extrema by weighted average     \n"
00983     "                                                                       \n"
00984     "Command line arguments after the above are taken to be input datasets. \n"
00985     "\n") ;
00986 
00987    printf("\n" MASTER_SHORTHELP_STRING ) ;
00988 
00989    exit(0) ;
00990 
00991 }
00992 
00993 
00994 /*---------------------------------------------------------------------------*/
00995 /*
00996    Read the arguments, load the global variables
00997 
00998 */
00999 
01000 int EX_read_opts( int argc , char * argv[] )
01001 {
01002    int nopt = 1 ;
01003    char message[80];        /* error message */
01004 
01005 
01006    while( nopt < argc ){
01007 
01008       /**** -prefix prefix ****/
01009       if( strcmp(argv[nopt],"-prefix") == 0 ||
01010           strcmp(argv[nopt],"-output") == 0   ){
01011          nopt++ ;
01012          if( nopt >= argc ){
01013             EX_error (" need argument after -prefix!");
01014          }
01015          EX_output_prefix = (char *) malloc (sizeof(char) * THD_MAX_PREFIX); 
01016          MCW_strncpy( EX_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
01017          continue ;
01018       }
01019 
01020 
01021       /**** -session directory ****/
01022       if( strcmp(argv[nopt],"-session") == 0 ){
01023          nopt++ ;
01024          if( nopt >= argc ){
01025             EX_error (" need argument after -session!"); 
01026          }
01027          EX_session = (char *) malloc (sizeof(char) * THD_MAX_NAME); 
01028          MCW_strncpy( EX_session , argv[nopt++] , THD_MAX_NAME ) ;
01029          continue ;
01030       }
01031 
01032 
01033       /**** -quiet ****/
01034       if( strcmp(argv[nopt],"-quiet") == 0 ){
01035          EX_quiet = 1;
01036          nopt++ ; continue ;
01037       }
01038 
01039 
01040       /**** -minima ****/
01041       if( strcmp(argv[nopt],"-minima") == 0 ){
01042          EX_maxima = 0;
01043          nopt++;  continue;
01044       }
01045 
01046 
01047       /**** -maxima ****/
01048       if( strcmp(argv[nopt],"-maxima") == 0 ){
01049          EX_maxima = 1;
01050          nopt++;  continue;
01051       }
01052 
01053 
01054       /**** -strict ****/
01055       if( strcmp(argv[nopt],"-strict") == 0 ){
01056          EX_strict = 1;
01057          nopt++;  continue;
01058       }
01059 
01060 
01061       /**** -partial ****/
01062       if( strcmp(argv[nopt],"-partial") == 0 ){
01063          EX_strict = 0;
01064          nopt++;  continue;
01065       }
01066 
01067 
01068       /**** -interior ****/
01069       if( strcmp(argv[nopt],"-interior") == 0 ){
01070          EX_interior = 1;
01071          nopt++;  continue;
01072       }
01073 
01074 
01075       /**** -closure ****/
01076       if( strcmp(argv[nopt],"-closure") == 0 ){
01077          EX_interior = 0;
01078          nopt++;  continue;
01079       }
01080 
01081 
01082       /**** -slice ****/
01083       if( strcmp(argv[nopt],"-slice") == 0 ){
01084          EX_slice = 1;
01085          nopt++;  continue;
01086       }
01087 
01088 
01089       /**** -volume ****/
01090       if( strcmp(argv[nopt],"-volume") == 0 ){
01091          EX_slice = 0;
01092          nopt++;  continue;
01093       }
01094 
01095 
01096       /**** -sort ****/
01097       if( strcmp(argv[nopt],"-sort") == 0 ){
01098          EX_sort = 1;
01099          nopt++;  continue;
01100       }
01101 
01102 
01103       /**** -mask_file mname ****/
01104       if( strcmp(argv[nopt],"-mask_file") == 0 ){
01105          nopt++ ;
01106          if( nopt >= argc ){
01107             EX_error (" need 1 argument after -mask_file"); 
01108          }
01109 
01110          EX_mask_filename = (char *) malloc (sizeof(char) * THD_MAX_NAME); 
01111          MCW_strncpy( EX_mask_filename , argv[nopt++] , THD_MAX_NAME ) ;
01112          continue;
01113       }
01114 
01115 
01116       /**** -mask_thr m ****/
01117       if( strcmp(argv[nopt],"-mask_thr") == 0 ){
01118          float fval;
01119          nopt++ ;
01120          if( nopt >= argc ){
01121             EX_error (" need 1 argument after -mask_thr"); 
01122          }
01123          sscanf (argv[nopt], "%f", &fval); 
01124          if (fval < 0.0){
01125             EX_error (" Require mask_thr >= 0.0 ");
01126          }
01127          EX_mask_thr = fval;
01128          nopt++;  continue;
01129 
01130       }
01131 
01132 
01133       /**** -data_thr d ****/
01134       if( strcmp(argv[nopt],"-data_thr") == 0 ){
01135          float fval;
01136          nopt++ ;
01137          if( nopt >= argc ){
01138             EX_error (" need 1 argument after -data_thr"); 
01139          }
01140          sscanf (argv[nopt], "%f", &fval); 
01141          if (fval < 0.0){
01142             EX_error (" Require data_thr >= 0.0 ");
01143          }
01144          EX_data_thr = fval;
01145          nopt++;  continue;
01146 
01147       }
01148       
01149 
01150       /**** -remove ****/
01151       if( strcmp(argv[nopt],"-remove") == 0 ){
01152          EX_merge = 1;
01153          nopt++;  continue;
01154       }
01155 
01156 
01157       /**** -average ****/
01158       if( strcmp(argv[nopt],"-average") == 0 ){
01159          EX_merge = 2;
01160          nopt++;  continue;
01161       }
01162 
01163 
01164       /**** -weight ****/
01165       if( strcmp(argv[nopt],"-weight") == 0 ){
01166          EX_merge = 3;
01167          nopt++;  continue;
01168       }
01169 
01170 
01171       /**** -sep_dist d ****/
01172       if( strcmp(argv[nopt],"-sep_dist") == 0 ){
01173          float fval;
01174          nopt++ ;
01175          if( nopt >= argc ){
01176             EX_error (" need 1 argument after -sep_dist"); 
01177          }
01178          sscanf (argv[nopt], "%f", &fval); 
01179          if (fval < 0.0){
01180             EX_error (" Require data_thr >= 0.0 ");
01181          }
01182          EX_sep_dist = fval;
01183          nopt++;  continue;
01184 
01185       }
01186 
01187 
01188      /*----- Invalid option -----*/
01189       if( argv[nopt][0] == '-' ){
01190         sprintf (message, " Unknown option: %s ", argv[nopt]);
01191         EX_error (message);
01192       }
01193 
01194 
01195       /*----- Remaining inputs should be parameter datasets -----*/
01196       break;
01197 
01198 
01199    }  /* end of loop over command line arguments */
01200 
01201    
01202    /*----- Check for valid input -----*/
01203    if ((EX_merge == 3) && (EX_maxima == 0))
01204      EX_error ("-weight option is not compatible with -minima option");
01205 
01206 
01207    /*----- Set flag for binary relation -----*/
01208    if ((EX_maxima == 1) && (EX_strict == 1))  EX_relation = EX_GT;
01209    if ((EX_maxima == 1) && (EX_strict == 0))  EX_relation = EX_GE;
01210    if ((EX_maxima == 0) && (EX_strict == 1))  EX_relation = EX_LT;
01211    if ((EX_maxima == 0) && (EX_strict == 0))  EX_relation = EX_LE;
01212 
01213 
01214    return (nopt) ;
01215 }
01216 
01217 
01218 /*---------------------------------------------------------------------------*/
01219 /*
01220   Routine to initialize the program: get all operator inputs; create mask
01221   for voxels above mask.
01222 */
01223 
01224 void * initialize_program (int argc, char * argv[], int * nopt)
01225 {
01226   const int MIN_NTHR = 10;    /* minimum number of voxels above threshold */
01227 
01228   int iv;                  /* index number of sub-brick */
01229   void * vfim = NULL;      /* sub-brick data pointer */
01230   float * ffim = NULL;     /* sub-brick data in floating point format */
01231   int ixyz, nthr;          /* voxel indices */
01232   int nx, ny, nz, nxy, nxyz = 0;     /* numbers of voxels */
01233   char message[80];        /* error message */
01234 
01235 
01236   /*----- Save command line for history notes -----*/
01237   commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
01238 
01239 
01240   /*----- Does user request help menu? -----*/
01241   if( argc < 2 || strcmp(argv[1],"-help") == 0 ) EX_Syntax() ;
01242 
01243   
01244   /*----- Add to program log -----*/
01245   AFNI_logger (PROGRAM_NAME,argc,argv); 
01246 
01247 
01248   /*----- Read input options -----*/
01249   *nopt = EX_read_opts( argc , argv ) ;
01250 
01251 
01252   /*----- Get coordinate directions -----*/
01253   THD_coorder_fill (my_getenv("AFNI_ORIENT"), &EX_cord);
01254 
01255 
01256   /*----- Open the mask dataset -----*/
01257   if (EX_mask_filename != NULL)
01258     {
01259       if (!EX_quiet)  printf ("Reading mask dataset: %s \n", EX_mask_filename);
01260       DOPEN (EX_dset, EX_mask_filename);
01261 
01262       if (EX_dset == NULL)
01263         {
01264           sprintf (message, "Cannot open mask dataset %s", EX_mask_filename); 
01265           EX_error (message);
01266         }
01267 
01268       if (DSET_NVALS(EX_dset) != 1)
01269         EX_error ("Must specify single sub-brick for mask data");
01270 
01271 
01272       /*----- Get dimensions of mask dataset -----*/
01273       nx   = DSET_NX(EX_dset);   
01274       ny   = DSET_NY(EX_dset);   
01275       nz   = DSET_NZ(EX_dset);
01276       nxy  = nx * ny;  nxyz = nx*ny*nz;
01277 
01278 
01279       /*----- Allocate memory for float data -----*/
01280       ffim = (float *) malloc (sizeof(float) * nxyz);   MTEST (ffim);
01281 
01282 
01283       /*----- Convert mask dataset sub-brick to floats (in ffim) -----*/
01284       iv = DSET_PRINCIPAL_VALUE (EX_dset);
01285       SUB_POINTER (EX_dset, iv, 0, vfim);
01286       EDIT_coerce_scale_type (nxyz, DSET_BRICK_FACTOR(EX_dset,iv),
01287                               DSET_BRICK_TYPE(EX_dset,iv), vfim,   /* input  */
01288                               MRI_float                   , ffim); /* output */
01289   
01290       
01291       /*----- Allocate memory for mask volume -----*/
01292       EX_mask = (byte *) malloc (sizeof(byte) * nxyz);
01293       MTEST (EX_mask);
01294       
01295       
01296       /*----- Create mask of voxels above mask threshold -----*/
01297       nthr = 0;
01298       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01299         if (fabs(ffim[ixyz]) >= EX_mask_thr)  
01300           { 
01301             EX_mask[ixyz] = 1;
01302             nthr++;
01303           }
01304         else
01305           EX_mask[ixyz] = 0;
01306 
01307       if (!EX_quiet)  
01308         printf ("Number of voxels above mask threshold = %d \n", nthr);
01309       if (nthr < MIN_NTHR)  
01310         {
01311           sprintf (message, 
01312                    "Only %d voxels above mask threshold.  Cannot continue.", 
01313                    nthr);
01314           EX_error (message);
01315         }
01316 
01317 
01318       /*----- Delete floating point sub-brick -----*/
01319       if (ffim != NULL) { free (ffim);   ffim = NULL; }
01320 
01321       /*----- Delete mask dataset -----*/
01322       THD_delete_3dim_dataset (EX_dset, False);  EX_dset = NULL ;
01323 
01324     }
01325 
01326 
01327   /*----- Open first input dataset -----*/
01328   {
01329     /*----- Check if this is an input option -----*/
01330     if (argv[*nopt][0] == '-')
01331       EX_error ("ALL input options must precede ALL input datasets");
01332 
01333     /*----- Open the input dataset -----*/
01334     if (!EX_quiet)  printf ("Reading input dataset: %s \n", argv[*nopt]);
01335     EX_dset = THD_open_dataset(argv[*nopt]);
01336 
01337     if (EX_dset == NULL)
01338       {
01339         sprintf (message, "Cannot open input dataset %s", argv[*nopt]); 
01340         EX_error (message);
01341       }
01342 
01343     /*----- Get dimensions of input dataset -----*/
01344     if (nxyz == 0)
01345       {
01346         nx   = DSET_NX(EX_dset);   
01347         ny   = DSET_NY(EX_dset);   
01348         nz   = DSET_NZ(EX_dset);
01349         nxy  = nx * ny;  nxyz = nx*ny*nz;
01350       }
01351     
01352     /*----- Create mask, if not already done -----*/
01353     if (EX_mask == NULL)
01354       {
01355         /*----- Allocate memory for mask volume -----*/
01356         EX_mask = (byte *) malloc (sizeof(byte) * nxyz);
01357         MTEST (EX_mask);
01358         for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01359           EX_mask[ixyz] = 1;
01360       }
01361 
01362   }
01363 
01364   
01365   /*----- Save dimensions of dataset for compatibility test -----*/
01366   EX_nx   = nx;     EX_ny   = ny;     EX_nz   = nz;
01367   EX_nxy  = nxy;    EX_nxyz = nxyz; 
01368 
01369 
01370   /*----- The offset array definitions were borrowed from plug_maxima -----*/
01371   EX_offset[ 0] = +1;     EX_offset[ 8] = nxy;       EX_offset[17] = -nxy;
01372   EX_offset[ 1] = -1;     EX_offset[ 9] = nxy+1;     EX_offset[18] = -nxy+1;
01373   EX_offset[ 2] =  nx;    EX_offset[10] = nxy-1;     EX_offset[19] = -nxy-1; 
01374   EX_offset[ 3] =  nx+1;  EX_offset[11] = nxy+nx;    EX_offset[20] = -nxy+nx;
01375   EX_offset[ 4] =  nx-1;  EX_offset[12] = nxy+nx+1;  EX_offset[21] = -nxy+nx+1;
01376   EX_offset[ 5] = -nx;    EX_offset[13] = nxy+nx-1;  EX_offset[22] = -nxy+nx-1;
01377   EX_offset[ 6] = -nx+1;  EX_offset[14] = nxy-nx;    EX_offset[23] = -nxy-nx;
01378   EX_offset[ 7] = -nx-1;  EX_offset[15] = nxy-nx+1;  EX_offset[24] = -nxy-nx+1;
01379                           EX_offset[16] = nxy-nx-1;  EX_offset[25] = -nxy-nx-1;
01380   
01381 }
01382 
01383 
01384 /*---------------------------------------------------------------------------*/
01385 /*
01386   Process all input datasets.
01387 */
01388 
01389 void process_all_datasets (int argc, char * argv[], int nopt)
01390 {
01391   THD_3dim_dataset * input_dset=NULL;   /* input dataset(s) */
01392 
01393   int iv;                        /* index number of sub-brick */
01394   void * vfim = NULL;            /* sub-brick data pointer */
01395   float * ffim = NULL;           /* sub-brick data in floating point format */
01396   int ibrick, nbricks;           /* sub-brick indices */
01397   char message[80];              /* error message */
01398   int nx, ny, nz, nxy, nxyz;     /* numbers of voxels */
01399   int kz, nfirst, nlast;         /* range of voxel indices */
01400 
01401 
01402   /*----- Initialize local variables -----*/
01403   nx = EX_nx;  ny = EX_ny;  nz = EX_nz;  
01404   nxy = nx * ny;  nxyz = nx*ny*nz;
01405   
01406 
01407   /*----- Allocate memory for float data -----*/
01408   ffim = (float *) malloc (sizeof(float) * EX_nxyz);   MTEST (ffim);
01409 
01410 
01411   /*----- Begin loop over input datasets -----*/
01412   EX_nbricks = 0;
01413   while (nopt < argc)
01414     {
01415       /*----- Check if this is an input option -----*/
01416       if (argv[nopt][0] == '-')
01417         EX_error ("ALL input options must precede ALL input datasets");
01418 
01419       /*----- Open the input dataset -----*/
01420       if (!EX_quiet)  printf ("Reading input dataset: %s \n", argv[nopt]);
01421       DOPEN (input_dset, argv[nopt]);
01422 
01423       if (input_dset == NULL)
01424         {
01425           sprintf (message, "Cannot open input dataset %s", argv[nopt]); 
01426           EX_error (message);
01427         }
01428 
01429       /*----- Test for dataset compatibility -----*/
01430       if ((EX_nx != DSET_NX(input_dset)) || (EX_ny != DSET_NY(input_dset))
01431           || (EX_nz != DSET_NZ(input_dset)))
01432         {
01433           sprintf (message, "Input dataset %s has incompatible dimensions",
01434                    argv[nopt]); 
01435           EX_error (message);
01436         }
01437         
01438      
01439       /*----- Get number of volumes specified for this dataset -----*/
01440       nbricks = DSET_NVALS(input_dset);
01441 
01442 
01443       /*----- Loop over sub-bricks selected from input dataset -----*/
01444       for (ibrick = 0;  ibrick < nbricks;  ibrick++)
01445         {
01446           if (!EX_quiet)  printf ("Reading volume #%d \n", EX_nbricks);
01447 
01448           SUB_POINTER (input_dset, ibrick, 0, vfim);
01449           EDIT_coerce_scale_type 
01450             (EX_nxyz, DSET_BRICK_FACTOR(input_dset,ibrick),
01451              DSET_BRICK_TYPE(input_dset,ibrick), vfim,   /* input  */
01452              MRI_float                         , ffim);  /* output */
01453           
01454           if (EX_slice)  /*----- Find extrema slice-by-slice -----*/
01455             {
01456               for (kz = 0;  kz < nz;  kz++)
01457                 {
01458                   nfirst = kz*nxy;  nlast = nfirst + nxy;
01459                   EX_head_extrema[EX_num_ll] 
01460                     = find_extrema (ffim, 8, nfirst, nlast);
01461                   EX_num_ll++;
01462                   if (EX_num_ll >= EX_MAX_LL)  
01463                     EX_error ("Exceeded Max. Number of Linked Lists");
01464                 }
01465             }
01466           else           /*----- Find extrema for entire volume -----*/
01467             {
01468               EX_head_extrema[EX_num_ll] = find_extrema (ffim, 26, 0, nxyz);
01469               EX_num_ll++;
01470               if (EX_num_ll >= EX_MAX_LL)  
01471                 EX_error ("Exceeded Max. Number of Linked Lists");
01472             }
01473 
01474           /*----- Increment count of sub-bricks -----*/
01475           EX_nbricks++;  
01476         }
01477 
01478 
01479       /*----- Delete input dataset -----*/
01480       THD_delete_3dim_dataset (input_dset, False);  input_dset = NULL ;
01481      
01482       nopt++;
01483     }
01484 
01485 
01486   /*----- Delete floating point sub-brick -----*/
01487   if (ffim != NULL) { free (ffim);   ffim = NULL; }
01488 
01489 
01490   if (!EX_quiet)  printf ("Number of volumes = %d \n", EX_nbricks);
01491   if (EX_nbricks < 1)  EX_error ("No input data?");
01492 
01493 
01494 }
01495 
01496 
01497 /*---------------------------------------------------------------------------*/
01498 /*
01499   Write extrema to output bucket dataset.
01500 */
01501 
01502 void write_bucket ()
01503 
01504 {
01505   THD_3dim_dataset * new_dset = NULL;      /* output bucket dataset */
01506   int ixyz, nxyz;                  /* voxel indices */
01507   int kz, nz;                      /* slice indices */
01508   extrema * head_extrema = NULL;   /* pointer to linked list of extrema */
01509   byte ** bar = NULL;              /* array of extrema sub-bricks */
01510   int nbricks;                     /* number of extrema sub-bricks */
01511   int ibrick;                      /* extrema sub-brick index */
01512   int ierror;                      /* number of errors in editing data */
01513 
01514   char message[80];           /* error message */
01515 
01516 
01517   /*----- Initialize local variables -----*/
01518   nz = EX_nz;  
01519   nxyz = EX_nxyz;
01520   nbricks = EX_nbricks;
01521   if (!EX_quiet) 
01522     printf ("\nOutput dataset will have %d sub-bricks\n", nbricks);
01523 
01524 
01525   /*-- Make an empty copy of input dataset, for eventual output --*/
01526   new_dset = EDIT_empty_copy (EX_dset);
01527 
01528 
01529   /*----- Record history of dataset -----*/
01530   tross_Copy_History (EX_dset, new_dset);
01531   if( commandline != NULL ) tross_Append_History( new_dset , commandline ) ;
01532 
01533 
01534   /*----- Modify some structural properties.  Note that the nbricks
01535           just make empty sub-bricks, without any data attached. -----*/
01536   ierror = EDIT_dset_items (new_dset,
01537                             ADN_prefix,          EX_output_prefix,
01538                             ADN_directory_name,  EX_session,
01539                             ADN_type,            HEAD_FUNC_TYPE,
01540                             ADN_func_type,       FUNC_BUCK_TYPE,
01541                             ADN_ntt,             0,               /* no time */
01542                             ADN_nvals,           nbricks,
01543                             ADN_malloc_type,     DATABLOCK_MEM_MALLOC ,  
01544                             ADN_none ) ;
01545   
01546   if( ierror > 0 )
01547     {
01548       sprintf (message, 
01549               " %d errors in attempting to create bucket dataset! ", 
01550               ierror);
01551       EX_error (message);
01552     }
01553   
01554   if (THD_is_file(DSET_HEADNAME(new_dset))) 
01555     {
01556       sprintf (message,
01557               " Output dataset file %s already exists--cannot continue! ",
01558               DSET_HEADNAME(new_dset));
01559       EX_error (message);
01560     }
01561 
01562 
01563   /*----- Allocate memory -----*/
01564   bar  = (byte **) malloc (sizeof(byte *) * nbricks);
01565   MTEST (bar);
01566   
01567 
01568   /*----- Save extrema into sub-bricks -----*/
01569   for (ibrick = 0;  ibrick < nbricks;  ibrick++)
01570     {
01571       /*----- allocate memory for output sub-brick -----*/
01572       bar[ibrick]  = (byte *) malloc (sizeof(byte) * nxyz);
01573       MTEST (bar[ibrick]);
01574 
01575       /*----- Save extrema into output sub-brick -----*/
01576       for (ixyz = 0;  ixyz < nxyz;  ixyz++)        
01577         bar[ibrick][ixyz] = 0;
01578       if (EX_slice)
01579         for (kz = 0;  kz < nz;  kz++)
01580           {
01581             head_extrema = EX_head_extrema[kz+ibrick*nz];
01582             save_all_extrema (head_extrema, bar[ibrick]);           
01583           }
01584       else
01585         {
01586           head_extrema = EX_head_extrema[ibrick];
01587           save_all_extrema (head_extrema, bar[ibrick]); 
01588         }
01589 
01590       /*----- attach bar[ib] to be sub-brick #ibrick -----*/
01591       EDIT_substitute_brick (new_dset, ibrick, MRI_byte, bar[ibrick]);
01592 
01593     }
01594 
01595 
01596   /*----- Output the extrema dataset -----*/
01597   if( !EX_quiet ) fprintf(stderr,"Computing sub-brick statistics\n") ;
01598   THD_load_statistics( new_dset ) ;
01599 
01600   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01601   if( !EX_quiet ) fprintf(stderr,"Wrote output dataset: %s\n", DSET_BRIKNAME(new_dset) );
01602   
01603 
01604   /*----- Deallocate memory for extrema dataset -----*/   
01605   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01606   
01607   
01608 }
01609 
01610 
01611 /*---------------------------------------------------------------------------*/
01612 /*
01613   Output the results.
01614 */
01615 
01616 void output_results ()
01617 {
01618   int islice,  ibrick;
01619   int nslices, nbricks;
01620 
01621     
01622   /*----- Initialize local variables -----*/
01623   nbricks = EX_nbricks;
01624   nslices = EX_nz;
01625 
01626 
01627   /*----- Write extrema to screen -----*/
01628   if (!EX_quiet)
01629     if (EX_slice)
01630       for (ibrick = 0;  ibrick < nbricks;  ibrick++)
01631         {
01632           for (islice = 0;  islice < nslices;  islice++)
01633             {
01634               print_all_extrema (ibrick, islice, 
01635                                  EX_head_extrema[islice + ibrick*nslices]);
01636             }
01637         }
01638     else
01639       for (ibrick = 0;  ibrick < nbricks;  ibrick++)
01640         {
01641           print_all_extrema (ibrick, -1, EX_head_extrema[ibrick]);
01642         }
01643 
01644 
01645   /*----- Write extrema to bucket dataset -----*/
01646   if (EX_output_prefix != NULL)
01647     write_bucket ();
01648 
01649 }
01650 
01651 
01652 /*---------------------------------------------------------------------------*/
01653 
01654 int main( int argc , char * argv[] )
01655 
01656 {
01657   int nopt;
01658   
01659 
01660   /*----- Identify software -----*/
01661 #if 0
01662   printf ("\n\n");
01663   printf ("Program:          %s \n", PROGRAM_NAME);
01664   printf ("Author:           %s \n", PROGRAM_AUTHOR); 
01665   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
01666   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
01667   printf ("\n");
01668 #else
01669   PRINT_VERSION("3dExtrema") ; mainENTRY("3dExtrema main") ;
01670 #endif
01671 
01672 
01673   /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
01674   machdep() ; 
01675   { int new_argc ; char ** new_argv ;
01676   addto_args( argc , argv , &new_argc , &new_argv ) ;
01677   if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
01678   }
01679   
01680 
01681   /*----- Initialize program:  get all operator inputs; 
01682     create mask for voxels above mask threshold -----*/
01683   initialize_program (argc, argv, &nopt);
01684 
01685 
01686   /*----- Process all datasets -----*/
01687   process_all_datasets (argc, argv, nopt);
01688 
01689 
01690   /*----- Output the results -----*/
01691   output_results ();
01692 
01693 
01694   exit(0) ;
01695 }
01696 
01697 
01698 /*---------------------------------------------------------------------------*/
01699 
01700