1 #include "mrilib.h"
2
3 /*------------------------------------------------------------------------*/
4
main(int argc,char * argv[])5 int main( int argc , char * argv[] )
6 {
7 float mrad=0.0f , fwhm=0.0f ;
8 int nrep=1 ;
9 char *prefix = "Polyfit" ;
10 char *resid = NULL ;
11 char *cfnam = NULL ;
12 int iarg , verb=0 , do_automask=0 , nord=3 , meth=2 , do_mclip=0 ;
13 THD_3dim_dataset *inset ;
14 MRI_IMAGE *imout , *imin ;
15 byte *mask=NULL ; int nvmask=0 , nmask=0 , do_mone=0 , do_byslice=0 ;
16 MRI_IMARR *exar=NULL ;
17 floatvec *fvit=NULL ; /* 26 Feb 2019 */
18
19 if( argc < 2 || strcasecmp(argv[1],"-help") == 0 ){
20 printf("\n"
21 "Usage: 3dPolyfit [options] dataset ~1~\n"
22 "\n"
23 "* Fits a polynomial in space to the input dataset and outputs that fitted dataset.\n"
24 "\n"
25 "* You can also add your own basis datasets to the fitting mix, using the\n"
26 " '-base' option.\n"
27 "\n"
28 "* You can get the fit coefficients using the '-1Dcoef' option.\n"
29 "\n"
30 "--------\n"
31 "Options: ~1~\n"
32 "--------\n"
33 "\n"
34 " -nord n = Maximum polynomial order (0..9) [default order=3]\n"
35 " [n=0 is the constant 1]\n"
36 " [n=-1 means only use volumes from '-base']\n"
37 "\n"
38 " -blur f = Gaussian blur input dataset (inside mask) with FWHM='f' (mm)\n"
39 "\n"
40 " -mrad r = Radius (voxels) of preliminary median filter of input\n"
41 " [default is no blurring of either type; you can]\n"
42 " [do both types (Gaussian and median), but why??]\n"
43 " [N.B.: median blur is slower than Gaussian]\n"
44 "\n"
45 " -prefix pp = Use 'pp' for prefix of output dataset (the fit).\n"
46 " [default prefix is 'Polyfit'; use NULL to skip this output]\n"
47 "\n"
48 " -resid rr = Use 'rr' for the prefix of the residual dataset.\n"
49 " [default is not to output residuals]\n"
50 "\n"
51 " -1Dcoef cc = Save coefficients of fit into text file cc.1D.\n"
52 " [default is not to save these coefficients]\n"
53 "\n"
54 " -automask = Create a mask (a la 3dAutomask)\n"
55 " -mask mset = Create a mask from nonzero voxels in 'mset'.\n"
56 " [default is not to use a mask, which is probably a bad idea]\n"
57 "\n"
58 " -mone = Scale the mean value of the fit (inside the mask) to 1.\n"
59 " [probably this option is not useful for anything]\n"
60 "\n"
61 " -mclip = Clip fit values outside the rectilinear box containing the\n"
62 " mask to the edge of that box, to avoid weird artifacts.\n"
63 "\n"
64 " -meth mm = Set 'mm' to 2 for least squares fit;\n"
65 " set it to 1 for L1 fit [default method=2]\n"
66 " [Note that L1 fitting is slower than L2 fitting!]\n"
67 "\n"
68 " -base bb = In addition to the polynomial fit, also use\n"
69 " the volumes in dataset 'bb' as extra basis functions.\n"
70 " [If you use a base dataset, then you can set nord]\n"
71 " [to -1, to skip using any spatial polynomial fit.]\n"
72 "\n"
73 " -verb = Print fun and useful progress reports :-)\n"
74 "\n"
75 "------\n"
76 "Notes: ~1~\n"
77 "------\n"
78 "* Output dataset is always stored in float format.\n"
79 "\n"
80 "* If the input dataset has more than 1 sub-brick, only sub-brick #0\n"
81 " is processed. To fit more than one volume, you'll have to use a script\n"
82 " to loop over the input sub-bricks, and then glue (3dTcat) the results\n"
83 " together to get a final result. A simple example:\n"
84 " #!/bin/tcsh\n"
85 " set base = model.nii\n"
86 " set dset = errts.nii\n"
87 " set nval = `3dnvals $dset`\n"
88 " @ vtop = $nval - 1\n"
89 " foreach vv ( `count 0 $vtop` )\n"
90 " 3dPolyfit -base \"$base\" -nord 0 -mask \"$base\" -1Dcoef QQ.$vv -prefix QQ.$vv.nii $dset\"[$vv]\"\n"
91 " end\n"
92 " 3dTcat -prefix QQall.nii QQ.0*.nii\n"
93 " 1dcat QQ.0*.1D > QQall.1D\n"
94 " \rm QQ.0*\n"
95 " exit 0\n"
96 "\n"
97 "* If the '-base' dataset has multiple sub-bricks, all of them are used.\n"
98 "\n"
99 "* You can use the '-base' option more than once, if desired or needed.\n"
100 "\n"
101 "* The original motivation for this program was to fit a spatial model\n"
102 " to a field map MRI, but that didn't turn out to be useful. Nevertheless,\n"
103 " I make this program available to someone who might find it beguiling.\n"
104 "\n"
105 "* If you really want, I could allow you to put sign constraints on the\n"
106 " fit coefficients (e.g., say that the coefficient for a given base volume\n"
107 " should be non-negative). But you'll have to beg for this.\n"
108 "\n"
109 "-- Emitted by RWCox\n"
110 ) ;
111 PRINT_COMPILE_DATE ; exit(0) ;
112 }
113
114 /*-- startup paperwork --*/
115
116 mainENTRY("3dPolyfit main"); machdep(); AFNI_logger("3dPolyfit",argc,argv);
117 PRINT_VERSION("3dPolyfit") ;
118
119 /*-- scan command line --*/
120
121 iarg = 1 ;
122 while( iarg < argc && argv[iarg][0] == '-' ){
123
124 if( strcasecmp(argv[iarg],"-base") == 0 ){
125 THD_3dim_dataset *bset ; int kk ; MRI_IMAGE *bim ;
126 if( ++iarg >= argc ) ERROR_exit("Need argument after '-base'") ;
127 bset = THD_open_dataset(argv[iarg]) ;
128 CHECK_OPEN_ERROR(bset,argv[iarg]) ;
129 DSET_load(bset) ; CHECK_LOAD_ERROR(bset) ;
130 if( exar == NULL ) INIT_IMARR(exar) ;
131 for( kk=0 ; kk < DSET_NVALS(bset) ; kk++ ){
132 bim = THD_extract_float_brick(kk,bset) ;
133 if( bim != NULL ) ADDTO_IMARR(exar,bim) ;
134 DSET_unload_one(bset,kk) ;
135 }
136 DSET_delete(bset) ;
137 iarg++ ; continue ;
138 }
139
140 if( strcasecmp(argv[iarg],"-verb") == 0 ){
141 verb++ ; iarg++ ; continue ;
142 }
143
144 if( strcasecmp(argv[iarg],"-hermite") == 0 ){ /* 25 Mar 2013 [New Year's Day] */
145 mri_polyfit_set_basis("hermite") ; /* HIDDEN */
146 iarg++ ; continue ;
147 }
148
149 if( strcasecmp(argv[iarg],"-byslice") == 0 ){ /* 25 Mar 2013 [New Year's Day] */
150 do_byslice++ ; iarg++ ; continue ; /* HIDDEN */
151 }
152
153 if( strcasecmp(argv[iarg],"-mask") == 0 ){
154 THD_3dim_dataset *mset ;
155 if( ++iarg >= argc ) ERROR_exit("Need argument after '-mask'") ;
156 if( mask != NULL || do_automask ) ERROR_exit("Can't have two mask inputs") ;
157 mset = THD_open_dataset(argv[iarg]) ;
158 CHECK_OPEN_ERROR(mset,argv[iarg]) ;
159 DSET_load(mset) ; CHECK_LOAD_ERROR(mset) ;
160 nvmask = DSET_NVOX(mset) ;
161 mask = THD_makemask( mset , 0 , 0.5f, 0.0f ) ; DSET_delete(mset) ;
162 if( mask == NULL ) ERROR_exit("Can't make mask from dataset '%s'",argv[iarg]) ;
163 nmask = THD_countmask( nvmask , mask ) ;
164 if( nmask < 99 ) ERROR_exit("Too few voxels in mask (%d)",nmask) ;
165 if( verb ) INFO_message("Number of voxels in mask = %d",nmask) ;
166 iarg++ ; continue ;
167 }
168
169 if( strcasecmp(argv[iarg],"-nord") == 0 ){
170 nord = (int)strtol( argv[++iarg], NULL , 10 ) ;
171 if( nord < -1 || nord > 9 )
172 ERROR_exit("Illegal value after -nord :(") ;
173 iarg++ ; continue ;
174 }
175
176 if( strcasecmp(argv[iarg],"-meth") == 0 ){
177 meth = (int)strtol( argv[++iarg], NULL , 10 ) ;
178 if( meth < 1 || meth > 2 )
179 ERROR_exit("Illegal value after -meth :(") ;
180 iarg++ ; continue ;
181 }
182
183 if( strncmp(argv[iarg],"-automask",5) == 0 ){
184 if( mask != NULL ) ERROR_exit("Can't use -mask and -automask together!") ;
185 do_automask++ ; iarg++ ; continue ;
186 }
187
188 if( strncmp(argv[iarg],"-mclip",5) == 0 ){
189 do_mclip++ ; iarg++ ; continue ;
190 }
191
192 if( strncmp(argv[iarg],"-mone",5) == 0 ){
193 do_mone++ ; iarg++ ; continue ;
194 }
195
196 if( strcasecmp(argv[iarg],"-mrad") == 0 ){
197 mrad = strtod( argv[++iarg] , NULL ) ; iarg++ ; continue ;
198 }
199
200 if( strcasecmp(argv[iarg],"-blur") == 0 ){
201 fwhm = strtod( argv[++iarg] , NULL ) ; iarg++ ; continue ;
202 }
203
204 if( strcasecmp(argv[iarg],"-prefix") == 0 ){
205 prefix = argv[++iarg] ;
206 if( !THD_filename_ok(prefix) )
207 ERROR_exit("Illegal value after -prefix :(");
208 if( strcasecmp(prefix,"NULL") == 0 ) prefix = NULL ;
209 iarg++ ; continue ;
210 }
211
212 if( strcasecmp(argv[iarg],"-resid") == 0 ){
213 resid = argv[++iarg] ;
214 if( !THD_filename_ok(resid) )
215 ERROR_exit("Illegal value after -resid :(");
216 if( strcasecmp(resid,"NULL") == 0 ) resid = NULL ;
217 iarg++ ; continue ;
218 }
219
220 if( strcasecmp(argv[iarg],"-1Dcoef") == 0 ){ /* 26 Feb 2019 */
221 cfnam = argv[++iarg] ;
222 if( !THD_filename_ok(cfnam) )
223 ERROR_exit("Illegal value after -1Dcoef :(");
224 if( strcasecmp(cfnam,"NULL") == 0 ) cfnam = NULL ;
225 iarg++ ; continue ;
226 }
227
228 ERROR_exit("Unknown option: %s\n",argv[iarg]);
229 }
230
231 /*--- check for blatant errors ---*/
232
233 if( iarg >= argc )
234 ERROR_exit("No input dataset name on command line?");
235
236 if( prefix == NULL && resid == NULL && cfnam == NULL )
237 ERROR_exit("-prefix and -resid and -1Dcoef are all NULL?!") ;
238
239 if( do_byslice && cfnam != NULL ){
240 WARNING_message("-byslice does not work with -1Dcoef option :(") ;
241 cfnam = NULL ;
242 }
243
244 if( nord < 0 && exar == NULL )
245 ERROR_exit("no polynomial fit AND no -base option ==> nothing to compute :(") ;
246
247 /*-- read input --*/
248
249 if( verb ) INFO_message("Load input dataset") ;
250
251 inset = THD_open_dataset( argv[iarg] ) ;
252 CHECK_OPEN_ERROR(inset,argv[iarg]) ;
253 DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
254 if( DSET_NVALS(inset) > 1 )
255 WARNING_message( "Only processing sub-brick #0 (out of %d)" , DSET_NVALS(inset) );
256
257 /* check input mask or create automask */
258
259 if( mask != NULL ){
260 if( nvmask != DSET_NVOX(inset) )
261 ERROR_exit("-mask and input datasets don't match in voxel counts :-(") ;
262 } else if( do_automask ){
263 THD_automask_verbose( (verb > 1) ) ;
264 THD_automask_extclip( 1 ) ;
265 mask = THD_automask( inset ) ; nvmask = DSET_NVOX(inset) ;
266 nmask = THD_countmask( nvmask , mask ) ;
267 if( nmask < 99 ) ERROR_exit("Too few voxels in automask (%d)",nmask) ;
268 if( verb ) ININFO_message("Number of voxels in automask = %d",nmask) ;
269 } else {
270 WARNING_message("3dPolyfit is running without a mask") ;
271 }
272
273 #undef GOOD
274 #define GOOD(i) (mask == NULL || mask[i])
275
276 /* check -base input datasets */
277
278 if( exar != NULL ){
279 int ii,kk , nvbad=0 , nvox=DSET_NVOX(inset),nm ; float *ex , exb ;
280 for( kk=0 ; kk < IMARR_COUNT(exar) ; kk++ ){
281 if( nvox != IMARR_SUBIM(exar,kk)->nvox ){
282 if( IMARR_SUBIM(exar,kk)->nvox != nvbad ){
283 ERROR_message("-base volume (%d voxels) doesn't match input dataset grid size (%d voxels)",
284 IMARR_SUBIM(exar,kk)->nvox , nvox ) ;
285 nvbad = IMARR_SUBIM(exar,kk)->nvox ;
286 }
287 }
288 }
289 if( nvbad != 0 ) ERROR_exit("Cannot continue :-(") ;
290
291 /* subtract mean from each base input, if is a constant polynomial in the fit */
292
293 if( nord >= 0 ){
294 if( verb ) INFO_message("subtracting spatial mean from '-base'") ;
295 for( kk=0 ; kk < IMARR_COUNT(exar) ; kk++ ){
296 exb = 0.0f ; ex = MRI_FLOAT_PTR(IMARR_SUBIM(exar,kk)) ;
297 for( nm=ii=0 ; ii < nvox ; ii++ ){ if( GOOD(ii) ){ exb += ex[ii]; nm++; } }
298 exb /= nm ;
299 for( ii=0 ; ii < nvox ; ii++ ) ex[ii] -= exb ;
300 }
301 }
302 }
303
304 /* if blurring, edit mask a little */
305
306 if( mask != NULL && (fwhm > 0.0f || mrad > 0.0f) ){
307 int ii ;
308 ii = THD_mask_remove_isolas( DSET_NX(inset),DSET_NY(inset),DSET_NZ(inset),mask ) ;
309 if( ii > 0 ){
310 nmask = THD_countmask( nvmask , mask ) ;
311 if( verb )
312 ININFO_message("Removed %d isola%s from mask, leaving %d voxels" ,
313 ii,(ii==1)?"\0":"s" , nmask ) ;
314 if( nmask < 99 )
315 ERROR_exit("Too few voxels left in mask after isola removal :-(") ;
316 }
317 }
318
319 /* convert input to float, which is simpler to deal with */
320
321 imin = THD_extract_float_brick(0,inset) ;
322 if( imin == NULL ) ERROR_exit("Can't extract input dataset brick?! :-(") ;
323 DSET_unload(inset) ;
324
325 if( verb ) INFO_message("Start fitting process") ;
326
327 /* do the Gaussian blurring */
328
329 if( fwhm > 0.0f ){
330 if( verb ) ININFO_message("Gaussian blur: FWHM=%g mm",fwhm) ;
331 imin->dx = fabsf(DSET_DX(inset)) ;
332 imin->dy = fabsf(DSET_DY(inset)) ;
333 imin->dz = fabsf(DSET_DZ(inset)) ;
334 mri_blur3D_addfwhm( imin , mask , fwhm ) ;
335 }
336
337 /* do the fitting */
338
339 mri_polyfit_verb(verb) ;
340 if( do_byslice )
341 imout = mri_polyfit_byslice( imin , nord , exar , mask , mrad , meth ) ;
342 else
343 imout = mri_polyfit ( imin , nord , exar , mask , mrad , meth ) ;
344
345 /* WTF? */
346
347 if( imout == NULL )
348 ERROR_exit("Can't compute polynomial fit :-( !?") ;
349 if( resid == NULL ) mri_free(imin) ;
350
351 if( ! do_byslice )
352 fvit = mri_polyfit_get_fitvec() ; /* get coefficients of fit [26 Feb 2019] */
353
354 /* scale the fit dataset? */
355
356 if( do_mone ){
357 float sum=0.0f ; int nsum=0 , ii,nvox ; float *par=MRI_FLOAT_PTR(imout) ;
358 nvox = imout->nvox ;
359 for( ii=0 ; ii < nvox ; ii++ ){
360 if( mask != NULL && mask[ii] == 0 ) continue ;
361 sum += par[ii] ; nsum++ ;
362 }
363 if( nsum > 0 && sum != 0.0f ){
364 sum = nsum / sum ;
365 if( verb ) ININFO_message("-mone: scaling fit by %g",sum) ;
366 for( ii=0 ; ii < nvox ; ii++ ) par[ii] *= sum ;
367 }
368 }
369
370 /* if there's a mask, clip values outside of its box */
371
372 #undef PF
373 #define PF(i,j,k) par[(i)+(j)*nx+(k)*nxy]
374 if( mask != NULL && do_mclip ){
375 int xm,xp,ym,yp,zm,zp , ii,jj,kk , nx,ny,nz,nxy ; float *par ;
376 MRI_IMAGE *bim = mri_empty_conforming( imout , MRI_byte ) ;
377 mri_fix_data_pointer(mask,bim) ;
378 if( verb ) ININFO_message("-mclip: polynomial fit to autobox of mask") ;
379 MRI_autobbox( bim , &xm,&xp , &ym,&yp , &zm,&zp ) ;
380 mri_clear_data_pointer(bim) ; mri_free(bim) ;
381 nx = imout->nx ; ny = imout->ny ; nz = imout->nz ; nxy = nx*ny ;
382 par = MRI_FLOAT_PTR(imout) ;
383 for( ii=0 ; ii < xm ; ii++ )
384 for( kk=0 ; kk < nz ; kk++ )
385 for( jj=0 ; jj < ny ; jj++ ) PF(ii,jj,kk) = PF(xm,jj,kk) ;
386 for( ii=xp+1 ; ii < nx ; ii++ )
387 for( kk=0 ; kk < nz ; kk++ )
388 for( jj=0 ; jj < ny ; jj++ ) PF(ii,jj,kk) = PF(xp,jj,kk) ;
389 for( jj=0 ; jj < ym ; jj++ )
390 for( kk=0 ; kk < nz ; kk++ )
391 for( ii=0 ; ii < nx ; ii++ ) PF(ii,jj,kk) = PF(ii,ym,kk) ;
392 for( jj=yp+1 ; jj < ny ; jj++ )
393 for( kk=0 ; kk < nz ; kk++ )
394 for( ii=0 ; ii < nx ; ii++ ) PF(ii,jj,kk) = PF(ii,yp,kk) ;
395 for( kk=0 ; kk < zm ; kk++ )
396 for( jj=0 ; jj < ny ; jj++ )
397 for( ii=0 ; ii < nx ; ii++ ) PF(ii,jj,kk) = PF(ii,jj,zm) ;
398 for( kk=zp+1 ; kk < nz ; kk++ )
399 for( jj=0 ; jj < ny ; jj++ )
400 for( ii=0 ; ii < nx ; ii++ ) PF(ii,jj,kk) = PF(ii,jj,zp) ;
401 }
402
403 if( mask != NULL ) free(mask) ;
404
405 /* write outputs */
406
407 if( prefix != NULL ){
408 THD_3dim_dataset *outset = EDIT_empty_copy( inset ) ;
409 EDIT_dset_items( outset ,
410 ADN_prefix , prefix ,
411 ADN_nvals , 1 ,
412 ADN_ntt , 0 ,
413 ADN_none ) ;
414 EDIT_substitute_brick( outset , 0 , MRI_float , MRI_FLOAT_PTR(imout) ) ;
415 tross_Copy_History( inset , outset ) ;
416 tross_Make_History( "3dPolyfit" , argc,argv , outset ) ;
417 DSET_write(outset) ;
418 WROTE_DSET(outset) ;
419 }
420
421 if( resid != NULL ){
422 THD_3dim_dataset *outset = EDIT_empty_copy( inset ) ;
423 float *inar=MRI_FLOAT_PTR(imin) , *outar=MRI_FLOAT_PTR(imout) ;
424 int nx,ny,nz , nxyz , kk ;
425 nx = imout->nx ; ny = imout->ny ; nz = imout->nz ; nxyz = nx*ny*nz ;
426 for( kk=0 ; kk < nxyz ; kk++ ) outar[kk] = inar[kk] - outar[kk] ;
427 mri_free(imin) ;
428 EDIT_dset_items( outset ,
429 ADN_prefix , resid ,
430 ADN_nvals , 1 ,
431 ADN_ntt , 0 ,
432 ADN_none ) ;
433 EDIT_substitute_brick( outset , 0 , MRI_float , MRI_FLOAT_PTR(imout) ) ;
434 tross_Copy_History( inset , outset ) ;
435 tross_Make_History( "3dPolyfit" , argc,argv , outset ) ;
436 DSET_write(outset) ;
437 WROTE_DSET(outset) ;
438 }
439
440 if( cfnam != NULL && fvit != NULL ){ /* won't work with '-byslice' */
441 char *qn ;
442 qn = STRING_HAS_SUFFIX(cfnam,".1D") ? cfnam : modify_afni_prefix(cfnam,NULL,".1D") ;
443 mri_write_floatvec( qn , fvit ) ;
444 }
445
446 exit(0) ;
447 }
448