00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011 #define ORCODE(aa) \
00012 ( (aa)=='R' ? ORI_R2L_TYPE : (aa)=='L' ? ORI_L2R_TYPE : \
00013 (aa)=='P' ? ORI_P2A_TYPE : (aa)=='A' ? ORI_A2P_TYPE : \
00014 (aa)=='I' ? ORI_I2S_TYPE : (aa)=='S' ? ORI_S2I_TYPE : ILLEGAL_TYPE )
00015
00016 #define OR3OK(x,y,z) ( ((x)&6) + ((y)&6) + ((z)&6) == 6 )
00017
00018
00019
00020 void Syntax(char * msg)
00021 {
00022 if( msg != NULL ){
00023 fprintf(stderr,"*** %s\n",msg) ; exit(1) ;
00024 }
00025
00026 printf(
00027 "Usage: 3dUndump [options] infile ...\n"
00028 "Assembles a 3D dataset from an ASCII list of coordinates and\n"
00029 "(optionally) values.\n"
00030 "\n"
00031 "Options:\n"
00032 " -prefix ppp = 'ppp' is the prefix for the output dataset\n"
00033 " [default = undump].\n"
00034 " -master mmm = 'mmm' is the master dataset, whose geometry\n"
00035 " *OR* will determine the geometry of the output.\n"
00036 " -dimen I J K = Sets the dimensions of the output dataset to\n"
00037 " be I by J by K voxels. (Each I, J, and K\n"
00038 " must be >= 2.) This option can be used to\n"
00039 " create a dataset of a specific size for test\n"
00040 " purposes, when no suitable master exists.\n"
00041 " ** N.B.: Exactly one of -master or -dimen must be given.\n"
00042 " -mask kkk = This option specifies a mask dataset 'kkk', which\n"
00043 " will control which voxels are allowed to get\n"
00044 " values set. If the mask is present, only\n"
00045 " voxels that are nonzero in the mask can be\n"
00046 " set in the new dataset.\n"
00047 " * A mask can be created with program 3dAutomask.\n"
00048 " * Combining a mask with sphere insertion makes\n"
00049 " a lot of sense (to me, at least).\n"
00050 " -datum type = 'type' determines the voxel data type of the\n"
00051 " output, which may be byte, short, or float\n"
00052 " [default = short].\n"
00053 " -dval vvv = 'vvv' is the default value stored in each\n"
00054 " input voxel that does not have a value\n"
00055 " supplied in the input file [default = 1].\n"
00056 " -fval fff = 'fff' is the fill value, used for each voxel\n"
00057 " in the output dataset that is NOT listed\n"
00058 " in the input file [default = 0].\n"
00059 " -ijk = Coordinates in the input file are (i,j,k) index\n"
00060 " *OR* triples, as might be output by 3dmaskdump.\n"
00061 " -xyz = Coordinates in the input file are (x,y,z)\n"
00062 " spatial coordinates, in mm. If neither\n"
00063 " -ijk or -xyz is given, the default is -ijk.\n"
00064 " ** N.B.: -xyz can only be used with -master. If -dimen\n"
00065 " is used to specify the size of the output dataset,\n"
00066 " (x,y,z) coordinates are not defined (until you\n"
00067 " use 3drefit to define the spatial structure).\n"
00068 " -srad rrr = Specifies that a sphere of radius 'rrr' will be\n"
00069 " filled about each input (x,y,z) or (i,j,k) voxel.\n"
00070 " If the radius is not given, or is 0, then each\n"
00071 " input data line sets the value in only one voxel.\n"
00072 " * If '-master' is used, then 'rrr' is in mm.\n"
00073 " * If '-dimen' is used, then 'rrr' is in voxels.\n"
00074 " -orient code = Specifies the coordinate order used by -xyz.\n"
00075 " The code must be 3 letters, one each from the pairs\n"
00076 " {R,L} {A,P} {I,S}. The first letter gives the\n"
00077 " orientation of the x-axis, the second the orientation\n"
00078 " of the y-axis, the third the z-axis:\n"
00079 " R = right-to-left L = left-to-right\n"
00080 " A = anterior-to-posterior P = posterior-to-anterior\n"
00081 " I = inferior-to-superior S = superior-to-inferior\n"
00082 " If -orient isn't used, then the coordinate order of the\n"
00083 " -master dataset is used to interpret (x,y,z) inputs.\n"
00084 " ** N.B.: If -dimen is used (which implies -ijk), then the\n"
00085 " only use of -orient is to specify the axes ordering\n"
00086 " of the output dataset. If -master is used instead,\n"
00087 " the output dataset's axes ordering is the same as the\n"
00088 " -master dataset's, regardless of -orient.\n"
00089 "\n"
00090 "Input File Format:\n"
00091 " The input file(s) are ASCII files, with one voxel specification per\n"
00092 " line. A voxel specification is 3 numbers (-ijk or -xyz coordinates),\n"
00093 " with an optional 4th number giving the voxel value. For example:\n"
00094 "\n"
00095 " 1 2 3 \n"
00096 " 3 2 1 5\n"
00097 " 5.3 6.2 3.7\n"
00098 " // this line illustrates a comment\n"
00099 "\n"
00100 " The first line puts a voxel (with value given by -dval) at point\n"
00101 " (1,2,3). The second line puts a voxel (with value 5) at point (3,2,1).\n"
00102 " The third line puts a voxel (with value given by -dval) at point\n"
00103 " (5.3,6.2,3.7). If -ijk is in effect, and fractional coordinates\n"
00104 " are given, they will be rounded to the nearest integers; for example,\n"
00105 " the third line would be equivalent to (i,j,k) = (5,6,4).\n"
00106 "\n"
00107 "Notes:\n"
00108 "* This program creates a 1 sub-brick file. You can 'glue' multiple\n"
00109 " files together using 3dbucket or 3dTcat to make multi-brick datasets.\n"
00110 "* If an input filename is '-', then stdin is used for input.\n"
00111 "* By default, the output dataset is of type '-fim', unless the -master\n"
00112 " dataset is an anat type. You can change the output type using 3drefit.\n"
00113 "* You could use program 1dcat to extract specific columns from a\n"
00114 " multi-column rectangular file (e.g., to get a specific sub-brick\n"
00115 " from the output of 3dmaskdump), and use the output of 1dcat as input\n"
00116 " to this program.\n"
00117 "* [19 Feb 2004] The -mask and -srad options were added this day.\n"
00118 " Also, a fifth value on an input line, if present, is taken as a\n"
00119 " sphere radius to be used for that input point only. Thus, input\n"
00120 " 3.3 4.4 5.5 6.6 7.7\n"
00121 " means to put the value 6.6 into a sphere of radius 7.7 mm centered\n"
00122 " about (x,y,z)=(3.3,4.4,5.5).\n"
00123 "\n"
00124 "-- RWCox -- October 2000\n"
00125 ) ;
00126
00127 exit(0) ;
00128 }
00129
00130
00131
00132 #define NBUF 1024
00133
00134 int main( int argc , char * argv[] )
00135 {
00136 int do_ijk=1 , dimen_ii=0 , dimen_jj=0 , dimen_kk=0 , datum=MRI_short ;
00137 THD_3dim_dataset *mset=NULL ;
00138 char *prefix="undump" , *orcode=NULL ;
00139 THD_coorder cord ;
00140 float dval_float=1.0 , fval_float=0.0 , *fbr=NULL ;
00141 short dval_short=1 , fval_short=0 , *sbr=NULL ;
00142 byte dval_byte =1 , fval_byte =0 , *bbr=NULL , *mmask=NULL ;
00143
00144 FILE *fp ;
00145 THD_3dim_dataset *dset , *maskset=NULL ;
00146 int iarg , ii,jj,kk,ll,ijk , nx,ny,nz , nxyz , nn ;
00147 float xx,yy,zz,vv=0.0 ;
00148 short sv=0 ;
00149 byte bv=0 ;
00150 char linbuf[NBUF] , *cp ;
00151
00152 float xxdown,xxup , yydown,yyup , zzdown,zzup ;
00153
00154 float srad=0.0 , vrad,rii,rjj,rkk,qii,qjj,qkk , dx,dy,dz ;
00155
00156
00157
00158 if( argc < 3 || strcmp(argv[1],"-help") == 0 ) Syntax(NULL) ;
00159
00160
00161
00162 mainENTRY("3dUndump main") ; machdep() ; PRINT_VERSION("3dUndump");
00163
00164 machdep() ;
00165 { int new_argc ; char ** new_argv ;
00166 addto_args( argc , argv , &new_argc , &new_argv ) ;
00167 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00168 }
00169
00170 AFNI_logger("3dUndump",argc,argv) ;
00171
00172
00173
00174 iarg = 1 ;
00175 while( iarg < argc && argv[iarg][0] == '-' ){
00176
00177 if( strcmp(argv[iarg],"-") == 0 ){
00178 break ;
00179 }
00180
00181
00182
00183 if( strcmp(argv[iarg],"-prefix") == 0 ){
00184 if( iarg+1 >= argc )
00185 Syntax("-prefix: no argument follows!?") ;
00186 else if( !THD_filename_ok(argv[++iarg]) )
00187 Syntax("-prefix: Illegal prefix given!") ;
00188 prefix = argv[iarg] ;
00189 iarg++ ; continue ;
00190 }
00191
00192
00193
00194 if( strcmp(argv[iarg],"-master") == 0 ){
00195 if( iarg+1 >= argc )
00196 Syntax("-master: no argument follows!?") ;
00197 else if( mset != NULL )
00198 Syntax("-master: can't have two -master options!") ;
00199 else if( dimen_ii > 0 )
00200 Syntax("-master: conflicts with previous -dimen!") ;
00201
00202 mset = THD_open_dataset( argv[++iarg] ) ;
00203 if( mset == NULL )
00204 Syntax("-master: can't open dataset" ) ;
00205
00206 iarg++ ; continue ;
00207 }
00208
00209
00210
00211 if( strcmp(argv[iarg],"-mask") == 0 ){
00212 if( iarg+1 >= argc )
00213 Syntax("-mask: no argument follows!?") ;
00214 else if( maskset != NULL )
00215 Syntax("-mask: can't have two -mask options!") ;
00216
00217 maskset = THD_open_dataset( argv[++iarg] ) ;
00218 if( maskset == NULL )
00219 Syntax("-mask: can't open dataset" ) ;
00220
00221 iarg++ ; continue ;
00222 }
00223
00224
00225
00226 if( strcmp(argv[iarg],"-dimen") == 0 ){
00227 if( iarg+3 >= argc )
00228 Syntax("-dimen: don't have 3 arguments following!?") ;
00229 else if( mset != NULL )
00230 Syntax("-dimen: conflicts with previous -master!") ;
00231 else if( dimen_ii > 0 )
00232 Syntax("-dimen: can't have two -dimen options!") ;
00233 dimen_ii = strtol(argv[++iarg],NULL,10) ;
00234 dimen_jj = strtol(argv[++iarg],NULL,10) ;
00235 dimen_kk = strtol(argv[++iarg],NULL,10) ;
00236 if( dimen_ii < 2 || dimen_jj < 2 || dimen_kk < 2 )
00237 Syntax("-dimen: values following are not all >= 2!") ;
00238
00239 iarg++ ; continue ;
00240 }
00241
00242
00243
00244 if( strcmp(argv[iarg],"-datum") == 0 ){
00245 if( ++iarg >= argc )
00246 Syntax("-datum: no argument follows?!") ;
00247
00248 if( strcmp(argv[iarg],"short") == 0 )
00249 datum = MRI_short ;
00250 else if( strcmp(argv[iarg],"float") == 0 )
00251 datum = MRI_float ;
00252 else if( strcmp(argv[iarg],"byte") == 0 )
00253 datum = MRI_byte ;
00254 else
00255 Syntax("-datum: illegal type given!") ;
00256
00257 iarg++ ; continue ;
00258 }
00259
00260
00261
00262 if( strcmp(argv[iarg],"-srad") == 0 ){
00263 if( iarg+1 >= argc )
00264 Syntax("-srad: no argument follows!?") ;
00265
00266 srad = strtod( argv[++iarg] , NULL ) ;
00267 if( srad <= 0.0 ){
00268 fprintf(stderr,"++ WARNING: -srad value of %g is ignored!\n",srad);
00269 srad = 0.0 ;
00270 }
00271 iarg++ ; continue ;
00272 }
00273
00274
00275
00276 if( strcmp(argv[iarg],"-dval") == 0 ){
00277 if( iarg+1 >= argc )
00278 Syntax("-dval: no argument follows!?") ;
00279
00280 dval_float = strtod( argv[++iarg] , NULL ) ;
00281 dval_short = (short) rint(dval_float) ;
00282 dval_byte = (byte) dval_short ;
00283 iarg++ ; continue ;
00284 }
00285
00286
00287
00288 if( strcmp(argv[iarg],"-fval") == 0 ){
00289 if( iarg+1 >= argc )
00290 Syntax("-fval: no argument follows!?") ;
00291
00292 fval_float = strtod( argv[++iarg] , NULL ) ;
00293 fval_short = (short) rint(fval_float) ;
00294 fval_byte = (byte) fval_short ;
00295 iarg++ ; continue ;
00296 }
00297
00298 if( strcmp(argv[iarg],"-ijk") == 0 ){
00299 do_ijk = 1 ;
00300 iarg++ ; continue ;
00301 }
00302
00303
00304
00305 if( strcmp(argv[iarg],"-xyz") == 0 ){
00306 do_ijk = 0 ;
00307 iarg++ ; continue ;
00308 }
00309
00310
00311
00312 if( strcmp(argv[iarg],"-orient") == 0 ){
00313 int xx,yy,zz ;
00314 if( iarg+1 >= argc )
00315 Syntax("-orient: no argument follows!?") ;
00316
00317 orcode = argv[++iarg] ;
00318 if( strlen(orcode) != 3 )
00319 Syntax("-orient: illegal argument follows") ;
00320
00321 xx = ORCODE(orcode[0]) ; yy = ORCODE(orcode[1]) ; zz = ORCODE(orcode[2]) ;
00322 if( xx < 0 || yy < 0 || zz < 0 || !OR3OK(xx,yy,zz) )
00323 Syntax("-orient: illegal argument follows") ;
00324
00325 iarg++ ; continue ;
00326 }
00327
00328
00329
00330 fprintf(stderr,"*** Unknown option: %s\n",argv[iarg]) ; exit(1) ;
00331
00332 }
00333
00334
00335
00336 if( iarg >= argc )
00337 Syntax("No input files on command line!?") ;
00338
00339 if( do_ijk == 0 && mset == NULL )
00340 Syntax("Can't use -xyz without -master also!") ;
00341
00342 if( mset == NULL && dimen_ii < 2 )
00343 Syntax("Must use exactly one of -master or -dimen options on command line");
00344
00345 if( (datum == MRI_short && dval_short == fval_short) ||
00346 (datum == MRI_float && dval_float == fval_float) ||
00347 (datum == MRI_byte && dval_byte == fval_byte ) ){
00348
00349 fprintf(stderr,"+++ Warning: -dval and -fval are the same!\n") ;
00350 }
00351
00352
00353
00354 if( mset != NULL && do_ijk == 0 && orcode == NULL ){
00355 orcode = malloc(4) ;
00356 orcode[0] = ORIENT_typestr[mset->daxes->xxorient][0] ;
00357 orcode[1] = ORIENT_typestr[mset->daxes->yyorient][0] ;
00358 orcode[2] = ORIENT_typestr[mset->daxes->zzorient][0] ;
00359 orcode[3] = '\0' ;
00360 }
00361
00362 THD_coorder_fill( orcode , &cord ) ;
00363
00364
00365
00366 if( mset != NULL ){
00367
00368 dset = EDIT_empty_copy( mset ) ;
00369 EDIT_dset_items( dset ,
00370 ADN_prefix , prefix ,
00371 ADN_datum_all , datum ,
00372 ADN_nvals , 1 ,
00373 ADN_ntt , 0 ,
00374 ADN_func_type , ISANAT(mset) ? mset->func_type
00375 : FUNC_FIM_TYPE ,
00376
00377 ADN_directory_name , "./" ,
00378 ADN_none ) ;
00379
00380 } else {
00381
00382 THD_ivec3 iv_nxyz , iv_xyzorient ;
00383 THD_fvec3 fv_xyzorg , fv_xyzdel ;
00384
00385 LOAD_IVEC3( iv_nxyz , dimen_ii , dimen_jj , dimen_kk ) ;
00386 LOAD_IVEC3( iv_xyzorient , cord.xxor , cord.yyor , cord.zzor ) ;
00387 LOAD_FVEC3( fv_xyzdel ,
00388 ORIENT_sign[iv_xyzorient.ijk[0]]=='+' ? 1.0 : -1.0 ,
00389 ORIENT_sign[iv_xyzorient.ijk[1]]=='+' ? 1.0 : -1.0 ,
00390 ORIENT_sign[iv_xyzorient.ijk[2]]=='+' ? 1.0 : -1.0 ) ;
00391 LOAD_FVEC3( fv_xyzorg ,
00392 ORIENT_sign[iv_xyzorient.ijk[0]]=='+' ? -0.5*dimen_ii : 0.5*dimen_ii,
00393 ORIENT_sign[iv_xyzorient.ijk[1]]=='+' ? -0.5*dimen_jj : 0.5*dimen_jj,
00394 ORIENT_sign[iv_xyzorient.ijk[2]]=='+' ? -0.5*dimen_kk : 0.5*dimen_kk ) ;
00395
00396 dset = EDIT_empty_copy( NULL ) ;
00397
00398 EDIT_dset_items( dset ,
00399 ADN_nxyz , iv_nxyz ,
00400 ADN_xyzdel , fv_xyzdel ,
00401 ADN_xyzorg , fv_xyzorg ,
00402 ADN_xyzorient , iv_xyzorient ,
00403 ADN_prefix , prefix ,
00404 ADN_datum_all , datum ,
00405 ADN_nvals , 1 ,
00406 ADN_ntt , 0 ,
00407 ADN_type , HEAD_FUNC_TYPE ,
00408 ADN_func_type , FUNC_FIM_TYPE ,
00409 ADN_none ) ;
00410 }
00411
00412 if( THD_is_file(DSET_HEADNAME(dset)) )
00413 Syntax("Output dataset already exists -- can't overwrite") ;
00414
00415
00416
00417 EDIT_substitute_brick( dset , 0 , datum , NULL ) ;
00418
00419 nx = DSET_NX(dset); ny = DSET_NY(dset); nz = DSET_NZ(dset); nxyz = nx*ny*nz;
00420
00421
00422
00423 if( maskset != NULL &&
00424 ( DSET_NX(maskset) != nx ||
00425 DSET_NY(maskset) != ny ||
00426 DSET_NZ(maskset) != nz ) )
00427 Syntax("-mask dataset doesn't match dimension of output dataset") ;
00428
00429 if( maskset != NULL ){
00430 mmask = THD_makemask( maskset , 0 , 1.0,-1.0 ) ;
00431 if( mmask == NULL ){
00432 fprintf(stderr,"++ WARNING: can't create mask for some reason!\n") ;
00433 } else {
00434 int nmask = THD_countmask( nxyz , mmask ) ;
00435 if( nmask == 0 ){
00436 fprintf(stderr,"++ WARNING: 0 voxels in mask -- ignoring it!\n") ;
00437 free((void *)mmask) ; mmask = NULL ;
00438 } else {
00439 fprintf(stderr,"++ %d voxels found in mask\n",nmask) ;
00440 }
00441 }
00442 DSET_delete(maskset) ;
00443 }
00444
00445
00446
00447 switch( datum ){
00448 case MRI_short:
00449 sbr = (short *) DSET_BRICK_ARRAY(dset,0) ;
00450 for( ii=0 ; ii < nxyz ; ii++ ) sbr[ii] = fval_short ;
00451 break ;
00452
00453 case MRI_float:
00454 fbr = (float *) DSET_BRICK_ARRAY(dset,0) ;
00455 for( ii=0 ; ii < nxyz ; ii++ ) fbr[ii] = fval_float ;
00456 break ;
00457
00458 case MRI_byte:
00459 bbr = (byte *) DSET_BRICK_ARRAY(dset,0) ;
00460 for( ii=0 ; ii < nxyz ; ii++ ) bbr[ii] = fval_byte ;
00461 break ;
00462 }
00463
00464
00465
00466 dx = fabs(dset->daxes->xxdel) ; if( dx <= 0.0 ) dx = 1.0 ;
00467 dy = fabs(dset->daxes->yydel) ; if( dy <= 0.0 ) dy = 1.0 ;
00468 dz = fabs(dset->daxes->zzdel) ; if( dz <= 0.0 ) dz = 1.0 ;
00469
00470 if( !do_ijk ){
00471 #ifndef EXTEND_BBOX
00472 xxdown = dset->daxes->xxmin - 0.501 * dx ;
00473 xxup = dset->daxes->xxmax + 0.501 * dx ;
00474 yydown = dset->daxes->yymin - 0.501 * dy ;
00475 yyup = dset->daxes->yymax + 0.501 * dy ;
00476 zzdown = dset->daxes->zzmin - 0.501 * dz ;
00477 zzup = dset->daxes->zzmax + 0.501 * dz ;
00478 #else
00479 xxdown = dset->daxes->xxmin ;
00480 xxup = dset->daxes->xxmax ;
00481 yydown = dset->daxes->yymin ;
00482 yyup = dset->daxes->yymax ;
00483 zzdown = dset->daxes->zzmin ;
00484 zzup = dset->daxes->zzmax ;
00485 #endif
00486 }
00487
00488
00489
00490 for( ; iarg < argc ; iarg++ ){
00491
00492
00493
00494 if( strcmp(argv[iarg],"-") == 0 ){
00495 fp = stdin ;
00496 } else {
00497 fp = fopen( argv[iarg] , "r" ) ;
00498 if( fp == NULL ){
00499 fprintf(stderr,
00500 "+++ Warning: can't open input file %s -- skipping it\n",
00501 argv[iarg]) ;
00502 continue ;
00503 }
00504 }
00505
00506
00507
00508 ll = 0 ;
00509 while(1){
00510 ll++ ;
00511 cp = fgets( linbuf , NBUF , fp ) ;
00512 if( cp == NULL ) break ;
00513 kk = strlen(linbuf) ;
00514 if( kk == 0 ) continue ;
00515
00516
00517
00518 for( ii=0 ; ii < kk && isspace(linbuf[ii]) ; ii++ ) ;
00519 if( ii == kk ) continue ;
00520 if( linbuf[ii] == '/' && linbuf[ii+1] == '/' ) continue ;
00521 if( linbuf[ii] == '#' ) continue ;
00522
00523
00524
00525 vv = dval_float ;
00526 vrad = srad ;
00527 nn = sscanf(linbuf+ii , "%f%f%f%f%f" , &xx,&yy,&zz,&vv,&vrad ) ;
00528 if( nn < 3 ){
00529 fprintf(stderr,"+++ Warning: file %s line %d: incomplete\n",argv[iarg],ll) ;
00530 continue ;
00531 }
00532
00533
00534
00535 if( do_ijk ){
00536
00537 ii = (int) rint(xx) ; jj = (int) rint(yy) ; kk = (int) rint(zz) ;
00538 if( ii < 0 || ii >= nx ){
00539 fprintf(stderr,
00540 "+++ Warning: file %s line %d: i index=%d is invalid\n",
00541 argv[iarg],ll,ii) ;
00542 continue ;
00543 }
00544 if( jj < 0 || jj >= ny ){
00545 fprintf(stderr,
00546 "+++ Warning: file %s line %d: j index=%d is invalid\n",
00547 argv[iarg],ll,jj) ;
00548 continue ;
00549 }
00550 if( kk < 0 || kk >= nz ){
00551 fprintf(stderr,
00552 "+++ Warning: file %s line %d: k index=%d is invalid\n",
00553 argv[iarg],ll,kk) ;
00554 continue ;
00555 }
00556
00557 } else {
00558
00559 THD_fvec3 mv , dv ;
00560 THD_ivec3 iv ;
00561
00562 THD_coorder_to_dicom( &cord , &xx,&yy,&zz ) ;
00563 LOAD_FVEC3( dv , xx,yy,zz ) ;
00564 mv = THD_dicomm_to_3dmm( dset , dv ) ;
00565
00566
00567
00568 if( mv.xyz[0] < xxdown || mv.xyz[0] > xxup ){
00569 fprintf(stderr,"+++ Warning: file %s line %d: x coord=%g is outside %g .. %g\n" ,
00570 argv[iarg],ll,mv.xyz[0] , xxdown,xxup ) ;
00571 continue ;
00572 }
00573 if( mv.xyz[1] < yydown || mv.xyz[1] > yyup ){
00574 fprintf(stderr,"+++ Warning: file %s line %d: y coord=%g is outside %g .. %g\n" ,
00575 argv[iarg],ll,mv.xyz[1] , yydown , yyup ) ;
00576 continue ;
00577 }
00578 if( mv.xyz[2] < zzdown || mv.xyz[2] > zzup ){
00579 fprintf(stderr,"+++ Warning: file %s line %d: z coord=%g is outside %g .. %g\n" ,
00580 argv[iarg],ll,mv.xyz[2] , zzdown , zzup ) ;
00581 continue ;
00582 }
00583
00584 iv = THD_3dmm_to_3dind( dset , mv ) ;
00585 ii = iv.ijk[0]; jj = iv.ijk[1]; kk = iv.ijk[2];
00586 }
00587
00588
00589
00590 ijk = ii + jj*nx + kk*nx*ny ;
00591 if( mmask == NULL || mmask[ijk] ){
00592 switch( datum ){
00593 case MRI_float:{
00594 if( fbr[ijk] != fval_float && fbr[ijk] != vv )
00595 fprintf(stderr,"Overwrite voxel %d %d %d\n",ii,jj,kk) ;
00596 fbr[ijk] = vv ;
00597 }
00598 break ;
00599 case MRI_short:{
00600 sv = SHORTIZE(vv) ;
00601 if( sbr[ijk] != fval_short && sbr[ijk] != sv )
00602 fprintf(stderr,"Overwrite voxel %d %d %d\n",ii,jj,kk) ;
00603 sbr[ijk] = sv ;
00604 }
00605 break ;
00606 case MRI_byte:{
00607 bv = BYTEIZE(vv) ;
00608 if( bbr[ijk] != fval_byte && bbr[ijk] != bv )
00609 fprintf(stderr,"Overwrite voxel %d %d %d\n",ii,jj,kk) ;
00610 bbr[ijk] = bv ;
00611 }
00612 break ;
00613 }
00614 }
00615
00616
00617
00618
00619
00620
00621 vrad *= 1.00001 ;
00622 rii = vrad/dx ; rjj = vrad/dy ; rkk = vrad/dz ;
00623 qii = rii*rii ; qjj = rjj*rjj ; qkk = rkk*rkk ;
00624
00625 if( rii >= 1.0 || rjj >= 1.0 || rkk >= 1.0 ){
00626 int aa,bb,cc , abot,atop,bbot,btop,cbot,ctop; float rr;
00627 abot = ii-(int)rint(rii) ; atop = ii+(int)rint(rii) ;
00628 if( abot < 0 ) abot = 0 ; if( atop >= nx ) atop = nx-1 ;
00629 bbot = jj-(int)rint(rjj) ; btop = jj+(int)rint(rjj) ;
00630 if( bbot < 0 ) bbot = 0 ; if( btop >= ny ) btop = ny-1 ;
00631 cbot = kk-(int)rint(rkk) ; ctop = kk+(int)rint(rkk) ;
00632 if( cbot < 0 ) cbot = 0 ; if( ctop >= nz ) ctop = nz-1 ;
00633 for( cc=cbot ; cc <= ctop ; cc++ ){
00634 for( bb=bbot ; bb <= btop ; bb++ ){
00635 for( aa=abot ; aa <= atop ; aa++ ){
00636 rr = (aa-ii)*(aa-ii)/qii
00637 + (bb-jj)*(bb-jj)/qjj + (cc-kk)*(cc-kk)/qkk ;
00638 if( rr <= 1.00001 ){
00639 ijk = aa + bb*nx + cc*nx*ny ;
00640 if( mmask == NULL || mmask[ijk] ){
00641 switch( datum ){
00642 case MRI_float: fbr[ijk] = vv ; break ;
00643 case MRI_short: sbr[ijk] = sv ; break ;
00644 case MRI_byte: bbr[ijk] = bv ; break ;
00645 }
00646 }
00647 }
00648 }
00649 }
00650 }
00651 }
00652
00653 }
00654
00655
00656
00657 if( fp != stdin ) fclose( fp ) ;
00658
00659 }
00660
00661 tross_Make_History( "3dUndump" , argc,argv , dset ) ;
00662 DSET_write(dset) ;
00663 fprintf(stderr,"+++ Wrote results to dataset %s\n",DSET_BRIKNAME(dset)) ;
00664 exit(0) ;
00665 }
