1 #include "mrilib.h"
2
3 #define NTYPE_SPHERE 1
4 #define NTYPE_RECT 2
5 #define NTYPE_NULL 666
6
7 #define BAILOUT 3
8 #define BLURMAX 0.05f
9 #define DFGOAL 0.2f
10 #define DFGOALC (1.0f-DFGOAL)
11
12 static int verb = 1 ;
13
14 void estimate_blur_map( MRI_IMARR *bmar , byte *mask , MCW_cluster *nbhd ,
15 float *fxar , float *fyar , float *fzar ) ;
16
17 #ifdef USE_OMP
18 # include <omp.h>
19 # include "mri_fwhm.c"
20 # include "mri_blur3d_variable.c"
21 #endif
22
23 /*----------------------------------------------------------------------------*/
24
main(int argc,char * argv[])25 int main( int argc , char *argv[] )
26 {
27 int iarg=1 , ii,nvox ;
28 THD_3dim_dataset *bmset=NULL, *inset=NULL, *outset=NULL , *mset=NULL ;
29 char *prefix="./blurto" ;
30 float fwhm_goal=0.0f , fwhm_subgoal ; int fwhm_2D=0 ;
31 MCW_cluster *nbhd=NULL ;
32 byte *mask=NULL ; int mask_nx=0,mask_ny=0,mask_nz=0 , automask=0 , nmask ;
33 int ntype=0 ; float na=0.0f,nb=0.0f,nc=0.0f ;
34 MRI_IMARR *bmar ; MRI_IMAGE *bmed , *bmim ; int ibm ;
35 MRI_IMARR *dsar ; MRI_IMAGE *dsim ; int ids ;
36 MRI_IMAGE *fxim=NULL , *fyim=NULL , *fzim=NULL ;
37 float *fxar=NULL , *fyar=NULL , *fzar=NULL ;
38 float dx,dy,dz=0.0f , hx,hy,hz=0.0f , qx,qy,qz=0.0f ;
39 float gx,gy,gz=0.0f , val , maxfxyz , maxfx,maxfy,maxfz ;
40 int nite , bmeqin=0 , maxite=0 , numfxyz , nd,nblur , xdone,ydone,zdone ;
41 int xstall , ystall , zstall ;
42 float bx,by,bz ;
43 float last_fwx=0.0f , last_fwy=0.0f , last_fwz=0.0f ;
44 float delt_fwx , delt_fwy , delt_fwz ;
45 char *bsave_prefix=NULL ;
46 THD_fvec3 fw ; int do_acf=0 ;
47 int nbail=0 , xalmost,yalmost,zalmost , xstopped=0,ystopped=0,zstopped=0 ;
48 float blurfac , blurmax=BLURMAX ;
49 float xrat,yrat,zrat , dmin ;
50 int temperize=0 , temper_fx=0, temper_fy=0, temper_fz=0 ;
51 float fx_tbot=0.0f,fx_ttop=0.0f,fy_tbot=0.0f,fy_ttop=0.0f,fz_tbot=0.0f,fz_ttop=0.0f ;
52 float fx_trat=0.0f,fy_trat=0.0f,fz_trat=0.0f ;
53 int bmall=0,do_unif=0 ; /* 11 Dec 2006 */
54 MRI_IMARR *imar ; MRI_IMAGE *imed, *imad ; float *imedar=NULL, *imadar=NULL, *bar ;
55
56 int corder_bm=-1 , corder_in=0 ; /* 04 Jun 2007: detrending */
57 MRI_IMARR *corder_inar=NULL ; MRI_IMAGE *corder_invv ;
58 float **corder_inref=NULL ;
59 int corder_inrefnum=0 ;
60
61 /*-- help the pitifully ignorant luser? --*/
62
63 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
64 printf(
65 "Usage: 3dBlurToFWHM [options]\n"
66 "\n"
67 "Blurs a 'master' dataset until it reaches a specified FWHM\n"
68 "smoothness (approximately). The same blurring schedule is\n"
69 "applied to the input dataset to produce the output. The goal\n"
70 "is to make the output dataset have the given smoothness, no\n"
71 "matter what smoothness it had on input (however, the program\n"
72 "cannot 'unsmooth' a dataset!). See below for the METHOD used.\n"
73 "\n"
74 "OPTIONS\n"
75 "-------\n"
76 " -input ddd = This required 'option' specifies the dataset\n"
77 " that will be smoothed and output.\n"
78 " -blurmaster bbb = This option specifies the dataset whose\n"
79 " whose smoothness controls the process.\n"
80 " **N.B.: If not given, the input dataset is used.\n"
81 " **N.B.: This should be one continuous run.\n"
82 " Do not input catenated runs!\n"
83 " -prefix ppp = Prefix for output dataset will be 'ppp'.\n"
84 " **N.B.: Output dataset is always in float format.\n"
85 " -mask mmm = Mask dataset, if desired. Blurring will\n"
86 " occur only within the mask. Voxels NOT in\n"
87 " the mask will be set to zero in the output.\n"
88 " -automask = Create an automask from the input dataset.\n"
89 " **N.B.: Not useful if the input dataset has been\n"
90 " detrended or otherwise regressed before input!\n"
91 " -FWHM f = Blur until the 3D FWHM is 'f'.\n"
92 " -FWHMxy f = Blur until the 2D (x,y)-plane FWHM is 'f'.\n"
93 " No blurring is done along the z-axis.\n"
94 " **N.B.: Note that you can't REDUCE the smoothness\n"
95 " of a dataset.\n"
96 " **N.B.: Here, 'x', 'y', and 'z' refer to the\n"
97 " grid/slice order as stored in the dataset,\n"
98 " not DICOM ordered coordinates!\n"
99 " **N.B.: With -FWHMxy, smoothing is done only in the\n"
100 " dataset xy-plane. With -FWHM, smoothing\n"
101 " is done in 3D.\n"
102 " **N.B.: The actual goal is reached when\n"
103 " -FWHM : cbrt(FWHMx*FWHMy*FWHMz) >= f\n"
104 " -FWHMxy: sqrt(FWHMx*FWHMy) >= f\n"
105 " That is, when the area or volume of a\n"
106 " 'resolution element' goes past a threshold.\n"
107 " -quiet Shut up the verbose progress reports.\n"
108 " **N.B.: This should be the first option, to stifle\n"
109 " any verbosity from the option processing code.\n"
110 "\n"
111 "FILE RECOMMENDATIONS for -blurmaster:\n"
112 "For FMRI statistical purposes, you DO NOT want the FWHM to reflect\n"
113 " the spatial structure of the underlying anatomy. Rather, you want\n"
114 " the FWHM to reflect the spatial structure of the noise. This means\n"
115 " that the -blurmaster dataset should not have anatomical structure. One\n"
116 " good form of input is the output of '3dDeconvolve -errts', which is\n"
117 " the residuals left over after the GLM fitted signal model is subtracted\n"
118 " out from each voxel's time series. You can also use the output of\n"
119 " '3dREMLfit -Rerrts' or '3dREMLfit -Rwherr' for this purpose.\n"
120 "You CAN give a multi-brick EPI dataset as the -blurmaster dataset; the\n"
121 " dataset will be detrended in time (like the -detrend option in 3dFWHMx)\n"
122 " which will tend to remove the spatial structure. This makes it\n"
123 " practicable to make the input and blurmaster datasets be the same,\n"
124 " without having to create a detrended or residual dataset beforehand.\n"
125 " Considering the accuracy of blurring estimates, this is probably good\n"
126 " enough for government work [that is an insider's joke :-]. \n"
127 " N.B.: Do not use catenated runs as blurmasters. There should\n"
128 " be no discontinuities in the time axis of blurmaster, which would\n"
129 " make the simple regression detrending do peculiar things.\n"
130 "\n"
131 "ALSO SEE:\n"
132 " * 3dFWHMx, which estimates smoothness globally\n"
133 " * 3dLocalstat -stat FWHM, which estimates smoothness locally\n"
134 " * This paper, which discusses the need for a fixed level of smoothness\n"
135 " when combining FMRI datasets from different scanner platforms:\n"
136 " Friedman L, Glover GH, Krenz D, Magnotta V; The FIRST BIRN. \n"
137 " Reducing inter-scanner variability of activation in a multicenter\n"
138 " fMRI study: role of smoothness equalization.\n"
139 " Neuroimage. 2006 Oct 1;32(4):1656-68.\n"
140 "\n"
141 "METHOD:\n"
142 "The blurring is done by a conservative finite difference approximation\n"
143 "to the diffusion equation:\n"
144 " du/dt = d/dx[ D_x(x,y,z) du/dx ] + d/dy[ D_y(x,y,z) du/dy ]\n"
145 " + d/dz[ D_z(x,y,z) du/dz ]\n"
146 " = div[ D(x,y,z) grad[u(x,y,z)] ]\n"
147 "where diffusion tensor D() is diagonal, Euler time-stepping is used, and\n"
148 "with Neumann (reflecting) boundary conditions at the edges of the mask\n"
149 "(which ensures that voxel data inside and outside the mask don't mix).\n"
150 "* At each pseudo-time step, the FWHM is estimated globally (like '3dFWHMx')\n"
151 " and locally (like '3dLocalstat -stat FWHM'). Voxels where the local FWHM\n"
152 " goes past the goal will not be smoothed any more (D gets set to zero).\n"
153 "* When the global smoothness estimate gets close to the goal, the blurring\n"
154 " rate (pseudo-time step) will be reduced, to avoid over-smoothing.\n"
155 "* When an individual direction's smoothness (e.g., FWHMz) goes past the goal,\n"
156 " all smoothing in that direction stops, but the other directions continue\n"
157 " to be smoothed until the overall resolution element goal is achieved.\n"
158 "* When the global FWHM estimate reaches the goal, the program is done.\n"
159 " It will also stop if progress stalls for some reason, or if the maximum\n"
160 " iteration count is reached (infinite loops being unpopular).\n"
161 "* The output dataset will NOT have exactly the smoothness you ask for, but\n"
162 " it will be close (fondly we do hope). In our Imperial experiments, the\n"
163 " results (measured via 3dFWHMx) are within 10%% of the goal (usually better).\n"
164 "* 2D blurring via -FWHMxy may increase the smoothness in the z-direction\n"
165 " reported by 3dFWHMx, even though there is no inter-slice processing.\n"
166 " At this moment, I'm not sure why. It may be an estimation artifact due\n"
167 " to increased correlation in the xy-plane that biases the variance estimates\n"
168 " used to calculate FWHMz.\n"
169 "\n"
170 "ADVANCED OPTIONS:\n"
171 " -maxite ccc = Set maximum number of iterations to 'ccc' [Default=variable].\n"
172 " -rate rrr = The value of 'rrr' should be a number between\n"
173 " 0.05 and 3.5, inclusive. It is a factor to change\n"
174 " the overall blurring rate (slower for rrr < 1) and thus\n"
175 " require more or less blurring steps. This option should only\n"
176 " be needed to slow down the program if the it over-smooths\n"
177 " significantly (e.g., it overshoots the desired FWHM in\n"
178 " Iteration #1 or #2). You can increase the speed by using\n"
179 " rrr > 1, but be careful and examine the output.\n"
180 " -nbhd nnn = As in 3dLocalstat, specifies the neighborhood\n"
181 " used to compute local smoothness.\n"
182 " [Default = 'SPHERE(-4)' in 3D, 'SPHERE(-6)' in 2D]\n"
183 " ** N.B.: For the 2D -FWHMxy, a 'SPHERE()' nbhd\n"
184 " is really a circle in the xy-plane.\n"
185 " ** N.B.: If you do NOT want to estimate local\n"
186 " smoothness, use '-nbhd NULL'.\n"
187 " -ACF or -acf = Use the 'ACF' method (from 3dFWHMx) to estimate\n"
188 " the global smoothness, rather than the 'classic'\n"
189 " Forman 1995 method. This option will be somewhat\n"
190 " slower. It will also set '-nbhd NULL', since there\n"
191 " is no local ACF estimation method implemented.\n"
192 " -bsave bbb = Save the local smoothness estimates at each iteration\n"
193 " with dataset prefix 'bbb' [for debugging purposes].\n"
194 " -bmall = Use all blurmaster sub-bricks.\n"
195 " [Default: a subset will be chosen, for speed]\n"
196 " -unif = Uniformize the voxel-wise MAD in the blurmaster AND\n"
197 " input datasets prior to blurring. Will be restored\n"
198 " in the output dataset.\n"
199 " -detrend = Detrend blurmaster dataset to order NT/30 before starting.\n"
200 " -nodetrend = Turn off detrending of blurmaster.\n"
201 " ** N.B.: '-detrend' is the new default [05 Jun 2007]!\n"
202 " -detin = Also detrend input before blurring it, then retrend\n"
203 " it afterwards. [Off by default]\n"
204 " -temper = Try harder to make the smoothness spatially uniform.\n"
205 "\n"
206 "-- Author: The Dreaded Emperor Zhark - Nov 2006\n"
207 ) ;
208 PRINT_COMPILE_DATE ; exit(0) ;
209 }
210
211 /*---- official startup ---*/
212
213 PRINT_VERSION("3dBlurToFWHM"); mainENTRY("3dBlurToFWHM main"); machdep();
214 AFNI_logger("3dBlurToFWHM",argc,argv);
215
216 /*---- loop over options ----*/
217
218 while( iarg < argc && argv[iarg][0] == '-' ){
219
220 if( strcmp(argv[iarg],"-nodetrend") == 0 ){ /* 05 Jun 2007 */
221 corder_bm = corder_in = 0 ;
222 iarg++ ; continue ;
223 }
224
225 if( strcmp(argv[iarg],"-detin") == 0 ){ /* 05 Jun 2007 */
226 corder_bm = corder_in = -1 ; do_unif = 0 ;
227 iarg++ ; continue ;
228 }
229
230 if( strcmp(argv[iarg],"-detrend") == 0 ){ /* 04 Jun 2007 */
231 corder_bm = -1 ; do_unif = 0 ;
232 iarg++ ; continue ;
233 }
234
235 if( strcmp(argv[iarg],"-rate") == 0 ){
236 if( ++iarg >= argc ) ERROR_exit("Need argument after '-rate'") ;
237 val = (float)strtod(argv[iarg],NULL) ;
238 if( val >= 0.05f && val <= 3.555f ) blurmax = BLURMAX * val ;
239 else ERROR_exit("Illegal value after '-rate': '%s'",argv[iarg]) ;
240 iarg++ ; continue ;
241 }
242
243 if( strncmp(argv[iarg],"-q",2) == 0 ){
244 verb = 0 ; iarg++ ; continue ;
245 }
246 if( strncmp(argv[iarg],"-verb",4) == 0 ){
247 verb = 1 ; iarg++ ; continue ;
248 }
249
250 if( strncmp(argv[iarg],"-bmal",5) == 0 ){
251 bmall = 1 ; iarg++ ; continue ;
252 }
253
254 if( strncmp(argv[iarg],"-unif",4) == 0 ){
255 do_unif = 1 ; iarg++ ; continue ;
256 }
257 if( strncmp(argv[iarg],"-nounif",6) == 0 ){
258 do_unif = 0 ; iarg++ ; continue ;
259 }
260
261 if( strncmp(argv[iarg],"-temper",6) == 0 ){
262 temperize = 1 ; iarg++ ; continue ;
263 }
264
265 if( strcmp(argv[iarg],"-maxite") == 0 ){
266 if( ++iarg >= argc ) ERROR_exit("Need argument after '-maxite'") ;
267 maxite = (int)strtod(argv[iarg],NULL) ;
268 if( maxite <= 0 ) maxite = 666 ;
269 iarg++ ; continue ;
270 }
271
272 if( strcmp(argv[iarg],"-input") == 0 || strcmp(argv[iarg],"-dset") == 0 ){
273 if( inset != NULL ) ERROR_exit("Can't have two -input options") ;
274 if( ++iarg >= argc ) ERROR_exit("Need argument after '-input'") ;
275 inset = THD_open_dataset( argv[iarg] );
276 CHECK_OPEN_ERROR(inset,argv[iarg]) ;
277 iarg++ ; continue ;
278 }
279
280 if( strncmp(argv[iarg],"-blurm",6) == 0 ){
281 if( bmset != NULL ) ERROR_exit("Can't have two -blurmaster options") ;
282 if( ++iarg >= argc ) ERROR_exit("Need argument after '-blurmaster'") ;
283 bmset = THD_open_dataset( argv[iarg] ) ;
284 CHECK_OPEN_ERROR(bmset,argv[iarg]) ;
285 iarg++ ; continue ;
286 }
287
288 if( strcmp(argv[iarg],"-prefix") == 0 ){
289 if( ++iarg >= argc ) ERROR_exit("Need argument after '-prefix'") ;
290 prefix = strdup(argv[iarg]) ;
291 if( !THD_filename_ok(prefix) ) ERROR_exit("Bad name after '-prefix'") ;
292 iarg++ ; continue ;
293 }
294
295 if( strcmp(argv[iarg],"-bsave") == 0 ){
296 if( ++iarg >= argc ) ERROR_exit("Need argument after '-bsave'") ;
297 bsave_prefix = strdup(argv[iarg]) ;
298 if( !THD_filename_ok(bsave_prefix) ) ERROR_exit("Bad name after '-bsave'") ;
299 iarg++ ; continue ;
300 }
301
302 if( strcmp(argv[iarg],"-mask") == 0 ){
303 if( ++iarg >= argc ) ERROR_exit("Need argument after '-mask'") ;
304 if( mask != NULL || automask ) ERROR_exit("Can't have two mask inputs") ;
305 mset = THD_open_dataset( argv[iarg] ) ;
306 CHECK_OPEN_ERROR(mset,argv[iarg]) ;
307 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
308 mask_nx = DSET_NX(mset); mask_ny = DSET_NY(mset); mask_nz = DSET_NZ(mset);
309 mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_unload(mset) ;
310 if( mask == NULL ) ERROR_exit("Can't make mask from dataset '%s'",argv[iarg]) ;
311 nmask = THD_countmask( mask_nx*mask_ny*mask_nz , mask ) ;
312 if( verb ) INFO_message("Number of voxels in mask = %d",nmask) ;
313 if( nmask < 333 ) ERROR_exit("Mask is too small to process") ;
314 iarg++ ; continue ;
315 }
316
317 if( strcmp(argv[iarg],"-automask") == 0 ){
318 if( mask != NULL ) ERROR_exit("Can't have -automask and -mask") ;
319 automask = 1 ; iarg++ ; continue ;
320 }
321
322 if( strcasecmp(argv[iarg],"-acf") == 0 ){ /* 30 Dec 2015 */
323 do_acf++ ; iarg++ ; continue ;
324 }
325
326 if( strcmp(argv[iarg],"-nbhd") == 0 ){
327 char *cpt ;
328 if( ntype > 0 ) ERROR_exit("Can't have 2 '-nbhd' options") ;
329 if( ++iarg >= argc ) ERROR_exit("Need argument after '-nbhd'") ;
330
331 cpt = argv[iarg] ;
332 if( strncasecmp(cpt,"SPHERE",6) == 0 ){
333 sscanf( cpt+7 , "%f" , &na ) ;
334 if( na == 0.0f ) ERROR_exit("Can't have a SPHERE of radius 0") ;
335 ntype = NTYPE_SPHERE ;
336 } else if( strncasecmp(cpt,"RECT",4) == 0 ){
337 sscanf( cpt+5 , "%f,%f,%f" , &na,&nb,&nc ) ;
338 if( na == 0.0f && nb == 0.0f && nc == 0.0f )
339 ERROR_exit("'RECT(0,0,0)' is not a legal neighborhood") ;
340 ntype = NTYPE_RECT ;
341 } else if( strcasecmp(cpt,"NULL") == 0 || *cpt == '0' ){
342 ntype = NTYPE_NULL ;
343 } else {
344 ERROR_exit("Unknown -nbhd shape: '%s'",cpt) ;
345 }
346 iarg++ ; continue ;
347 }
348
349 if( strcmp(argv[iarg],"-FWHM") == 0 || strcmp(argv[iarg],"-FHWM") == 0 ){
350 if( ++iarg >= argc ) ERROR_exit("Need argument after '%s'",argv[iarg-1]);
351 val = (float)strtod(argv[iarg],NULL) ;
352 if( val <= 0.0f ) ERROR_exit("Illegal value after '%s': '%s'",
353 argv[iarg-1],argv[iarg]) ;
354 fwhm_goal = val ; fwhm_2D = 0 ;
355 iarg++ ; continue ;
356 }
357
358 if( strcmp(argv[iarg],"-FWHMxy") == 0 || strcmp(argv[iarg],"-FHWMxy") == 0 ){
359 if( ++iarg >= argc ) ERROR_exit("Need argument after '%s'",argv[iarg-1]);
360 val = (float)strtod(argv[iarg],NULL) ;
361 if( val <= 0.0f ) ERROR_exit("Illegal value after '%s': '%s'",
362 argv[iarg-1],argv[iarg]) ;
363 fwhm_goal = val ; fwhm_2D = 1 ;
364 iarg++ ; continue ;
365 }
366
367 ERROR_exit("Unknown option '%s'",argv[iarg]) ;
368
369 } /*--- end of loop over options ---*/
370
371 /*----- check for stupid inputs, load datasets, et cetera -----*/
372
373 if( do_acf ) ntype = NTYPE_NULL ; /* 30 Dec 2015 */
374
375 if( fwhm_goal == 0.0f )
376 ERROR_exit("No -FWHM option given! What do you want?") ;
377
378 fwhm_subgoal = 0.95f*fwhm_goal ; /* stop one axis if it gets 95% of the way */
379 /* while the other axes need to catch up */
380 blurfac = blurmax ;
381
382 if( inset == NULL ){
383 if( iarg >= argc ) ERROR_exit("No input dataset on command line?") ;
384 inset = THD_open_dataset( argv[iarg] ) ;
385 CHECK_OPEN_ERROR(inset,argv[iarg]) ;
386 }
387 nvox = DSET_NVOX(inset) ;
388 dx = fabs(DSET_DX(inset)) ; if( dx == 0.0f ) dx = 1.0f ;
389 dy = fabs(DSET_DY(inset)) ; if( dy == 0.0f ) dy = 1.0f ;
390 dz = fabs(DSET_DZ(inset)) ; if( dz == 0.0f ) dz = 1.0f ;
391 dmin = MIN(dx,dy) ; if( !fwhm_2D ) dmin = MIN(dmin,dz) ;
392 xrat = dmin/dx ; xrat = xrat*xrat ;
393 yrat = dmin/dy ; yrat = yrat*yrat ;
394 zrat = dmin/dz ; zrat = zrat*zrat ;
395
396 if( bmset == NULL ){
397 bmset = inset ; bmeqin = 1 ;
398 if( verb ) INFO_message("Using input dataset as blurmaster") ;
399 }
400
401 if( DSET_NX(inset) != DSET_NX(bmset) ||
402 DSET_NY(inset) != DSET_NY(bmset) ||
403 DSET_NZ(inset) != DSET_NZ(bmset) )
404 ERROR_exit("-blurmaster dataset grid doesn't match input dataset") ;
405
406 if( DSET_NZ(inset) == 1 && !fwhm_2D ){
407 WARNING_message("Dataset is 2D ==> switching from -FWHM to -FWHMxy") ;
408 fwhm_2D = 1 ;
409 }
410
411 if( maxite <= 0 ){
412 if( fwhm_2D ) maxite = (int)( 66.6 * fwhm_goal / sqrt(dx*dy) ) ;
413 else maxite = (int)( 66.6 * fwhm_goal / cbrt(dx*dy*dz) ) ;
414 if( maxite < 66 ) maxite = 66 ;
415 INFO_message("Max number iterations set to %d",maxite) ;
416 }
417
418 /*--- deal with mask or automask ---*/
419
420 if( mask != NULL ){
421 if( mask_nx != DSET_NX(inset) ||
422 mask_ny != DSET_NY(inset) ||
423 mask_nz != DSET_NZ(inset) )
424
425 ERROR_exit("-mask dataset grid doesn't match input dataset") ;
426
427 } else if( automask ){
428 mask = THD_automask( inset ) ;
429 if( mask == NULL )
430 ERROR_message("Can't create -automask from input dataset?") ;
431 nmask = THD_countmask( DSET_NVOX(inset) , mask ) ;
432 if( verb ) INFO_message("Number of voxels in automask = %d (out of %d total)",
433 nmask,DSET_NVOX(inset));
434 if( nmask < 333 ) ERROR_exit("Automask is too small to process") ;
435
436 } else {
437 mask = (byte *)malloc(sizeof(byte)*nvox) ; nmask = nvox ;
438 memset(mask,1,sizeof(byte)*nvox) ;
439 if( verb ) INFO_message("No mask ==> processing all %d voxels",nvox);
440 }
441
442 /*---- create neighborhood -----*/
443
444 if( ntype == 0 ){
445 ntype = NTYPE_SPHERE ; na = (fwhm_2D) ? -6.0f : -4.0f ;
446 if( verb ) INFO_message("Using default neighborhood 'SPHERE(%.1f)'",na);
447 }
448
449 switch( ntype ){
450 default:
451 ERROR_exit("WTF? ntype=%d",ntype) ; /* should not happen */
452
453 case NTYPE_NULL: nbhd = NULL ; break ;
454
455 case NTYPE_SPHERE:{
456 float ddx , ddy , ddz ;
457 if( na < 0.0f ){ ddx = ddy = ddz = 1.0f ; na = -na ; }
458 else { ddx = dx ; ddy = dy ; ddz = dz ; }
459 if( fwhm_2D ) ddz = 1.e+10 ; /* 2D only */
460 nbhd = MCW_spheremask( ddx,ddy,ddz , na ) ;
461 }
462 break ;
463
464 case NTYPE_RECT:{
465 float ddx , ddy , ddz ;
466 if( na < 0.0f ){ ddx = 1.0f; na = -na; } else ddx = dx ;
467 if( nb < 0.0f ){ ddy = 1.0f; nb = -nb; } else ddy = dy ;
468 if( nc < 0.0f ){ ddz = 1.0f; nc = -nc; } else ddz = dz ;
469 nbhd = MCW_rectmask( ddx,ddy,ddz , na,nb,nc ) ;
470 }
471 break ;
472 }
473
474 if( nbhd != NULL ){
475 if( verb || nbhd->num_pt < 19)
476 INFO_message("Neighborhood comprises %d voxels",nbhd->num_pt) ;
477 if( nbhd->num_pt < 19 )
478 ERROR_exit("FWHM estimation requires neighborhood of at least 19 voxels!") ;
479 } else {
480 if( bsave_prefix != NULL ){
481 WARNING_message("Option '-bsave' means nothing with '-nbhd NULL'") ;
482 bsave_prefix = NULL ;
483 }
484 if( temperize ){
485 WARNING_message("Option '-temper' means nothing with '-nbhd NULL'") ;
486 temperize = 0 ;
487 }
488 }
489
490 /*--- process blurmaster dataset to produce blur master image array ---*/
491
492 DSET_load(bmset) ; CHECK_LOAD_ERROR(bmset) ;
493 INIT_IMARR(bmar) ;
494
495 if( corder_bm ) do_unif = 0 ;
496 if( corder_bm < 0 ) corder_bm = DSET_NVALS(bmset) / 30 ;
497 #if 0
498 if( corder_bm > 0 && 2*corder_bm+3 >= DSET_NVALS(bmset) )
499 ERROR_exit("-corder %d is too big for blurmaster dataset",corder_bm) ;
500 #endif
501
502 if( DSET_NVALS(bmset) == 1 ){
503
504 bmed = THD_median_brick(bmset) ; /* scaled to floats, if need be */
505 if( mri_allzero(bmed) )
506 ERROR_exit("blurmaster dataset has 1 sub-brick, which is all zero!") ;
507 ADDTO_IMARR(bmar,bmed) ;
508 if( do_unif )
509 WARNING_message("Can't apply -unif: only 1 sub-brick in blurmaster dataset") ;
510 else if( corder_bm > 0 ){
511 WARNING_message("Can't apply -detrend: only 1 sub-brick in blurmaster dataset") ;
512 corder_bm = 0 ;
513 }
514
515 } else {
516
517 MRI_IMAGE *smed ;
518 float *mar , *sar ;
519 int ntouse , nvb=DSET_NVALS(bmset) , ibot , idel , nzs ;
520 int *az=(int *)malloc(sizeof(int)*nvb) , naz=0 , jbm , jj ;
521
522 for( jj=0 ; jj < nvb ; jj++ ){
523 az[jj] = mri_allzero(DSET_BRICK(bmset,jj)); if( az[jj] ) naz++;
524 }
525 if( corder_bm > 0 && naz > 0 ){ /* 04 Nov 2008 */
526 corder_bm = 0 ;
527 WARNING_message("blurmaster detrend disabled: %d all zero sub-bricks found",naz) ;
528 }
529
530 if( corder_bm > 0 ){ /*--- 04 Jun 2007 ---*/
531 int nref=2*corder_bm+3, iv,kk ;
532 float **ref , tm,fac,fq ;
533 THD_3dim_dataset *newset ;
534 INFO_message("detrending blurmaster: %d ref funcs, %d time points",nref,nvb) ;
535
536 ref = THD_build_trigref( corder_bm , nvb ) ;
537 if( ref == NULL ) ERROR_exit("THD_build_trigref failed!") ;
538
539 newset = THD_detrend_dataset( bmset, nref,ref , 2,1 , mask , NULL ) ;
540 if( newset == NULL ) ERROR_exit("detrending failed!?") ;
541
542 if( !bmeqin ) DSET_delete(bmset) ;
543 bmset = newset ;
544 ININFO_message("detrending of blurmaster complete") ;
545
546 for( jj=0 ; jj < nref ; jj++ ) free(ref[jj]) ;
547 free(ref) ;
548 }
549
550 imar = THD_medmad_bricks(bmset) ;
551 bmed = IMARR_SUBIM(imar,0) ; /* 11 Dec 2006: -do_unif: */
552 smed = IMARR_SUBIM(imar,1) ; /* normalize each voxel */
553 mar = MRI_FLOAT_PTR(bmed) ; /* by 1/MAD, to allow for */
554 sar = MRI_FLOAT_PTR(smed) ; /* spatial variability */
555
556 for( nzs=ii=0 ; ii < nvox ; ii++ ) /* 10 May 2007 */
557 if( mask[ii] && sar[ii] == 0.0f ){ mask[ii] = 0; nzs++ ; }
558 if( nzs > 0 )
559 INFO_message("Removed %d voxels with 0 MAD in blurmaster from mask",nzs) ;
560
561 if( do_unif )
562 for( ii=0 ; ii < nvox ; ii++ ) if( sar[ii] != 0.0f ) sar[ii] = 1.0/sar[ii] ;
563
564 if( nbhd != NULL ){
565 ntouse = (int)ceil(9999.9999/nbhd->num_pt) ;
566 ntouse = MAX(2,ntouse) ;
567 ntouse = MIN(nvb,ntouse) ;
568 } else {
569 ntouse = MIN(nvb,32) ;
570 }
571 if( bmall ) ntouse = nvb ;
572 idel = nvb / ntouse ;
573 ibot = (nvb-1 - idel*(ntouse-1)) / 2 ;
574 if( verb ) INFO_message("Using blurmaster sub-bricks [%d..%d(%d)]",
575 ibot,ibot+(ntouse-1)*idel,idel);
576 for( ibm=ibot ; ibm < nvb && IMARR_COUNT(bmar) < ntouse ; ibm+=idel ){
577 jbm = ibm ;
578 if( az[jbm] ){
579 if( idel == 1 ){
580 WARNING_message("skip blurmaster sub-brick #%d because it is all zero",ibm) ;
581 continue ;
582 }
583 for( jbm=ibm-1 ; jbm >= 0 && jbm > ibm-idel && az[jbm] ; jbm-- ) ; /*nada*/
584 if( jbm < 0 || jbm == ibm-idel || az[jbm] ){
585 WARNING_message("skip blurmaster sub-brick #%d because it is all zero",ibm) ;
586 continue ;
587 }
588 if( verb )
589 INFO_message("replace blurmaster all zero sub-brick #%d with #%d",ibm,jbm) ;
590 }
591 bmim = mri_scale_to_float(DSET_BRICK_FACTOR(bmset,jbm),DSET_BRICK(bmset,jbm));
592 bar = MRI_FLOAT_PTR(bmim) ;
593 if( do_unif )
594 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] = (bar[ii]-mar[ii])*sar[ii] ;
595 else if( corder_bm == 0 )
596 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] = bar[ii]-mar[ii] ;
597 ADDTO_IMARR(bmar,bmim) ;
598 if( !bmeqin ) DSET_unload_one(bmset,jbm) ;
599 }
600 DESTROY_IMARR(imar) ; /* median and MAD bricks */
601 free(az) ;
602 if( IMARR_COUNT(bmar) == 0 )
603 ERROR_exit("Can't create blurmaster subset collection? Help me, Mr. Wizard!") ;
604 }
605
606 { MRI_IMAGE *bbim = IMARR_SUBIM(bmar,0) ; /* 30 Dec 2015 */
607 bbim->dx = fabsf(DSET_DX(bmset)) ;
608 bbim->dy = fabsf(DSET_DY(bmset)) ;
609 bbim->dz = fabsf(DSET_DZ(bmset)) ;
610 }
611
612 if( !bmeqin ) DSET_unload(bmset) ;
613 for( ibm=0 ; ibm < IMARR_COUNT(bmar) ; ibm++ ){
614 IMARR_SUBIM(bmar,ibm)->dx = dx ;
615 IMARR_SUBIM(bmar,ibm)->dy = dy ;
616 IMARR_SUBIM(bmar,ibm)->dz = dz ;
617 }
618
619 /*--- pre-process input dataset ---*/
620
621 #if 1
622 if( corder_in < 0 ) corder_in = DSET_NVALS(inset) / 30 ;
623 #else
624 corder_in = 0 ;
625 #endif
626
627 DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
628
629 outset = EDIT_empty_copy( inset ) ; /* moved here 04 Jun 2007 */
630 EDIT_dset_items( outset , ADN_prefix , prefix , ADN_none ) ;
631 EDIT_dset_items( outset , ADN_brick_fac , NULL , ADN_none ) ; /* 11 Sep 2007 */
632 tross_Copy_History( inset , outset ) ;
633 tross_Make_History( "3dBlurToFWHM" , argc,argv , outset ) ;
634
635 if( corder_in > 0 ){ /* 04 Jun 2007 -- detrend input (but don't scale) */
636
637 int jj,iv,kk , nvi = DSET_NVALS(inset) , nzs ;
638 float tm,fac,fq , *vv,*ww , vmed ;
639 THD_3dim_dataset *newset ;
640
641 corder_inrefnum = 2*corder_in+3 ;
642 INFO_message("detrending input: %d ref funcs, %d time points",corder_inrefnum,nvi) ;
643
644 corder_inref = THD_build_trigref( corder_in , nvi ) ;
645 if( corder_inref == NULL ) ERROR_exit("THD_build_trigref failed!") ;
646
647 newset = THD_detrend_dataset( inset, corder_inrefnum,corder_inref ,
648 2,0 , mask , &corder_inar ) ;
649 if( newset == NULL ) ERROR_exit("detrending failed!?") ;
650 corder_invv = IMARR_SUBIM(corder_inar,corder_inrefnum) ;
651 vv = MRI_FLOAT_PTR(corder_invv) ;
652 ww = (float *)malloc(sizeof(float)*nvox) ;
653 for( ii=jj=0 ; ii < nvox ; ii++ ) if( mask[ii] ) ww[jj++] = vv[ii] ;
654 vmed = 0.05 * qmed_float(jj,ww) ; free(ww) ;
655 for( nzs=ii=0 ; ii < nvox ; ii++ )
656 if( mask[ii] && vv[ii] < vmed ){ mask[ii] = 0; nzs++ ; }
657 if( nzs > 0 )
658 INFO_message("Removed %d voxels with small input stdev from mask",nzs) ;
659
660 DSET_delete(inset) ; inset = newset ;
661 ININFO_message("detrending of input complete") ;
662
663 } else if( do_unif ){
664
665 int nvi = DSET_NVALS(inset) ;
666 if( nvi < 3 ){
667 WARNING_message(
668 "Can't apply -unif: only %d sub-brick%s in input dataset" ,
669 nvi , (nvi==1) ? "\0" : "s" ) ;
670 do_unif = 0 ;
671 } else {
672 imar = THD_medmad_bricks(inset) ;
673 imed = IMARR_SUBIM(imar,0) ; imedar = MRI_FLOAT_PTR(imed) ;
674 imad = IMARR_SUBIM(imar,1) ; imadar = MRI_FLOAT_PTR(imad) ;
675 for( ii=0 ; ii < nvox ; ii++ )
676 if( imadar[ii] != 0.0f ) imadar[ii] = 1.0f/imadar[ii] ;
677 }
678 }
679 INIT_IMARR(dsar) ;
680 for( ids=0 ; ids < DSET_NVALS(inset) ; ids++ ){
681 dsim = mri_scale_to_float(DSET_BRICK_FACTOR(inset,ids),DSET_BRICK(inset,ids));
682 if( do_unif ){
683 bar = MRI_FLOAT_PTR(dsim) ;
684 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] = (bar[ii]-imedar[ii])*imadar[ii] ;
685 }
686 ADDTO_IMARR(dsar,dsim) ; DSET_unload_one(inset,ids) ;
687 }
688 if( do_unif ){
689 for( ii=0 ; ii < nvox ; ii++ )
690 if( imadar[ii] != 0.0f ) imadar[ii] = 1.0f/imadar[ii] ;
691 }
692 DSET_unload(inset) ;
693 for( ids=0 ; ids < IMARR_COUNT(dsar) ; ids++ ){
694 IMARR_SUBIM(dsar,ids)->dx = dx ;
695 IMARR_SUBIM(dsar,ids)->dy = dy ;
696 IMARR_SUBIM(dsar,ids)->dz = dz ;
697 }
698
699 /*--- make arrays for FWHM estimates and then diffusion parameters ---*/
700
701 numfxyz = 0 ;
702 fxim = mri_new_conforming( IMARR_SUBIM(bmar,0) , MRI_float ) ;
703 fxar = MRI_FLOAT_PTR(fxim) ;
704 gx = fwhm_goal - 0.011f*dx ; numfxyz++ ;
705 hx = 0.8f*gx ; qx = 1.0f/(gx-hx) ;
706
707 fyim = mri_new_conforming( IMARR_SUBIM(bmar,0) , MRI_float ) ;
708 fyar = MRI_FLOAT_PTR(fyim) ;
709 gy = fwhm_goal - 0.011f*dy ; numfxyz++ ;
710 hy = 0.8f*gy ; qy = 1.0f/(gy-hy) ;
711
712 if( !fwhm_2D ){
713 fzim = mri_new_conforming( IMARR_SUBIM(bmar,0) , MRI_float ) ;
714 fzar = MRI_FLOAT_PTR(fzim) ;
715 gz = fwhm_goal - 0.011f*dz ; numfxyz++ ;
716 hz = 0.8f*gz ; qz = 1.0f/(gz-hz) ;
717 }
718 maxfxyz = blurfac / numfxyz ; /* maximum lambda */
719 #if 1
720 if( do_acf ) fwhm_goal *= 0.985f ; /* fudge factor */
721 #endif
722
723 /*----------- Do the work:
724 estimate smoothness of blur master
725 - exit loop when all done
726 set up incremental blurring parameters
727 blur the master and the input
728 loop back to top --------------*/
729
730 for( nite=1 ; nite <= maxite ; nite++ ){
731
732 /*--- global smoothness estimation into bx,by,bz ---*/
733
734 if( !do_acf ){ /* the old way */
735 fw = mriarr_estimate_FWHM_1dif( bmar , mask , do_unif ) ;
736 UNLOAD_FVEC3(fw,bx,by,bz) ;
737 } else { /* the ACF way [30 Dec 2015] */
738 bx = mriarr_estimate_FWHM_acf( bmar, mask, do_unif, 3.0f*fwhm_goal ) ;
739 by = bz = bx ; /* spherical assumption */
740 }
741
742 if( bx <= 0.0f ) bx = dx ; /* should not happen */
743 if( by <= 0.0f ) by = dy ;
744 if( bz <= 0.0f ) bz = dz ;
745
746 /*--- test for progress towards our goal ---*/
747
748 if( fwhm_2D ){ /*-- 2D work --*/
749 val = (float)sqrt(bx*by) ;
750 if( verb )
751 INFO_message("-- Iteration #%d: 2D FWHMx=%.4f FWHMy=%.4f sqrt()=%.4f",
752 nite,bx,by,val) ;
753 if( val >= fwhm_goal ){
754 if( verb ){
755 if( !do_acf )
756 INFO_message("** Passes 2D threshold sqrt(FWHMx*FWHMy) ==> done!");
757 else
758 INFO_message("** Passes 2D threshold for ACF FWHM ==> done!");
759 }
760 break;
761 }
762 if( nite > 3 && bx < last_fwx && by < last_fwy ){ /* negative progress? */
763 nbail++ ;
764 if( nbail == BAILOUT ){
765 if( verb ) INFO_message("** Bailing out for sluggish progress!") ;
766 break ;
767 } else {
768 if( verb ) INFO_message("** Progress is slow for some reason?!") ;
769 }
770 } else nbail = 0 ;
771 xdone = (bx >= fwhm_goal) ; /* is x done? */
772 ydone = (by >= fwhm_goal) ; /* is y done? */
773 zdone = 1 ;
774 xalmost = (!xdone && bx >= DFGOALC*fwhm_goal) ; /* is x 80% done? */
775 yalmost = (!ydone && by >= DFGOALC*fwhm_goal) ; /* is y 80% done? */
776 zalmost = 0 ;
777 xstall = (!xdone && (bx < last_fwx && nbail==BAILOUT)) ; /* is x going backwards? */
778 ystall = (!ydone && (by < last_fwy && nbail==BAILOUT)) ; /* is y going backwards? */
779 zstall = 1 ;
780 } else { /*--- 3D work ---*/
781 val = (float)cbrt(bx*by*bz) ;
782 if( verb ){
783 if( !do_acf )
784 INFO_message("-- Iteration #%d: 3D FWHMx=%.4f FWHMy=%.4f FWHMz=%.4f cbrt()=%.4f",
785 nite,bx,by,bz,val) ;
786 else
787 INFO_message("-- Iteration #%d: 3D ACF FWHM=%.4f",nite,val) ;
788 }
789 if( val >= fwhm_goal ){
790 if( verb ){
791 if( !do_acf ){
792 INFO_message("** Passes 3D threshold cbrt(FWHMx*FWHMy*FWHMz) ==> done!");
793 } else {
794 INFO_message("** Passes 3D threshold for ACF FWHM ==> done!") ;
795 }
796 }
797 break;
798 }
799 if( nite > 3 && bx < last_fwx && by < last_fwy && bz < last_fwz ){ /* negative progress? */
800 nbail++ ;
801 if( nbail == BAILOUT ){
802 if( verb ) INFO_message("** Bailing out for sluggish progress!") ;
803 break ;
804 } else {
805 if( verb ) INFO_message("** Progress is slow for some reason?!") ;
806 }
807 } else nbail = 0 ;
808 xdone = (bx >= fwhm_goal) ; /* is x done? */
809 ydone = (by >= fwhm_goal) ; /* is y done? */
810 zdone = (bz >= fwhm_goal) ; /* is z done? */
811 xalmost = (!xdone && bx >= DFGOALC*fwhm_goal) ; /* is x 80% done? */
812 yalmost = (!ydone && by >= DFGOALC*fwhm_goal) ; /* is y 80% done? */
813 zalmost = (!zdone && bz >= DFGOALC*fwhm_goal) ; /* is z 80% done? */
814 xstall = (!xdone && (bx < last_fwx && nbail==BAILOUT)) ; /* is x going backwards? */
815 ystall = (!ydone && (by < last_fwy && nbail==BAILOUT)) ; /* is y going backwards? */
816 zstall = (!zdone && (bz < last_fwz && nbail==BAILOUT)) ; /* is z going backwards? */
817 }
818
819 /* these cases should not happen, but just in case ... */
820 if( xdone && ydone && zdone ){
821 if( verb ) INFO_message("** All axes done!") ;
822 break ;
823 }
824 if( xstall && ystall && zstall ){
825 if( verb ) INFO_message("** All axes stalled ==> quitting in disgust") ;
826 break ;
827 }
828
829 /** if one axis is very close but others still need work, stop that one now **/
830
831 if( fwhm_2D && !xdone && !ydone ){
832 if( by > 1.01f*bx && by >= fwhm_subgoal ) ydone = 1; /* stop y */
833 else if( bx > 1.01f*by && bx >= fwhm_subgoal ) xdone = 1; /* stop y */
834 } else if( !fwhm_2D && !xdone && !ydone && !zdone ){
835 if( bz > 1.01f*by && bz > 1.01f*bx && bz >= fwhm_subgoal ) zdone = 1; /* stop z */
836 else if( by > 1.01f*bx && by > 1.01f*bz && by >= fwhm_subgoal ) ydone = 1; /* stop y */
837 else if( bx > 1.01f*by && bx > 1.01f*bz && bx >= fwhm_subgoal ) xdone = 1; /* stop x */
838 }
839
840 /** if global blurring is going too fast
841 (more than 20% progress per iteration), slow it down **/
842
843 if( nite > 1 ){
844 float bfac=1.0f , dfg=DFGOAL*fwhm_goal ;
845 delt_fwx = bx-last_fwx; delt_fwy = by-last_fwy; delt_fwz = bz-last_fwz;
846 if( !xdone && !xstall && delt_fwx > 0.0f && (val=dfg/delt_fwx) < 1.0f )
847 bfac = MIN( bfac , val ) ;
848 if( !ydone && !ystall && delt_fwy > 0.0f && (val=dfg/delt_fwy) < 1.0f )
849 bfac = MIN( bfac , val ) ;
850 if( !zdone && !zstall && delt_fwz > 0.0f && (val=dfg/delt_fwz) < 1.0f )
851 bfac = MIN( bfac , val ) ;
852 bfac = MAX(0.123f,bfac) ; blurfac *= bfac ;
853 maxfxyz = blurfac / numfxyz ;
854 if( bfac < 0.999f )
855 if( verb )
856 ININFO_message(" Slowing overall blur rate by factor of %.3f",bfac) ;
857 }
858
859 /** if near the goal in any direction, slow that direction's blurring down **/
860
861 maxfx = maxfxyz * xrat ;
862 if( xdone || xstall ){
863 maxfx = 0.0f; xstopped = 1;
864 if( verb ) ININFO_message(" x-blur is stopped");
865 } else if( xalmost ){
866 val = (nite==1) ? 0.123f : sqrt((fwhm_goal-bx)/(DFGOAL*fwhm_goal)) ;
867 if( val >= 0.0f && val < 0.999f ){
868 val = MAX(0.09f,val) ;
869 if( verb ) ININFO_message(" Slowing x-blur rate by factor of %.3f",val) ;
870 maxfx *= val ;
871 }
872 }
873
874 maxfy = maxfxyz * yrat ;
875 if( ydone || ystall ){
876 maxfy = 0.0f; ystopped = 1;
877 if( verb ) ININFO_message(" y-blur is stopped");
878 } else if( yalmost ){
879 val = (nite==1) ? 0.123f : sqrt((fwhm_goal-by)/(DFGOAL*fwhm_goal)) ;
880 if( val >= 0.0f && val < 0.999f ){
881 val = MAX(0.09f,val) ;
882 if( verb ) ININFO_message(" Slowing y-blur rate by factor of %.3f",val) ;
883 maxfy *= val ;
884 }
885 }
886
887 maxfz = maxfxyz * zrat ;
888 if( zdone || zstall ){
889 maxfz = 0.0f ; zstopped = 1;
890 if( !fwhm_2D && verb ) ININFO_message(" z-blur is stopped");
891 } else if( zalmost ){
892 val = (nite==1) ? 0.123f : sqrt((fwhm_goal-bz)/(DFGOAL*fwhm_goal)) ;
893 if( val >= 0.0f && val < 0.999f ){
894 val = MAX(0.09f,val) ;
895 if( verb ) ININFO_message(" Slowing z-blur rate by factor of %.3f",val) ;
896 maxfz *= val ;
897 }
898 }
899
900 last_fwx = bx; last_fwy = by; last_fwz = bz;
901
902 /*--- blur map estimation (within each nbhd) ---*/
903
904 if( nbhd != NULL ){
905 if( verb ) ININFO_message(" Estimate smoothness in each voxel's nbhd") ;
906 estimate_blur_map( bmar , mask,nbhd , fxar,fyar,fzar ) ;
907
908 /*--- Save blur maps for debuggers? ---*/
909
910 if( bsave_prefix != NULL ){
911 char bp[1024]; THD_3dim_dataset *bset; float *qqq;
912 { float *bbb=MRI_FLOAT_PTR( IMARR_SUBIM(bmar,0) ) ;
913 qqq = malloc(sizeof(float)*nvox); memcpy(qqq,bbb,sizeof(float)*nvox);
914 bset = EDIT_empty_copy(inset); sprintf(bp,"%s_%04d_bm",bsave_prefix,nite);
915 EDIT_dset_items( bset, ADN_prefix,bp, ADN_nvals,1, ADN_ntt,0, ADN_none );
916 EDIT_substitute_brick( bset , 0 , MRI_float , qqq ) ;
917 EDIT_BRICK_LABEL( bset , 0 , "Bmar0" ) ;
918 DSET_write(bset) ; WROTE_DSET(bset) ; DSET_delete(bset) ;
919 }
920 if( fxar != NULL ){
921 qqq = malloc(sizeof(float)*nvox); memcpy(qqq,fxar,sizeof(float)*nvox);
922 bset = EDIT_empty_copy(inset); sprintf(bp,"%s_%04dx",bsave_prefix,nite);
923 EDIT_dset_items( bset, ADN_prefix,bp, ADN_nvals,1, ADN_ntt,0, ADN_none );
924 EDIT_substitute_brick( bset , 0 , MRI_float , qqq ) ;
925 EDIT_BRICK_LABEL( bset , 0 , "FWHMx" ) ;
926 DSET_write(bset) ; WROTE_DSET(bset) ; DSET_delete(bset) ;
927 }
928 if( fyar != NULL ){
929 qqq = malloc(sizeof(float)*nvox); memcpy(qqq,fyar,sizeof(float)*nvox);
930 bset = EDIT_empty_copy(inset); sprintf(bp,"%s_%04dy",bsave_prefix,nite);
931 EDIT_dset_items( bset, ADN_prefix,bp, ADN_nvals,1, ADN_ntt,0, ADN_none );
932 EDIT_substitute_brick( bset , 0 , MRI_float , qqq ) ;
933 EDIT_BRICK_LABEL( bset , 0 , "FWHMy" ) ;
934 DSET_write(bset) ; WROTE_DSET(bset) ; DSET_delete(bset) ;
935 }
936 if( fzar != NULL ){
937 qqq = malloc(sizeof(float)*nvox); memcpy(qqq,fzar,sizeof(float)*nvox);
938 bset = EDIT_empty_copy(inset); sprintf(bp,"%s_%04dz",bsave_prefix,nite);
939 EDIT_dset_items( bset, ADN_prefix,bp, ADN_nvals,1, ADN_ntt,0, ADN_none );
940 EDIT_substitute_brick( bset , 0 , MRI_float , qqq ) ;
941 EDIT_BRICK_LABEL( bset , 0 , "FWHMz" ) ;
942 DSET_write(bset) ; WROTE_DSET(bset) ; DSET_delete(bset) ;
943 }
944 }
945
946 /*--- compute local blurring parameters where needed ---*/
947
948 /* temper the blurring rate depending on the local smoothness? */
949
950 #undef TB
951 #undef TT
952 #undef TR
953 #define TB 0.8f
954 #define TT 0.2f
955 #define TR 0.05f
956 if( temperize ){
957 float_pair qp; MRI_IMAGE *qim; float *qar; int ng;
958 temper_fx = temper_fy = temper_fz = 0 ;
959 if( fxim != NULL && !xstopped ){
960 for( ii=ng=0 ; ii < nvox ; ii++ )
961 if( mask[ii] == 1 && fxar[ii] > 0.0f && fxar[ii] < gx ) ng++ ;
962 if( ng >= 666 ){
963 qim = mri_new_vol( ng,1,1 , MRI_float ); qar = MRI_FLOAT_PTR(qim);
964 for( ii=ng=0 ; ii < nvox ; ii++ )
965 if( mask[ii] == 1 && fxar[ii] > 0.0f && fxar[ii] < gx )
966 qar[ng++] = fxar[ii] ;
967 qp = mri_twoquantiles( qim , 0.15f , 0.85f ) ;
968 fx_tbot = qp.a ; fx_ttop = qp.b ; mri_free(qim) ;
969 fx_trat = fx_ttop - fx_tbot ;
970 temper_fx = (fx_trat > 0.0f) && (fx_trat > TR*gx) ;
971 if( temper_fx ) fx_trat = TB / fx_trat ;
972 if( verb ) ININFO_message(" temper x: %.4f .. %.4f %s",
973 fx_tbot,fx_ttop,temper_fx?"ON":"OFF") ;
974 }
975 }
976 if( fyim != NULL && !ystopped ){
977 for( ii=ng=0 ; ii < nvox ; ii++ )
978 if( mask[ii] == 1 && fyar[ii] > 0.0f && fyar[ii] < gy ) ng++ ;
979 if( ng >= 666 ){
980 qim = mri_new_vol( ng,1,1 , MRI_float ); qar = MRI_FLOAT_PTR(qim);
981 for( ii=ng=0 ; ii < nvox ; ii++ )
982 if( mask[ii] == 1 && fyar[ii] > 0.0f && fyar[ii] < gy )
983 qar[ng++] = fyar[ii] ;
984 qp = mri_twoquantiles( qim , 0.15f , 0.85f ) ;
985 fy_tbot = qp.a ; fy_ttop = qp.b ; mri_free(qim) ;
986 fy_trat = fy_ttop - fy_tbot ;
987 temper_fy = (fy_trat > 0.0f) && (fy_trat > TR*gy) ;
988 if( verb ) ININFO_message(" temper y: %.4f .. %.4f %s",
989 fy_tbot,fy_ttop,temper_fy?"ON":"OFF") ;
990 }
991 }
992 if( fzim != NULL && !zstopped ){
993 for( ii=ng=0 ; ii < nvox ; ii++ )
994 if( mask[ii] == 1 && fzar[ii] > 0.0f && fzar[ii] < gz ) ng++ ;
995 if( ng >= 666 ){
996 qim = mri_new_vol( ng,1,1 , MRI_float ); qar = MRI_FLOAT_PTR(qim);
997 for( ii=ng=0 ; ii < nvox ; ii++ )
998 if( mask[ii] == 1 && fzar[ii] > 0.0f && fzar[ii] < gz )
999 qar[ng++] = fzar[ii] ;
1000 qp = mri_twoquantiles( qim , 0.15f , 0.85f ) ;
1001 fz_tbot = qp.a ; fz_ttop = qp.b ; mri_free(qim) ;
1002 fz_trat = fz_ttop - fz_tbot ;
1003 temper_fz = (fz_trat > 0.0f) && (fz_trat > TR*gz) ;
1004 if( verb ) ININFO_message(" temper z: %.4f .. %.4f %s",
1005 fz_tbot,fz_ttop,temper_fz?"ON":"OFF") ;
1006 }
1007 }
1008 }
1009
1010 nblur = 0 ;
1011 for( ii=0 ; ii < nvox ; ii++ ){
1012 if( mask[ii] == 1 ){
1013 nd = 0 ;
1014 if( fxar != NULL ){
1015 if( fxar[ii] <= 0.0f || fxar[ii] >= gx )
1016 fxar[ii] = 0.0f ;
1017 else if( !temper_fx )
1018 fxar[ii] = (fxar[ii] <= hx) ? maxfx : maxfx*qx*(gx-fxar[ii]) ;
1019 else if( fxar[ii] <= fx_tbot )
1020 fxar[ii] = maxfx ;
1021 else if( fxar[ii] >= fx_ttop )
1022 fxar[ii] = TT*maxfx ;
1023 else
1024 fxar[ii] = (TT+fx_trat*(fx_ttop-fxar[ii]))*maxfx ;
1025
1026 if( fxar[ii] > 0.0f ) nd++ ;
1027 }
1028 if( fyar != NULL ){
1029 if( fyar[ii] <= 0.0f || fyar[ii] >= gy )
1030 fyar[ii] = 0.0f ;
1031 else if( !temper_fy )
1032 fyar[ii] = (fyar[ii] <= hy) ? maxfy : maxfy*qy*(gy-fyar[ii]) ;
1033 else if( fyar[ii] <= fy_tbot )
1034 fyar[ii] = maxfy ;
1035 else if( fyar[ii] >= fy_ttop )
1036 fyar[ii] = TT*maxfy ;
1037 else
1038 fyar[ii] = (TT+fy_trat*(fy_ttop-fyar[ii]))*maxfy ;
1039
1040 if( fyar[ii] > 0.0f ) nd++ ;
1041 }
1042 if( fzar != NULL ){
1043 if( fzar[ii] <= 0.0f || fzar[ii] >= gz )
1044 fzar[ii] = 0.0f ;
1045 else if( !temper_fz )
1046 fzar[ii] = (fzar[ii] <= hz) ? maxfz : maxfz*qz*(gz-fzar[ii]) ;
1047 else if( fzar[ii] <= fz_tbot )
1048 fzar[ii] = maxfz ;
1049 else if( fzar[ii] >= fz_ttop )
1050 fzar[ii] = TT*maxfz ;
1051 else
1052 fzar[ii] = (TT+fz_trat*(fz_ttop-fzar[ii]))*maxfz ;
1053
1054 if( fzar[ii] > 0.0f ) nd++ ;
1055 }
1056 if( nd == 0 ) mask[ii] = 2 ; /* turn off future diffusion here */
1057 else nblur++ ;
1058 }
1059 }
1060 if( nblur == 0 ) break ; /* no blurring ==> done?! */
1061
1062 } else { /* no nbhd ==> blur all non-mask points equally */
1063 nblur = 0 ;
1064 for( ii=0 ; ii < nvox ; ii++ ){
1065 if( fxar != NULL ) fxar[ii] = (mask[ii]) ? maxfx : 0.0f ;
1066 if( fyar != NULL ) fyar[ii] = (mask[ii]) ? maxfy : 0.0f ;
1067 if( fzar != NULL ) fzar[ii] = (mask[ii]) ? maxfz : 0.0f ;
1068 if( mask[ii] ) nblur++ ;
1069 }
1070 }
1071
1072 /*--- blur the master and the input ---*/
1073
1074 if( verb ) ININFO_message(" Blurring %d voxels in master",nblur) ;
1075 AFNI_OMP_START ;
1076 #pragma omp parallel if( IMARR_COUNT(bmar) > 2 )
1077 { int ii ;
1078 #pragma omp for
1079 for( ii=0 ; ii < IMARR_COUNT(bmar) ; ii++ )
1080 mri_blur3D_variable( IMARR_SUBIM(bmar,ii) , mask , fxim,fyim,fzim ) ;
1081 }
1082 AFNI_OMP_END ;
1083
1084 if( verb ) ININFO_message(" Blurring input") ;
1085 AFNI_OMP_START ;
1086 #pragma omp parallel if( IMARR_COUNT(dsar) > 2 )
1087 { int ii ;
1088 #pragma omp for
1089 for( ii=0 ; ii < IMARR_COUNT(dsar) ; ii++ )
1090 mri_blur3D_variable( IMARR_SUBIM(dsar,ii) , mask , fxim,fyim,fzim ) ;
1091 }
1092 AFNI_OMP_END ;
1093
1094 } /*-- loop back to estimate smoothness, check done-ness, smooth, etc --*/
1095
1096 /*----- tossolo some trashola ---*/
1097
1098 DESTROY_IMARR(bmar) ;
1099
1100 if( fxim != NULL ) mri_free(fxim) ;
1101 if( fyim != NULL ) mri_free(fyim) ;
1102 if( fzim != NULL ) mri_free(fzim) ;
1103
1104 /*--- put blurred data into the output dataset ---*/
1105
1106 for( ids=0 ; ids < DSET_NVALS(outset) ; ids++ ){
1107 bar = MRI_FLOAT_PTR( IMARR_SUBIM(dsar,ids) ) ;
1108 if( do_unif ){
1109 for( ii=0 ; ii < nvox ; ii++ )
1110 if( mask[ii] ) bar[ii] = bar[ii]*imadar[ii] + imedar[ii] ;
1111 }
1112 EDIT_substitute_brick( outset , ids , MRI_float , bar ) ;
1113 }
1114 #if 1
1115 if( corder_inrefnum > 0 ){
1116 INFO_message("re-trending output dataset") ;
1117 THD_retrend_dataset( outset , corder_inrefnum,corder_inref ,
1118 0 , mask , corder_inar ) ;
1119 }
1120 #endif
1121 DSET_write(outset) ;
1122 WROTE_DSET(outset) ;
1123
1124 /*--- run 3dFWHMx for fun and profit(?) ---*/
1125
1126 if( verb ){
1127 char *buf , *pg ;
1128 pg = THD_find_executable("3dFWHMx") ;
1129 buf = malloc( ((pg != NULL) ? strlen(pg)+999 : 999 ) ) ;
1130 sprintf(buf,"%s", ((pg != NULL) ? pg : "3dFWHMx") ) ;
1131 if( !do_acf ) sprintf(buf+strlen(buf)," -arith") ;
1132 else sprintf(buf+strlen(buf)," -acf %s.1D %.3g",prefix,3.0*fwhm_goal ) ;
1133 if( do_unif )
1134 sprintf(buf+strlen(buf)," -unif") ;
1135 else if( corder_bm > 0 )
1136 sprintf(buf+strlen(buf)," -detrend") ;
1137 if( automask )
1138 sprintf(buf+strlen(buf)," -automask") ;
1139 else if( mset != NULL )
1140 sprintf(buf+strlen(buf)," -mask %s",DSET_BRIKNAME(mset)) ;
1141 sprintf(buf+strlen(buf)," -input %s",DSET_BRIKNAME(outset)) ;
1142 if( pg != NULL ){
1143 INFO_message("-------------------------------------------------------") ;
1144 ININFO_message("Checking results by running command below:") ;
1145 ININFO_message(" %s",buf) ;
1146 (void)system(buf) ;
1147 INFO_message("-------------------------------------------------------") ;
1148 } else {
1149 INFO_message("-------------------------------------------------------") ;
1150 ININFO_message("To check results, run the command below:") ;
1151 ININFO_message(" %s",buf) ;
1152 INFO_message("-------------------------------------------------------") ;
1153 }
1154 free(buf) ;
1155 }
1156 exit(0) ;
1157 }
1158
1159 /*----------------------------------------------------------------------------*/
1160
1161 #undef ASSIF
1162 #define ASSIF(f,i,v) if(f != NULL)f[i]=v
1163
estimate_blur_map(MRI_IMARR * bmar,byte * mask,MCW_cluster * nbhd,float * fxar,float * fyar,float * fzar)1164 void estimate_blur_map( MRI_IMARR *bmar , byte *mask , MCW_cluster *nbhd ,
1165 float *fxar , float *fyar , float *fzar )
1166 {
1167 int ibm , ii , nvox , nar , nx,nxy , xx,yy,zz , pp ;
1168 THD_fvec3 fw ;
1169 float fxx , fyy , fzz ;
1170 int nfxx , nfyy , nfzz;
1171 float vxx , vyy , vzz;
1172
1173 ENTRY("estimate_blur_map") ;
1174
1175 nvox = IMARR_SUBIM(bmar,0)->nvox ;
1176 nar = IMARR_COUNT(bmar) ;
1177 nx = IMARR_SUBIM(bmar,0)->nx ;
1178 nxy = nx * IMARR_SUBIM(bmar,0)->ny ;
1179
1180 for( ii=0 ; ii < nvox ; ii++ ){
1181 vxx = vyy = vzz = 0.0f ;
1182 if( mask[ii] == 1 ){
1183 nfxx = nfyy = nfzz = 0 ;
1184 fxx = fyy = fzz = 0.0f ;
1185 IJK_TO_THREE(ii,xx,yy,zz,nx,nxy) ;
1186 for( ibm=0 ; ibm < nar ; ibm++ ){
1187 fw = mri_nstat_fwhmxyz( xx,yy,zz , IMARR_SUBIM(bmar,ibm),mask,nbhd );
1188 if( fxar != NULL && fw.xyz[0] > 0.0f ){ nfxx++; fxx += fw.xyz[0]; }
1189 if( fyar != NULL && fw.xyz[1] > 0.0f ){ nfyy++; fyy += fw.xyz[1]; }
1190 if( fzar != NULL && fw.xyz[2] > 0.0f ){ nfzz++; fzz += fw.xyz[2]; }
1191 }
1192 if( nfxx > 0 ) vxx = fxx / nfxx ;
1193 if( nfyy > 0 ) vyy = fyy / nfyy ;
1194 if( nfzz > 0 ) vzz = fzz / nfzz ;
1195 }
1196 ASSIF(fxar,ii,vxx); ASSIF(fyar,ii,vyy); ASSIF(fzar,ii,vzz);
1197 }
1198
1199 EXRETURN;
1200 }
1201