1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8 #include "uthash.h"
9
10 typedef struct {
11 int id; /* keep it named 'id' to facilitate use of convenience
12 macros in uthash . */
13 UT_hash_handle hh; /* keep name for same reason */
14 int index;
15 } INT_HASH_DATUM;
16
17 /*---------------------------------------------------------------------*/
18 /*! Make a byte mask from mask dataset:
19 miv = sub-brick of input
20 if( mask_bot <= mask_top ) then
21 only nonzero values in this range will be used
22 else
23 all nonzero values in the mask will be used
24 The input dataset should be byte-, short-, or float-valued.
25
26 The output is a byte array with 1s in "hit" locations and 0s in
27 other locations. The number of bytes is DSET_NVOX(mask_dset).
28 This array should be free()-d someday. If NULL is returned,
29 some grotesque error transpired.
30 -----------------------------------------------------------------------*/
31
THD_makemask(THD_3dim_dataset * mask_dset,int miv,float mask_bot,float mask_top)32 byte * THD_makemask( THD_3dim_dataset *mask_dset ,
33 int miv , float mask_bot , float mask_top )
34 {
35 float maxval ; /* for computing limits for an empty mask */
36 byte *mmm = NULL ;
37 int nvox , ii ;
38 int empty = 0 ; /* do we return an empty mask */
39
40 if( !ISVALID_DSET(mask_dset) ||
41 miv < 0 ||
42 miv >= DSET_NVALS(mask_dset) ) return NULL ;
43
44 nvox = DSET_NVOX(mask_dset) ;
45
46 DSET_load(mask_dset) ; if( !DSET_LOADED(mask_dset) ) return NULL ;
47
48 mmm = (byte *) calloc( sizeof(byte) * nvox , 1 ) ;
49
50 switch( DSET_BRICK_TYPE(mask_dset,miv) ){
51 default:
52 WARNING_message("makemask: bad brick type %d",
53 DSET_BRICK_TYPE(mask_dset,miv));
54 free(mmm) ; DSET_unload(mask_dset) ; return NULL ;
55
56 case MRI_short:{
57 short mbot , mtop ;
58 short *mar = (short *) DSET_ARRAY(mask_dset,miv) ;
59 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
60 if( mfac == 0.0 ) mfac = 1.0 ;
61 if( mask_bot <= mask_top ){
62 /* maybe this mask is empty, allow for rounding */
63 maxval = MRI_TYPE_maxval[MRI_short] + 0.5 ;
64 if( mask_bot/mfac >= maxval || mask_top/mfac <= -maxval ) empty=1;
65
66 mbot = SHORTIZE(mask_bot/mfac) ;
67 mtop = SHORTIZE(mask_top/mfac) ;
68 } else {
69 mbot = (short) -MRI_TYPE_maxval[MRI_short] ;
70 mtop = (short) MRI_TYPE_maxval[MRI_short] ;
71 }
72 if( !empty ) /* 6 Jun 2007 */
73 for( ii=0 ; ii < nvox ; ii++ )
74 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 )
75 mmm[ii]=1;
76 }
77 break ;
78
79 case MRI_byte:{
80 byte mbot , mtop ;
81 byte *mar = (byte *) DSET_ARRAY(mask_dset,miv) ;
82 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
83 if( mfac == 0.0 ) mfac = 1.0 ;
84 if( mask_bot <= mask_top && mask_top > 0.0 ){
85 /* maybe this mask is empty, allow for rounding */
86 /* (top <= 0 is flag for full mask) */
87 maxval = MRI_TYPE_maxval[MRI_byte] + 0.5 ;
88 if( mask_bot/mfac >= maxval ) empty = 1;
89
90 mbot = BYTEIZE(mask_bot/mfac) ;
91 mtop = BYTEIZE(mask_top/mfac) ;
92 } else {
93 mbot = 0 ;
94 // [PT: Dec 22, 2020] Thanks for the fix here, C Rorden!
95 mtop = (byte) MRI_TYPE_maxval[MRI_byte] ;
96 }
97 if( !empty ) /* 6 Jun 2007 */
98 for( ii=0 ; ii < nvox ; ii++ )
99 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 )
100 mmm[ii]=1;
101 }
102 break ;
103
104 case MRI_float:{
105 float mbot , mtop ;
106 float *mar = (float *) DSET_ARRAY(mask_dset,miv) ;
107 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
108 if( mfac == 0.0 ) mfac = 1.0 ;
109 if( mask_bot <= mask_top ){
110 mbot = (float) (mask_bot/mfac) ;
111 mtop = (float) (mask_top/mfac) ;
112 } else {
113 mbot = -WAY_BIG ;
114 mtop = WAY_BIG ;
115 }
116 for( ii=0 ; ii < nvox ; ii++ )
117 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 ) mmm[ii]=1;
118 }
119 break ;
120 }
121
122 return mmm ;
123 }
124
125 /*----------------------------------------------------------------------------*/
126 /*! Remove isolated voxels from a byte mask [21 May 2009 - RWCox]. */
127
THD_mask_remove_isolas(int nx,int ny,int nz,byte * mmm)128 int THD_mask_remove_isolas( int nx, int ny, int nz , byte *mmm )
129 {
130 int ii,jj,kk , ip,im , jp,jm , kp,km , qq,nxy , niso=0 ;
131
132 if( nx < 1 || ny < 1 || nz < 1 || mmm == NULL ) return 0 ;
133 nxy = nx*ny ;
134
135 for( qq=kk=0 ; kk < nz ; kk++ ){ /* qq = voxel index in 1D array */
136 km = kk-1 ; kp = kk+1 ;
137 for( jj=0 ; jj < ny ; jj++ ){
138 jm = jj-1 ; jp = jj+1 ;
139 for( ii=0 ; ii < nx ; ii++,qq++ ){
140 if( !mmm[qq] ) continue ; /* already 0 */
141 im = ii-1 ; ip = ii+1 ;
142 if( im >= 0 && mmm[qq-1] ) continue ; /* -x nbhr */
143 if( ip < nx && mmm[qq+1] ) continue ; /* +x */
144 if( jm >= 0 && mmm[qq-nx] ) continue ; /* -y */
145 if( jp < ny && mmm[qq+nx] ) continue ; /* +y */
146 if( km >= 0 && mmm[qq-nxy] ) continue ; /* -z */
147 if( kp < nz && mmm[qq+nxy] ) continue ; /* +z */
148 mmm[qq] = 0 ; niso++ ;
149 }}}
150 return niso ;
151 }
152
153 /*----------------------------------------------------------------------------*/
154 /*!
155 Similar to THD_makemask except that it turns the dset itself to mask values
156 returns (-1) if it fails, number of non-zero voxels if OK
157 */
158
THD_makedsetmask(THD_3dim_dataset * mask_dset,int miv,float mask_bot,float mask_top,byte * cmask)159 int THD_makedsetmask( THD_3dim_dataset *mask_dset ,
160 int miv , float mask_bot , float mask_top,
161 byte *cmask )
162 {
163 float maxval ; /* for computing limits for an empty mask */
164 int nvox , ii, nonzero=-1 , empty = 0 ;
165
166 if( !ISVALID_DSET(mask_dset) ||
167 miv < 0 ||
168 miv >= DSET_NVALS(mask_dset) ) return (-1) ;
169
170 nvox = DSET_NVOX(mask_dset) ;
171
172 DSET_mallocize(mask_dset); /* do this or else it could be a read only dset! */
173 DSET_load(mask_dset) ; if( !DSET_LOADED(mask_dset) ) return (-1) ;
174
175 nonzero = 0;
176 switch( DSET_BRICK_TYPE(mask_dset,miv) ){
177 default:
178 DSET_unload(mask_dset) ; return (-1) ;
179
180 case MRI_short:{
181 short mbot , mtop ;
182 short *mar = (short *) DSET_ARRAY(mask_dset,miv) ;
183 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
184 if( mfac == 0.0 ) mfac = 1.0 ;
185 if( mask_bot <= mask_top ){
186 /* maybe this mask is empty, allow for rounding */
187 maxval = MRI_TYPE_maxval[MRI_short] + 0.5 ;
188 if( mask_bot/mfac >= maxval || mask_top/mfac <= -maxval ) empty=1;
189
190 mbot = SHORTIZE(mask_bot/mfac) ;
191 mtop = SHORTIZE(mask_top/mfac) ;
192 } else {
193 mbot = (short) -MRI_TYPE_maxval[MRI_short] ;
194 mtop = (short) MRI_TYPE_maxval[MRI_short] ;
195 }
196 if (empty) { /* if empty, clear result 6 Jun 2007 */
197 for( ii=0 ; ii < nvox ; ii++ ) mar[ii] = 0;
198 } else if (cmask) {
199 for( ii=0 ; ii < nvox ; ii++ )
200 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0
201 && cmask[ii]) { mar[ii]=1; ++nonzero; }
202 else { mar[ii] = 0; }
203 } else {
204 for( ii=0 ; ii < nvox ; ii++ )
205 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 )
206 { mar[ii]=1; ++nonzero; }
207 else { mar[ii] = 0; }
208 }
209 }
210 break ;
211
212 case MRI_byte:{
213 byte mbot , mtop ;
214 byte *mar = (byte *) DSET_ARRAY(mask_dset,miv) ;
215 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
216 if( mfac == 0.0 ) mfac = 1.0 ;
217 if( mask_bot <= mask_top && mask_top > 0.0 ){
218 /* maybe this mask is empty, allow for rounding */
219 /* (top <= 0 is flag for full mask) */
220 maxval = MRI_TYPE_maxval[MRI_byte] + 0.5 ;
221 if( mask_bot/mfac >= maxval ) empty = 1;
222
223 mbot = BYTEIZE(mask_bot/mfac) ;
224 mtop = BYTEIZE(mask_top/mfac) ;
225 } else {
226 mbot = 0 ;
227 // [PT: Dec 22, 2020] Thanks for the fix here, C Rorden!
228 mtop = (byte) MRI_TYPE_maxval[MRI_byte] ;
229 }
230 if (empty) { /* if empty, clear result 6 Jun 2007 */
231 for( ii=0 ; ii < nvox ; ii++ ) mar[ii] = 0;
232 } else if (cmask) {
233 for( ii=0 ; ii < nvox ; ii++ )
234 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0
235 && cmask[ii]){ mar[ii]=1; ++nonzero; }
236 else { mar[ii] = 0; }
237 } else {
238 for( ii=0 ; ii < nvox ; ii++ )
239 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 )
240 { mar[ii]=1; ++nonzero; }
241 else { mar[ii] = 0; }
242 }
243 }
244 break ;
245
246 case MRI_float:{
247 float mbot , mtop ;
248 float *mar = (float *) DSET_ARRAY(mask_dset,miv) ;
249 float mfac = DSET_BRICK_FACTOR(mask_dset,miv) ;
250 if( mfac == 0.0 ) mfac = 1.0 ;
251 if( mask_bot <= mask_top ){
252 mbot = (float) (mask_bot/mfac) ;
253 mtop = (float) (mask_top/mfac) ;
254 } else {
255 mbot = -WAY_BIG ;
256 mtop = WAY_BIG ;
257 }
258 if (cmask) {
259 for( ii=0 ; ii < nvox ; ii++ )
260 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0
261 && cmask[ii]) { mar[ii]=1; ++nonzero; }
262 else { mar[ii] = 0; }
263 } else {
264 for( ii=0 ; ii < nvox ; ii++ )
265 if( mar[ii] >= mbot && mar[ii] <= mtop && mar[ii] != 0 )
266 { mar[ii]=1; ++nonzero; }
267 else { mar[ii] = 0; }
268 }
269 }
270 break ;
271 }
272
273 /* remove any scaling factor ZSS April 24 06*/
274 EDIT_BRICK_FACTOR(mask_dset,miv , 0.0);
275
276 return (nonzero) ;
277 }
278
279
280 /*---------------------------------------------------------------------------*/
281 /*! Convert an entire dataset of MRI_byte, and all values to binary
282 (set or not).
283
284 return 1 on failure, 0 on success 3 Jun, 2014 [rickr]
285 */
286
THD_dset_to_mask(THD_3dim_dataset * dset,float mask_bot,float mask_top)287 int THD_dset_to_mask(THD_3dim_dataset * dset, float mask_bot, float mask_top)
288 {
289 byte * bvol = NULL;
290 int ivol;
291
292 ENTRY("THD_dset_to_mask");
293
294 if( !ISVALID_DSET(dset) ) {
295 ERROR_message("dset_to_mask: dset not valid");
296 RETURN(1);
297 }
298
299 DSET_mallocize(dset); DSET_load(dset) ;
300 if( !DSET_LOADED(dset) ) {
301 ERROR_message("dset_to_mask: dset not loaded");
302 RETURN(1);
303 }
304
305 for( ivol = 0; ivol < DSET_NVALS(dset); ivol++ ) {
306
307 bvol = THD_makemask(dset, ivol, mask_bot, mask_top);
308 if( !bvol ) {
309 ERROR_message("dset_to_mask: failed to mask vol %d", ivol);
310 RETURN(1);
311 }
312
313 EDIT_substitute_brick(dset, ivol, MRI_byte, bvol);
314 EDIT_BRICK_FACTOR(dset, ivol, 0.0);
315 }
316
317 RETURN(0);
318 }
319
320 /*
321 Zero out voxels vv in dset where cmask[vv]=0
322 Returns the number of voxels edited in dset (across all sub-bricks)
323 -1 if dset was null
324 */
THD_applydsetmask(THD_3dim_dataset * dset,byte * cmask)325 int THD_applydsetmask( THD_3dim_dataset *dset , byte *cmask )
326 {
327 int ss, ii, jj, kk, vv, nedited = -1 ;
328
329 ENTRY("THD_applydsetmask");
330
331 if (!dset) RETURN(nedited);
332
333 if (!cmask) RETURN(0);
334
335 DSET_mallocize(dset); DSET_load(dset);
336 for (ss=0; ss<DSET_NVALS(dset); ++ss) {
337 switch (DSET_BRICK_TYPE(dset,ss)) {
338 case MRI_byte:
339 { byte *bv = (byte *)DSET_ARRAY(dset,ss) ;
340 vv=0;
341 for (kk=0; kk<DSET_NZ(dset); ++kk) {
342 for (jj=0; jj<DSET_NY(dset); ++jj) {
343 for (ii=0; ii<DSET_NX(dset); ++ii) {
344 if (!cmask[vv]) {
345 bv[vv] = 0;
346 ++nedited;
347 }
348 ++vv;
349 } } }
350 }
351 break;
352 case MRI_short:
353 { short *sv = (short *)DSET_ARRAY(dset,ss) ;
354 vv=0;
355 for (kk=0; kk<DSET_NZ(dset); ++kk) {
356 for (jj=0; jj<DSET_NY(dset); ++jj) {
357 for (ii=0; ii<DSET_NX(dset); ++ii) {
358 if (!cmask[vv]) {
359 sv[vv] = 0;
360 ++nedited;
361 }
362 ++vv;
363 } } }
364 }
365 break;
366 case MRI_float:
367 { float *fv = (float *)DSET_ARRAY(dset,ss) ;
368 vv=0;
369 for (kk=0; kk<DSET_NZ(dset); ++kk) {
370 for (jj=0; jj<DSET_NY(dset); ++jj) {
371 for (ii=0; ii<DSET_NX(dset); ++ii) {
372 if (!cmask[vv]) {
373 fv[vv] = 0;
374 ++nedited;
375 }
376 ++vv;
377 } } }
378 }
379 break;
380 case MRI_complex:
381 { complex *cv = (complex *)DSET_ARRAY(dset,ss) ;
382 vv=0;
383 for (kk=0; kk<DSET_NZ(dset); ++kk) {
384 for (jj=0; jj<DSET_NY(dset); ++jj) {
385 for (ii=0; ii<DSET_NX(dset); ++ii) {
386 if (!cmask[vv]) {
387 cv[vv].i = cv[vv].r = 0.0;
388 ++nedited;
389 }
390 ++vv;
391 } } }
392 }
393 break;
394 default:
395 ERROR_message(
396 "THD_applydsetmask: Dset type %d for subbrick %d not supported\n",
397 DSET_BRICK_TYPE(dset,ss), ss);
398 break;
399 }
400 }
401
402 RETURN(nedited);
403 }
404
405 /*----------------------------------------------------------------------------*/
406 extern int * UniqueInt (int *y, int ysz, int *kunq, int Sorted );
407
is_integral_sub_brick(THD_3dim_dataset * dset,int isb,int check_values)408 int is_integral_sub_brick ( THD_3dim_dataset *dset, int isb, int check_values)
409 {
410 float mfac = 0.0;
411 void *vv=NULL;
412
413 if( !ISVALID_DSET(dset) ||
414 isb < 0 ||
415 isb >= DSET_NVALS(dset) ) {
416
417 fprintf(stderr,"** Bad dset or sub-brick index.\n");
418 return (0) ;
419
420 }
421 if( !DSET_LOADED(dset) ) DSET_load(dset);
422
423 switch( DSET_BRICK_TYPE(dset,isb) ){
424 case MRI_short:
425 case MRI_byte:
426 if (check_values) {
427 mfac = DSET_BRICK_FACTOR(dset,isb) ;
428 if (mfac != 0.0f && mfac != 1.0f) return(0);
429 }
430 break;
431 case MRI_double:
432 case MRI_complex:
433 case MRI_float:
434 vv = (void *)DSET_ARRAY(dset,isb);
435 mfac = DSET_BRICK_FACTOR(dset,isb) ;
436 if (mfac != 0.0f && mfac != 1.0f) return(0);
437 if (!vv) {
438 fprintf(stderr,"** NULL array!\n");
439 return(0);
440 }
441 return(is_integral_data(DSET_NVOX(dset),
442 DSET_BRICK_TYPE(dset,isb),
443 DSET_ARRAY(dset,isb) ) );
444 break;
445 default:
446 return(0);
447 }
448
449 return(1);
450 }
451
is_integral_dset(THD_3dim_dataset * dset,int check_values)452 int is_integral_dset ( THD_3dim_dataset *dset, int check_values)
453 {
454 int i=0;
455
456 if( !ISVALID_DSET(dset) ) return(0);
457 for (i=0; i<DSET_NVALS(dset); ++i) {
458 if (!is_integral_sub_brick(dset, i, check_values)) return(0);
459 }
460 return(1);
461 }
462
463 /*!
464 Returns a list of the unique values in a dataset.
465 */
466
THD_unique_vals(THD_3dim_dataset * mask_dset,int miv,int * n_unique,byte * cmask)467 int *THD_unique_vals( THD_3dim_dataset *mask_dset ,
468 int miv,
469 int *n_unique ,
470 byte *cmask)
471 {
472 int nvox , ii, *unq = NULL, *vals=NULL;
473
474 *n_unique = 0;
475 unq = NULL ;
476
477 if( !ISVALID_DSET(mask_dset) ||
478 miv < 0 ||
479 miv >= DSET_NVALS(mask_dset) ) {
480
481 fprintf(stderr,"** Bad mask_dset or sub-brick index.\n");
482 return (unq) ;
483
484 }
485 nvox = DSET_NVOX(mask_dset) ;
486
487 DSET_load(mask_dset) ; if( !DSET_LOADED(mask_dset) ) return (unq) ;
488
489 if (!is_integral_sub_brick (mask_dset, miv, 0)) {
490 fprintf(stderr,"** Sub-brick %d of %s is not of an integral data type.\n",
491 miv, DSET_PREFIX(mask_dset) ? DSET_PREFIX(mask_dset):"NULL");
492 return (unq) ;
493 }
494
495 vals = (int *)malloc(sizeof(int)*nvox);
496 if (!vals) {
497 fprintf(stderr,"** Failed to allocate.\n");
498 return (unq) ;
499 }
500
501 switch( DSET_BRICK_TYPE(mask_dset,miv) ){
502 default:
503 fprintf(stderr,"** Bad dset type for unique operation.\n"
504 "Only integral valued dsets are allowed.\n");
505 DSET_unload(mask_dset) ; if (vals) free(vals); return (unq) ;
506
507 case MRI_short:{
508 short *mar = (short *) DSET_ARRAY(mask_dset,miv) ;
509 if (cmask) {
510 for( ii=0 ; ii < nvox ; ii++ )
511 if (cmask[ii]) vals[ii] = (int)(mar[ii]); else vals[ii] = 0;
512 } else {
513 for( ii=0 ; ii < nvox ; ii++ )
514 vals[ii] = (int)(mar[ii]);
515 }
516
517 }
518 break ;
519
520 case MRI_byte:{
521 byte *mar = (byte *) DSET_ARRAY(mask_dset,miv) ;
522 if (cmask) {
523 for( ii=0 ; ii < nvox ; ii++ )
524 if (cmask[ii]) vals[ii] = (int)(mar[ii]); else vals[ii] = 0;
525 } else {
526 for( ii=0 ; ii < nvox ; ii++ )
527 vals[ii] = (int)(mar[ii]);
528 }
529
530 }
531 break ;
532
533 case MRI_float:{
534 float *mar = (float *) DSET_ARRAY(mask_dset,miv) ;
535 if (cmask) {
536 for( ii=0 ; ii < nvox ; ii++ )
537 if (cmask[ii]) vals[ii] = (int)(mar[ii]); else vals[ii] = 0;
538 } else {
539 for( ii=0 ; ii < nvox ; ii++ )
540 vals[ii] = (int)(mar[ii]);
541 }
542
543 }
544 break ;
545 }
546
547 /* unique */
548 unq = UniqueInt (vals, nvox, n_unique, 0 );
549
550 free(vals); vals = NULL;
551
552 return (unq) ;
553 }
554
555 /*----------------------------------------------------------------------------*/
556 /* returns an nvox int array which represents
557 the rank of the voxel value in mask_dset
558 */
559
THD_unique_rank(THD_3dim_dataset * mask_dset,int miv,byte * cmask,char * mapname,int ** unqp,int * N_unq)560 int *THD_unique_rank( THD_3dim_dataset *mask_dset ,
561 int miv,
562 byte *cmask,
563 char *mapname,
564 int **unqp, int *N_unq)
565 {
566 int nvox , ii, *unq = NULL, *vals=NULL, imax=0;
567 INT_HASH_DATUM *rmap=NULL, *hd=NULL;
568 int n_unique, r;
569 FILE *fout=NULL;
570
571 n_unique = 0;
572 unq = NULL ;
573
574 if (unqp && *unqp!=NULL) {
575 fprintf(stderr,"** unqp (%p) not initialized properly to NULL", *unqp);
576 return (vals) ;
577 }
578
579 if( !ISVALID_DSET(mask_dset) ||
580 miv < 0 ||
581 miv >= DSET_NVALS(mask_dset) ) {
582
583 fprintf(stderr,"** Bad mask_dset or sub-brick index.\n");
584 return (vals) ;
585
586 }
587 nvox = DSET_NVOX(mask_dset) ;
588
589 DSET_load(mask_dset) ; if( !DSET_LOADED(mask_dset) ) return (vals) ;
590
591 if (!is_integral_sub_brick (mask_dset, miv, 0)) {
592 fprintf(stderr,"** Sub-brick %d of %s is not integral valued.\n",
593 miv, DSET_PREFIX(mask_dset) ? DSET_PREFIX(mask_dset):"NULL");
594 return (vals) ;
595 }
596
597
598 vals = (int *)malloc(sizeof(int)*nvox);
599 if (!vals) {
600 fprintf(stderr,"** Failed to allocate.\n");
601 return (vals) ;
602 }
603
604 switch( DSET_BRICK_TYPE(mask_dset,miv) ){
605 default:
606 fprintf( stderr,
607 "** Bad dset type for unique operation.\n"
608 "Only Byte, Short and float dsets are allowed.\n");
609 DSET_unload(mask_dset) ;
610 if (vals) free(vals); vals = NULL; return (vals) ;
611
612 case MRI_short:{
613 short *mar = (short *) DSET_ARRAY(mask_dset,miv) ;
614 if (cmask) {
615 for( ii=0 ; ii < nvox ; ii++ )
616 if (cmask[ii]) vals[ii] = (int)(mar[ii]); else vals[ii] = 0;
617 } else {
618 for( ii=0 ; ii < nvox ; ii++ )
619 vals[ii] = (int)(mar[ii]);
620 }
621
622 }
623 break ;
624
625 case MRI_byte:{
626 byte *mar = (byte *) DSET_ARRAY(mask_dset,miv) ;
627 if (cmask) {
628 for( ii=0 ; ii < nvox ; ii++ )
629 if (cmask[ii]) vals[ii] = (int)(mar[ii]); else vals[ii] = 0;
630 } else {
631 for( ii=0 ; ii < nvox ; ii++ )
632 vals[ii] = (int)(mar[ii]);
633 }
634
635 }
636 break ;
637
638 case MRI_float:{ /* not an integral type but we store ints (from NIFTI)
639 as floats */
640 float *mar = (float *) DSET_ARRAY(mask_dset,miv) ;
641 if (cmask) {
642 for( ii=0 ; ii < nvox ; ii++ )
643 if (cmask[ii]) vals[ii] = (int)(mar[ii]); else vals[ii] = 0;
644 } else {
645 for( ii=0 ; ii < nvox ; ii++ )
646 vals[ii] = (int)(mar[ii]);
647 }
648
649 }
650 break ;
651 }
652
653 /* unique */
654 unq = UniqueInt (vals, nvox, &n_unique, 0 );
655 /* fprintf(stderr,"-- Have %d unique values\n", n_unique); */
656 if (!unq) {
657 fprintf(stderr,"** Failed to create unique list\n");
658 free(vals); return (NULL);
659 }
660
661 if (mapname && mapname[0]) {
662 /*fprintf(stderr,"-- Writing mapping to >>%s<<\n", mapname);*/
663 if ((fout = fopen(mapname,"w"))) {
664 fprintf(fout, "#Rank Map (%d unique values)\n", n_unique);
665 fprintf(fout, "#Col. 0: Rank\n");
666 fprintf(fout, "#Col. 1: Input Dset Value\n");
667 }
668 }
669 /* now replace by rank */
670 #if 0
671 for (r=0; r<n_unique; ++r) {
672 /* fprintf(stderr,"-- Doing %d ...\n", unq[r]); */
673 if (cmask) {
674 for (ii=0; ii<nvox; ii++) {
675 if (cmask[ii]) {
676 if (vals[ii] == unq[r]) vals[ii] = r;
677 } else vals[ii] = 0;
678 }
679 } else {
680 for( ii=0 ; ii < nvox ; ii++ ) if (vals[ii] == unq[r]) vals[ii] = r;
681 }
682 }
683 #else /* faster approach */
684 imax=0;
685 for (r=0; r<n_unique; ++r) {
686 if (imax < unq[r]) imax = unq[r];
687 if (fout) fprintf(fout, "%d %d\n", r, unq[r]);
688 hd = (INT_HASH_DATUM*)calloc(1,sizeof(INT_HASH_DATUM));
689 hd->id = unq[r];
690 hd->index = r;
691 HASH_ADD_INT(rmap, id, hd);
692 }
693 for (ii=0; ii<nvox; ii++)
694 if (!cmask || cmask[ii]) {
695 HASH_FIND_INT(rmap,&(vals[ii]),hd);
696 if (hd) vals[ii] = hd->index;
697 else {
698 fprintf(stderr,
699 "** Failed to find key %d inhash table\n",
700 vals[ii]);
701 free(vals);
702 while (rmap) { hd=rmap; HASH_DEL(rmap,hd); free(hd); }
703 return (NULL);
704 }
705 }
706 /* destroy hash */
707 while (rmap) {
708 hd=rmap;
709 HASH_DEL(rmap,hd);
710 if (hd) free(hd);
711 }
712
713 #endif
714
715 if (!unqp) free(unq); else *unqp = unq; unq = NULL;
716 if (N_unq) *N_unq = n_unique;
717 if (fout) fclose(fout); fout = NULL;
718
719 return (vals) ;
720 }
721
722
723 /*----------------------------------------------------------------------------*/
724 /* Same as THD_unique_rank but replaces values in mask_dset with rank */
725
THD_unique_rank_edit(THD_3dim_dataset * mask_dset,int miv,byte * cmask,char * mapname,int ** unqp,int * N_unq)726 int THD_unique_rank_edit( THD_3dim_dataset *mask_dset ,
727 int miv,
728 byte *cmask,
729 char *mapname, int **unqp, int *N_unq)
730 {
731 int *vals=NULL, nvox, mxval, ii;
732
733 if (!(vals = THD_unique_rank(mask_dset, miv, cmask, mapname, unqp, N_unq))) {
734 fprintf(stderr,"** Failed to uniquate\n");
735 return (0);
736 }
737
738 mxval = -1;
739 nvox = DSET_NVOX(mask_dset) ;
740 for( ii=0 ; ii < nvox ; ii++ ) { if (vals[ii] > mxval) mxval = vals[ii]; }
741 /* fprintf (stderr,"-- Have maxval of %d\n", mxval); */
742
743 switch( DSET_BRICK_TYPE(mask_dset,miv) ){
744 default:
745 fprintf(stderr,"** Bad dset type for unique operation.\n"
746 "Should have been stopped a while ago.\n");
747 if (vals) free(vals); vals = NULL; return (0) ;
748
749 case MRI_short:{
750 short *mar = (short *) DSET_ARRAY(mask_dset,miv) ;
751 if (mxval > MRI_TYPE_maxval[MRI_short]) {
752 fprintf(stderr,
753 "** Have too many unique values (%d) for "
754 "datatype short (limit %f)!\n",
755 mxval, MRI_TYPE_maxval[MRI_short]);
756 if (vals) free(vals); vals = NULL; return (0) ;
757 }
758 EDIT_BRICK_FACTOR(mask_dset,miv,0.0);
759 for( ii=0 ; ii < nvox ; ii++ )
760 mar[ii] = (short)(vals[ii]);
761 }
762 break ;
763
764 case MRI_byte:{
765 byte *mar = (byte *) DSET_ARRAY(mask_dset,miv) ;
766 if (mxval > MRI_TYPE_maxval[MRI_byte]) {
767 fprintf(stderr,
768 "** Have too many unique values (%d) for "
769 "datatype byte (limit %f)!\n",
770 mxval, MRI_TYPE_maxval[MRI_byte]);
771 if (vals) free(vals); vals = NULL; return (0) ;
772 }
773 EDIT_BRICK_FACTOR(mask_dset,miv,0.0);
774 for( ii=0 ; ii < nvox ; ii++ )
775 mar[ii] = (byte)(vals[ii]);
776 }
777 break ;
778
779 case MRI_float:{
780 float *mar = (float *) DSET_ARRAY(mask_dset,miv) ;
781 EDIT_BRICK_FACTOR(mask_dset,miv,0.0);
782 for( ii=0 ; ii < nvox ; ii++ )
783 mar[ii] = (float)(vals[ii]);
784 }
785 break ;
786 }
787
788 return (1);
789
790 }
791
792 /*---------------------------------------------------------------------*/
793 /*! Count the number of nonzero voxels in a mask.
794 -----------------------------------------------------------------------*/
795
THD_countmask(int nvox,byte * mmm)796 int THD_countmask( int nvox , byte *mmm )
797 {
798 int ii,mc ;
799
800 if( nvox <= 0 || mmm == NULL ) return 0 ;
801
802 for( ii=mc=0 ; ii < nvox ; ii++ ) if( mmm[ii] ) mc++ ;
803
804 return mc ;
805 }
806
807 /*---------------------------------------------------------------------*/
808
THD_parse_boxball(int * boxball_num,float ** boxball_dat,char ** argv)809 int THD_parse_boxball( int *boxball_num , float **boxball_dat , char **argv )
810 {
811 int bnum , narg=0 ; float *bdat ;
812
813 if( boxball_num == NULL || boxball_dat == NULL || argv == NULL ) return 0 ;
814
815 bnum = *boxball_num ; if( bnum < 0 ) bnum = 0 ;
816 bdat = *boxball_dat ;
817
818 if( strcmp(argv[narg]+2,"box") == 0 ){
819 float xbot,xtop , ybot,ytop , zbot,ztop , btyp ; int nn ;
820 char code = *(argv[narg]+1) ; /* should be 'x', 'd' , 'n', or 'i' */
821 switch( code ){
822 case 'x': btyp = BOX_XYZ ; break ;
823 case 'd': btyp = BOX_DIC ; break ;
824 case 'n': btyp = BOX_NEU ; break ;
825 case 'i': btyp = BOX_IJK ; break ;
826 default: WARNING_message("Unknown 'box' option %s\n",argv[narg]) ; return 0 ;
827 }
828 nn = sscanf( argv[narg+1] , "%f:%f" , &xbot , &xtop ) ;
829 if( nn < 1 ){
830 WARNING_message("Can't decode %s after %s\n",argv[narg+1],argv[narg]); return 0 ;
831 }
832 else if( nn == 1 ) xtop=xbot ;
833 nn = sscanf( argv[narg+2] , "%f:%f" , &ybot , &ytop ) ;
834 if( nn < 1 ){
835 WARNING_message("Can't decode %s after %s\n",argv[narg+2],argv[narg]); return 0 ;
836 }
837 else if( nn == 1 ) ytop=ybot ;
838 nn = sscanf( argv[narg+3] , "%f:%f" , &zbot , &ztop ) ;
839 if( nn < 1 ){
840 WARNING_message("Can't decode %s after %s\n",argv[narg+3],argv[narg]); return 0 ;
841 }
842 else if( nn == 1 ) ztop=zbot ;
843 bdat = (float *) realloc( bdat , sizeof(float)*BOXLEN*(bnum+1) ) ;
844 bdat[0+BOXLEN*bnum] = btyp ;
845 bdat[1+BOXLEN*bnum] = xbot ;
846 bdat[2+BOXLEN*bnum] = xtop ;
847 bdat[3+BOXLEN*bnum] = ybot ;
848 bdat[4+BOXLEN*bnum] = ytop ;
849 bdat[5+BOXLEN*bnum] = zbot ;
850 bdat[6+BOXLEN*bnum] = ztop ;
851 bnum++ ; narg = 4 ;
852
853 } else if( strcmp(argv[narg]+2,"ball") == 0 ){
854 float xcen,ycen,zcen,rad , btyp ;
855 char code = *(argv[narg]+1) ; /* should be 'x', 'd' , or 'n' */
856 switch( code ){
857 case 'x': btyp = BALL_XYZ ; break ;
858 case 'd': btyp = BALL_DIC ; break ;
859 case 'n': btyp = BALL_NEU ; break ;
860 default: WARNING_message("Unknown 'ball' option %s",argv[narg]) ; return 0 ;
861 }
862 xcen = strtod( argv[narg+1] , NULL ) ;
863 ycen = strtod( argv[narg+2] , NULL ) ;
864 zcen = strtod( argv[narg+3] , NULL ) ;
865 rad = strtod( argv[narg+4] , NULL ) ;
866 if( rad <= 0.0f ){
867 WARNING_message("%s radius=%s !?",argv[narg],argv[narg+4]) ; rad = 0.0f;
868 }
869
870 bdat = (float *) realloc( bdat , sizeof(float)*BOXLEN*(bnum+1) ) ;
871 bdat[0+BOXLEN*bnum] = btyp ;
872 bdat[1+BOXLEN*bnum] = xcen ;
873 bdat[2+BOXLEN*bnum] = ycen ;
874 bdat[3+BOXLEN*bnum] = zcen ;
875 bdat[4+BOXLEN*bnum] = rad ;
876 bnum++ ; narg = 5 ;
877 }
878
879 *boxball_num = bnum ; *boxball_dat = bdat ; return narg ;
880 }
881
882 /*----------------------------------------------------------------------------*/
883
THD_boxballmask(THD_3dim_dataset * dset,int boxball_num,float * boxball_dat)884 byte * THD_boxballmask( THD_3dim_dataset *dset ,
885 int boxball_num , float *boxball_dat )
886 {
887 int nx,ny,nz , nxy,nxyz , ii,jj,kk ;
888 byte *bmask ;
889 int bb, ibot,itop, jbot,jtop, kbot,ktop , btyp ;
890 float xbot,xtop, ybot,ytop, zbot,ztop ;
891 float xcen,ycen,zcen , icen,jcen,kcen ;
892 float xmin,xmax , ymin,ymax , zmin,zmax , rad,dist , xx,yy,zz ;
893 THD_fvec3 dv,xv ;
894
895 ENTRY("THD_boxballmask") ;
896
897 if( !ISVALID_DSET(dset) || boxball_num <= 0 || boxball_dat == NULL ) RETURN(NULL) ;
898
899 xmin=dset->daxes->xxmin ; xmax=dset->daxes->xxmax ;
900 ymin=dset->daxes->yymin ; ymax=dset->daxes->yymax ;
901 zmin=dset->daxes->zzmin ; zmax=dset->daxes->zzmax ;
902
903 nx = DSET_NX(dset) ;
904 ny = DSET_NY(dset) ; nxy = nx*ny ;
905 nz = DSET_NZ(dset) ; nxyz = nxy*nz ;
906
907 bmask = (byte *)calloc(sizeof(byte),nxyz) ;
908
909 for( bb=0 ; bb < boxball_num ; bb++ ){
910
911 btyp = boxball_dat[0+BOXLEN*bb] ;
912
913 if( btyp < BALL_XYZ ){ /*---- box ----*/
914
915 xbot = boxball_dat[1+BOXLEN*bb]; xtop = boxball_dat[2+BOXLEN*bb];
916 ybot = boxball_dat[3+BOXLEN*bb]; ytop = boxball_dat[4+BOXLEN*bb];
917 zbot = boxball_dat[5+BOXLEN*bb]; ztop = boxball_dat[6+BOXLEN*bb];
918
919 if( btyp != BOX_IJK ){ /* convert coords to indexes */
920
921 if( btyp == BOX_NEU ){ /* coords from Neuroscience to DICOM */
922 xbot = -xbot; xtop = -xtop; ybot = -ybot; ytop = -ytop; btyp = BOX_DIC;
923 }
924 if( btyp == BOX_DIC ){ /* coords from DICOM to dataset */
925 LOAD_FVEC3(dv,xbot,ybot,zbot) ;
926 xv = THD_dicomm_to_3dmm( dset , dv ) ;
927 UNLOAD_FVEC3(xv,xbot,ybot,zbot) ;
928 LOAD_FVEC3(dv,xtop,ytop,ztop) ;
929 xv = THD_dicomm_to_3dmm( dset , dv ) ;
930 UNLOAD_FVEC3(xv,xtop,ytop,ztop) ;
931 }
932 if( xbot < xmin && xtop < xmin ) continue ; /* skip box if outside dataset */
933 if( xbot > xmax && xtop > xmax ) continue ;
934 if( ybot < ymin && ytop < ymin ) continue ;
935 if( ybot > ymax && ytop > ymax ) continue ;
936 if( zbot < zmin && ztop < zmin ) continue ;
937 if( zbot > zmax && ztop > zmax ) continue ;
938 LOAD_FVEC3(dv,xbot,ybot,zbot) ;
939 xv = THD_3dmm_to_3dfind( dset , dv ) ; /* coords from dataset to index */
940 UNLOAD_FVEC3(xv,xbot,ybot,zbot) ;
941 LOAD_FVEC3(dv,xtop,ytop,ztop) ;
942 xv = THD_3dmm_to_3dfind( dset , dv ) ;
943 UNLOAD_FVEC3(xv,xtop,ytop,ztop) ;
944 }
945 ibot = rint(xbot) ; jbot = rint(ybot) ; kbot = rint(zbot) ; /* round */
946 itop = rint(xtop) ; jtop = rint(ytop) ; ktop = rint(ztop) ;
947 if( ibot > itop ){ btyp = ibot; ibot = itop; itop = btyp; } /* flip? */
948 if( jbot > jtop ){ btyp = jbot; jbot = jtop; jtop = btyp; }
949 if( kbot > ktop ){ btyp = kbot; kbot = ktop; ktop = btyp; }
950
951 /* skip box if outside dataset */
952 if ( itop < 0 || ibot >= nx ) continue;
953 if ( jtop < 0 || jbot >= ny ) continue;
954 if ( ktop < 0 || kbot >= nz ) continue;
955
956 /* constrain values to dataset dimensions */
957 if ( ibot < 0 ) ibot = 0; if ( itop >= nx ) itop = nx-1;
958 if ( jbot < 0 ) jbot = 0; if ( jtop >= ny ) jtop = ny-1;
959 if ( kbot < 0 ) kbot = 0; if ( ktop >= nz ) ktop = nz-1;
960
961 for( kk=kbot ; kk <= ktop ; kk++ )
962 for( jj=jbot ; jj <= jtop ; jj++ )
963 for( ii=ibot ; ii <= itop ; ii++ ) bmask[ii+jj*nx+kk*nxy] = 1 ;
964
965 } else { /*---- ball ----*/
966
967 xcen = boxball_dat[1+BOXLEN*bb] ; ycen = boxball_dat[2+BOXLEN*bb] ;
968 zcen = boxball_dat[3+BOXLEN*bb] ; rad = boxball_dat[4+BOXLEN*bb] ;
969
970 /* convert center coords to dataset indexes */
971
972 if( btyp == BALL_NEU ){ /* coords from Neuroscience to DICOM */
973 xcen = -xcen; ycen = -ycen; btyp = BALL_DIC;
974 }
975 if( btyp == BALL_DIC ){ /* coords from DICOM to dataset */
976 LOAD_FVEC3(dv,xcen,ycen,zcen) ;
977 xv = THD_dicomm_to_3dmm( dset , dv ) ;
978 UNLOAD_FVEC3(xv,xcen,ycen,zcen) ;
979 }
980 if( xcen < xmin || xcen > xmax ) continue ; /* skip ball if outside */
981 if( ycen < ymin || ycen > ymax ) continue ;
982 if( zcen < zmin || zcen > zmax ) continue ;
983 LOAD_FVEC3(dv,xcen,ycen,zcen) ;
984 xv = THD_3dmm_to_3dfind( dset , dv ) ; /* coords from dataset to index */
985 UNLOAD_FVEC3(xv,icen,jcen,kcen) ;
986
987 ibot = rint(icen-rad) ; itop = rint(icen+rad) ; /* box around ball */
988 jbot = rint(jcen-rad) ; jtop = rint(jcen+rad) ;
989 kbot = rint(kcen-rad) ; ktop = rint(kcen+rad) ;
990
991 rad = rad*rad ;
992
993 for( kk=kbot ; kk <= ktop ; kk++ ){
994 for( jj=jbot ; jj <= jtop ; jj++ ){
995 for( ii=ibot ; ii <= itop ; ii++ ){
996 LOAD_FVEC3( dv , ii,jj,kk ) ; /* convert to xyz coords */
997 xv = THD_3dfind_to_3dmm( dset , dv ) ; /* then test distance^2 */
998 UNLOAD_FVEC3( xv , xx,yy,zz ) ; /* xyz of ball center. */
999 dist = SQR(xx-xcen) + SQR(yy-ycen) + SQR(zz-zcen) ;
1000 if( dist <= rad ) bmask[ii+jj*nx+kk*nxy] = 1 ;
1001 }}}
1002 }
1003
1004 } /*----- end of loop over box/ball list -----*/
1005
1006 RETURN(bmask) ;
1007 }
1008
1009 /*------- functions moved from vol2surf.c ------------- 13 Nov 2006 [rickr] */
1010
1011 /*----------------------------------------------------------------------
1012 * thd_mask_from_brick - create a mask from a sub-brick and threshold
1013 *
1014 * return the number of set voxels in the mask
1015 *----------------------------------------------------------------------
1016 */
1017
thd_mask_from_brick(THD_3dim_dataset * dset,int volume,float thresh,byte ** mask,int absolute)1018 int thd_mask_from_brick(THD_3dim_dataset * dset, int volume, float thresh,
1019 byte ** mask, int absolute)
1020 {
1021 float factor;
1022 byte * tmask;
1023 int nvox, type, c, size = 0;
1024
1025 ENTRY("thd_mask_from_brick");
1026
1027 if ( mask ) *mask = NULL; /* to be sure */
1028
1029 if ( !ISVALID_DSET(dset) || ! mask || volume < 0 )
1030 RETURN(-1);
1031
1032 if ( volume >= DSET_NVALS(dset) )
1033 {
1034 fprintf(stderr,"** tmfb: sub-brick %d out-of-range\n", volume);
1035 RETURN(-1);
1036 }
1037
1038 if( !DSET_LOADED(dset) ) DSET_load(dset);
1039 nvox = DSET_NVOX(dset);
1040 type = DSET_BRICK_TYPE(dset, volume);
1041
1042 if ( type != MRI_byte && type != MRI_short &&
1043 type != MRI_int && type != MRI_float )
1044 {
1045 fprintf(stderr,"** tmfb: invalid dataset type %s, sorry...\n",
1046 MRI_type_name[type]);
1047 RETURN(-1);
1048 }
1049
1050 tmask = (byte *)calloc(nvox, sizeof(byte));
1051 if ( ! tmask )
1052 {
1053 fprintf(stderr,"** tmfb: failed to allocate mask of %d bytes\n", nvox);
1054 RETURN(-1);
1055 }
1056
1057 factor = DSET_BRICK_FACTOR(dset, volume);
1058
1059 /* cheat: adjust threshold, not data */
1060 if ( factor != 0.0 ) thresh /= factor;
1061
1062 switch( DSET_BRICK_TYPE(dset, volume) )
1063 {
1064 case MRI_byte:
1065 {
1066 if (thresh <= (float)MRI_maxbyte) { /* ZSS: Oct 2011
1067 Without this test, a high threshold value might end up
1068 equal to MRI_maxbyte when BYTEIZED below, resulting in
1069 the highest voxel making it to the mask no matter how
1070 much higher the threshold is set. */
1071 byte * dp = DSET_ARRAY(dset, volume);
1072 byte thr = BYTEIZE(thresh + 0.99999); /* ceiling */
1073 for ( c = 0; c < nvox; c++ )
1074 if ( dp[c] != 0 && ( dp[c] >= thr ) )
1075 {
1076 size++;
1077 tmask[c] = 1;
1078 }
1079 }
1080 }
1081 break;
1082
1083 case MRI_short:
1084 {
1085 if (thresh <= (float)MRI_maxshort) { /* ZSS: Oct 2011 */
1086 short * dp = DSET_ARRAY(dset, volume);
1087 short thr = SHORTIZE(thresh + 0.99999); /* ceiling */
1088 for ( c = 0; c < nvox; c++, dp++ )
1089 if ( *dp != 0 && ( *dp >= thr || (absolute && *dp <= -thr) ) )
1090 {
1091 size++;
1092 tmask[c] = 1;
1093 }
1094 }
1095 }
1096 break;
1097
1098 case MRI_int:
1099 {
1100 int * dp = DSET_ARRAY(dset, volume);
1101 int thr = (int)(thresh + 0.99999); /* ceiling */
1102 for ( c = 0; c < nvox; c++, dp++ )
1103 if ( *dp != 0 && ( *dp >= thr || (absolute && *dp <= -thr) ) )
1104 {
1105 size++;
1106 tmask[c] = 1;
1107 }
1108 }
1109 break;
1110
1111 case MRI_float:
1112 {
1113 float * dp = DSET_ARRAY(dset, volume);
1114 for ( c = 0; c < nvox; c++, dp++ )
1115 if (*dp != 0 && (*dp >= thresh || (absolute && *dp <= -thresh)))
1116 {
1117 size++;
1118 tmask[c] = 1;
1119 }
1120 }
1121 break;
1122
1123 default: /* let's be sure */
1124 {
1125 fprintf(stderr,"** tmfb: invalid dataset type, sorry...\n");
1126 free(tmask);
1127 }
1128 break;
1129 }
1130
1131 *mask = tmask;
1132
1133 RETURN(size);
1134 }
1135
1136 /*----------------------------------------------------------------------
1137 * thd_multi_mask_from_brick - create a valued mask from a sub-brick
1138 *
1139 * return 0 on success, else failure 10 Nov 2006 [rickr]
1140 *----------------------------------------------------------------------
1141 */
1142
thd_multi_mask_from_brick(THD_3dim_dataset * dset,int volume,byte ** mask)1143 int thd_multi_mask_from_brick(THD_3dim_dataset * dset, int volume, byte ** mask)
1144 {
1145 float factor;
1146 byte * tmask;
1147 int nvox, type, c;
1148
1149 ENTRY("thd_multi_mask_from_brick");
1150
1151 if ( mask ) *mask = NULL; /* to be sure */
1152
1153 if ( !ISVALID_DSET(dset) || ! mask || volume < 0 )
1154 RETURN(-1);
1155
1156 if ( volume >= DSET_NVALS(dset) )
1157 {
1158 fprintf(stderr,"** tmmfb: sub-brick %d out-of-range\n", volume);
1159 RETURN(-1);
1160 }
1161
1162 if( !DSET_LOADED(dset) ) DSET_load(dset);
1163 nvox = DSET_NVOX(dset);
1164 type = DSET_BRICK_TYPE(dset, volume);
1165
1166 if ( type != MRI_byte && type != MRI_short &&
1167 type != MRI_int && type != MRI_float )
1168 {
1169 fprintf(stderr,"** tmmfb: invalid dataset type %s, sorry...\n",
1170 MRI_type_name[type]);
1171 RETURN(-1);
1172 }
1173
1174 tmask = (byte *)calloc(nvox, sizeof(byte));
1175 if ( ! tmask )
1176 {
1177 fprintf(stderr,"** tmmfb: failed to allocate mask of %d bytes\n", nvox);
1178 RETURN(-1);
1179 }
1180
1181 factor = DSET_BRICK_FACTOR(dset, volume);
1182 if( factor == 1.0 ) factor = 0.0;
1183
1184 switch( DSET_BRICK_TYPE(dset, volume) )
1185 {
1186 case MRI_byte:
1187 {
1188 byte * dp = DSET_ARRAY(dset, volume);
1189 if( factor )
1190 for ( c = 0; c < nvox; c++ )
1191 tmask[c] = (byte)(int)rint(dp[c]*factor);
1192 else
1193 for ( c = 0; c < nvox; c++ )
1194 tmask[c] = dp[c];
1195 }
1196 break;
1197
1198 case MRI_short:
1199 {
1200 short * dp = DSET_ARRAY(dset, volume);
1201 if( factor )
1202 for ( c = 0; c < nvox; c++ )
1203 tmask[c] = (byte)(int)rint(dp[c]*factor);
1204 else
1205 for ( c = 0; c < nvox; c++ )
1206 tmask[c] = (byte)dp[c];
1207 }
1208 break;
1209
1210 case MRI_int:
1211 {
1212 int * dp = DSET_ARRAY(dset, volume);
1213 if( factor )
1214 for ( c = 0; c < nvox; c++ )
1215 tmask[c] = (byte)(int)rint(dp[c]*factor);
1216 else
1217 for ( c = 0; c < nvox; c++ )
1218 tmask[c] = (byte)dp[c];
1219 }
1220 break;
1221
1222 case MRI_float:
1223 {
1224 float * dp = DSET_ARRAY(dset, volume);
1225 if( factor )
1226 for ( c = 0; c < nvox; c++ )
1227 tmask[c] = (byte)(int)rint(dp[c]*factor);
1228 else
1229 for ( c = 0; c < nvox; c++ )
1230 tmask[c] = (byte)dp[c];
1231 }
1232 break;
1233
1234 default: /* let's be sure */
1235 {
1236 fprintf(stderr,"** tmmfb: invalid dataset type, sorry...\n");
1237 free(tmask);
1238 }
1239 break;
1240 }
1241
1242 *mask = tmask;
1243
1244 RETURN(0);
1245 }
1246
1247 /****************************************************************************
1248 ** The functions below are for converting a byte-valued 0/1 mask to/from **
1249 ** an ASCII representation. The ASCII representation is formed like so: **
1250 ** 1. convert it to binary = 8 bits stored in each byte, rather than 1 **
1251 ** - this takes the mask from nvox bytes to 1+(nvox-1)/8 bytes **
1252 ** - this operation is done in function mask_binarize() [below] **
1253 ** 2. compress the binarized array with zlib **
1254 ** - this step is done in function array_to_zzb64(), which uses **
1255 ** function zz_compress_all() [in zfun.c] **
1256 ** 3. express the compressed array into Base64 notation (in ASCII) **
1257 ** - this step is also done in function array_to_zzb64(), which **
1258 ** uses function B64_to_base64() [in niml/niml_b64.c] **
1259 ** 4. attach at the end a string to indicate the number of voxels **
1260 ** in the mask. **
1261 ** + The above steps are done in function mask_to_b64string() [below]. **
1262 ** + The inverse is done in function mask_from_b64string() [below]. **
1263 ** + Function mask_b64string_nvox() [below] can be used to get the voxel **
1264 ** count from the end of the string, which can be used to check if a **
1265 ** mask is compatible with a given dataset for which it is intended. **
1266 ****************************************************************************
1267 * See program 3dMaskToASCII.c for sample usage of these functions. *
1268 *****************************************************************************/
1269
1270 /*-------------------------------------------------------------------------*/
1271 /*! Convert a byte-value 0/1 mask to an ASCII string in Base64. */
1272
mask_to_b64string(int nvox,byte * mful)1273 char * mask_to_b64string( int nvox , byte *mful )
1274 {
1275 byte *mbin ; char *str ; int nstr ;
1276
1277 if( nvox < 1 || mful == NULL ) return NULL ; /* bad inputs */
1278
1279 mbin = mask_binarize( nvox , mful ) ;
1280 str = array_to_zzb64( 1+(nvox-1)/8 , (char *)mbin , 72 ) ; free(mbin) ;
1281 if( str == NULL ) return NULL ; /* should never happen */
1282
1283 nstr = strlen(str) ;
1284 str = (char *)realloc( str , sizeof(char)*(nstr+16) ) ;
1285 sprintf( str+nstr-1 , "===%d" , nvox ) ; /* -1 to erase last linefeed */
1286
1287 return str ;
1288 }
1289
1290 /*-------------------------------------------------------------------------*/
1291 /*! Convert an ASCII string in Base64 to a byte-valued 0/1 mask. */
1292
mask_from_b64string(char * str,int * nvox)1293 byte * mask_from_b64string( char *str , int *nvox )
1294 {
1295 byte *mful , *mbin=NULL ; int nvvv , nstr , ii,ibot , nbin ;
1296
1297 if( str == NULL || nvox == NULL ) return NULL ; /* bad inputs */
1298
1299 nvvv = mask_b64string_nvox(str) ;
1300 if( nvvv <= 0 ) return NULL ; /* WTF? */
1301
1302 /* decode string to binarized array */
1303
1304 nbin = zzb64_to_array( str , (char **)&mbin ) ;
1305 if( nbin <= 0 || mbin == NULL ) return NULL ; /* bad decode */
1306
1307 /* decode binarized array to byte mask */
1308
1309 mful = mask_unbinarize( nvvv , mbin ) ; free(mbin) ;
1310
1311 *nvox = nvvv ; return mful ;
1312 }
1313
1314 /*-------------------------------------------------------------------------*/
1315
mask_b64string_nvox(char * str)1316 int mask_b64string_nvox( char *str )
1317 {
1318 int nstr , ii , ibot ;
1319
1320 if( str == NULL ) return 0 ;
1321 nstr = strlen(str) ; if( nstr < 7 ) return 0 ; /* too short */
1322
1323 /* find the last '=' at the end of the string */
1324
1325 ibot = nstr-16 ; if( ibot < 3 ) ibot = 3 ;
1326 for( ii=nstr-1 ; ii > ibot && str[ii] != '=' ; ii-- ) ; /*nada*/
1327 if( str[ii] != '=' ) return 0 ; /* badly formatted */
1328
1329 ibot = (int)strtod(str+ii+1,NULL) ; /* number of voxels */
1330 return ibot ;
1331 }
1332
1333 /*-------------------------------------------------------------------------*/
1334 /* Input: 0/1 array of bytes, nvox long.
1335 Output: compressed by a factor of 8: 1+(nvox-1)/8 bytes long.
1336 *//*-----------------------------------------------------------------------*/
1337
1338 static byte binar[8] = { 1 , 1<<1 , 1<<2 , 1<<3 , 1<<4 , 1<<5 , 1<<6 , 1<<7 } ;
1339
mask_binarize(int nvox,byte * mful)1340 byte * mask_binarize( int nvox , byte *mful )
1341 {
1342 register byte *mbin ; register int ii ;
1343
1344 if( nvox < 1 || mful == NULL ) return NULL ;
1345
1346 mbin = (byte *)calloc(sizeof(byte),1+(nvox-1)/8) ;
1347
1348 for( ii=0 ; ii < nvox ; ii++ )
1349 if( mful[ii] != 0 ) mbin[ii>>3] |= binar[ii&0x7] ;
1350
1351 return mbin ;
1352 }
1353
1354 /*-------------------------------------------------------------------------*/
1355
mask_unbinarize(int nvox,byte * mbin)1356 byte * mask_unbinarize( int nvox , byte *mbin )
1357 {
1358 register byte *mful ; register int ii ;
1359
1360 if( nvox < 1 || mbin == NULL ) return NULL ;
1361
1362 mful = (byte *)calloc(sizeof(byte),nvox) ;
1363
1364 for( ii=0 ; ii < nvox ; ii++ )
1365 mful[ii] = ( mbin[ii>>3] & binar[ii&0x7] ) != 0 ;
1366
1367 return mful ;
1368 }
1369
1370 /*===========================================================================*/
1371 /*! Create a binary byte-valued mask from an input string:
1372 - a dataset filename
1373 - a Base64 mask string
1374 - filename with data containing a Base64 mask string
1375 - future editions?
1376 *//*-------------------------------------------------------------------------*/
1377
THD_create_mask_from_string(char * str)1378 bytevec * THD_create_mask_from_string( char *str ) /* Jul 2010 */
1379 {
1380 bytevec *bvec=NULL ; int nstr ; char *buf=NULL ;
1381 int ferr=0;
1382
1383 ENTRY("THD_create_mask") ;
1384
1385 if( str == NULL || *str == '\0' ) RETURN(NULL) ;
1386
1387 nstr = strlen(str) ;
1388 bvec = (bytevec *)malloc(sizeof(bytevec)) ;
1389
1390 /* try to read it as a dataset */
1391
1392 if( nstr < THD_MAX_NAME ){
1393 THD_3dim_dataset *dset = THD_open_one_dataset(str) ;
1394 if( dset != NULL ){
1395 bvec->nar = DSET_NVOX(dset) ;
1396 bvec->ar = THD_makemask( dset , 0 , 1.0f,0.0f ) ;
1397 DSET_delete(dset) ;
1398 if( bvec->ar == NULL ){
1399 ERROR_message("Can't make mask from dataset '%s'",str) ;
1400 free(bvec) ; bvec = NULL ;
1401 }
1402 RETURN(bvec) ;
1403 }
1404
1405 ferr = 1; /* string is short, but failed to open as dataset */
1406 }
1407
1408 /* if str is a filename, read that file;
1409 otherwise, use the string itself to find the mask */
1410
1411 if( THD_is_file(str) ){
1412 buf = AFNI_suck_file(str) ;
1413 if( buf != NULL ) nstr = strlen(buf) ;
1414 } else {
1415 buf = str ;
1416 }
1417
1418 /* try to read buf as a Base64 mask string */
1419
1420 if( strrchr(buf,'=') != NULL ){
1421 int nvox ;
1422 bvec->ar = mask_from_b64string( buf , &nvox ) ;
1423 if( bvec->ar != NULL ){
1424 bvec->nar = nvox ;
1425 } else {
1426 /* might be a non-existent file 14 Jan 2014 [rickr] */
1427 if( ferr ) ERROR_message("Failed to open mask from '%s'", str);
1428 else ERROR_message("Can't make mask from string '%.16s' %s",
1429 buf,(nstr<=16)?" ":"...") ;
1430 free(bvec) ; bvec = NULL ;
1431 }
1432 } else {
1433 if( ferr ) ERROR_message("Failed to open mask '%s'", str);
1434 else ERROR_message("Don't understand mask string '%.16s'",
1435 buf,(nstr<=16)?" ":"...") ;
1436 free(bvec) ; bvec = NULL ;
1437 }
1438
1439 if( buf != str && buf != NULL ) free(buf) ;
1440 RETURN(bvec) ;
1441 }
1442