1 #include "mrilib.h"
2 #include "mri_rbfinterp.c"
3
4 /**** Adapted from 3dUndump.c ****/
5
Syntax(void)6 void Syntax(void)
7 {
8 printf(
9 "Usage: 3dRBFdset [options]\n"
10 "Uses Radial Basis Functions to fit the values given at a set\n"
11 "of points and then creates a 3D dataset by evaluating the\n"
12 "fit at each point on a 3D grid. Output dataset is stored\n"
13 "in floating point format.\n"
14 "\n"
15 "You might ask 'What use is this program?' The answer is that\n"
16 "it is purely for testing the RBF expansion code, in preparation\n"
17 "for using those functions for something glorious.\n"
18 "\n"
19 "Options:\n"
20 " -prefix ppp = 'ppp' is the prefix for the output dataset\n"
21 " [default = rbf].\n"
22 " -master mmm = 'mmm' is the master dataset, whose geometry\n"
23 " will determine the geometry of the output.\n"
24 " ** This is a MANDATORY option.\n"
25 " -mask MMM = This option specifies a mask dataset 'MMM', which\n"
26 " will control which voxels are allowed to get\n"
27 " values set. If the mask is present, only\n"
28 " voxels that are nonzero in the mask will be\n"
29 " nonzero in the new dataset.\n"
30 " -knotx kfile = 3 column .1D file that lists the 'knot' coordinates\n"
31 " (center of each radial basis function).\n"
32 " ** This is a MANDATORY option.\n"
33 " ** There must be at least 5 knot locations.\n"
34 " -vals vfile = .1D file that lists the values to be fitted at\n"
35 " each fit location.\n"
36 " ** This is a MANDATORY option.\n"
37 " ** One sub-brick will be created for each column\n"
38 " in 'vfile'.\n"
39 " ** 'vfile' must have the same number of values per\n"
40 " column that the '-knotx' file does.\n"
41 " -radius rr = Set the RBF radius to 'rr' mm.\n"
42 " ** Default = compute from knot distribution.\n"
43 " -nolinear = Don't use a global linear polynomial.\n"
44 "\n"
45 "-- RWCox -- February 2009\n"
46 ) ;
47
48 PRINT_COMPILE_DATE ; exit(0) ;
49 }
50
51 /*---------------------------------------------------------------------------*/
52
ijk_to_xyz(THD_3dim_dataset * dset,int ijk,float * x,float * y,float * z)53 static void ijk_to_xyz( THD_3dim_dataset *dset , int ijk ,
54 float *x , float *y , float *z )
55 {
56 THD_fvec3 xyz_mm , xyz_dicom ;
57 THD_ivec3 ivec ;
58
59 ivec.ijk[0] = DSET_index_to_ix(dset,ijk) ;
60 ivec.ijk[1] = DSET_index_to_jy(dset,ijk) ;
61 ivec.ijk[2] = DSET_index_to_kz(dset,ijk) ;
62 xyz_mm = THD_3dind_to_3dmm ( dset , ivec ) ;
63 xyz_dicom = THD_3dmm_to_dicomm( dset , xyz_mm ) ;
64 UNLOAD_FVEC3( xyz_dicom , *x , *y , *z ) ;
65 return ;
66 }
67
68 /*---------------------------------------------------------------------------*/
69
main(int argc,char * argv[])70 int main( int argc , char * argv[] )
71 {
72 char *prefix="rbf" ;
73 MRI_IMAGE *kfim=NULL , *vfim=NULL ;
74 float *kfar=NULL , *vfar=NULL ;
75 int kfnx , vfnx,vfny , uselin=1 ;
76 float rad=0.0f ;
77 byte *mmask=NULL ;
78 THD_3dim_dataset *dset , *maskset=NULL , *mset=NULL ;
79 int iarg , ii,ijk , nx,ny,nz , nxyz,ival,nmask ;
80 float xx,yy,zz , *vv , *dar ;
81 RBF_knots *rbk ;
82 RBF_evalues *rbe ;
83 RBF_evalgrid *rbg ;
84
85 /*------------ help? ------------*/
86
87 if( argc < 3 || strcmp(argv[1],"-help") == 0 ) Syntax() ;
88
89 mainENTRY("3dRBFdset main") ; machdep() ; PRINT_VERSION("3dRBFdset");
90 machdep() ; AFNI_logger("3dRBFdset",argc,argv) ;
91
92 /*------------ command line options ------------*/
93
94 iarg = 1 ;
95 while( iarg < argc && argv[iarg][0] == '-' ){
96
97 /*-----*/
98
99 if( strncmp(argv[iarg],"-rad",4) == 0 ){
100 if( iarg+1 >= argc )
101 ERROR_exit("%s: no argument follows!?",argv[iarg]) ;
102 rad = (float)strtod(argv[++iarg],NULL) ;
103 iarg++ ; continue ;
104 }
105
106 /*-----*/
107
108 if( strncmp(argv[iarg],"-nolinear",5) == 0 ){
109 uselin = 0 ; iarg++ ; continue ;
110 }
111
112 /*-----*/
113
114 if( strncmp(argv[iarg],"-knotx",5) == 0 ){
115 if( kfim != NULL )
116 ERROR_exit("Can't have 2 %s options!",argv[iarg]) ;
117 if( iarg+1 >= argc )
118 ERROR_exit("%s: no argument follows!?",argv[iarg]) ;
119 kfim = mri_read_1D( argv[++iarg] ) ;
120 if( kfim == NULL )
121 ERROR_exit("%s %s: can't read file!?",argv[iarg-1],argv[iarg]) ;
122 if( kfim->ny != 3 )
123 ERROR_exit("%s %s: file must have exactly 3 columns!",argv[iarg-1],argv[iarg]) ;
124 if( kfim->nx < 5 )
125 ERROR_exit("%s %s: file must have 5 or more rows!",argv[iarg-1],argv[iarg]) ;
126 kfnx = kfim->nx ; kfar = MRI_FLOAT_PTR(kfim) ;
127 iarg++ ; continue ;
128 }
129
130 /*-----*/
131
132 if( strncmp(argv[iarg],"-vals",4) == 0 ){
133 if( vfim != NULL )
134 ERROR_exit("Can't have 2 %s options!",argv[iarg]) ;
135 if( iarg+1 >= argc )
136 ERROR_exit("%s: no argument follows!?",argv[iarg]) ;
137 vfim = mri_read_1D( argv[++iarg] ) ;
138 if( vfim == NULL )
139 ERROR_exit("%s %s: can't read file!?",argv[iarg-1],argv[iarg]) ;
140 vfnx = vfim->nx ; vfny = vfim->ny ; vfar = MRI_FLOAT_PTR(vfim) ;
141 iarg++ ; continue ;
142 }
143
144 /*-----*/
145
146 if( strcmp(argv[iarg],"-prefix") == 0 ){
147 if( iarg+1 >= argc )
148 ERROR_exit("-prefix: no argument follows!?") ;
149 else if( !THD_filename_ok(argv[++iarg]) )
150 ERROR_exit("-prefix: Illegal prefix given!") ;
151 prefix = argv[iarg] ;
152 iarg++ ; continue ;
153 }
154
155 /*-----*/
156
157 if( strcmp(argv[iarg],"-master") == 0 ){
158 if( iarg+1 >= argc )
159 ERROR_exit("-master: no argument follows!?") ;
160 else if( mset != NULL )
161 ERROR_exit("-master: can't have two -master options!") ;
162 mset = THD_open_dataset( argv[++iarg] ) ;
163 if( mset == NULL )
164 ERROR_exit("-master: can't open dataset" ) ;
165 iarg++ ; continue ;
166 }
167
168 /*-----*/
169
170 if( strcmp(argv[iarg],"-mask") == 0 ){
171 if( iarg+1 >= argc )
172 ERROR_exit("-mask: no argument follows!?") ;
173 else if( maskset != NULL )
174 ERROR_exit("-mask: can't have two -mask options!") ;
175 maskset = THD_open_dataset( argv[++iarg] ) ;
176 if( maskset == NULL )
177 ERROR_exit("-mask: can't open dataset" ) ;
178 iarg++ ; continue ;
179 }
180
181 /*-----*/
182
183 ERROR_exit("Unknown option: %s",argv[iarg]) ;
184
185 } /* end of loop over command line options */
186
187 /*----- check for inconsistencies or criminal negligence -----*/
188
189 if( mset == NULL )
190 ERROR_exit("3dRBFdset: -master is a mandatory option!") ;
191 if( kfar == NULL )
192 ERROR_exit("3dRBFdset: -knotx is a mandatory option!") ;
193 if( vfar == NULL )
194 ERROR_exit("3dRBFdset: -vals is a mandatory option!") ;
195
196 if( kfnx != vfnx )
197 ERROR_exit("-vals file must have same number lines as -knotx file!") ;
198
199 /*------- make empty dataset -------*/
200
201 dset = EDIT_empty_copy( mset ) ;
202 EDIT_dset_items( dset ,
203 ADN_prefix , prefix ,
204 ADN_datum_all , MRI_float ,
205 ADN_nvals , vfny ,
206 ADN_ntt , 0 ,
207 ADN_func_type , ISANAT(mset) ? ANAT_BUCK_TYPE
208 : FUNC_BUCK_TYPE ,
209 ADN_directory_name , "./" ,
210 ADN_none ) ;
211
212 if( THD_deathcon() && THD_is_file(DSET_HEADNAME(dset)) )
213 ERROR_exit("Output dataset already exists -- can't overwrite") ;
214
215 nx = DSET_NX(dset); ny = DSET_NY(dset); nz = DSET_NZ(dset); nxyz = nx*ny*nz;
216
217 /*--- check and make mask if desired ---*/
218
219 if( maskset != NULL ){
220 if( DSET_NX(maskset) != nx || DSET_NY(maskset) != ny || DSET_NZ(maskset) != nz )
221 ERROR_exit("-mask dataset doesn't match dimension of output dataset") ;
222 mmask = THD_makemask( maskset , 0 , 1.0,-1.0 ) ;
223 DSET_delete(maskset) ;
224 if( mmask == NULL ){
225 WARNING_message("Can't create mask for some unknowable reason!") ;
226 } else {
227 nmask = THD_countmask( nxyz , mmask ) ;
228 if( nmask == 0 ){
229 WARNING_message("0 voxels in mask -- ignoring it!") ;
230 free((void *)mmask) ; mmask = NULL ;
231 } else {
232 INFO_message("%d voxels found in mask",nmask) ;
233 }
234 }
235 }
236 if( mmask == NULL ){
237 mmask = (byte *)malloc(sizeof(byte)*nxyz) ;
238 memset( mmask , 1 , sizeof(byte)*nxyz ) ;
239 nmask = nxyz ;
240 }
241
242 /*----- setup for RBF interpolation at the given set of points -----*/
243
244 RBF_set_verbosity( 2 ) ;
245
246 rbk = RBF_setup_knots( kfnx , rad,uselin ,
247 kfar+0*kfnx , kfar+1*kfnx , kfar+2*kfnx ) ;
248 if( rbk == NULL )
249 ERROR_exit("Can't setup RBF interpolation for some reason!") ;
250
251 /*-- setup output grid points --*/
252
253 MAKE_RBF_evalgrid(rbg,nmask) ;
254 for( ii=ijk=0 ; ijk < nxyz ; ijk++ ){
255 if( mmask[ijk] ){
256 ijk_to_xyz( dset , ijk , &xx,&yy,&zz ) ;
257 rbg->xpt[ii] = xx ;
258 rbg->ypt[ii] = yy ;
259 rbg->zpt[ii] = zz ; ii++ ;
260 }
261 }
262 RBF_setup_kranges( rbk , rbg ) ; /* 20 Feb 2009 */
263
264 /*-- setup space for evaluation --*/
265
266 MAKE_RBF_evalues(rbe,kfnx) ;
267 vv = (float *)malloc(sizeof(float)*nmask) ;
268
269 /*-- loop over columns in vfar and compute results --*/
270
271 for( ival=0 ; ival < vfny ; ival++ ){
272 memcpy( rbe->val , vfar+ival*kfnx , sizeof(float)*kfnx ) ;
273 rbe->code = 0 ; ii = RBF_setup_evalues( rbk , rbe ) ;
274 if( ii == 0 ) ERROR_exit("Can't compute knot coefficients!?") ;
275 ii = RBF_evaluate( rbk , rbe , rbg , vv ) ;
276 if( ii == 0 ) ERROR_exit("Can't compute RBF expansion!?") ;
277 EDIT_substitute_brick( dset , ival , MRI_float , NULL ) ;
278 dar = DSET_BRICK_ARRAY(dset,ival) ;
279 for( ii=ijk=0 ; ijk < nxyz ; ijk++ ){
280 if( mmask[ijk] ) dar[ijk] = vv[ii++] ;
281 }
282 }
283
284 free(vv) ; free(mmask) ;
285 DESTROY_RBF_evalgrid(rbg); DESTROY_RBF_evalues(rbe); DESTROY_RBF_knots(rbk);
286 mri_free(vfim) ; mri_free(kfim) ;
287
288 dset_floatscan(dset) ;
289 tross_Make_History( "3dRBFdset" , argc,argv , dset ) ;
290 DSET_write(dset) ;
291 WROTE_DSET(dset) ;
292 exit(0) ;
293 }
294