00001 #include "mrilib.h"
00002
00003
00004
00005 typedef struct {
00006 int npts ;
00007 int * nstart ;
00008 int * nwt ;
00009 float ** wt ;
00010 } intermap ;
00011
00012
00013
00014
00015
00016
00017 intermap * INTERP_setup_linear( int nin, float * xin, int nout, float * xout )
00018 {
00019 intermap * imap ;
00020 int ii,ib ;
00021 float xx ;
00022
00023 if( nin < 2 || nout < 2 || xin == NULL || xout == NULL ) return NULL ;
00024
00025 imap = (intermap *) malloc(sizeof(intermap)) ;
00026
00027 imap->npts = nout ;
00028 imap->nstart = (int *) calloc(sizeof(int) ,nout) ;
00029 imap->nwt = (int *) calloc(sizeof(int) ,nout) ;
00030 imap->wt = (float **) calloc(sizeof(float *),nout) ;
00031
00032 for( ii=0 ; ii < nout ; ii++ ){
00033 xx = xout[ii] ;
00034
00035 if( xx < xin[0] ){
00036 if( xin[0]-xx <= 0.5*(xin[1]-xin[0]) ){
00037 imap->nstart[ii] = 0 ;
00038 imap->nwt[ii] = 1 ;
00039 imap->wt[ii] = (float *) malloc(sizeof(float)) ;
00040 imap->wt[ii][0] = 1.0 ;
00041 }
00042 continue ;
00043 } else if ( xx > xin[nin-1] ){
00044 if( xx-xin[nin-1] <= 0.5*(xin[nin-1]-xin[nin-2]) ){
00045 imap->nstart[ii] = nin-1 ;
00046 imap->nwt[ii] = 1 ;
00047 imap->wt[ii] = (float *) malloc(sizeof(float)) ;
00048 imap->wt[ii][0] = 1.0 ;
00049 }
00050 continue ;
00051 }
00052
00053
00054
00055
00056 for( ib=0 ; ib < nin-1 ; ib++ )
00057 if( xx >= xin[ib] && xx <= xin[ib+1] ) break ;
00058 if( ib == nin ) continue ;
00059
00060
00061
00062 imap->nstart[ii] = ib ;
00063 imap->nwt[ii] = 2 ;
00064 imap->wt[ii] = (float *) malloc(sizeof(float)*2) ;
00065 imap->wt[ii][1] = (xx-xin[ib]) / (xin[ib+1]-xin[ib]) ;
00066 imap->wt[ii][0] = 1.0 - imap->wt[ii][1] ;
00067 }
00068
00069 return imap ;
00070 }
00071
00072
00073
00074 void INTERP_destroy( intermap * imap )
00075 {
00076 int ii ;
00077
00078 if( imap == NULL ) return ;
00079
00080 for( ii=0 ; ii < imap->npts ; ii++ )
00081 if( imap->wt[ii] != NULL ) free(imap->wt[ii]) ;
00082
00083 free(imap->wt) ; free(imap->nwt) ; free(imap->nstart) ;
00084 free(imap) ; return ;
00085 }
00086
00087
00088
00089
00090
00091
00092 void INTERP_evaluate( intermap * imap , float * vin , float * vout )
00093 {
00094 int ii , jj,ib,nj , nout ;
00095
00096 if( imap == NULL || vin == NULL || vout == NULL ) return ;
00097
00098 nout = imap->npts ;
00099 for( ii=0 ; ii < nout ; ii++ ){
00100 vout[ii] = 0.0 ;
00101 ib = imap->nstart[ii] ;
00102 nj = imap->nwt[ii] ;
00103 for( jj=0 ; jj < nj ; jj++ )
00104 vout[ii] += imap->wt[ii][jj] * vin[ib+jj] ;
00105 }
00106 return ;
00107 }
00108
00109
00110
00111 int main( int argc , char * argv[] )
00112 {
00113 THD_3dim_dataset * dset , * nset ;
00114 int iarg=1 ;
00115 float new_dz=0.0 , old_dz ;
00116 int new_nz=0 , old_nz ;
00117 float new_zsize=0.0 , old_zsize ;
00118 char * prefix = "regrid" ;
00119 float old_zcen , new_zcen , zshift ;
00120 int bkind , nx,ny,nxy , iv,ii,kk , verb=0 ;
00121 byte *bar , *nbar ;
00122 short *sar , *nsar ;
00123 float *far , *nfar , *ofz , *nfz ;
00124 complex *car , *ncar ;
00125 void * nvar ;
00126 float * zin , * zout ;
00127 intermap * imap ;
00128
00129 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00130 printf("Usage: 3dZregrid [option] dataset\n"
00131 "Alters the input dataset's slice thickness and/or number.\n"
00132 "\n"
00133 "OPTIONS:\n"
00134 " -dz D = sets slice thickness to D mm\n"
00135 " -nz N = sets slice count to N\n"
00136 " -zsize Z = sets thickness of dataset (center-to-center of\n"
00137 " first and last slices) to Z mm\n"
00138 " -prefix P = write result in dataset with prefix P\n"
00139 " -verb = write progress reports to stderr\n"
00140 "\n"
00141 "At least one of '-dz', '-nz', or '-zsize' must be given.\n"
00142 "On the other hand, using all 3 is over-specification.\n"
00143 "The following combinations make sense:\n"
00144 " -dz only ==> N stays fixed from input dataset\n"
00145 " and then is like setting Z = N*D\n"
00146 " -dz and -nz together ==> like setting Z = N*D\n"
00147 " -dz and -zsize together ==> like setting N = Z/D\n"
00148 " -nz only ==> D stays fixed from input dataset\n"
00149 " and then is like setting Z = N*D\n"
00150 " -zsize only ==> D stays fixed from input dataset\n"
00151 " and then is like setting N = Z/D\n"
00152 " -nsize and -zsize together ==> like setting D = Z/N\n"
00153 "\n"
00154 "NOTES:\n"
00155 " * If the input is a 3D+time dataset with slice-dependent time\n"
00156 " offsets, the output will have its time offsets cleared.\n"
00157 " It probably makes sense to do 3dTshift BEFORE using this\n"
00158 " program in such a case.\n"
00159 " * The output of this program is centered around the same\n"
00160 " location as the input dataset. Slices outside the\n"
00161 " original volume (e.g., when Z is increased) will be\n"
00162 " zero. This is NOT the same as using 3dZeropad, which\n"
00163 " only adds zeros, and does not interpolate to a new grid.\n"
00164 " * Linear interpolation is used between slices. However,\n"
00165 " new slice positions outside the old volume but within\n"
00166 " 0.5 old slice thicknesses will get a copy of the last slice.\n"
00167 " New slices outside this buffer zone will be all zeros.\n"
00168 "\n"
00169 "EXAMPLE:\n"
00170 " You have two 3D anatomical datasets from the same subject that\n"
00171 " need to be registered. Unfortunately, the first one has slice\n"
00172 " thickness 1.2 mm and the second 1.3 mm. Assuming they have\n"
00173 " the same number of slices, then do something like\n"
00174 " 3dZregrid -dz 1.2 -prefix ElvisZZ Elvis2+orig\n"
00175 " 3dvolreg -base Elvis1+orig -prefix Elvis2reg ElvisZZ+orig\n"
00176 ) ;
00177 exit(0) ;
00178 }
00179
00180 mainENTRY("3dZregrid main"); machdep(); PRINT_VERSION("3dZregrid") ;
00181 AFNI_logger("3dZregrid",argc,argv) ;
00182
00183
00184
00185 while( iarg < argc && argv[iarg][0] == '-' ){
00186
00187 if( strncmp(argv[iarg],"-verb",5) == 0 ){
00188 verb++ ; iarg++ ; continue ;
00189 }
00190
00191 if( strcmp(argv[iarg],"-dz") == 0 ){
00192 new_dz = strtod( argv[++iarg] , NULL ) ;
00193 if( new_dz <= 0.0 ){ fprintf(stderr,"** ILLEGAL -dz value!\n"); exit(1); }
00194 iarg++ ; continue ;
00195 }
00196
00197 if( strcmp(argv[iarg],"-nz") == 0 ){
00198 new_nz = strtod( argv[++iarg] , NULL ) ;
00199 if( new_nz <= 2 ){ fprintf(stderr,"** ILLEGAL -nz value!\n"); exit(1); }
00200 iarg++ ; continue ;
00201 }
00202
00203 if( strcmp(argv[iarg],"-zsize") == 0 ){
00204 new_zsize = strtod( argv[++iarg] , NULL ) ;
00205 if( new_zsize <= 0.0 ){ fprintf(stderr,"** ILLEGAL -zsize value!\n"); exit(1); }
00206 iarg++ ; continue ;
00207 }
00208
00209 if( strcmp(argv[iarg],"-prefix") == 0 ){
00210 prefix = argv[++iarg] ;
00211 if( !THD_filename_ok(prefix) ){ fprintf(stderr,"** ILLEGAL -prefix value!\n"); exit(1); }
00212 iarg++ ; continue ;
00213 }
00214
00215 fprintf(stderr,"** ILLEGAL option: %s\n",argv[iarg]) ; exit(1) ;
00216 }
00217
00218
00219
00220 if( verb ) fprintf(stderr,"++ Reading input dataset %s\n",argv[iarg]) ;
00221 dset = THD_open_dataset( argv[iarg] ) ;
00222 if( dset == NULL ){
00223 fprintf(stderr,"** CANNOT open dataset %s\n",argv[iarg]) ; exit(1) ;
00224 }
00225 DSET_load(dset) ;
00226 if( !DSET_LOADED(dset) ){
00227 fprintf(stderr,"** CANNOT load dataset %s\n",argv[iarg]) ; exit(1) ;
00228 }
00229
00230
00231
00232 if( new_dz == 0.0 && new_nz == 0 && new_zsize == 0.0 ){
00233 fprintf(stderr,"** AT LEAST one of -dz, -nz, -zsize must be given!\n"); exit(1);
00234 }
00235
00236
00237
00238 if( new_dz > 0.0 ){
00239 if( new_nz > 0 ){
00240 new_zsize = new_dz * new_nz ;
00241 } else if( new_zsize > 0.0 ){
00242 new_nz = rint( new_zsize / new_dz ) ;
00243 } else {
00244 new_nz = DSET_NZ(dset) ;
00245 new_zsize = new_dz * new_nz ;
00246 }
00247 } else if( new_nz > 0 ){
00248
00249 if( new_zsize > 0.0 ){
00250 new_dz = new_zsize / new_nz ;
00251 } else {
00252 new_dz = fabs(DSET_DZ(dset)) ;
00253 new_zsize = new_nz * new_dz ;
00254 }
00255 } else {
00256 new_dz = fabs(DSET_DZ(dset)) ;
00257 new_nz = rint( new_zsize / new_dz ) ;
00258 }
00259
00260
00261
00262 if( new_nz < 2 ){
00263 fprintf(stderr,"** ILLEGAL number of new slices computed = %d\n",new_nz) ;
00264 exit(1) ;
00265 }
00266
00267
00268
00269 old_nz = DSET_NZ(dset) ;
00270 old_dz = fabs(DSET_DZ(dset)) ;
00271 old_zsize = old_nz * old_dz ;
00272
00273 nset = EDIT_empty_copy( dset ) ;
00274 EDIT_dset_items( nset , ADN_prefix , prefix , ADN_none ) ;
00275 if( THD_is_file( DSET_HEADNAME(nset) ) ){
00276 fprintf(stderr,"** Output file %s exists -- cannot overwrite!\n",
00277 DSET_HEADNAME(nset) ) ;
00278 exit(1) ;
00279 }
00280
00281 tross_Copy_History( dset , nset ) ;
00282 tross_Make_History( "3dZregrid" , argc,argv , nset ) ;
00283
00284
00285
00286 if( new_nz != old_nz ){
00287 THD_ivec3 iv_nxyz ;
00288 LOAD_IVEC3( iv_nxyz , DSET_NX(dset) , DSET_NY(dset) , new_nz ) ;
00289 EDIT_dset_items( nset , ADN_nxyz , iv_nxyz , ADN_none ) ;
00290 if( verb )
00291 fprintf(stderr,"++ Adjusting slice count from %d to %d\n",old_nz,new_nz) ;
00292 }
00293
00294 if( new_dz != old_dz ){
00295 THD_fvec3 fv_dxyz ; float dz=new_dz ;
00296 if( DSET_DZ(dset) < 0.0 ) dz = -dz ;
00297 LOAD_FVEC3( fv_dxyz , DSET_DX(dset) , DSET_DY(dset) , dz ) ;
00298 EDIT_dset_items( nset , ADN_xyzdel , fv_dxyz , ADN_none ) ;
00299 if( verb )
00300 fprintf(stderr,"++ Adjusting slice thickness from %6.3f to %6.3f\n",old_dz,new_dz) ;
00301 }
00302
00303 old_zcen = dset->daxes->zzorg + 0.5*(dset->daxes->nzz-1)*dset->daxes->zzdel ;
00304 new_zcen = nset->daxes->zzorg + 0.5*(nset->daxes->nzz-1)*nset->daxes->zzdel ;
00305 zshift = new_zcen - old_zcen ;
00306 if( fabs(zshift) > 0.01*new_dz ){
00307 THD_fvec3 fv_xyzorg ; float zz ;
00308 zz = nset->daxes->zzorg - zshift ;
00309 LOAD_FVEC3( fv_xyzorg , nset->daxes->xxorg , nset->daxes->yyorg , zz ) ;
00310 EDIT_dset_items( nset , ADN_xyzorg , fv_xyzorg , ADN_none ) ;
00311 }
00312
00313 if( nset->taxis != NULL && nset->taxis->nsl > 0 ){
00314 EDIT_dset_items( nset , ADN_nsl , 0 , ADN_none ) ;
00315 fprintf(stderr,
00316 "++ WARNING - slice-dependent time shifts have been removed!\n") ;
00317 }
00318
00319
00320
00321 zin = (float *) malloc(sizeof(float)*old_nz) ;
00322 zout = (float *) malloc(sizeof(float)*new_nz) ;
00323 for( kk=0 ; kk < old_nz ; kk++ ) zin[kk] = (kk-0.5*old_nz) * old_dz ;
00324 for( kk=0 ; kk < new_nz ; kk++ ) zout[kk] = (kk-0.5*new_nz) * new_dz ;
00325
00326 imap =INTERP_setup_linear( old_nz,zin , new_nz,zout ) ;
00327
00328 free(zin) ; free(zout) ;
00329
00330
00331
00332 nx = DSET_NX(nset) ; ny = DSET_NY(nset) ; nxy = nx*ny ;
00333
00334 ofz = (float *) malloc(sizeof(float)*old_nz*2) ;
00335 nfz = (float *) malloc(sizeof(float)*new_nz*2) ;
00336
00337 if( verb ) fprintf(stderr,"++ Sub-bricks:") ;
00338
00339 for( iv=0 ; iv < DSET_NVALS(nset) ; iv++ ){
00340
00341 if( verb ) fprintf(stderr,"%d",iv) ;
00342
00343 bkind = DSET_BRICK_TYPE(dset,iv) ;
00344 switch( bkind ){
00345 default:
00346 fprintf(stderr,"** ILLEGAL brick type %s at sub-brick %d\n",
00347 MRI_type_name[bkind] , iv ) ;
00348 exit(1) ;
00349
00350 case MRI_byte: bar = DSET_ARRAY(dset,iv) ;
00351 nbar = (byte *)malloc(sizeof(byte)*nxy*new_nz) ;
00352 nvar = (void *) nbar ;
00353 break ;
00354
00355 case MRI_short: sar = DSET_ARRAY(dset,iv) ;
00356 nsar = (short *)malloc(sizeof(short)*nxy*new_nz) ;
00357 nvar = (void *) nsar ;
00358 break ;
00359
00360 case MRI_float: far = DSET_ARRAY(dset,iv) ;
00361 nfar = (float *)malloc(sizeof(float)*nxy*new_nz) ;
00362 nvar = (void *) nfar ;
00363 break ;
00364
00365 case MRI_complex: car = DSET_ARRAY(dset,iv) ;
00366 ncar = (complex *)malloc(sizeof(complex)*nxy*new_nz) ;
00367 nvar = (void *) ncar ;
00368 break ;
00369 }
00370
00371
00372
00373 for( ii=0 ; ii < nxy ; ii++ ){
00374
00375
00376
00377 switch( bkind ){
00378 case MRI_byte:
00379 for( kk=0 ; kk < old_nz ; kk++ ) ofz[kk] = bar[ii+kk*nxy] ;
00380 break ;
00381
00382 case MRI_short:
00383 for( kk=0 ; kk < old_nz ; kk++ ) ofz[kk] = sar[ii+kk*nxy] ;
00384 break ;
00385
00386 case MRI_float:
00387 for( kk=0 ; kk < old_nz ; kk++ ) ofz[kk] = far[ii+kk*nxy] ;
00388 break ;
00389
00390 case MRI_complex:
00391 for( kk=0 ; kk < old_nz ; kk++ ){
00392 ofz[kk] = car[ii+kk*nxy].r ;
00393 ofz[kk+old_nz] = car[ii+kk*nxy].i ;
00394 }
00395 break ;
00396 }
00397
00398
00399
00400 INTERP_evaluate( imap , ofz , nfz ) ;
00401 if( bkind == MRI_complex )
00402 INTERP_evaluate( imap , ofz+old_nz , nfz+new_nz ) ;
00403
00404
00405
00406 switch( bkind ){
00407 case MRI_byte:
00408 for( kk=0 ; kk < new_nz ; kk++ ) nbar[ii+kk*nxy] = BYTEIZE(nfz[kk]) ;
00409 break ;
00410
00411 case MRI_short:
00412 for( kk=0 ; kk < new_nz ; kk++ ) nsar[ii+kk*nxy] = SHORTIZE(nfz[kk]) ;
00413 break ;
00414
00415 case MRI_float:
00416 for( kk=0 ; kk < new_nz ; kk++ ) nfar[ii+kk*nxy] = nfz[kk] ;
00417 break ;
00418
00419 case MRI_complex:
00420 for( kk=0 ; kk < new_nz ; kk++ ){
00421 ncar[ii+kk*nxy].r = nfz[kk] ;
00422 ncar[ii+kk*nxy].i = nfz[kk+new_nz] ;
00423 }
00424 break ;
00425 }
00426 }
00427
00428 if( verb ) fprintf(stderr,".") ;
00429
00430 EDIT_substitute_brick( nset , iv , bkind , nvar ) ;
00431 DSET_unload_one( dset , iv ) ;
00432
00433 }
00434
00435 DSET_delete(dset) ; INTERP_destroy(imap) ; free(nfz) ; free(ofz) ;
00436
00437 DSET_write(nset) ;
00438 if( verb ) fprintf(stderr,"\n++ Output dataset = %s\n",DSET_BRIKNAME(nset)) ;
00439 exit(0) ;
00440 }
