00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009 #if 0
00010 # define DB(sss,str) printf("%s %s\n",sss,str)
00011 #else
00012 # define DB(sss,str)
00013 #endif
00014
00015 static EDIT_options PRED_edopt ;
00016
00017
00018
00019 typedef struct {
00020 float xbot , xtop , ybot , ytop , zbot , ztop ;
00021 } dset_range ;
00022
00023 dset_range PR_get_range( THD_3dim_dataset * ) ;
00024 void PR_one_slice( int , MRI_IMAGE * , MRI_IMAGE * ) ;
00025 float PR_type_scale( int , MRI_IMAGE * ) ;
00026
00027 #define PRED_err(str) \
00028 do{ fprintf(stderr,"\n*** %s\n",str) ; exit(-1) ; } while(1)
00029
00030
00031
00032 void Syntax(void)
00033 {
00034 printf(
00035 "Projection along cardinal axes from a 3D dataset\n"
00036 "Usage: 3dproject [editing options]\n"
00037 " [-sum|-max|-amax|-smax] [-output root] [-nsize] [-mirror]\n"
00038 " [-RL {all | x1 x2}] [-AP {all | y1 y2}] [-IS {all | z1 z2}]\n"
00039 " [-ALL] dataset\n"
00040 "\n"
00041 "Program to produce orthogonal projections from a 3D dataset.\n"
00042 " -sum ==> Add the dataset voxels along the projection direction\n"
00043 " -max ==> Take the maximum of the voxels [the default is -sum]\n"
00044 " -amax ==> Take the absolute maximum of the voxels\n"
00045 " -smax ==> Take the signed maximum of the voxels; for example,\n"
00046 " -max ==> -7 and 2 go to 2 as the projected value\n"
00047 " -amax ==> -7 and 2 go to 7 as the projected value\n"
00048 " -smax ==> -7 and 2 go to -7 as the projected value\n"
00049 " -first x ==> Take the first value greater than x\n"
00050 " -nsize ==> Scale the output images up to 'normal' sizes\n"
00051 " (e.g., 64x64, 128x128, or 256x256)\n"
00052 " This option only applies to byte or short datasets.\n"
00053 " -mirror ==> The radiologists' and AFNI convention is to display\n"
00054 " axial and coronal images with the subject's left on\n"
00055 " the right of the image; the use of this option will\n"
00056 " mirror the axial and coronal projections so that\n"
00057 " left is left and right is right.\n"
00058 "\n"
00059 " -output root ==> Output projections will named\n"
00060 " root.sag, root.cor, and root.axi\n"
00061 " [the default root is 'proj']\n"
00062 "\n"
00063 " -RL all ==> Project in the Right-to-Left direction along\n"
00064 " all the data (produces root.sag)\n"
00065 " -RL x1 x2 ==> Project in the Right-to-Left direction from\n"
00066 " x-coordinate x1 to x2 (mm)\n"
00067 " [negative x is Right, positive x is Left]\n"
00068 " [OR, you may use something like -RL 10R 20L\n"
00069 " to project from x=-10 mm to x=+20 mm ]\n"
00070 "\n"
00071 " -AP all ==> Project in the Anterior-to-Posterior direction along\n"
00072 " all the data (produces root.cor)\n"
00073 " -AP y1 y2 ==> Project in the Anterior-to-Posterior direction from\n"
00074 " y-coordinate y1 to y2 (mm)\n"
00075 " [negative y is Anterior, positive y is Posterior]\n"
00076 " [OR, you may use something like -AP 10A 20P\n"
00077 " to project from y=-10 mm to y=+20 mm ]\n"
00078 "\n"
00079 " -IS all ==> Project in the Inferior-to-Superior direction along\n"
00080 " all the data (produces root.axi)\n"
00081 " -IS y1 y2 ==> Project in the Inferior-to-Superior direction from\n"
00082 " z-coordinate z1 to z2 (mm)\n"
00083 " [negative z is Inferior, positive z is Superior]\n"
00084 " [OR, you may use something like -IS 10I 20S\n"
00085 " to project from z=-10 mm to z=+20 mm ]\n"
00086 "\n"
00087 " -ALL ==> Equivalent to '-RL all -AP all -IS all'\n"
00088 "\n"
00089 "* NOTE that a projection direction will not be used if the bounds aren't\n"
00090 " given for that direction; thus, at least one of -RL, -AP, or -IS must\n"
00091 " be used, or nothing will be computed!\n"
00092 "* NOTE that in the directions transverse to the projection direction,\n"
00093 " all the data is used; that is, '-RL -5 5' will produce a full sagittal\n"
00094 " image summed over a 10 mm slice, irrespective of the -IS or -AP extents.\n"
00095 "* NOTE that the [editing options] are the same as in 3dmerge.\n"
00096 " In particular, the '-1thtoin' option can be used to project the\n"
00097 " threshold data (if available).\n"
00098 ) ;
00099
00100 printf("\n" MASTER_SHORTHELP_STRING ) ;
00101
00102 exit(0) ;
00103 }
00104
00105 #define PROJ_SUM 0
00106 #define PROJ_MMAX 1
00107 #define PROJ_AMAX 2
00108 #define PROJ_SMAX 3
00109 #define PROJ_FIRST 4
00110
00111 #define MIRR_NO 0
00112 #define MIRR_YES 1
00113
00114 #define BIGG 9999999.99
00115
00116 static float first_thresh=0.0 ;
00117
00118 int main( int argc , char * argv[] )
00119 {
00120 THD_3dim_dataset * dset ;
00121 THD_dataxes * daxes ;
00122 FD_brick ** brarr , * baxi , * bsag , * bcor ;
00123
00124 int iarg , ii ;
00125 Boolean ok ;
00126 MRI_IMAGE * pim , * flim , * slim ;
00127 float * flar ;
00128 dset_range dr ;
00129 float val , fimfac ;
00130 int ival,ityp , kk ;
00131
00132 float xbot=BIGG,xtop=BIGG , ybot=BIGG,ytop=BIGG , zbot=BIGG,ztop=BIGG ;
00133 int xgood=0 , ygood=0 , zgood=0 ,
00134 proj_code=PROJ_SUM , mirror_code=MIRR_NO , nsize=0 ;
00135 char root[THD_MAX_NAME] = "proj." ;
00136 char fname[THD_MAX_NAME] ;
00137
00138 int ixbot,ixtop , jybot,jytop , kzbot,kztop ;
00139 THD_fvec3 fv ;
00140 THD_ivec3 iv ;
00141
00142
00143
00144 if( argc < 2 || strncmp(argv[1],"-help",6) == 0 ) Syntax() ;
00145
00146 INIT_EDOPT( &PRED_edopt ) ;
00147 iarg = 1 ;
00148 while( iarg < argc && argv[iarg][0] == '-' ){
00149
00150 DB("new arg:",argv[iarg]) ;
00151
00152
00153
00154 ii = EDIT_check_argv( argc , argv , iarg , &PRED_edopt ) ;
00155 if( ii > 0 ){
00156 iarg += ii ;
00157 continue ;
00158 }
00159
00160
00161
00162 if( strncmp(argv[iarg],"-sum",6) == 0 ){
00163 proj_code = PROJ_SUM ;
00164 iarg++ ; continue ;
00165 }
00166
00167 if( strncmp(argv[iarg],"-max",6) == 0 ){
00168 proj_code = PROJ_MMAX ;
00169 iarg++ ; continue ;
00170 }
00171
00172 if( strncmp(argv[iarg],"-amax",6) == 0 ){
00173 proj_code = PROJ_AMAX ;
00174 iarg++ ; continue ;
00175 }
00176
00177 if( strncmp(argv[iarg],"-smax",6) == 0 ){
00178 proj_code = PROJ_SMAX ;
00179 iarg++ ; continue ;
00180 }
00181
00182 if( strcmp(argv[iarg],"-first") == 0 ){
00183 proj_code = PROJ_FIRST ;
00184 first_thresh = strtod( argv[++iarg] , NULL ) ;
00185 iarg++ ; continue ;
00186 }
00187
00188
00189
00190 if( strncmp(argv[iarg],"-mirror",6) == 0 ){
00191 mirror_code = MIRR_YES ;
00192 iarg++ ; continue ;
00193 }
00194
00195
00196
00197 if( strncmp(argv[iarg],"-nsize",6) == 0 ){
00198 nsize = 1 ;
00199 iarg++ ; continue ;
00200 }
00201
00202
00203
00204 if( strncmp(argv[iarg],"-output",6) == 0 ||
00205 strncmp(argv[iarg],"-root",6) == 0 ){
00206
00207 if( iarg+1 >= argc ){
00208 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
00209 exit(-1) ;
00210 }
00211
00212 MCW_strncpy( root , argv[++iarg] , THD_MAX_NAME-1 ) ;
00213 ii = strlen(root) ;
00214 if( ii == 0 ){
00215 fprintf(stderr,"\n*** illegal rootname!\n") ; exit(-1) ;
00216 }
00217 if( root[ii-1] != '.' ){ root[ii] = '.' ; root[ii+1] = '\0' ; }
00218 iarg++ ; continue ;
00219 }
00220
00221
00222
00223 if( strncmp(argv[iarg],"-ALL",6) == 0 ||
00224 strncmp(argv[iarg],"-all",6) == 0 ){
00225
00226 xgood = ygood = zgood = 1 ;
00227 xbot = ybot = zbot = -BIGG ;
00228 xtop = ytop = ztop = BIGG ;
00229 iarg++ ; continue ;
00230 }
00231
00232
00233
00234 if( strncmp(argv[iarg],"-RL",6) == 0 ||
00235 strncmp(argv[iarg],"-LR",6) == 0 ||
00236 strncmp(argv[iarg],"-rl",6) == 0 ||
00237 strncmp(argv[iarg],"-lr",6) == 0 ||
00238 strncmp(argv[iarg],"-sag",6)== 0 ){
00239
00240 char * cerr ; float tf ;
00241
00242 xgood = 1 ;
00243
00244 if( iarg+1 >= argc ){
00245 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
00246 exit(-1) ;
00247 }
00248
00249 if( strncmp(argv[iarg+1],"all",6) == 0 ||
00250 strncmp(argv[iarg+1],"ALL",6) == 0 ){
00251
00252 xbot = -BIGG;
00253 xtop = BIGG ;
00254 iarg += 2 ; continue ;
00255 }
00256
00257 if( iarg+2 >= argc ){
00258 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
00259 exit(-1) ;
00260 }
00261
00262 xbot = strtod( argv[iarg+1] , &cerr ) ;
00263 if( cerr == argv[iarg+1] ){
00264 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
00265 argv[iarg],argv[iarg+1] ) ; exit(-1) ;
00266 }
00267 if( *cerr == 'R' && xbot > 0.0 ) xbot = -xbot ;
00268
00269 xtop = strtod( argv[iarg+2] , &cerr ) ;
00270 if( cerr == argv[iarg+2] ){
00271 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
00272 argv[iarg],argv[iarg+2] ) ; exit(-1) ;
00273 }
00274 if( *cerr == 'R' && xtop > 0.0 ) xtop = -xtop ;
00275
00276 if( xbot > xtop ){ tf = xbot ; xbot = xtop ; xtop = tf ; }
00277 iarg +=3 ; continue ;
00278 }
00279
00280
00281
00282 if( strncmp(argv[iarg],"-AP",6) == 0 ||
00283 strncmp(argv[iarg],"-PA",6) == 0 ||
00284 strncmp(argv[iarg],"-ap",6) == 0 ||
00285 strncmp(argv[iarg],"-pa",6) == 0 ||
00286 strncmp(argv[iarg],"-cor",6)== 0 ){
00287
00288 char * cerr ; float tf ;
00289
00290 ygood = 1 ;
00291
00292 if( iarg+1 >= argc ){
00293 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
00294 exit(-1) ;
00295 }
00296
00297 if( strncmp(argv[iarg+1],"all",6) == 0 ||
00298 strncmp(argv[iarg+1],"ALL",6) == 0 ){
00299
00300 ybot = -BIGG ;
00301 ytop = BIGG ;
00302 iarg += 2 ; continue ;
00303 }
00304
00305 if( iarg+2 >= argc ){
00306 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
00307 exit(-1) ;
00308 }
00309
00310 ybot = strtod( argv[iarg+1] , &cerr ) ;
00311 if( cerr == argv[iarg+1] ){
00312 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
00313 argv[iarg],argv[iarg+1] ) ; exit(-1) ;
00314 }
00315 if( *cerr == 'A' && ybot > 0.0 ) ybot = -ybot ;
00316
00317 ytop = strtod( argv[iarg+2] , &cerr ) ;
00318 if( cerr == argv[iarg+2] ){
00319 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
00320 argv[iarg],argv[iarg+2] ) ; exit(-1) ;
00321 }
00322 if( *cerr == 'A' && ytop > 0.0 ) ytop = -ytop ;
00323
00324 if( ybot > ytop ){ tf = ybot ; ybot = ytop ; ytop = tf ; }
00325 iarg +=3 ; continue ;
00326 }
00327
00328
00329
00330 if( strncmp(argv[iarg],"-IS",6) == 0 ||
00331 strncmp(argv[iarg],"-SI",6) == 0 ||
00332 strncmp(argv[iarg],"-is",6) == 0 ||
00333 strncmp(argv[iarg],"-si",6) == 0 ||
00334 strncmp(argv[iarg],"-axi",6)== 0 ){
00335
00336 char * cerr ; float tf ;
00337
00338 zgood = 1 ;
00339
00340 if( iarg+1 >= argc ){
00341 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
00342 exit(-1) ;
00343 }
00344
00345 if( strncmp(argv[iarg+1],"all",6) == 0 ||
00346 strncmp(argv[iarg+1],"ALL",6) == 0 ){
00347
00348 zbot = -BIGG ;
00349 ztop = BIGG ;
00350 iarg += 2 ; continue ;
00351 }
00352
00353 if( iarg+2 >= argc ){
00354 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
00355 exit(-1) ;
00356 }
00357
00358 zbot = strtod( argv[iarg+1] , &cerr ) ;
00359 if( cerr == argv[iarg+1] ){
00360 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
00361 argv[iarg],argv[iarg+1] ) ; exit(-1) ;
00362 }
00363 if( *cerr == 'I' && zbot > 0.0 ) zbot = -zbot ;
00364
00365 ztop = strtod( argv[iarg+2] , &cerr ) ;
00366 if( cerr == argv[iarg+2] ){
00367 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
00368 argv[iarg],argv[iarg+2] ) ; exit(-1) ;
00369 }
00370 if( *cerr == 'I' && ztop > 0.0 ) ztop = -ztop ;
00371
00372 if( zbot > ztop ){ tf = zbot ; zbot = ztop ; ztop = tf ; }
00373 iarg +=3 ; continue ;
00374 }
00375
00376
00377
00378 fprintf(stderr,"\n*** Unknown option: %s\n",argv[iarg]) ;
00379 exit(-1) ;
00380 }
00381
00382 if( ! xgood && ! ygood && ! zgood ){
00383 fprintf(stderr,"\n*** No projections ordered!?\n") ; exit(-1) ;
00384 }
00385
00386
00387
00388 dset = THD_open_dataset( argv[iarg] ) ;
00389 if( dset == NULL ){
00390 fprintf(stderr,"\n*** Can't open dataset file %s\n",argv[iarg]) ;
00391 exit(-1) ;
00392 }
00393 if( DSET_NUM_TIMES(dset) > 1 ){
00394 fprintf(stderr,"\n*** Can't project time-dependent dataset!\n") ;
00395 exit(1) ;
00396 }
00397 EDIT_one_dataset( dset, &PRED_edopt ) ;
00398 daxes = dset->daxes ;
00399 brarr = THD_setup_bricks( dset ) ;
00400 baxi = brarr[0] ; bsag = brarr[1] ; bcor = brarr[2] ;
00401
00402
00403
00404 dr = PR_get_range( dset ) ;
00405
00406 if( xgood ){
00407 if( xbot < dr.xbot ) xbot = dr.xbot ;
00408 if( xtop > dr.xtop ) xtop = dr.xtop ;
00409
00410 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(xbot,0,0) ) ;
00411 iv = THD_3dmm_to_3dind ( dset , fv ) ;
00412 iv = THD_3dind_to_fdind( bsag , iv ) ; ixbot = iv.ijk[2] ;
00413
00414 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(xtop,0,0) ) ;
00415 iv = THD_3dmm_to_3dind ( dset , fv ) ;
00416 iv = THD_3dind_to_fdind( bsag , iv ) ; ixtop = iv.ijk[2] ;
00417
00418 if( ixbot > ixtop ) { ii = ixbot ; ixbot = ixtop ; ixtop = ii ; }
00419
00420 if( ixbot < 0 ) ixbot = 0 ;
00421 if( ixtop >= bsag->n3 ) ixtop = bsag->n3 - 1 ;
00422 }
00423
00424 if( ygood ){
00425 if( ybot < dr.ybot ) ybot = dr.ybot ;
00426 if( ytop > dr.ytop ) ytop = dr.ytop ;
00427
00428 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(0,ybot,0) ) ;
00429 iv = THD_3dmm_to_3dind ( dset , fv ) ;
00430 iv = THD_3dind_to_fdind( bcor , iv ) ; jybot = iv.ijk[2] ;
00431
00432 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(0,ytop,0) ) ;
00433 iv = THD_3dmm_to_3dind ( dset , fv ) ;
00434 iv = THD_3dind_to_fdind( bcor , iv ) ; jytop = iv.ijk[2] ;
00435
00436 if( jybot > jytop ) { ii = jybot ; jybot = jytop ; jytop = ii ; }
00437
00438 if( jybot < 0 ) jybot = 0 ;
00439 if( jytop >= bcor->n3 ) jytop = bcor->n3 - 1 ;
00440 }
00441
00442 if( zgood ){
00443 if( zbot < dr.zbot ) zbot = dr.zbot ;
00444 if( ztop > dr.ztop ) ztop = dr.ztop ;
00445
00446 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(0,0,zbot) ) ;
00447 iv = THD_3dmm_to_3dind ( dset , fv ) ;
00448 iv = THD_3dind_to_fdind( baxi , iv ) ; kzbot = iv.ijk[2] ;
00449
00450 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(0,0,ztop) ) ;
00451 iv = THD_3dmm_to_3dind ( dset , fv ) ;
00452 iv = THD_3dind_to_fdind( baxi , iv ) ; kztop = iv.ijk[2] ;
00453
00454 if( kzbot > kztop ) { ii = kzbot ; kzbot = kztop ; kztop = ii ; }
00455
00456 if( kzbot < 0 ) kzbot = 0 ;
00457 if( kztop >= baxi->n3 ) kztop = baxi->n3 - 1 ;
00458 }
00459
00460 ival = DSET_PRINCIPAL_VALUE(dset) ;
00461 ityp = DSET_BRICK_TYPE(dset,ival) ;
00462 fimfac = DSET_BRICK_FACTOR(dset,ival) ;
00463
00464
00465
00466 if( xgood ){
00467 int n1 = bsag->n1 , n2 = bsag->n2 ;
00468 int ss , npix ;
00469 float fmax , fmin , scl ;
00470
00471
00472
00473 npix = n1*n2 ;
00474 flim = mri_new( n1 , n2 , MRI_float ) ;
00475 flar = mri_data_pointer( flim ) ;
00476 for( ii=0 ; ii < npix ; ii++ ) flar[ii] = 0.0 ;
00477
00478
00479
00480 for( ss=ixbot ; ss <= ixtop ; ss++ ){
00481 slim = FD_brick_to_mri( ss , ival , bsag ) ;
00482 if( slim->kind != MRI_float ){
00483 pim = mri_to_float( slim ) ;
00484 mri_free( slim ) ; slim = pim ;
00485 }
00486 PR_one_slice( proj_code , slim , flim ) ;
00487 mri_free( slim ) ;
00488 }
00489
00490
00491
00492 if( fimfac != 0.0 && fimfac != 1.0 ){
00493 slim = mri_scale_to_float( 1.0/fimfac , flim ) ;
00494 mri_free(flim) ; flim = slim ;
00495 }
00496
00497 scl = PR_type_scale( ityp , flim ) ;
00498 pim = mri_to_mri_scl( ityp , scl , flim ) ; mri_free( flim ) ;
00499 if( nsize ){
00500 slim = mri_nsize( pim ) ;
00501 if( slim != NULL && slim != pim ) { mri_free(pim) ; pim = slim ; }
00502 }
00503
00504 if( scl != 1.0 )
00505 printf("Sagittal projection pixels scaled by %g to avoid overflow!\n",
00506 scl ) ;
00507
00508 fmax = mri_max(pim) ; fmin = mri_min(pim) ;
00509 printf("Sagittal projection min = %g max = %g\n",fmin,fmax) ;
00510
00511 strcpy(fname,root) ; strcat(fname,"sag") ;
00512 mri_write( fname, pim ) ;
00513 mri_free(pim) ;
00514 }
00515
00516
00517
00518 if( ygood ){
00519 int n1 = bcor->n1 , n2 = bcor->n2 ;
00520 int ss , npix ;
00521 float fmax , fmin , scl ;
00522
00523
00524
00525 npix = n1*n2 ;
00526 flim = mri_new( n1 , n2 , MRI_float ) ;
00527 flar = mri_data_pointer( flim ) ;
00528 for( ii=0 ; ii < npix ; ii++ ) flar[ii] = 0.0 ;
00529
00530
00531
00532 for( ss=jybot ; ss <= jytop ; ss++ ){
00533 slim = FD_brick_to_mri( ss , ival , bcor ) ;
00534 if( slim->kind != MRI_float ){
00535 pim = mri_to_float( slim ) ;
00536 mri_free( slim ) ; slim = pim ;
00537 }
00538 PR_one_slice( proj_code , slim , flim ) ;
00539 mri_free( slim ) ;
00540 }
00541
00542
00543
00544 if( fimfac != 0.0 && fimfac != 1.0 ){
00545 slim = mri_scale_to_float( 1.0/fimfac , flim ) ;
00546 mri_free(flim) ; flim = slim ;
00547 }
00548
00549 scl = PR_type_scale( ityp , flim ) ;
00550 pim = mri_to_mri_scl( ityp , scl , flim ) ; mri_free( flim ) ;
00551 if( nsize ){
00552 slim = mri_nsize( pim ) ;
00553 if( slim != NULL && slim != pim ) { mri_free(pim) ; pim = slim ; }
00554 }
00555
00556 if( scl != 1.0 )
00557 printf("Coronal projection pixels scaled by %g to avoid overflow!\n",
00558 scl ) ;
00559
00560 fmax = mri_max(pim) ; fmin = mri_min(pim) ;
00561 printf("Coronal projection min = %g max = %g\n",fmin,fmax) ;
00562
00563 if( mirror_code == MIRR_YES ){
00564 slim = mri_flippo( MRI_ROT_0 , TRUE , pim ) ;
00565 mri_free(pim) ; pim = slim ;
00566 }
00567 strcpy(fname,root) ; strcat(fname,"cor") ;
00568 mri_write( fname, pim ) ;
00569 mri_free(pim) ;
00570 }
00571
00572
00573
00574 if( zgood ){
00575 int n1 = baxi->n1 , n2 = baxi->n2 ;
00576 int ss , npix ;
00577 float fmax , fmin , scl ;
00578
00579
00580
00581 npix = n1*n2 ;
00582 flim = mri_new( n1 , n2 , MRI_float ) ;
00583 flar = mri_data_pointer( flim ) ;
00584 for( ii=0 ; ii < npix ; ii++ ) flar[ii] = 0.0 ;
00585
00586
00587
00588 for( ss=kzbot ; ss <= kztop ; ss++ ){
00589 slim = FD_brick_to_mri( ss , ival , baxi ) ;
00590 if( slim->kind != MRI_float ){
00591 pim = mri_to_float( slim ) ;
00592 mri_free( slim ) ; slim = pim ;
00593 }
00594 PR_one_slice( proj_code , slim , flim ) ;
00595 mri_free( slim ) ;
00596 }
00597
00598
00599
00600 if( fimfac != 0.0 && fimfac != 1.0 ){
00601 slim = mri_scale_to_float( 1.0/fimfac , flim ) ;
00602 mri_free(flim) ; flim = slim ;
00603 }
00604
00605 scl = PR_type_scale( ityp , flim ) ;
00606 pim = mri_to_mri_scl( ityp , scl , flim ) ; mri_free( flim ) ;
00607 if( nsize ){
00608 slim = mri_nsize( pim ) ;
00609 if( slim != NULL && slim != pim ) { mri_free(pim) ; pim = slim ; }
00610 }
00611
00612 if( scl != 1.0 )
00613 printf("Axial projection pixels scaled by %g to avoid overflow!\n",
00614 scl ) ;
00615
00616 fmax = mri_max(pim) ; fmin = mri_min(pim) ;
00617 printf("Axial projection min = %g max = %g\n",fmin,fmax) ;
00618
00619 if( mirror_code == MIRR_YES ){
00620 slim = mri_flippo( MRI_ROT_0 , TRUE , pim ) ;
00621 mri_free(pim) ; pim = slim ;
00622 }
00623 strcpy(fname,root) ; strcat(fname,"axi") ;
00624 mri_write( fname, pim ) ;
00625 mri_free(pim) ;
00626 }
00627
00628 exit(0) ;
00629 }
00630
00631
00632
00633 dset_range PR_get_range( THD_3dim_dataset * dset )
00634 {
00635 THD_dataxes * daxes ;
00636 THD_fvec3 fv1 , fv2 , fv3 ;
00637 THD_ivec3 iv ;
00638 float tf ;
00639 dset_range dr ;
00640
00641 #define FSWAP(x,y) (tf=(x),(x)=(y),(y)=tf)
00642
00643 daxes = dset->daxes ;
00644
00645 LOAD_FVEC3(fv1 , daxes->xxorg , daxes->yyorg , daxes->zzorg) ;
00646 fv1 = THD_3dmm_to_dicomm( dset , fv1 ) ;
00647
00648 LOAD_FVEC3(fv2 , daxes->xxorg + (daxes->nxx-1)*daxes->xxdel ,
00649 daxes->yyorg + (daxes->nyy-1)*daxes->yydel ,
00650 daxes->zzorg + (daxes->nzz-1)*daxes->zzdel ) ;
00651 fv2 = THD_3dmm_to_dicomm( dset , fv2 ) ;
00652
00653 if( fv1.xyz[0] > fv2.xyz[0] ) FSWAP( fv1.xyz[0] , fv2.xyz[0] ) ;
00654 if( fv1.xyz[1] > fv2.xyz[1] ) FSWAP( fv1.xyz[1] , fv2.xyz[1] ) ;
00655 if( fv1.xyz[2] > fv2.xyz[2] ) FSWAP( fv1.xyz[2] , fv2.xyz[2] ) ;
00656
00657 dr.xbot = fv1.xyz[0] ; dr.xtop = fv2.xyz[0] ;
00658 dr.ybot = fv1.xyz[1] ; dr.ytop = fv2.xyz[1] ;
00659 dr.zbot = fv1.xyz[2] ; dr.ztop = fv2.xyz[2] ;
00660
00661 return dr ;
00662 }
00663
00664
00665
00666
00667
00668
00669
00670
00671 void PR_one_slice( int proj_code , MRI_IMAGE * slim , MRI_IMAGE * prim )
00672 {
00673 float * slar = MRI_FLOAT_PTR(slim) ;
00674 float * flar = MRI_FLOAT_PTR(prim) ;
00675 int ii , npix = slim->nvox ;
00676
00677 switch( proj_code ){
00678 default:
00679 case PROJ_SUM:
00680 for( ii=0 ; ii < npix ; ii++ ) flar[ii] += slar[ii] ;
00681 break ;
00682
00683 case PROJ_FIRST:
00684 for( ii=0 ; ii < npix ; ii++ )
00685 if( flar[ii] == 0.0 && slar[ii] > first_thresh )
00686 flar[ii] = slar[ii] ;
00687 break ;
00688
00689 case PROJ_MMAX:
00690 for( ii=0 ; ii < npix ; ii++ )
00691 if( flar[ii] < slar[ii] ) flar[ii] = slar[ii] ;
00692 break ;
00693
00694 case PROJ_AMAX:{
00695 int ss ;
00696 for( ii=0 ; ii < npix ; ii++ ){
00697 ss = abs(slar[ii]) ;
00698 if( flar[ii] < ss ) flar[ii] = ss ;
00699 }
00700 }
00701 break ;
00702
00703 case PROJ_SMAX:{
00704 int ss ;
00705 for( ii=0 ; ii < npix ; ii++ ){
00706 ss = abs(slar[ii]) ;
00707 if( fabs(flar[ii]) < ss ) flar[ii] = slar[ii] ;
00708 }
00709 }
00710 break ;
00711 }
00712 return ;
00713 }
00714
00715
00716
00717
00718
00719
00720 float PR_type_scale( int itype , MRI_IMAGE * prim )
00721 {
00722 float ptop , scl ;
00723
00724 if( ! MRI_IS_INT_TYPE(itype) ) return 1.0 ;
00725
00726 ptop = mri_maxabs( prim ) ;
00727 scl = 1.0 ;
00728
00729 switch( itype ){
00730
00731 default: return scl ;
00732
00733 case MRI_short:
00734 while( ptop > 32767.0 ){
00735 scl /= 10.0 ; ptop /= 10.0 ;
00736 }
00737 return scl ;
00738
00739 case MRI_byte:
00740 while( ptop > 255.0 ){
00741 scl /= 10.0 ; ptop /= 10.0 ;
00742 }
00743 return scl ;
00744
00745 case MRI_int:
00746 while( ptop > 2147483647.0 ){
00747 scl /= 10.0 ; ptop /= 10.0 ;
00748 }
00749 return scl ;
00750 }
00751 }