1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 #if 0
10 # define DB(sss,str) printf("%s %s\n",sss,str)
11 #else
12 # define DB(sss,str)
13 #endif
14
15 static EDIT_options PRED_edopt ;
16
17 /*------------------------------ prototypes ------------------------------*/
18
19 typedef struct {
20 float xbot , xtop , ybot , ytop , zbot , ztop ;
21 } dset_range ;
22
23 dset_range PR_get_range( THD_3dim_dataset * ) ;
24 void PR_one_slice( int , MRI_IMAGE * , MRI_IMAGE * ) ;
25 float PR_type_scale( int , MRI_IMAGE * ) ;
26
27 #define PRED_err(str) \
28 do{ fprintf(stderr,"\n*** %s\n",str) ; exit(-1) ; } while(1)
29
30 /*-------------------------------------------------------------------------*/
31
Syntax(void)32 void Syntax(void)
33 {
34 printf(
35 "Projection along cardinal axes from a 3D dataset\n"
36 "Usage: 3dproject [editing options]\n"
37 " [-sum|-max|-amax|-smax] [-output root] [-nsize] [-mirror]\n"
38 " [-RL {all | x1 x2}] [-AP {all | y1 y2}] [-IS {all | z1 z2}]\n"
39 " [-ALL] dataset\n"
40 "\n"
41 "Program to produce orthogonal projections from a 3D dataset.\n"
42 " -sum ==> Add the dataset voxels along the projection direction\n"
43 " -max ==> Take the maximum of the voxels [the default is -sum]\n"
44 " -amax ==> Take the absolute maximum of the voxels\n"
45 " -smax ==> Take the signed maximum of the voxels; for example,\n"
46 " -max ==> -7 and 2 go to 2 as the projected value\n"
47 " -amax ==> -7 and 2 go to 7 as the projected value\n"
48 " -smax ==> -7 and 2 go to -7 as the projected value\n"
49 " -first x ==> Take the first value greater than x\n"
50 " -nsize ==> Scale the output images up to 'normal' sizes\n"
51 " (e.g., 64x64, 128x128, or 256x256)\n"
52 " This option only applies to byte or short datasets.\n"
53 " -mirror ==> The radiologists' and AFNI convention is to display\n"
54 " axial and coronal images with the subject's left on\n"
55 " the right of the image; the use of this option will\n"
56 " mirror the axial and coronal projections so that\n"
57 " left is left and right is right.\n"
58 "\n"
59 " -output root ==> Output projections will named\n"
60 " root.sag, root.cor, and root.axi\n"
61 " [the default root is 'proj']\n"
62 "\n"
63 " -RL all ==> Project in the Right-to-Left direction along\n"
64 " all the data (produces root.sag)\n"
65 " -RL x1 x2 ==> Project in the Right-to-Left direction from\n"
66 " x-coordinate x1 to x2 (mm)\n"
67 " [negative x is Right, positive x is Left]\n"
68 " [OR, you may use something like -RL 10R 20L\n"
69 " to project from x=-10 mm to x=+20 mm ]\n"
70 "\n"
71 " -AP all ==> Project in the Anterior-to-Posterior direction along\n"
72 " all the data (produces root.cor)\n"
73 " -AP y1 y2 ==> Project in the Anterior-to-Posterior direction from\n"
74 " y-coordinate y1 to y2 (mm)\n"
75 " [negative y is Anterior, positive y is Posterior]\n"
76 " [OR, you may use something like -AP 10A 20P\n"
77 " to project from y=-10 mm to y=+20 mm ]\n"
78 "\n"
79 " -IS all ==> Project in the Inferior-to-Superior direction along\n"
80 " all the data (produces root.axi)\n"
81 " -IS y1 y2 ==> Project in the Inferior-to-Superior direction from\n"
82 " z-coordinate z1 to z2 (mm)\n"
83 " [negative z is Inferior, positive z is Superior]\n"
84 " [OR, you may use something like -IS 10I 20S\n"
85 " to project from z=-10 mm to z=+20 mm ]\n"
86 "\n"
87 " -ALL ==> Equivalent to '-RL all -AP all -IS all'\n"
88 "\n"
89 "* NOTE that a projection direction will not be used if the bounds aren't\n"
90 " given for that direction; thus, at least one of -RL, -AP, or -IS must\n"
91 " be used, or nothing will be computed!\n"
92 "* NOTE that in the directions transverse to the projection direction,\n"
93 " all the data is used; that is, '-RL -5 5' will produce a full sagittal\n"
94 " image summed over a 10 mm slice, irrespective of the -IS or -AP extents.\n"
95 "* NOTE that the [editing options] are the same as in 3dmerge.\n"
96 " In particular, the '-1thtoin' option can be used to project the\n"
97 " threshold data (if available).\n"
98 ) ;
99
100 printf("\n" MASTER_SHORTHELP_STRING ) ;
101
102 PRINT_COMPILE_DATE ; exit(0) ;
103 }
104
105 #define PROJ_SUM 0
106 #define PROJ_MMAX 1
107 #define PROJ_AMAX 2
108 #define PROJ_SMAX 3
109 #define PROJ_FIRST 4 /* 02 Nov 2000 */
110
111 #define MIRR_NO 0
112 #define MIRR_YES 1
113
114 #define BIGG 9999999.99
115
116 static float first_thresh=0.0 ;
117
main(int argc,char * argv[])118 int main( int argc , char * argv[] )
119 {
120 THD_3dim_dataset * dset ;
121 THD_dataxes * daxes ;
122 FD_brick ** brarr , * baxi , * bsag , * bcor ;
123
124 int iarg , ii ;
125 MRI_IMAGE * pim , * flim , * slim ;
126 float * flar ;
127 dset_range dr ;
128 float val , fimfac ;
129 int ival,ityp , kk ;
130
131 float xbot=BIGG,xtop=BIGG , ybot=BIGG,ytop=BIGG , zbot=BIGG,ztop=BIGG ;
132 int xgood=0 , ygood=0 , zgood=0 ,
133 proj_code=PROJ_SUM , mirror_code=MIRR_NO , nsize=0 ;
134 char root[THD_MAX_NAME] = "proj." ;
135 char fname[THD_MAX_NAME] ;
136
137 int ixbot=0,ixtop=0 , jybot=0,jytop=0 , kzbot=0,kztop=0 ;
138 THD_fvec3 fv ;
139 THD_ivec3 iv ;
140
141 /*--- read command line arguments ---*/
142
143 WARNING_message("This program (3dproject) is old, not maintained, and probably useless!") ;
144
145 if( argc < 2 || strncmp(argv[1],"-help",6) == 0 ) Syntax() ;
146
147 INIT_EDOPT( &PRED_edopt ) ;
148 iarg = 1 ;
149 while( iarg < argc && argv[iarg][0] == '-' ){
150
151 DB("new arg:",argv[iarg]) ;
152
153 /**** check for editing option ****/
154
155 ii = EDIT_check_argv( argc , argv , iarg , &PRED_edopt ) ;
156 if( ii > 0 ){
157 iarg += ii ;
158 continue ;
159 }
160
161 /**** -sum or -max ****/
162
163 if( strncmp(argv[iarg],"-sum",6) == 0 ){
164 proj_code = PROJ_SUM ;
165 iarg++ ; continue ;
166 }
167
168 if( strncmp(argv[iarg],"-max",6) == 0 ){
169 proj_code = PROJ_MMAX ;
170 iarg++ ; continue ;
171 }
172
173 if( strncmp(argv[iarg],"-amax",6) == 0 ){
174 proj_code = PROJ_AMAX ;
175 iarg++ ; continue ;
176 }
177
178 if( strncmp(argv[iarg],"-smax",6) == 0 ){
179 proj_code = PROJ_SMAX ;
180 iarg++ ; continue ;
181 }
182
183 if( strcmp(argv[iarg],"-first") == 0 ){ /* 02 Nov 2000 */
184 proj_code = PROJ_FIRST ;
185 first_thresh = strtod( argv[++iarg] , NULL ) ;
186 iarg++ ; continue ;
187 }
188
189 /**** -mirror ****/
190
191 if( strncmp(argv[iarg],"-mirror",6) == 0 ){
192 mirror_code = MIRR_YES ;
193 iarg++ ; continue ;
194 }
195
196 /**** -nsize ****/
197
198 if( strncmp(argv[iarg],"-nsize",6) == 0 ){
199 nsize = 1 ;
200 iarg++ ; continue ;
201 }
202
203 /**** -output root ****/
204
205 if( strncmp(argv[iarg],"-output",6) == 0 ||
206 strncmp(argv[iarg],"-root",6) == 0 ){
207
208 if( iarg+1 >= argc ){
209 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
210 exit(-1) ;
211 }
212
213 MCW_strncpy( root , argv[++iarg] , THD_MAX_NAME-1 ) ;
214 ii = strlen(root) ;
215 if( ii == 0 ){
216 fprintf(stderr,"\n*** illegal rootname!\n") ; exit(-1) ;
217 }
218 if( root[ii-1] != '.' ){ root[ii] = '.' ; root[ii+1] = '\0' ; }
219 iarg++ ; continue ;
220 }
221
222 /**** -ALL ****/
223
224 if( strncmp(argv[iarg],"-ALL",6) == 0 ||
225 strncmp(argv[iarg],"-all",6) == 0 ){
226
227 xgood = ygood = zgood = 1 ;
228 xbot = ybot = zbot = -BIGG ;
229 xtop = ytop = ztop = BIGG ;
230 iarg++ ; continue ;
231 }
232
233 /**** -RL {all | x1 x2} ****/
234
235 if( strncmp(argv[iarg],"-RL",6) == 0 ||
236 strncmp(argv[iarg],"-LR",6) == 0 ||
237 strncmp(argv[iarg],"-rl",6) == 0 ||
238 strncmp(argv[iarg],"-lr",6) == 0 ||
239 strncmp(argv[iarg],"-sag",6)== 0 ){
240
241 char * cerr ; float tf ;
242
243 xgood = 1 ; /* mark for x projection */
244
245 if( iarg+1 >= argc ){
246 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
247 exit(-1) ;
248 }
249
250 if( strncmp(argv[iarg+1],"all",6) == 0 ||
251 strncmp(argv[iarg+1],"ALL",6) == 0 ){
252
253 xbot = -BIGG;
254 xtop = BIGG ;
255 iarg += 2 ; continue ;
256 }
257
258 if( iarg+2 >= argc ){
259 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
260 exit(-1) ;
261 }
262
263 xbot = strtod( argv[iarg+1] , &cerr ) ;
264 if( cerr == argv[iarg+1] ){
265 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
266 argv[iarg],argv[iarg+1] ) ; exit(-1) ;
267 }
268 if( *cerr == 'R' && xbot > 0.0 ) xbot = -xbot ;
269
270 xtop = strtod( argv[iarg+2] , &cerr ) ;
271 if( cerr == argv[iarg+2] ){
272 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
273 argv[iarg],argv[iarg+2] ) ; exit(-1) ;
274 }
275 if( *cerr == 'R' && xtop > 0.0 ) xtop = -xtop ;
276
277 if( xbot > xtop ){ tf = xbot ; xbot = xtop ; xtop = tf ; }
278 iarg +=3 ; continue ;
279 }
280
281 /**** -AP {all | y1 y2} ****/
282
283 if( strncmp(argv[iarg],"-AP",6) == 0 ||
284 strncmp(argv[iarg],"-PA",6) == 0 ||
285 strncmp(argv[iarg],"-ap",6) == 0 ||
286 strncmp(argv[iarg],"-pa",6) == 0 ||
287 strncmp(argv[iarg],"-cor",6)== 0 ){
288
289 char * cerr ; float tf ;
290
291 ygood = 1 ; /* mark for y projection */
292
293 if( iarg+1 >= argc ){
294 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
295 exit(-1) ;
296 }
297
298 if( strncmp(argv[iarg+1],"all",6) == 0 ||
299 strncmp(argv[iarg+1],"ALL",6) == 0 ){
300
301 ybot = -BIGG ;
302 ytop = BIGG ;
303 iarg += 2 ; continue ;
304 }
305
306 if( iarg+2 >= argc ){
307 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
308 exit(-1) ;
309 }
310
311 ybot = strtod( argv[iarg+1] , &cerr ) ;
312 if( cerr == argv[iarg+1] ){
313 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
314 argv[iarg],argv[iarg+1] ) ; exit(-1) ;
315 }
316 if( *cerr == 'A' && ybot > 0.0 ) ybot = -ybot ;
317
318 ytop = strtod( argv[iarg+2] , &cerr ) ;
319 if( cerr == argv[iarg+2] ){
320 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
321 argv[iarg],argv[iarg+2] ) ; exit(-1) ;
322 }
323 if( *cerr == 'A' && ytop > 0.0 ) ytop = -ytop ;
324
325 if( ybot > ytop ){ tf = ybot ; ybot = ytop ; ytop = tf ; }
326 iarg +=3 ; continue ;
327 }
328
329 /**** -IS {all | z1 z2} ****/
330
331 if( strncmp(argv[iarg],"-IS",6) == 0 ||
332 strncmp(argv[iarg],"-SI",6) == 0 ||
333 strncmp(argv[iarg],"-is",6) == 0 ||
334 strncmp(argv[iarg],"-si",6) == 0 ||
335 strncmp(argv[iarg],"-axi",6)== 0 ){
336
337 char * cerr ; float tf ;
338
339 zgood = 1 ; /* mark for y projection */
340
341 if( iarg+1 >= argc ){
342 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
343 exit(-1) ;
344 }
345
346 if( strncmp(argv[iarg+1],"all",6) == 0 ||
347 strncmp(argv[iarg+1],"ALL",6) == 0 ){
348
349 zbot = -BIGG ;
350 ztop = BIGG ;
351 iarg += 2 ; continue ;
352 }
353
354 if( iarg+2 >= argc ){
355 fprintf(stderr,"\n*** no argument after option %s\n",argv[iarg]) ;
356 exit(-1) ;
357 }
358
359 zbot = strtod( argv[iarg+1] , &cerr ) ;
360 if( cerr == argv[iarg+1] ){
361 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
362 argv[iarg],argv[iarg+1] ) ; exit(-1) ;
363 }
364 if( *cerr == 'I' && zbot > 0.0 ) zbot = -zbot ;
365
366 ztop = strtod( argv[iarg+2] , &cerr ) ;
367 if( cerr == argv[iarg+2] ){
368 fprintf(stderr,"\n*** illegal argument after %s: %s\n",
369 argv[iarg],argv[iarg+2] ) ; exit(-1) ;
370 }
371 if( *cerr == 'I' && ztop > 0.0 ) ztop = -ztop ;
372
373 if( zbot > ztop ){ tf = zbot ; zbot = ztop ; ztop = tf ; }
374 iarg +=3 ; continue ;
375 }
376
377 /**** unknown option ****/
378
379 fprintf(stderr,"\n*** Unknown option: %s\n",argv[iarg]) ;
380 exit(-1) ;
381 } /* end of loop over input options */
382
383 if( ! xgood && ! ygood && ! zgood ){
384 fprintf(stderr,"\n*** No projections ordered!?\n") ; exit(-1) ;
385 }
386
387 /*--- open dataset and set up to extract data slices ---*/
388
389 dset = THD_open_dataset( argv[iarg] ) ;
390 if( dset == NULL ){
391 fprintf(stderr,"\n*** Can't open dataset file %s\n",argv[iarg]) ;
392 exit(-1) ;
393 }
394 if( DSET_NUM_TIMES(dset) > 1 ){
395 fprintf(stderr,"\n*** Can't project time-dependent dataset!\n") ;
396 exit(1) ;
397 }
398 EDIT_one_dataset( dset, &PRED_edopt ) ;
399 daxes = dset->daxes ;
400 brarr = THD_setup_bricks( dset ) ;
401 baxi = brarr[0] ; bsag = brarr[1] ; bcor = brarr[2] ;
402
403 /*--- determine index range for each direction ---*/
404
405 dr = PR_get_range( dset ) ;
406
407 if( xgood ){
408 if( xbot < dr.xbot ) xbot = dr.xbot ;
409 if( xtop > dr.xtop ) xtop = dr.xtop ;
410
411 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(xbot,0,0) ) ;
412 iv = THD_3dmm_to_3dind ( dset , fv ) ;
413 iv = THD_3dind_to_fdind( bsag , iv ) ; ixbot = iv.ijk[2] ;
414
415 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(xtop,0,0) ) ;
416 iv = THD_3dmm_to_3dind ( dset , fv ) ;
417 iv = THD_3dind_to_fdind( bsag , iv ) ; ixtop = iv.ijk[2] ;
418
419 if( ixbot > ixtop ) { ii = ixbot ; ixbot = ixtop ; ixtop = ii ; }
420
421 if( ixbot < 0 ) ixbot = 0 ;
422 if( ixtop >= bsag->n3 ) ixtop = bsag->n3 - 1 ;
423 }
424
425 if( ygood ){
426 if( ybot < dr.ybot ) ybot = dr.ybot ;
427 if( ytop > dr.ytop ) ytop = dr.ytop ;
428
429 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(0,ybot,0) ) ;
430 iv = THD_3dmm_to_3dind ( dset , fv ) ;
431 iv = THD_3dind_to_fdind( bcor , iv ) ; jybot = iv.ijk[2] ;
432
433 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(0,ytop,0) ) ;
434 iv = THD_3dmm_to_3dind ( dset , fv ) ;
435 iv = THD_3dind_to_fdind( bcor , iv ) ; jytop = iv.ijk[2] ;
436
437 if( jybot > jytop ) { ii = jybot ; jybot = jytop ; jytop = ii ; }
438
439 if( jybot < 0 ) jybot = 0 ;
440 if( jytop >= bcor->n3 ) jytop = bcor->n3 - 1 ;
441 }
442
443 if( zgood ){
444 if( zbot < dr.zbot ) zbot = dr.zbot ;
445 if( ztop > dr.ztop ) ztop = dr.ztop ;
446
447 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(0,0,zbot) ) ;
448 iv = THD_3dmm_to_3dind ( dset , fv ) ;
449 iv = THD_3dind_to_fdind( baxi , iv ) ; kzbot = iv.ijk[2] ;
450
451 fv = THD_dicomm_to_3dmm( dset , TEMP_FVEC3(0,0,ztop) ) ;
452 iv = THD_3dmm_to_3dind ( dset , fv ) ;
453 iv = THD_3dind_to_fdind( baxi , iv ) ; kztop = iv.ijk[2] ;
454
455 if( kzbot > kztop ) { ii = kzbot ; kzbot = kztop ; kztop = ii ; }
456
457 if( kzbot < 0 ) kzbot = 0 ;
458 if( kztop >= baxi->n3 ) kztop = baxi->n3 - 1 ;
459 }
460
461 ival = DSET_PRINCIPAL_VALUE(dset) ; /* index to project */
462 ityp = DSET_BRICK_TYPE(dset,ival) ; /* type of this data */
463 fimfac = DSET_BRICK_FACTOR(dset,ival) ; /* scale factor of this data */
464
465 /*--- project the x direction, if desired ---*/
466
467 if( xgood ){
468 int n1 = bsag->n1 , n2 = bsag->n2 ;
469 int ss , npix ;
470 float fmax , fmin , scl ;
471
472 /*-- set up --*/
473
474 npix = n1*n2 ;
475 flim = mri_new( n1 , n2 , MRI_float ) ;
476 flar = mri_data_pointer( flim ) ;
477 for( ii=0 ; ii < npix ; ii++ ) flar[ii] = 0.0 ;
478
479 /*-- actually project --*/
480
481 for( ss=ixbot ; ss <= ixtop ; ss++ ){
482 slim = FD_brick_to_mri( ss , ival , bsag ) ;
483 if( slim->kind != MRI_float ){
484 pim = mri_to_float( slim ) ;
485 mri_free( slim ) ; slim = pim ;
486 }
487 PR_one_slice( proj_code , slim , flim ) ;
488 mri_free( slim ) ;
489 }
490
491 /*-- form output --*/
492
493 if( fimfac != 0.0 && fimfac != 1.0 ){
494 slim = mri_scale_to_float( 1.0/fimfac , flim ) ;
495 mri_free(flim) ; flim = slim ;
496 }
497
498 scl = PR_type_scale( ityp , flim ) ;
499 pim = mri_to_mri_scl( ityp , scl , flim ) ; mri_free( flim ) ;
500 if( nsize ){
501 slim = mri_nsize( pim ) ;
502 if( slim != NULL && slim != pim ) { mri_free(pim) ; pim = slim ; }
503 }
504
505 if( scl != 1.0 )
506 printf("Sagittal projection pixels scaled by %g to avoid overflow!\n",
507 scl ) ;
508
509 fmax = mri_max(pim) ; fmin = mri_min(pim) ;
510 printf("Sagittal projection min = %g max = %g\n",fmin,fmax) ;
511
512 strcpy(fname,root) ; strcat(fname,"sag") ;
513 mri_write( fname, pim ) ;
514 mri_free(pim) ;
515 }
516
517 /*--- project the y direction, if desired ---*/
518
519 if( ygood ){
520 int n1 = bcor->n1 , n2 = bcor->n2 ;
521 int ss , npix ;
522 float fmax , fmin , scl ;
523
524 /*-- set up --*/
525
526 npix = n1*n2 ;
527 flim = mri_new( n1 , n2 , MRI_float ) ;
528 flar = mri_data_pointer( flim ) ;
529 for( ii=0 ; ii < npix ; ii++ ) flar[ii] = 0.0 ;
530
531 /*-- actually project --*/
532
533 for( ss=jybot ; ss <= jytop ; ss++ ){
534 slim = FD_brick_to_mri( ss , ival , bcor ) ;
535 if( slim->kind != MRI_float ){
536 pim = mri_to_float( slim ) ;
537 mri_free( slim ) ; slim = pim ;
538 }
539 PR_one_slice( proj_code , slim , flim ) ;
540 mri_free( slim ) ;
541 }
542
543 /*-- form output --*/
544
545 if( fimfac != 0.0 && fimfac != 1.0 ){
546 slim = mri_scale_to_float( 1.0/fimfac , flim ) ;
547 mri_free(flim) ; flim = slim ;
548 }
549
550 scl = PR_type_scale( ityp , flim ) ;
551 pim = mri_to_mri_scl( ityp , scl , flim ) ; mri_free( flim ) ;
552 if( nsize ){
553 slim = mri_nsize( pim ) ;
554 if( slim != NULL && slim != pim ) { mri_free(pim) ; pim = slim ; }
555 }
556
557 if( scl != 1.0 )
558 printf("Coronal projection pixels scaled by %g to avoid overflow!\n",
559 scl ) ;
560
561 fmax = mri_max(pim) ; fmin = mri_min(pim) ;
562 printf("Coronal projection min = %g max = %g\n",fmin,fmax) ;
563
564 if( mirror_code == MIRR_YES ){
565 slim = mri_flippo( MRI_ROT_0 , TRUE , pim ) ;
566 mri_free(pim) ; pim = slim ;
567 }
568 strcpy(fname,root) ; strcat(fname,"cor") ;
569 mri_write( fname, pim ) ;
570 mri_free(pim) ;
571 }
572
573 /*--- project the z direction, if desired ---*/
574
575 if( zgood ){
576 int n1 = baxi->n1 , n2 = baxi->n2 ;
577 int ss , npix ;
578 float fmax , fmin , scl ;
579
580 /*-- set up --*/
581
582 npix = n1*n2 ;
583 flim = mri_new( n1 , n2 , MRI_float ) ;
584 flar = mri_data_pointer( flim ) ;
585 for( ii=0 ; ii < npix ; ii++ ) flar[ii] = 0.0 ;
586
587 /*-- actually project --*/
588
589 for( ss=kzbot ; ss <= kztop ; ss++ ){
590 slim = FD_brick_to_mri( ss , ival , baxi ) ;
591 if( slim->kind != MRI_float ){
592 pim = mri_to_float( slim ) ;
593 mri_free( slim ) ; slim = pim ;
594 }
595 PR_one_slice( proj_code , slim , flim ) ;
596 mri_free( slim ) ;
597 }
598
599 /*-- form output --*/
600
601 if( fimfac != 0.0 && fimfac != 1.0 ){
602 slim = mri_scale_to_float( 1.0/fimfac , flim ) ;
603 mri_free(flim) ; flim = slim ;
604 }
605
606 scl = PR_type_scale( ityp , flim ) ;
607 pim = mri_to_mri_scl( ityp , scl , flim ) ; mri_free( flim ) ;
608 if( nsize ){
609 slim = mri_nsize( pim ) ;
610 if( slim != NULL && slim != pim ) { mri_free(pim) ; pim = slim ; }
611 }
612
613 if( scl != 1.0 )
614 printf("Axial projection pixels scaled by %g to avoid overflow!\n",
615 scl ) ;
616
617 fmax = mri_max(pim) ; fmin = mri_min(pim) ;
618 printf("Axial projection min = %g max = %g\n",fmin,fmax) ;
619
620 if( mirror_code == MIRR_YES ){
621 slim = mri_flippo( MRI_ROT_0 , TRUE , pim ) ;
622 mri_free(pim) ; pim = slim ;
623 }
624 strcpy(fname,root) ; strcat(fname,"axi") ;
625 mri_write( fname, pim ) ;
626 mri_free(pim) ;
627 }
628
629 exit(0) ;
630 }
631
632 /*------------------------------------------------------------------*/
633
PR_get_range(THD_3dim_dataset * dset)634 dset_range PR_get_range( THD_3dim_dataset * dset )
635 {
636 THD_dataxes * daxes ;
637 THD_fvec3 fv1 , fv2 , fv3 ;
638 THD_ivec3 iv ;
639 float tf ;
640 dset_range dr ;
641
642 #define FSWAP(x,y) (tf=(x),(x)=(y),(y)=tf)
643
644 daxes = dset->daxes ;
645
646 LOAD_FVEC3(fv1 , daxes->xxorg , daxes->yyorg , daxes->zzorg) ;
647 fv1 = THD_3dmm_to_dicomm( dset , fv1 ) ;
648
649 LOAD_FVEC3(fv2 , daxes->xxorg + (daxes->nxx-1)*daxes->xxdel ,
650 daxes->yyorg + (daxes->nyy-1)*daxes->yydel ,
651 daxes->zzorg + (daxes->nzz-1)*daxes->zzdel ) ;
652 fv2 = THD_3dmm_to_dicomm( dset , fv2 ) ;
653
654 if( fv1.xyz[0] > fv2.xyz[0] ) FSWAP( fv1.xyz[0] , fv2.xyz[0] ) ;
655 if( fv1.xyz[1] > fv2.xyz[1] ) FSWAP( fv1.xyz[1] , fv2.xyz[1] ) ;
656 if( fv1.xyz[2] > fv2.xyz[2] ) FSWAP( fv1.xyz[2] , fv2.xyz[2] ) ;
657
658 dr.xbot = fv1.xyz[0] ; dr.xtop = fv2.xyz[0] ;
659 dr.ybot = fv1.xyz[1] ; dr.ytop = fv2.xyz[1] ;
660 dr.zbot = fv1.xyz[2] ; dr.ztop = fv2.xyz[2] ;
661
662 return dr ;
663 }
664
665 /**********************************************************************
666 Add one slice to the projection.
667 slim = input image from brick (both images are MRI_float)
668 prim = accumulating projection image
669 proj_code = one of the PROJ_* values
670 **********************************************************************/
671
PR_one_slice(int proj_code,MRI_IMAGE * slim,MRI_IMAGE * prim)672 void PR_one_slice( int proj_code , MRI_IMAGE * slim , MRI_IMAGE * prim )
673 {
674 float * slar = MRI_FLOAT_PTR(slim) ;
675 float * flar = MRI_FLOAT_PTR(prim) ;
676 int ii , npix = slim->nvox ;
677
678 switch( proj_code ){
679 default:
680 case PROJ_SUM:
681 for( ii=0 ; ii < npix ; ii++ ) flar[ii] += slar[ii] ;
682 break ;
683
684 case PROJ_FIRST: /* 02 Nov 2000 */
685 for( ii=0 ; ii < npix ; ii++ )
686 if( flar[ii] == 0.0 && slar[ii] > first_thresh )
687 flar[ii] = slar[ii] ;
688 break ;
689
690 case PROJ_MMAX:
691 for( ii=0 ; ii < npix ; ii++ )
692 if( flar[ii] < slar[ii] ) flar[ii] = slar[ii] ;
693 break ;
694
695 case PROJ_AMAX:{
696 int ss ;
697 for( ii=0 ; ii < npix ; ii++ ){
698 ss = abs(slar[ii]) ;
699 if( flar[ii] < ss ) flar[ii] = ss ;
700 }
701 }
702 break ;
703
704 case PROJ_SMAX:{
705 int ss ;
706 for( ii=0 ; ii < npix ; ii++ ){
707 ss = abs(slar[ii]) ;
708 if( fabs(flar[ii]) < ss ) flar[ii] = slar[ii] ;
709 }
710 }
711 break ;
712 }
713 return ;
714 }
715
716 /*--------------------------------------------------------------
717 Find the scale factor needed to convert the MRI_float input
718 image to the desired output type.
719 ----------------------------------------------------------------*/
720
PR_type_scale(int itype,MRI_IMAGE * prim)721 float PR_type_scale( int itype , MRI_IMAGE * prim )
722 {
723 float ptop , scl ;
724
725 if( ! MRI_IS_INT_TYPE(itype) ) return 1.0 ;
726
727 ptop = mri_maxabs( prim ) ;
728 scl = 1.0 ;
729
730 switch( itype ){
731
732 default: return scl ;
733
734 case MRI_short:
735 while( ptop > 32767.0 ){
736 scl /= 10.0 ; ptop /= 10.0 ;
737 }
738 return scl ;
739
740 case MRI_byte:
741 while( ptop > 255.0 ){
742 scl /= 10.0 ; ptop /= 10.0 ;
743 }
744 return scl ;
745
746 case MRI_int:
747 while( ptop > 2147483647.0 ){
748 scl /= 10.0 ; ptop /= 10.0 ;
749 }
750 return scl ;
751 }
752 }
753