1 /***** This code is part of the AFNI software package, which is *****
2 ***** partly in the public domain and partly covered by the GPL. *****
3 ***** See https://afni.nimh.nih.gov/afni for more information. *****/
4
5 #include "mrilib.h"
6
7 static void vstep_print(void) ; /* prototype */
8
main(int argc,char * argv[])9 int main( int argc , char *argv[] )
10 {
11 int vstep=0 , ii,nvox , ntin , ntout , do_one=0 , nup=-1 ;
12 THD_3dim_dataset *inset=NULL , *outset ;
13 char *prefix="Upsam", *dsetname=NULL ;
14 int verb=0 , iarg=1, datum = MRI_float;
15 float *ivec , *ovec , trin , trout, *fac=NULL, *ofac=NULL,
16 top=0.0, maxtop=0.0;
17
18 /*------- help the pitifully ignorant user? -------*/
19
20 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
21 printf(
22 "Usage: 3dUpsample [options] n dataset\n"
23 "\n"
24 "* Upsamples a 3D+time dataset, in the time direction,\n"
25 " by a factor of 'n'.\n"
26 "* The value of 'n' must be between 2 and 320 (inclusive).\n"
27 "* The output dataset is in float format by default.\n"
28 "\n"
29 "Options:\n"
30 "--------\n"
31 " -1 or -one = Use linear interpolation. Otherwise,\n"
32 " or -linear 7th order polynomial interpolation is used.\n"
33 "\n"
34 " -prefix pp = Define the prefix name of the output dataset.\n"
35 " [default prefix is 'Upsam']\n"
36 "\n"
37 " -verb = Be eloquently and mellifluosly verbose.\n"
38 "\n"
39 " -n n = An alternate way to specify n\n"
40 " -input dataset = An alternate way to specify dataset\n"
41 "\n"
42 " -datum ddd = Use datatype ddd at output. Choose from\n"
43 " float (default), short, byte.\n"
44 "Example:\n"
45 "--------\n"
46 " 3dUpsample -prefix LongFred 5 Fred+orig\n"
47 "\n"
48 "Nota Bene:\n"
49 "----------\n"
50 "* You should not use this for files that were 3dTcat-ed across\n"
51 " imaging run boundaries, since that will result in interpolating\n"
52 " between non-contiguous time samples!\n"
53 "* If the input has M time points, the output will have n*M time\n"
54 " points. The last n-1 of them will be past the end of the original\n"
55 " time series.\n"
56 "* This program gobbles up memory and diskspace as a function of n.\n"
57 " You can reduce output file size with -datum option.\n"
58 "\n"
59 "--- RW Cox - April 2008\n"
60 ) ;
61 PRINT_COMPILE_DATE ; exit(0) ;
62 }
63
64 mainENTRY("3dUpsample"); machdep();
65 PRINT_VERSION("3dUpsample"); AUTHOR("RWCox") ;
66 AFNI_logger("3dUpsample",argc,argv);
67
68 /*------- read command line args -------*/
69
70 datum = MRI_float;
71 iarg = 1 ;
72 while( iarg < argc && argv[iarg][0] == '-' ){
73
74 if( strncasecmp(argv[iarg],"-prefix",5) == 0 ){
75 if( ++iarg >= argc )
76 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
77 prefix = argv[iarg] ;
78 if( !THD_filename_ok(prefix) )
79 ERROR_exit("Illegal string after -prefix: '%s'",prefix) ;
80 iarg++ ; continue ;
81 }
82
83 if( strncasecmp(argv[iarg],"-one",4) == 0 ||
84 strcmp (argv[iarg],"-1" ) == 0 ||
85 strncasecmp(argv[iarg],"-lin",4) == 0 ){
86 do_one = 1 ; iarg++ ; continue ;
87 }
88
89 if( strncasecmp(argv[iarg],"-verb",3) == 0 ){
90 verb = 1 ; iarg++ ; continue ;
91 }
92
93 if( strcasecmp(argv[iarg],"-n") == 0 ){
94 if( ++iarg >= argc )
95 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
96 nup = (int)strtod(argv[iarg],NULL) ;
97 if( nup < 2 || nup > 320 )
98 ERROR_exit("3dUpsample rate '%d' is outside range 2..320",nup) ;
99 iarg++ ; continue ;
100 }
101
102 if( strcasecmp(argv[iarg],"-input") == 0 ){
103 if( ++iarg >= argc )
104 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
105 dsetname = argv[iarg];
106 iarg++ ; continue ;
107 }
108
109 if( strcasecmp(argv[iarg],"-datum") == 0 ){
110 if( ++iarg >= argc )
111 ERROR_exit("Need argument after '%s'",argv[iarg-1]);
112
113 if( strcmp(argv[iarg],"short") == 0 ){
114 datum = MRI_short ;
115 } else if( strcmp(argv[iarg],"float") == 0 ){
116 datum = MRI_float ;
117 } else if( strcmp(argv[iarg],"byte") == 0 ){
118 datum = MRI_byte ;
119 } else {
120 ERROR_message("-datum of type '%s' not supported in 3dUpsample!\n",
121 argv[iarg] ) ;
122 exit(1) ;
123 }
124
125 iarg++ ; continue ;
126 }
127
128 ERROR_message("Unknown argument on command line: '%s'",argv[iarg]) ;
129 suggest_best_prog_option(argv[0], argv[iarg]);
130 exit (1);
131 }
132
133 /*------- check options for completeness and consistency -----*/
134
135 if (nup == -1) {
136 if( iarg+1 >= argc )
137 ERROR_exit("need 'n' and 'dataset' on command line!") ;
138
139 nup = (int)strtod(argv[iarg++],NULL) ;
140 if( nup < 2 || nup > 320 )
141 ERROR_exit("3dUpsample rate '%d' is outside range 2..320",nup) ;
142 }
143 if (!dsetname) {
144 if( iarg >= argc )
145 ERROR_exit("need 'dataset' on command line!") ;
146 dsetname = argv[iarg];
147 }
148
149 inset = THD_open_dataset(dsetname) ;
150 if( !ISVALID_DSET(inset) )
151 ERROR_exit("3dUpsample can't open dataset '%s'", dsetname) ;
152 ntin = DSET_NVALS(inset) ; trin = DSET_TR(inset) ;
153 if( ntin < 2 )
154 ERROR_exit("dataset '%s' has only 1 value per voxel?!",dsetname) ;
155
156 nvox = DSET_NVOX(inset) ;
157
158 if( verb ) INFO_message("loading input dataset into memory") ;
159
160 DSET_load(inset) ; CHECK_LOAD_ERROR(inset) ;
161
162
163 /*------ create output dataset ------*/
164
165 ntout = ntin * nup ; trout = trin / nup ;
166
167 /* scaling factor for output */
168 fac = NULL; maxtop = 0.0;
169 if (MRI_IS_INT_TYPE(datum)) {
170 fac = (float *)calloc(DSET_NVALS(inset), sizeof(float));
171 ofac = (float *)calloc(ntout, sizeof(float));
172 for (ii=0; ii<DSET_NVALS(inset); ++ii) {
173 top = MCW_vol_amax( DSET_NVOX(inset),1,1 ,
174 DSET_BRICK_TYPE(inset,ii),
175 DSET_BRICK_ARRAY(inset,ii) ) ;
176 if (DSET_BRICK_FACTOR(inset, ii))
177 top = top * DSET_BRICK_FACTOR(inset,ii);
178 fac[ii] = (top > MRI_TYPE_maxval[datum]) ?
179 top/MRI_TYPE_maxval[datum] : 0.0 ;
180 if (top > maxtop) maxtop = top;
181 }
182 if (storage_mode_from_filename(prefix) != STORAGE_BY_BRICK) {
183 fac[0] = (maxtop > MRI_TYPE_maxval[datum]) ?
184 maxtop/MRI_TYPE_maxval[datum] : 0.0 ;
185 for (ii=0; ii<ntout; ++ii)
186 ofac[ii] = fac[0];
187 if (verb) INFO_message("Forcing global scaling, Max = %f, fac = %f\n",
188 maxtop, fac[0]);
189 } else {
190 if (verb) INFO_message("Reusing scaling factors of input dset\n");
191 upsample_1( nup, DSET_NVALS(inset), fac, ofac);
192 }
193 }
194 free(fac); fac = NULL;
195 outset = EDIT_empty_copy(inset) ;
196 EDIT_dset_items( outset ,
197 ADN_nvals , ntout ,
198 ADN_ntt , DSET_NUM_TIMES(inset) > 1 ? ntout : 0 ,
199 ADN_datum_all , datum ,
200 ADN_brick_fac , ofac ,
201 ADN_prefix , prefix ,
202 ADN_none ) ;
203 tross_Copy_History( inset , outset ) ;
204 tross_Make_History( "3dUpsample" , argc,argv , outset ) ;
205 free(ofac); ofac = NULL;
206
207 if( outset->taxis != NULL ){
208 outset->taxis->ttdel /= nup ;
209 outset->taxis->ttdur /= nup ;
210 if( outset->taxis->toff_sl != NULL ){
211 for( ii=0 ; ii < outset->taxis->nsl ; ii++ )
212 outset->taxis->toff_sl[ii] /= nup ;
213 }
214 }
215
216 for( ii=0 ; ii < ntout ; ii++ ){ /* create empty bricks to be filled below */
217 EDIT_substitute_brick( outset , ii , datum , NULL ) ;
218 }
219
220 /*------- loop over voxels and process them one at a time ---------*/
221
222 if( verb )
223 INFO_message("Upsampling time series from %d to %d: %s interpolation",
224 ntin , ntout , (do_one) ? "linear" : "heptic" ) ;
225
226 if( verb && nvox > 499 ) vstep = nvox / 50 ;
227 if( vstep > 0 ) fprintf(stderr,"++ voxel loop: ") ;
228
229 ivec = (float *)malloc(sizeof(float)*ntin) ;
230 ovec = (float *)malloc(sizeof(float)*ntout) ;
231
232 for( ii=0 ; ii < nvox ; ii++ ){
233
234 if( vstep > 0 && ii%vstep==vstep-1 ) vstep_print() ;
235
236 THD_extract_array( ii , inset , 0 , ivec ) ;
237
238 if( do_one ) upsample_1( nup , ntin , ivec , ovec ) ;
239 else upsample_7( nup , ntin , ivec , ovec ) ;
240
241 THD_insert_series( ii , outset , ntout , MRI_float , ovec ,
242 datum==MRI_float ? 1:0 ) ;
243 } /* end of loop over voxels */
244
245 if( vstep > 0 ) fprintf(stderr," Done!\n") ;
246
247 /*----- clean up and go away -----*/
248
249 DSET_write(outset) ;
250 if( verb ) WROTE_DSET(outset) ;
251 if( verb ) INFO_message("Total CPU time = %.1f s",COX_cpu_time()) ;
252 exit(0);
253 }
254
255 /*---------------------------------------------------------------------------*/
256
vstep_print(void)257 static void vstep_print(void)
258 {
259 static int nn=0 ;
260 static char xx[10] = "0123456789" ;
261 fprintf(stderr , "%c" , xx[nn%10] ) ;
262 if( nn%10 == 9) fprintf(stderr,".") ;
263 nn++ ;
264 }
265