1 #define MAIN
2
3 #include "mrilib.h"
4 #include "afni.h"
5 #include <stdio.h>
6 #include <stdlib.h>
7
8 #define SEG_UNKNOWN 0
9 #define SEG_CSF 1
10 #define SEG_GM 2
11 #define SEG_WM 4
12
SegValToSegName(int s)13 const char *SegValToSegName(int s)
14 {
15 switch (s) {
16 case SEG_UNKNOWN:
17 RETURN("Unknown");
18 case SEG_CSF:
19 RETURN("CSF");
20 case SEG_WM:
21 RETURN("WM");
22 case SEG_GM:
23 RETURN("GM");
24 default:
25 RETURN("Smoking?");
26 }
27 }
28
29 #define SEG_METH_UNKNOWN 0
30 #define SEG_METH_DIST 1
31 #define SEG_METH_T 2
32
SegMethToMethName(int s)33 const char *SegMethToMethName(int s)
34 {
35 switch (s) {
36 case SEG_METH_UNKNOWN:
37 RETURN("Unknown");
38 case SEG_METH_DIST:
39 RETURN("Distance");
40 case SEG_METH_T:
41 RETURN("T");
42 default:
43 RETURN("Smoking?");
44 }
45 }
46
47 static int vn=0 ;
48
vstep_print(void)49 static void vstep_print(void)
50 {
51 static char xx[10] = "0123456789" ;
52 fprintf(stderr , "%c" , xx[vn%10] ) ;
53 if( vn%10 == 9) fprintf(stderr,".") ;
54 vn++ ;
55 }
56
Segment_usage(void)57 void Segment_usage(void)
58 {
59 int i = 0;
60
61 ENTRY("Segment_usage");
62
63 printf( "Usage: Segment -anat ANAT -csf CSF -wm WM -gm GM -mask mset\n"
64 ".\n"
65 "-debug DEBUG\n"
66 "-prefix PREFIX\n"
67 "-sphere_hood RAD default is 4.0\n"
68 "-T_met\n"
69 "-Dist_met\n"
70 "-max_loops\n"
71 "-use_global_estimates\n"
72 "\n");
73 EXRETURN;
74 }
75
76 typedef struct {
77 char *aset_name;
78 char *mset_name;
79 char *csfinit_name;
80 char *wminit_name;
81 char *gminit_name;
82 char *prefix;
83 THD_3dim_dataset *aset;
84 THD_3dim_dataset *mset;
85 THD_3dim_dataset *csfinit;
86 THD_3dim_dataset *wminit;
87 THD_3dim_dataset *gminit;
88 THD_3dim_dataset *oset;
89 int debug;
90 int idbg, jdbg, kdbg;
91 float na;
92 int DistMetric;
93 int N_loopmax;
94 int UseGlobal;
95 } SEG_OPTS;
96
Segment_ParseInput(char * argv[],int argc)97 SEG_OPTS *Segment_ParseInput (char *argv[], int argc)
98 {
99 static char FuncName[]={"Segment_ParseInput"};
100 SEG_OPTS *Opt=NULL;
101 int kar, i, ind, exists;
102 char *outname, cview[10];
103 int brk = 0;
104
105 ENTRY("Segment_ParseInput");
106
107 Opt = (SEG_OPTS *)malloc(sizeof(SEG_OPTS));
108
109 kar = 1;
110 Opt->aset_name = NULL;
111 Opt->mset_name = NULL;
112 Opt->csfinit_name = NULL;
113 Opt->wminit_name = NULL;
114 Opt->gminit_name = NULL;
115 Opt->prefix = NULL;
116 Opt->aset = NULL;
117 Opt->mset = NULL;
118 Opt->csfinit = NULL;
119 Opt->wminit = NULL;
120 Opt->gminit = NULL;
121 Opt->oset = NULL;
122 Opt->debug = 0;
123 Opt->idbg = Opt->kdbg = Opt->jdbg = -1;
124 Opt->na = 4.0;
125 Opt->DistMetric = SEG_METH_T;
126 Opt->N_loopmax = 10;
127 Opt->UseGlobal = 0; /* Uses global estimates of class medians,
128 to help when local samples of a particular
129 class are rare. Does not seem to help
130 much. The most critical remains a good
131 initialization with the startup masks*/
132 brk = 0;
133 while (kar < argc) { /* loop accross command ine options */
134 /*fprintf(stdout, "%s verbose: Parsing command line...\n", FuncName);*/
135 if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
136 Segment_usage();
137 exit (0);
138 }
139
140 #ifdef USE_TRACING
141 if( strncmp(argv[kar],"-trace",5) == 0 ){
142 DBG_trace = 1 ;
143 brk = 1 ;
144 }
145 if( strncmp(argv[kar],"-TRACE",5) == 0 ){
146 DBG_trace = 2 ;
147 brk = 1 ;
148 }
149 #endif
150
151 if (!brk && (strcmp(argv[kar], "-debug") == 0)) {
152 kar ++;
153 if (kar >= argc) {
154 fprintf (stderr, "need argument after -debug \n");
155 exit (1);
156 }
157 Opt->debug = atoi(argv[kar]);
158 brk = 1;
159 }
160
161 if (!brk && (strcmp(argv[kar], "-max_loops") == 0)) {
162 kar ++;
163 if (kar >= argc) {
164 fprintf (stderr, "need argument after -max_loops \n");
165 exit (1);
166 }
167 Opt->N_loopmax = atoi(argv[kar]);
168 brk = 1;
169 }
170
171 if (!brk && (strcmp(argv[kar], "-sphere_hood") == 0)) {
172 kar ++;
173 if (kar >= argc) {
174 fprintf (stderr, "need argument after -sphere_hood \n");
175 exit (1);
176 }
177 Opt->na = atof(argv[kar]);
178 brk = 1;
179 }
180
181 if( strcmp(argv[kar],"-T_met") == 0 ){
182 Opt->DistMetric = SEG_METH_T ;
183 brk = 1;
184 }
185
186 if( strcmp(argv[kar],"-Dist_met") == 0 ){
187 Opt->DistMetric = SEG_METH_DIST ;
188 brk = 1;
189 }
190
191 if( strcmp(argv[kar],"-use_global_estimates") == 0 ){
192 Opt->UseGlobal = 1 ;
193 brk = 1;
194 }
195
196 if (!brk && (strcmp(argv[kar], "-vox_debug") == 0)) {
197 kar ++;
198 if (kar+2 >= argc) {
199 fprintf (stderr, "need 3 arguments after -vox_debug \n");
200 exit (1);
201 }
202 Opt->idbg = atoi(argv[kar]); ++kar;
203 Opt->jdbg = atoi(argv[kar]); ++kar;
204 Opt->kdbg = atoi(argv[kar]);
205 brk = 1;
206 }
207
208 if (!brk && (strcmp(argv[kar], "-anat") == 0)) {
209 kar ++;
210 if (kar >= argc) {
211 fprintf (stderr, "need argument after -anat \n");
212 exit (1);
213 }
214 Opt->aset_name = argv[kar];
215 brk = 1;
216 }
217
218 if (!brk && (strcmp(argv[kar], "-csf") == 0)) {
219 kar ++;
220 if (kar >= argc) {
221 fprintf (stderr, "need argument after -csf \n");
222 exit (1);
223 }
224 Opt->csfinit_name = argv[kar];
225 brk = 1;
226 }
227
228 if (!brk && (strcmp(argv[kar], "-wm") == 0)) {
229 kar ++;
230 if (kar >= argc) {
231 fprintf (stderr, "need argument after -wm \n");
232 exit (1);
233 }
234 Opt->wminit_name = argv[kar];
235 brk = 1;
236 }
237
238 if (!brk && (strcmp(argv[kar], "-mask") == 0)) {
239 kar ++;
240 if (kar >= argc) {
241 fprintf (stderr, "need argument after -mask \n");
242 exit (1);
243 }
244 Opt->mset_name = argv[kar];
245 brk = 1;
246 }
247
248 if (!brk && (strcmp(argv[kar], "-gm") == 0)) {
249 kar ++;
250 if (kar >= argc) {
251 fprintf (stderr, "need argument after -gm \n");
252 exit (1);
253 }
254 Opt->gminit_name = argv[kar];
255 brk = 1;
256 }
257
258 if (!brk && (strcmp(argv[kar], "-prefix") == 0)) {
259 kar ++;
260 if (kar >= argc) {
261 fprintf (stderr, "need argument after -prefix \n");
262 exit (1);
263 }
264 Opt->prefix = (argv[kar]);
265 brk = 1;
266 }
267
268 if (!brk) {
269 fprintf (stderr,"Option %s not understood. \n"
270 "Try -help for usage\n", argv[kar]);
271 exit (1);
272 } else {
273 brk = 0;
274 kar ++;
275 }
276
277 }
278 if ( !Opt->gminit_name ||
279 !Opt->csfinit_name ||
280 !Opt->wminit_name ||
281 !Opt->aset_name) {
282 ERROR_exit("Missing input");
283 }
284
285 if (!Opt->prefix) Opt->prefix = "SegOut";
286
287 RETURN(Opt);
288 }
289
290 /*!
291 based on qmedmad_float, sped up for 3dSegment
292 */
293
qmedmad_float_3dSeg(int n,float * ar,float * med,float * mad)294 void qmedmad_float_3dSeg( int n, float *ar, float *med, float *mad )
295 {
296 float me=0.0f , ma=0.0f , *q ;
297 register int ii ;
298
299 if( med == NULL && mad == NULL || n <= 0 || ar == NULL ) return ;
300
301 #if 0 /* Protectionism */
302 q = (float *)malloc(sizeof(float)*n) ;
303 memcpy(q,ar,sizeof(float)*n) ; /* duplicate input array */
304 #else
305 q = ar;
306 #endif
307 me = qmed_float( n , q ) ; /* compute median */
308
309 if( mad != NULL && n > 1 ){
310 for( ii=0 ; ii < n ; ii++ ) /* subtract off median */
311 q[ii] = fabsf(q[ii]-me) ; /* (absolute deviation) */
312 ma = qmed_float( n , q ) ; /* MAD = median absolute deviation */
313 }
314
315 #if 0
316 free(q) ; /* 05 Nov 2001 - oops */
317 #endif
318
319 if( med != NULL ) *med = me ; /* 19 Aug 2005 - assign only */
320 if( mad != NULL ) *mad = ma ; /* if not NULL */
321 return ;
322 }
323
324
325 /*!
326 Based on mri_nstat, spedup for 3dSegment
327 */
mri_nstat_3dSeg(MRI_IMAGE * im,int * code_vec,int N_code,float * val_vec)328 int mri_nstat_3dSeg( MRI_IMAGE *im, int *code_vec, int N_code, float *val_vec )
329 {
330 MRI_IMAGE *fim ;
331 float *far , outval=0.0f, med = 0.0f, mad = 0.0f;
332 float var = -1.0f , cvar = -1.0f, sig = -1.0f, mean = -1.0f;
333 int npt, code , medmaded = 0;
334 int ic = 0;
335
336 if( im == NULL || im->nvox == 0 ) return outval ;
337
338 /* convert input to float format, if not already there */
339
340 if( im->kind != MRI_float ) fim = mri_to_float(im) ;
341 else fim = im ;
342 far = MRI_FLOAT_PTR(fim) ; /* array of values to statisticate */
343 npt = fim->nvox ; /* number of values */
344
345 for (ic=0; ic < N_code; ++ic) {
346 code = code_vec[ic];
347 outval = 0.0f;
348 switch( code ){
349
350 case NSTAT_NUM: outval = (float)npt ; break ; /* quite easy */
351
352 default:
353 case NSTAT_MEAN:{
354 register int ii ;
355 for( ii=0 ; ii < npt ; ii++ ) outval += far[ii] ;
356 outval /= npt ;
357 }
358 break ;
359 case NSTAT_NZNUM:{
360 register int ii ;
361 for( ii=0 ; ii < npt ; ii++ ) if (far[ii] != 0.0f) outval += 1 ;
362 }
363 break ;
364 case NSTAT_FNZNUM:{
365 register int ii ;
366 for( ii=0 ; ii < npt ; ii++ ) if (far[ii] != 0.0f) outval += 1 ;
367 outval /= npt;
368 }
369 break ;
370
371 case NSTAT_SIGMA: /* these 3 need the mean and variance sums */
372 case NSTAT_CVAR:
373 case NSTAT_VAR:{
374 register float mm,vv ; register int ii ;
375 if( npt == 1 ) break ; /* will return 0.0 */
376 for( mm=0.0,ii=0 ; ii < npt ; ii++ ) mm += far[ii] ;
377 mm /= npt ; mean = mm;
378 for( vv=0.0,ii=0 ; ii < npt ; ii++ ) vv += (far[ii]-mm)*(far[ii]-mm) ;
379 vv /= (npt-1) ;
380 if( code == NSTAT_SIGMA ) sig = outval = sqrt(vv) ;
381 else if( code == NSTAT_VAR ) var = outval = vv ;
382 else if( mm != 0.0f ) cvar = outval = sqrt(vv) / fabsf(mm) ;
383 }
384 break ;
385
386 case NSTAT_MEDIAN:
387 if (!medmaded) {
388 qmedmad_float_3dSeg( npt , far , &med , &mad ) ;
389 medmaded = 1;
390 }
391 outval = med;
392 break ;
393
394 case NSTAT_MAD:
395 if (!medmaded) {
396 qmedmad_float_3dSeg( npt , far , &med , &mad ) ;
397 medmaded = 1;
398 }
399 outval = mad;
400 break ;
401 case NSTAT_KURT:
402 outval = 0;
403 ERROR_message("Not supported");
404 break;
405 case NSTAT_P2SKEW:
406 /* Pearson's second skewness coefficient */
407 if (!medmaded) {
408 qmedmad_float_3dSeg( npt , far , &med , &mad ) ;
409 medmaded = 1;
410 }
411 if (sig < 0 ) {
412 register float mm,vv ; register int ii ;
413 if( npt == 1 ) break ; /* will return 0.0 */
414 for( mm=0.0,ii=0 ; ii < npt ; ii++ ) mm += far[ii] ;
415 mm /= npt ; mean = mm;
416 for( vv=0.0,ii=0 ; ii < npt ; ii++ )
417 vv += (far[ii]-mm)*(far[ii]-mm) ;
418 vv /= (npt-1) ;
419 sig = sqrt(vv) ;
420 }
421 if (sig) outval = 3.0 * (mean - med) / sig;
422 else outval = 0.0;
423 break ;
424
425 case NSTAT_MAX:{
426 register int ii ;
427 outval = far[0] ;
428 for( ii=1 ; ii < npt ; ii++ ) if( far[ii] > outval ) outval = far[ii] ;
429 }
430 break ;
431
432 case NSTAT_MIN:{
433 register int ii ;
434 outval = far[0] ;
435 for( ii=1 ; ii < npt ; ii++ ) if( far[ii] < outval ) outval = far[ii] ;
436 }
437 break ;
438
439 case NSTAT_ABSMAX:{
440 register int ii ; register float vv ;
441 outval = fabsf(far[0]) ;
442 for( ii=1 ; ii < npt ; ii++ ){
443 vv = fabsf(far[ii]) ; if( vv > outval ) outval = vv ;
444 }
445 }
446 break ;
447 }
448
449 val_vec[ic] = outval;
450 }
451
452 /* cleanup and exit */
453 if( fim != im ) mri_free(fim) ;
454 return (1) ;
455 }
456
457
Segment(SEG_OPTS * Opt)458 int Segment(SEG_OPTS *Opt)
459 {
460 MCW_cluster *nbhd=NULL ;
461 float dx , dy , dz ;
462 int nxyz, nx, ny, nz, vstep, ijk, ii, jj, kk;
463 byte *vval=NULL, *csfmask=NULL, *gmmask=NULL, *wmmask=NULL, *mmask=NULL;
464 float *aval = NULL;
465 MRI_IMAGE *nbim_csf=NULL, *nbim_wm=NULL, *nbim_gm=NULL;
466 float dist_csf, dist_wm, dist_gm, min_dist,
467 med_csf, med_wm, med_gm, mad_csf, mad_wm, mad_gm,
468 med_gm_global, mad_gm_global, med_wm_global, mad_wm_global,
469 med_csf_global, mad_csf_global;
470 float dist_mad_rat, dist_mad_rat_csf, dist_mad_rat_wm, dist_mad_rat_gm ;
471 float min_T, T_csf, T_gm, T_wm, vval_candidate;
472 byte bb=0;
473 int ShowVoxSketchy=0, many = 0, far = 0,
474 n_csf, n_gm, n_wm, n_m, iloop, iii, nnn;
475 int N_code =0, code_vec[20], nl_csf=0, nl_gm = 0, nl_wm = 0;
476 float val_vec[20], Bcsf=0.0f, Bwm=0.0f, Bgm = 0.0f;
477 float *buff=NULL, *buffc=NULL, sfac=0.0,
478 Exp_gm_frac=0.0, Exp_wm_frac=0.0, Exp_csf_frac=0.0;
479 ENTRY("Segment");
480
481 /* prep output */
482 Opt->oset = EDIT_empty_copy( Opt->aset ) ;
483 EDIT_dset_items( Opt->oset ,
484 ADN_nvals , 5 ,
485 ADN_datum_all , MRI_byte ,
486 ADN_brick_fac , NULL ,
487 ADN_prefix , Opt->prefix ,
488 ADN_none ) ;
489 nx = DSET_NX(Opt->aset);
490 ny = DSET_NY(Opt->aset);
491 nz = DSET_NZ(Opt->aset);
492 nxyz = DSET_NVOX(Opt->aset);
493
494 csfmask = THD_makemask( Opt->csfinit, 0 , 1.0 , 0.0 );
495 gmmask = THD_makemask( Opt->gminit , 0 , 1.0 , 0.0 );
496 wmmask = THD_makemask( Opt->wminit , 0 , 1.0 , 0.0 );
497 if (Opt->mset) mmask = THD_makemask( Opt->mset , 0 , 1.0 , 0.0 );
498 else mmask = THD_makemask( Opt->aset , 0 , 1.0 , 0.0 );
499 aval = (float *)malloc(sizeof(float)*nxyz);
500 vval = (byte *)calloc(nxyz, sizeof(byte));
501 if (!vval || !csfmask || !gmmask || !wmmask || !aval || !mmask) {
502 ERROR_exit("Failed to allocate for vval or masks");
503 }
504
505 EDIT_coerce_scale_type( DSET_NVOX(Opt->aset) ,
506 DSET_BRICK_FACTOR(Opt->aset,0) ,
507 DSET_BRICK_TYPE(Opt->aset,0),
508 DSET_ARRAY(Opt->aset, 0) , /* input */
509 MRI_float , aval ) ; /* output */
510
511 if( Opt->na < 0.0f ){ dx = dy = dz = 1.0f ; Opt->na = -Opt->na ; }
512 else { dx = fabsf(DSET_DX(Opt->aset)) ;
513 dy = fabsf(DSET_DY(Opt->aset)) ;
514 dz = fabsf(DSET_DZ(Opt->aset)) ; }
515
516 nbhd = MCW_spheremask( dx,dy,dz , Opt->na ) ;
517 if (Opt->debug) {
518 fprintf(stderr,"nbhd: %p\n"
519 "%d voxels.\n",
520 nbhd, nbhd->num_pt);
521 }
522
523
524 iloop = 0;
525 do {
526 n_csf = THD_countmask(nxyz, csfmask);
527 n_gm = THD_countmask(nxyz, gmmask );
528 n_wm = THD_countmask(nxyz, wmmask );
529 n_m = THD_countmask(nxyz, mmask );
530
531 if (Opt->debug) {
532 fprintf(stderr,"Pass %d, have:\n"
533 "%d voxels in csf,\n"
534 "%d voxels in gm,\n"
535 "%d voxels in wm\n"
536 "%.2f%% gm/wm\n"
537 "%d voxels in mask\n",
538 iloop, n_csf,
539 n_gm,
540 n_wm,
541 (float)n_gm/(float)n_wm*100.0f,
542 n_m);
543 }
544
545 vstep = (Opt->debug && nxyz > 99999) ? nxyz/50 : 0 ;
546 if( vstep ) fprintf(stderr,"++ voxel loop:") ;
547
548 SetSearchAboutMaskedVoxel(1); /* search the neighborhood of a voxel that
549 is itself not in the mask.
550 The default behaviour for
551 THD_get_dset_nbhd is to return a NULL
552 if the voxel ii, jj, kk itself is not
553 in the mask */
554 /* get global stats for 3 classes */
555 if (!buff) buff = (float *)calloc(sizeof(float), nxyz);
556 if (!buff) {
557 fprintf( stderr,"Failed to allocate for buff. Misere!");
558 return(0);
559 }
560 EDIT_coerce_scale_type( nxyz ,
561 DSET_BRICK_FACTOR(Opt->aset,0) ,
562 DSET_BRICK_TYPE(Opt->aset,0),
563 DSET_ARRAY(Opt->aset, 0) , /* input */
564 MRI_float , buff ) ; /* output */
565
566 /* get csf global stats */
567 med_csf_global=-1;
568 mad_csf_global=-1;
569 buffc = (float *)calloc(sizeof(float), n_csf);
570 if (!buffc) {
571 fprintf( stderr,"Failed to allocate for buffc. Misere!");
572 return(0);
573 }
574 nnn = 0;
575 for (iii=0; iii<nxyz;++iii) {
576 if (csfmask[iii]) { buffc[nnn] = buff[iii]; ++nnn; }
577 }
578 qmedmad_float_3dSeg( nnn , buffc , &med_csf_global , &mad_csf_global ) ;
579 free(buffc); buffc = NULL;
580
581 /* get wm global stats */
582 med_wm_global=-1;
583 mad_wm_global=-1;
584 buffc = (float *)calloc(sizeof(float), n_wm);
585 if (!buffc) {
586 fprintf( stderr,"Failed to allocate for buffc. Misere!");
587 return(0);
588 }
589 nnn = 0;
590 for (iii=0; iii<nxyz;++iii) {
591 if (wmmask[iii]) { buffc[nnn] = buff[iii]; ++nnn; }
592 }
593 qmedmad_float_3dSeg( nnn , buffc , &med_wm_global , &mad_wm_global ) ;
594 free(buffc); buffc = NULL;
595
596 /* get gm global stats */
597 med_gm_global=-1;
598 mad_gm_global=-1;
599 buffc = (float *)calloc(sizeof(float), n_gm);
600 if (!buffc) {
601 fprintf( stderr,"Failed to allocate for buffc. Misere!");
602 return(0);
603 }
604 nnn = 0;
605 for (iii=0; iii<nxyz;++iii) {
606 if (gmmask[iii]) { buffc[nnn] = buff[iii]; ++nnn; }
607 }
608 qmedmad_float_3dSeg( nnn , buffc , &med_gm_global , &mad_gm_global ) ;
609 free(buffc); buffc = NULL;
610
611 free(buff); buff=NULL;
612
613 if (Opt->debug) {
614 fprintf(stderr,"Global stats\n"
615 "CSF: (med,mad) = (%5.2f,%5.2f)\n"
616 " WM: (med,mad) = (%5.2f,%5.2f)\n"
617 " GM: (med,mad) = (%5.2f,%5.2f)\n"
618 ,
619 med_csf_global , mad_csf_global,
620 med_wm_global , mad_wm_global,
621 med_gm_global , mad_gm_global);
622 }
623
624 /* for each voxel in aset */
625 far = 0; many = 0;
626 code_vec[0] = NSTAT_MEDIAN;
627 code_vec[1] = NSTAT_MAD;
628 code_vec[2] = NSTAT_NUM;
629 N_code = 3;
630 for( ijk=kk=0 ; kk < nz ; kk++ ){
631 for( jj=0 ; jj < ny ; jj++ ){
632 for( ii=0 ; ii < nx ; ii++,ijk++ ){
633 if( vstep && ijk%vstep==vstep-1 ) vstep_print() ;
634 if ( mmask[ijk] &&
635 ( Opt->debug < 3 ||
636 ( Opt->debug > 2 && Opt->idbg == ii &&
637 Opt->jdbg == jj && Opt->kdbg == kk) ) ) {
638 nbim_csf = THD_get_dset_nbhd( Opt->aset , 0 ,
639 csfmask , ii,jj,kk , nbhd ) ;
640 nbim_wm = THD_get_dset_nbhd( Opt->aset , 0 ,
641 wmmask , ii,jj,kk , nbhd ) ;
642 nbim_gm = THD_get_dset_nbhd( Opt->aset , 0 ,
643 gmmask , ii,jj,kk , nbhd ) ;
644
645 dist_csf = dist_wm = dist_gm = -1.0;
646 if (nbim_csf) {
647 if (!mri_nstat_3dSeg( nbim_csf, code_vec, N_code, val_vec )) {
648 fprintf( stderr,
649 "Failed in mri_nstat_3dSeg @csf[%d %d %d]!\n",
650 ii, jj, kk);
651 }
652 med_csf = val_vec[0] ;
653 mad_csf = val_vec[1] ;
654 nl_csf = val_vec[2];
655 } else {
656 if (Opt->debug > 2)
657 fprintf(stderr,"NULL nbim_csf @[%d %d %d]\n", ii, jj, kk);
658 med_csf = -1.0;
659 mad_csf = -1.0;
660 }
661 if (nbim_wm) {
662 if (!mri_nstat_3dSeg( nbim_wm, code_vec, N_code, val_vec )) {
663 fprintf(stderr,"Failed in mri_nstat_3dSeg @wm[%d %d %d]!\n", ii, jj, kk);
664 }
665 med_wm = val_vec[0] ;
666 mad_wm = val_vec[1] ;
667 nl_wm = val_vec[2];
668 } else {
669 if (Opt->debug > 2)
670 fprintf(stderr,"NULL nbim_wm @[%d %d %d]\n", ii, jj, kk);
671 med_wm = -1.0;
672 mad_wm = -1.0;
673 }
674 if (nbim_gm) {
675 if (!mri_nstat_3dSeg( nbim_gm, code_vec, N_code, val_vec )) {
676 fprintf(stderr,"Failed in mri_nstat_3dSeg @gm[%d %d %d]!\n", ii, jj, kk);
677 }
678 med_gm = val_vec[0] ;
679 mad_gm = val_vec[1] ;
680 nl_gm = val_vec[2];
681 } else {
682 if (Opt->debug > 2)
683 fprintf(stderr,"NULL nbim_gm @[%d %d %d]\n", ii, jj, kk);
684 med_gm = -1.0;
685 mad_wm = -1.0;
686 }
687 if (Opt->UseGlobal) { /* adjust median values by giving weight
688 to global stats */
689 if (nl_csf > nl_wm > nl_gm) {
690 /* shading compensation comes from csf sample */
691 sfac = med_csf / med_csf_global;
692 } else if (nl_wm > nl_csf > nl_gm) {
693 /* shading compensation comes from wm sample */
694 sfac = med_wm / med_wm_global;
695 } else {
696 /* shading compensation comes from gm sample */
697 sfac = med_gm / med_gm_global;
698 }
699 } else {
700 sfac = 1.0f;
701 }
702
703 /* Now use previous calculations to get the distance
704 One of the classes, will get no correction because
705 sfac was derived from its class*/
706 if (nbim_csf) {
707 if (Opt->UseGlobal) {
708 Exp_csf_frac = 0.2; /* Expected fraction of csf
709 in neighborhood */
710 Bcsf = (float)nl_csf/((float)nbhd->num_pt*Exp_csf_frac);
711 if (Bcsf > 1.0) Bcsf = 1.0; /* only use B when too few
712 samples present*/
713 med_csf = med_csf*Bcsf +
714 med_csf_global*sfac*(1-Bcsf);
715 mad_csf = mad_csf*Bcsf +
716 mad_csf_global*sfac*(1-Bcsf);
717 }
718 dist_csf = med_csf - aval[ijk];
719 if (dist_csf < 0.0) dist_csf = -dist_csf;
720 }
721 if (nbim_gm) {
722 if (Opt->UseGlobal) {
723 Exp_gm_frac = 0.4;/* Expected fraction of gm
724 in neighborhood */
725 Bgm = (float)nl_gm/((float)nbhd->num_pt*Exp_gm_frac);
726 if (Bgm > 1.0) Bgm = 1.0; /* only use B when too few
727 samples present*/
728 med_gm = med_gm*Bgm +
729 med_gm_global*sfac*(1-Bgm);
730 mad_gm = mad_gm*Bgm +
731 mad_gm_global*sfac*(1-Bgm);
732 }
733 dist_gm = med_gm - aval[ijk];
734 if (dist_gm < 0.0) dist_gm = -dist_gm;
735 }
736 if (nbim_wm) {
737 if (Opt->UseGlobal) {
738 Exp_wm_frac = 0.4;/* Expected fraction of wm
739 in neighborhood */
740 Bwm = (float)nl_wm/((float)nbhd->num_pt*Exp_wm_frac);
741 if (Bwm > 1.0) Bwm = 1.0; /* only use B when too few
742 samples present*/
743 med_wm = med_wm*Bwm +
744 med_wm_global*sfac*(1-Bwm);
745 mad_wm = mad_wm*Bwm +
746 mad_wm_global*sfac*(1-Bwm);
747 }
748 dist_wm = med_wm - aval[ijk];
749 if (dist_wm < 0.0) dist_wm = -dist_wm;
750 }
751
752
753 if (Opt->DistMetric == SEG_METH_DIST) {/* distance based, what is
754 the closest to
755 Opt->aset[ijk] ? */
756 {
757 vval_candidate = SEG_UNKNOWN;
758 min_dist = 1e30;
759 dist_mad_rat = dist_mad_rat_csf =
760 dist_mad_rat_wm = dist_mad_rat_gm = 1e30;
761 if (dist_csf >= 0.0) {
762 dist_mad_rat_csf = dist_csf / mad_csf;
763 if (dist_csf < min_dist) {
764 vval_candidate = SEG_CSF;
765 min_dist = dist_csf;
766 dist_mad_rat = dist_mad_rat_csf; }
767 }
768 if (dist_gm >= 0.0) {
769 dist_mad_rat_gm = dist_gm / mad_gm;
770 if (dist_gm < min_dist) {
771 vval_candidate = SEG_GM ;
772 min_dist = dist_gm;
773 dist_mad_rat = dist_mad_rat_gm; }
774 }
775 if (dist_wm >= 0.0) {
776 dist_mad_rat_wm = dist_wm / mad_wm;
777 if (dist_wm < min_dist) {
778 vval_candidate = SEG_WM ;
779 min_dist = dist_wm;
780 dist_mad_rat = dist_mad_rat_wm; }
781 }
782 }
783 /* some checks to see if decision was somewhat shaky */
784 ShowVoxSketchy = 0;
785 if (dist_mad_rat > 2.0) { /* This decision is shaky,
786 value quite far from median */
787 ShowVoxSketchy = 1;
788 ++far;
789 } else { /* decision is OK, but
790 do we have other possible ones? */
791 #if 1
792 if (vval_candidate == SEG_CSF) {
793 if (dist_mad_rat_wm < 2.0 || dist_mad_rat_gm < 2.0) {
794 ShowVoxSketchy = 2;
795 ++many;
796 }
797 }else if (vval_candidate == SEG_GM) {
798 if (dist_mad_rat_csf < 2.0 || dist_mad_rat_wm < 2.0) {
799 ShowVoxSketchy = 2;
800 ++many;
801 }
802 }else if (vval_candidate == SEG_WM) {
803 if (dist_mad_rat_csf < 2.0 || dist_mad_rat_gm < 2.0) {
804 ShowVoxSketchy = 2;
805 ++many;
806 }
807 }
808 #else /* An alternative, but gives the same results. */
809 if (vval_candidate == SEG_CSF) {
810 if ( dist_mad_rat_wm < dist_mad_rat ||
811 dist_mad_rat_gm < dist_mad_rat) {
812 ShowVoxSketchy = 2;
813 ++many;
814 }
815 }else if (vval_candidate == SEG_GM) {
816 if ( dist_mad_rat_csf < dist_mad_rat ||
817 dist_mad_rat_wm < dist_mad_rat) {
818 ShowVoxSketchy = 2;
819 ++many;
820 }
821 }else if (vval_candidate == SEG_WM) {
822 if ( dist_mad_rat_csf < dist_mad_rat ||
823 dist_mad_rat_gm < dist_mad_rat) {
824 ShowVoxSketchy = 2;
825 ++many;
826 }
827 }
828 #endif
829 }
830
831 if (!ShowVoxSketchy) {
832 /* congrats, get that voxel out of the mask */
833 vval[ijk] = vval_candidate;
834 mmask[ijk] = 0;
835 } else {
836 if (iloop == Opt->N_loopmax) {
837 /* take what you can get */
838 vval[ijk] = vval_candidate;
839 /* don't set mmask[ijk] = 0;,
840 leave mmask's contents to flag sketchy voxels
841 mmask will reflect the choice made for sketchy
842 voxels*/
843 mmask[ijk] = vval_candidate;
844 }
845 }
846
847 if ( (Opt->debug > 1 && ShowVoxSketchy) ||
848 (Opt->idbg == ii && Opt->jdbg == jj && Opt->kdbg == kk)){
849 if (ShowVoxSketchy == 1)
850 fprintf(stdout,"\nSketchy, nothing close ");
851 else if (ShowVoxSketchy == 2)
852 fprintf(stdout,"\nSketchy, many close ");
853 else fprintf(stdout,"\nDebug ");
854 fprintf(stdout,"for voxel [%d %d %d], pass %d\n"
855 "aval = %.2f, sfac = %.4f\n"
856 "med:mad_csf = %.2f:%.2f, dist_csf = %.2f, "
857 "dist_mad_rat_csf = %.2f, Bcsf = %.4f"
858 "(%d) voxels in hood\n"
859 "med:mad_gm = %.2f:%.2f, dist_gm = %.2f, "
860 "dist_mad_rat_gm = %.2f, Bgm = %.4f"
861 "(%d) voxels in hood\n"
862 "med:mad_wm = %.2f:%.2f, dist_wm = %.2f, "
863 "dist_mad_rat_wm = %.2f, Bwm = %.4f"
864 "(%d) voxels in hood\n"
865 "Voxel set to %s\n",
866 ii, jj, kk, iloop,
867 aval[ijk], sfac,
868 med_csf, mad_csf, dist_csf,
869 (dist_mad_rat_csf < 5000 ?
870 dist_mad_rat_csf : -1.0), Bcsf,
871 ((nbim_csf) ? nbim_csf->nvox : 0),
872 med_gm , mad_gm, dist_gm,
873 (dist_mad_rat_gm < 5000 ?
874 dist_mad_rat_gm : -1.0), Bgm,
875 ((nbim_gm) ? nbim_gm->nvox : 0),
876 med_wm , mad_wm, dist_wm,
877 (dist_mad_rat_wm < 5000 ?
878 dist_mad_rat_wm : -1.0), Bwm,
879 ((nbim_wm) ? nbim_wm->nvox : 0),
880 SegValToSegName(vval[ijk]));
881 }
882 } else if (Opt->DistMetric == SEG_METH_T) {
883 /* T based, what is Opt->aset[ijk] closest to,
884 based on the T value ? (Forgive me Stat Gods)*/
885 vval[ijk] = SEG_UNKNOWN;
886 min_dist = 1e30; min_T = 1e30;
887 T_csf = T_wm = T_gm = 1e30;
888 if (dist_csf >= 0.0) {
889 T_csf = dist_csf / ( mad_csf / sqrt (nbim_csf->nvox) );
890 if (T_csf < min_T) { vval[ijk] = SEG_CSF; min_T = T_csf; }
891 }
892 if (dist_gm >= 0.0) {
893 T_gm = dist_gm / ( mad_gm / sqrt (nbim_gm->nvox) );
894 if (T_gm < min_T) { vval[ijk] = SEG_GM ; min_T = T_gm; }
895 }
896 if (dist_wm >= 0.0) {
897 T_wm = dist_wm / ( mad_wm / sqrt (nbim_wm->nvox) );
898 if (T_wm < min_T) { vval[ijk] = SEG_WM ; min_T = T_wm; }
899 }
900 if ((Opt->idbg == ii && Opt->jdbg == jj && Opt->kdbg == kk)) {
901 fprintf(stdout,"\nDebug ");
902 fprintf(stdout,"for voxel [%d %d %d]\n"
903 "aval = %.2f\n"
904 "med:mad_csf = %.2f:%.2f, dist_csf = %.2f, "
905 "T_csf = %.2f, (%d) voxels in hood\n"
906 "med:mad_gm = %.2f:%.2f, dist_gm = %.2f, "
907 "T_gm = %.2f, (%d) voxels in hood\n"
908 "med:mad_wm = %.2f:%.2f, dist_wm = %.2f, "
909 "T_wm = %.2f, (%d) voxels in hood\n"
910 "Voxel set to %s\n",
911 ii, jj, kk,
912 aval[ijk],
913 med_csf, mad_csf, dist_csf,
914 (T_csf < 5000 ? T_csf : -1.0),
915 ((nbim_csf) ? nbim_csf->nvox : 0),
916 med_gm , mad_gm, dist_gm,
917 (T_gm < 5000 ? T_gm : -1.0),
918 ((nbim_gm) ? nbim_gm->nvox : 0),
919 med_wm , mad_wm, dist_wm,
920 (T_wm < 5000 ? T_wm : -1.0),
921 ((nbim_wm) ? nbim_wm->nvox : 0),
922 SegValToSegName(vval[ijk]));
923 }
924 }
925 mri_free(nbim_csf) ;
926 mri_free(nbim_wm) ;
927 mri_free(nbim_gm) ;
928 }
929 }}}
930
931 if (Opt->debug) {
932 if (Opt->DistMetric == SEG_METH_DIST) {
933 fprintf(stderr,
934 "\nPass %d:\n"
935 "Out of %d voxels in mask\n"
936 " %d (%.2f%%) were within 2 MAD of only one class\n"
937 " %d (%.2f%%) were far\n"
938 " %d (%.2f%%) had many close options\n",
939 iloop, n_m,
940 n_m - (far + many),(float)(n_m - (far + many))/(float)n_m*100.0 ,
941 far ,(float)far /(float)n_m*100.0 ,
942 many ,(float)many /(float)n_m*100.0 );
943 } else if (Opt->DistMetric == SEG_METH_T) {
944
945 }
946 }
947 if (iloop == 0) {
948 /* reset masks from initial guess to voxels that were a sure thing */
949 for (ijk=0; ijk<nxyz; ++ijk) {
950 csfmask[ijk] = wmmask[ijk] = gmmask[ijk] = SEG_UNKNOWN;
951 }
952 }
953 /* Update masks to reflect current result */
954 for (ijk=0; ijk<nxyz; ++ijk) {
955 if (vval[ijk] == SEG_CSF) csfmask[ijk] = SEG_CSF;
956 else if (vval[ijk] == SEG_GM) gmmask[ijk] = SEG_GM;
957 else if (vval[ijk] == SEG_WM) wmmask[ijk] = SEG_WM;
958 }
959 ++iloop;
960 } while (iloop <= Opt->N_loopmax && Opt->DistMetric == SEG_METH_DIST);
961
962 /* report final count */
963 if (Opt->debug) {
964 n_csf = THD_countmask(nxyz, csfmask);
965 n_gm = THD_countmask(nxyz, gmmask );
966 n_wm = THD_countmask(nxyz, wmmask );
967 n_m = THD_countmask(nxyz, mmask );
968
969 fprintf(stderr,"Final result:\n"
970 "%d voxels in csf,\n"
971 "%d voxels in gm,\n"
972 "%d voxels in wm\n"
973 "%.2f%% gm/wm\n"
974 "%d voxels were still sketchy but decided upon anyway\n",
975 n_csf,
976 n_gm,
977 n_wm,
978 (float)n_gm/(float)n_wm*100.0f,
979 n_m);
980 }
981
982 EDIT_substitute_brick(Opt->oset, 0, MRI_byte, vval );
983 EDIT_dset_items (Opt->oset, ADN_brick_label_one + 0, "All_Masks", ADN_none);
984 EDIT_substitute_brick(Opt->oset, 1, MRI_byte, wmmask );
985 EDIT_dset_items (Opt->oset, ADN_brick_label_one + 1, "White", ADN_none);
986 EDIT_substitute_brick(Opt->oset, 2, MRI_byte, gmmask );
987 EDIT_dset_items (Opt->oset, ADN_brick_label_one + 2, "Gray", ADN_none);
988 EDIT_substitute_brick(Opt->oset, 3, MRI_byte, csfmask );
989 EDIT_dset_items (Opt->oset, ADN_brick_label_one + 3, "CSF", ADN_none);
990 EDIT_substitute_brick(Opt->oset, 4, MRI_byte, mmask );
991 EDIT_dset_items (Opt->oset, ADN_brick_label_one + 4, "Sketcht", ADN_none);
992
993 KILL_CLUSTER(nbhd); nbhd = NULL;
994 free(aval); aval = NULL;
995 if (0) { /* free no more, now being placed in oset ... */
996 free(csfmask); free(gmmask); free(wmmask);
997 csfmask = NULL; gmmask = NULL; wmmask = NULL;
998 free(mmask); mmask = NULL;
999 }
1000
1001 RETURN(1);
1002 }
1003
main(int argc,char ** argv)1004 int main(int argc, char **argv)
1005 {
1006 SEG_OPTS *Opt=NULL;
1007
1008 mainENTRY("Segment");
1009
1010 Opt = Segment_ParseInput (argv, argc);
1011
1012 /* load the input data */
1013 Opt->aset = THD_open_dataset( Opt->aset_name );
1014 if( !ISVALID_DSET(Opt->aset) ){
1015 fprintf(stderr,"**ERROR: can't open dataset %s\n",Opt->aset_name) ;
1016 exit(1);
1017 }
1018 /*
1019 Opt->mset = THD_open_dataset( Opt->mset_name );
1020 if( !ISVALID_DSET(Opt->mset) ){
1021 fprintf(stderr,"**ERROR: can't open dataset %s\n",Opt->mset_name) ;
1022 exit(1);
1023 }
1024 */
1025 Opt->wminit = THD_open_dataset( Opt->wminit_name );
1026 if( !ISVALID_DSET(Opt->wminit) ){
1027 fprintf(stderr,"**ERROR: can't open dataset %s\n",Opt->wminit_name) ;
1028 exit(1);
1029 }
1030 Opt->gminit = THD_open_dataset( Opt->gminit_name );
1031 if( !ISVALID_DSET(Opt->gminit) ){
1032 fprintf(stderr,"**ERROR: can't open dataset %s\n",Opt->gminit_name) ;
1033 exit(1);
1034 }
1035 Opt->csfinit = THD_open_dataset( Opt->csfinit_name );
1036 if( !ISVALID_DSET(Opt->gminit) ){
1037 fprintf(stderr,"**ERROR: can't open dataset %s\n",Opt->csfinit_name) ;
1038 exit(1);
1039 }
1040
1041 DSET_mallocize(Opt->aset) ; DSET_load(Opt->aset);
1042 DSET_mallocize(Opt->gminit) ; DSET_load(Opt->gminit);
1043 DSET_mallocize(Opt->wminit) ; DSET_load(Opt->wminit);
1044 DSET_mallocize(Opt->csfinit); DSET_load(Opt->csfinit);
1045
1046 /* check on sizes */
1047 if ( DSET_NX(Opt->aset) != DSET_NX(Opt->gminit)
1048 || DSET_NX(Opt->aset) != DSET_NX(Opt->wminit)
1049 || DSET_NX(Opt->aset) != DSET_NX(Opt->gminit)
1050 || DSET_NY(Opt->aset) != DSET_NY(Opt->gminit)
1051 || DSET_NY(Opt->aset) != DSET_NY(Opt->wminit)
1052 || DSET_NY(Opt->aset) != DSET_NY(Opt->gminit)
1053 || DSET_NZ(Opt->aset) != DSET_NZ(Opt->gminit)
1054 || DSET_NZ(Opt->aset) != DSET_NZ(Opt->wminit)
1055 || DSET_NZ(Opt->aset) != DSET_NZ(Opt->gminit) ) {
1056 ERROR_exit("All input data must have same grid");
1057 }
1058
1059 if (Opt->mset) { DSET_mallocize(Opt->mset); DSET_load(Opt->mset); }
1060 if (!Segment(Opt)) {
1061 ERROR_exit("Failed in Segment");
1062 }
1063
1064 /* write output */
1065 DSET_write(Opt->oset);
1066
1067 /* all done, free */
1068 DSET_delete(Opt->aset); Opt->aset = NULL;
1069 DSET_delete(Opt->mset); Opt->mset = NULL;
1070 DSET_delete(Opt->oset); Opt->oset = NULL;
1071 DSET_delete(Opt->gminit); Opt->gminit = NULL;
1072 DSET_delete(Opt->wminit); Opt->wminit = NULL;
1073 DSET_delete(Opt->csfinit); Opt->csfinit = NULL;
1074 free(Opt); Opt = NULL;
1075
1076 PRINT_COMPILE_DATE ; exit(0);
1077 }
1078