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
9 /*--------------------------------------------------------------------
10 routine to edit an input dataset in place according to inputs
11 in "edopt" (see editvol.h).
12
13 Feb 1996: This routine is much more complex now due to the need to deal
14 with byte, short, float, or complex data in sub-bricks.
15
16 30 Nov 1997: added ability to edit a given sub-brick, using the
17 edopt->iv_fim entry
18
19 17 June 1998: Modifications for erosion and dilation of clusters.
20 ----------------------------------------------------------------------*/
21
EDIT_one_dataset(THD_3dim_dataset * dset,EDIT_options * edopt)22 void EDIT_one_dataset( THD_3dim_dataset *dset , EDIT_options *edopt )
23 {
24 int edit_thtoin = edopt->thtoin ; /* copy into local variables */
25 int edit_noneg = edopt->noneg ; /* for historical reasons */
26 int edit_abs = edopt->abss ;
27 float edit_clip_bot = edopt->clip_bot ; /* Nov 1995: changed to floats */
28 float edit_clip_top = edopt->clip_top ;
29 float edit_thresh = edopt->thresh ;
30 float edit_thbot = edopt->thbot ; /* 26 Dec 2007 */
31 int edit_clust = edopt->edit_clust ; /* 10 Sept 1996 */
32 float clust_rmm = edopt->clust_rmm ;
33 float clust_vmul = edopt->clust_vmul ;
34 float erode_pv = edopt->erode_pv; /* 17 June 1998 */
35 int dilate = edopt->dilate; /* 17 June 1998 */
36 int filter_opt = edopt->filter_opt; /* 11 Sept 1996 */
37 float filter_rmm = edopt->filter_rmm; /* 11 Sept 1996 */
38 int thrfilter_opt = edopt->thrfilter_opt; /* 1 Oct 1996 */
39 float thrfilter_rmm = edopt->thrfilter_rmm; /* 1 Oct 1996 */
40 float edit_blur = edopt->blur ;
41 float edit_blurx = edopt->blurx ; /* 17 Jun 2019 */
42 float edit_blury = edopt->blury ;
43 float edit_blurz = edopt->blurz ;
44 float edit_thrblur = edopt->thrblur; /* 4 Oct 1996 */
45 int edit_scale = edopt->scale ;
46 float edit_mult = edopt->mult ;
47 int edit_zvol = edopt->do_zvol ;
48 int edit_ivfim = edopt->iv_fim ; /* 30 Nov 1997 */
49 int edit_ivthr = edopt->iv_thr ; /* 30 Nov 1997 */
50 int verbose = edopt->verbose ; /* 01 Nov 1999 */
51 int fake_dxyz = edopt->fake_dxyz ; /* 11 Sep 2000 */
52 int rank = edopt->rank; /* 13 Nov 2007 */
53 int edit_clip_unscaled = edopt->clip_unscaled ; /* 09 Aug 1996 */
54
55 THD_dataxes *daxes ;
56 short *sfim = NULL , *sthr = NULL ;
57 float *ffim = NULL , *fthr = NULL ;
58 complex *cfim = NULL ;
59 byte *bfim = NULL , *bthr = NULL ;
60 void *vfim = NULL , *vthr = NULL ;
61 int nx,ny,nz,nxy,nxyz , jj,kk , ptmin , iclu,nclu , fim_max ;
62 int iv_fim=0, iv_thr=0, fim_type=0, thr_type=0, icl, ijk ;
63 register int ii ;
64 float dx,dy,dz , dxyz , rmm,vmul , val , vvv ;
65 MCW_cluster_array *clar ;
66 MCW_cluster *blur=NULL ;
67 int fimtype , thrtype ;
68 float fimfac=0.0, thrfac=0.0;
69 int blur3d=0 ; /* 17 Jun 2019 */
70
71 /** get the data from this dataset **/
72
73 ENTRY("EDIT_one_dataset") ;
74
75 THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ;
76 THD_load_datablock( dset->dblk ) ;
77
78 if( !DSET_LOADED(dset) ){
79 fprintf(stderr,
80 "\n*** Cannot read data brick for dataset %s\a\n",
81 DSET_BRIKNAME(dset) ) ;
82 EXRETURN ;
83 }
84
85 /** load the data sub-brick indexes (iv_*) and check types for legality **/
86
87 if( ISANAT(dset) ){
88 if( edit_ivfim >= 0 && edit_ivfim < DSET_NVALS(dset) ) /* 30 Nov 1997 */
89 iv_fim = edit_ivfim ;
90 else
91 iv_fim = ANAT_ival_zero[dset->func_type] ;
92
93 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: fim index = %d\n",iv_fim) ;
94
95 fim_type = DSET_BRICK_TYPE(dset,iv_fim) ;
96 fimfac = DSET_BRICK_FACTOR(dset,iv_fim) ;
97 if( edit_ivthr >= 0 && edit_ivthr < DSET_NVALS(dset) ){ /* 26 Dec 2007 */
98 iv_thr = edit_ivthr ;
99 thr_type = DSET_BRICK_TYPE(dset,iv_thr) ;
100 thrfac = DSET_BRICK_FACTOR(dset,iv_thr) ;
101 if( thrfac == 0.0f ) thrfac = 1.0f ;
102 } else {
103 iv_thr = -1 ; /* old code didn't allow thresholding */
104 thr_type = ILLEGAL_TYPE ; /* of an anatomical type dataset! */
105 }
106
107 if( !AFNI_GOOD_DTYPE(fim_type) || fim_type == MRI_rgb ){
108 fprintf(stderr,"\n*** Illegal anatomy data type %s in dataset %s\a\n" ,
109 MRI_type_name[fim_type] ,
110 dset->dblk->diskptr->brick_name ) ;
111 EXRETURN ;
112 }
113
114 #ifdef AFNI_DEBUG
115 { char str[256] ;
116 sprintf(str,"Anat dset: iv=%d type=%s fac=%g",iv_fim,MRI_TYPE_name[fim_type],fimfac) ;
117 STATUS(str) ; }
118 #endif
119
120 }
121
122 if( ISFUNC(dset) ){
123 if( edit_ivfim >= 0 && edit_ivfim < DSET_NVALS(dset) ) /* 30 Nov 1997 */
124 iv_fim = edit_ivfim ;
125 else
126 iv_fim = FUNC_ival_fim[dset->func_type] ;
127
128 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: fim index = %d\n",iv_fim) ;
129
130 fim_type = DSET_BRICK_TYPE(dset,iv_fim) ;
131 fimfac = DSET_BRICK_FACTOR(dset,iv_fim) ;
132
133 if( edit_ivthr >= 0 && edit_ivthr < DSET_NVALS(dset) ) /* 30 Nov 1997 */
134 iv_thr = edit_ivthr ;
135 else
136 iv_thr = FUNC_ival_thr[dset->func_type] ;
137
138 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: thr index = %d\n",iv_thr) ;
139
140 if( iv_thr < 0 ){
141 thr_type = ILLEGAL_TYPE ;
142 thrfac = 0.0 ;
143 } else {
144 thr_type = DSET_BRICK_TYPE(dset,iv_thr) ;
145 thrfac = DSET_BRICK_FACTOR(dset,iv_thr) ;
146 if( thrfac == 0.0 ){
147 switch( thr_type ){
148 case MRI_short: thrfac = 1.0/FUNC_scale_short[dset->func_type]; break;
149 case MRI_byte : thrfac = 1.0/FUNC_scale_byte [dset->func_type]; break;
150 }
151 }
152 }
153
154 if( !AFNI_GOOD_FUNC_DTYPE(fim_type) || fim_type == MRI_rgb ){
155 fprintf(stderr,"\n*** Illegal functional data type %s in dataset %s\a\n" ,
156 MRI_type_name[fim_type], dset->dblk->diskptr->brick_name ) ;
157 EXRETURN ;
158 }
159
160 if( thr_type >= 0 && (!AFNI_GOOD_FUNC_DTYPE(thr_type) || fim_type == MRI_rgb) ){
161 fprintf(stderr,"\n*** Illegal threshold data type %s in dataset %s\a\n" ,
162 MRI_type_name[fim_type] , dset->dblk->diskptr->brick_name ) ;
163 EXRETURN ;
164 }
165
166 #ifdef AFNI_DEBUG
167 { char str[256] ;
168 sprintf(str,"Func dset: iv_fim=%d type=%s fac=%g",iv_fim,MRI_TYPE_name[fim_type],fimfac) ;
169 STATUS(str) ;
170 if( iv_thr >= 0 ){
171 sprintf(str,"Func dset: iv_thr=%d type=%s fac=%g",iv_thr,MRI_TYPE_name[thr_type],thrfac) ;
172 STATUS(str) ; } }
173 #endif
174
175 }
176
177 /** load the pointers to the sub-bricks **/
178
179 vfim = DSET_ARRAY(dset,iv_fim) ;
180 switch( fim_type ){
181 default:
182 fprintf(stderr,"\n*** Illegal data type in dataset %s\a\n",
183 dset->dblk->diskptr->brick_name ) ;
184 EXRETURN ;
185
186 case MRI_short: sfim = (short *) vfim ; break ;
187 case MRI_float: ffim = (float *) vfim ; break ;
188 case MRI_byte: bfim = (byte *) vfim ; break ;
189 case MRI_complex: cfim = (complex *) vfim ; break ;
190 }
191
192 if( iv_thr >= 0 ){
193 vthr = DSET_ARRAY(dset,iv_thr) ;
194 switch( thr_type ){
195 default:
196 fprintf(stderr,"\n*** Illegal thresh data type in dataset %s\a\n",
197 dset->dblk->diskptr->brick_name ) ;
198 EXRETURN ;
199
200 case MRI_short: sthr = (short *)vthr ; break ;
201 case MRI_float: fthr = (float *)vthr ; break ;
202 case MRI_byte: bthr = (byte *) vthr ; break ;
203 }
204 }
205
206 /** load the grid parameters **/
207
208 daxes = dset->daxes ;
209 nx = daxes->nxx ; dx = fabs(daxes->xxdel) ;
210 ny = daxes->nyy ; dy = fabs(daxes->yydel) ;
211 nz = daxes->nzz ; dz = fabs(daxes->zzdel) ;
212
213 if( fake_dxyz ) dx = dy = dz = 1.0 ; /* 11 Sep 2000 */
214
215 nxy = nx * ny ; nxyz = nxy * nz ; dxyz = dx*dy*dz ;
216
217 /*----- copy threshold over intensity? -----*/
218
219 STATUS("dataset loaded") ;
220
221 if( edit_thtoin && iv_thr >= 0 ){
222 float new_fimfac=0.0, scaling=0.0;
223
224 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: copy thr over fim\n") ;
225
226 /****
227 Find scaling factors for various conversions (0 --> no scaling)
228 scaling = factor to actually scale data by when copying to new brick
229 new_fimfac = factor to later scale data by when converting to floats
230 ****/
231
232 if( edit_thtoin == 2 ){
233 new_fimfac = scaling = 0.0 ; /** -2thtoin --> no scaling **/
234 } else {
235 switch( thr_type ){
236
237 /** threshold datum is shorts **/
238
239 case MRI_short:{
240 switch( fim_type ){
241 case MRI_short: /* fim datum is shorts --> no new scaling needed */
242 new_fimfac = thrfac ;
243 scaling = 0.0 ;
244 break ;
245
246 case MRI_float: /* fim datum is floats --> will be scaled properly */
247 new_fimfac = 0.0 ;
248 scaling = thrfac ;
249 break ;
250
251 case MRI_byte: /* fim datum is bytes */
252 new_fimfac = 1.0 / FUNC_scale_byte[dset->func_type] ;
253 scaling = thrfac * FUNC_scale_byte[dset->func_type] ;
254 break ;
255 }
256 }
257 break ;
258
259 /** threshold datum is bytes **/
260
261 case MRI_byte:{
262 switch( fim_type ){
263 case MRI_short: /* fim datum is shorts */
264 new_fimfac = 1.0 / FUNC_scale_short[dset->func_type] ;
265 scaling = thrfac * FUNC_scale_short[dset->func_type] ;
266 break ;
267
268 case MRI_float: /* fim datum is floats */
269 new_fimfac = 0.0 ;
270 scaling = thrfac ;
271 break ;
272
273 case MRI_byte: /* fim datum is bytes */
274 new_fimfac = thrfac ;
275 scaling = 0.0 ;
276 break ;
277 }
278 }
279 break ;
280
281 /** threshold datum is floats **/
282
283 case MRI_float:{
284 switch( fim_type ){
285 case MRI_short: /* fim datum is shorts */
286 new_fimfac = 1.0 / FUNC_scale_short[dset->func_type] ;
287 scaling = FUNC_scale_short[dset->func_type] ;
288 break ;
289
290 case MRI_float: /* fim datum is floats --> no scaling needed */
291 new_fimfac = 0.0 ;
292 scaling = 0.0 ;
293 break ;
294
295 case MRI_byte: /* fim datum is bytes */
296 new_fimfac = 1.0 / FUNC_scale_byte[dset->func_type] ;
297 scaling = FUNC_scale_byte[dset->func_type] ;
298 break ;
299 }
300 }
301 break ;
302 }
303 }
304
305 #ifdef AFNI_DEBUG
306 { char str[256] ;
307 sprintf(str,"thtoin: scaling=%f new_fimfac=%f input=%s output=%s",
308 scaling,new_fimfac,MRI_TYPE_name[thr_type],MRI_TYPE_name[fim_type]) ;
309 STATUS(str) ; }
310 #endif
311
312 /** have scaling factors, so use them **/
313
314 EDIT_coerce_scale_type( nxyz , scaling ,
315 thr_type , vthr , fim_type , vfim ) ;
316
317 DSET_BRICK_FACTOR(dset,iv_fim) = fimfac = new_fimfac ;
318 } /* end -1thtoin */
319
320 /*----- non-negative? -----*/
321
322 if( edit_noneg ){ /* meaningless for byte and complex */
323 STATUS("noneg") ;
324
325 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: remove negative values\n") ;
326
327 switch( fim_type ){
328 case MRI_short:
329 for( ii=0 ; ii < nxyz ; ii++ ) if( sfim[ii] < 0 ) sfim[ii] = 0 ;
330 break ;
331
332 case MRI_float:
333 for( ii=0 ; ii < nxyz ; ii++ ) if( ffim[ii] < 0 ) ffim[ii] = 0 ;
334 break ;
335
336 default:
337 STATUS("noneg applied to meaningless type: will be ignored") ;
338 }
339 }
340
341 /*----- absolute? -----*/
342
343 if( edit_abs ){ /* meaningless for byte */
344 STATUS("abs") ;
345
346 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: take absolute value\n") ;
347
348 switch( fim_type ){
349 case MRI_short:
350 for( ii=0 ; ii < nxyz ; ii++ ) if( sfim[ii] < 0 ) sfim[ii] = -sfim[ii] ;
351 break ;
352
353 case MRI_float:
354 for( ii=0 ; ii < nxyz ; ii++ ) if( ffim[ii] < 0 ) ffim[ii] = -ffim[ii] ;
355 break ;
356
357 case MRI_complex:
358 for( ii=0 ; ii < nxyz ; ii++ ){
359 cfim[ii].r = CABS(cfim[ii]) ; cfim[ii].i = 0.0 ;
360 }
361 break ;
362
363 default:
364 STATUS("abs applied to meaningless type: will be ignored") ;
365 }
366 }
367
368 /*----- clip? -----*/
369
370 if( edit_clip_bot < edit_clip_top ){
371
372 if( verbose )
373 fprintf(stderr,"--- EDIT_one_dataset: clip fim values %f %f\n",
374 edit_clip_bot, edit_clip_top) ;
375
376 switch( fim_type ){
377 case MRI_short:{
378 int top , bot ;
379 float ftop,fbot ;
380 if( fimfac > 0.0 && ! edit_clip_unscaled ){
381 ftop = edit_clip_top / fimfac ;
382 fbot = edit_clip_bot / fimfac ;
383 } else {
384 ftop = edit_clip_top ;
385 fbot = edit_clip_bot ;
386 }
387
388 top = rint(ftop) ; /* this code was modifed 28 Sep 1998 */
389 if( top >= MRI_maxshort ) top = MRI_maxshort + 1 ;
390 if( top <= -MRI_maxshort ) top = -(MRI_maxshort + 1) ;
391
392 bot = rint(fbot) ;
393 if( bot >= MRI_maxshort ) bot = MRI_maxshort + 1 ;
394 if( bot <= -MRI_maxshort ) bot = -(MRI_maxshort + 1) ;
395
396 #ifdef AFNI_DEBUG
397 { char str[256] ;
398 sprintf(str,"clipping short from %d to %d",bot,top) ;
399 STATUS(str) ; }
400 #endif
401 for( ii=0 ; ii < nxyz ; ii++ )
402 if( sfim[ii] > bot && sfim[ii] < top ) sfim[ii] = 0 ;
403 }
404 break ;
405
406 case MRI_byte:{
407 int top , bot ;
408 float ftop,fbot ;
409 if( fimfac > 0.0 && ! edit_clip_unscaled ){
410 ftop = edit_clip_top / fimfac ;
411 fbot = edit_clip_bot / fimfac ;
412 } else {
413 ftop = edit_clip_top ;
414 fbot = edit_clip_bot ;
415 }
416
417 top = rint(ftop) ;
418 if( top >= MRI_maxbyte ) top = MRI_maxbyte + 1 ;
419 if( top <= -MRI_maxbyte ) top = -(MRI_maxbyte + 1) ;
420
421 bot = rint(fbot) ;
422 if( bot >= MRI_maxbyte ) bot = MRI_maxbyte + 1 ;
423 if( bot <= -MRI_maxbyte ) bot = -(MRI_maxbyte + 1) ;
424
425 if( bot < 0 ) bot = 0 ;
426 if( top < bot ) top = bot ;
427 #ifdef AFNI_DEBUG
428 { char str[256] ;
429 sprintf(str,"clipping byte from %d to %d",bot,top) ;
430 STATUS(str) ; }
431 #endif
432 for( ii=0 ; ii < nxyz ; ii++ )
433 if( bfim[ii] > bot && bfim[ii] < top ) bfim[ii] = 0 ;
434 }
435 break ;
436
437 case MRI_float:{
438 float top , bot ;
439 if( fimfac > 0.0 && ! edit_clip_unscaled ){
440 top = edit_clip_top / fimfac ;
441 bot = edit_clip_bot / fimfac ;
442 } else {
443 top = edit_clip_top ;
444 bot = edit_clip_bot ;
445 }
446 #ifdef AFNI_DEBUG
447 { char str[256] ;
448 sprintf(str,"clipping float from %g to %g",bot,top) ;
449 STATUS(str) ; }
450 #endif
451 for( ii=0 ; ii < nxyz ; ii++ )
452 if( ffim[ii] > bot && ffim[ii] < top ) ffim[ii] = 0.0 ;
453 }
454 break ;
455
456 case MRI_complex:{
457 float val ;
458 float top , bot ;
459 if( fimfac > 0.0 && ! edit_clip_unscaled ){
460 top = edit_clip_top / fimfac ;
461 bot = edit_clip_bot / fimfac ;
462 } else {
463 top = edit_clip_top ;
464 bot = edit_clip_bot ;
465 }
466 #ifdef AFNI_DEBUG
467 { char str[256] ;
468 sprintf(str,"clipping complex from %g to %g",bot,top) ;
469 STATUS(str) ; }
470 #endif
471 for( ii=0 ; ii < nxyz ; ii++ ){
472 val = CABS(cfim[ii]) ;
473 if( val > bot && val < top ) cfim[ii].r = cfim[ii].i = 0.0 ;
474 }
475 }
476 break ;
477 }
478 }
479
480 /*----- apply threshold? -----*/
481
482 if( iv_thr >= 0 && (edit_thresh > 0.0 || edit_thresh > edit_thbot) ){
483
484 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: apply threshold\n") ;
485
486 switch( thr_type ){
487
488 /** threshold datum is shorts **/
489
490 case MRI_short:{
491 float thrplu , thrmin ;
492 float fplu = edit_thresh / thrfac ;
493 float fmin = edit_thbot / thrfac ;
494 thrplu = THRESH_SHORTIZE(fplu) ; thrmin = THRESH_SHORTIZE(fmin) ;
495 switch( fim_type ){
496 case MRI_short: /* fim datum is shorts */
497 for( ii=0 ; ii < nxyz ; ii++ )
498 if( sthr[ii] < thrplu && sthr[ii] > thrmin ){ sfim[ii] = 0; }
499 break ;
500
501 case MRI_byte: /* fim datum is bytes */
502 for( ii=0 ; ii < nxyz ; ii++ )
503 if( sthr[ii] < thrplu && sthr[ii] > thrmin ){ bfim[ii] = 0; }
504 break ;
505
506 case MRI_float: /* fim datum is floats */
507 for( ii=0 ; ii < nxyz ; ii++ )
508 if( sthr[ii] < thrplu && sthr[ii] > thrmin ){ ffim[ii] = 0.0; }
509 break ;
510
511 case MRI_complex: /* fim datum is complex */
512 for( ii=0 ; ii < nxyz ; ii++ )
513 if( sthr[ii] < thrplu && sthr[ii] > thrmin ){
514 cfim[ii].r = cfim[ii].i = 0.0 ;
515 }
516 break ;
517 }
518 }
519 break ;
520
521 /** threshold datum is bytes **/
522
523 case MRI_byte:{
524 byte thrplu , thrmin ;
525 float fplu = edit_thresh / thrfac ;
526 float fmin = edit_thbot / thrfac ;
527 thrplu = BYTEIZE(fplu) ; thrmin = BYTEIZE(fmin) ;
528 switch( fim_type ){
529 case MRI_short: /* fim datum is shorts */
530 for( ii=0 ; ii < nxyz ; ii++ )
531 if( bthr[ii] < thrplu && (thrmin == 0 || bthr[ii] > thrmin) ){
532 sfim[ii] = 0;
533 }
534 break ;
535
536 case MRI_byte: /* fim datum is bytes */
537 for( ii=0 ; ii < nxyz ; ii++ )
538 if( bthr[ii] < thrplu && (thrmin == 0 || bthr[ii] > thrmin) ){
539 bfim[ii] = 0;
540 }
541 break ;
542
543 case MRI_float: /* fim datum is floats */
544 for( ii=0 ; ii < nxyz ; ii++ )
545 if( bthr[ii] < thrplu && (thrmin == 0 || bthr[ii] > thrmin) ){
546 ffim[ii] = 0.0;
547 }
548 break ;
549
550 case MRI_complex: /* fim datum is complex */
551 for( ii=0 ; ii < nxyz ; ii++ )
552 if( bthr[ii] < thrplu && (thrmin == 0 || bthr[ii] > thrmin) ){
553 cfim[ii].r = cfim[ii].i = 0.0 ;
554 }
555 break ;
556 }
557 }
558 break ;
559
560 /** threshold datum is floats **/
561
562 case MRI_float:{
563 float thrplu , thrmin ;
564 thrplu = edit_thresh ; if( thrfac > 0.0 ) thrplu /= thrfac ;
565 thrmin = edit_thbot ; if( thrfac > 0.0 ) thrmin /= thrfac ;
566 switch( fim_type ){
567 case MRI_short: /* fim datum is shorts */
568 for( ii=0 ; ii < nxyz ; ii++ )
569 if( fthr[ii] < thrplu && fthr[ii] > thrmin ){ sfim[ii] = 0 ; }
570 break ;
571
572 case MRI_byte: /* fim datum is bytes */
573 for( ii=0 ; ii < nxyz ; ii++ )
574 if( fthr[ii] < thrplu && fthr[ii] > thrmin ){ bfim[ii] = 0 ; }
575 break ;
576
577 case MRI_float: /* fim datum is floats */
578 for( ii=0 ; ii < nxyz ; ii++ )
579 if( fthr[ii] < thrplu && fthr[ii] > thrmin ){ ffim[ii] = 0.0 ; }
580 break ;
581
582 case MRI_complex: /* fim datum is complex */
583 for( ii=0 ; ii < nxyz ; ii++ )
584 if( fthr[ii] < thrplu && fthr[ii] > thrmin ){
585 cfim[ii].r = cfim[ii].i = 0.0 ;
586 }
587 break ;
588 }
589 }
590 break ;
591 }
592 }
593
594 /*----- blur? -----*/
595
596 if( AFNI_yesenv("AFNI_BLUR_FFT") ) EDIT_blur_allow_fir(0) ; /* 04 Oct 2005 */
597
598 blur3d = ( edit_blurx > 0.0f || edit_blury > 0.0f || edit_blurz > 0.0f ) ;
599
600 if( blur3d ){ /* 17 Jun 2019 */
601
602 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: blurring_3D fim\n") ;
603 EDIT_blur_volume_3d( nx,ny,nz, dx,dy,dz , fim_type,vfim ,
604 edit_blurx,edit_blury,edit_blurz ) ;
605
606 } else if( edit_blur > 0.0 ){
607
608 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: blurring fim\n") ;
609 EDIT_blur_volume( nx,ny,nz, dx,dy,dz , fim_type,vfim , edit_blur ) ;
610
611 }
612
613 /*----- threshold blur? -----*/ /* 4 Oct 1996 */
614 if(( edit_thrblur > 0.0) && (vthr != NULL) ){
615
616 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: blurring threshold\n") ;
617
618 EDIT_blur_volume( nx,ny,nz, dx,dy,dz , thr_type,vthr , edit_thrblur ) ;
619 }
620
621
622 /*----- zvol? -----*/
623
624 if( edit_zvol ){
625 THD_ivec3 iv1 , iv2 ;
626 int ix1,ix2 , jy1,jy2 , kz1,kz2 , jj,kk ;
627
628 iv1 = THD_3dmm_to_3dind(dset,TEMP_FVEC3(edopt->zv_x1,edopt->zv_y1,edopt->zv_z1));
629 iv2 = THD_3dmm_to_3dind(dset,TEMP_FVEC3(edopt->zv_x2,edopt->zv_y2,edopt->zv_z2));
630
631 ix1 = iv1.ijk[0] ; ix2 = iv2.ijk[0] ;
632 jy1 = iv1.ijk[1] ; jy2 = iv2.ijk[1] ;
633 kz1 = iv1.ijk[2] ; kz2 = iv2.ijk[2] ;
634
635 if( ix1 > ix2 ){ ii=ix1 ; ix1=ix2 ; ix2=ii ; }
636 if( jy1 > jy2 ){ ii=jy1 ; jy1=jy2 ; jy2=ii ; }
637 if( kz1 > kz2 ){ ii=kz1 ; kz1=kz2 ; kz2=ii ; }
638
639 #ifdef AFNI_DEBUG
640 { char str[256] ;
641 sprintf(str,"edit_zvol: x1=%g x2=%g y1=%g y2=%g z1=%g z2=%g",
642 edopt->zv_x1,edopt->zv_x2,edopt->zv_y1,edopt->zv_y2,edopt->zv_z1,edopt->zv_z2) ;
643 STATUS(str) ;
644 sprintf(str," : ix1=%d ix2=%d jy1=%d jy2=%d kz1=%d kz2=%d",
645 ix1,ix2,jy1,jy2,kz1,kz2) ;
646 STATUS(str) ; }
647 #endif
648
649 if( verbose )
650 fprintf(stderr,"--- EDIT_one_dataset: zeroing indexes [%d,%d]x[%d,%d]x[%d,%d]\n",
651 ix1,ix2,jy1,jy2,kz1,kz2 ) ;
652
653 for( kk=kz1 ; kk <= kz2 ; kk++ ){
654 for( jj=jy1 ; jj <= jy2 ; jj++ ){
655 switch( fim_type ){
656 case MRI_short:
657 for( ii=ix1 ; ii <= ix2 ; ii++ ) sfim[ii+jj*nx+kk*nxy] = 0 ;
658 break ;
659
660 case MRI_byte:
661 for( ii=ix1 ; ii <= ix2 ; ii++ ) bfim[ii+jj*nx+kk*nxy] = 0 ;
662 break ;
663
664 case MRI_float:
665 for( ii=ix1 ; ii <= ix2 ; ii++ ) ffim[ii+jj*nx+kk*nxy] = 0 ;
666 break ;
667
668 case MRI_complex:
669 for( ii=ix1 ; ii <= ix2 ; ii++ )
670 cfim[ii+jj*nx+kk*nxy].r = cfim[ii+jj*nx+kk*nxy].i = 0 ;
671 break ;
672 }
673 }
674 }
675 }
676
677 /*----- form clusters? -----*/
678
679 rmm = clust_rmm ;
680 vmul = clust_vmul ;
681
682 if( rmm >= 0.0 ){ /* do clustering? */
683
684 MCW_cluster_array *clbig ;
685 MCW_cluster *cl ;
686
687 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: clustering with rmm=%g vmul=%g\n",
688 rmm,vmul ) ;
689
690 /*----- Erosion and dilation of clusters -----*/ /* 17 June 1998 */
691 if (erode_pv > 0.0)
692 MCW_erode_clusters (nx, ny, nz, dx, dy, dz, fim_type, vfim, rmm,
693 erode_pv, dilate);
694
695
696 STATUS("clustering") ;
697
698 if( vmul >= 0.0 )
699 ptmin = (int)( vmul / dxyz + 0.99 ) ;
700 else
701 ptmin = (int) fabs(vmul) ; /* 30 Apr 2002 */
702
703 vmul = MAX(1,ptmin) * dxyz ; /* for use below */
704
705 /* ZSS March 03 2010, changed from MCW_find_clusters
706 to NIH_find_clusters to accommodate isomode */
707 clar = NIH_find_clusters( nx,ny,nz , dx,dy,dz ,
708 fim_type,vfim , rmm, edopt->isomode) ;
709 nclu = 0 ;
710
711 if( clar != NULL ){
712 INIT_CLARR(clbig) ;
713 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
714 cl = clar->clar[iclu] ;
715 if( cl->num_pt >= ptmin ){ /* big enough */
716 ADDTO_CLARR(clbig,cl) ; /* copy pointer */
717 clar->clar[iclu] = NULL ; /* null out original */
718 nclu++ ;
719 }
720 }
721 DESTROY_CLARR(clar) ;
722 clar = clbig ;
723 if( nclu == 0 || clar == NULL || clar->num_clu == 0 ){
724 printf("*** NO CLUSTERS FOUND ***\n") ;
725 if( clar != NULL ) DESTROY_CLARR(clar) ;
726 EXRETURN ;
727 }
728 SORT_CLARR(clar) ;
729 }
730
731 if( nclu == 0 ){ /* no data left */
732 STATUS("no data left after cluster edit!") ;
733 DESTROY_CLARR(clar) ;
734 EXRETURN ;
735 }
736
737 #ifdef AFNI_DEBUG
738 { char str[256] ;
739 sprintf(str,"number clusters = %d",nclu) ; STATUS(str) ; }
740 #endif
741
742 /*----- edit clusters? -----*/ /* 10 Sept 1996 */
743 if (edit_clust > ECFLAG_SAME)
744 EDIT_cluster_array (clar, edit_clust, dxyz, vmul);
745 if (edit_clust == ECFLAG_SIZE ||
746 edit_clust == ECFLAG_ORDER ||
747 edit_clust == ECFLAG_DEPTH )
748 DSET_BRICK_FACTOR(dset,iv_fim) = 0.0;
749
750 for( iclu=0 ; iclu < clar->num_clu ; iclu++ )
751 if( clar->clar[iclu] != NULL && clar->clar[iclu]->num_pt > 0 ){
752 MCW_cluster_to_vol( nx,ny,nz , fim_type, vfim , clar->clar[iclu] ) ;
753 } else {
754 }
755
756 /* now deal with ECFLAG_DEPTH */
757
758 if (edit_clust == ECFLAG_DEPTH) {
759 int iwarn=0;
760 byte *mask=NULL;
761 short *depth=NULL;
762 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ) {
763 if ( clar->clar[iclu] && clar->clar[iclu]->num_pt > 0 ){
764 /* allocate if 1st time */
765 if (!mask) mask = (byte *)calloc(nxyz, sizeof(byte));
766 if (!depth) depth = (short*)calloc(nxyz, sizeof(short));
767 if (!mask || !depth) {
768 ERROR_message("NULL mask and/or depth");
769 EXRETURN ;
770 }
771 /* start clean */
772 memset(mask, 0, sizeof(byte)*nxyz);
773 memset(depth, 0, sizeof(short)*nxyz);
774 /* fill mask */
775 #if 0
776 switch( fim_type ){
777 case MRI_short: /* fim datum is shorts */
778
779 for( ii=0 ; ii < nxyz ; ii++ ) {
780 if( sfim[ii] == (short)clar->clar[iclu]->mag[0] )
781 mask[ii] = 1 ;
782 }
783 break ;
784
785 case MRI_byte: /* fim datum is bytes */
786 for( ii=0 ; ii < nxyz ; ii++ ) {
787 if( bfim[ii] == (byte)clar->clar[iclu]->mag[0] )
788 mask[ii] = 1 ;
789 }
790 break ;
791
792 case MRI_float: /* fim datum is floats */
793 for( ii=0 ; ii < nxyz ; ii++ ) {
794 if( ffim[ii] == (float)clar->clar[iclu]->mag[0] )
795 mask[ii] = 1 ;
796 }
797 break ;
798
799 case MRI_complex: /* fim datum is complex */
800 for( ii=0 ; ii < nxyz ; ii++ ) {
801 if( cfim[ii].r == (float)clar->clar[iclu]->mag[0] )
802 mask[ii] = 1 ;
803 }
804 break ;
805
806 default:
807 ERROR_message("Bad type");
808 break;
809
810 }/* end of switch */
811 #endif
812 for( icl=0 ; icl < clar->clar[iclu]->num_pt ; icl++ ) {
813 ijk = THREE_TO_IJK (clar->clar[iclu]->i[icl],
814 clar->clar[iclu]->j[icl],
815 clar->clar[iclu]->k[icl],
816 nx, nxy);
817 mask[ijk] = 1 ;
818 }
819 if (!THD_mask_depth ( nx, ny, nz, mask, 1, depth,2)) {
820 ERROR_message("Failed in mask depth");
821 break;
822 }
823 /* Now put depth values back in fim */
824 switch( fim_type ){
825 case MRI_short:
826 for( ii=0 ; ii < nxyz ; ii++ ) {
827 if (mask[ii]) sfim[ii] = depth[ii];
828 }
829 break;
830 case MRI_byte:
831 for( ii=0 ; ii < nxyz ; ii++ ) {
832 if (mask[ii]) {
833 if (depth[ii] > MRI_maxbyte) {
834 if (!iwarn) {
835 WARNING_message("Depth values exceed byte range,"
836 " clipping at %d\n", MRI_maxbyte);
837 ++iwarn;
838 }
839 depth[ii] = MRI_maxbyte;
840 }
841 bfim[ii] = (byte)depth[ii];
842 }
843 }
844 break;
845 case MRI_float:
846 for( ii=0 ; ii < nxyz ; ii++ ) {
847 if (mask[ii]) ffim[ii] = (float)depth[ii];
848 }
849 break;
850 case MRI_complex:
851 for( ii=0 ; ii < nxyz ; ii++ ) {
852 if (mask[ii]) cfim[ii].r = cfim[ii].i =(float)depth[ii];
853 }
854 break;
855 } /* end of second switch */
856 } /* end of is non empty cluster */
857 }/* end of loop over clusters */
858 if (mask) free(mask); mask=NULL;
859 if (depth) free(depth); depth=NULL;
860 } /* end of ORDER trip */
861
862 DESTROY_CLARR(clar) ;
863 }
864
865
866 /*----- filter? -----*/ /* 11 Sept 1996 */
867 if (filter_opt > FCFLAG_NONE){
868
869 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: filtering fim\n") ;
870
871 EDIT_filter_volume (nx, ny, nz, dx, dy, dz, fim_type, vfim,
872 filter_opt, filter_rmm ,
873 edopt->fmask,edopt->fmclip , edopt->fexpr );
874 }
875
876
877 /*----- threshold filter? -----*/ /* 1 Oct 1996 */
878 if ((thrfilter_opt > FCFLAG_NONE) && (vthr != NULL)){
879
880 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: filtering thr\n") ;
881
882 EDIT_filter_volume (nx, ny, nz, dx, dy, dz, thr_type, vthr,
883 thrfilter_opt, thrfilter_rmm ,
884 edopt->fmask,edopt->fmclip , edopt->fexpr );
885 }
886
887
888 /*----- scale? -----*/
889
890 #ifdef ALLOW_SCALE_TO_MAX
891 if( edit_scale ){
892 STATUS("scale") ;
893
894 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: scaling fim to max\n") ;
895
896 MCW_scale_to_max( nx,ny,nz , fim_type , vfim ) ;
897 }
898 #endif
899
900 /*----- mult? -----*/
901 /*--- correction for scaling of short and byte bricks (13 Sept. 1996) ---*/
902
903 if( edit_mult != 0.0 ){
904 STATUS("mult") ;
905
906 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: multiplying fim\n") ;
907
908 switch( fim_type ){
909 case MRI_short:
910 if (fimfac > 0)
911 DSET_BRICK_FACTOR(dset,iv_fim) =
912 DSET_BRICK_FACTOR(dset,iv_fim) * edit_mult ;
913 else
914 for( ii=0 ; ii < nxyz ; ii++ ) sfim[ii] *= edit_mult ;
915 break ;
916
917 case MRI_byte :
918 if (fimfac > 0)
919 DSET_BRICK_FACTOR(dset,iv_fim) =
920 DSET_BRICK_FACTOR(dset,iv_fim) * edit_mult ;
921 else
922 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] *= edit_mult ;
923 break ;
924
925 case MRI_float: for( ii=0 ; ii < nxyz ; ii++ ) ffim[ii] *= edit_mult ;
926 break ;
927
928 case MRI_complex: for( ii=0 ; ii < nxyz ; ii++ )
929 cfim[ii].r *= edit_mult , cfim[ii].i *= edit_mult ;
930 break ;
931 }
932 }
933
934 /*----- 17 Sep 1998: conversion to z-score? -----*/
935
936 if( edopt->zscore ){ /* 07 Jun 1999! How did this get lost? */
937 int kv = DSET_BRICK_STATCODE(dset,iv_fim) ;
938 float par[2] ;
939
940 if( FUNC_IS_STAT(kv) && kv != FUNC_ZT_TYPE ){
941
942 #if 0
943 fprintf(stderr," -1zscore: converting\n") ;
944 #endif
945
946 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: converting to zscore\n") ;
947
948 EDIT_zscore_vol( nxyz , fim_type , fimfac , vfim ,
949 kv , DSET_BRICK_STATAUX(dset,iv_fim) ) ;
950
951 if( ISBUCKET(dset) || dset->dblk->brick_statcode != NULL ){
952
953 #if 0
954 fprintf(stderr," -1zscore: bucketing\n") ;
955 #endif
956
957 par[0] = FUNC_ZT_TYPE ;
958 par[1] = 0 ;
959 EDIT_dset_items( dset , ADN_brick_stataux_one+iv_fim,par , ADN_none ) ;
960
961 } else if( ISFUNC(dset) &&
962 FUNC_IS_STAT(dset->func_type) &&
963 iv_fim == FUNC_ival_thr[dset->func_type] ){
964
965 #if 0
966 fprintf(stderr," -1zscore: retyping\n") ;
967 #endif
968
969 dset->func_type = FUNC_ZT_TYPE ;
970 dset->stat_aux[0] = 0.0 ;
971
972 } else {
973 fprintf(stderr,"*** -1zscore error: non-bucket & non-func!\n") ;
974 }
975
976 if( fim_type == MRI_short )
977 DSET_BRICK_FACTOR(dset,iv_fim) = 1.0 / FUNC_ZT_SCALE_SHORT ;
978 }
979 }
980
981 /* turn dset to rank */
982 if ( edopt->rank) {
983 if( verbose )
984 fprintf(stderr,"--- EDIT_one_dataset: Converting to rank value.\n");
985 if (!(THD_unique_rank_edit(dset ,
986 iv_fim,
987 edopt->fmask,
988 edopt->rankmapname, NULL, NULL))) {
989 fprintf(stderr,"*** Ranking error.\n");
990 }
991 }
992 /*------ DONE! -----*/
993
994 EXRETURN ;
995 }
996