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 #define PROGRAM_NAME "3dExtrema"
00026 #define PROGRAM_AUTHOR "B. Douglas Ward"
00027 #define PROGRAM_INITIAL "12 April 2001"
00028 #define PROGRAM_LATEST "15 August 2001"
00029
00030
00031
00032
00033
00034
00035
00036 #include "mrilib.h"
00037
00038
00039
00040
00041
00042
00043
00044 struct extrema;
00045
00046 typedef struct extrema
00047 {
00048 float intensity;
00049 float * centroid;
00050 int count;
00051 float sum;
00052 float nearest_dist;
00053 struct extrema * nearest_extrema;
00054 struct extrema * next_extrema;
00055 } extrema;
00056
00057
00058
00059
00060 #define EX_MAX_LL 5000
00061 #define EX_MAX_NV 26
00062 #define EX_DIMENSION 3
00063
00064 #define EX_GT 1
00065 #define EX_GE 2
00066 #define EX_LT 3
00067 #define EX_LE 4
00068
00069 static THD_coorder EX_cord;
00070
00071 static int EX_nx, EX_ny, EX_nz,
00072 EX_nxy, EX_nxyz;
00073
00074 static int EX_quiet = 0;
00075 static int EX_relation = 1;
00076 static int EX_maxima = 1;
00077 static int EX_strict = 1;
00078 static int EX_interior = 1;
00079 static int EX_slice = 1;
00080 static int EX_sort = 1;
00081 static int EX_merge = 1;
00082 static float EX_sep_dist = 0.0;
00083
00084 static float EX_mask_thr = 1.0;
00085 static float EX_data_thr = 0.0;
00086 static int EX_nbricks = 0;
00087
00088 static extrema * EX_head_extrema[EX_MAX_LL];
00089 static int EX_num_ll = 0;
00090
00091 static int EX_offset[EX_MAX_NV];
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;
00098
00099 static char * commandline = NULL;
00100
00101 static THD_3dim_dataset * EX_dset = NULL;
00102
00103
00104
00105
00106
00107 #define MTEST(ptr) \
00108 if((ptr)==NULL) \
00109 ( printf ("Cannot allocate memory \n"), exit(1) )
00110
00111
00112
00113
00114
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
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
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
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
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
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
00267
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
00291
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
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
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
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
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
00396
00397
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
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
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
00444 delete_extrema (aextrema);
00445
00446 return (head_extrema);
00447 }
00448
00449
00450
00451
00452
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
00463 abextrema = initialize_extrema ();
00464
00465
00466
00467 na = aextrema->count;
00468 nb = bextrema->count;
00469 abextrema->count = na + nb;
00470
00471
00472
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
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
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
00504 abextrema = initialize_extrema ();
00505
00506
00507
00508 abextrema->count = aextrema->count + bextrema->count;
00509
00510
00511
00512 suma = aextrema->sum;
00513 sumb = bextrema->sum;
00514 abextrema->sum = suma + sumb;
00515
00516
00517
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
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
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:
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:
00566
00567
00568 abextrema = average_extrema (aextrema, bextrema);
00569
00570
00571 head_extrema = remove_extrema (aextrema, head_extrema);
00572 head_extrema = remove_extrema (bextrema, head_extrema);
00573
00574
00575
00576 head_extrema = add_extrema (abextrema, head_extrema);
00577
00578 break;
00579
00580
00581 case 3:
00582
00583
00584 abextrema = weight_extrema (aextrema, bextrema);
00585
00586
00587 head_extrema = remove_extrema (aextrema, head_extrema);
00588 head_extrema = remove_extrema (bextrema, head_extrema);
00589
00590
00591
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
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
00620 guard = initialize_extrema();
00621 guard->next_extrema = head_extrema;
00622 ip = guard;
00623
00624 while (ip->next_extrema != NULL)
00625 {
00626
00627 i = ip->next_extrema;
00628 mp = ip;
00629 m = i;
00630 jp = i;
00631
00632
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
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
00655 ip = i;
00656
00657 }
00658
00659
00660
00661 head_extrema = guard->next_extrema;
00662 delete_extrema (guard);
00663
00664 return (head_extrema);
00665 }
00666
00667
00668
00669
00670
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
00688 if (EX_sort) head_extrema = sort_extrema (head_extrema);
00689
00690
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
00706
00707
00708
00709
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
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
00755
00756
00757 extrema * find_extrema
00758 (
00759 float * fvol,
00760 int num_nv,
00761 int nfirst,
00762 int nlast
00763 )
00764
00765 {
00766 int ix, jy, kz, ixyz;
00767 int it, ijk;
00768 int nx, ny, nz, nxy, nxyz;
00769 float fval;
00770 float * location = NULL;
00771 extrema * head_extrema = NULL;
00772 int passed_test;
00773
00774
00775
00776 nx = EX_nx; ny = EX_ny; nz = EX_nz;
00777 nxy = nx * ny; nxyz = nx*ny*nz;
00778
00779
00780
00781 for (ixyz = nfirst; ixyz < nlast; ixyz++)
00782 {
00783
00784
00785 if (!EX_mask[ixyz]) continue;
00786
00787
00788 fval = fvol[ixyz];
00789 if (fabs(fval) < EX_data_thr) continue;
00790
00791
00792 it = 0;
00793 passed_test = 1;
00794 while ((it < num_nv) && passed_test)
00795 {
00796 ijk = ixyz + EX_offset[it];
00797
00798
00799 if ((ijk < nfirst) || (ijk >= nlast))
00800 {
00801 if (EX_interior) passed_test = 0;
00802 }
00803
00804
00805 else if (!EX_mask[ijk])
00806 {
00807 if (EX_interior) passed_test = 0;
00808 }
00809
00810
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
00833 if (passed_test)
00834 {
00835 location = to_3dmm (ixyz);
00836 head_extrema = new_extrema (location, fval, head_extrema);
00837 }
00838 }
00839
00840
00841
00842 if (EX_sep_dist > 0.0) head_extrema = agglomerate_extrema (head_extrema,
00843 EX_sep_dist);
00844
00845
00846
00847 if (EX_sort) head_extrema = sort_extrema (head_extrema);
00848
00849
00850 return (head_extrema);
00851 }
00852
00853
00854
00855
00856
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
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
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
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
00997
00998
00999
01000 int EX_read_opts( int argc , char * argv[] )
01001 {
01002 int nopt = 1 ;
01003 char message[80];
01004
01005
01006 while( nopt < argc ){
01007
01008
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
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
01034 if( strcmp(argv[nopt],"-quiet") == 0 ){
01035 EX_quiet = 1;
01036 nopt++ ; continue ;
01037 }
01038
01039
01040
01041 if( strcmp(argv[nopt],"-minima") == 0 ){
01042 EX_maxima = 0;
01043 nopt++; continue;
01044 }
01045
01046
01047
01048 if( strcmp(argv[nopt],"-maxima") == 0 ){
01049 EX_maxima = 1;
01050 nopt++; continue;
01051 }
01052
01053
01054
01055 if( strcmp(argv[nopt],"-strict") == 0 ){
01056 EX_strict = 1;
01057 nopt++; continue;
01058 }
01059
01060
01061
01062 if( strcmp(argv[nopt],"-partial") == 0 ){
01063 EX_strict = 0;
01064 nopt++; continue;
01065 }
01066
01067
01068
01069 if( strcmp(argv[nopt],"-interior") == 0 ){
01070 EX_interior = 1;
01071 nopt++; continue;
01072 }
01073
01074
01075
01076 if( strcmp(argv[nopt],"-closure") == 0 ){
01077 EX_interior = 0;
01078 nopt++; continue;
01079 }
01080
01081
01082
01083 if( strcmp(argv[nopt],"-slice") == 0 ){
01084 EX_slice = 1;
01085 nopt++; continue;
01086 }
01087
01088
01089
01090 if( strcmp(argv[nopt],"-volume") == 0 ){
01091 EX_slice = 0;
01092 nopt++; continue;
01093 }
01094
01095
01096
01097 if( strcmp(argv[nopt],"-sort") == 0 ){
01098 EX_sort = 1;
01099 nopt++; continue;
01100 }
01101
01102
01103
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
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
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
01151 if( strcmp(argv[nopt],"-remove") == 0 ){
01152 EX_merge = 1;
01153 nopt++; continue;
01154 }
01155
01156
01157
01158 if( strcmp(argv[nopt],"-average") == 0 ){
01159 EX_merge = 2;
01160 nopt++; continue;
01161 }
01162
01163
01164
01165 if( strcmp(argv[nopt],"-weight") == 0 ){
01166 EX_merge = 3;
01167 nopt++; continue;
01168 }
01169
01170
01171
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
01189 if( argv[nopt][0] == '-' ){
01190 sprintf (message, " Unknown option: %s ", argv[nopt]);
01191 EX_error (message);
01192 }
01193
01194
01195
01196 break;
01197
01198
01199 }
01200
01201
01202
01203 if ((EX_merge == 3) && (EX_maxima == 0))
01204 EX_error ("-weight option is not compatible with -minima option");
01205
01206
01207
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
01221
01222
01223
01224 void * initialize_program (int argc, char * argv[], int * nopt)
01225 {
01226 const int MIN_NTHR = 10;
01227
01228 int iv;
01229 void * vfim = NULL;
01230 float * ffim = NULL;
01231 int ixyz, nthr;
01232 int nx, ny, nz, nxy, nxyz = 0;
01233 char message[80];
01234
01235
01236
01237 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
01238
01239
01240
01241 if( argc < 2 || strcmp(argv[1],"-help") == 0 ) EX_Syntax() ;
01242
01243
01244
01245 AFNI_logger (PROGRAM_NAME,argc,argv);
01246
01247
01248
01249 *nopt = EX_read_opts( argc , argv ) ;
01250
01251
01252
01253 THD_coorder_fill (my_getenv("AFNI_ORIENT"), &EX_cord);
01254
01255
01256
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
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
01280 ffim = (float *) malloc (sizeof(float) * nxyz); MTEST (ffim);
01281
01282
01283
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,
01288 MRI_float , ffim);
01289
01290
01291
01292 EX_mask = (byte *) malloc (sizeof(byte) * nxyz);
01293 MTEST (EX_mask);
01294
01295
01296
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
01319 if (ffim != NULL) { free (ffim); ffim = NULL; }
01320
01321
01322 THD_delete_3dim_dataset (EX_dset, False); EX_dset = NULL ;
01323
01324 }
01325
01326
01327
01328 {
01329
01330 if (argv[*nopt][0] == '-')
01331 EX_error ("ALL input options must precede ALL input datasets");
01332
01333
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
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
01353 if (EX_mask == NULL)
01354 {
01355
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
01366 EX_nx = nx; EX_ny = ny; EX_nz = nz;
01367 EX_nxy = nxy; EX_nxyz = nxyz;
01368
01369
01370
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
01387
01388
01389 void process_all_datasets (int argc, char * argv[], int nopt)
01390 {
01391 THD_3dim_dataset * input_dset=NULL;
01392
01393 int iv;
01394 void * vfim = NULL;
01395 float * ffim = NULL;
01396 int ibrick, nbricks;
01397 char message[80];
01398 int nx, ny, nz, nxy, nxyz;
01399 int kz, nfirst, nlast;
01400
01401
01402
01403 nx = EX_nx; ny = EX_ny; nz = EX_nz;
01404 nxy = nx * ny; nxyz = nx*ny*nz;
01405
01406
01407
01408 ffim = (float *) malloc (sizeof(float) * EX_nxyz); MTEST (ffim);
01409
01410
01411
01412 EX_nbricks = 0;
01413 while (nopt < argc)
01414 {
01415
01416 if (argv[nopt][0] == '-')
01417 EX_error ("ALL input options must precede ALL input datasets");
01418
01419
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
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
01440 nbricks = DSET_NVALS(input_dset);
01441
01442
01443
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,
01452 MRI_float , ffim);
01453
01454 if (EX_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
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
01475 EX_nbricks++;
01476 }
01477
01478
01479
01480 THD_delete_3dim_dataset (input_dset, False); input_dset = NULL ;
01481
01482 nopt++;
01483 }
01484
01485
01486
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
01500
01501
01502 void write_bucket ()
01503
01504 {
01505 THD_3dim_dataset * new_dset = NULL;
01506 int ixyz, nxyz;
01507 int kz, nz;
01508 extrema * head_extrema = NULL;
01509 byte ** bar = NULL;
01510 int nbricks;
01511 int ibrick;
01512 int ierror;
01513
01514 char message[80];
01515
01516
01517
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
01526 new_dset = EDIT_empty_copy (EX_dset);
01527
01528
01529
01530 tross_Copy_History (EX_dset, new_dset);
01531 if( commandline != NULL ) tross_Append_History( new_dset , commandline ) ;
01532
01533
01534
01535
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,
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
01564 bar = (byte **) malloc (sizeof(byte *) * nbricks);
01565 MTEST (bar);
01566
01567
01568
01569 for (ibrick = 0; ibrick < nbricks; ibrick++)
01570 {
01571
01572 bar[ibrick] = (byte *) malloc (sizeof(byte) * nxyz);
01573 MTEST (bar[ibrick]);
01574
01575
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
01591 EDIT_substitute_brick (new_dset, ibrick, MRI_byte, bar[ibrick]);
01592
01593 }
01594
01595
01596
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
01605 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01606
01607
01608 }
01609
01610
01611
01612
01613
01614
01615
01616 void output_results ()
01617 {
01618 int islice, ibrick;
01619 int nslices, nbricks;
01620
01621
01622
01623 nbricks = EX_nbricks;
01624 nslices = EX_nz;
01625
01626
01627
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
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
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
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
01682
01683 initialize_program (argc, argv, &nopt);
01684
01685
01686
01687 process_all_datasets (argc, argv, nopt);
01688
01689
01690
01691 output_results ();
01692
01693
01694 exit(0) ;
01695 }
01696
01697
01698
01699
01700